In the following code, with an elasticity equation, i am trying to impose a zero normal displacement on a dirichlet boundary using lagrange multiplier. Can someone let me know what i should fix, probably in the mixed bilinear form (a_N in my code) to achieve this?
from dolfin import *
# dimensions and material parameters
L = 1; W = 0.2
mu = 1; lambda_ = 1.25
# mesh generation
nx = 5; ny = 2; nz = 2
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(L, W, W), nx, ny, nz)
class Dirichlet(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Neumann(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)
# markers for lagrange multiplier's mesh
facets = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Dirichlet().mark(facets, 1)
# Domain markers
domain = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
# two meshes
mesh_dmin = MeshView.create(domain, 0)
mesh_bdry = MeshView.create(facets, 1)
Velem = VectorElement("CG", mesh_dmin.ufl_cell(), 1, dim=3)
U = FunctionSpace(mesh_dmin, Velem)
Qelem = FiniteElement("CG", mesh_bdry.ufl_cell(), 1)
Q = FunctionSpace(mesh_bdry, Qelem)
# mixed space
W = MixedFunctionSpace(U, Q)
# markers for Dirichlet and Neumann BCs
facet_bdry = MeshFunction("size_t", mesh_dmin, mesh_dmin.geometry().dim()-1, 0)
Dirichlet().mark(facet_bdry, 1)
Neumann().mark(facet_bdry, 2)
# trial and test functions
(u, q) = TrialFunctions(W)
(u_t, q_t) = TestFunctions(W)
def epsilon(u):
return sym(grad(u))
def sigma(u):
return 2 * mu * epsilon(u) + lambda_ * tr(epsilon(u)) * Identity(len(u))
# measures
dx_ = Measure("dx", domain=W.sub_space(0).mesh())
ds_d = Measure("ds", domain=W.sub_space(0).mesh(), subdomain_data=facet_bdry, subdomain_id=1)
ds_n = Measure("ds", domain=W.sub_space(0).mesh(), subdomain_data=facet_bdry, subdomain_id=2)
# measure for LM mesh
dx_i = Measure("dx", domain=W.sub_space(1).mesh())
# first bilinear form w/o LM
a_M = inner(sigma(u), epsilon(u_t)) * dx_
g_L = Constant((1e0, 0, 0))
L_M = inner(g_L, u_t) * ds_n
w_ = Function(W)
(u_, q_) = w_.split()
## to solve elasticity problem with LM for imposing zero normal displacement on zero boundary
# dirichlet BCs
bc_Ly = DirichletBC(W.sub_space(0).sub(1), Constant(0), facet_bdry, 1)
bc_Lz = DirichletBC(W.sub_space(0).sub(2), Constant(0), facet_bdry, 1)
bcs = [bc_Ly, bc_Lz]
# mixed bilinear form --> !NOT WORKING!
a_N = q_t * inner(u, FacetNormal(mesh_bdry)) * dx_i + q * inner(u_t, FacetNormal(mesh_bdry)) * dx_i
a = a_M + a_N
problem = LinearVariationalProblem(a, L_M, w_, bcs)
nsolver = LinearVariationalSolver(problem)
nsolver.solve()