Rigid Plate Boundary Condition

Hello Everyone,

I pretty much want to do exactly what “SM” was asking in this post: https://fenicsproject.org/qa/12513/how-can-a-rigid-plate-constraint-be-applied-to-a-boundary/ but is unfortunatly without an answer.
I want to verify my Poromechanic Code with the 2D Mandel Problem where the Area is a rectangle with lenght a, and b.

I have found the paper of Castelletto et. al. (https://onlinelibrary.wiley.com/doi/pdf/10.1002/nag.2400) where one Boundary Condition (page 23) is:

u(x,y=b) * n = constant in space

Which means, that the normal Displacement at the Top of the rectangle ist constant in space.
My Question is: how can I implement this in my Code?

Much Thanks in advance for you help!

Hi,
consider

DirichletBC(V.sub(1), Constant(1.), top)

where top is your function for the top boundary and you can change the constant to be any value you like.
Best

1 Like

Hello bleyerj,

thank you for your answer!
But isn’t the displacement on the top then not only constant in space, but also in time?

Hi,
you didn’t specify that you wanted a time-dependent value. Then consider using an Expression such as:

Uimp = Expression("t", t=0, degree=0)
bc = DirichletBC(V.sub(1), Uimp, top)

and later in your code you can assign a specified value for the time t with

Uimp.t = t
1 Like

Thank you, I will try this!

I’m sorry, but I think that isn’t really what I want.

I wasn’t sure myself what exactly the BC meant and because of this, I didn’t asked the right question.
What I really want is: to constrain all degrees of freedom pertaining to u_y on the top boundary to have the
same value, which I don’t know before. The values of the displacement can vary for each time-step.

Is this possible? For Example to declare one node at the top as a “master node/element” and all the other nodes/elements at the top boundary as slave-nodes/Elements that have to have the same displacement as the “master-Element/node”?

OK,
for this you will have to use lagrange multipliers and impose the constraint directly in the variational formulation. However, these Lagrange multipliers shloud live only on the corresponding boundary. This is not possible in the current versions of FEniCS but should be possible in C. Daversin mixed-dimensional branch which will soon be merged in the next version.

2 Likes

Regarding the mixed-dimensional branch, you can either use the dedicated branches (installation from source), or the Docker container named ceciledc/fenics_mixed_dimensional including a 2D-1D demo using a Lagrange multiplier.

1 Like

I’m using linux on ubuntu, can I still use the the Docker container that you linked ( https://hub.docker.com/repository/docker/ceciledc/fenics_mixed_dimensional) or do I have to use one of the dedicated branches from here: https://bitbucket.org/fenics-project/dolfin/branch/cecile/mixed-dimensional ?
Sorry, I’m still new to fenics and Linux.

Hi,

In ubuntu it is also quite convinient to use Docker to isolate different develop environment. That holds for exploring (different versions of) FEniCS as well.

1 Like

Yes, you can use the Docker container on Ubuntu. I would recommend using this option, since all you need is already included and installed in the container.
To do so, just install Docker on your ubuntu (Docker documentation) and run the Docker container ceciledc/fenics_mixed_dimensional

1 Like

Thank you everyone for your help!

I will try to solve this with Docker as you suggested and your Code. :slight_smile:

I downloaded Docker and installed the https://hub.docker.com/repository/docker/ceciledc/fenics_mixed_dimensional

But I don’t really understand how the demo there can help me for my Problem, I guess I don’t understand it well enough. How would you use a lagrange multiplier to implement the Constrain, that the displacement on the top is constant?

And if I have my lagrange multiplier that only lives on the top-boundary, how can I access only the displacements on the top for implementation in my equation?
I guess it has to look something like

k = 10**6 # lagrange multipliers
eqn = … + k*(u(1,0) - u(everything else on the top boundary)) * dS(top)

But I really don’t know how to convert this to a form that I can give to fenics.
I’m sorry for the many questions but this is my first time implementing code and I cant fine a solution on my own.

Having the displacement constant on the top would correspond to a null gradient condition for the displacement on this region.

To have your Lagrange multiplier living only on the top-boundary, you need to define the corresponding submesh and functionspace. Note that if your constrain is a zero gradient, your Lagrange multiplier (and then the function space it belongs to) will be vectorial.

As you might have seen in the demo, you can define such a submesh by marking the appropriate mesh entities and using MeshView

marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
class TopBoundary(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[1], 1) # if top is y = 1 for example

# Marking facets that are on TopBoundary and creating the corresponding submesh
TopBoundary().mark(marker,1)
submesh_top = MeshView.create(marker, 1)

Then you can define the function space for your Lagrange multiplier as

LM = VectorFunctionSpace(submesh_top, ... )

and add this function space in the definition of your MixedFunctionSpace.

When writing your variational form, you need to define the integration domain, especially for the terms involving the Lagrange multiplier that you will integrate only on the top boundary. It might look like that :

dL = Measure("dx", domain=submesh_top)
a = .... + inner(grad(u),e)*dL

You can also have a look to this post involving Lagrange multiplier(s).

1 Like

Thank you so much for your detailed explanation!!!

I tried to implement, what you wrote, and used

Lag = inner(l,v)*dL + inner(u,e)*dL

for the implementation of the lagrangian multiplier like it was done here: Stress continuity using mixed-dimensional branch

But when I try to run it, I get the error, and I can’t figure out what I did wrong.

> Traceback (most recent call last):
>   File "mandel_v2.py", line 142, in <module>
>     solve(l == r, w, bc_ges)    
>   File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/solving.py", line 233, in solve
>     _solve_varproblem(*args, **kwargs)
>   File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/solving.py", line 269, in _solve_varproblem
>     form_compiler_parameters=form_compiler_parameters)
>   File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/problem.py", line 59, in __init__
>     a = Form(a, form_compiler_parameters=form_compiler_parameters)
>   File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/form.py", line 23, in __init__
>     self.subdomains, = list(sd.values())  # Assuming single domain
> ValueError: too many values to unpack (expected 1)

I tried to create a “minimal” example by removing the inital conditions and the plotting etc.

from dolfin import *
import matplotlib.pyplot as plt

# Elasticity parameters 

E= 100*10**6		# E in Pa
nu = 0.25 		# Poissonnumber
k0= 1*10**(-14)
visk = 1 * 10**(-5)
k = k0/visk		# permeability	in m^4/Ns
b = 1           # Biot Koefficent
M = 10**14      # Biot Modul
F = 2        # Force at the top of the plate
a = 1           # length of the plate
G = Constant(E/(2*(1 + nu))) 

T = 1           # final time
num_steps = 1000    # number of time steps
dt = T / num_steps # time step size

# Create mesh and define function space
nx0 = 0
ny0 = 0
nx1 = 1
ny1 = 1
mesh = RectangleMesh(Point(nx0, ny0), Point(nx1, ny1), 40,40)

# Mark boundary subdomians
class rechts(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0],1,tol)
class links(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0],0,tol)
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1],1,tol)
class unten(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1],0,tol)

