Mesh, Submesh , Torus, Boundary

mesh
3

Hi @Noya,

You may want to try creating your mesh using an external meshing tool, e.g. Gmsh. The advantage of this approach is that you can mark the subdomains and boundaries using Gmsh, save the markers to files, and load them into FEniCS MeshFunctions, which can be used to define integration measures and boundary conditions. Using your current code (i.e. with BoxMesh), you obtain a mesh that doesn’t respect the internal boundary of the torus, i.e. you will have elements that pass through the torus boundary.

Below is a code using the Gmsh Python API to generate the mesh, and meshio to convert to XDMF files. Running this code produces three files: mesh.xdmf contains only the mesh, domains.xdmf contains markers which identify the domains \Omega and T, and boundaries.xdmf contains markers that identify the boundaries \partial \Omega and \partial T.

import gmsh
import meshio
import numpy as np

gmsh.initialize()

R = 0.5  # Major radius
r = 0.1  # Minor radius
L = 2    # Box dimension

hmin = r/2   # Min mesh size
hmax = L/10  # Max mesh size

# Create geometry
gmsh.model.occ.addTorus(0,0,0,R,r,tag=1,angle=np.pi)
gmsh.model.occ.addTorus(0,0,0,R,r,tag=2,angle=np.pi)
gmsh.model.occ.rotate([(3,2)],0,0,0,0,0,1,np.pi)
gmsh.model.occ.addBox(-L/2,-L/2,-L/2,L,L,L,tag=3)
gmsh.model.occ.fragment([(3,3)], [(3,1),(3,2)])
gmsh.model.occ.synchronize()

# Create mesh fields to control element size
d = gmsh.model.mesh.field.add("MathEval")
gmsh.model.mesh.field.setString(d, "F", "Sqrt((Sqrt(x^2+y^2)-{R:f})^2+z^2)".format(R=R))
t = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(t, "InField", d)
gmsh.model.mesh.field.setNumber(t, "DistMin", r)
gmsh.model.mesh.field.setNumber(t, "DistMax", 2*r)
gmsh.model.mesh.field.setNumber(t, "SizeMin", hmin)
gmsh.model.mesh.field.setNumber(t, "SizeMax", hmax)
gmsh.model.mesh.field.setAsBackgroundMesh(t)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)

# Generate physical groups to mark domains and boundaries
gmsh.model.addPhysicalGroup(3, [1,2], 1)         # Torus
gmsh.model.addPhysicalGroup(3, [3], 2)           # Surroundings
gmsh.model.addPhysicalGroup(2, [7,8], 3)         # Internal boundary
gmsh.model.addPhysicalGroup(2,[1,2,3,4,5,6], 4)  # External boundary

# Generate mesh and save
gmsh.model.mesh.generate()
gmsh.write("torus_mesh.msh")
gmsh.finalize()

# Convert mesh to XDMF for use in FEniCS
msh = meshio.read("torus_mesh.msh")
tet_data = msh.cell_data_dict["gmsh:physical"]["tetra"]
tri_data = msh.cell_data_dict["gmsh:physical"]["triangle"]

meshio.write("mesh.xdmf",
    meshio.Mesh(points=msh.points,
        cells={"tetra": msh.cells_dict["tetra"]}
    )
)
meshio.write("boundaries.xdmf",
    meshio.Mesh(points=msh.points,
        cells={"triangle": msh.cells_dict["triangle"]},
        cell_data={"bnd_marker": [msh.cell_data_dict["gmsh:physical"]["triangle"]]}
    )
)
meshio.write("domains.xdmf",
     meshio.Mesh(
        points=msh.points,
        cells={"tetra": msh.cells_dict["tetra"]},
        cell_data={"dom_marker": [msh.get_cell_data("gmsh:physical","tetra")]}
    )
)

You can then load the markers contained in domains.xdmf and boundaries.xdmf, e.g.:

from fenics import*
import numpy as np

# -------------------------------------------------
#Mesh Omega
N = 50
mesh = Mesh()
with XDMFFile("mesh.xdmf") as mshfile:
    mshfile.read(mesh)

# -------------------------------------------------
# Meshfunction over domains
mft = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
with XDMFFile("domains.xdmf") as domfile:
    domfile.read(mft)
File("mft.pvd") << mft
dxT = Measure('dx', domain=mesh, subdomain_data=mft)

submesh_torus = SubMesh(mesh, mft, 1)
submesh_OT = SubMesh(mesh, mft, 2)

# -------------------------------------------------

V = FiniteElement("N1curl", mesh.ufl_cell(), 1)
S = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Z = MixedElement([V, S])
W = FunctionSpace(mesh, Z)
W1 = FunctionSpace(submesh_OT, Z)
T = FunctionSpace(submesh_torus, S)

# -------------------------------------------------
# Meshfunction over facets
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("boundaries.xdmf") as bndfile:
    bndfile.read(mvc, "bnd_marker")
mff = MeshFunction('size_t', mesh, mvc)
for i in range(mff.size()):
    if mff[i] not in (3,4):
        mff[i] = 0
File("mff.pvd") << mff
dS = Measure('dS', domain=mesh, subdomain_data=mff)

#Boundary \Omega
w0 = Function(W)
bc0 = DirichletBC(W, w0, mff, 4)
#Boundary T
w2 = Function(T)
bc2 = DirichletBC(W, w0, mff, 3)

bc = [bc0, bc2]

With a suitable set of filters in Paraview, I obtain the following for mff.pvd:
mff.pvd
and the following for mft.pvd:
mft.pvd

2 Likes