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

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 = Poisson_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) # Plot solution plot (u);

Code: Define Poisson problem

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)