# Problems encountered in solving thermoelastic problems

Hi all, recently I’m trying to learn some code to solve the thermoplastic problem with the following code:

``````import numpy as np
from arguments import args
from utils import walltime
import fenics as fe
import os
import glob

fe.parameters["form_compiler"]["representation"] = 'uflacs'
import warnings

@walltime
def simulation():
ambient_T = args.T_ambient
rho = args.rho
Cp = args.c_p
k = args.kappa_T
h = args.h_conv
eta = args.power_fraction
r = args.r_beam
P = args.power
EPS = 1e-8

E = fe.Constant(70e3)
nu = fe.Constant(0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
sig0 = fe.Constant(250.)
Et = E/100.  ###  切向模量
H = E*Et/(E-Et)  ###  塑性模量，又称硬化系数

x0 = 0.2*args.domain_x
y0 = 0.5*args.domain_y

total_t = 1200*1e-6
vel = 0.6*args.domain_x/total_t
dt = 1e-6
ts = np.arange(0., total_t + dt, dt)
print(f"total time steps = {len(ts)}")

ele_size = 0.01

mesh = fe.BoxMesh(fe.Point(0., 0., 0.), fe.Point(args.domain_x, args.domain_y, args.domain_z),
round(args.domain_x/ele_size), round(args.domain_y/ele_size), round(args.domain_z/ele_size))
mesh_file = fe.File(f'data/vtk/pbf/mesh/mesh.pvd')
mesh_file << mesh

# Define bottom surface
class Bottom(fe.SubDomain):
def inside(self, x, on_boundary):
# The condition for a point x to be on bottom side is that x[2] < EPS
return on_boundary and x[2] < EPS

# Define top surface
class Top(fe.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[2] > args.domain_z - EPS

# Define the other four surfaces
class SurroundingSurfaces(fe.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (x[0] < EPS or x[0] > args.domain_x - EPS or x[1] < EPS or x[1] > args.domain_y - EPS)

bottom = Bottom()
top = Top()
surrounding_surfaces = SurroundingSurfaces()
boundaries = fe.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
bottom.mark(boundaries, 1)
top.mark(boundaries, 2)
surrounding_surfaces.mark(boundaries, 3)
ds = fe.Measure('ds')(subdomain_data=boundaries)
normal = fe.FacetNormal(mesh)

deg_u = 2
deg_stress = 2
# TODO: can we do a tensor element?
U = fe.VectorFunctionSpace(mesh, "CG", deg_u)
W = fe.FunctionSpace(mesh, We)
# W = fe.FunctionSpace(mesh, We)
W0 = fe.FunctionSpace(mesh, W0e)

sig = fe.Function(W)
sig_old = fe.Function(W)
n_elas = fe.Function(W)
beta = fe.Function(W0)
# name = "a b" is bad (having a space)!
p = fe.Function(W0, name="Cumulative plastic strain")
u = fe.Function(U, name="Total displacement")
du = fe.Function(U, name="Iteration correction")
Du = fe.Function(U, name="Current increment")
v = fe.TrialFunction(U)
u_ = fe.TestFunction(U)

u_bcs = [fe.DirichletBC(U, fe.Constant((0., 0., 0.)), bottom)]

def eps(v):
return e

def sigma(eps_el):
return lmbda*fe.tr(eps_el)*fe.Identity(3) + 2*mu*eps_el

def as_3D_tensor(X):
return fe.as_tensor([[X[0], X[1], X[2]],
[X[3], X[4], X[5]],
[X[6], X[7], X[8]]])

def as_long_vector(X):
return fe.as_vector([X[0, 0], X[0, 1], X[0, 2], X[1, 0], X[1, 1], X[1, 2], X[2, 0], X[2, 1], X[2, 2]])

ppos = lambda x: (x+abs(x))/2.

def thermal_strain(dT):
alpha_V = 1e-5
return alpha_V*dT*fe.Identity(3)

def proj_sig(deps, dT, sig_old, p_old):
sig_n = as_3D_tensor(sig_old)   ### 将旧的应力场转为三维张量

# sig_elas = sig_n + sigma(deps)

d_eps_T = thermal_strain(dT)   ### 计算温度变化dT引起的热应变
sig_elas = sig_n + sigma(deps - d_eps_T)   ### 计算当前弹性应力张量

s = fe.dev(sig_elas)   ### 计算弹性应力张量的偏量部分
sig_eq = fe.sqrt(3/2.*fe.inner(s, s))   ### 计算等效应力
f_elas = sig_eq - sig0 - H*p_old   ###   计算弹性力，用于评估是否需要激活塑性本构。>0,应力状态超过屈服应力；
dp = ppos(f_elas)/(3*mu+H)   ###  计算塑性应变增量
n_elas = s/sig_eq*ppos(f_elas)/f_elas   ###   计算弹性流动方向张量
beta = 3*mu*dp/sig_eq   ###  计算硬化参数
new_sig = sig_elas-beta*s   ###   更新应力张量

return as_long_vector(new_sig), as_long_vector(n_elas), beta, dp

def sigma_tang(e):
N_elas = as_3D_tensor(n_elas)
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*fe.inner(N_elas, e)*N_elas-2*mu*beta*fe.dev(e)

# TODO
a_Newton = fe.inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -fe.inner(as_3D_tensor(sig), eps(u_))*dxm + 1e-10*fe.dot(normal, u_)*ds(2)

def local_project(v, V, u=None):
dv = fe.TrialFunction(V)
v_ = fe.TestFunction(V)
a_proj = fe.inner(dv, v_)*dxm
b_proj = fe.inner(v, v_)*dxm
solver = fe.LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = fe.Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)

def global_project(v, V, u=None):
'''
Not working.
'''
if u is None:
return fe.project(v, V)
else:
u.assign(fe.project(v, V))

P0 = fe.FunctionSpace(mesh, "DG", 0)
p_avg = fe.Function(P0, name="Plastic strain")

V = fe.FunctionSpace(mesh, 'CG', 1)
T_crt = fe.interpolate(fe.Constant(ambient_T), V)
T_pre = fe.interpolate(fe.Constant(ambient_T), V)
T_old = fe.interpolate(fe.Constant(ambient_T), V)
T_old.rename('T', 'T')
v = fe.TestFunction(V)

# If theta = 0., we recover implicit Eulear; if theta = 1., we recover explicit Euler; theta = 0.5 seems to be a good choice.
theta = 0.5
T_rhs = theta*T_pre + (1 - theta)*T_crt
T_bcs = [fe.DirichletBC(V, fe.Constant(ambient_T), bottom)]

class LaserExpression(fe.UserExpression):
def __init__(self, t):
super(LaserExpression, self).__init__()
self.t = t

def eval(self, values, x):
t = self.t
values[0] = 2*P*eta/(np.pi*r**2) * np.exp(-2*((x[0] - x0 - vel*t)**2 + (x[1] - y0)**2) / r**2)

def value_shape(self):
return ()

q_laser = LaserExpression(None)
q_convection = h * (T_rhs - ambient_T)

q_top = q_convection + q_laser
q_surr = q_convection

res_T = rho*Cp/dt*(T_crt - T_pre) * v * fe.dx + k * fe.dot(fe.grad(T_rhs), fe.grad(v)) * fe.dx \
- q_top * v * ds(2) - q_surr * v * ds(3)

files_vtk = glob.glob('data/vtk/pbf/sols' + f"/*")
files_xdmf = glob.glob('data/xdmf/pbf/' + f"/*")
for f in files_vtk + files_xdmf:
os.remove(f)

T_vtk_file = fe.File(f'data/vtk/pbf/sols/T.pvd')
T_vtk_file << T_old
u_vtk_file = fe.File(f'data/vtk/pbf/sols/u.pvd')
u_vtk_file << u
p_vtk_file = fe.File(f'data/vtk/pbf/sols/p.pvd')
p_vtk_file << p_avg

file_results = fe.XDMFFile('data/xdmf/pbf/u.xdmf')
file_results.parameters["functions_share_mesh"] = True
file_results.write(T_old, 0)
file_results.write(u, 0)
file_results.write(p_avg, 0)

plastic_inverval = 10

for i in range(len(ts) - 1):
# for i in range(100):

print(f"step {i + 1}, time = {ts[i + 1]}")
q_laser.t = theta*ts[i] + (1 - theta)*ts[i + 1]
solver_parameters = {'newton_solver': {'maximum_iterations': 20, 'linear_solver': 'mumps'}}
fe.solve(res_T == 0, T_crt, T_bcs, solver_parameters=solver_parameters)
T_pre.assign(T_crt)
print(f"min T = {np.min(np.array(T_pre.vector()))}")
print(f"max T = {np.max(np.array(T_pre.vector()))}")

if (i + 1) % plastic_inverval == 0:

T_crt_array = np.array(T_crt.vector())
T_crt_array = np.where(T_crt_array < args.T_ambient, args.T_ambient, T_crt_array)
T_crt_array = np.where(T_crt_array > args.T_melt, args.T_melt, T_crt_array)
T_old_array = np.array(T_old.vector())
T_old_array = np.where(T_old_array < args.T_ambient, args.T_ambient, T_old_array)
T_old_array = np.where(T_old_array > args.T_melt, args.T_melt, T_old_array)
dT = fe.Function(V)
dT.vector()[:] = T_crt_array - T_old_array

Nitermax, tol = 200, 1e-8  # parameters of the Newton-Raphson procedure
Nincr = 20
A, Res = fe.assemble_system(a_Newton, res, u_bcs)
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(fe.Constant((0., 0., 0.)))

niter = 0
nRes = 1.
tol = 1e-8
# while nRes/nRes0 > tol and niter < Nitermax:
# while niter < Nitermax:
while nRes > tol:
print("\n")
fe.solve(A, du.vector(), Res, "mumps")
Du.assign(Du+du)
deps = eps(Du)

print(f"du norm = {np.linalg.norm(np.array(du.vector()))}")

sig_, n_elas_, beta_, dp_ = proj_sig(deps, dT, sig_old, p)

local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
local_project(beta_, W0, beta)

print(f"sig dof = {np.array(sig.vector()).shape}")
print(f"sig norm = {np.linalg.norm(np.array(sig.vector()))}")
print(f"n_elas norm = {np.linalg.norm(np.array(n_elas.vector()))}")
print(f"beta norm = {np.linalg.norm(np.array(beta.vector()))}")

A, Res = fe.assemble_system(a_Newton, res, u_bcs)
nRes = Res.norm("l2")
print(f"Residual: {nRes}")
niter += 1

u.assign(u + Du)
sig_old.assign(sig)
p.assign(p + local_project(dp_, W0))
p_avg.assign(fe.project(p, P0))

T_old.assign(T_crt)

T_vtk_file << T_old
u_vtk_file << u
p_vtk_file << p_avg

file_results.write(T_old, i)
file_results.write(u, i)
file_results.write(p_avg, i)

if __name__ == '__main__':
simulation()
``````

However, the following error occurs in the runtime：

``````Moving new file over differing existing file:
Backup file exists, overwriting.
------------------- Start compiler output ------------------------
/tmp/tmpx7xoc3ez/ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16.cpp: In member function 'virtual void ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16_cell_integral_main_otherwise::tabulate_tensor(double*, const double* const*, const double*, int) const':
/tmp/tmpx7xoc3ez/ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16.cpp:124:29: error: 'FE17_C0_Q4' was not declared in this scope
124 |             TF0[iq] = fw0 * FE17_C0_Q4[0][iq][iq];
|                             ^~~~~~~~~~
/tmp/tmpx7xoc3ez/ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16.cpp:125:34: error: 'FE17_C0_Q4' was not declared in this scope
125 |         BF0[iq][iq] += TF0[iq] * FE17_C0_Q4[0][iq][iq];
|                                  ^~~~~~~~~~

-------------------  End compiler output  ------------------------
Compilation failed! Sources, command, and errors have been written to: /home/ljr/plasticity-main/src/jitfailure-ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16

---------------------------------------------------------------------------
DijitsoError                              Traceback (most recent call last)
Cell In[8], line 320
315             file_results.write(p_avg, i)
319 if __name__ == '__main__':
--> 320     simulation()

File ~/plasticity-main/src/utils.py:9, in walltime.<locals>.wrapper(*list_args, **keyword_wargs)
7 def wrapper(*list_args, **keyword_wargs):
8     start_time = time.time()
----> 9     func(*list_args, **keyword_wargs)
10     end_time = time.time()
11     time_elapsed = end_time - start_time

Cell In[8], line 287, in simulation()
283 print(f"du norm = {np.linalg.norm(np.array(du.vector()))}")
285 sig_, n_elas_, beta_, dp_ = proj_sig(deps, dT, sig_old, p)
--> 287 local_project(sig_, W, sig)
288 local_project(n_elas_, W, n_elas)
289 local_project(beta_, W0, beta)

Cell In[8], line 164, in simulation.<locals>.local_project(v, V, u)
162 a_proj = fe.inner(dv, v_)*dxm
163 b_proj = fe.inner(v, v_)*dxm
--> 164 solver = fe.LocalSolver(a_proj, b_proj)
165 solver.factorize()
166 if u is None:

File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/fem/solvers.py:39, in LocalSolver.__init__(self, a, L, solver_type)
36 self.L_ufl = L
38 # Wrap as DOLFIN forms
---> 39 a = Form(a)
40 if L is None:
41     # Initialize C++ base class
42     cpp.fem.LocalSolver.__init__(self, a, solver_type)

File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/fem/form.py:43, in Form.__init__(self, form, **kwargs)
40 if form_compiler_parameters is None:
41     form_compiler_parameters = {"external_include_dirs": dolfin_pc["include_dirs"]}
---> 43 ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
44                    mpi_comm=mesh.mpi_comm())
45 ufc_form = cpp.fem.make_ufc_form(ufc_form[0])
47 function_spaces = kwargs.get("function_spaces")

File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/jit/jit.py:47, in mpi_jit_decorator.<locals>.mpi_jit(*args, **kwargs)
45 # Just call JIT compiler when running in serial
46 if MPI.size(mpi_comm) == 1:
---> 47     return local_jit(*args, **kwargs)
49 # Default status (0 == ok, 1 == fail)
50 status = 0

File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/jit/jit.py:97, in ffc_jit(ufl_form, form_compiler_parameters)
95 p.update(dict(parameters["form_compiler"]))
96 p.update(form_compiler_parameters or {})
---> 97 return ffc.jit(ufl_form, parameters=p)

File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/ffc/jitcompiler.py:217, in jit(ufl_object, parameters, indirect)
214 kind, module_name = compute_jit_prefix(ufl_object, parameters)
216 # Inspect cache and generate+build if necessary
--> 217 module = jit_build(ufl_object, module_name, parameters)
219 # Raise exception on failure to build or import module
220 if module is None:
221     # TODO: To communicate directory name here, need dijitso params to call
222     #fail_dir = dijitso.cache.create_fail_dir_path(signature, dijitso_cache_params)

File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/ffc/jitcompiler.py:130, in jit_build(ufl_object, module_name, parameters)
123 params = dijitso.validate_params({
124     "cache": cache_params,
125     "build": build_params,
126     "generator": parameters,  # ffc parameters, just passed on to jit_generate
127 })
129 # Carry out jit compilation, calling jit_generate only if needed
--> 130 module, signature = dijitso.jit(jitable=ufl_object,
131                                 name=module_name,
132                                 params=params,
133                                 generate=jit_generate)
134 return module

File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dijitso/jit.py:216, in jit(jitable, name, params, generate, send, receive, wait)
214 if err_info:
215     # TODO: Parse output to find error(s) for better error messages
--> 216     raise DijitsoError("Dijitso JIT compilation failed, see '%s' for details"
217                        % err_info['fail_dir'], err_info)
219 # Return built library and its signature
220 return lib, signature

DijitsoError: Dijitso JIT compilation failed, see '/home/ljr/plasticity-main/src/jitfailure-ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16' for details
``````

Quadrature elements have always been fishy in the legacy FEniCS, yielding strange errors such as this one. I suggest switching to FEniCSx (dolfinx v0.8.0) in which everything works better

Thanks for your suggestion, I tried to use fenicsx to solve the thermoelasticity problem with the following code:

``````import numpy as np
import ufl
import dolfinx
from dolfinx import fem
from mpi4py import MPI
from petsc4py import PETSc
import glob
import os
import sys
import basix
from pprint import pprint
from arguments import args
from utils import walltime

comm = MPI.COMM_WORLD

def mpi_print(msg):
if comm.rank == 0:
print(f"Rank {comm.rank} print: {msg}")
sys.stdout.flush()

@walltime
def simulation():
case_name = 'mechanical'

if comm.rank == 0:
files_vtk = glob.glob(f'data/vtk/{case_name}/sols' + f"/*")
files_xdmf = glob.glob(f'data/xdmf/{case_name}/' + f"/*")
for f in files_vtk + files_xdmf:
os.remove(f)

ambient_T = args.T_ambient
rho = args.rho
Cp = args.c_p
k = args.kappa_T
h = args.h_conv
eta = args.power_fraction
r = args.r_beam
P = args.power

x0 = 0.2*args.domain_x
y0 = 0.5*args.domain_y

simulation_t = 2400*1e-6
total_t = 1200*1e-6
vel = 0.6*args.domain_x/total_t
dt = 1e-6
ts = np.arange(0., simulation_t + dt, dt)
mpi_print(f"total time steps = {len(ts)}")
ele_size = 0.01

Nx, Ny, Nz = round(args.domain_x/ele_size), round(args.domain_y/ele_size), round(args.domain_z/ele_size)

mpi_print(f"Nx = {Nx}, Ny = {Ny}, Nz = {Nz}")

mesh = dolfinx.mesh.create_box(MPI.COMM_WORLD, [np.array([0., 0., 0.]), np.array([args.domain_x, args.domain_y, args.domain_z])],
[Nx, Ny, Nz], cell_type=dolfinx.mesh.CellType.hexahedron)  # cell_type=mesh.CellType.hexahedron/tetrahedron

mesh_vtk_file = dolfinx.io.VTKFile(MPI.COMM_WORLD, f'data/vtk/{case_name}/mesh/mesh.pvd', 'w')
mesh_vtk_file.write_mesh(mesh)

# pprint(dir(mesh.geometry))
print(f"Total number of local mesh vertices {len(mesh.geometry.x)}" )

def bottom(x):
return np.isclose(x[2], 0.)

def top(x):
return np.isclose(x[2], args.domain_z)

fdim = mesh.topology.dim - 1
bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom)
top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top)

