Why are boundary and surface markers not carried over to the refined mesh?

1

Hi! I am having trouble in carrying over the boundary and the surface tags to a refined mesh. I tried some of the suggestions in previous posts but couldn’t make them work. In the following MWE, I created and meshed a disk with the boundary tagged 2 and the surface of the disk tagged 1. Then I did

mesh1 = refine(mesh)

But the refined mesh doesn’t preserve the tags of the coarser mesh.

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp

# 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}

import numpy as np
import meshio
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
import numpy as np

        
class boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary        


class disk(SubDomain):
    def inside(self, x, on_boundary):
        return True

C = Circle(Point(0, 0), sqrt(3))

mesh = generate_mesh(C, 100)

boundary_nodes = BoundaryMesh(mesh, "exterior",  True)
boundary_coordinates = boundary_nodes.coordinates()
print(len(boundary_coordinates))

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
boundary().mark(boundary_markers, 2)
disk().mark(surface_markers, 1)

mesh1 = refine(mesh)
#bmarkers = adapt(mesh1,boundary_markers)
#smarkers = adapt(mesh1,surface_markers)
boundary_nodes = BoundaryMesh(mesh1, "exterior",  True)
boundary_coordinates = boundary_nodes.coordinates()
print(len(boundary_coordinates))

ds = Measure('ds', subdomain_data=boundary_markers)
dx = Measure('dx', subdomain_data=surface_markers)
n = FacetNormal(mesh1)

W1 = FunctionSpace(mesh1, "Lagrange", 1)
n = FacetNormal(mesh1)
G , mu = 1, 0.1
u_D=Constant(0.0)

bc = DirichletBC(W1, u_D, boundary_markers, 2)

# Define variational problem
u = TrialFunction(W1)
v = TestFunction(W1)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = -f*v*dx 

# Compute solution
u = Function(W1)
solve(a == L, u, bc)

Error:

Error:   Unable to create Dirichlet boundary condition.
*** Reason:  User MeshFunction and FunctionSpace meshes are different.

I also tried (but didn’t work!)

bmarkers = adapt(mesh1,boundary_markers)
smarkers = adapt(mesh1,surface_markers)

hoping that boundary and surface markers of the coarser mesh will be carried over to the refined one. Also according to python.org, adapt(X,Y)" means “continue by using X according to protocol Y” (Reference: https://www.python.org/dev/peps/pep-0246/)

So is it that every time I refine a mesh, I have to redefine the markers? Is there a way to preserve the tags of the coarser mesh?

2

Please consider the following example:


from dolfin import *
from mshr import *
import numpy as np
parameters["refinement_algorithm"] = "plaza_with_parent_facets"
        
class boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary        


class disk(SubDomain):
    def inside(self, x, on_boundary):
        return True

C = Circle(Point(0, 0), sqrt(3))

mesh = generate_mesh(C, 100)

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)

boundary().mark(boundary_markers, 2)
disk().mark(surface_markers, 1)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)
print(assemble(1*dx(1)), assemble(1*ds(2)))



mesh1 = adapt(mesh)
bmarkers = adapt(boundary_markers, mesh1)
smarkers = adapt(surface_markers, mesh1)

ds1 = Measure('ds', domain=mesh1, subdomain_data=bmarkers)
dx1 = Measure('dx', domain=mesh1, subdomain_data=smarkers)
print(assemble(1*dx1(1)), assemble(1*ds1(2)))

File("ref_marker.pvd") << smarkers
File("marker.pvd") << surface_markers

Note that you are not referring to the FEniCS (dolfin) adapt function with your documentation. adapt is a function supplied by dolfin.
Please also reduce your example next time, by removing all unnecessary imports and code called after you receive an error message.

1 Like
3

Thanks @dokken ! I will try it out. And my apologies for not trimming down the minimal code. But doesn’t an MWE require a minimal code that can be copied and executed directly?

Just to make sure I supply the minimal code next time, what do you mean by “removing all unnecessary imports and code called after you receive an error message”? Does it mean removing lines like

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp

# 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}

and

bmarkers = adapt(mesh1,boundary_markers)
smarkers = adapt(mesh1,surface_markers)
4

For instance consider the code I supplied above. It only has three import statements, it only has the functionality required to check that the mesh and markers have been refined.

In the code you supplied you import things from dolfin in three different ways:

  1. from dolfin import *
  2. from fenics import *
  3. from dolfin import Mesh, ....

You should only use one of these, and remove imports such as

  1. from __future__ import print_function (not used in the code)
  2. import matplotlib.pyplot as plt (not used in the code)
  3. import meshio (not used in the code)

You also define a BoundaryMesh that is unused, the FacetNormal is defined twice, etc.

You example above, with the error message you supplied should stop after:

i.e. not include

(and therefore definitions of G, mu, u_D can be removed as well.

1 Like
5

Got it! Will do that from next time on. Thanks for the explanation.

6

Hi @dokken !

The code worked. But I have a question:

What does parameters["refinement_algorithm"] = "plaza_with_parent_facets" do? In the link below from another post, by “refine the facet markers appropriately”, do you mean preserving facet markers of the original coarse mesh? I don’t understand what refining a maker means?

7

Plaza with parent facets means that the plaza algorithms stores the parent of a facet when it is refined. Imagine the outer boundary of your initial mesh. When you adapt/refine the mesh, each existing facet is divided into 2. you Need to keep this information to transfer markers from a coarse mesh to a finer one.

1 Like
8

Got it! Thanks for the explanation.