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

variational formulation
1

This is continuation of the problem described here with an additional constraint (tagged as New Constraint below). The problem is as follows:

- \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(1,0) = 0
u_x(0,1)=0

along with the constraints
(a) 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)
and
(b)u_y(x,1) - u_y(x,0)=u_y(1,1) - u_y(0,0)
\underbrace{\Rightarrow u_y(x,1) - u_y(x,0)=u_y(1,1) }_{\textrm{New Constraint}}

where \mathbf{u}=(u_x,u_y) is the displacement vector. Please note that unlike u_x(1,1) = 0.05 in Constraint (a), u_y(1,1) is unknown in Constraint (b)

I am attempting this problem by defining four scalar Lagrange multipliers on the right face \Gamma_R=(1,y) , left face \Gamma_L=(0,y), bottom face \Gamma_B=(x,0) and top face \Gamma_T=(x,1). While the code runs, the solution does not seem to be correct. 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 

#  Dirac Delta
class Delta(UserExpression):
    def __init__(self, eps, x0, **kwargs):
        self.eps = eps
        self.x0 = x0
        UserExpression.__init__(self, **kwargs) 
    def eval(self, values, x):
        eps = self.eps
        values[1] = eps/(np.linalg.norm(x-self.x0)**2 + eps**2)
        values[0] = 0.0

    def value_shape(self): return (2, )

# Origin Point x=(0,0)
class Originpt(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0.0) and near(x[1], 0.0)

# Bottom Right point x=(1,0)
class BRpt(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],1.0) and near(x[1], 0.0)

# Left Top point x=(0,1)
class LTpt(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0.0) and near(x[1], 1.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(40, 40)

origin = Originpt()
cp = CPpt()                                                                 # Origin point
br = BRpt()                                                                 # Bottom Right point
lt = LTpt()                                                                   # Left Top point

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)

# file = File("subdomains.pvd")
# file << facets
###########################################
# Define Function Spaces

Wu = VectorElement("CG", mesh.ufl_cell(), 2)          # for u
Wlam_R = FiniteElement("R", mesh_R.ufl_cell(), 0)     # for lam on right face
Wlam_L = FiniteElement("R", mesh_L.ufl_cell(), 0)     # for lam on left face
Wlam_T = FiniteElement("R", mesh_T.ufl_cell(), 0)     # for lam on top face
Wlam_B = FiniteElement("R", mesh_B.ufl_cell(), 0)     # for lam on bottom face

V_u = FunctionSpace(mesh,Wu) 
V_lam_R = FunctionSpace(mesh_R, Wlam_R)
V_lam_L = FunctionSpace(mesh_L, Wlam_L)
V_lam_T = FunctionSpace(mesh_T, Wlam_T)
V_lam_B = FunctionSpace(mesh_B, Wlam_B)

V = MixedFunctionSpace(V_u, V_lam_L, V_lam_R, V_lam_B, V_lam_T)

u, lam_L, lam_R, lam_B, lam_T= TrialFunctions(V)
uh, lam_Lh, lam_Rh, lam_Bh, lam_Th = TestFunctions(V)

# Measures
ds_L = Measure("dx", domain = V.sub_space(1).mesh())
ds_R = Measure("dx", domain = V.sub_space(2).mesh())
ds_B = Measure("dx", domain = V.sub_space(3).mesh())
ds_T = Measure("dx", domain = V.sub_space(4).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_br = DirichletBC(V.sub_space(0).sub(1), Constant(0.0), br, method="pointwise")
bc_lt = DirichletBC(V.sub_space(0).sub(0), Constant(0.0), lt, method="pointwise")

bcs=[bcOrigin, bc_br, bc_lt, bc_cp]
##############################################
# Problem Formulation

delta = Delta(eps=1E-14, x0=np.array([1.0, 1.0]), degree=5)
u_cp = inner(u, delta)                  # u_cp =u_y(1,1)
uh_cp = inner(uh, delta)              # uh_cp =uh_y(1,1)

a = inner(sigma(u), eps(uh))*dx

# Constraint (a)
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 

# Constraint (b)
a += inner(lam_T, uh[1])* ds_T + inner(lam_Th, u[1])* ds_T \
    - inner(lam_B, uh[1])* ds_B - inner(lam_Bh, u[1])* ds_B \
        - inner(lam_T, uh_cp)* ds_T - inner(lam_Th, u_cp)* ds_T

L =  inner(lam_Rh, stretch)* ds_R

usol = Function(V)

solve(a == L, usol, bcs, solver_parameters={"linear_solver": "lu"})
disp, lamL, lamR, lamB, lamT = usol.sub(0), usol.sub(1), usol.sub(2), usol.sub(3), usol.sub(4)

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)
2

Dear @Kamalendu_Ghosh,

I am currently working on creating a framework for multi point constraints indolfinx. I can let you know when this is made public.

The approach I will be using is elimination of degrees of freedom with a master-slave concept, thus, for each degree of freedom, one would have an equation such as:
u_x(1,y_i) = u_x(0, y_i) - u_x(1,1)
and
u_y(x_i,1) = u_y(x_i,0) - u_y(1,1)
for each 1,\dots, N-1 for a UnitSquareMesh(N, N).

I have not tried to use the approach of Lagrange multipliers for this, so I can only give you some pointer:

  1. Start with a simpler problem (First do pure periodic conditions with Lagrange multipliers, verify that it works as it should)
  2. Try to have an equality constraint u_y(x_i,1) = u_y(x_i,0) and see if you can get this to work.

If you get 1 and 2 to work, it should be possible to solve your full problem.

2 Likes
3

Hi @dokken,

Thanks for your reply. It would be great when you make the MPC using the master-slave concept public. That is precisely what I am looking for.

In the meantime, I will try to implement the problem using Lagrange multipliers. Thanks for your suggestions regarding the same.

4

Hi again,
My repo is now available at: https://github.com/jorgensd/dolfinx_mpc

1 Like
5

Thanks @dokken for sharing it.