Fem-fenics

Package for solving Partial Differential Equations based on Fenics.

Introduction

Fem-Fenics Octave is a package for solving partial differential equations. Obviously, Fem-fenics is not the only extra package for Octave with this purpose. For example, Bim_package uses finite volumes to solve diffusion-advection-reaction equations, while secs1d/2d/3d [1] are suited for the resolution of the drift-diffusion system. Furthermore, to use profitably the software, you can integrate it with msh [2] for the generation of the mesh and with fpl [3] for the post-processing of data. The objective of Fem-fenics is to be a generic library of finite elements, thereby allowing the user to resolve any type of pde, choosing also the most appropriate Finite Element space for any specific problem.

As the name suggests, the Fem-fenics pkg is a wrapper for FEniCS [4] functions and classes. Thus, ideally the Fem-fenics final goal would be to be able to reproduce all the features available in FEniCS, simplifying them where it is possible or using the Octave function whenever available (like the "\" for the resolution of a linear system or the odepkg [5] for the resolution of a time dependent problem).

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}$

A complete description of the problem is avilable on the Fenics website.

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) © Copyright 2011, The FEniCS Project

Mixed Formulation for Poisson Equation

A complete description of the problem is avilable on the Fenics website.

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() © Copyright 2011, The FEniCS Project

Hyperelasticity

This time we compare the code with the c++ version of DOLFIN. A complete description of the problem is avilable on the Fenics website. The final solution will look like this

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) © Copyright 2011, The FEniCS Project

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; } © Copyright 2011, The FEniCS Project

Incompressible Navier-Stokes equations

A complete description of the problem is avilable on the Fenics website.

Code: Define HyperElasticity problem with fem-fenics

pkg load fem-fenics msh import_ufl_Problem ("TentativeVelocity"); import_ufl_Problem ("VelocityUpdate"); import_ufl_Problem ("PressureUpdate"); # We can either load the mesh from the file as in Dolfin but # we can also use the msh pkg to generate the L-shape domain name = [tmpnam ".geo"]; fid = fopen (name, "w"); fputs (fid,"Point (1) = {0, 0, 0, 0.1};\n"); fputs (fid,"Point (2) = {1, 0, 0, 0.1};\n"); fputs (fid,"Point (3) = {1, 0.5, 0, 0.1};\n"); fputs (fid,"Point (4) = {0.5, 0.5, 0, 0.1};\n"); fputs (fid,"Point (5) = {0.5, 1, 0, 0.1};\n"); fputs (fid,"Point (6) = {0, 1, 0,0.1};\n"); fputs (fid,"Line (1) = {5, 6};\n"); fputs (fid,"Line (2) = {2, 3};\n"); fputs (fid,"Line(3) = {6,1,2};\n"); fputs (fid,"Line(4) = {5,4,3};\n"); fputs (fid,"Line Loop(7) = {3,2,-4,1};\n"); fputs (fid,"Plane Surface(8) = {7};\n"); fclose (fid); msho = msh2m_gmsh (canonicalize_file_name (name)(1:end-4),... "scale", 1,"clscale", .2); unlink (canonicalize_file_name (name)); mesh = Mesh (msho); # Define function spaces (P2-P1). From ufl file # V = VectorElement("CG", triangle, 2) # Q = FiniteElement("CG", triangle, 1) V = FunctionSpace ('VelocityUpdate', mesh); Q = FunctionSpace ('PressureUpdate', mesh); # Define trial and test functions. From ufl file # u = TrialFunction(V) # p = TrialFunction(Q) # v = TestFunction(V) # q = TestFunction(Q) # Set parameter values. From ufl file # nu = 0.01 dt = 0.01; T = 3.; # Define boundary conditions noslip = DirichletBC (V, @(x,y) [0; 0], [3, 4]); outflow = DirichletBC (Q, @(x,y) 0, 2); # Create functions u0 = Expression ('u0', @(x,y) [0; 0]); # Define coefficients k = Constant ('k', dt); f = Constant ('f', [0; 0]); # Tentative velocity step. From ufl file # eq = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx \ # + nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx a1 = BilinearForm ('TentativeVelocity', V, k); # Pressure update. From ufl file # a = inner(grad(p), grad(q))*dx # L = -(1/k)*div(u1)*q*dx a2 = BilinearForm ('PressureUpdate', Q); # Velocity update # a = inner(u, v)*dx # L = inner(u1, v)*dx - k*inner(grad(p1), v)*dx a3 = BilinearForm ('VelocityUpdate', V); # Assemble matrices A1 = assemble (a1, noslip); A3 = assemble (a3, noslip); # Time-stepping t = dt; i = 0; while t < T # Update pressure boundary condition inflow = DirichletBC (Q, @(x,y) sin(3.0*t), 1); # Compute tentative velocity step "Computing tentative velocity" L1 = LinearForm ('TentativeVelocity', V, k, u0, f); b1 = assemble (L1, noslip); utmp = A1 \ b1; u1 = Function ('u1', V, utmp); # Pressure correction "Computing pressure correction" L2 = LinearForm ('PressureUpdate', Q, u1, k); [A2, b2] = assemble_system (a2, L2, inflow, outflow); ptmp = A2 \ b2; p1 = Function ('p1', Q, ptmp); # Velocity correction "Computing velocity correction" L3 = LinearForm ('VelocityUpdate', V, k, u1, p1); b3 = assemble (L3, noslip); ut = A3 \ b3; u1 = Function ('u0', V, ut); # Plot solution plot (p1); plot (u1); # Save to file save (p1, sprintf ("p_%3.3d", ++i)); save (u1, sprintf ("u_%3.3d", i)); # Move to next time step u0 = u1; t += dt end

