G[i] is the shear modulus of the subdomain i.

```
from __future__ import print_function
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
Nx = 10
Ny = 10
Lx = Nx
Ly = Ny
mesh = RectangleMesh(Point(0., 0.), Point(Lx, Ly), Nx, Ny, "crossed")
class Inclusion(SubDomain):
def inside(self, x, on_boundary):
return (between(x[1], (2, 6)) and between(x[0], (2, 6)))
inclusion = Inclusion()
subdomains = MeshFunction("size_t", mesh, 2)
subdomains.set_all(0)
inclusion.mark(subdomains, 1)
nu = Constant(0.3)
G = [Constant(1), Constant(10)]
def eps(v):
return sym(grad(v))
def sigma(v,i):
lmbda = 2*G[i]*nu/(1-2*nu)
return lmbda*tr(eps(v))*Identity(2) + 2.0*G[i]*eps(v)
f = Constant((0,0))
V = VectorFunctionSpace(mesh, 'Lagrange', 1)
du = TrialFunction(V)
u_ = TestFunction(V)
dx = Measure("dx")(subdomain_data=subdomains)
a = inner(sigma(du,0), eps(u_))*dx(0) + inner(sigma(du,1), eps(u_))*dx(1)
l = inner(f, u_)*dx(0) + inner(f, u_)*dx(1)
def bottom(x, on_boundary): return near(x[1],0.)
def top(x, on_boundary): return near(x[1],Ly)
def left(x, on_boundary): return near(x[0],0.)
def right(x, on_boundary): return near(x[0],Lx)
dirichelt_BCs = [DirichletBC(V, Constant((0.,0.)), bottom),
DirichletBC(V, Constant((1.,0.)), top),
DirichletBC(V.sub(1), Constant(0.), left),
DirichletBC(V.sub(1), Constant(0.), right)]
u = Function(V, name="Displacement")
solve(a == l, u, dirichelt_BCs)
T = TensorFunctionSpace(mesh, "Lagrange", 1)
stress = project(sigma(u), T)
[s11, s12, s21, s22] = stress.split(True)
```