marker = MeshFunction('size_t',mesh,mesh.topology().dim()-1 ,0)
TopBoundary().mark(marker,1)
submesh_top = MeshView.create(marker,1)

right = rechts()
left = links()
bottom = unten()

# Function spaces
LM = VectorElement('Lagrange', mesh.ufl_cell(), 2)      # Lagrange multiplier
P2 = VectorElement('Lagrange', mesh.ufl_cell(), 2)      # für die Verschiebung quadratisch
P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)      # für den Druck linear

# Create mixed space
TH = MixedElement([P2, P1, LM])
V = FunctionSpace(mesh, TH)

# Define functions
(u, p, l) = TrialFunctions(V)		# Verschiebung, Druck, Lagrangsch-Multiplyer
(v, r, e)  = TestFunctions(V)           
du = Function(V)        # Hilfsgröße zum splitten
u0, p0, l0 = split(du)

# Define Constants
n   = FacetNormal(mesh)
k   = Constant(k)			# Permeabilität 
nu  = Constant(nu)			# Querkontraktionszahl
b = Constant(b)
M = Constant(M)
F = Constant(F)
a = Constant(a)
mu = Constant(E/(2*(1 + nu)))
lam = Constant(E*nu/( (1+nu)*(1-2*nu)))

# Fluidflux
q = 0   # -q_vektor*n= q

# Kinematics
d = u.geometric_dimension()
I = Identity(d)             # Identity tensor

# linearer Straintensor
def epsilon(u):
#    return 0.5*(nabla_grad(u) + nabla_grad(u).T) 
    return sym(nabla_grad(u))
# Spur des Spannungstensors
def tr_epsilon(u):
    return tr(epsilon(u))    

def n_por(u0):      # Porositiy
    return 0.2/(1+tr_epsilon(u0)) 

