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

mesh
1

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:

re_5_slice_mesh_resol_32
re_5_slice_mesh_resol_321273×831 119 KB

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

re_10_slice_mesh_resol_32
re_10_slice_mesh_resol_321273×831 119 KB

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

re_10_slice_mesh_resol_64
re_10_slice_mesh_resol_641273×831 124 KB

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)
subdomains.set_all(mark["generic"])

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"])

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

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
3

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