With Dolfin version 2019.1.0 I in the Docker environment, I experience problems with correctly assigning mesh vertices to a sub-domain. It seems the edges get assigned all-right. But whether I use dolfin.RectangleMesh or mshr.Rectangle, a lot of boundary vertices clearly matching my criteria do not get assigned the marked correctly.
The MWEs below should produce a rectangular domain and assign a MeshFunction value of 1 to all elements on the outer boundaries of the domain while the inside should be assigned 0. Both vertices and edges seem to be affected, but vertices show this problem much more often (you’ll barely spot the affected edges at first glance).
Here’s my MWE for dolfin.RectangleMesh:
import mshr
import dolfin as do
import fenics as fe
import numpy as np
import matplotlib.pyplot as plt
print(do.__version__)
X, Y = 7e-3, 3.5e-3
mesh = do.RectangleMesh(do.Point(0,0), do.Point(X,Y), 20, 10)
do.plot(mesh)
FE = do.FiniteElement("P", mesh.ufl_cell(), 1)
FS = do.FunctionSpace(mesh, FE)
sub_domains = do.MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
sub_domains.set_all(0)
class thePEC(do.SubDomain):
def inside(self, x, on_boundary):
val = do.near(x[0], 0.) | do.near(x[0], X) | do.near(x[1], 0.) | do.near(x[1], Y) # & on_boundary
return val
outerwall = thePEC()
outerwall.mark(sub_domains, 1)
bcs = [
do.DirichletBC(FS, do.Constant(0.0), sub_domains, 1),
]
do.File("bcs2.pvd") << sub_domains
Here’s my MWE for mshr.Rectangle:
import mshr
import dolfin as do
import fenics as fe
import numpy as np
import matplotlib.pyplot as plt
print(do.__version__)
X, Y = 7e-3, 3.5e-3
dom = mshr.Rectangle(do.Point(0,0), do.Point(X, Y))
mesh = mshr.generate_mesh(dom, 20)
do.plot(mesh)
FE = do.FiniteElement("P", mesh.ufl_cell(), 1)
FS = do.FunctionSpace(mesh, FE)
sub_domains = do.MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
sub_domains.set_all(0)
class thePEC(do.SubDomain):
def inside(self, x, on_boundary):
val = do.near(x[0], 0.) | do.near(x[0], X) | do.near(x[1], 0.) | do.near(x[1], Y) # & on_boundary
return val
outerwall = thePEC()
outerwall.mark(sub_domains, 1)
bcs = [
do.DirichletBC(FS, do.Constant(0.0), sub_domains, 1),
]
do.File("bcs1.pvd") << sub_domains
Checking the output files (bcs1.pvd and bcs2.pvd) against each other e.g. using “Representation: Points” in Paraview reveals the problem.
It seems when using a circle things work better…
Any ideas?
