Dear all,
Thinking in a more complex problem, I’m trying to solve the elasticity equation on two coupled domains \Omega_l=[0,1]\times[0,1] and \Omega_r=[1,2]\times[0,1] (see figure). The problem is as follows:
\begin{alignat}{2}
-\text{div}(\sigma(\boldsymbol{u}_l)) &= \boldsymbol{0} \quad &&\text{in } \Omega_l \\
-\text{div}(\sigma(\boldsymbol{u}_r)) &= \boldsymbol{0} \quad &&\text{in } \Omega_r
\end{alignat}
subject to boundary and coupling conditions
\begin{alignat}{2}
\boldsymbol{u_l} = \boldsymbol{u_r} &= \boldsymbol{0} \quad &&\text{in } \Gamma_D \\
\sigma(\boldsymbol{u_l})\boldsymbol{n}_l &= \boldsymbol{g} \quad &&\text{in } \Gamma_N \\
\boldsymbol{u_l} &= \boldsymbol{u_r} \quad &&\text{in } \Gamma_I \\
\sigma(\boldsymbol{u_l})\boldsymbol{n}_l &= -\sigma(\boldsymbol{u_r})\boldsymbol{n}_r \quad &&\text{in } \Gamma_I
\end{alignat}
where \boldsymbol{u}_l and \boldsymbol{u}_r are the displacement fields on \Omega_l and \Omega_r, \boldsymbol{n}_l and \boldsymbol{n}_r outward normal vectors on the boundaries of \Omega_l and \Omega_r, and \sigma the stress tensor.
I managed to impose the displacement continuity in \Gamma_I using Lagrange multipliers. However, I’m still not sure about the stress continuity. How can I manage the stress continuity in \Gamma_I?
The below code is a MWE of my attempt.
from dolfin import *
class Dirichlet(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 2) \
or near(x[1], 0) \
or near(x[1], 1)
class Neumann(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0)
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1)
# Mesh generation
n = 20
mesh = RectangleMesh(Point(0,0,0),Point(2,1,0),2*n,n)
# Domain markers
domain = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
for c in cells(mesh):
domain[c] = c.midpoint().x() > 1.0
# Interface markers
facets = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Interface().mark(facets,1)
# Meshes
mesh_L = MeshView.create(domain, 0)
mesh_R = MeshView.create(domain, 1)
mesh_I = MeshView.create(facets, 1)
# Function spaces
FE = VectorElement("CG", mesh_I.ufl_cell(), 1, dim=2) # Lagrange multiplier
VE = VectorElement("CG", mesh.ufl_cell(), 1, dim=2) # Displacement
V_L = FunctionSpace(mesh_L,VE) # Displacement on left domain
V_R = FunctionSpace(mesh_R,VE) # Displacement on right domain
V_I = FunctionSpace(mesh_I,FE) # Lagrange multiplier on interface
W = MixedFunctionSpace(V_L, V_R, V_I)
# Markers for Dirichlet and Neuman bcs
ff_L = MeshFunction("size_t", mesh_L, mesh_L.geometry().dim()-1, 0)
Dirichlet().mark(ff_L, 1)
Neumann().mark(ff_L, 2)
Interface().mark(ff_L, 3)
ff_R = MeshFunction("size_t", mesh_R, mesh_R.geometry().dim()-1, 0)
Dirichlet().mark(ff_R, 1)
Interface().mark(ff_R, 3)
# Boundary conditions
bc_L = DirichletBC(W.sub_space(0), Constant((0,0)), ff_L, 1)
bc_R = DirichletBC(W.sub_space(1), Constant((0,0)), ff_R, 1)
bcs = [bc_L, bc_R]
# Elasticity parameters
rho = Constant(2200)
lmbda = Constant(7.428e+04)
mu = Constant(7.428e+04)
# Neumann condition
g = Constant((1e+5, 0))
# Measures
dx_L = Measure("dx", domain=W.sub_space(0).mesh())
ds_L = Measure("ds", domain=W.sub_space(0).mesh(), subdomain_data=ff_L)
dx_R = Measure("dx", domain=W.sub_space(1).mesh())
ds_R = Measure("ds", domain=W.sub_space(1).mesh(), subdomain_data=ff_R)
ds_I = Measure("dx", domain=W.sub_space(2).mesh())
# Function and test functions
(ul, ur, l) = TrialFunctions(W)
(wl, wr, e) = TestFunctions(W)
# Stress tensor
def sigma(r):
return 2.0*mu*sym(grad(r)) + lmbda*tr(sym(grad(r)))*Identity(len(r))
# Variational problem
a_left = inner(sigma(ul),sym(grad(wl)))*dx_L \
+ inner(sigma(ul)*FacetNormal(mesh_L), wl)*ds_L(3)
L_left = inner(g, wl)*ds_L(2)
a_right = inner(sigma(ur),sym(grad(wr)))*dx_R \
- inner(sigma(ul)*FacetNormal(mesh_R), wr)*ds_R(3)
L_right = inner(Constant((0,0)), wr)*dx_R
a = a_left + a_right \
+ inner(l, wr)*ds_I + inner(ur-ul, e)*ds_I
L = L_left + L_right
# Output files
pvda_file = File("output/ul.pvd")
pvdb_file = File("output/ur.pvd")
# Solve
u_mixed = Function(W)
solve(a == L, u_mixed, bcs, solver_parameters={"linear_solver":"direct"})
ul, ur, uc = u_mixed.sub(0), u_mixed.sub(1), u_mixed.sub(2)
# Save solution
pvda_file << ul
pvdb_file << ur
PS: sorry for the number of edits. It took me a lot to figure out how to write the equations.
.