Local Refinement of SubDomain

[uday-rade]

Hi Everyone,

I want to refine the mesh in the solid domain of the fluid structure interaction problem. For that, I am using a refine() function and a boolean mesh function. However, I believe I am losing markers, and I need to use adapt() because my relative and absolute residuals have become nan. But I’m not sure how or where to use adapt(). Please help me figure out how to locally refine it.

from dolfin import *
from ufl import indices, Jacobian, Min
from mshr import *

# Parameters defining domain geometry:
SOLID_LEFT = 0.45
SOLID_RIGHT = 0.55
SOLID_TOP = 0.5
OMEGA_H = 0.75
OMEGA_W = 1.0

# Define the mshr geometrical primitives for this domain:
r_Omega = Rectangle(Point(0,0),Point(OMEGA_W,OMEGA_H))
r_Omega_s = Rectangle(Point(SOLID_LEFT,0),
                      Point(SOLID_RIGHT,SOLID_TOP))


SOLID_FLAG = 1
r_Omega.set_subdomain(SOLID_FLAG,r_Omega_s)

# Define subdomains for use in boundary condition definitions:
class Walls(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary 
                and ((x[1] < DOLFIN_EPS) 
                     or (x[1] > (OMEGA_H - DOLFIN_EPS))))
class Inflow(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and (x[0] < DOLFIN_EPS))
class Outflow(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and (x[0] > (OMEGA_W - DOLFIN_EPS)))
class SolidDomainClosure(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] > SOLID_LEFT - DOLFIN_EPS
                and x[0] < SOLID_RIGHT + DOLFIN_EPS
                and x[1] < SOLID_TOP + DOLFIN_EPS)

mesh = generate_mesh(r_Omega, 70)
d = mesh.geometry().dim()
DomainToRefine= CompiledSubDomain("x[0] > 0.45 - DOLFIN_EPS and x[0] < 0.55 + DOLFIN_EPS and x[1] < 0.5 + DOLFIN_EPS")

 r_markers = MeshFunction("bool", mesh, d, False)
DomainToRefine.mark(r_markers, True)
refinedMesh = refine(mesh,r_markers)

mesh=refinedMesh
markers = MeshFunction('size_t', mesh, d, mesh.domains())

bdry_markers = MeshFunction('size_t', mesh, d-1,0)
OUTFLOW_FLAG = 2
Outflow().mark(bdry_markers,OUTFLOW_FLAG)

dx = dx(metadata={'quadrature_degree': 2},
        subdomain_data=markers)
ds = ds(metadata={'quadrature_degree': 2},
        subdomain_data=bdry_markers)

cell = mesh.ufl_cell()
Ve = VectorElement("CG", cell, 1)
Qe = FiniteElement("CG", cell, 1)
VQe = MixedElement((Ve,Qe))
W = FunctionSpace(mesh,VQe)
V = FunctionSpace(mesh,Ve)

[dokken]

