Hi all,
I'm trying to mark surfaces of a box for FEM simulations with DG elements.
Something however appears to go wrong. Below, you can find a minimal example to reproduce the problem.
Running the code with
myside = CompiledSubDomain("near(x[0], X)", X = 0.0) (left picture)
and
myside = CompiledSubDomain("near(x[0], X)", X = 1.0) (right picture)
results in two different pictures, which (in my opinion) should be identical (up to the orientation of the x-axis). I don't think the problem is in de marking, since when I plot(boundaries) I do get identical pictures as expected.
Experimenting further with other sides of the box results in similar results.
When considering the side x[2] == 3., only one node appears to be affected...
Am I doing something wrong? Is there a way to circumvent/overcome this?
Thanks in advance
Steven
from dolfin import *
mesh = BoxMesh(0., 0., 0., 1., 2., 3., 2, 4, 6)
# Define the boundaries
myside = CompiledSubDomain("near(x[0], X)", X = 0.0)
#myside = CompiledSubDomain("near(x[0], X)", X = 1.0)
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
myside.mark(boundaries, 1)
ds = Measure("ds")[boundaries]
# Define functionspaces
U = FunctionSpace(mesh, "DG", 1)
# Variational formulation
u = TrialFunction(U)
q = TestFunction(U)
F = inner(u, q)*dx - inner(1., q)*ds(1)
a, L = lhs(F), rhs(F)
# solve
u = Function(U)
solve(a == L, u)
File("solu.pvd")<<u