mpi_print(f"bottom_facets.shape = {bottom_facets.shape}")

marked_facets = np.hstack([bottom_facets, top_facets])
marked_values = np.hstack([np.full(len(bottom_facets), 1, dtype=np.int32), np.full(len(top_facets), 2, dtype=np.int32)])
sorted_facets = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(mesh, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

deg_u = 1
deg_stress = 2
degree_T = 1

# "quadrature_degree": 2 means that use 8 integrations ponits for a hexahedron element
normal = ufl.FacetNormal(mesh)

P0_ele = ufl.FiniteElement("DG", mesh.ufl_cell(), 0)
P0 = fem.FunctionSpace(mesh, P0_ele)
p_avg = fem.Function(P0, name="Plastic_strain")
strain_xx = fem.Function(P0, name="strain_xx")
stress_xx = fem.Function(P0, name="stress_xx")
phase_avg = fem.Function(P0, name="phase")

V_ele = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), degree=degree_T)
V = fem.FunctionSpace(mesh, V_ele)

U_ele = ufl.VectorElement("Lagrange", mesh.ufl_cell(), degree=deg_u)
U = fem.FunctionSpace(mesh, U_ele)

# W_ele = ufl.TensorElement("DG", mesh.ufl_cell(), 0)
# W = fem.FunctionSpace(mesh, W_ele)
# W0_ele = ufl.FiniteElement("DG", mesh.ufl_cell(), 0)
# W0 = fem.FunctionSpace(mesh, W0_ele)

W = fem.FunctionSpace(mesh, W_ele)
W0 = fem.FunctionSpace(mesh, W0_ele)

def ini_T(x):
return np.full(x.shape[1], ambient_T)

dT = fem.Function(V)

T_crt = fem.Function(V)
T_crt.interpolate(ini_T)
T_pre = fem.Function(V)
T_pre.interpolate(ini_T)
T_old = fem.Function(V, name='T')
T_old.interpolate(ini_T)

T_trial = ufl.TrialFunction(V)
T_test = ufl.TestFunction(V)

phase = fem.Function(V, name='phase')
alpha_V = fem.Function(V)
E = fem.Function(V)

# alpha_V.x.array[:] = 1e-5
# E.x.array[:] = args.Young_mod

nu = 0.3
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
sig0 = 250.
Et = E/100.
H = E*Et/(E-Et)

sig = fem.Function(W)
# Something like "Cumulative plastic strain" may cause an error due to the space - probably a bug of dolfinx
cumulative_p = fem.Function(W0, name="Cumulative_plastic_strain")
u = fem.Function(U, name="Total_displacement")
du = fem.Function(U, name="Iteration_correction")
Du = fem.Function(U, name="Current_increment")
v = ufl.TrialFunction(U)
u_ = ufl.TestFunction(U)

