 # Inner and outer boundaries of spherical annulus getting marked incorrectly; behavior seems to depend on mesh resolution

Dear Community

I am trying to create a 3D spherical anulus and marking the inner and outer boundaries as shown in the code below. I want to check if the facets have been marked correctly.

For an external radius of Re=5 and mesh resolution of 32, it is clearly seen that the boundaries have been marked correctly:

However, setting Re=10 and using a mesh resolution of 32 results in an incorrect identification of the boundaries as shown below:

Increasing the mesh resolution to 64, for Re=10, however, results in the boundaries being marked correctly:

I would like to understand this behavior, and learn if it is possible for the boundaries to get marked correctly irrespective of the mesh resolution.

Thank You
Warm Regards

Code segment

``````import numpy as np
from dolfin import *
import numpy as np
from mshr import *

Re = 10. # parameter under discussion
Ri = 1.

bound = Sphere(Point(0., 0., 0.), Re)
ball = Sphere(Point(0., 0., 0.), Ri)
domain = bound - ball

mesh_resolution = 32 # parameter under discussion
mesh = generate_mesh(domain, mesh_resolution)

# Making a mark dictionary
# Note: the values should be UNIQUE identifiers.
mark = {"generic": 0,
"inner": 1,
"outer": 2}

subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
subdomains.set_all(mark["generic"])

r_mid = Ri + (Re - Ri) / 2

class inner_boun(SubDomain):
def inside(self, x, on_boundary):
return sqrt(x**2 + x**2 + x**2) < r_mid and on_boundary

class outer_boun(SubDomain):
def inside(self, x, on_boundary):
return sqrt(x**2 + x**2 + x**2) > r_mid and on_boundary

## Marking boundaries
in_boun = inner_boun()
in_boun.mark(subdomains, mark["inner"])
out_boun = outer_boun()
out_boun.mark(subdomains, mark["outer"])

geom_name="geometry.pvd"
File(geom_name) << subdomains
``````

I would consider moving to gmsh instead of mshr, as mshr is no longer maintained (Different behaviour of a mesh coming from generate_mesh() resp. UnitSquareMesh()? - #5 by dokken)
I think then your geometry is much cleaner, especially for coarser meshes.

1 Like

Thank you, I will go through the documentation of gmsh and try to implement the mesh generation using that.