How to correctly implement non homogeneous neumann condtion on a hole?

variational formulation
1

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

2

Edit:

  • \S=S
  • At the end of the post: \int_{\partial S} v.n_x~d\sigma (I forgot the “” for all my integral)
3

Hi,

Could you please supply the value that you are using for comparasion. It is true that using ds will put your Neumann condition on the full boundary. However, even if you have a non-homogenous Neumann-boundary condition on parts of the boundary and a homogenous one on the other part, you still have to constrain your solution. The code will in both cases (using a subdomain and not using a subdomain), give you a solution, but it is not the solution you are looking for.

I’ve reimplemented your problem, using Lagrange multipliers (Using the strategy described at: https://fenicsproject.org/docs/dolfin/2018.1.0/python/demos/neumann-poisson/demo_neumann-poisson.py.html). Does this code give you your reference values?
This code was run using docker:

 docker run -ti -v $(pwd):/home/fenics/shared quay.io/fenicsproject/stable:latest


from dolfin import *

# Generate Mesh using MSHR
from mshr import *
r0 = Rectangle(Point(0.,0.), Point(1, 1))
r1 = Rectangle(Point(0.25,0.25), Point(0.75, 0.75))
domain = r0 - r1
mesh = generate_mesh(domain, 10)

class InnerBoundary(SubDomain):
    """
    The inner boundaries of the mesh
    """
    def inside(self, x, on_boundary):
        return (x[0] > DOLFIN_EPS and x[0] < 1-DOLFIN_EPS and
                x[1] > DOLFIN_EPS and x[1] < 1-DOLFIN_EPS and
                on_boundary)



# Define facet function and integral measure for the inner boundary
facet_func = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
inner_b = InnerBoundary()
inner_b.mark(facet_func, 2)
dInner = Measure("ds", subdomain_data=facet_func, subdomain_id=2)

# Define mixed space for Poisson problem with neumann conditions,
# see: https://fenicsproject.org/docs/dolfin/2018.1.0/python/demos/neumann-poisson/demo_neumann-poisson.py.html
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
R = FiniteElement("Real", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, P1 * R)

# Define variational form
(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)

f = Constant(0)
n = FacetNormal(mesh)
a = inner(grad(u), grad(v))*dx + (c*v + u*d)*dx
L = inner(f,v)*dx - n[0]*v*dInner

uh = Function(W)
solve(a == L, uh)
(u, c) = uh.split(deepcopy=True)

J = assemble(Constant(1/0.75)*u*dInner)
print(J)
File("tmp/u.pvd").write(u)
6

So, I finally implement the code given by dokken (with my periodic condition and with some new scale for my geometry):

from future import print_function
from dolfin import *
from fenics import *
from dolfin.cpp.mesh import *
from mshr import*

r = Rectangle(Point(0,0), Point(0.1, 0.1))
r1 = Rectangle(Point(0.025,0.025), Point(0.075, 0.075))
domain = r-r1
mesh = generate_mesh(domain, 50)

class PeriodicBoundary(SubDomain):

def inside(self, x, on_boundary):
    return bool((near(x[0], 0) or near(x[1], 0)) and 
            (not ((near(x[0], 0) and near(x[1], 0.1)) or 
                    (near(x[0], 0.1) and near(x[1], 0)))) and on_boundary)

def map(self, x, y):
    if near(x[0], 0.1) and near(x[1], 0.1):
        y[0] = x[0] - 0.1
        y[1] = x[1] - 0.1
    elif near(x[0], 0.1):
        y[0] = x[0] - 0.1
        y[1] = x[1]
    else:  
        y[0] = x[0]
        y[1] = x[1] - 0.1

pbc = PeriodicBoundary()

class InnerBoundary(SubDomain):

def inside(self, x, on_boundary):
    return (x[0] > DOLFIN_EPS and x[0] < 0.1-DOLFIN_EPS and
            x[1] > DOLFIN_EPS and x[1] < 0.1-DOLFIN_EPS and
            on_boundary)

facet_func = MeshFunction(“size_t”, mesh, mesh.topology().dim()-1)
inner_b = InnerBoundary()
inner_b.mark(facet_func, 2)
dInner = Measure(“ds”, subdomain_data=facet_func, subdomain_id=2)

P1 = FiniteElement(“Lagrange”, mesh.ufl_cell(), 1)
R = FiniteElement(“Real”, mesh.ufl_cell(), 0)
TH = P1*R
W = FunctionSpace(mesh, TH, constrained_domain= pbc)

(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)

f = Constant(0)
n = FacetNormal(mesh)
a = inner(grad(u), grad(v))dx + (cv + u*d)*dx
L = inner(f,v)*dx - n[0]vdInner

uh = Function(W)
solve(a == L, uh)
(u, c) = uh.split()

J = assemble((Constant(1/0.0075)un[0])*dInner)
print(J)

Edit: It’s working!!!

(I made a muistake in my scaling of the geometry and in my new InnerBoundary Class)
We obtain -0.22855 for the integral \frac{1}{0.0075}\int_{\partial [0.025,0.075]^2} bn_x~d\sigma VS -0.2299 in COMSOL.

CapturePbFEniCS2
CapturePbFEniCS2.jpg3200×1080 449 KB

Thank you for your help