Hi. I am creating a linear elasticity model. I need to work whit anisotropy materials, then i use a tensor like:
from mshr import*
from fenics import*
from ufl import nabla_div
import numpy as np
from dolfin import *
mesh = Mesh("Brillador.xml")
#plot(mesh)
cd=MeshFunction('size_t',mesh,"Brillador_physical_region.xml")
#plot(cd)
alfa=30
a=alfa*np.pi/180
G=5
nu=0.3
R=2.25
Ex=30
Ey=Ex/R
a11 = np.sin(a)**4/Ey + np.cos(a)**4/Ex + np.sin(2*a)**2*0.25*(1/G-2*nu/Ey)
a12 = np.sin(2*a)**2*0.25*(Ey**-1+Ex**-1-G**-1)-nu*Ey**-1*(np.cos(a)**4+np.sin(a)**4)
a16 = np.sin(2*a)*((np.sin(a)**2/Ey-np.cos(a)**2/Ex)+(0.5*G**-1-nu*Ex**-1)*np.cos(2*a))
a22 = np.cos(a)**4/Ey+np.sin(a)**4/Ex+np.sin(2*a)**2/4*(G**-1-2*nu/Ey)
a26 = np.sin(2*a)*((np.cos(a)**2/Ey-np.sin(a)**2/Ex)-(0.5/G-nu/Ey)*np.cos(2*a))
a66 = np.sin(2*a)**2*(1/Ey+1/Ex+2*nu/Ey)+np.cos(2*a)**2/G
S = as_matrix([[a11,a12,a16],[a12,a22,a26],[a16,a26,a66]])
C = inv(S)
x_max = 969.3292
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 sigma(v):
return voigt2stress(dot(C, strain2voigt(eps(v))))
This example is and homogeneous model, but, i need to create a heterogeneous model in which y need to assign a diferent tensor to diferents physical volumes, is this possible? In others models, i have assign diferents elastic propierties to the pyshical groups like:
Er=30e3
Ev=15e3
class Young(UserExpression): # UserExpression instead of Expression
def __init__(self, markers, **kwargs):
super().__init__(**kwargs) # This part is new!
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 1:
values[0] = Er
elif self.markers[cell.index] ==2:
values[0] = Ev
else:
values[0] = Er
E=Young(cd,degree=0)
But i dont know how use it to assign another tensor like āSā as i show above.
I hope you can help me