# Fem-fenics

Package for solving Partial Differential Equations based on Fenics.

## Tutorials

### Poisson Equation

Here is a first example for the solution of the Poisson equation. The equation being solved is

 Code: Define Poisson problem with fem-fenics  pkg load fem-fenics msh import_ufl_Problem ('Poisson') # Create mesh and define function space x = y = linspace (0, 1, 33); mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4)); # File Poisson.ufl # element = FiniteElement("Lagrange", triangle, 1) V = FunctionSpace('Poisson', mesh); # Define boundary condition bc = DirichletBC(V, @(x, y) 0.0, [2;4]); # Define variational problem # File Poisson.ufl # u = TrialFunction(element) # v = TestFunction(element) # f = Coefficient(element) # g = Coefficient(element) f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02)); g = Expression ('g', @(x,y) sin (5.0 * x)); # File Poisson.ufl # a = inner(grad(u), grad(v))*dx # L = f*v*dx + g*v*ds a = BilinearForm ('Poisson', V); L = LinearForm ('Poisson', V, f, g); # Compute solution [A, b] = assemble_system (a, L, bc); sol = A \ b; u = Function ('u', V, sol); # Save solution in VTK format save (u, 'poisson') # Plot solution plot (u); 
 Code: Define Poisson problem with fenics python  from dolfin import * # Create mesh and define function space mesh = UnitSquareMesh(32, 32) V = FunctionSpace(mesh, "Lagrange", 1) # Define Dirichlet boundary (x = 0 or x = 1) def boundary(x): return x < DOLFIN_EPS or x > 1.0 - DOLFIN_EPS # Define boundary condition u0 = Constant(0.0) bc = DirichletBC(V, u0, boundary) # Define variational problem u = TrialFunction(V) v = TestFunction(V) f = Expression("10*exp(-(pow(x - 0.5, 2) + pow(x - 0.5, 2)) / 0.02)") g = Expression("sin(5*x)") a = inner(grad(u), grad(v))*dx L = f*v*dx + g*v*ds # Compute solution u = Function(V) (A, b) = assemble_system (a, L, bc); solve(A, u.vector(), b, "gmres", "default") # Save solution in VTK format file = File("poisson.pvd") file << u # Plot solution plot(u, interactive=True) 

### Mixed Formulation for Poisson Equation

 Code: Define MixedPoisson problem with fem-fenics  pkg load fem-fenics msh import_ufl_Problem ('MixedPoisson') # Create mesh x = y = linspace (0, 1, 33); mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4)); # File MixedPoisson.ufl # BDM = FiniteElement("BDM", triangle, 1) # DG = FiniteElement("DG", triangle, 0) # W = BDM * DG V = FunctionSpace('MixedPoisson', mesh); # Define trial and test function # File MixedPoisson.ufl # (sigma, u) = TrialFunctions(W) # (tau, v) = TestFunctions(W) # CG = FiniteElement("CG", triangle, 1) # f = Coefficient(CG) f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02)); # Define variational form # File MixedPoisson.ufl # a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx # L = - f*v*dx a = BilinearForm ('MixedPoisson', V); L = LinearForm ('MixedPoisson', V, f); # Define essential boundary bc1 = DirichletBC (SubSpace (V, 1), @(x,y) [0; -sin(5.0*x)], 1); bc2 = DirichletBC (SubSpace (V, 1), @(x,y) [0; sin(5.0*x)], 3); # Compute solution [A, b] = assemble_system (a, L, bc1, bc2); sol = A \ b; func = Function ('func', V, sol); sigma = Function ('sigma', func, 1); u = Function ('u', func, 2); # Plot solution plot (sigma); plot (u); 
 Code: Define MixedPoisson problem with fenics python  from dolfin import * # Create mesh mesh = UnitSquareMesh(32, 32) # Define function spaces and mixed (product) space BDM = FunctionSpace(mesh, "BDM", 1) DG = FunctionSpace(mesh, "DG", 0) W = BDM * DG # Define trial and test functions (sigma, u) = TrialFunctions(W) (tau, v) = TestFunctions(W) f = Expression ("10*exp(-(pow(x - 0.5, 2) + pow(x - 0.5, 2)) / 0.02)") # Define variational form a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx L = - f*v*dx # Define function G such that G \cdot n = g class BoundarySource(Expression): def __init__(self, mesh): self.mesh = mesh def eval_cell(self, values, x, ufc_cell): cell = Cell(self.mesh, ufc_cell.index) n = cell.normal(ufc_cell.local_facet) g = sin(5*x) values = g*n values = g*n def value_shape(self): return (2,) G = BoundarySource(mesh) # Define essential boundary def boundary(x): return x < DOLFIN_EPS or x > 1.0 - DOLFIN_EPS bc = DirichletBC(W.sub(0), G, boundary) # Compute solution w = Function(W) solve(a == L, w, bc) (sigma, u) = w.split() # Plot sigma and u plot(sigma) plot(u) interactive() 

### Mixed Formulation for Poisson Equation

This time we compare the code with the c++ version of DOLFIN.

 Code: HyperElasticity Problem: the ufl file  # Function spaces element = VectorElement("Lagrange", tetrahedron, 1) # Trial and test functions du = TrialFunction(element) # Incremental displacement v = TestFunction(element) # Test function # Functions u = Coefficient(element) # Displacement from previous iteration B = Coefficient(element) # Body force per unit volume T = Coefficient(element) # Traction force on the boundary # Kinematics I = Identity(element.cell().d) # Identity tensor F = I + grad(u) # Deformation gradient C = F.T*F # Right Cauchy-Green tensor # Invariants of deformation tensors Ic = tr(C) J = det(F) # Elasticity parameters mu = Constant(tetrahedron) lmbda = Constant(tetrahedron) # Stored strain energy density (compressible neo-Hookean model) psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2 # Total potential energy Pi = psi*dx - inner(B, u)*dx - inner(T, u)*ds # First variation of Pi (directional derivative about u in the direction of v) F = derivative(Pi, u, v) # Compute Jacobian of F J = derivative(F, u, du) 
 Code: Define HyperElasticity problem with fem-fenics
 Code: Define HyperElasticity problem with fenics c++