1D boundary integrals and marking
I'm trying to implement a Robin boundary condition for a toy-problem in 1D.
I have a 1D domain and thus two boundaries:
mesh = IntervalMesh(100, 0., 1.) class boundaryLeft(SubDomain): def inside (self, x, on_boundary): return on_boundary and near(x, 0.) class boundaryRight(SubDomain): def inside (self, x, on_boundary): return on_boundary and near(x,1.) boundaryLeft = boundaryLeft(); boundaryRight = boundaryRight(); def MarkBoundaries(mesh): sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim()- 1) boundaryLeft.mark(sub_domains,1) boundaryRight.mark(sub_domains,2) ds = Measure('ds', domain=mesh, subdomain_data=sub_domains)
The Left Boundary is marked as '1', and the Right boundary is marked as '2'
Now, the problem arises when I try to construct a weak form (this is just a general idea, not the full code):
p = TrialFunction(V) q = TestFunction(V) a = -inner(grad(p), grad(q))*dx + (q*p)*ds
As you can see, I have a boundary term "q*p*ds" here. If I use ds(1) or ds(2) or (ds(1) + ds(2)) instead of ds, it seems like the boundary term just doesn't appear in the weak form, i.e. the multiplication by ds(1),ds(2) just zeroes the boundary term.
UPD: I also tried VertexFunction, FacetFunction, etc. However, these don't work as well. Even if I just write:
and then use ds(0) instead of ds, it zeroes the boundary term.
My question is:
Why do I get different answers (or different weak forms) if I use (ds(1) + ds(2)) or ds(0) instead of ds in my weak formulation ?
It seems like I do something wrong with the definition of MeshFunction, but I can't explain why and how to avoid this.
PS. MWE: https://gist.github.com/Corwinpro/0c32046e56023e1e1dd26dc4b57360bb