This page was generated from unit-6.3-plasticity/plasticity.ipynb.
6.3 Elasto-Plasticity¶
This tutorial considers Hencky-type [1] plasticity with isotropic hardening in 2d (plane strain). It showcases the classes IntegrationRuleSpace, NewtonCF and MinimizationCF.
We begin with a classical formulation with an explicit yield surface that is probably more familiar to engineers. In a second stage, we consider a (constrained) minimization formulation which leads to a more streamlined implementation. Both formulations correspond to the isotropic hardening model discussed in [2].
The essential model ingredients are given in terms of the dissipation potential \(\Phi\) and the stored energy density contribution \(\Psi\), see eg, [3]. For convenience, we first introduce the strain energy density \(\Psi^\text{e}\)
\begin{align*} \Psi^\text{e}(\epsilon) = \frac{1}{2}\| \epsilon\|^2_{M} = \frac{E}{2(1+\nu)(1-2\nu)}\left((1-2\nu)\epsilon:\epsilon + \nu \,\text{tr}(\epsilon)^2\right) = \frac{1}{2} \epsilon : \mathbf{C} : \epsilon. \end{align*}
Moreover, we consider the full strain \(\epsilon\) to be the sum of the elastic \(\epsilon^e\) and the plastic strain contribution \(p\). Conversely, the elastic strain is the given as \begin{equation} \epsilon^\text{e} = \epsilon - p ,\end{equation} which is the strain quantity that will be actually employed as argument to \(\Psi^\text{e}\).
The isotropic hardening variable will be denoted by \(\alpha\) throughout this tutorial. The hardening potential density is simply chosen as \(\Psi^\text{h} = \frac{1}{2}\alpha^2\) such that the total energy density reads \begin{equation} \Psi(\epsilon, p, \alpha) = \Psi^\text{e}(\epsilon-p) + \Psi^\text{h}(\alpha) = (\epsilon - p) : \mathbf{C} : (\epsilon - p) + \frac{1}{2} \alpha^2. \end{equation}
The yield function \(f\) that represents the elastic limit is given as
\begin{equation} f(\sigma, \beta) = \left|\mathrm{dev} \sigma\right| - \sqrt{2/3}\sigma_Y(1 - H \beta), \end{equation}
where we note that the hardening parameter \(H\) enters in \(f\) rather than in \(\Psi\) for consistency with [2]. Deviating from [2], we scale the yield stress \(\sigma_Y\) by \(\sqrt{2/3}\) to obtain the von-Mises yield surface radius. We remark that by virtue of \(\mathrm{dev}\sigma\) appearing in \(f\), the plastic strain (rate) \(p\) (\(\dot p\)) will be deviatoric.
The boundary value problem solved in this tutorial is in plain strains. However, this does not mean that the plastic strains are plain and neither are the elastic ones. Only their sum has vanishing out-of-plane components. Therefore, for simplicity, the material model is implemented in 3d whereas displacements are in 2d.
Formulation with explicit yield surface¶
Based on the principle of maximum dissipation we can abstractly formulate the dissipation potential as
:nbsphinx-math:`begin{align*} Phi(dot{p}, dot{alpha}) =
- sup_{sigma} sup_{beta} inf_{Lambda geq 0}
: sigma : dot{p} + beta dot{alpha} - Lambda f(sigma, beta)
end{align*}`
The stationary conditions of the above functional read as
- :nbsphinx-math:`begin{align*}
&dot{p} = Lambda frac{partial left|mathrm{dev} sigmaright|}{partial sigma} \ &dot{alpha} = Lambda,sqrt{2/3}sigma_YH \ &{Lambda geq 0 : land : (left|mathrm{dev} sigmaright| - sqrt{2/3}sigma_Y(1 - Hbeta)) leq 0 : land : Lambda(left|mathrm{dev} sigmaright| - sqrt{2/3}sigma_Y(1 - Hbeta)) = 0}
end{align*}`
which are regarded as evolution equations for \(p\) and \(\alpha\). The actual values of \(\sigma\) and \(\beta\) are obtained through the constitutive relations provided by the framework of generalized standard materials, we have
\begin{align*} &\frac{\partial\Psi((\epsilon - p), \alpha)}{\partial p} + \frac{\partial\Phi}{\partial \dot{p}} = 0 \quad \Rightarrow \quad \sigma = \mathbf{C} : (\epsilon - p) \\ &\frac{\partial\Psi((\epsilon - p), \alpha)}{\partial \alpha} + \frac{\partial\Phi}{\partial \dot{\alpha}} = 0 \quad \Rightarrow \quad \beta = - \alpha . \end{align*}
The time-discrete form is simply obtained by replacing the time rates with the approximations \(\dot p \approx \frac{p_{k+1} - p_k}{\Delta t}\) and \(\dot \alpha \approx \frac{\alpha_{k+1} - \alpha_k}{\Delta t}\), whereby the value of \(\Delta t\) has no effect for the present rate-independent model.
Constrained minimization formulation¶
where \(\Psi\) denotes the augmented stored energy density and \(\|\bullet \|_{\varepsilon}\) indicates a perturbed norm in order to avoid divisions by zero in the evaluation of the time-discrete evolution equation derived from an incremental variational principle.
Notes in the spatial discretization¶
The state \(u\) is spatially discretized by continuous Lagrange elements, the internal states \(p\) and \(\alpha\) as well as the multiplier \(\Lambda\) reside at quadrature points, for which NGSolve has the notion of an integration rule space. Thus, the problem above can be decomposed into a “global” problem for the displacements and local problems for the internal variables. The latter can be solved individually for each quadrature point. However, there is a coupling between
the local problems and displacements, which has to be accounted for by what is called “algorithmically consistent linearization”. The latter is automatically obtained by using an appropriate BilinearForm object representing the boundary value problem and the internal evolution problems as demonstrated for both implementations.
References 1. W. Han and B.D. Reddy: Mathematical Theory and Numerical Analysis, Springer 1999 2. C. Carstensen, Domain Decomposition for a Non-smooth Convex Minimization Problem and its Application in Plasticity, 1997 3. K. Hackl, F.D. Fischer, On the relation between the principle of maximum dissipation and inelastic evolution given by dissipation potentials, PRSA, 2008
Implementation¶
General setup¶
[1]:
import numpy as np
from ngsolve import *
from ngsolve.comp import IntegrationRuleSpace
from ngsolve.fem import MinimizationCF, NewtonCF
from netgen.geom2d import CSG2d, Circle, Rectangle
from ngsolve.webgui import Draw
SetNumThreads(4)
The problem domain is the classical plate with hole. Thanks to symmetry of the problem to be solved, we actually use only one quarter of the domain together with additional symmetry boundary conditions. Parameters for the geometry and the material employed are taken from [4].
[2]:
# polynomial order for geometry and displacements
order = 3
geo = CSG2d()
circle = Circle(center=(100,100), radius=10.0, bc="curve").Maxh(1)
rect = Rectangle(pmin=(0,100), pmax=(100,200), bottom="bottom", left="left", top="top", right="right")
geo.Add(rect-circle)
mesh = Mesh(geo.GenerateMesh(maxh=5))
mesh.Curve(order)
Draw(mesh)
# points of interest for post-processing
node_A = mesh(100,200)
node_B = mesh(0,200)
Strain energy density and generic helper functions
Model parameters
[3]:
# Shear modulus
mu = 80193.8
# Bulk modulus
kappa = 164206
# Young's modulus
Emod = Parameter((9 * kappa * mu) / (3 * kappa + mu))
# Poisson's ratio
nu = Parameter((3 * kappa - 2 * mu)/(2 * (3 * kappa + mu)))
print("E =", Emod, "\nnu =", nu)
# Hardening parameter
H = Parameter(10)
# Yield stress
sigma_Y = Parameter(450)
# Perturbation parameter for tensor norms: this value is decisive for the accuracy of the result!
# for higher values of H is must be considerably smaller
norm_pert = Parameter(1e-16)
E = ParameterCF, val = 206900
nu = ParameterCF, val = 0.29
[4]:
def strain(u):
_strain = Sym(Grad(u))
return CoefficientFunction((_strain[0,0], _strain[0,1], 0,
_strain[1,0], _strain[1,1], 0,
0 , 0, 0), dims=(3, 3))
def elastic_strain_energy(eps, Emod, nu):
return 1/2 * Emod / ((1 + nu) * (1 - 2 * nu)) * ((1 - 2 * nu) * InnerProduct(eps, eps) + nu * Trace(eps)**2)
def tensor_norm(a, pert=0):
return sqrt(InnerProduct(a, a) + pert)
def dev(a):
return a - Trace(a) * Id(3) / 3
Common FE spaces¶
We employ a \(H^1\) space for the displacements and a “integration rule space” (IR space) as basis space for all internal variables.
[5]:
fes_u = VectorH1(mesh, order=order, dirichletx="right", dirichlety="bottom")
fes_ir = IntegrationRuleSpace(mesh, order=order-1)
Scalar internal states can be defined based on fes_ir as needed.
The “measure” corresponding to fes_ir:
[6]:
# extract integration rules from IntegrationRuleSpace
irs_dx = dx(intrules=fes_ir.GetIntegrationRules())
Load definitions¶
[7]:
force = 450
# For loadsteps
loadfactor = Parameter(1)
Some data structures for postprocessing¶
[8]:
# for postprocessing
drawfes = MatrixValued(L2(mesh, order=order-1), dim=3, symmetric=True, deviatoric=False)
pd, qd = drawfes.TnT()
pdraw = GridFunction(drawfes)
ad = BilinearForm(InnerProduct(pd, qd)*irs_dx, symmetric=True).Assemble()
invad = ad.mat.Inverse()
drawafes = L2(mesh, order=order-1)
pda, qda = drawafes.TnT()
adraw = GridFunction(drawafes)
ada = BilinearForm(InnerProduct(pda, qda)*irs_dx, symmetric=True).Assemble()
invada = ada.mat.Inverse()
Implementation with explicit yield surface¶
In this approach we employ a generic symmetric-matrix-valued space for \(p\).
[9]:
fes_p = MatrixValued(fes_ir, symmetric=True, deviatoric=False, dim=3)
fes_int = fes_p * fes_ir * fes_ir # p x alpha x Lambda
# We need a coupled space as well.
fes = fes_u * fes_int
u, p, alpha, Lambda = fes.TrialFunction()
# GridFunction for full solution
gfsol = GridFunction(fes)
gfu = gfsol.components[0]
# GridFunction for internal states only
gfint = GridFunction(fes_int)
gfp, gfalpha, gfLambda = gfint.components
# For history states
fes_hist = fes_p * fes_ir
gfhist = GridFunction(fes_hist)
gfp_k, gfalpha_k = gfhist.components
# "Trial" states for plastic evolution (these are not trial functions in FE parlance)
gftrial = GridFunction(fes_hist)
gfsigma_trial, gfbeta_trial = gftrial.components
# Time increment
Delta_t = Parameter(1.0)
# Not really needed for a rate-independent model. Value could be used for scaling the equations though.
# Short-hand for storing history variables. Note that the history value of Lambda is not relevant and thus omitted.
def store_internal():
gfp_k.Interpolate(gfp)
gfalpha_k.Interpolate(gfalpha)
Next, we setup functions for the energy density \(\Psi\) and the dissipation potential \(\Phi\)
[10]:
# Energy density
def Psi(strain, p, alpha):
strain_energy = elastic_strain_energy(strain-p, Emod, nu)
hardening = 1/2 * alpha**2
return strain_energy + hardening
# The yield condition: result > 0 ==> plastic evolution
def yield_condition(sigma, beta, pert=0):
return tensor_norm(dev(sigma), pert=pert) - sqrt(2/3) * sigma_Y * (1 - H*beta)
# The objective function defining the dissipation potential
def Phi_obj(sigma, beta, p_dot, alpha_dot, Lambda):
return InnerProduct(sigma, p_dot) + InnerProduct(beta, alpha_dot) - \
Lambda * yield_condition(sigma, beta, pert=norm_pert)
# The evolution equations
def evolution_eqs(strain, p, alpha, Lambda, evol_tol):
p_dot = (p - gfp_k) / Delta_t
alpha_dot = (alpha - gfalpha_k) / Delta_t
sigma = -Psi(strain, p, alpha).Diff(p)
beta = -Psi(strain, p, alpha).Diff(alpha)
_Phi_obj = Phi_obj(sigma, beta, p_dot, alpha_dot, Lambda)
dsigma = _Phi_obj.Diff(sigma)
# only active when there is hardening
dbeta = IfPos(H - 1e-16, _Phi_obj.Diff(beta), alpha_dot)
dLambda = _Phi_obj.Diff(Lambda)
# Only solve the evolution equations if the elastic trial state is beyond the yield surface.
# This essentially leads to the fulfillment of Lambda >= 0
return IfPos(yield_condition(gfsigma_trial, gfbeta_trial) - evol_tol,
CoefficientFunction((dsigma, dbeta, dLambda)),
CoefficientFunction((p_dot, alpha_dot, Lambda - gfLambda)))
Setup the handling of the evolution of internal variables. The evolution equations are a system of nonlinear equations per quadrature point. NewtonCF solves such equations with given starting points, tolerance and maximum interations. On failure, it returns NaN.
[11]:
# Parameters for evolution
tol = 1e-10
maxiter = 20
evol_tol = Parameter(tol)
eqs = evolution_eqs(strain(gfu), p, alpha, Lambda, evol_tol).Compile(realcompile=False)
evolution = NewtonCF(eqs, gfsol.components[1:], tol=tol, maxiter=maxiter)
The “trial state” that indicates whether there will be a plastic evolution or is computed via an exact (interpolating) projection.
[12]:
# We employ a BilinearForm to evaluate the weak form of the projection problem.
# Thus, we will only call L_trial.Apply
L_trial = BilinearForm(fes_hist, nonassemble=True)
p_trial, alpha_trial = fes_hist.TrialFunction()
p_test, alpha_test = fes_hist.TestFunction()
Psi_trial = Psi(strain(gfu), p_trial, alpha_trial)
L_trial += InnerProduct(-Psi_trial.Diff(p_trial), p_test) * irs_dx
L_trial += InnerProduct(-Psi_trial.Diff(alpha_trial), alpha_test) * irs_dx
# The mass matrix of the projection.
M_trial = BilinearForm(fes_hist, symmetric=True, diagonal=True)
for _trial, _test in zip(*fes_hist.TnT()):
M_trial += InnerProduct(_trial, _test).Compile() * irs_dx
M_trial.Assemble()
M_trial_inv = M_trial.mat.Inverse()
# RHS vector
vec_trial = gftrial.vec.CreateVector()
def compute_trial_state():
L_trial.Apply(gfhist.vec, vec_trial)
gftrial.vec.data = M_trial_inv * vec_trial
The actual solution of the equations happens during an again exact (interpolating) projection that sets gfint.
[13]:
# A Linearform for which the evluation triggers the solution of the evolution equations
testfunc = CoefficientFunction(tuple(fes_int.TestFunction()))
L_evol = LinearForm(fes_int)
L_evol += InnerProduct(CacheCF(evolution), testfunc) * irs_dx
# Mass-Matrix for the internal variable space as generic interpolation tool
M_evol = BilinearForm(fes_int, symmetric=True, diagonal=True)
for _trial, _test in zip(*fes_int.TnT()):
M_evol += InnerProduct(_trial, _test).Compile() * irs_dx
M_evol.Assemble()
M_evol_inv = M_evol.mat.Inverse()
def evolve_internal():
compute_trial_state()
L_evol.Assemble()
gfint.vec.data = M_evol_inv * L_evol.vec
for i, comp in enumerate(gfint.components):
gfsol.components[i+1].vec.data = gfint.components[i].vec
What remains is to setup the BilinearForm that describes the global boundary value problem.
[14]:
testfunc = fes.TestFunction()
u_test = testfunc[0]
test_int = CoefficientFunction(tuple(testfunc[1:]))
# Create a strain variable for partial differentiation
strain_var = strain(u)
internal_virtual_work_density = InnerProduct(Psi(strain_var, p, alpha).Diff(strain_var), strain(u_test))
# This "constraint" ensures the consistent linearization
evolution_constraint = InnerProduct(evolution_eqs(strain_var, p, alpha, Lambda, evol_tol), test_int)
loading_potential = -u[1] * force * loadfactor * ds("top")
realcompile = False
# Possible optimization: use static condensation for p
a = BilinearForm(fes, symmetric=False)
a += internal_virtual_work_density.Compile(realcompile=realcompile) * irs_dx
a += evolution_constraint.Compile(realcompile=realcompile) * irs_dx
a += Variation(loading_potential).Compile(realcompile=realcompile)
The nonlinear PDE is solved with a Newton-Raphson scheme. Thereby, in each Newton iteration we solve a linearized PDE for the displacements and a nonlinear evolution problem at each quadrature point via ``NewtonCF``
In order to run the problem, we need to define vectors that store respective solution data, define loadsteps and setup a loop that iterates through the latter.
[15]:
# A set of functions/vectors for keeping data
w = gfsol.vec.CreateVector()
rhs = gfsol.vec.CreateVector()
# Postprocess internal states
fd = LinearForm(InnerProduct(gfp, qd)*irs_dx)
fda = LinearForm(InnerProduct(gfalpha_k, qda)*irs_dx)
# Load steps (chosen based on experience)
loadsteps = [0.1,0.3,0.5,0.7,0.8,0.9,0.95,1.0]
# Set solution to zero initially
gfsol.vec[:] = 0
gfint.vec[:] = 0
gfhist.vec[:] = 0
gftrial.vec[:] = 0
# Iterate through load steps
for ls in loadsteps:
loadfactor.Set(ls)
# Update old solution at time t = t_k
store_internal()
with TaskManager():
for i in range(20):
a.AssembleLinearization(gfsol.vec)
a.Apply(gfsol.vec, rhs)
# Static condensation could be employed to speedup solving
w.data = a.mat.Inverse(freedofs=fes.FreeDofs(False), inverse="umfpack") * rhs
gfsol.vec.data -= w
# Solve the evolution equastion (quadrature-point-wise)
evolve_internal()
if np.isnan(np.max(np.abs(gfp.vec.FV().NumPy()))):
raise Exception("Evolution solver failed")
# NOTE: this problem is not a minimization problem. Therefore, we cannot just sum up (W * R).
err = np.max(np.abs(w.FV().NumPy() * rhs.FV().NumPy()))
print("step ", i, "err = ", err)
# Check convergence
if err < 1e-6: break
print("force = ", ls * force, ", uy_A =", gfu(node_A)[1], ", ux_B =", gfu(node_B)[0],\
", int u2 =", Integrate(gfu[1] * ds("top"),mesh))
called base class apply, type = N5ngfem29SymMatrixDifferentialOperatorE
called base class apply, type = N5ngfem29SymMatrixDifferentialOperatorE
called base class apply, type = N5ngfem29SymMatrixDifferentialOperatorE
called base class apply trans, type = called base class apply trans, type = called base class apply trans, type = N5ngfem29SymMatrixDifferentialOperatorEN5ngfem29SymMatrixDifferentialOperatorE
N5ngfem29SymMatrixDifferentialOperatorEcalled base class apply trans, type = N5ngfem29SymMatrixDifferentialOperatorE
step 0 err = 4.712439334558208
step 1 err = 6.245282381131554e-28
force = 45.0 , uy_A = 0.020951379066401273 , ux_B = 0.0076757893577595404 , int u2 = 2.040344506156137
step 0 err = 18.84975733823284
step 1 err = 1.728968421404716e-27
force = 135.0 , uy_A = 0.06285413719920381 , ux_B = 0.023027368073278622 , int u2 = 6.121033518468413
step 0 err = 18.849757338232656
step 1 err = 7.12307056823448e-07
force = 225.0 , uy_A = 0.1047574352802679 , ux_B = 0.03837870025764183 , int u2 = 10.201748629630139
step 0 err = 18.850532723835535
step 1 err = 3.547156586217235e-06
step 2 err = 4.487561979041226e-17
force = 315.0 , uy_A = 0.14666332319174966 , ux_B = 0.05372884322602869 , int u2 = 14.282591447766409
step 0 err = 4.7129242441516395
step 1 err = 4.1926122760202725e-07
force = 360.0 , uy_A = 0.16761767283367765 , ux_B = 0.06140327364897293 , int u2 = 16.323081973835155
step 0 err = 4.713583054844594
step 1 err = 4.936185159622969e-06
step 2 err = 1.7895139415362223e-18
force = 405.0 , uy_A = 0.18857831811672796 , ux_B = 0.06907507603344919 , int u2 = 18.363849243671936
step 0 err = 1.1789000392728395
step 1 err = 8.764763248513859e-07
force = 427.5 , uy_A = 0.1990627297411155 , ux_B = 0.07290934219032726 , int u2 = 19.384411829814205
step 0 err = 1.179300535818519
step 1 err = 1.155307828057331e-06
step 2 err = 2.0308123074112682e-19
force = 450.0 , uy_A = 0.2095516123099226 , ux_B = 0.0767418892060738 , int u2 = 20.405179684012705
[16]:
fd.Assemble()
pdraw.vec.data = invad * fd.vec
Draw(Norm(pdraw), mesh, "p")
Draw(gfu,mesh, "u")
if H.Get() > 1e-16:
fda.Assemble()
adraw.vec.data = invada * fda.vec
Draw(Norm(adraw), mesh, "alpha")
[17]:
# confirm that the trace vanishes
p1 = np.array(pdraw(mesh(90, 100, 0))).reshape((3,3))
np.einsum('ij,ij', p1, np.eye(3))
[17]:
2.879912020664621e-20
[ ]:
Implementation based on constrained minimization¶
We again start with the (coupled) function spaces. This time, we won’t need a TrialFunction “variable” for \(\alpha\) but only a GridFunction to store its value. Also, there is no need for \(\Lambda\) in this formulation.
Moreover, we use a deviatoric space for \(p\).
[18]:
fes_p = MatrixValued(fes_ir, symmetric=True, deviatoric=True, dim=3)
fes = fes_u * fes_p
u, p = fes.TrialFunction()
# GridFunction for solution
gfsol = GridFunction(fes)
gfu, gfp = gfsol.components
# Save previous solution
gfsol_k = GridFunction(fes)
gfu_k, gfp_k = gfsol_k.components
# alpha_0 = 0
alpha_k = GridFunction(fes_ir)
The energy desnity \(\Psi\) and the dissipation potential \(\Phi\)
[19]:
# The full energy density (incl. hardening)
def Psi(p, strain, p_k, alpha_k):
strain_energy = elastic_strain_energy(strain-p, Emod, nu)
alpha = alpha_k + IfPos(H - 1e-16, sqrt(2/3) * sigma_Y * H * tensor_norm(p - p_k, pert=norm_pert), 0)
hardening = 1/2 * alpha**2
return strain_energy + hardening
# Dissipation potential: For rate-independent plasticity (homogeneous of degree 1 in the rate of p)
# this is the dissipation! Also, because of rate independence, the parameter Delta_t is not needed for
# the evolution.
def Phi(p, p_k):
return sqrt(2/3) * sigma_Y * tensor_norm(p - p_k, pert=norm_pert)
Finally, we have all ingredients for the variational formulation:
Find for given \((u^k,p^k)\) a solution \((u^{k+1},p^{k+1})\) such that \begin{align*} \Psi(\epsilon(u^{k+1}), p^{k+1},p^{k}) + \Phi(p^{k+1},p^k)-\int f^{k+1}\cdot u^{k+1}\,dx \rightarrow \min! \end{align*}
[20]:
realcompile = False
# Possible optimization: use static condensation for p
a = BilinearForm(fes, symmetric=True)
a += Variation( Psi(p, strain(u), gfp_k, alpha_k) * irs_dx ).Compile(realcompile=realcompile)
a += Variation( Phi(p, gfp_k) * irs_dx ).Compile(realcompile=realcompile)
a += Variation( -u[1] * force * loadfactor * ds("top") ).Compile(realcompile=realcompile)
Since Phi(p, gfp_k) is included above, we automatically obtain an “algorithmically consistent linearization”, correspondig to the local evolution
[21]:
evolution_objective = (Psi(p, strain(gfu), gfp_k, alpha_k) + Phi(p,gfp_k)).Compile(realcompile=False)
which will be used to solve for p in a nested way, ie. for each “guess” for u^{k+1}.
In order to run the problem, we need to define vectors that store respective solution data, define loadsteps and setup a loop that iterates through the latter.
The nonlinear PDE is solved with a Newton-Raphson scheme. Thereby, in each Newton iteration we solve a linearized PDE for the displacements and a nonlinear evolution problem at each quadrature point via ``MinimizationCF``
[22]:
# a set of vectors for keeping data
vec_k = gfsol.vec.CreateVector()
w = gfsol.vec.CreateVector()
rhs = gfsol.vec.CreateVector()
# postprocess internal states
fd = LinearForm(InnerProduct(gfp, qd)*irs_dx)
fda = LinearForm(InnerProduct(gfalpha_k, qda)*irs_dx)
# load steps (chosen based on experience)
loadsteps = [0.1,0.3,0.5,0.7,0.8,0.9,0.95,1.0]
# set solution to zero initially
gfsol.vec[:] = 0
alpha_k.vec[:] = 0
# evolution params (default values of MinimizationCF at the time of writing)
tol = 1e-6
maxiter = 20
# iterate through load steps
for ls in loadsteps:
loadfactor.Set(ls)
# update old solution at time t = t_k
gfsol_k.vec.data = gfsol.vec
with TaskManager():
for i in range(20):
energy_k = a.Energy(gfsol.vec)
#print ("energy(t_k) = ", energy_k)
a.AssembleLinearization(gfsol.vec)
a.Apply(gfsol.vec, rhs)
# static condensation could be employed to speedup solving
w.data = a.mat.Inverse(freedofs=fes.FreeDofs(False), inverse="sparsecholesky") * rhs
# linesearch ( with damping)
vec_k.data = gfsol.vec
scale = 1
while scale > 1e-7:
gfsol.vec.data -= scale * w
energy1 = a.Energy(gfsol.vec)
# Here we solve the evolution equation as each quadrature point through evaluation of
# the MinimizationCF. This function takes the "objective", the initial guess (most
# likely the previous solution), tolerances and the maximum number of iterations.
gfp.Interpolate(MinimizationCF(evolution_objective, gfp, tol=tol, maxiter=maxiter))
energy2 = a.Energy(gfsol.vec)
if energy2 < energy_k + 1e-12:
break
scale *= 0.5
gfsol.vec.data = vec_k
err = sqrt(InnerProduct(w, rhs))
print("step ", i, "err = ", err)
# check convergence
if err < 1e-6: break
alpha_k.Interpolate(alpha_k + sqrt(2/3) * sigma_Y * H * sqrt(InnerProduct(gfp - gfp_k, gfp - gfp_k)))
print("force = ", ls * force, ", uy_A =", gfu(node_A)[1], ", ux_B =", gfu(node_B)[0],\
", int u2 =", Integrate(gfu[1] * ds("top"),mesh))
step 0 err = 9.582056225662322
step 1 err = 1.5780944678840155e-07
force = 45.0 , uy_A = 0.02095144754202526 , ux_B = 0.007675845674702156 , int u2 = 2.040351173243636
step 0 err = 19.164151419493592
step 1 err = 4.6989951225375645e-06
step 2 err = 5.430056955037999e-14
force = 135.0 , uy_A = 0.06285441912483467 , ux_B = 0.023027598311892786 , int u2 = 6.1210609217462
step 0 err = 19.16427609990064
step 1 err = 0.0022109828051319487
step 2 err = 5.613339911029854e-09
force = 225.0 , uy_A = 0.10475810223204271 , ux_B = 0.03837924063529988 , int u2 = 10.201813369883073
step 0 err = 19.164462947492712
step 1 err = 0.008483893342100467
step 2 err = 6.256877107602741e-08
force = 315.0 , uy_A = 0.1466646154772785 , ux_B = 0.05372987454255943 , int u2 = 14.282716517935217
step 0 err = 9.583036183681568
step 1 err = 0.006165763008128386
step 2 err = 1.4664681909784846e-08
force = 360.0 , uy_A = 0.16761977271708897 , ux_B = 0.06140492220520205 , int u2 = 16.32328447096966
step 0 err = 9.583299064510136
step 1 err = 0.010614623935780314
step 2 err = 4.210088375059896e-08
force = 405.0 , uy_A = 0.1885813183739071 , ux_B = 0.06907759651425152 , int u2 = 18.364147590546313
step 0 err = 4.793816022500711
step 1 err = 0.007956947826377103
step 2 err = 2.45808619655413e-08
force = 427.5 , uy_A = 0.19906676318716343 , ux_B = 0.07291289320540822 , int u2 = 19.38482281942709
step 0 err = 4.794109890790945
step 1 err = 0.010136194946113881
step 2 err = 3.9074072898010314e-08
force = 450.0 , uy_A = 0.20955665403136225 , ux_B = 0.07674678471387063 , int u2 = 20.405718414146126
[23]:
fd.Assemble()
pdraw.vec.data = invad * fd.vec
Draw(Norm(pdraw), mesh, "p")
Draw(gfu,mesh, "u")
if H.Get() > 1e-16:
fda.Assemble()
adraw.vec.data = invada * fda.vec
Draw(Norm(adraw), mesh, "alpha")
[24]:
# confirm that the trace vanishes
p1 = np.array(pdraw(mesh(90, 100, 0))).reshape((3,3))
np.einsum('ij,ij', p1, np.eye(3))
[24]:
1.6940658945086007e-21
[ ]: