Add in situ stress in linear elasticity model

General
1

I need to add a in situ stress condition in mi sigma function wich is depth dependent like:

t=np.array([[22.26,4.29],[4.29,13.8],[0,0,0])
T=as_tensor(t)
y_insitu=849.9 #depth stress in situ medition
y_max=1179.735 #upper limit

#Dpth depend function 

stress_in = Expression(("(y_max-x[1])/(y_max-y_insitu)",0),degree=1,y_max=y_max,y_insitu=y_insitu,

#Functions

def eps(v):
  return sym(grad(v))
def sigma(v):
  return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) + T*stress_in

The error said: "ufl.log.UFLExceptio: cant add expressions with diferen shapes.
How can i fixed it?

The complete code is:

from dolfin import*
from mshr import*
from fenics import*
from ufl import nabla_div
import numpy as np

#Malla

from dolfin import *
mesh = Mesh("fantasma1.xml")
##plot(mesh)
#cd=MeshFunction('size_t',mesh,"Brillador_physical_region.xml")
#plot(cd)

#Inputs
E = Constant(46.9e3)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda =E*nu/(1+nu)/(1-2*nu)
x_max = 1669.4431
y_r=1078.926
t=np.array([[22.26,4.29],[4.29,13.8],[0,0,0])
T=as_tensor(t)
y_insitu=849.9 #depth stress in situ medition
y_max=1179.735 #upper limit

#Dpth depend function 

stress_in = Expression(("(y_max-x[1])/(y_max-y_insitu)",0),degree=1,y_max=y_max,y_insitu=y_insitu,

#Functions

def eps(v):
  return sym(grad(v))
def sigma(v):
  return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) + T*stress_in
rho_g = 0.027
f = Constant((0, -rho_g))
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
#Dirichlet BC

def bottom(x, on_boundary):
  return near(x[1], 0.,tol)
def left(x, on_boundary):
  return near(x[0], 0.,tol)
def right(x, on_boundary):
  return near(x[0], x_max,tol)  
bc3 = DirichletBC(V.sub(0), Constant(0.), right)
bc2 = DirichletBC(V.sub(0), Constant(0.), le

ft)
bc1 = DirichletBC(V.sub(1), Constant(0.), bottom)
bc=[bc1,bc2,bc3]

#Weak formulation

a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx #+inner(rigth_function, u_)*ds(2)
u = Function(V, name="Displacement")
solve(a == l, u, bc)

#Principal Stress

Vsig = TensorFunctionSpace(mesh, "DG",degree=1)
sig = Function(Vsig, name="Cauchy_Stress")
sig.assign(project(-sigma(u), Vsig)); 

sx = sig.sub(0)
sy = sig.sub(3)
txy = sig.sub(1)

def sigma1(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
    return((sigma_x+sigma_y)/2 + sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2))

def sigma3(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
    return((sigma_x+sigma_y)/2 - sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2));

S1 = sigma1(sx,sy,txy);
S3 = sigma3(sx,sy,txy)

Vp = TensorFunctionSpace(mesh, "DG",degree=1)
sp = Function(Vsig, name="Principal_Stress")
sp.assign(project(as_tensor(((S1,0.),(0.,S3))), Vp));

file_results = XDMFFile("fantasma1_Prueba5.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
file_results.write(sig, 0.)
file_results.write(sp,0.)