Mesh generation in FEniCS

FEniCS comes with built-in mesh generation that allows relatively complex domains to be defined and meshed using simple Python code. The meshing capabilities of FEniCS are handled by the (optional) FEniCS component mshr and the actual meshing is performed by CGAL or Tetgen. The interface is inspired by that of Netgen in that domains are defined in terms of Constructive Solid Geometry (CSG).

For example, here’s how to create a mesh of the Death Star in FEniCS:

from fenics import *
from mshr import *
from math import pi, sin, cos, sqrt

# Parameters
R = 1.1
r = 0.4
t = 10
x = R*cos(float(t) / 180 * pi)
y = 0
z = R*sin(t)

# Create geometry
s1 = Sphere(Point(0, 0, 0), 1)
s2 = Sphere(Point(x, y, z), r)
b1 = Box(Point(-2, -2, -0.03), Point(2, 2, 0.03))
geometry = s1 - s2 - b1

# Create mesh
mesh = generate_mesh(geometry, 32)

Death Star

mshr is the work by my student Benjamin Kehlet who has done a great job with designing an intuitive interface and handling corner cases, cleanups of geometries and various other issues that must be handled to deliver a smooth user experience. A feature that he recently added is the ability to define subdomains as part of domains and meshes. The following code taken from the upcoming FEniCS tutorial illustrates how to define subdomains in a mesh, here for a 2D cross-section of an iron cylinder with a coppper wire wound 10 times around the cylinder:

from fenics import *
from mshr import *
from math import sin, cos, pi

a = 1.0   # inner radius of iron cylinder
b = 1.2   # outer radius of iron cylinder
c_1 = 0.8 # radius for inner circle of copper wires
c_2 = 1.4 # radius for outer circle of copper wires
r = 0.1   # radius of copper wires
R = 5.0   # radius of domain
n = 10    # number of windings

# FIXME: Use 'domain' instead of 'geometry' in other examples

# Define geometry for background
domain = Circle(Point(0, 0), R)

# Define geometry for iron cylinder
cylinder = Circle(Point(0, 0), b) - Circle(Point(0, 0), a)

# Define geometry for wires (N = North (up), S = South (down))
angles_N = [i*2*pi/n for i in range(n)]
angles_S = [(i + 0.5)*2*pi/n for i in range(n)]
wires_N = [Circle(Point(c_1*cos(v), c_1*sin(v)), r) for v in angles_N]
wires_S = [Circle(Point(c_2*cos(v), c_2*sin(v)), r) for v in angles_S]

# Set subdomain for iron cylinder
domain.set_subdomain(1, cylinder)

# Set subdomains for wires
for (i, wire) in enumerate(wires_N):
    domain.set_subdomain(2 + i, wire)
for (i, wire) in enumerate(wires_S):
    domain.set_subdomain(2 + n + i, wire)

# Create mesh
mesh = generate_mesh(domain, 32)

The following plot shows a detail of the mesh. The circles representing cross-sections of the copper wire are clearly visible, as is the outline of the iron cylinder.

Magnetostatics mesh

The subdomain information is carried over to the FEniCS mesh and can be used to define variational problems, set boundary conditions etc. Here’s a code snippet illustrating how to extract the subdomain data and use it for defining measures for the different subdomains:

# Define subdomain markers and integration measure
markers = MeshFunction('size_t', mesh, 2, mesh.domains())
dx = Measure('dx', domain=mesh, subdomain_data=markers)

# Define current densities
J_N = Constant(1.0)
J_S = Constant(-1.0)

# Define magnetic permeability
class Permeability(Expression):
    def __init__(self, mesh, **kwargs):
        self.markers = markers
    def eval_cell(self, values, x, ufc_cell):
        if markers[ufc_cell.index] == 0:
            values[0] = 4*pi*1e-7 # vacuum
        elif markers[ufc_cell.index] == 1:
            values[0] = 1e-5      # iron (should really be 2.5e-1)
        else:
            values[0] = -6.4e-6   # copper

mu = Permeability(mesh, degree=1)

# Define variational problem
A_z = TrialFunction(V)
v = TestFunction(V)
a = (1 / mu)*dot(grad(A_z), grad(v))*dx
L_N = sum(J_N*v*dx(i) for i in range(2, 2 + n))
L_S = sum(J_S*v*dx(i) for i in range(2 + n, 2 + 2*n))
L = L_N + L_S

Finally, calling `solve(a == L)` in FEniCS gives us the solution. The plots below show the (z-component of the) magnetic vector potential and the magnetic field.

Magnetostatics potential

Magnetostatics magnetic field

Recent Posts

Recent Comments

Meta

logg Written by:

Anders Logg is Professor of Computational Mathematics at Chalmers University of Technology. His research interests are adaptive finite element methods, high-level automating software systems for solution of PDE, domain-specific languages and compilers in scientific computing, and applications in biomedicine, general relativity, and architecture. Logg is head of the unit for Computational Mathematics and excellence profile leader within Chalmers Area of Advance the Built Environment. He is director of the Swedish Network for Mathematics in Industry and member of the EMS Applied Mathematics Committee and the Swedish National Committee for Mathematics. He is co-founder and a core developer of the FEniCS Project, a leading open-source software for automated solution of PDE. Logg works part-time as scientific advisor to Fraunhofer-Chalmers Centre and Simula Research Laboratory.

2 Comments

  1. Dante
    June 8, 2016
    Reply

    It looks great.
    I’ve already try generate some 2D geometry and I have one question:
    How I can generate a structured mesh?
    For example, I want create 2D rectangular mesh with rectangular holes.

    • June 8, 2016
      Reply

      There’s no problem doing rectangular holes but I don’t think it’s possible (either because it’s not available in the meshing backends or we don’t have an interface for it) to guarantee that the mesh is structured.

Leave a Reply

Your email address will not be published. Required fields are marked *