Code: Define NS problem with fenics python

from dolfin import * # Load mesh from file mesh = Mesh("lshape.xml.gz") # Define function spaces (P2-P1) V = VectorFunctionSpace(mesh, "CG", 2) Q = FunctionSpace(mesh, "CG", 1) # Define trial and test functions u = TrialFunction(V) p = TrialFunction(Q) v = TestFunction(V) q = TestFunction(Q) # Set parameter values dt = 0.01 T = 3 nu = 0.01 # Define time-dependent pressure boundary condition p_in = Expression("sin(3.0*t)", t=0.0) # Define boundary conditions noslip = DirichletBC(V, (0, 0), "on_boundary && \ (x[0] < DOLFIN_EPS | x[1] < DOLFIN_EPS | \ (x[0] > 0.5 - DOLFIN_EPS && x[1] > 0.5 - DOLFIN_EPS))") inflow = DirichletBC(Q, p_in, "x[1] > 1.0 - DOLFIN_EPS") outflow = DirichletBC(Q, 0, "x[0] > 1.0 - DOLFIN_EPS") bcu = [noslip] bcp = [inflow, outflow] # Create functions u0 = Function(V) u1 = Function(V) p1 = Function(Q) # Define coefficients k = Constant(dt) f = Constant((0, 0)) # Tentative velocity step F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx \ + nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx a1 = lhs(F1) L1 = rhs(F1) # Pressure update a2 = inner(grad(p), grad(q))*dx L2 = -(1/k)*div(u1)*q*dx # Velocity update a3 = inner(u, v)*dx L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx # Assemble matrices A1 = assemble(a1) A2 = assemble(a2) A3 = assemble(a3) # Use amg preconditioner if available prec = "amg" if has_krylov_solver_preconditioner("amg") else "default" # Create files for storing solution ufile = File("results/velocity.pvd") pfile = File("results/pressure.pvd") # Time-stepping t = dt while t < T + DOLFIN_EPS: # Update pressure boundary condition p_in.t = t # Compute tentative velocity step begin("Computing tentative velocity") b1 = assemble(L1) [bc.apply(A1, b1) for bc in bcu] solve(A1, u1.vector(), b1, "gmres", "default") end() # Pressure correction begin("Computing pressure correction") b2 = assemble(L2) [bc.apply(A2, b2) for bc in bcp] solve(A2, p1.vector(), b2, "gmres", prec) end() # Velocity correction begin("Computing velocity correction") b3 = assemble(L3) [bc.apply(A3, b3) for bc in bcu] solve(A3, u1.vector(), b3, "gmres", "default") end() # Plot solution plot(p1, title="Pressure", rescale=True) plot(u1, title="Velocity", rescale=True) # Save to file ufile << u1 pfile << p1 # Move to next time step u0.assign(u1) t += dt print "t =", t # Hold plot interactive() © Copyright 2011, The FEniCS Project