# linear solid strain
def sigma_s(u):
    return lam*tr_epsilon(u)*Identity(d) + 2*mu*epsilon(u)

# strain from both solid and fluid
def sigma(u,p):
	return (1-n_por(u0)) * sigma_s(u) -n_por(u0)*(p-p0)*I

# Measures
ds_test = ds(domain=mesh, subdomain_data=marker)
dL = Measure('dx', domain=submesh_top)

# Strss-Vector
T1 = Expression(("0", "-F"),degree =2, F=F)

Imp =inner(sigma(u,p),nabla_grad(v))*dx -inner(T1,v)*ds_test(1)  
Fl = - k*inner(nabla_grad((p-p0)),nabla_grad(r))*dx -(tr_epsilon(u) - tr_epsilon(u0))/dt *r *dx
Lag = inner(l,v)*dL + inner(u,e)*dL
eqn = Imp + Fl + Lag

# free-slip BC
bcu_B  = DirichletBC(V.sub(0).sub(1), Constant((0)), bottom)     
bcu_L  = DirichletBC(V.sub(0).sub(0), Constant((0)), left)   
bcp_R  = DirichletBC(V.sub(1), Constant((0)), right)          

bc_ges = [bcu_L, bcu_B, bcp_R]

l = lhs(eqn)
r = rhs(eqn)

# Solve variational problem 
# Step in time
t = 0
      
for n in range(num_steps):
    # Update current time
    t += dt
    w = Function(V)

    solve(l == r, w, bc_ges)    
    (u,p,l) = w.split(True)

    assign(du.sub(0), u)

Hi Martin,

It seems that the definition of the mixed function space is wrong. Try this:

# Function spaces
LM = VectorElement('Lagrange', submesh_top.ufl_cell(), 2)      # Lagrange multiplier
P2 = VectorElement('Lagrange', mesh.ufl_cell(), 2)      # für die Verschiebung quadratisch
P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)      # für den Druck linear

# Create mixed space
V = MixedFunctionSpace(FunctionSpace(mesh, P2),
                       FunctionSpace(mesh, P1),
                       FunctionSpace(submesh_top, LM))

# Define functions
(u, p, l) = TrialFunctions(V)		# Verschiebung, Druck, Lagrangsch-Multiplyer
(v, r, e)  = TestFunctions(V)           
du = Function(V)        # Hilfsgröße zum splitten
u0, p0, l0 = du.sub(0), du.sub(1), du.sub(2)
.
.
.

# free-slip BC
bcu_B  = DirichletBC(V.sub_space(0).sub(1), Constant((0)), bottom)
bcu_L  = DirichletBC(V.sub_space(0).sub(0), Constant((0)), left)
bcp_R  = DirichletBC(V.sub_space(1), Constant((0)), right)

If you separate properly your lhs from your rhs and then call solve your code should work. Apparently the functions lhs and rhs don’t work with the mixed dimensional branch (or maybe yes @cdaversin?

1 Like

Hello,

I have implemented your suggesting (thank you for them!) in my code with the change of

LM = FiniteElement('Lagrange', submesh_top.ufl_cell(), 1)      # Lagrange multiplier

because I thought, when the displacement on the top is supposed to be constant, the divergenz is supposed to be zero, and then I would have: scalar*scalar

l = inner(sigma(u,p),nabla_grad(v))*dx - k*inner(nabla_grad((p-p0)),nabla_grad(r))*dx -\
        tr_epsilon(u)/dt *r *dx +\
        nabla_div(u)*e*dL 
    r = -inner(T1,v)*ds_test(1) - tr_epsilon(u0)/dt *r *dx

Is there an error in my lhs and rhs? Because, when I run the file I get the error:
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs ('v_1^0',)
But as far as I see it, there are only constants on my rhs and only unknowns on my lhs.

I am not sure if it helps, but this is my Code with the implemented changes that gives the error message I wrote in the last post:

from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div

# Elasticity parameters 

E= 100*10**6		# E in Pa
nu = 0.25 		# Querkontraktionszahl
k0= 1*10**(-14)
visk = 1 * 10**(-5)
k = k0/visk		# Permeabilität	in m^4/Ns
b = 1           # Biot Koefficent
M = 10**14      # Biot Modul
F = 2        # Kraft welche auf die Platte wirkt 
a = 1           # Länge des Plattenabschnitts
G = Constant(E/(2*(1 + nu))) # Schubmodul

T = 1           # final time
num_steps = 1000    # number of time steps
dt = T / num_steps # time step size

# Create mesh and define function space
nx0 = 0
ny0 = 0
nx1 = 1
ny1 = 1
mesh = RectangleMesh(Point(nx0, ny0), Point(nx1, ny1), 40,40)

# Mark boundary subdomians
class rechts(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0],1,tol)
class links(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0],0,tol)
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1],1,tol)
class unten(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1],0,tol)

