Dear Community.
In my previous models I used im_data to calculate new or store previous quantities at the im_nodes (internal or boundary ones) for assembly purposes or visualization,
but I’ve never considered before the usage of add_im_variable() or add_internal_im_variable().
In which cases it is advisable to introduce im variables into a FEM model.
Thank you.
Lorenzo
these are so called internal variables and can optionally be condensed out from the overall system. If you define them as “internal”, then GetFEM will solve for those first and update the stiffness matrix of the remaining variables accordingly so that the linearization is consistent.
Basically im_variables allow you to write equations that are defined at each Gauss point, with something like
md.add_im_variable("var", mimd)
md.add_nonlinear_term(mim, "(.... some... equation)*Test_var")
Hello Kostas,
Are there examples using im_variables and im_internal_variables?
I’m a little bit confused about defining an expression at each Gauss points using im_internal_variables instead of just defining a macro with the expression and putting it in an assembly term.
Thank you,
Lorenzo
start with cheking these two tests in C++:
https://git.savannah.nongnu.org/cgit/getfem.git/tree/tests/test_internal_variables.cc
and
https://git.savannah.nongnu.org/cgit/getfem.git/tree/tests/test_condensation.cc
the use of these functions in python is totally equivalent
Typically one would use an internal variable for the plastic multiplier in a model for plasticity with nonlinear hardening. Then the yield equation is linearized and solved at the Gauss point. I can find some example if you are interested in further details.
Thank you Kostas.
I will check the provided examples.
Since I’m very interested in understanding the topic and potential applications, if you can provide further details and examples it would be highly appreciated.
Have a nice day.
Lorenzo
a bit experimental script, not working 100% correct but you can consider for inspiration
import getfem as gf
import numpy as np
import os, sys, subprocess
import time
gf.util_trace_level(1)
# Input data
L = 2*26.667 # block length
H = 2*6.413 # block height
dH = 0.018*H # height reduction at the center of the block
N_L = 20 # number of elements in block length direction
N_H = 10 # number of elements in block height direction
E = 210e3 # Young's modulus
nu = 0.3 # Poisson's ratio
#plast_law = "Ramberg Osgood"
#pl_sigma_0 = 0.45e3 # Initial yield stress
#pl_n = 4.
#pl_alpha = 3./7.
#plast_law = "Linear hardening"
#pl_sigma_0 = 0.45e3 # Initial yield stress
#pl_H = 0.12924e3 # Hardening/softening tangent modulus
plast_law = "Linear and exponential hardening"
pl_sigma_0 = 0.45e3 # Initial yield stress
pl_sigma_inf = 0.715e3
pl_H = 0.12924e3 # Hardening/softening tangent modulus
pl_delta = 16.93
steps = 200
numerical_continuation = False
dirichlet_str = "[{0}*{1}*x,0]".format("{0}",14./L) # x stretching with 0.2 the stretching ratio
#------------------------------------
geotrans = "GT_QK(2,2)" # geometric transformation
disp_fem_order = 2 # displacements finite element order
press_fem_order = -1 # pressure finite element order
#integration_degree = 3 # 4 gauss points per quad
integration_degree = 5 # 9 gauss points per quad
#integration_degree = 7 # 16 gauss points per quad
#integration_degree = 9 # 25 gauss points per quad
integration_degree_press = 3
rigid_grip = False
#------------------------------------
tee = subprocess.Popen(["tee", "./finite_strain_plasticity_dksi.log"],
stdin=subprocess.PIPE)
sys.stdout.flush()
os.dup2(tee.stdin.fileno(), sys.stdout.fileno())
sys.stderr.flush()
os.dup2(tee.stdin.fileno(), sys.stderr.fileno())
mixed_up = (press_fem_order >= 0)
# auxiliary constants
L_RG = 3
R_RG = 4
B_RG = 5
T_RG = 6
#DIR_RG = 7
xmin = 0.
ymin = 0.
dx = L/2
dy = H/2
m = gf.Mesh("import", "structured",
f"GT='{geotrans}';ORG=[{xmin},{ymin}];SIZES=[{dx},{dy}];NSUBDIV=[{N_L},{N_H}]")
N = m.dim()
m.set_region(L_RG, m.outer_faces_with_direction([-1.,0.], np.pi/36))
m.set_region(R_RG, m.outer_faces_with_direction([1.,0.], np.pi/36))
m.set_region(B_RG, m.outer_faces_with_direction([0.,-1.], np.pi/36))
m.set_region(T_RG, m.outer_faces_with_direction([0.,1.], np.pi/36))
#m.merge_region(DIR_RG, L_RG)
#m.merge_region(DIR_RG, R_RG)
if dH > 0:
pts = m.pts()
dy = 0.2 # fine element size ratio w.r.t. uniform mesh size
y0 = 0.04 # ratio of the finer meshed area w.r.t. the total length
x0 = y0/dy
b2 = (1-dy)/(x0-1)**2
b1 = dy - 2*b2*x0
b0 = 1 - b1 -b2
for i in range(pts.shape[1]):
x = pts[0,i]
y = pts[1,i]
t = 2*abs(x)/L
if t < x0:
x *= dy;
else:
x = (b0 + b1*t + b2*t**2) * np.sign(x)*L/2
pts[0,i] = x
pts[1,i] -= (y*dH)/(2*H) * (1 + np.cos(2.*np.pi*x/L))
m.set_pts(pts)
m.export_to_vtu("mesh.vtu")
# FEM
mfN = gf.MeshFem(m, N)
mfN.set_classical_fem(disp_fem_order)
#if disp_fem_order == 2:
# mfN.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(2)"))
#else:
# mfN.set_classical_fem(disp_fem_order)
xdofsL = mfN.basic_dof_on_region(L_RG)[0::N]
ydofsB = mfN.basic_dof_on_region(B_RG)[1::N]
kept_dofs = np.setdiff1d(np.arange(mfN.nb_basic_dof()),
np.union1d(xdofsL,ydofsB))
mfu = gf.MeshFem("partial", mfN, kept_dofs)
if mixed_up:
mfp = gf.MeshFem(m, 1)
mfp.set_classical_fem(press_fem_order)
if not rigid_grip:
mfdir = gf.MeshFem(m, 1)
if disp_fem_order == 2:
mfdir.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(2)"))
else:
mfdir.set_classical_fem(disp_fem_order)
else:
mfdir = mfN
mfout1 = gf.MeshFem(m)
mfout1.set_classical_discontinuous_fem(disp_fem_order-1)
mfout2 = gf.MeshFem(m)
mfout2.set_classical_discontinuous_fem(disp_fem_order)
mim = gf.MeshIm(m, integration_degree)
mimd1 = gf.MeshImData(mim, -1)
mimd4 = gf.MeshImData(mim, -1, [4])
if mixed_up:
mimpress = gf.MeshIm(m, integration_degree_press)
# Model
md = gf.Model("real")
md.add_fem_variable("u", mfu) # displacements
if mixed_up:
md.add_fem_variable("p", mfp) # pressure
md.add_internal_im_variable("dksi", mimd1) # plastic multiplier increase in the current step
md.add_im_data("ksi0", mimd1) # plastic multiplier value at previous time step
md.add_im_data("invCp0vec", mimd4) # Components 11, 22, 33 and 12 of the plastic part
md.set_variable("invCp0vec", # of the inverse right Cauchy Green tensor at the
np.tile([1,1,1,0], # previous step.
mimd4.nbpts())) # Component 21 is equal to 12 and is omitted.
md.add_initialized_data("gamma", [0.]) # Current load level (in percent)
md.add_initialized_data("K", E/(3.*(1.-2.*nu))) # Bulk modulus
md.add_initialized_data("mu", E/(2*(1+nu))) # Shear modulus
md.add_macro("F", "Id(2)+Grad_u")
md.add_macro("F3d", "Id(3)+[0,0,0;0,0,0;0,0,1/X(2)]*u(2)+[1,0;0,1;0,0]*Grad_u*[1,0,0;0,1,0]")
md.add_macro("J", "Det(F)*(1+u(2)/X(2))")
#md.add_macro("J", "Det(F3d)")
md.add_macro("invCp0", "[[[1,0,0],[0,0,0],[0,0,0]],"+\
" [[0,0,0],[0,1,0],[0,0,0]],"+\
" [[0,0,0],[0,0,0],[0,0,1]],"+\
" [[0,1,0],[1,0,0],[0,0,0]]].invCp0vec") #Vec4ToMat3x3
md.add_macro("devlogbetr", "Deviator(Logm(F3d*invCp0*F3d'))")
md.add_macro("devlogbe", "(1-2*dksi)*devlogbetr")
md.add_macro("tauD2d", "mu*pow(J,-2/3)*"+\
"([[[[1,0],[0,0]] , [[0,1],[0,0]] , [[0,0],[0,0]]],"+\
" [[[0,0],[1,0]] , [[0,0],[0,1]] , [[0,0],[0,0]]],"+\
" [[[0,0],[0,0]] , [[0,0],[0,0]] , [[0,0],[0,0]]]]:devlogbe)")
if mixed_up:
md.add_macro("tauH", "-p*J")
else:
md.add_macro("tauH", "K*log(J)")
md.add_macro("tauD33", "mu*devlogbe(3,3)")
md.add_nonlinear_term\
(mim, "2*pi*X(2)*(((tauH*Id(2)+tauD2d)*Inv(F')):Grad_Test_u+(tauH+tauD33)/(X(2)+u(2))*Test_u(2))")
if mixed_up:
md.add_nonlinear_term(mimpress, "(p/K+log(J)/J)*Test_p")
if plast_law == "Ramberg Osgood":
yield_lim_str = "{coef}*pow(ksi+1e-12,{exponent})"\
.format(coef=pow(E/pl_alpha * pow(2./3., (pl_n+1)/2.) * pow(pl_sigma_0, pl_n-1),
1./pl_n),
exponent=1./pl_n)
elif plast_law == "Linear hardening":
yield_lim_str = "{A}+{B}*ksi".format(A=np.sqrt(2./3.)*pl_sigma_0, B=2./3.*pl_H)
elif plast_law == "Linear and exponential hardening":
yield_lim_str = "{A}+{B}*(1-exp(-{C}*ksi))+{D}*ksi"\
.format(A=np.sqrt(2./3.)*pl_sigma_0, B=np.sqrt(2./3.)*(pl_sigma_inf-pl_sigma_0),
C=np.sqrt(2./3.)*pl_delta, D=2./3.*pl_H)
else:
raise NameError(plast_law)
md.add_macro("ksi", "ksi0+dksi*Norm(devlogbetr)")
md.add_macro("sigmaD", "mu*pow(J,-5/3)*devlogbe")
md.add_macro("sigmaDtr", "mu*pow(J,-5/3)*devlogbetr")
md.add_nonlinear_term\
(mim, f"(Norm(sigmaD)-min({yield_lim_str},Norm(sigmaDtr)+1e-12*dksi))*Test_dksi")
if rigid_grip:
md.add_fem_data("dirichlet_data", mfN);
ibdir = md.add_Dirichlet_condition_with_multipliers(mim, "u", mfdir, R_RG, "dirichlet_data")
else:
dirichlet_str = dirichlet_str[dirichlet_str.find("[")+1:dirichlet_str.find(",")]
md.add_fem_data("dirichlet_data", mfdir);
ibdir = md.add_normal_Dirichlet_condition_with_multipliers(mim, "u", mfdir, R_RG, "dirichlet_data")
# ibdir = md.add_Dirichlet_condition_with_penalization(mim, "u", 1e6, DIR_RG, "dirichlet_data")
mfdirdata = md.mesh_fem_of_variable("dirichlet_data")
# Numerical continuation
if numerical_continuation:
md.add_fem_data("dirichlet_data_init", mfdirdata);
md.add_fem_data("dirichlet_data_final", mfdirdata);
md.set_variable("dirichlet_data_init",
mfdirdata.eval(dirichlet_str.format(0),globals(),locals()).flatten("F"))
md.set_variable("dirichlet_data_final",
mfdirdata.eval(dirichlet_str.format(1),globals(),locals()).flatten("F"))
scfac = 1e0/md.interval_of_variable("u")[1]
S = gf.ContStruct\
(md, "gamma", "dirichlet_data_init", "dirichlet_data_final",
"dirichlet_data", scfac,
"variable_name", "u",
"h_init", 1e-1,
"h_max", 1e0,
"h_min", 1e-5,
"h_dec", 1./np.sqrt(10.),
"h_inc", np.sqrt(10.),
"lsolver", "mumps",
"max_iter", 12,
"thr_iter", 6,
"max_res", 1e-8,
"max_diff", 1e-8,
"min_cos", 0.995,
"singularities", 0, "very_noisy") #, "non-smooth"
print(f"Displacement dofs: {mfu.nbdof()}\nTotal model dofs: {md.nbdof()}")
dirmultname = md.mult_varname_Dirichlet(ibdir)
mfdirmult = md.mesh_fem_of_variable(dirmultname)
sl_shrinked = gf.Slice(("explode", 1.), m, 4)
sl_tracked = gf.Slice("points", m, [[0],[0]])
IU = md.interval_of_variable("u")
IDIR = md.interval_of_variable(dirmultname)
IU = range(IU[0],IU[0]+IU[1])
IDIR = range(IDIR[0],IDIR[0]+IDIR[1])
if mixed_up:
IP = md.interval_of_variable("p")
IP = range(IP[0],IP[0]+IP[1])
md.add_macro("invCp", "(Inv(F3d)*Expm(-dksi*devlogbetr)*(F3d))*invCp0"\
"*(Inv(F3d)*Expm(-dksi*devlogbetr)*(F3d))'")
md.add_macro("tau", "tauH*Id(3)+mu*pow(J,-2/3)*devlogbe")
starttime_overall = time.process_time()
with open("./finite_strain_plasticity_dksi_axisym_forces.dat", "w") as f:
for step in range(steps+1):
if md.variable("gamma") >= 1. or md.variable("gamma") < 0. or\
os.path.exists("./abort"):
break # simulation completed (or returned to zero load)
print(f"SOLVING STEP {step}")
starttime = time.process_time()
if step == 0 or not numerical_continuation:
loadfactor = float(step)/float(steps)
md.set_variable("gamma", [loadfactor])
md.set_variable("dirichlet_data",
mfdirdata.eval(dirichlet_str.format(loadfactor), globals(), locals()).flatten("F"))
md.solve("noisy", "max_iter", 25, "max_res", 1e-10,
"lsearch", "simplest", "alpha max ratio", 10, "alpha min", 0.2, "alpha mult", 0.6,
"alpha threshold res", 1e9)
if numerical_continuation:
init_dir = 1.
gamma = 0.
x = md.from_variables()
[tx, tgamma, h] = S.init_Moore_Penrose_continuation(x, gamma, init_dir)
hpre = h
limit_points = []
gammas = [0]
progress = [0]
else:
x = md.from_variables()
[x, gamma, tx, tgamma, h, hpre, sing_label] = \
S.Moore_Penrose_continuation(x, gamma, tx, tgamma, h)
if h == 0:
break # Continuation has failed
elif sing_label == "limit point":
limit_points.append(step)
md.to_variables(x)
print(f"STEP {step} COMPLETED IN {time.process_time()-starttime:f} SEC")
output = (mfout1, md.local_projection(mim, "sqrt(1.5)*Norm(sigmaD)", mfout1), "Von Mises Stress")
for i,j in [[1,1],[2,2],[3,3],[1,2],[2,1]]:
output += (mfout1, md.local_projection(mim, f"tau({i},{j})", mfout1),
f"Kirchhoff stress {i},{j}")
output += (mfout1, md.local_projection(mim, f"tau({i},{j})/J", mfout1),
f"Cauchy Stress {i},{j}")
output += (mfout2, md.local_projection(mim, f"devlogbetr({i},{j})", mfout2),
f"devlogbetr {i},{j}")
for i,j in [[1,1],[2,2],[1,2],[2,1]]:
output += (mfout2, md.local_projection(mim, f"Logm(Left_Cauchy_Green(F3d))({i},{j})", mfout2),
f"Log Left Cauchy Green {i},{j}")
output += (mfout2, md.local_projection(mim, "Norm(devlogbetr)", mfout2), "normdevbe")
RHS = md.rhs()
output += (mfu, RHS[IU], "residual u",
mfdirmult, RHS[IDIR], "residual dirichlet")
if mixed_up:
output = (mfp, RHS[IP], "residual p") + output
if numerical_continuation:
output += (mfu, tx[IU], "tangent u",
mfdirmult, tx[IDIR], "tangent dirichlet")
if mixed_up:
output += (mfp, tx[IP], "tangent p") + output
if mixed_up:
output = (mfp, md.variable("p"), "Pressure") + output
dksi = md.local_projection(mim, "dksi*Norm(devlogbetr)", mfout2)
if numerical_continuation:
output = (mfout2, dksi/hpre, "dksi over h") + output
output += (mfout2, dksi, "dksi")
output = (mfu, md.variable("u"), "Displacements",
mfout2, md.local_projection(mim, "dksi", mfout2), "dksi over normdevlogbetr",
mfdirmult, -md.variable(dirmultname), "Reaction stresses") + output
# update ksi0 after all stress calculations are done with the old ksi0 plus dksi
cc = 0.99
md.set_variable("dksi", cc*md.variable("dksi"))
md.set_variable("ksi0", md.interpolation("ksi", mimd1, -1))
md.set_variable("invCp0vec",
md.interpolation("[[[1,0,0,0] ,[0,0,0,0.5],[0,0,0,0]],"+\
" [[0,0,0,0.5],[0,1,0,0] ,[0,0,0,0]],"+\
" [[0,0,0,0] ,[0,0,0,0] ,[0,0,1,0]]]:invCp", mimd4, -1))
md.set_variable("dksi", (1./cc-1.)*md.variable("dksi"))
# output of updated state variables
output += (mfout2, md.local_projection(mim, "ksi", mfout2), "ksi")
for i,j in [[1,1],[2,2],[3,3],[1,2]]:
output += (mfout2, md.local_projection(mim, f"Inv(invCp0)({i},{j})", mfout2),
"Cp {0}{1}".format(i,j))
for i in range(4):
output += (mfout2, md.local_projection(mim, f"invCp0vec({i+1})", mfout2),
"invCp " + ['11','22','33','12'][i])
output += (mfout2, md.local_projection(mim, "Det(Inv(invCp0))", mfout2), "detCp")
sl_shrinked.export_to_vtu(f"finite_strain_plasticity_dksi_sl{step}.vtu", *output)
mfout2.export_to_vtu(f"finite_strain_plasticity_dksi_{step}.vtu", *output)
if numerical_continuation and step > 0:
progress.append(progress[-1]+hpre)
gammas.append(md.variable("gamma"))
elif step > 0:
progress.append(float(step)/float(steps))
ksi_center = gf.compute_interpolate_on(mfout2, md.local_projection(mim, "ksi0", mfout2),
sl_tracked)[0]
DIRMULT = -md.variable(dirmultname)
DIRMULT = np.reshape(DIRMULT,[1,DIRMULT.size])
dfR = gf.asm_boundary_source(R_RG, mim, mfdir, mfdirmult, DIRMULT)
if not rigid_grip:
f.write((f"step={step} gamma={md.variable("gamma")[0]:e} fR=({dfR.sum():e},0) "
f"sR=({dfR.sum()/H:e},0) progress={progress[-1]:e} ksi_center={ksi_center:e}\n"))
else:
f.write(("step=%i gamma=%e fR=(%e,%e) sR=(%e,%e) progress=%e ksi_center=%e\n") %
(step, md.variable("gamma"),
dfR[0::N].sum(), dfR[1::N].sum(),
dfR[0::N].sum()/H, dfR[1::N].sum()/H,
progress[-1], ksi_center))
f.flush()
print(f"OVERALL SOLUTION TIME IS {time.process_time()-starttime_overall:f} SEC")
if numerical_continuation:
from subprocess import call
progress = gammas[-2]/progress[-2] * np.array(progress)
progress[-1] = gammas[-1]
for postfix in ["", "sl"]:
with open(f"./postproc{postfix}.pvd", "w") as f:
f.write('<?xml version="1.0"?>\n')
f.write('<VTKFile type="Collection" version="0.1"\n')
f.write(' byte_order="LittleEndian">\n')
f.write(' <Collection>\n')
for step in range(len(progress)):
f.write(f' <DataSet timestep="{progress[step]:f}" group="" part="0" ')
f.write(f'file="finite_strain_plasticity_dksi_{postfix}{step}.vtu"/>\n')
f.write(' </Collection>\n')
f.write('</VTKFile>')
Hello Kostas,
Thank you for indicating these examples. I tried to metabolize them but I’ve still some doubts.
The C++ guide of add_internal_im_variable() states that the im_variable will be “statically” condensed out during the linearization of the problem.
My questions are:
Thank you,
Lorenzo
The static condensation code, which I added a couple of years ago, lacks proper documentation documentation indeed.
I will try to describe the most important aspects here.
Imagine that u_R and u_I are respectively your “retained” and “internal” variables. Assume also that you are solving for a residulal vector
If you define u_I with add_im_variable, then the system R=0 is solved as normally. Internal variables in this case are treated as normal variables. The only slight difference is that the tangent matrix of the system will turn out to have many disconnected regions, because there is no direct coupling between the dofs in u_I, in contrast to the dofs for field variables contained in u_R. The tangent system in this case is:
with K_{II} being diagonal or block diagonal.
If you define u_I with add_internal_im_variable, then when you assemble your RHS and your tangent matrix, only the parts for the retained variable u_R will be assembled, let’s call them R'_R and K'_ {RR} respectively. The tangent system in this case is
with
and
In this case, the method model.tangent_matrix() will just give you K'_{RR} and model.rhs() will give you R'_R.
So if you are performing time discretization, you should check if this treatment is consistent with your time discretization. E.g. if your time discretization is reduced to fulfilling a system like R(u_R,u_I)=0 at every time step, then everything should be fine.