mpi_print(f"facet_tag.dim = {facet_tag.dim}")

bottom_dofs_u = fem.locate_dofs_topological(U, facet_tag.dim, bottom_facets)
bcs_u = [fem.dirichletbc(PETSc.ScalarType((0., 0., 0.)), bottom_dofs_u, U)]

def eps(v):
return e

def sigma(eps_el):
return lmbda*ufl.tr(eps_el)*ufl.Identity(3) + 2*mu*eps_el

deps = eps(Du)

def thermal_strain():
# alpha_V = 1e-5
return alpha_V*dT*ufl.Identity(3)

ppos = lambda x: (x + abs(x))/2.
heaviside = lambda x: ufl.conditional(ufl.gt(x, 0.), 1., 0.)

def proj_sig():
EPS = 1e-10
d_eps_T = thermal_strain()
sig_elas = sig + sigma(deps - d_eps_T)
s = ufl.dev(sig_elas)
sig_eq = ufl.sqrt(3/2.*ufl.inner(s, s))
f_elas = sig_eq - sig0 - H*cumulative_p
dp = ppos(f_elas)/(3*mu+H)
# Prevent divided by zero error
# didn't consider this, and can cause nan error in the solver.
n_elas = s/(sig_eq + EPS)*heaviside(f_elas)
beta = 3*mu*dp/(sig_eq + EPS)
new_sig = sig_elas - beta*s

return new_sig, n_elas, beta, dp

def sigma_tang(e):
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*ufl.inner(n_elas, e)*n_elas  -2*mu*beta*ufl.dev(e)

# If theta = 0., we recover implicit Eulear; if theta = 1., we recover explicit Euler; theta = 0.5 seems to be a good choice.
theta = 0.5
T_rhs = theta*T_pre + (1 - theta)*T_trial

bottom_dofs_T = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.indices[facet_tag.values==1])
bcs_T = [fem.dirichletbc(PETSc.ScalarType(ambient_T), bottom_dofs_T, V)]

x = ufl.SpatialCoordinate(mesh)
crt_time = fem.Constant(mesh, PETSc.ScalarType(0.))

q_laser = 2*P*eta/(np.pi*r**2) * ufl.exp(-2*((x[0] - x0 - vel*crt_time)**2 + (x[1] - y0)**2) / r**2) * ufl.conditional(ufl.gt(crt_time, total_t), 0., 1.)
# q_laser = 2*P*eta/(np.pi*r**2) * ufl.exp(-2*((x[0] - x0 - vel*crt_time)**2 + (x[1] - y0)**2) / r**2) * heaviside(total_t - crt_time.value)

q_convection = h * (T_rhs - ambient_T)
res_T = rho*Cp/dt*(T_trial - T_pre) * T_test * dxm + k * ufl.dot(ufl.grad(T_rhs), ufl.grad(T_test)) * dxm \
- q_laser * T_test * ds(2) - q_convection * T_test * ds

new_sig, n_elas, beta, dp = proj_sig()

# ufl diff might be used to automate the computation of tangent stiffness tensor
res_u_lhs = ufl.inner(eps(v), sigma_tang(eps(u_)))*dxm
res_u_rhs = -ufl.inner(new_sig, eps(u_))*dxm

problem_T = fem.petsc.LinearProblem(ufl.lhs(res_T), ufl.rhs(res_T), bcs=bcs_T, petsc_options={"ksp_type": "bicg", "pc_type": "jacobi"})

problem_u = fem.petsc.LinearProblem(res_u_lhs, res_u_rhs, bcs=bcs_u, petsc_options={"ksp_type": "bicg", "pc_type": "jacobi"})

def l2_projection(v, V):
dv = ufl.TrialFunction(V)
v_ = ufl.TestFunction(V)
a_proj = ufl.inner(dv, v_)*dxm
b_proj = ufl.inner(v, v_)*dxm
problem = fem.petsc.LinearProblem(a_proj, b_proj, petsc_options={"ksp_type": "bicg", "pc_type": "jacobi"})
u = problem.solve()
return u

'''
See https://github.com/FEniCS/dolfinx/issues/2243
'''
u = fem.Function(V)
map_c = mesh.topology.index_map(mesh.topology.dim)
num_cells = map_c.size_local + map_c.num_ghosts
cells = np.arange(0, num_cells, dtype=np.int32)
e_eval = e_expr.eval(cells)

