Multi Point Constraint using Mixed Dimensional Branch (1/2)

variational formulation
1

Hello All,

I am trying to solve a constrained problem in linear elasticity in the domain \Omega=(0,1)^2 \subset \mathbb{R}^2. The problem described below:

- \textrm{div} ( \mathbf{\sigma}(\mathbf{u}))=\mathbf{0} in \Omega
subject to the boundary conditions:

u_x (0,0)= 0
u_x (1,1)= 0.05
u_y = 0 in (x,0)

along with the constraint
u_x(1,y) - u_x(0,y)=u_x(1,1) - u_x(0,0)
\Rightarrow u_x(1,y) - u_x(0,y)=u_x(1,1) ,

where \mathbf{u}=(u_x,u_y) is the displacement vector.

I am attempting this problem by defining two scalar Lagrange multipliers on the right face \Gamma_R=(1,y) and left face \Gamma_L=(0,y). But I am getting an “Assertion error”. Please help.

Following is the MWE:

from dolfin import *
## Definitions
# strain
def eps(v):
    return sym(grad(v))
# stress
def sigma(v):
    lmbda = E*nu/(1+nu)/(1-2*nu)
    mu = E/2/(1+nu)
    return lmbda*tr(eps(v)) * Identity(2) + 2*mu*(eps(v))

## Classes 
# Origin Point
class Originpt(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0.0) and near(x[1], 0.0)
# control point x=(1,1)
class CPpt(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],1.0) and near(x[1], 1.0)
######################################################
# Material Properties 
E = 50e3
nu = 0.25 
######################################################
#  Mesh and subdomains

mesh = UnitSquareMesh(20, 20)

origin = Originpt()
cp = CPpt()                                         

bottom_fc = CompiledSubDomain("near(x[1],0.0) && on_boundary")            # bottom face
top_fc = CompiledSubDomain("near(x[1],1.0) && on_boundary")                 # top face
right_fc = CompiledSubDomain("near(x[0],1.0) && on_boundary")              # right face
left_fc = CompiledSubDomain("near(x[0],0.0) && on_boundary")                # left face        

facets = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

facets.set_all(0)
bottom_fc.mark(facets, 1)
top_fc.mark(facets, 2)
left_fc.mark(facets, 3)
right_fc.mark(facets, 4)

mesh_B=MeshView.create(facets,1)
mesh_T=MeshView.create(facets,2)
mesh_L=MeshView.create(facets,3)
mesh_R=MeshView.create(facets,4)

###################################################
# Define Function Spaces

Wu = VectorElement("CG", mesh.ufl_cell(), 2)              # for u
Wlam_R = FiniteElement("R", mesh_R.ufl_cell(), 0)     # for Lagrange Multiplier on Right face
Wlam_L = FiniteElement("R", mesh_L.ufl_cell(), 0)     # for Lagrange Multiplier on Left face
V_u = FunctionSpace(mesh,Wu) 
V_lam_R = FunctionSpace(mesh_R, Wlam_R)
V_lam_L = FunctionSpace(mesh_L, Wlam_L)

V = MixedFunctionSpace(V_u, V_lam_L, V_lam_R)

u, lam_L, lam_R = TrialFunctions(V)
uh, lam_Lh, lam_Rh = TestFunctions(V)

# Measures
ds_L = Measure("dx", domain = V.sub_space(1).mesh())
ds_R = Measure("dx", domain = V.sub_space(2).mesh())

###################################################
#Define Dirichlet Boundary
stretch = Constant(0.05)

bcOrigin = DirichletBC(V.sub_space(0), Constant((0.0, 0.0)), origin, method="pointwise")
bc_cp = DirichletBC(V.sub_space(0).sub(0), stretch, cp, method="pointwise")
bc_bottom = DirichletBC(V.sub_space(0).sub(1), Constant(0.0), bottom_fc)

bcs=[bcOrigin, bc_bottom, bc_cp]
###################################################
# Problem Formulation

a = inner(sigma(u), eps(uh))*dx
a += inner(lam_R, uh[0])* ds_R + inner(lam_Rh, u[0])* ds_R \
    - inner(lam_R, uh[0])* ds_L - inner(lam_Rh, u[0])* ds_L \
        - inner(lam_Rh, stretch)* ds_R

f = Constant((0, 0))
L = dot(f, uh)*dx 

usol = Function(V)
solve(a == L, usol, bcs)

Attached is the error message

Error Message
Error Message1018Ă—267 13.5 KB

2

Hello,

The failing assertion means that the size (number of blocks) in your bilinear form is not coherent with the size of your solution vector. I see two mistakes in your bilinear form :
(1) The term inner(lam_Rh, stretch)* ds_R doesn’t involve trial functions, so it should be in the linear form instead.
(2) The term - inner(lam_R, uh[0])* ds_L - inner(lam_Rh, u[0])* ds_L integrates lam_R and lam_Rh on V_lam_L (by definition of ds_L), but these are not defined on this domain. I guess this should be replaced by : - inner(lam_L, uh[0])* ds_L - inner(lam_Lh, u[0])* ds_L ?

Note : there is ongoing work on Multi Point Constraint methods by @dokken in Dolfin-x that could be of interest.

1 Like
3

@cdaversin Thanks for your reply. Both the points that you raised were valid. I have incorporated the changes and the code seems to work now.

Updated part:

a = inner(sigma(u), eps(uh))*dx
a += inner(lam_R, uh[0])* ds_R + inner(lam_Rh, u[0])* ds_R \
    - inner(lam_L, uh[0])* ds_L - inner(lam_Lh, u[0])* ds_L 

L = + inner(lam_Rh, stretch)* ds_R

usol = Function(V)
solve(a == L, usol, bcs, solver_parameters={"linear_solver": "lu"})

disp, lamL, lamR = usol.sub(0), usol.sub(1), usol.sub(2)
wfil = XDMFFile(comm,'Solution.xdmf')
wfil.parameters["flush_output"] = True
wfil.parameters["functions_share_mesh"] = True
wfil.parameters["rewrite_function_mesh"] = False
wfil.write(disp)
1 Like