Mixed-space formulation on distinct domains

Question

Hi,

I am trying to impose a Dirichlet boundary condition on apart $S$ of a domain's boundary using a Lagrange multiplier. Therefore I need to work with two unknown fields $(u,\lambda)$ ($u$ is a displacement field defined on the entire domain and $\lambda$ a Lagrange multiplier associated with the kinematic constraint defined on $S$), the energy for the variational problem is : $$\varepsilon(u) + \int_S \lambda \,(u-\overline{u}) \, \mathrm{d}S$$
where $\varepsilon(u)$ is an elastic energy (defined on the entire domain) and $\overline{u}$ the prescribed displacement.

I imagined a solution using a MixedFunctionSpace made of two FunctionSpaces : V defined over the entire mesh and V0 defined on a part of the boundary (corresponding to $S$) built with BoundaryMesh. Here is a 2-d test case :

mesh = UnitSquareMesh(10, 10)
TOL = 0.00001  

class FaceS(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0],1., TOL)

faceS = FaceS()

# V is the FunctionSpace for u
V = VectorFunctionSpace(mesh, "Lagrange", 1)

bmesh = BoundaryMesh(mesh,'exterior') 
S_bmesh = SubMesh(bmesh,faceS) 
# V0 is the FunctionSpace for \lambda
V0 = FunctionSpace(S_bmesh, "DG",0)

W = MixedFunctionSpace([V,V0]) 

It seems impossible to build W:

ufl.log.UFLException: Sub elements must live on the same domain (for now).

Has someone alternative ideas to define such a variational problem involving a Lagrange multiplier defined on a boundary ? Many thanks in advance !

Claire

Comments

I found this old question : https://fenicsproject.org/qa/3076/mixedfunctionspace-with-boundary-function-spaces?show=3076#q3076 strongly related to this topic..

So I have the impression that such a formulation cannot be solved within FEniCS yet.

Answer 1

Merging FunctionSpace that are defined on different meshes is not valid in FEniCS.

However, the integral can be set on different domains. Hence there might be some trick handling your issue.

You may define DirichletBC on a SubDomain having inner nodes, this might cancel out all the undesired DoFs of lamba.

Answer 2

I found a reasonably satisfying solution: define $\lambda$ on entire domain and impose its value on inner dofs by adding a term to the energy:

$$\epsilon(u) + \int_S \lambda (u-\overline{u}) \mathrm{d} S + \int_\Omega \delta \lambda^ 2\mathrm{d}V$$

where $\delta$ is a scalar field equal to $0.$ on $S$ :

mesh = UnitSquareMesh(10, 10)

TOL = 0.00001  
class FaceT(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0],1., TOL)

class FaceS(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0],0., TOL)

faceS = FaceS()
faceTd = FaceT()

V = VectorFunctionSpace(mesh, "Lagrange", 1)
V0 = FunctionSpace(S_bmesh, "DG",0)
W = MixedFunctionSpace([V,V0])

bc1 = DirichletBC(W.sub(0), (0.0,0.0), faceT)
bcs = [bc1]

subdomains = MeshFunction("size_t", mesh, 1)
subdomains.set_all(0)
faceS.mark(subdomains, 1)
my_ds = Measure("ds", domain=mesh, subdomain_data=subdomains)

Exp = Expression('x[0] > 0.1 ? 1. : 0.0001')
project_x = Constant((1.,0.))

w  = Function(W)
(u,p) = split(w)

d = u.geometric_dimension()  
FF = Identity(d) + grad(u)                  
C = FF.T*FF                        

Ic = tr(C)
JJ  = det(FF)

# Elasticity parameters
mu, Kpar, Jm = 1., 10., 100.0

ux = dot(project_x,u)

psi = mu/2.*(Ic-2 - 2*ln(JJ))+Kpar/2.*(ln(JJ))**2
E = psi*dx + inner(p,(ux+0.1))*my_ds(1)+ Exp*inner(p,p)*dx

v  = TestFunction(W)
du = TrialFunction(W) 

J = derivative(derivative(E,w,v),w,du)
F = derivative(E,w,v)  

problem = NonlinearVariationalProblem(F, w, bcs, J=J)
solver  = NonlinearVariationalSolver(problem)
solver.parameters.nonlinear_solver = "snes"
solver.parameters.snes_solver.linear_solver = "umfpack"