1D boundary integrals and marking


253
views
3
11 months ago by
Hello everyone,

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], 0.)
class boundaryRight(SubDomain):
	def inside (self, x, on_boundary):
		return on_boundary and near(x[0],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:

sub_domains.set_all(0)


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.

Thank you!

Cheers, Petr.

PS. MWE: https://gist.github.com/Corwinpro/0c32046e56023e1e1dd26dc4b57360bb
Community: FEniCS Project
Answer: ds should be returned in def Markboundaries(). Otherwise, it is just a local instance.
written 11 months ago by Corwinpro  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »