Hello everyone,
i am running follwing code.
from dolfin import *
import matplotlib.pyplot as plt
import os
import dolfin as dlf
import numpy as np
import math
import ufl as uf
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
# traction
class Traction(UserExpression):
def __init__(self):
super().__init__(self)
self.t = 0.0
def eval(self, values, x):
values[0] = 0*self.t
values[1] = 0.0
values[2] = 0.0
def value_shape(self):
return (3,)
# Kinematics
def pk1Stress(u,mu,nu):
I = Identity(V.mesh().geometry().dim()) # Identity tensor
F = I + grad(u)
J = det(F) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
I1 = tr(C) # Invariants of deformation tensors
C2 = C*C
IS1 = tr(C2)
I2= (1/2)*(I1**(2)-IS1)
DI2C= I1*I-C
DJDC =(J/2)*inv(C)
pk2 = mu*((1/J**(2))*DI2C-(2*I2/J**(3))*DJDC+2*J**((4*nu-1)/(1-2*nu))*DJDC) # second PK stress
return pk2
def geometry_3d():
mesh = UnitCubeMesh(6, 6, 6)
boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
x0 = AutoSubDomain(lambda x: near(x[0], 0))
x1 = AutoSubDomain(lambda x: near(x[0], 1))
y0 = AutoSubDomain(lambda x: near(x[1], 0))
z0 = AutoSubDomain(lambda x: near(x[2], 0))
x0.mark(boundary_parts, 1)
y0.mark(boundary_parts, 2)
z0.mark(boundary_parts, 3)
x1.mark(boundary_parts, 4)
return boundary_parts
# Create mesh and define function space ============================================
facet_function = geometry_3d()
mesh = facet_function.mesh()
gdim = mesh.geometry().dim()
dx = Measure("dx")
ds = Measure("ds", subdomain_data=facet_function, subdomain_id=4)
print('Number of nodes: ',mesh.num_vertices())
print('Number of cells: ',mesh.num_cells())
# Limit quadrature degree
dx = dx(degree=4)
ds = ds(degree=4)
# Create function space
element_v = VectorElement("P", mesh.ufl_cell(), 1) #
V = FunctionSpace(mesh, element_v )
#V = VectorFunctionSpace(mesh, "Lagrange", 1)
# Define test and trial functions
du = TrialFunction(V)
_u = TestFunctions(V)
u = Function(V)
# Create tensor function spaces for stress and stretch output
W_DFnStress = TensorFunctionSpace(mesh, "DG", degree=0)
defGrad = Function(W_DFnStress, name='F')
PK1_stress = Function(W_DFnStress, name='PK1')
# Displacement from previous iteration
b = Constant((0.0,0.0, 0.0)) # Body force per unit mass
h = Traction() # Traction force on the boundary
# Define Dirichlet boundary
bc0 = DirichletBC(V.sub(0), Constant(0.), facet_function, 1)
bc1 = DirichletBC(V.sub(1), Constant(0.), facet_function, 2)
bc2 = DirichletBC(V.sub(2), Constant(0.), facet_function, 3)
tDirBC = Expression(('1.0*time_'),time_=0.0 , degree=0)
bc3 = DirichletBC(V.sub(0), tDirBC, facet_function, 4)
bcs = [bc0,bc1,bc2,bc3]
# material parameters
mu = 1.1712
nu = 0.25
# # Total potential energy
pkstrs = pk1Stress(u,mu,nu)
I = Identity(V.mesh().geometry().dim())
dgF = I + grad(u)
F = inner(dot(dgF,pkstrs),grad(_u))*dx - dot(b, _u)*dx - dot(h, _u)*ds
J = derivative(F , u, du)
# Create nonlinear variational problem and solve
problem = NonlinearVariationalProblem(F, u, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.parameters['newton_solver']['linear_solver'] = 'mumps'
# Time stepping parameters
dt = 0.1
t, T = 0.0, 10*dt
# Save solution in VTK format
file_results = XDMFFile("./Blatz/TestUniaxialLoading/Uniaxial.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
stretch_vec = []
stress_vec=[]
stress_ana=[]
while t <= T:
print('time: ', t)
# Increase traction
h.t = t
tDirBC.time_ = t
# solve and save disp
solver.solve()
# Extract solution components
u_plot = u
u_plot.rename("u", "displacement")
# get stretch at a point for plotting
point = (0.5,0.5,0.5)
DF = I + grad(u_plot)
defGrad.assign(project(DF, W_DFnStress))
stretch_vec.append(defGrad(point)[0])
# get stress at a point for plotting
PK1_s= pk1Stress(u_plot,mu,nu)
PK1_stress.assign(project(PK1_s, W_DFnStress))
stress_vec.append(PK1_stress(point)[0])
# save xdmf file
file_results.write(u_plot, t)
file_results.write(defGrad, t)
file_results.write(PK1_stress, t)
# time increment
t += float(dt)
# get analytical solution
stretch_vec = np.array(stretch_vec)
stress_vec = np.array(stress_vec)
# plot results
f = plt.figure(figsize=(12,6))
plt.plot(stretch_vec, stress_vec,'r-')
plt.show()
I am getting this error. can someone help me with error?
Number of nodes: 343
Number of cells: 1296
Traceback (most recent call last):
File "/home/willkin/Desktop/APS3_Fenics/blatz.py", line 110, in <module>
F = inner(dot(dgF,pkstrs),grad(_u))*dx - dot(b, _u)*dx - dot(h, _u)*ds
File "/home/willkin/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/operators.py", line 380, in grad
f = as_ufl(f)
File "/home/willkin/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/constantvalue.py", line 471, in as_ufl
raise UFLValueError("Invalid type conversion: %s can not be converted"
ufl.log.UFLValueError: Invalid type conversion: (Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 0, None), MultiIndex((FixedIndex(0),))), Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 0, None), MultiIndex((FixedIndex(1),))), Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 0, None), MultiIndex((FixedIndex(2),)))) can not be converted to any UFL type.
Thanks in advance.