BoundaryMesh(mesh, "interior", True) is returning an empty set

1

In the below MWE, I have an annulus with inner radius 1 and outer radius 2. I introduced a mesh on the annulus. I want to extract the coordinates of the nodes that lie on the inner circle. But

BoundaryMesh(mesh, "interior", True)

is returning an empty set. According to me, the above command must return the nodal coordinates on the inner circle.

MWE:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
from dolfin import *
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File

C1 = Circle(Point(0, 0),1)
C2 = Circle(Point(0, 0),2)

C_1 = C2 - C1

class b1(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 1, eps=1e-3) and on_boundary
class b2(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 4, eps=1e-3) and on_boundary

class s1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 > 1 and x[0]**2 + x[1]**2 < 4

mesh = generate_mesh(C_1, 5)

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
b1().mark(boundary_markers, 1)
b2().mark(boundary_markers, 2)
s1().mark(surface_markers, 5)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)

inner_nodes = BoundaryMesh(mesh, "interior", True)
boundary_coordinates = inner_nodes.coordinates()
print(len(boundary_coordinates))
print(boundary_coordinates)

I also tried the commands BoundaryMesh(mesh, "local", True). This returns all the coordinates of the nodes both on the inner and the outer boundaries of the annulus. So how can I extract the nodal coordinates of just the nodes on the inner circle?

Any suggestions would be very helpful.

2

I would strongly suggest using an exterior mesh generator such as gmsh to generate your meshes. Then you can tag interior boundaries in the mesh generation, and import them to dolfin as xml or XDMF.

I would also suggest you have a look at the documentation, as from it is clear that “interior” is not what you want to use: https://fenicsproject.org/olddocs/dolfin/latest/cpp/d3/dcc/classdolfin_1_1BoundaryMesh.html#a864bd8b53c0588ba2a077feb4d500973

( std::string ) The type of BoundaryMesh, which can be “exterior”, “interior” or “local”. “exterior” is the globally external boundary, “interior” is the inter-process mesh and “local” is the boundary of the local (this process) mesh.

3

Hi @dokken ,

I always use gmsh for creating the mesh. I created an in-script mesh just for the MWE to demonstrate the problem.

I usually import the msh file from gmsh to dolfin and create xdmf for the surface (using surface markers) and the boundaries (using boundary markers). But I am not sure how to extract the boundary nodes from the xdmf itself.

I tried to do it using a vertex_to_dof_map in the following example

Consider the same annulus created in gmsh, with the gmsh script:

SetFactory("OpenCASCADE");
//+
Circle(1) = {0, 0, 0, 2, 0, 2*Pi};
//+
Circle(2) = {0, 0, 0, 1, 0, 2*Pi};
//+
Curve Loop(1) = {1};
//+
Curve Loop(2) = {2};
//+
Plane Surface(1) = {1, 2};
//+
Curve Loop(3) = {1};
//+
Physical Curve("C1", 1) = {2};
//+
Physical Curve("C2", 2) = {1};
//+
Physical Surface("S1", 3) = {1};

The Python script for extracting the inner boundary coordinates is:

from dolfin import *     
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp
import numpy as np
import matplotlib.pyplot as plt

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

msh = meshio.read("bndycoord.msh")  # Inner circle has 12 nodes. Outer circle has 23 nodes. Interior has 24 nodes
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)

mvc1 = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc1)
sur_mesh = cpp.mesh.MeshFunctionSizet(mesh, mvc1)

mvc2 = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc2)
bnd_mesh = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
###############################################################
File('boundary1.pvd')<<bnd_mesh

V = FunctionSpace(mesh,'P',1)

vx = project(Expression("x[0]",degree=1),V)
vy = project(Expression("x[1]",degree=1),V)
vz = project(Expression("x[2]",degree=1),V)

v2d = vertex_to_dof_map(V)
dofs =[]
cx=[]
cy=[]
cz=[]
for facet in facets(mesh):
    if bnd_mesh[facet.index()]==3: # The surface has the physical tag 3
        vertices = facet.entities(0)
        for vertex in vertices:
            dofs.append(v2d[vertex])

unique_dofs = np.array(list(set(dofs)),dtype=np.int32)
boundary_coords = V.tabulate_dof_coordinates()[unique_dofs]

for i, dof in enumerate(unique_dofs):
    print(vx.vector()[dof][0],vy.vector()[dof][0],vz.vector()[dof][0])
    
print('number of boundary coordinates =', len(vx.vector())) # Gives 59: all the mesh nodes

What should be different here in order to find just the node coordinates on the inner circle?

Also, suppose I have a solution u on the annulus. Is it possible to find the values of u at the nodes on the inner circle without explicitly plugging in the coordinates of the nodes on the inner circle?

4

Consider the following (using the mesh that you supplied above):

from dolfin import *
import meshio
import numpy as np


def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={
                           cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    return out_mesh


# Inner circle has 12 nodes. Outer circle has 23 nodes. Interior has 24 nodes
msh = meshio.read("annulus.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mf.array()[mf.array() > 1e3] = 0
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)

File('boundary.pvd') << mf

V = FunctionSpace(mesh, 'P', 1)

vx = project(Expression("x[0]", degree=1), V)
vy = project(Expression("x[1]", degree=1), V)
vz = project(Expression("x[2]", degree=1), V)

v2d = vertex_to_dof_map(V)
dofs = []
cx = []
cy = []
cz = []
for facet in facets(mesh):
    if mf[facet.index()] == 1:  # The inner circle has boundary 1
        vertices = facet.entities(0)
        for vertex in vertices:
            dofs.append(v2d[vertex])

unique_dofs = np.array(list(set(dofs)), dtype=np.int32)
boundary_coords = V.tabulate_dof_coordinates()[unique_dofs]
u = interpolate(Expression("x[0]", degree=1), V)
print(boundary_coords, boundary_coords.shape)
print(u.vector().get_local()[unique_dofs])
1 Like