with u.vector.localForm() as u_local:
u_local.setBlockSize(u.function_space.dofmap.bs)

return u

def update_modului():
# 0: powder, 1: liquid, 2: solid
T_array = T_crt.x.array

powder_to_liquid = (phase.x.array == 0) & (T_array > args.T_melt)
liquid_to_solid = (phase.x.array == 1) & (T_array < args.T_melt)

phase.x.array[powder_to_liquid] = 1
phase.x.array[liquid_to_solid] = 2

# print(f"number of powder = {np.sum(phase.x.array == 0)}, liquid = {np.sum(phase.x.array == 1)}, solid = {np.sum(phase.x.array == 2)}")

E.x.array[(phase.x.array == 0) | (phase.x.array == 1)]  = 1e-2*args.Young_mod
E.x.array[phase.x.array == 2] = args.Young_mod

alpha_V.x.array[(phase.x.array == 0) | (phase.x.array == 1)] = 0.

alpha_V.x.array[phase.x.array == 2] = args.alpha_V

def write_sol(file, step):
file.write_function(T_old, step)
file.write_function(u, step)
file.write_function(p_avg, step)
file.write_function(strain_xx, step)
file.write_function(stress_xx, step)
file.write_function(phase_avg, step)

vtk_file = dolfinx.io.VTKFile(mesh.comm, f'data/vtk/{case_name}/sols/u.pvd', 'w')
xdmf_file = dolfinx.io.XDMFFile(mesh.comm, f'data/xdmf/{case_name}/u.xdmf', 'w')

xdmf_file.write_mesh(mesh)

write_sol(vtk_file, 0)
write_sol(xdmf_file, 0)

plastic_inverval = 20

for i in range(len(ts) - 1):
# for i in range(20):

mpi_print(f"step {i + 1}/{len(ts) - 1}, time = {ts[i + 1]}")
crt_time.value = theta*ts[i] + (1 - theta)*ts[i + 1]

update_modului()

T_crt = problem_T.solve()

T_pre.x.array[:] = T_crt.x.array

# print(f"min T = {np.min(np.array(T_pre.x.array))}")
# print(f"max T = {np.max(np.array(T_pre.x.array))}\n")

if (i + 1) % plastic_inverval == 0:

T_crt_array = np.array(T_crt.x.array)
T_crt_array = np.where(T_crt_array < args.T_ambient, args.T_ambient, T_crt_array)
T_crt_array = np.where(T_crt_array > args.T_melt, args.T_melt, T_crt_array)
T_old_array = np.array(T_old.x.array)
T_old_array = np.where(T_old_array < args.T_ambient, args.T_ambient, T_old_array)
T_old_array = np.where(T_old_array > args.T_melt, args.T_melt, T_old_array)
dT.x.array[:] = T_crt_array - T_old_array

Du.x.array[:] = 0.

niter = 0
nRes = 1.
tol = 1e-8

while nRes > tol or niter < 2:
mpi_print(f"At iteration step {niter + 1}")
du = problem_u.solve()
Du.x.array[:] = Du.x.array + du.x.array
mpi_print(f"du norm = {np.linalg.norm(du.x.array)}")

