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

$-\mathrm {div} \ (\nabla u(x,y)))=1\qquad {\mbox{ in }}\Omega$

$u(x,y)=0\qquad {\mbox{ on }}\Gamma _{D}$

$(\nabla u(x,y))\cdot \mathbf {n} =0\qquad {\mbox{ on }}\Gamma _{N}$

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[0] < DOLFIN_EPS or x[0] > 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] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)") g = Expression("sin(5*x[0])") 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] - 0.5, 2) + pow(x[1] - 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[0]) values[0] = g*n[0] values[1] = g*n[1] def value_shape(self): return (2,) G = BoundarySource(mesh) # Define essential boundary def boundary(x): return x[1] < DOLFIN_EPS or x[1] > 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()

Hyperelasticity

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

pkg load fem-fenics msh problem = 'HyperElasticity'; import_ufl_Problem (problem); Rotation = @(x,y,z) ... [0; ... 0.5*(0.5 + (y - 0.5)*cos(pi/3) - (z-0.5)*sin(pi/3) - y);... 0.5*(0.5 + (y - 0.5)*sin(pi/3) + (z-0.5)*cos(pi/3) - z)]; #Create mesh and define function space x = y = z = linspace (0, 1, 17); mshd = Mesh (msh3m_structured_mesh (x, y, z, 1, 1:6)); V = FunctionSpace (problem, mshd); # Create Dirichlet boundary conditions bcl = DirichletBC (V, @(x,y,z) [0; 0; 0], 1); bcr = DirichletBC (V, Rotation, 2); bcs = {bcl, bcr}; # Define source and boundary traction functions B = Constant ('B', [0.0; -0.5; 0.0]); T = Constant ('T', [0.1; 0.0; 0.0]); # Set material parameters E = 10.0; nu = 0.3; mu = Constant ('mu', E./(2*(1+nu))); lmbda = Constant ('lmbda', E*nu./((1+nu)*(1-2*nu))); u = Expression ('u', @(x,y,z) [0; 0; 0]); # Create (linear) form defining (nonlinear) variational problem L = ResidualForm (problem, V, mu, lmbda, B, T, u); # Solve nonlinear variational problem F(u; v) = 0 u0 = assemble (L, bcs{:}); # Create function for the resolution of the NL problem # function [y, jac] = f (problem, xx, V, bc1, bc2, B, T, mu, lmbda) # u = Function ('u', V, xx); # a = JacobianForm (problem, V, mu, lmbda, u); # L = ResidualForm (problem, V, mu, lmbda, B, T, u); # if (nargout == 1) # [y, xx] = assemble (L, xx, bc1, bc2); # elseif (nargout == 2) # [jac, y, xx] = assemble_system (a, L, xx, bc1, bc2); # endif # endfunction fs = @(xx) f (problem, xx, V, bcl, bcr, B, T, mu, lmbda); [x, fval, info] = fsolve (fs, u0, optimset ("jacobian", "on")); func = Function ('u', V, x); # Save solution in VTK format save (func, 'displacement'); # Plot solution plot (func);

Code: Define HyperElasticity problem with fenics c++

#include <dolfin.h> #include "HyperElasticity.h" using namespace dolfin; // Sub domain for clamp at left end class Left : public SubDomain { bool inside(const Array<double>& x, bool on_boundary) const { return (std::abs(x[0]) < DOLFIN_EPS) && on_boundary; } }; // Sub domain for rotation at right end class Right : public SubDomain { bool inside(const Array<double>& x, bool on_boundary) const { return (std::abs(x[0] - 1.0) < DOLFIN_EPS) && on_boundary; } }; // Dirichlet boundary condition for clamp at left end class Clamp : public Expression { public: Clamp() : Expression(3) {} void eval(Array<double>& values, const Array<double>& x) const { values[0] = 0.0; values[1] = 0.0; values[2] = 0.0; } }; // Dirichlet boundary condition for rotation at right end class Rotation : public Expression { public: Rotation() : Expression(3) {} void eval(Array<double>& values, const Array<double>& x) const { const double scale = 0.5; // Center of rotation const double y0 = 0.5; const double z0 = 0.5; // Large angle of rotation (60 degrees) double theta = 1.04719755; // New coordinates double y = y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta); double z = z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta); // Rotate at right end values[0] = 0.0; values[1] = scale*(y - x[1]); values[2] = scale*(z - x[2]); } }; int main() { // Create mesh and define function space UnitCubeMesh mesh (16, 16, 16); HyperElasticity::FunctionSpace V(mesh); // Define Dirichlet boundaries Left left; Right right; // Define Dirichlet boundary functions Clamp c; Rotation r; // Create Dirichlet boundary conditions DirichletBC bcl(V, c, left); DirichletBC bcr(V, r, right); std::vector<const BoundaryCondition*> bcs; bcs.push_back(&bcl); bcs.push_back(&bcr); // Define source and boundary traction functions Constant B(0.0, -0.5, 0.0); Constant T(0.1, 0.0, 0.0); // Define solution function Function u(V); // Set material parameters const double E = 10.0; const double nu = 0.3; Constant mu(E/(2*(1 + nu))); Constant lambda(E*nu/((1 + nu)*(1 - 2*nu))); // Create (linear) form defining (nonlinear) variational problem HyperElasticity::ResidualForm F(V); F.mu = mu; F.lmbda = lambda; F.B = B; F.T = T; F.u = u; // Create jacobian dF = F' (for use in nonlinear solver). HyperElasticity::JacobianForm J(V, V); J.mu = mu; J.lmbda = lambda; J.u = u; // Solve nonlinear variational problem F(u; v) = 0 solve(F == 0, u, bcs, J); // Save solution in VTK format File file("displacement.pvd"); file << u; // Plot solution plot(u); interactive(); return 0; }