First please use 3x` encapsulation to format the code, ie

```python
#add code here
```

Secondly, consider

Adapt for mesh and MeshFunctions

The issue is that when refining the mesh, new facets are introduced into the mesh (that had no previous parent facet, as they are internal to a cell, these facets get large numbers as values. I would like to note that your example is a very good one (good way of presenting a minimal working example), but I have simplified it slightly for easier visual interpretation: from dolfin import * parameters["refinement_algorithm"] = "plaza_with_parent_facets" class Top(SubDomain): def inside(self,x,…

[dokken]

You have not defined boundary markers (bdry_markers), so I cannot really help you any further with your post.
Also I don’t see the point of creating markers, as the original mesh as there are no markers of dimension mesh.geometry().dim().

[uday-rade]

Hi Dokken,
Thank you for your response!

I have updated the code with some other subdomains if you looking for them.

Outflow() is needed to apply stable Neumann BC i.e. zero traction BC at the outlet. That is why I am marking Outflow(). Please let me know if you need further information.

[dokken]

@uday-rade

bdry_markers

is still not defined.

[uday-rade]

Sorry, just updated!

[dokken]

I am not sure what your issue is, as the following code runs as expected:

from dolfin import *
from ufl import indices, Jacobian, Min
from mshr import *

# Parameters defining domain geometry:
SOLID_LEFT = 0.45
SOLID_RIGHT = 0.55
SOLID_TOP = 0.5
OMEGA_H = 0.75
OMEGA_W = 1.0

# Define the mshr geometrical primitives for this domain:
r_Omega = Rectangle(Point(0, 0), Point(OMEGA_W, OMEGA_H))
r_Omega_s = Rectangle(Point(SOLID_LEFT, 0), Point(SOLID_RIGHT, SOLID_TOP))


SOLID_FLAG = 1
r_Omega.set_subdomain(SOLID_FLAG, r_Omega_s)


# Define subdomains for use in boundary condition definitions:
class Walls(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((x[1] < DOLFIN_EPS) or (x[1] > (OMEGA_H - DOLFIN_EPS)))


class Inflow(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (x[0] < DOLFIN_EPS)


class Outflow(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (x[0] > (OMEGA_W - DOLFIN_EPS))


class SolidDomainClosure(SubDomain):
    def inside(self, x, on_boundary):
        return (
            x[0] > SOLID_LEFT - DOLFIN_EPS
            and x[0] < SOLID_RIGHT + DOLFIN_EPS
            and x[1] < SOLID_TOP + DOLFIN_EPS
        )


outflow = Outflow()
OUTFLOW_FLAG = 2


mesh = generate_mesh(r_Omega, 70)
d = mesh.topology().dim()
# Generate initial outflow marker
boundary_marker = MeshFunction("size_t", mesh, d - 1, 0)
outflow.mark(boundary_marker, OUTFLOW_FLAG)
File("initial_marker.pvd") << boundary_marker

# Refine mesh
DomainToRefine = CompiledSubDomain(
    "x[0] > 0.45 - DOLFIN_EPS and x[0] < 0.55 + DOLFIN_EPS and x[1] < 0.5 + DOLFIN_EPS"
)
r_markers = MeshFunction("bool", mesh, d, False)
DomainToRefine.mark(r_markers, True)
refinedMesh = refine(mesh, r_markers)
mesh = refinedMesh


# Create new boundary marker
boundary_marker = MeshFunction("size_t", mesh, d - 1, 0)
outflow.mark(boundary_marker, OUTFLOW_FLAG)
File("new_marker.pvd") << boundary_marker

yielding

[uday-rade]

The issue is that if I don’t locally refine the mesh (uniform refinement), my code works great. However, if locally refine it (mesh=refinedMesh), my relative and absolute residuals attain NaN values. I emailed you; please check it.
Thank you!

[uday-rade]

Actually, the issue is with losing domain markers after refinement. The r_omega_s domain with marker ‘SOLID_FLAG’ is lost after refinement. Please have a look at the code and images.
Thank you!

from dolfin import *
from ufl import indices, Jacobian, Min
from mshr import *
from matplotlib import pyplot as plt
####### Domain and mesh setup #######

# Parameters defining domain geometry:
SOLID_LEFT = 0.45
SOLID_RIGHT = 0.55
SOLID_TOP = 0.5
OMEGA_H = 0.75
OMEGA_W = 1.0
REF_MARGIN = 0.1

# Define the mshr geometrical primitives for this domain:
r_Omega = Rectangle(Point(0,0),Point(OMEGA_W,OMEGA_H))
r_Omega_s = Rectangle(Point(SOLID_LEFT,0),
                      Point(SOLID_RIGHT,SOLID_TOP))

# Enumerate subdomain markers
FLUID_FLAG = 0
SOLID_FLAG = 1
# Zero is the default flag, and does not need to be
# explicitly set for the fluid subdomain.
r_Omega.set_subdomain(SOLID_FLAG,r_Omega_s)

# Parameters defining refinement level:
N = 64

# Generate mesh of Omega, which will have a fitted
# subdomain for Omega_s.
mesh = generate_mesh(r_Omega, N)

# Mesh-derived quantities:
d = mesh.geometry().dim()
markers = MeshFunction('size_t', mesh, d, mesh.domains())
File("initial_marker.pvd") << markers
# Refine mesh
DomainToRefine = CompiledSubDomain(
    "x[0] > 0.45 - DOLFIN_EPS and x[0] < 0.55 + DOLFIN_EPS and x[1] < 0.5 + DOLFIN_EPS"
)
r_markers = MeshFunction("bool", mesh, d, False)
DomainToRefine.mark(r_markers, True)
refinedMesh = refine(mesh, r_markers)
mesh = refinedMesh
markers = MeshFunction('size_t', mesh, d, mesh.domains())
File("after_refinement_marker.pvd") << markers

[uday-rade]

Once mesh is refined with refine(), just add following line after mesh=refinedMesh

markers = adapt(markers,mesh)

So, the complete code will look like

from dolfin import *
from ufl import indices, Jacobian, Min
from mshr import *
from matplotlib import pyplot as plt
####### Domain and mesh setup #######

# Parameters defining domain geometry:
SOLID_LEFT = 0.45
SOLID_RIGHT = 0.55
SOLID_TOP = 0.5
OMEGA_H = 0.75
OMEGA_W = 1.0
REF_MARGIN = 0.1

# Define the mshr geometrical primitives for this domain:
r_Omega = Rectangle(Point(0,0),Point(OMEGA_W,OMEGA_H))
r_Omega_s = Rectangle(Point(SOLID_LEFT,0),
                      Point(SOLID_RIGHT,SOLID_TOP))

# Enumerate subdomain markers
FLUID_FLAG = 0
SOLID_FLAG = 1
# Zero is the default flag, and does not need to be
# explicitly set for the fluid subdomain.
r_Omega.set_subdomain(SOLID_FLAG,r_Omega_s)

# Parameters defining refinement level:
N = 64

# Generate mesh of Omega, which will have a fitted
# subdomain for Omega_s.
mesh = generate_mesh(r_Omega, N)

# Mesh-derived quantities:
d = mesh.geometry().dim()
markers = MeshFunction('size_t', mesh, d, mesh.domains())
File("initial_marker.pvd") << markers
# Refine mesh
DomainToRefine = CompiledSubDomain(
    "x[0] > 0.45 - DOLFIN_EPS and x[0] < 0.55 + DOLFIN_EPS and x[1] < 0.5 + DOLFIN_EPS"
)
r_markers = MeshFunction("bool", mesh, d, False)
DomainToRefine.mark(r_markers, True)
refinedMesh = refine(mesh, r_markers)
mesh = refinedMesh
markers = adapt(markers,mesh)