Thanks a ton ! It is working now.
How to change the limits of the colorbar in such plots ? I haven’t used matplotlib much before and very new to it.
Thanks again.
Matplotlib has an extensive documentation with plenty of examples, see https://matplotlib.org/3.2.1/api/_as_gen/matplotlib.pyplot.colorbar.html
1 Like
Hi,
I modified the Elasto-plastic demo to solve a thermomechanical problem in a staggered manner. I need to reset the equivalent plastic strain, p=0, once the temperature exceeds a certain value(e.g. 1514 Kelvin )
I have tried using an Expression that depends on the temperature field T
# thermal part
scalar_elem = FiniteElement('P', triangle, 2)
V_T = FunctionSpace(mesh, scalar_elem)
T = Function(V_T , name='T')
T_ann = Expression('temp > 1514 ? 0 :1', temp= T, degree=1)
in order to manipulate p :
# mechanical part
key = Function(W0) # to control p
for (i, t) in enumerate(load_steps):
....
p.assign(p+local_project(dp_, W0))
local_project(T_ann, V_T, key)
# reassignment
p.assign(p*key)
and I get the following error because of p.assign(p*key)
RuntimeError: Expected only linear combinations of Functions in the same FunctionSpaces
is there any way to do such manipulation correctly?
Thanks in advance
Ali
Please supply a minimal example as described at Read before posting: How do I get my question answered? as it for instance is not clear what kind of variable key is.
Thank you for your reply.
because the Function spaces of the thermal and mechanical problems are different, I thought of defining a new parameter key using the quadrature space function W0 and then locally project the Expression using the local_project function. At the key should represent T_ann on the Quadrature spaces.
I would like to provide MWE but with plasticity demo it will be a long code. I will post it once it is ready.
The code does not Need to solve the whole problem, only Give the error. Thus, if the product consist of two functions from different spaces, just define the two functions and do the assign operation .
Hi,
If p and key are both defined on the quadrature function space W0 you can perform the multiplication on the array of dofs directly:
p.vector().set_local(p.vector().get_local()*key.vector().get_local())
2 Likes
Thank you so much. This works now.
Hi @bleyerj thanks for the notification. I have a question that does this implementation https://comet-fenics.readthedocs.io/en/latest/demo/plasticity_mfront/plasticity_mfront.py.html has ability to take care for the calculation based upon load history i.e. when after the loading the object is going back to unloading or the direction of load is reversed? Thank you!
Yes, just change the load stepping and it should work.
Best
Hi @bleyerj I am having a similar problem with the Newton-Raphson method not converging when removing the forcing function and setting the DirichletBC to zero for one side. I am trying to introduce a thermal strain and it does not converge, ||residual||,L2=0, when I use an input thermal field, but when I put a constant thermal strain it starts to solve but the residual is quite large, so Im guessing it has something to do with the calculation of the displacement increment and adjusting incremental values during iterations but Im just not sure what.
I added these two strain functions to the example code, and the only other change is removing the forcing function. Where dT is the thermal field from solving the heat conduction equation, and by replacing with a constant there is some action by the solver.
def eps(v):
e = sym(nabla_grad(v))
return as_tensor([[e[0, 0], e[0, 1], 0],
[e[0, 1], e[1, 1], 0],
[0, 0, 0]])
def eps_el(v,dT):
e = sym(nabla_grad(v))
eT=alphaExp*(dT-u_0)
return as_tensor([[e[0, 0]-eT, e[0, 1], 0],
[e[0, 1], e[1, 1]-eT, 0],
[0, 0, 0]])
A bit more detail:
metadata = {"quadrature_degree": 2, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
eq = inner(eps(vu), sigma_tang(eps(u_)))*dxm +inner(eps(u_), as_3D_tensor(sig))*dxm
a_Newton = lhs(eq)
res = rhs(eq)
itermax, tol = 15, 1e-8 # parameters of the Newton-Raphson procedure
Nincr = 20
A, Res = assemble_system(a_Newton, res,bcu4)
bcu4.homogenize()
nRes0 = Res.norm("l2")
nRes = nRes0 +DOLFIN_EPS
Du.interpolate(Constant((0, 0)))
#NewtonSolver(res, Du, bcu4, a_Newton)
print("Increment:", str(n+1))
niter = 0
print(" Residual:", nRes)
niter += 1
try:
#nRes/(nRes0)
while niter < Nitermax: #nRes/(nRes0) > tol and
solve(A, du.vector(), Res, "gmres")
Du.assign(Du+du)
deps = eps_el(Du,u)
sig_, n_elas_,beta_, dp_ = proj_sig(u_n,deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
A, Res = assemble_system(a_Newton, res, bcu4)
nRes = Res.norm("l2")
print(" Residual:", nRes)
niter += 1
if nRes < tol:
break
except ZeroDivisionError as error:
print('Zero Residual...')Hi,
sorry but it is hard to help without a complete MWE which I can run. Are you sure that your thermal strain is non-zero initially ? The residual should have a non-zero L2 norm at the beginning of the computation for a non-zero thermal strain.
It’s quite long, Im more than happy to share. I can see the thermal field when outputting the results, so I know there is some thermal strain. The error I get is
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
But here is the code…
"""
2D SIMULATION
Heating of a rectangular block from top surface and subsequent equilibration.
Solve the following PDE for u (temperature field)
Solving of quasi-static thermo-elasto-plastic PDE with input temperature field.
"""
from fenics import *
from petsc4py import PETSc
import os
import argparse
from ufl import nabla_div
import numpy as np
import time
import sys
from matplotlib import pyplot as plt
sys.getrecursionlimit()
sys.setrecursionlimit(1000000)
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
#choose scanning type
parser = argparse.ArgumentParser()
parser.add_argument("-s", "--spot-source", action="store_true", default=False)
args = parser.parse_args()
output_folder = './runs'
file_results = XDMFFile(output_folder+"/2Dthermoplasticity_results.xdmf")
#test change
num_steps = 50 # number of time steps
end_wait = 0 # number of end wait time steps (laser off)
dt = 1e-3 # time step (in second)
t_final = num_steps * dt # final time (in second)
printfreq = 1
set_log_level(20) # CRITICAL=50, ERROR=40, WARNING=30, INFO=20, PROGRESS=16, TRACE=13, DBG=10
#room temperature thermal properties Cu
k = 385 # thermal conductivity, in W/m-K
c = 0.385 # heat capacity, in J/g-K
rho = 8.93e6 # density, in g/m^3
alpha = k / (c*rho) # thermal diffusivity, in m^2/s
u_0 = 273.15 # initial temperature, in K
#room temperature mechanical properties
nu = Constant(0.3) #Poisson's Ratio
E = Constant(500e6) #Elastic Modulus in Pascal
H=0.1*E #Hardening Modulus
mu = E/2/(1+nu) #Lame Constant
lambda_ = E*nu/(1+nu)/(1-2*nu) #Lame Constant
alphaExp = Constant(.01) #Thermal Expansion m/K?
sig0 = Constant(25e5) #yield strength
#Laser parameters
scan_power = 2000 # laser power, in W
scan_vel = 0.3 #speed in m/s
scan_tlim = 1 # time period of illumination, in second
sigma_ = 2e-3 #1/4 beam diameter
line_length = 0.1 #length of laser line
# Create mesh and define function space or import mesh
Nx, Ny,= 50, 50 #number of nodes
L = 0.05 # sample length, in m
mesh = RectangleMesh(Point(-L, -L/2), Point(L, 0), Nx, Ny, "crossed")
if False:
plot(mesh)
plt.savefig('mesh.png')
# Define boundary conditions
class Left(SubDomain):
def inside(self,x, on_boundary):
return on_boundary and near(x[0], -L)
class Right(SubDomain):
def inside(self,x, on_boundary):
return on_boundary and near(x[0], L)
class Top(SubDomain):
def inside(self,x, on_boundary):
return on_boundary and near(x[1], 0)
class Bottom(SubDomain):
def inside(self,x, on_boundary):
return on_boundary and near(x[1], -L/2)
# Add boundaries
left = Left()
right = Right()
top=Top()
bottom=Bottom()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
top.mark(boundaries,3)
bottom.mark(boundaries,4)
#Define function spaces for both problems
V = FunctionSpace(mesh, 'P', 1)
Vu = VectorFunctionSpace(mesh, 'CG', 2)
# Function space for properties defined on Gaussian Quadrature points
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=2,dim=4, quad_scheme='default')
W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=2, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)
sig = Function(W, name="Stress field on Gauss points")
sig_old = Function(W)
n_elas = Function(W)
beta = Function(W0)
p = Function(W0, name="Cumulative plastic strain")
uE = Function(Vu, name="Total displacement")
du = Function(Vu, name="Iteration correction")
Du = Function(Vu, name="Current increment")
vu = TrialFunction(Vu)
u_ = TestFunction(Vu)
u = TrialFunction(V)
v = TestFunction(V)
# Define initial values for thermal problems
u_n = interpolate(Constant(u_0), V)
#Potential boundary conditions for both problems
bcT1 = DirichletBC(V, Constant(u_0), boundaries, 1)
bcT2 = DirichletBC(V, Constant(u_0), boundaries, 2)
bcT3 = DirichletBC(V, Constant(u_0), boundaries, 3)
bcT4 = DirichletBC(V, Constant(u_0), boundaries, 4)
bcu1 = DirichletBC(Vu, Constant((0., 0.)), boundaries, 1)
bcu2 = DirichletBC(Vu, Constant((0., 0.)), boundaries, 2)
bcu3 = DirichletBC(Vu, Constant((0., 0.)), boundaries, 3)
bcu4 = DirichletBC(Vu, Constant((0., 0.)), boundaries, 4)
# total strain
def eps(v):
e = sym(nabla_grad(v))
return as_tensor([[e[0, 0], e[0, 1], 0],
[e[0, 1], e[1, 1], 0],
[0, 0, 0]])
def eps_el(v,dT):
e = sym(nabla_grad(v))
eT=alphaExp*(500-u_0)
return as_tensor([[e[0, 0]-eT, e[0, 1], 0],
[e[0, 1], e[1, 1]-eT, 0],
[0, 0, 0]])
#assumes elastic strain only
##def sigma(eps_el,dT):
## return (lambda_*tr(eps_el)- alphaExp*(3*lambda_+2*mu)*(dT-u_0))*Identity(3) + 2.0*mu*eps_el
def sigma(eps_el):
return (lambda_*tr(eps_el))*Identity(3) + 2*mu*eps_el
def as_3D_tensor(X):
return as_tensor([[X[0], X[3], 0],
[X[3], X[1], 0],
[0, 0, X[2]]])
def sigma_tang(e):
N_elas = as_3D_tensor(n_elas)
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*inner(N_elas, e)*N_elas-2*mu*beta*dev(e)
#return mapping procedure
ppos = lambda x: (x+abs(x))/2.
def proj_sig(u,deps, old_sig, old_p): #instead of strain put displacement then calculate here so we can sub thermal strain
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + sigma(deps)
s = dev(sig_elas)
sig_eq = sqrt(3/2.*inner(s, s))
f_elas = sig_eq - sig0 - H*old_p
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_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1]]), \
as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1]]), \
beta, dp
#define laser function and parameters
x_start=-3*L/4
if args.spot_source:
print("use moving spot source for heating")
output_folder = output_folder + '_spot_source'
spot_intensity=scan_power/(2*pi*sigma_*sigma_) / (c*rho)
g = Expression('(t<=scan_tlim)?(spot_intensity*exp(-(x[0]-xcenter)*(x[0]-xcenter)/(2*sigma*sigma) - (x[1]-ycenter)*(x[1]-ycenter)/(2*sigma*sigma))):0',
degree=2, ycenter=0,xcenter=x_start,scan_tlim=scan_tlim,sigma=sigma_,spot_intensity=spot_intensity,t=0)
else:
print("use moving line_source for heating")
output_folder = output_folder + '_line_source'
scan_intensity = scan_power / (line_length) / (c * rho)
g = Expression('(fabs(x[1])<=line_length/2 & t<=scan_tlim)?(scan_intensity/(sqrt(2*M_PI)*sigma)*exp(-(x[0]-xcenter)*(x[0]-xcenter)/(2*sigma*sigma))):0',
degree=2, t=0, scan_tlim=scan_tlim,scan_intensity=scan_intensity,xcenter=x_start,sigma=sigma_,line_length=line_length)
####
def NewtonSolver(G, u, bcs, dG):
tol = 1e-10
error = 1.
it = 0
bcs.apply(u.vector())
bcs.homogenize()
while error > tol and it < 10:
# Assemble
dg, g = assemble_system(dG, G, bcs)
# PETSc solver
dGp = as_backend_type(dg).mat()
x, b = dGp.getVecs()
b.setValues(range(g.size()), g)
ksp = PETSc.KSP().create()
ksp.setType('gmres')
ksp.setOperators(dGp)
ksp.solve(b, x)
error = g.norm('l2')
print('Iteration:', str(it), '; Error: %3.3e' % error)
if error < tol:
break
# Update
u.vector()[:] = u.vector()[:] - x
it += 1
####
#thermal problem
ds = Measure("ds", domain = mesh, subdomain_data = boundaries)
F = u*v*dx + dt*alpha*dot(grad(u), grad(v))*dx - (u_n)*v*dx - dt*(g)*v*ds(3)
a, L = lhs(F), rhs(F)
u = Function(V,name="Temperature")
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
#local solver for quadrature space
def local_project(v, V, u=None):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dxm
b_proj = inner(v, v_)*dxm
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return
metadata = {"quadrature_degree": 2, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
eq = inner(eps(vu), sigma_tang(eps(u_)))*dxm +inner(eps(u_), as_3D_tensor(sig))*dxm
a_Newton = lhs(eq)
res = rhs(eq)
t = 0
start = time.time()
#compute solution
for n in range(num_steps+end_wait):
#update position
if n<=num_steps:
#update current position (replace with text file locations)
g.xcenter=x_start+n*dt*scan_vel
else:
g.spot_intensity=0
# Compute thermal solution with boundary conditions
solve(a==L, u)
# Update previous solution
u_n.assign(u)
##################################SOLID SOLUTION##############################
#solid problem
Nitermax, tol = 15, 1e-8 # parameters of the Newton-Raphson procedure
Nincr = 20
A, Res = assemble_system(a_Newton, res,bcu4)
bcu4.homogenize()
nRes0 = Res.norm("l2")
nRes = nRes0 +DOLFIN_EPS
Du.interpolate(Constant((0, 0)))
#NewtonSolver(res, Du, bcu4, a_Newton)
print("Increment:", str(n+1))
niter = 0
print(" Residual:", nRes)
niter += 1
try:
#nRes/(nRes0)
while niter < Nitermax: #nRes/(nRes0) > tol and
solve(A, du.vector(), Res, "gmres")
Du.assign(Du+du)
deps = eps_el(Du,u)
sig_, n_elas_,beta_, dp_ = proj_sig(u_n,deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
A, Res = assemble_system(a_Newton, res, bcu4)
nRes = Res.norm("l2")
print(" Residual:", nRes)
niter += 1
if nRes < tol:
break
except ZeroDivisionError as error:
print('Zero Residual...')
uE.assign(uE+Du) #need to change u here to something else
sig_old.assign(sig)
p.assign(p+local_project(dp_, W0))
file_results.parameters["flush_output"] = True
file_results.parameters['rewrite_function_mesh'] = False
file_results.parameters["functions_share_mesh"] = True
file_results.write(u,t)
file_results.write(uE,t)
p_avg.assign(project(p, P0))
file_results.write(p_avg, t)
#save solution
if n % printfreq == 0 or n == num_steps-1:
print("n = %d / %d"%(n, num_steps))
# Save to file
file_results.parameters["flush_output"] = True
file_results.parameters['rewrite_function_mesh'] = False
file_results.parameters["functions_share_mesh"] = True
file_results.write(u,t)
file_results.write(uE,t)
p_avg.assign(project(p, P0))
file_results.write(p_avg, t)
# Update current time
t += dt
g.t = t
file_results.close()
end = time.time()
print(end - start)Thanks but I suggest you first remove as much as possible from the script so that it will be easier to debug. I understand that the plasticity script will not necessarily be minimal but at least remove unnecessary stuff such as moving spot sources etc…
Also, do you succeed in solving a linear problem without plasticity ?
I have successfully solved a thermoelastic problem using the following relations:
# Define strain and stress
def epsilon(u):
return sym(nabla_grad(u))
def sigma(v, dT):
return (lambda_*tr(epsilon(v))- alphaExp*(3*lambda_+2*mu_)*(dT-u_0))*Identity(3) + 2.0*mu_*epsilon(v)
FE = inner(sigma(du,u), epsilon(u_))*dx
aE, LE = lhs(FE), rhs(FE)Here is a stationary source:
"""
2D SIMULATION
Heating of a rectangular block from top surface and subsequent equilibration.
Solve the following PDE for u (temperature field)
Solving of quasi-static thermo-elasto-plastic PDE with input temperature field.
"""
from fenics import *
from petsc4py import PETSc
import os
import argparse
from ufl import nabla_div
import numpy as np
import time
import sys
from matplotlib import pyplot as plt
sys.getrecursionlimit()
sys.setrecursionlimit(1000000)
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
output_folder = './runs'
file_results = XDMFFile(output_folder+"/2Dthermoplasticity_results.xdmf")
#test change
num_steps = 50 # number of time steps
end_wait = 0 # number of end wait time steps (laser off)
dt = 1e-3 # time step (in second)
t_final = num_steps * dt # final time (in second)
printfreq = 1
set_log_level(20) # CRITICAL=50, ERROR=40, WARNING=30, INFO=20, PROGRESS=16, TRACE=13, DBG=10
#room temperature thermal properties Cu
k = 385 # thermal conductivity, in W/m-K
c = 0.385 # heat capacity, in J/g-K
rho = 8.93e6 # density, in g/m^3
alpha = k / (c*rho) # thermal diffusivity, in m^2/s
u_0 = 273.15 # initial temperature, in K
#room temperature mechanical properties
nu = Constant(0.3) #Poisson's Ratio
E = Constant(500e6) #Elastic Modulus in Pascal
H=0.1*E #Hardening Modulus
mu = E/2/(1+nu) #Lame Constant
lambda_ = E*nu/(1+nu)/(1-2*nu) #Lame Constant
alphaExp = Constant(.01) #Thermal Expansion m/K?
sig0 = Constant(25e5) #yield strength
#Laser parameters
scan_power = 2000 # laser power, in W
scan_vel = 0.3 #speed in m/s
scan_tlim = 1 # time period of illumination, in second
sigma_ = 2e-3 #1/4 beam diameter
line_length = 0.1 #length of laser line
# Create mesh and define function space or import mesh
Nx, Ny,= 50, 50 #number of nodes
L = 0.05 # sample length, in m
mesh = RectangleMesh(Point(-L, -L/2), Point(L, 0), Nx, Ny, "crossed")
if False:
plot(mesh)
plt.savefig('mesh.png')
# Define boundary conditions
class Left(SubDomain):
def inside(self,x, on_boundary):
return on_boundary and near(x[0], -L)
class Right(SubDomain):
def inside(self,x, on_boundary):
return on_boundary and near(x[0], L)
class Top(SubDomain):
def inside(self,x, on_boundary):
return on_boundary and near(x[1], 0)
class Bottom(SubDomain):
def inside(self,x, on_boundary):
return on_boundary and near(x[1], -L/2)
# Add boundaries
left = Left()
right = Right()
top=Top()
bottom=Bottom()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
top.mark(boundaries,3)
bottom.mark(boundaries,4)
#Define function spaces for both problems
V = FunctionSpace(mesh, 'P', 1)
Vu = VectorFunctionSpace(mesh, 'CG', 2)
# Function space for properties defined on Gaussian Quadrature points
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=2,dim=4, quad_scheme='default')
W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=2, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)
sig = Function(W, name="Stress field on Gauss points")
sig_old = Function(W)
n_elas = Function(W)
beta = Function(W0)
p = Function(W0, name="Cumulative plastic strain")
uE = Function(Vu, name="Total displacement")
du = Function(Vu, name="Iteration correction")
Du = Function(Vu, name="Current increment")
vu = TrialFunction(Vu)
u_ = TestFunction(Vu)
u = TrialFunction(V)
v = TestFunction(V)
# Define initial values for thermal problems
u_n = interpolate(Constant(u_0), V)
#Potential boundary conditions for both problems
bcT1 = DirichletBC(V, Constant(u_0), boundaries, 1)
bcT2 = DirichletBC(V, Constant(u_0), boundaries, 2)
bcT3 = DirichletBC(V, Constant(u_0), boundaries, 3)
bcT4 = DirichletBC(V, Constant(u_0), boundaries, 4)
bcu1 = DirichletBC(Vu, Constant((0., 0.)), boundaries, 1)
bcu2 = DirichletBC(Vu, Constant((0., 0.)), boundaries, 2)
bcu3 = DirichletBC(Vu, Constant((0., 0.)), boundaries, 3)
bcu4 = DirichletBC(Vu, Constant((0., 0.)), boundaries, 4)
# total strain
def eps(v):
e = sym(nabla_grad(v))
return as_tensor([[e[0, 0], e[0, 1], 0],
[e[0, 1], e[1, 1], 0],
[0, 0, 0]])
def eps_el(v,dT):
e = sym(nabla_grad(v))
eT=alphaExp*(500-u_0)
return as_tensor([[e[0, 0]-eT, e[0, 1], 0],
[e[0, 1], e[1, 1]-eT, 0],
[0, 0, 0]])
#assumes elastic strain only
##def sigma(eps_el,dT):
## return (lambda_*tr(eps_el)- alphaExp*(3*lambda_+2*mu)*(dT-u_0))*Identity(3) + 2.0*mu*eps_el
def sigma(eps_el):
return (lambda_*tr(eps_el))*Identity(3) + 2*mu*eps_el
def as_3D_tensor(X):
return as_tensor([[X[0], X[3], 0],
[X[3], X[1], 0],
[0, 0, X[2]]])
def sigma_tang(e):
N_elas = as_3D_tensor(n_elas)
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*inner(N_elas, e)*N_elas-2*mu*beta*dev(e)
#return mapping procedure
ppos = lambda x: (x+abs(x))/2.
def proj_sig(u,deps, old_sig, old_p): #instead of strain put displacement then calculate here so we can sub thermal strain
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + sigma(deps)
s = dev(sig_elas)
sig_eq = sqrt(3/2.*inner(s, s))
f_elas = sig_eq - sig0 - H*old_p
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_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1]]), \
as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1]]), \
beta, dp
#define laser function and parameters
x_start=-3*L/4
output_folder = output_folder + '_line_source'
scan_intensity = scan_power / (line_length) / (c * rho)
g = Expression('(fabs(x[1])<=line_length/2 & t<=scan_tlim)?(scan_intensity/(sqrt(2*M_PI)*sigma)*exp(-(x[0]-xcenter)*(x[0]-xcenter)/(2*sigma*sigma))):0',
degree=2, t=0, scan_tlim=scan_tlim,scan_intensity=scan_intensity,xcenter=x_start,sigma=sigma_,line_length=line_length)
#thermal problem
ds = Measure("ds", domain = mesh, subdomain_data = boundaries)
F = u*v*dx + dt*alpha*dot(grad(u), grad(v))*dx - (u_n)*v*dx - dt*(g)*v*ds(3)
a, L = lhs(F), rhs(F)
u = Function(V,name="Temperature")
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
#local solver for quadrature space
def local_project(v, V, u=None):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dxm
b_proj = inner(v, v_)*dxm
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return
metadata = {"quadrature_degree": 2, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
eq = inner(eps(vu), sigma_tang(eps(u_)))*dxm +inner(eps(u_), as_3D_tensor(sig))*dxm
a_Newton = lhs(eq)
res = rhs(eq)
t = 0
start = time.time()
#compute solution
for n in range(num_steps+end_wait):
#update position
if n<=num_steps:
#update current position (replace with text file locations)
g.xcenter=x_start
else:
g.spot_intensity=0
# Compute thermal solution with boundary conditions
solve(a==L, u)
# Update previous solution
u_n.assign(u)
##################################SOLID SOLUTION##############################
#solid problem
Nitermax, tol = 15, 1e-8 # parameters of the Newton-Raphson procedure
Nincr = 20
A, Res = assemble_system(a_Newton, res,bcu4)
bcu4.homogenize()
nRes0 = Res.norm("l2")
nRes = nRes0 +DOLFIN_EPS
Du.interpolate(Constant((0, 0)))
#NewtonSolver(res, Du, bcu4, a_Newton)
print("Increment:", str(n+1))
niter = 0
print(" Residual:", nRes)
niter += 1
try:
#nRes/(nRes0)
while niter < Nitermax: #nRes/(nRes0) > tol and
solve(A, du.vector(), Res, "gmres")
Du.assign(Du+du)
deps = eps_el(Du,u)
sig_, n_elas_,beta_, dp_ = proj_sig(u_n,deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
A, Res = assemble_system(a_Newton, res, bcu4)
nRes = Res.norm("l2")
print(" Residual:", nRes)
niter += 1
if nRes < tol:
break
except ZeroDivisionError as error:
print('Zero Residual...')
uE.assign(uE+Du) #need to change u here to something else
sig_old.assign(sig)
p.assign(p+local_project(dp_, W0))
file_results.parameters["flush_output"] = True
file_results.parameters['rewrite_function_mesh'] = False
file_results.parameters["functions_share_mesh"] = True
file_results.write(u,t)
file_results.write(uE,t)
p_avg.assign(project(p, P0))
file_results.write(p_avg, t)
#save solution
if n % printfreq == 0 or n == num_steps-1:
print("n = %d / %d"%(n, num_steps))
# Save to file
file_results.parameters["flush_output"] = True
file_results.parameters['rewrite_function_mesh'] = False
file_results.parameters["functions_share_mesh"] = True
file_results.write(u,t)
file_results.write(uE,t)
p_avg.assign(project(p, P0))
file_results.write(p_avg, t)
# Update current time
t += dt
g.t = t
file_results.close()
end = time.time()
print(end - start)Okay so somehow it pushed the temperature field when using the linear solver for temperature.
# Compute thermal solution with boundary conditions
# Setup iterative solver
problem = LinearVariationalProblem(a, L, u)
solver = LinearVariationalSolver(problem)
solver.parameters['linear_solver'] = 'cg'
solver.parameters['preconditioner'] = 'hypre_euclid'
solver.parameters["krylov_solver"]["relative_tolerance"] = 5e-4
solver.parameters["krylov_solver"]["monitor_convergence"] = True
#info(solver.parameters,True)
solver.solve()
# Update previous solution
u_n.assign(u)
Though the residual seems to slowly grow for the solid steps, I think indicating I am not correctly implementing the increment conditions…
Guys, in the original tutorial (Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation) it stated in a note that “We could have used conditionals to write more explicitly the difference between elastic and plastic evolution.”
How can I explicitly use conditionals to branch the elastic/plastic situations?
I tried something like this: flag_elastic = conditional(lt(sig_eq_trial, sig0), True, False)
But it seems that the variable flag_elastic is always true no matter how sig_eq_trial changes…