Doubts about the demo_overlapping-regions.py example

1

Esteemed,

I would like to know what exactly fenics does when the line below of code is called:

# Define regions as tuples of subdomain labels
DL, DM, DR = (1,2), (2,), (2,3) # ***

How does fenics do this “overlap” over domains?

I see that there is a spacing between the domains in the example, and if it were really considered an overlap of the domains, would there be any difference? As follows:

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] <= 0.25+DOLFIN_EPS)

class Mid(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] >= 0.2-DOLFIN_EPS) and (x[0] <= 0.6+DOLFIN_EPS)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] >= 0.55-DOLFIN_EPS)

Follow the link code below:
https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/undocumented/overlapping-regions/

Could the code authors help me?

2

FEniCS is open source. You are welcome to read the implementation details. Many of the CutFEM papers by Hansbo and Logg provide even greater detail. @dokken may know more literature sources.

3

You cannot mark an overlap in a single mesh function (as each cell only can contain one marker).
You need to assemble each domain, and the add together the matrices, as shown in the example below:

import matplotlib.pyplot as plt
from dolfin import *
from pprint import pprint


class Left(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] <= 2/3+DOLFIN_EPS)


class Right(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] >= 1/3-DOLFIN_EPS)


# Define mesh and subdomains
mesh = UnitSquareMesh(3, 1)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
Left().mark(domains, 1)
domains2 = MeshFunction("size_t", mesh, mesh.topology().dim())
domains2.set_all(0)
Right().mark(domains2, 2)

# Define new measures associated with the interior domains
dx = Measure("dx", domain=mesh, subdomain_data=domains, subdomain_id=1)
dx2 = Measure("dx", domain=mesh, subdomain_data=domains2, subdomain_id=2)
V = FunctionSpace(mesh, "CG", 1)
u, v = TrialFunction(V), TestFunction(V)
a1 = u*v*dx
A1 = assemble(u*v*dx)
a2 = u*v*dx2
A2 = assemble(a2)
A12 = A1 + A2
print("-"*5, "A1", "-"*5)
pprint(A1.array())
print("-"*5, "A2", "-"*5)
pprint(A2.array())
print("-"*5, "A1+A2", "-"*5)
pprint(A12.array())
plot(domains)
plt.savefig("domain.png")
plot(domains2)
plt.savefig("domain2.png")
1 Like