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.

Could you please help?

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)

r_mid = Ri + (Re - Ri) / 2

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

class outer_boun(SubDomain):
    def inside(self, x, on_boundary):
        return sqrt(x[0]**2 + x[1]**2 + x[2]**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"])

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.