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)
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.
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 = 4*pi*1e-7 # vacuum elif markers[ufc_cell.index] == 1: values = 1e-5 # iron (should really be 2.5e-1) else: values = -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.
Dante June 8, 2016
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.
logg June 8, 2016 — Post Author
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.