marker = MeshFunction('size_t',mesh,mesh.topology().dim()-1 ,0)
TopBoundary().mark(marker,1)
submesh_top = MeshView.create(marker,1)

right = rechts()
left = links()
bottom = unten()

left.mark(marker,2)
right.mark(marker,3)
bottom.mark(marker,4)

# Function spaces
LM = FiniteElement('Lagrange', submesh_top.ufl_cell(), 1)      # Lagrange multiplier
#LM = VectorElement('Lagrange', submesh_top.ufl_cell(), 2)      # Lagrange multiplier
P2 = VectorElement('Lagrange', mesh.ufl_cell(), 2)      # für die Verschiebung quadratisch
P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)      # für den Druck linear

# Create mixed space
#TH = MixedElement([P2, P1, LM])
#V = FunctionSpace(mesh, TH)
V = MixedFunctionSpace(FunctionSpace(mesh, P2),
                       FunctionSpace(mesh, P1),
                       FunctionSpace(submesh_top, LM))

# Define functions
(u, p, l) = TrialFunctions(V)		# Verschiebung, Druck, Lagrangsch-Multiplyer
(v, r, e)  = TestFunctions(V)           
du = Function(V)        # Hilfsgröße zum splitten
u0, p0, l0 = du.sub(0), du.sub(1), du.sub(2)
#u0, p0, l0 = split(du) # AttributeError: 'Function' object has no attribute '_ufl_function_space'

# Define Constants
n   = FacetNormal(mesh)
k   = Constant(k)			# Permeabilität 
nu  = Constant(nu)			# Querkontraktionszahl
b = Constant(b)
M = Constant(M)
F = Constant(F)
a = Constant(a)
mu = Constant(E/(2*(1 + nu)))
lam = Constant(E*nu/( (1+nu)*(1-2*nu)))

# Fluidflux
q = 0   # -q_vektor*n= q

# Kinematics
d = u.geometric_dimension()
I = Identity(d)             # Identity tensor

# linearer Straintensor
def epsilon(u):
#    return 0.5*(nabla_grad(u) + nabla_grad(u).T) 
    return sym(nabla_grad(u))
# Spur des Spannungstensors
def tr_epsilon(u):
    return tr(epsilon(u))    

#def n_por(u0):      # Porositiy
#    return 0.2/(1+tr_epsilon(u0)) 
n_por = 0.1

# linear solid strain
def sigma_s(u):
    return lam*tr_epsilon(u)*Identity(d) + 2*mu*epsilon(u)

# strain from both solid and fluid
def sigma(u,p):
	return (1-n_por) * sigma_s(u) -n_por*(p-p0)*I

# Measures
ds_test = ds(domain=mesh, subdomain_data=marker)
dL = Measure('dx', domain=submesh_top)

# Strss-Vector
T1 = Expression(("0", "-F"),degree =2, F=F)

Imp =inner(sigma(u,p),nabla_grad(v))*dx -inner(T1,v)*ds_test(1)  
Fl = - k*inner(nabla_grad((p-p0)),nabla_grad(r))*dx -(tr_epsilon(u) - tr_epsilon(u0))/dt *r *dx
Lag =  inner(nabla_div(u),e)*dL #inner(l,v)*dL + inner(u,e)*dL
eqn = Imp + Fl + Lag

l = inner(sigma(u,p),nabla_grad(v))*dx - k*inner(nabla_grad((p-p0)),nabla_grad(r))*dx -\
    tr_epsilon(u)/dt *r *dx +\
    inner(nabla_div(u),e)*dL
r = -inner(T1,v)*ds_test(1) - tr_epsilon(u0)/dt *r *dx

# free-slip BC
bcu_B  = DirichletBC(V.sub_space(0).sub(1), Constant((0)), bottom)     
bcu_L  = DirichletBC(V.sub_space(0).sub(0), Constant((0)), left)      
bcp_R  = DirichletBC(V.sub_space(1), Constant((0)), right)          

