Sorry Dokken,
As you have linked me to how to overcome my problem, I have developed the code below
from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
#***********************************************************#
# Create the mesh
Nx, Ny = 4 , 4
L, H = 1, 1
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, 'left')
#************************************************************#
E0 = Constant(2.1e11)
nu0 = Constant(0.3)
beta0 = E0/(1-(nu0*nu0))
C0 = as_matrix([[beta0, beta0*nu0, 0],[beta0*nu0, beta0, 0],[0, 0, beta0*(1-nu0)/2]])
#*************************************************************#
E1 = Constant(2.1e9)
nu1 = Constant(0.3)
beta1 = E1/(1-(nu1*nu1))
C1 = as_matrix([[beta1, beta1*nu1, 0],[beta1*nu1, beta1, 0],[0, 0, beta1*(1-nu1)/2]])
#*************************************************************#
def eps(v):
return sym(grad(v))
def strain2voigt(e):
"""e is a 2nd-order tensor, returns its Voigt vectorial representation"""
return as_vector([e[0,0],e[1,1],2*e[0,1]])
def voigt2stress(s):
"""
s is a stress-like vector (no 2 factor on last component)
returns its tensorial representation
"""
return as_tensor([[s[0],s[2]],[s[2],s[1]]])
def sigma0(v):
return voigt2stress(dot(C0,strain2voigt(eps(v))))
#*************************************************************#
def sigma1(v):
return voigt2stress(dot(C1,strain2voigt(eps(v))))
#*************************************************************#
# Create classes for defining parts of the boundaries
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
#*************************************************************#
# Create classes for defining parts of the interior of the domain
tol = 1E-14
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.5 + tol
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.5 - tol
#*************************************************************#
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
Omega_0 = Omega_0()
Omega_1 = Omega_1()
#*************************************************************#
# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
Omega_1.mark(domains, 1)
#*************************************************************#
# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
#*************************************************************#
# Define function space
V = VectorFunctionSpace(mesh, 'CG', 1)
#*************************************************************#
# Define variational problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = Function(V, name = 'Displacement')
#*************************************************************#
# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
#*************************************************************#
# Define Dirichlet boundary conditions at bottom boundary
bcs = DirichletBC(V, Constant((0., 0.)), boundaries, 4)
#*************************************************************#
a = inner(sigma0(du), eps(u_))*dx(0) + inner(sigma1(du), eps(u_))*dx(1)
#*************************************************************#
# uniform traction on top boundary
T = Constant((0, -1e7))
l = dot(T, u_)*ds(2)
#*************************************************************#
solve(a == l, u, bcs)
and it is working successfully based on what I have seen here:
