Practical Lab WS 23/24 Practical Lab Numerical Simulation

Under direction of

Dr. Martin Lenz

This lab is not targeted at students from numerical simulation, but at students from analysis and other related fields, who are interested in an approach to develop finite element models in a syntax that is very general and close to the variational formulation of a PDE.

In the first part of the lab, we will get to know the toolbox FEniCS(x) based on introductory examples describing some of the most important features. In the second part, participants are invited to develop finite element models for PDEs of their own personal interest, and we will discuss ideas, problems and possible solutions.

To get an idea of how the toolbox works, the following code demonstrates how to solve linearized elasticity problems by just defining the energy functional E and letting FEniCS do all the work…

Time: by arrangement

Please register via eCampus

from fenics import *
from ufl.algorithms import replace

# Elasticity parameters
la = 1
mu = 1
g = Constant ((0, -1))

# Mesh and finite element space
mesh = RectangleMesh (Point (0, 0), Point (1, 1), 100, 100)

U = VectorFunctionSpace (mesh, 'CG', 1)

u = Function (U, name='displacement')
phi = TestFunction (U)
ut = TrialFunction (U)

# Boundary conditions
fixed = CompiledSubDomain ("on_boundary && near (x[1], 0)") # Dirichlet
force = CompiledSubDomain ("on_boundary && near (x[1], 1)") # Neumann

bc = DirichletBC (U, Constant ((0, 0)), fixed)

# Introduce a boundary measure ds(1) on the Neumann boundary
subdomains = MeshFunction ("size_t", mesh, mesh.topology().dim() - 1)
force.mark (subdomains, 1)
ds = Measure ('ds', domain=mesh, subdomain_data=subdomains)

# Energy
def epsilon (u):
    return 0.5 * (grad(u) + grad(u).T)
def sigma (u):
    return la * div (u) * Identity (2) + 2 * mu * epsilon (u)
def edensity (u):
    return 0.5 * inner (sigma (u), epsilon (u))

E = edensity (u) * dx - dot (g, u) * ds(1)

# Automatic differentiation to compute the necessary condition
duE = replace (derivative (E, u, phi), {u: ut})

# Solve the linear system of equations
solve (lhs (duE) == rhs (duE), u, bc)