# nRes1 = np.sqrt(mesh.comm.allreduce(np.sum(problem_u.b.array**2), op=MPI.SUM))
nRes = problem_u.b.norm(1)
mpi_print(f"b norm: {nRes}\n")
niter += 1

u.x.array[:] = u.x.array + Du.x.array

mpi_print(f"sig dof = {sig.x.array.shape}")
mpi_print(f"sig norm = {np.linalg.norm(sig.x.array)}")

cumulative_p.x.array[:] = cumulative_p.x.array + quad_interpolation(dp, W0).x.array

# Remark: Can we do interpolation here?
# p_avg.interpolate(fem.Expression(cumulative_p, P0.element.interpolation_points))
p_avg.x.array[:] = l2_projection(cumulative_p, P0).x.array
stress_xx.x.array[:] = l2_projection(sig[0, 0], P0).x.array
phase_avg.x.array[:] = l2_projection(phase, P0).x.array

T_old.x.array[:] = T_crt.x.array

write_sol(vtk_file, i + 1)
write_sol(xdmf_file, i + 1)

if __name__ == '__main__':
simulation()
``````

Where the code for “arguments” is as follows:

``````##### import numpy as onp
import argparse
import sys
import numpy as onp
import matplotlib.pyplot as plt

# Set numpy printing format
onp.random.seed(0)                  #指定随机数生成时所用算法开始的整数值
onp.set_printoptions(threshold=sys.maxsize, linewidth=1000, suppress=True)     #控制输出方式：显示所有数组元素、每行字符的数目为1000、小数以科学计数法的形式输出
onp.set_printoptions(precision=10)                                             #控制输出方式：输出结果的精度（小数点后的位数）为10

# np.set_printoptions(threshold=sys.maxsize, linewidth=1000, suppress=True)
# np.set_printoptions(precision=5)

# Manage arguments
parser = argparse.ArgumentParser()                     #创建一个解析器，ArgumentParser对象包含将命令行解析成Python数据类型的全部信息

# parser.add_argument('--domain_x', type=float, help='Unit: mm', default=0.05)
# parser.add_argument('--domain_y', type=float, help='Unit: mm', default=0.02)
# parser.add_argument('--domain_z', type=float, help='Unit: mm', default=0.005)

# parser.add_argument('--laser_vel', type=float, help='Unit: mm/s', default=500.)
parser.add_argument('--write_sol_interval', type=int, help='interval of writing solutions to file', default=500)

args = parser.parse_args(args=[])

# Latex style plot              #更新全局绘图参数
plt.rcParams.update({
"text.latex.preamble": r"\usepackage{amsmath}",
"text.usetex": True,
"font.family": "sans-serif",
"font.sans-serif": ["Helvetica"]})
``````

However, the following error occurs during runtime:

``````usage: ipykernel_launcher.py [-h] [--num_oris NUM_ORIS]
[--num_grains NUM_GRAINS] [--dim DIM]
[--domain_x DOMAIN_X] [--domain_y DOMAIN_Y]
[--domain_z DOMAIN_Z] [--dt DT] [--T_melt T_MELT]
[--T_ambient T_AMBIENT] [--rho RHO] [--c_p C_P]
[--alpha_V ALPHA_V] [--Young_mod YOUNG_MOD]
[--power POWER] [--power_fraction POWER_FRACTION]
[--r_beam R_BEAM] [--emissivity EMISSIVITY]
[--SB_constant SB_CONSTANT] [--h_conv H_CONV]
[--kappa_T KAPPA_T]
[--write_sol_interval WRITE_SOL_INTERVAL]
ipykernel_launcher.py: error: unrecognized arguments: -f /home/ljr/.local/share/jupyter/runtime/kernel-418ab754-bea6-4cd3-9a12-d50bacc0d21e.json

An exception has occurred, use %tb to see the full traceback.

SystemExit: 2

WARNING:py.warnings:/home/ljr/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/IPython/core/interactiveshell.py:3561: UserWarning: To exit: use 'exit', 'quit', or Ctrl-D.
warn("To exit: use 'exit', 'quit', or Ctrl-D.", stacklevel=1)
``````

Sorry but this is far from a MWE, please make one and identify exactly on a much simpler version what is your error

You would have to switch the form_compiler representation back to “quadrature” if you want to use quadrature elements in legacy dolfin.