Hello everyone,
First of all, I’m new to FEniCS (coming from COMSOL).
I consider a perforated domain [0,1]^2\setminus [0.25,0.75]^2=\Omega \setminus \S and I have some difficulty when I want to implement the non homogeneous Neuman condtion \partial_n u1=-n_x only on the boundary of the little square S (where n_x is the x component of the normal)
In fact, I do not really know what represents ds and dS (internal and external facet).
For that, I do not know how to get my flow condition on \partial S but also how to perform the line integration \frac{1}{|\Omega \setminus \S |} int_{\partial S} b_1n_x d\sigma
I tried two things:
i) Weak formulation using « ds » and without declaring a « Subdomain » class for the obstacle
r=Rectangle(Point(0.,0.), Point(1, 1))
r1= Rectangle(Point(0.25,0.25), Point(0.75, 0.75))
domain = r-r1)
u1 = TrialFunction(V)
v1 = TestFunction(V)
n = FacetNormal(mesh)
f = Constant(0.0)
a1 = inner(grad(u1), grad(v1))dx
L1 = fv1*dx-n[0]v1ds
integral = (1/0.75)assemble(u1ds)
With this technique, I find the good value for \frac{1}{|=\Omega \setminus \S |} int_{\partial S} b_1.n_x (comparison with COMSOL). But I know that this code does not represent my real mathematical problem. Indeed, this system has an infinity of solution. So, it needs a constraint (pointwise constraint by lagrange multiplier or a zero mean value over the domain). Yet, the code doesn’t need it. That give me the impression that « ds » put the flux on \partial S but also on \partial \Omega.
So I tried the second method: create a class(SubDomain) for the obstacle and news measures “ds()”
class gauche(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.25) and between(x[1],(0.25,0.75))
class droite(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.75) and between(x[1],(0.25,0.75))
class bas(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.25) and between(x[0],(0.25,0.75))
class haut(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.75) and between(x[0],(0.25,0.75))
gauche = gauche()
haut = haut()
droite = droite()
bas = bas()
boundary_parts = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1)
gauche.mark(boundary_parts, 1)
haut.mark(boundary_parts, 2)
droite.mark(boundary_parts, 3)
bas.mark(boundary_parts, 4)
ds = Measure(“ds”)[boundary_parts]
L2 = fv2dx - n[0]v2ds(1)-n[0]v2ds(2)-n[0]v2ds(3)-n[0]v2ds(4)
integral = (1/0.75)assemble(u1ds(1))
integral1 = (1/0.75)assemble(u1ds(2))
integral2 = (1/0.75)assemble(u1ds(3))
integral3 = (1/0.75)assemble(u1ds(4))
Here, I don’t obtain the good result.
Can you tell me how to write correctly my weak formulation under FEniCS i.e how to had int_{\partial_S} v.n_x in my linear form L(v)?
Thank you