Hello, I am working on a problem in linearized elasticity in curvilinear coordinates.
When I implement the weak formulation of the problem, some troubles occur when I insert tensors in the definition of the weak formulation.
I enclose a MWE here below.
File CS1_fenics.txt
(('0.0','0.0','1.0/(Re+x[2])'),
('0.0','0.0','0.0'),
('1.0/(Re+x[2])','0.0','0.0'))
File CS2_fenics.txt
(('0.0','0.0','0.0'),
('0.0','0.0','0.0'),
('0.0','0.0','0.0'))
File CS3_fenics.txt
(('-Re-x[2]','0.0','0.0'),
('0.0','0.0','0.0'),
('0.0','0.0','0.0'))
File a_fenics.txt
(('1.0/pow(Re+x[2],4.0)*(lambda_+mu*2.0)','lambda_*1.0/pow(Re+x[2],2.0)','lambda_*1.0/pow(Re+x[2],2.0)','0.0','0.0','0.0'),
('lambda_*1.0/pow(Re+x[2],2.0)','lambda_+mu*2.0','lambda_','0.0','0.0','0.0'),
('lambda_*1.0/pow(Re+x[2],2.0)','lambda_','lambda_+mu*2.0','0.0','0.0','0.0'),
('0.0','0.0','0.0','mu*2.0','0.0','0.0'),
('0.0','0.0','0.0','0.0','mu*1.0/pow(Re+x[2],2.0)*2.0','0.0'),
('0.0','0.0','0.0','0.0','0.0','mu*1.0/pow(Re+x[2],2.0)*2.0'))
and finally, the Python script:
from ast import Expression
from fenics import *
import math
import sympy as sym
import os
## 3D box mesh -----------------------------------
epsilon = 0.01
H=2 # height of the cylinder
#-----------------------
# 3D Geometry
#-----------------------
mesh3 = BoxMesh(Point(0, 0, -epsilon), Point(math.pi, H, epsilon), 50, 50, 20)
# #---------------------
# # Function spaces
# # --------------------
V = VectorFunctionSpace(mesh3, "CG", degree=1,dim=3)
W=FunctionSpace(mesh3, "CG",1)
u_tr = TrialFunction(V)
u_test = TestFunction(V)
mu, Re, lambda_ = 1.0, 1.0, 1.0
# First Christoffel symbol
with open('./CS1_fenics.txt') as file:
expression = file.read()
CS1 = Expression(eval(expression), mu=mu, Re=Re, lambda_=lambda_, degree=3)
# Second Christoffel symbol
with open('./CS2_fenics.txt') as file:
expression = file.read()
CS2 = Expression(eval(expression), mu=mu, Re=Re, lambda_=lambda_, degree=3)
# Third Christoffel symbol
with open('./CS3_fenics.txt') as file:
expression = file.read()
CS3 = Expression(eval(expression), mu=mu, Re=Re, lambda_=lambda_, degree=3)
# Three-dimensional fourth order elasticity tensor (A)
with open('./a_fenics.txt') as a_fenics:
expression = a_fenics.read()
A = Expression(eval(expression), mu=mu, Re=Re, lambda_=lambda_, degree=3)
def epsV(u):
return Expression((
('Dx(u0,0)-CS1[0][0]*u0-CS2[0][0]*u1-CS3[0][0]*u2',
'Dx(u1,1)-CS1[1][1]*u0-CS2[1][1]*u1-CS3[1][1]*u1',
'Dx(u2,2)-CS1[2][2]*u0-CS2[2][2]*u1-CS3[2][2]*u2',
'Dx(u1,2)-2.*CS1[1][2]*u0-2.*CS2[1][2]*u1-2.*CS3[1][2]*u2+Dx(u2,1)',
'Dx(u0,2)-2.*CS1[0][2]*u0-2.*CS2[0][2]*u1-2.*CS3[0][2]*u2+Dx(u2,0)',
'Dx(u0,1)-2.*CS1[0][1]*u0-2.*CS2[0][1]*u1-2.*CS3[0][1]*u2+Dx(u1,0)')
),u0=u[0],u1=u[1],u2=u[2],CS1=CS1,CS2=CS2,CS3=CS3,degree=2)
# #----------------------
# # 3D applied body force
# #----------------------
f = Constant((0, 0, -1))
# vtkfile3 = File('./box_mesh.pvd')
# vtkfile3<<mesh3
# # --------------------
# # Weak form
# # --------------------
a = inner(epsV(u_tr), A*epsV(u_test))*dx
l = dot(f, u_test)*dx
When I run the code, it seems that Fenics cannot access the matrices of expressions CS1, CS2 and CS3.
Would anyone know how to solve this issue without having to declare the matrices entries one by one?
Apart from having a correct code, I would like to avoid having to expand all the sums in the matrix-vector products.
Thanks in advance.
but, for sure, the problem is linear.