bc_ges = [bcu_L, bcu_B, bcp_R]

#l = lhs(eqn)
#r = rhs(eqn)

# Solve variational problem 
# Step in time
t = 0
      
for n in range(num_steps):
    # Update current time
    t += dt
    w = Function(V)

    solve(l == r, w, bc_ges)    
    (u,p,l) = w.split(True)

    assign(du.sub(0), u)

Hello,

I can indeed reproduce the error with your code, and I see two issues in your implementation :

(1) The variational formulation corresponding to your left hand side is named l, as the trial function used for your Lagrange multiplier. You should use another variable name.

(2) If you decompose the term k*inner(nabla_grad((p-p0)),nabla_grad(r))*dx in your left hand side, you actually end up with a term -k*inner(nabla_grad((p0)),nabla_grad(r))*dx which doesn’t involve the trial function and should then go in the right hand side. This is probably what causes the ArityMismatch error. Note : same when using p - p0 in sigma(u,p)

Also, I notice that you do assign(du.sub(0), u) but du is never updated. And p0 doesn’t seem to be updated either.

Hope that helps

1 Like

Hello cdaversin,

Thank you for your answer.

You were right, I wasn’t paying attention with the names for my trial function and lhs. I decomposed my functions further and now it is running without any errors!!!
Although I still have an error, because my pressure is zero (10^(-12)) for all the time steps, which I don’t really understand.

I used the assign(du.sub(0), u) to updated my u^n to my u^n+1 at the end of the time step.

And I don’t update my p0 because I just wanted it as an inital Condition. I define the function p0 and in the equations with p I replaced p with p-p0, and set p0=0 after the first time step so that it is only an initial condition. (At least that was my Plan)

Here is my updated Code:

from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div

# Elasticity parameters 

E= 100*10**6		# E in Pa
nu = 0.25 		# Querkontraktionszahl
k0= 1*10**(-14)
visk = 1 * 10**(-5)
k = k0/visk		# Permeabilität	in m^4/Ns
b = 1           # Biot Koefficent
M = 10**14      # Biot Modul
F = 2        # Force on the plate in N
a = 1           # length of the plate
G = Constant(E/(2*(1 + nu))) # shear modulus
K = Constant((1+nu)/(3*(1-nu))*2*G*(1-nu)/(1-2*nu))
B_big = Constant(b*M/(K+pow(b,2)*M))
nu_u = Constant((3*nu+b*B_big*(1-2*nu))/(3-b*B_big*(1-2*nu)))

T = 1           # final time
num_steps = 1000    # number of time steps
dt = T / num_steps # time step size

# Create mesh and define function space
nx0 = 0
ny0 = 0
nx1 = 1
ny1 = 1
mesh = RectangleMesh(Point(nx0, ny0), Point(nx1, ny1), 20,20)

# Mark boundary subdomians
class rechts(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0],1,tol)
class links(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0],0,tol)
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1],1,tol)
class unten(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1],0,tol)

marker = MeshFunction('size_t',mesh,mesh.topology().dim()-1 ,0)
TopBoundary().mark(marker,1)
submesh_top = MeshView.create(marker,1)

right = rechts()
left = links()
bottom = unten()

left.mark(marker,2)
right.mark(marker,3)
bottom.mark(marker,4)

# Function spaces
LM = FiniteElement('Lagrange', submesh_top.ufl_cell(), 1)      # Lagrange multiplier
#LM = VectorElement('Lagrange', submesh_top.ufl_cell(), 2)      # Lagrange multiplier
P2 = VectorElement('Lagrange', mesh.ufl_cell(), 2)      # für die Verschiebung quadratisch
P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)      # für den Druck linear

# Create mixed space
#TH = MixedElement([P2, P1, LM])
#V = FunctionSpace(mesh, TH)
V = MixedFunctionSpace(FunctionSpace(mesh, P2),
                       FunctionSpace(mesh, P1),
                       FunctionSpace(submesh_top, LM))

# Define functions
(u, p, l) = TrialFunctions(V)		# Verschiebung, Druck, Lagrangsch-Multiplyer
(v, r, e)  = TestFunctions(V)           
du = Function(V)        # Hilfsgröße zum splitten
u0, p0, l0 = du.sub(0), du.sub(1), du.sub(2)

