1 Introduction
Contents
1. dx IL Setup problem and compute solution problem VariationalProblem a L bcs problem parameters linear_solver direct U problem solve Split the mixed solution using deep copy needed for further computation on coefficient vector u p U split True u_plot plot u p_plot plot p print Norm of velocity coefficient vector 4e u vector norm 12 print Norm of pressure coefficient vector 4e p vector norm 12 Alternative assembly process with more than one Dianuchilet boundary condation samply wsie A b assemble_system a L bcs Or boe oLlormimna s assemble a assemble L or condition in bcs condition apply A b hh oO ff HH KH FH Solve the manually assembled problem U2 Function w solve A U2 vector b Compute the difference should be zero of velocities u2 p2 U2 split True diff u vector u2 vector print Norm of velocity difference 4e diff norm 12 UD 12 ROLAND HERZOG Extract the boundary mesh boundary_mesh BoundaryMesh mesh When the boundary mesh is created we also get two maps attached to the boundary mesh as mesh data named cell map and vertex map They are mesh functions which map the cells and vertices of the boundary mesh to their counterparts in the original mesh cell_map boundary_mesh data mesh_function cell map Now we can create a plottable version of the subd2. differentiating L This is achieved by the statement a derivative L u du in the code below which solves the semilinear problem Aut w f in u 0 onr INTRODUCTION TO FENICS 13 from dolfin import HecCucave mesh wand detane mune talon space mesh UnitSquare 32 32 V FunctionSpace mesh Lagrange 2 Define Dirichlet boundary everywhere class DirichletBoundary SubDomain def inside self x on_boundary return on_boundary Define boundary condition u0 Constant 0 0 be DirichletBC V u0 DirichletBoundary Define the source and solution functions f Expression 500 0 exp pow x 0 0 5 2 pow x SOR SR 2 70 02 2 u Function V Define the variational problem which represents the Newton step v TestFunction V du TrialFunction V Wem mee det Omer ine mon lh eat hem mon inean esnia L inner grad v grad u dx v ux 3 dx v f dx since the bilinear form can be created automatically i OS We ClewiWeNelye Os lle weSacliell Ik fi Ww a derivative L u du Setup the nonlinear problem L u 0 problem VariationalProblem a L bc nonlinear True and solve it NOTE that u problem solve DOES NOT WORK no progress in iterations problem solve u Note that the high level command problem solve u calls the Newton solver and the assembly of the Jacobian and the residual is triggered automatically by the solver The following code shows how to achieve
3. 00 1 00 1 50 2 00 In the present situation we follow one of the examples that came with DOLFIN and consider the following boundary conditions for the Stokes system u 0O onT 9 no slip boundary u uog onl inflow boundary p 0 onT outflow boundary These three parts of the boundary are specified by the values 0 1 2 of a mesh function In this setting the value of the pressure is uniquely defined i e the FE stiffness matrix is invertible from dolfin import Load mesh and subdomains mesh Mesh dolfin 2 xml gz subdomains MeshFunction uint mesh subdomains xml gz Define function spaces V VectorFunctionSpace mesh Lagrange 2 Q FunctionSpace mesh Lagrange 1 W V Q No slip boundary condition for velocity bc_noslip Constant 0 0 bcO DirichletBC W sub 0 bc_noslip subdomains 0 Inflow boundary condition for velocity bc_inflow Expression sin x 1 pi 0 0 bel DirichletBC W sub 0 bc_inflow subdomains 1 Boundary condition for pressure at outflow bc_zero Constant 0 0 INTRODUCTION TO FEniCS 11 bc2 DirichletBC W sub 1 bc_zero subdomains 2 Collect boundary conditions bes bcO Del bc2 Define variational problem v q TestFunctions Ww u p TrialFunctions W f Constant 0 0 mu Constant 1 0e 3 n FacetNormal mesh a muxinner grad v grad u div v p q div u dx p dot v n ds inner v f
4. INTRODUCTION TO THE FINITE ELEMENT SOFTWARE FEniCS ROLAND HERZOG 1 Introduction FENICS is free software for automated solution of differential equations As a particular feature it allows the formulation of the weak variational forms in a rather natural way For instance the bilinear and linear forms associated to the Poisson problem Find ueV such that a u v L v forallv eV with aluo f Vv Vuda bw vfds can be formulated as a inner grad v grad u dx L v i dx The language used here is the so called Unified Form Language UFL which is an extension of the Python scripting language 1 1 Components in the FEniCS Project The interaction of components within the FENICS project is reflected in Figure 1 1 The core library is DOLFIN which is a library written in C with interfaces to C and Python as a mod ule It provides functionalities such as mesh handling refinement coarsening smoothing partitioning mesh generation currently only for simple shapes assembly of linear and bilinear in general multilinear forms input output routines Linear algebra and visualization is handled by external software The typical usage of the library within a C program is like this demo ufl demo h C compiler demo cpp parei demo In this case the form compiler FFC or its alternative SYF1 SFC generates a header file demo h which defines a number of proble
5. ce tern class Source public Expression publac Source Expression void eval Array lt double gt amp values const Array lt double gt amp x const double dx x 0 0 5 double dy x 1 0 5 values 0 500 0 exp dx dx dy dy 0 02 Jr r3 Subomain for Dirichlet boundary condition class DirichletBoundary public SubDomain bool inside const Array lt double gt amp x bool on_boundary const return x 0 lt DOLFIN_EPS or x 0 gt 1 0 DOLFIN_EPS int main al Create mesh and function space UnitSquare mesh 32 32 Example_Poisson FunctionSpace V mesh Define boundary condition Constant u0 0 0 DirichletBoundary boundary DirichletBC bc V u0 boundary INTRODUCTION TO FENICS 5 Define variational problem Example_Poisson BilinearForm a V V Example_Poisson LinearForm L V Source f aoe wt cae Setup problem and compute solution VariationalProblem problem a L bc Function u V problem parameters linear_solver iterative problem solve u Save solution in VTK format File solution_file poisson_solution pvd solution_file lt lt u Save the mesh also File mesh_file poisson_mesh xml mesh_file lt lt mesh Save the parameter file File parameter_file poisson_parameters xml parameter_file lt lt problem parameters Plot the mesh and the solution using viper plot mesh plot u r
6. ctor valued and compound finite elements velocity pressure e implement different boundary conditions for different components of a system e export finite element matrices e import finite element meshes We consider the Stokes system in its strong form pAut Vp f inQ divu 0 inQ where u represents the velocity of a fluid with dynamic viscosity u and p is its pressure Suitable boundary conditions are for example u ug on I which have to satisfy the compatibility condition Se uo n 0 Note that the pressure is only determined up to a constant factor in this setting We take as a first example Q 0 1 with uo 1 0 on the upper edge uo 0 0 elsewhere and f 0 0 This problem is known as a lid driven cavity flow The weak formulation of the Stokes problem is Find u p wo V x Q such that u vo Vud divvpde py ads f v far Q Q r ea Q q div u dz 0 Q for all v V v H Q v OonT and all q Q q P Q Jaade 0 L Q R The problem has a saddle point structure from dolfin import Create mesh and define function spaces Taylor Hood FE mesh UnitSquare 32 32 V VectorFunctionSpace mesh Lagrange 2 8 ROLAND HERZOG Q W FunctionSpace mesh Lagrange 1 V Q Define Dirichlet boundary everywhere class DirichletBoundary SubDomain def inside self x on_boundary return on_boundary Define Dirichlet boundary condition every
7. de after the script has been executed Changing the Boundary Conditions We now change the problem by modify ing the boundary conditions We use u sin a z2 on Tp Ou On au g with a x 1 and g x 10 x on the complement I y Moreover we replace Au by div AVu with A x 6 a The variational problem then becomes Find u ug V such that yea vuae avuds vfdz vgds Q TN Q TN for all v V v H Q v 0 on Tp Here are the relevant changes to the Python code above Define boundary condition u0 Expression sin x 1 pi be DirichletBC V u0 DirichletBoundary Define variational problem TrialFunction V TestFunction V Expression 500 0 exp pow x 0 0 5 2 pow x 1 8 D 02 h a E H INTRODUCTION TO FENICS T g Expression 10 0 x 0 A as_matrix 2 0 1 0 1 0 2 0 alpha Constant 1 0 a inner grad v A grad u dx alpha v u ds Suggestions Change the order of the finite element Change the setting to 3D Try in Python help u princ unr mCO 5 0 5 3 ce MOE u clememeOQs wu vector array 0 10 help V print V V dim type v help mesh print mesh mesh coordinates 0 10 mesh cells mesh hmax mesh num_entities 0 in the interactive Python shell Try also the command line tool ufl convert 2 2 The Stokes System Purpose of this example e introduce ve
8. eturn Q The we invoke the build tool CMake to obtain the executable Example_Poisson run it and visualize the solution with ParaView cmake amp amp make Example_Poisson paraview data poisson_solution pvd Here is a file with the same functionality in Python from dolfin import Create mesh and define function space mesh UnitSquare 32 32 V FunctionSpace mesh Lagrange 1 Define boundary condition x 0 or x 1 u0 Constant 0 0 be DirichletBC V u0 x 0 lt DOLFIN_EPS x 0 gt 1 0 DOLFIN_EPS Define variational problem u TrialFunction V 6 ROLAND HERZOG v TestFunction V f Expression 500 0 exp pow x 0 0 5 2 pow x 1 0 55 2 0 02 1 a inner grad u grad v dx L f v dx Setup problem and compute solution problem VariationalProblem a L bc problem parameters linear_solver iterative u problem solve Save solution in VTK format solution_file File poisson_solution pvd solution_file lt lt u Save the mesh also mesh_file File poisson_mesh xml mesh_file lt lt mesh Save the parameter file parameter_file File poisson_parameters xml parameter_file lt lt problem parameters Plot the mesh and the solution using viper plot mesh plot u interactive True Execute it by saying python i Example_Poisson py where the switch i stands for interactive mo
9. ghtly 2 1 The Poisson Equation Purpose of this example e solve Poisson equation with Lagrange FEs of various degrees present C and Python interfaces visualize the results implement various choices of boundary conditions and non constant coef ficients We consider the Poisson equation in its variational form i e find u V such that ve vude ofa Q Q for all v V v H Q v 0 on Tp where Ip is the Dirichlet part of the boundary I 0Q In its strong form this problem corresponds to Au f in Q with boundary conditions u 0 on Ip and 0u 0n 0 on the complement Ty The corresponding variational description in UFL looks like this Example_Poisson uf1 element FiniteElement Lagrange triangle 1 v TestFunction element u TrialFunction element f Coefficient element 4 ROLAND HERZOG a inner grad v grad u dx v f dx ll We compile it using FFC and obtain a header file Example_Poisson h ffc l1 dolfin Example_Poisson ufl The following main program makes use of the general DOLFIN classes Function SubDomain UnitSquare Constant DirichletBC VariationalProblem and File It also uses the classes FunctionSpace BilinearForm and LinearForm which have been defined in the header file Example_Poisson h automatically generated by FFC include lt dolfin h gt include Example_Poisson h using namespace dolfin Sour
10. ini elements of arbitrary degree INTRODUCTION TO FENICS 3 Brezzi Douglas Fortin Marini elements of arbitrary degree N d lec elements 1st kind of arbitrary degree N d lec elements 2nd kind of arbitrary degree e H curl conforming FE spaces composed of N d lec elements 1st kind of arbitrary degree N d lec elements 2nd kind of arbitrary degree e L conforming FE spaces composed of discontinuous Lagrange elements of arbitrary degree Crouzeix Raviart elements degree one All of these can be used on simplicial meshes in 1D 2D and 3D intervals triangles tetrahedra but currently not on other geometries rectangles cubes prisms etc The variational form language is described in more detail in Section 3 It supports in particular e scalar vector valued and matrix valued finite elements e integrals over cells interior facets and exterior facets but not over arbitrary lines or surfaces e jumps and averages over facets which can be used for error estimation and in discontinuous Galerkin DG approaches The visualization of results is accomplished by external packages such as ParaView MayaVi or Viper which is also part of the FENICS project 2 Examples The following examples conform to FFC version 0 9 4 ffc v and the DOLFIN core library version 0 9 9 which are installed on the compute servers of the faculty Note that the interfaces to the DOLFIN library are still changing sli
11. m specific classes that may be instantiated by the user in her main program and passed to the DOLFIN library These provide the functionality needed in the innermost loop of FE assembly rou tines e g e evaluation of basis functions and their derivatives e application of local degrees of freedom to a function Date April 1 2011 2 ROLAND HERZOG e local to global mapping of the degrees of freedom e evaluation of the local element tensor e g the local stiffness matrix The code which implements this functionality is automatically generated by the form compiler from the high level description in the uf1 file By contrast from within a Python program the form compiler is called automatically whenever nec essary just in time Applications FEniCS Apps Application Interfaces SEER DOLFIN Core components Instant External libraries PETSe UBLAS UMFRACK NumPy Scorch FIGURE 1 1 Interaction of components in FENICS UFL Uni fied Form Language FFC Unified Form Compiler UFC Uni fied Form Assembly Code Highlighted are the components we refer to this text 1 2 Some Features and Limitations FENICS currently supports the follow ing types of finite element function spaces see Section 4 3 e H conforming FE spaces composed of Lagrange elements of arbitrary degree e H div conforming FE spaces composed of Raviart Thomas elements of arbitrary degree Brezzi Douglas Mar
12. omains function only vertex and cell valued meshfunctions can be plotted in viper boundary_subdomains MeshFunction uint boundary_mesh 1 boundary_subdomains values subdomains values cell_map values Plot the mesh and boundary mesh mesh_plot plot mesh title Mesh boundary_mesh_plot plot boundary_mesh title Boundary Mesh boundary_subdomains_plot plot boundary_subdomains title Boundary Subdomains Save the boundary subdomains plot as a file boundary_subdomains_plot write_png stokes_subdomains png uD Save the boundary mesh this includes the cell map and vertex map mesh data mesh_file File stokes_boundary_mesh xml mesh_file lt lt boundary_mesh 2 3 A Nonlinear Equation Nonlinear equations can be solved through the use of Newton s method We denote the nonlinear residual by F u Every Newton step requires the solution of J u du F u where J u F u is the Jacobian of F In the context of PDEs F u is a vector whose components are generated from a linear form L u v as v ranges over all test functions F w L vi u And thus the Jacobian matrix corresponds to a bilinear form a in the following sense J u du a vi du u Recall that test functions are always the first argument in a bilinear form and trial functions are second It is a particular feature of FENICS that the Jacobian can be generated automatically from
13. rder of declaration not in the order in which they appear in the form This means that in a bilinear form one should always define the test function first v BasisFunction Lagrange triangle 2 u BasisFunction Lagrange triangle 2 a inner grad v grad u However TestFunction and TrialFunction are special instances of BasisFunction with the property that a TestFunction will always be the first argument in a form and a TrialFunction will always be the second argument Coefficients are declared by the keyword Function or the special cases Constant VectorConstant or TensorConstant in UFL In Python they can also be generated by the Expression constructor In UFL one can make extensive use of index notation and tensor algebraic concepts but we don t go into details here We only mention the following operators which act on tensors of appropriate ranks transpose tr trace dot collapse the last axis of the left and the first axis of the right argument inner collapse all axes outer outer product cross defined only for two vectors of length det determinant dev deviatoric part of a matrix sym and skew symmetric and skey symmetric parts of a matrix cofac cofactor of a matrix inv inverse of a matrix Forms defining PDEs need spatial derivatives of basis functions One may use f f Dx v 2 or equivalently v dx 2 to denote 0v 0x2 for example Useful abbreviations e
14. s are currently not supported assemble a assemble L mfile File A_before_bc m mfile lt lt A mfile File b_before_bc m mfile lt lt b H HHH OO OH HHH i Apply the Dirichlet boundary conditions and export the linear system again bc apply A b mfile File A_after_bc m mfile lt lt A mfile File b_after_bc m mfile lt lt b Solve the manually assembled problem U2 Function w solve A U2 vector b Compute the difference should be zero of velocities u2 p2 U2 split True diff u vector u2 vector print Norm of velocity difference 4 4e diff norm 12 Go interactive interactive Note that we have not enforced the condition Jo pdx 0 and thus the system matrix has a kernel of dimension one This is handled by the linear algebra sub system Changing the Boundary Conditions We now consider a different possibility of specifying Dirichlet boundary conditions Instead of using boolean indicator func tions like DirichletBoundary above one may use a MeshFunction to distinguish various parts of the boundary In general a MeshFunction is a function which is defined on the entities of a given mesh e g on its cells facets or vertices It can be used for many purposes such as to represent subdomains of different material to mark cells for refinement or to identify different types of boundary conditions 10 ROLAND HERZOG 0 00 0 5
15. the same with lower level code We manually create the problem class SemilinearEquation derived from the NonlinearProblem class Note that we need to implement the evaluation routine of the residual F v u and of the Jacobian J v du u These consist in calling the assemble routine and then applying the Dirichlet boundary conditions for the Newton step We also need to manually create the Newton solver instance We now solve the same problem again but we manually create the NonlinearProblem instance and the solver See also demo pde cahn hilliard python demo py class SemilinearEquation NonlinearProblem 14 ROLAND HERZOG The constructor is defined to take all arguments we need der Minit Csel ma Ih pe NonlinearProblem __init__ self Saala JL Ib self a a self bc bec self reset_sparsity True evaluation of residual F x def F self b x assemble self L tensor b self bc apply b evaluation of Jacobian F x def J self A x assemble self a tensor A reset_sparsity self reset_sparsity self bc apply A self reset sparsity False Define another solution function u2 Function V v TestFunction V du2 TrialFunction V This means that we also need to define the variational problem again L2 inner grad v grad u2 dx v u2 3 dx v f dx a2 derivative L2 u2 du2 Create nonlinear problem and Newton solver problem2 SemilinearEquation a2 L2 bc solver New
16. tonSolver lu solver parameters convergence_criterion incremental solver parameters relative_tolerance 1e 6 Solve the problem solver solve problem2 u2 vector Compare the solutions diff u vector u2 vector print Norm of solution difference 4e diff norm 12 up Suggestions Modify the problem to include non homogeneous boundary condi tions u g on and mixed boundary conditions 2 4 Further Reading e FENICS web site e FENICS tutorial e UFC specifications INTRODUCTION TO FENICS 15 3 Overview of Unified Form Language UFL According to the UFL Specification and User Manual UFL is designed to express variational forms of the following kind Ne a c A U1 Up W1 Wn J f Tp U1 Ur W1 Wn dz k 1 Qk Ne e 5 Tp v1 Ur W1 Wn ds k 1 Ok ni i 5 I Tj U1 5 Ur W1 Wn dS k 1 Ts where v U are form arguments basis functions and w1 Wp are form coefficients Moreover dz ds and dS refer to integration over a domain Q its boundary experior facets OQ x or its interior facets T respectively The number r of basis functions represents the arity of the form For instance bilinear forms have arity r 2 linear forms have arity r 1 Basis functions are declared by the keyword BasisFunction in the UFL Note When defining a multivariate form the arguments v v are assigned in their o
17. where uO Expression x 1 gt 1 0 str DOLFIN_EPS 0 be DirichletBC W sub 0 u0 DirichletBoundary Define variational problem v q TestFunctions W same as vq TestFunction W v q split vq u p TrialFunctions CW f Constant 0 0 0 0 mu Constant 1 0e 1 n FacetNormal mesh a muxinner grad v grad u div v p q div u dx p dot v n ds L inner v f dx Setup problem and compute solution problem VariationalProblem a L bc problem parameters linear_solver direct U problem solve Split the mixed solution using a deep copy copy data instead of sharing needed to compute norms below u p U split True The average value of the pressure is currently random since the condition int p dx 0 has not yet been enforced We adjust the average of p now as a post processing step int_p p dx average_p assemble int_p mesh mesh p_array p vector array average_p p vector p_array Save solution in VTK format ufile_pvd File stokes_velocity pvd ufile_pvd lt lt u pfile_pvd File stokes_pressure pvd pfile_pvd lt lt p Plot solution plot u title Velocity INTRODUCTION TO FENICS 9 plot sqrt u 0 2 u 1 2 title Absolute Velocity plot p title Pressure Alternatively we can assemble the linear system manwanmiy pn hiie horens BING POOR Sallda Isis Cre exportan e ale M File
18. xist for often used compound derivatives such as grad div and curl Jump and average operators for interior facets are also defined which are needed to solve problems with discontinuous Galerkin DG function spaces UFL also supports differentiation with respect to user defined variables see the demo file HyperElasticity ufl for an example and differentiation with respect to basis functions The latter can be used for instance to automatically generate linearized equations e g for Newton iterations 16 ROLAND HERZOG A Setup Environment Currently FENICS is installed on the compute servers in the following versions e FFC Version 0 9 4 e DOLFIN Version 0 9 9 Here is how to setup the environment Set environment variables source usr local FEniCS 0 9 9 setup_fenics References CHEMNITZ UNIVERSITY OF TECHNOLOGY FACULTY OF MATHEMATICS D 09107 CHEMNITZ GER MANY E mail address roland herzog mathematik tu chemnitz de URL http www tu chemnitz de herzog
Download Pdf Manuals
Related Search