# Define Constants
n   = FacetNormal(mesh)
k   = Constant(k)			# Permeabilität 
nu  = Constant(nu)			# Querkontraktionszahl
b = Constant(b)
M = Constant(M)
F = Constant(F)
a = Constant(a)
mu = Constant(E/(2*(1 + nu)))
lam = Constant(E*nu/( (1+nu)*(1-2*nu)))

# Fluidflux
q = 0   # -q_vektor*n= q

# Kinematics
d = u.geometric_dimension()
I = Identity(d)             # Identity tensor

# linearer Straintensor
def epsilon(u):
#    return 0.5*(nabla_grad(u) + nabla_grad(u).T) 
    return sym(nabla_grad(u))
# Spur des Spannungstensors
def tr_epsilon(u):
    return tr(epsilon(u))    

#def n_por(u0):      # Porositiy
#    return 0.2/(1+tr_epsilon(u0)) 
n_por = 0.1

# linear solid strain
def sigma_s(u):
    return lam*tr_epsilon(u)*Identity(d) + 2*mu*epsilon(u)

# strain from both solid and fluid
def sigma(u,p):
	return (1-n_por) * sigma_s(u) -n_por*(p)*I      # Anteil: n_por*(-p0)*I wird bei der rhs hinzugefügt

# Measures
ds_test = ds(domain=mesh, subdomain_data=marker)
dL = Measure('dx', domain=submesh_top)

# Stress-Vector
T1 = Expression(("0", "-F"),degree =2, F=F)

#Imp =inner(sigma(u,p),nabla_grad(v))*dx -inner(T1,v)*ds_test(1)  
#Fl = - k*inner(nabla_grad((p-p0)),nabla_grad(r))*dx -(tr_epsilon(u) - tr_epsilon(u0))/dt *r *dx
#Lag =  inner(nabla_div(u),e)*dL #inner(l,v)*dL + inner(u,e)*dL
#eqn = Imp + Fl + Lag

L_eqn1 = inner(sigma(u,p),nabla_grad(v))*dx                          # momentum balance
R_eqn1 = inner(T1,v)*ds_test(1) - inner(n_por*p0*I,nabla_grad(v))*dx # momentum balance 
L_eqn2 = - k*inner(nabla_grad((p)),nabla_grad(r))*dx - tr_epsilon(u)/dt *r *dx # Fluidmass balance with Darcy
R_eqn2 = - tr_epsilon(u0)/dt *r *dx - k*inner(nabla_grad((p0)),nabla_grad(r))*dx  
L_eqn3 = inner(nabla_div(u),e)*dL     # Lagrangian multiplier 

L_eqn = L_eqn1 + L_eqn2 + L_eqn3
R_eqn = R_eqn1 + R_eqn2

# free-slip BC
bcu_B  = DirichletBC(V.sub_space(0).sub(1), Constant((0)), bottom)     
bcu_L  = DirichletBC(V.sub_space(0).sub(0), Constant((0)), left)      
bcp_R  = DirichletBC(V.sub_space(1), Constant((0)), right)          

bc_ges = [bcu_L, bcu_B, bcp_R]

# Initial Conditions:
e_p0 = Expression('1/(3*a) *B_big*(1+nu_u)*F', degree =2, F=F, a=a, B_big=B_big,nu_u=nu_u)
e_u0 = Expression(('(F*nu_u/(2*G))*x[0]/a' , '-F*(1-nu_u)/(2*G)*(x[1]/a)'), degree =2, F=F, a=a, G=G ,nu_u=nu_u)

p_0 = interpolate(e_p0, V.sub_space(1))
u_0 = interpolate(e_u0, V.sub_space(0))

# Solve variational problem 
# Step in time
t = 0

vtkfile_p = File('Ergebnisse/p_Mandel_v10.pvd')
vtkfile_u= File('Ergebnisse/u_Mandel_v10.pvd')
      
for n in range(num_steps):
    if n == 0:
        # Assign initial Conditions to the functions
        assign(p0, p_0)
        assign(u0, u_0)
    else:
        p0 = 0

    # Update current time
    t += dt
    w = Function(V)

    solve(L_eqn == R_eqn, w, bc_ges)    
#    (u,p,l) = w.split(True)
    u,p,l = w.sub(0), w.sub(1), w.sub(2)
    
    u.rename('u','u')       # that Paraview doesnt view each new time step as a  
    p.rename('p','p')       # new variable
    
    vtkfile_u << (u, t)
    vtkfile_p << (p, t)

    assign(du.sub(0), u)  

`