[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] [getfem-commits] branch master updated: Remove erroneou
From: |
Konstantinos Poulios |
Subject: |
[Getfem-commits] [getfem-commits] branch master updated: Remove erroneous J^-2/3 factor from finite strain plasticity examples |
Date: |
Sat, 16 Sep 2023 16:18:57 -0400 |
This is an automated email from the git hooks/post-receive script.
logari81 pushed a commit to branch master
in repository getfem.
The following commit(s) were added to refs/heads/master by this push:
new 1910cb6e Remove erroneous J^-2/3 factor from finite strain plasticity
examples
1910cb6e is described below
commit 1910cb6ee6f72ffe2f819bdb054097aa6f830898
Author: Konstantinos Poulios <logari81@gmail.com>
AuthorDate: Sat Sep 16 22:18:46 2023 +0200
Remove erroneous J^-2/3 factor from finite strain plasticity examples
+ coding style improvements
---
...ticity_fin_strain_lin_hardening_axisymmetric.py | 32 +++++------
...ticity_fin_strain_lin_hardening_plane_strain.py | 64 ++++++++--------------
...ty_finite_strain_linear_hardening_tension_3D.py | 34 ++++++------
3 files changed, 55 insertions(+), 75 deletions(-)
diff --git
a/contrib/continuum_mechanics/plasticity_fin_strain_lin_hardening_axisymmetric.py
b/contrib/continuum_mechanics/plasticity_fin_strain_lin_hardening_axisymmetric.py
index 7969b238..f24d1db5 100644
---
a/contrib/continuum_mechanics/plasticity_fin_strain_lin_hardening_axisymmetric.py
+++
b/contrib/continuum_mechanics/plasticity_fin_strain_lin_hardening_axisymmetric.py
@@ -2,7 +2,7 @@
# -*- coding: UTF8 -*-
# Python GetFEM interface
#
-# Copyright (C) 2020-2020 Konstantinos Poulios.
+# Copyright (C) 2020-2023 Konstantinos Poulios.
#
# This file is a part of GetFEM
#
@@ -68,7 +68,7 @@ resultspath = "./results"
if not os.path.exists(resultspath):
os.makedirs(resultspath)
-tee = subprocess.Popen(["tee", "%s/tension_axisymmetric.log" % resultspath],
+tee = subprocess.Popen(["tee", f"{resultspath}/tension_axisymmetric.log"],
stdin=subprocess.PIPE)
sys.stdout.flush()
os.dup2(tee.stdin.fileno(), sys.stdout.fileno())
@@ -86,8 +86,7 @@ ymin = 0.
dx = L/2
dy = H/2
mesh = gf.Mesh("import", "structured",
- "GT='%s';ORG=[%f,%f];SIZES=[%f,%f];NSUBDIV=[%i,%i]"
- % (geotrans, xmin, ymin, dx, dy, N_L, N_H))
+
f"GT='{geotrans}';ORG=[{xmin},{ymin}];SIZES=[{dx},{dy}];NSUBDIV=[{N_L},{N_H}]")
N = mesh.dim()
outer_faces = mesh.outer_faces()
@@ -119,7 +118,7 @@ if dH > 0:
pts[1,i] -= (y*dH)/(2*H) * (1 + np.cos(2.*np.pi*x/L))
mesh.set_pts(pts)
-mesh.export_to_vtu("%s/mesh.vtu" % resultspath)
+mesh.export_to_vtu(f"{resultspath}/mesh.vtu")
# FEM
mfN = gf.MeshFem(mesh, N)
@@ -171,12 +170,12 @@ md.add_macro("invCp0", "[[[1,0,0],[0,0,0],[0,0,0]],"+\
md.add_macro("devlogbetr", "Deviator(Logm(F3d*invCp0*F3d'))")
md.add_macro("Y0","{A}+{B}*gamma0".format(A=np.sqrt(2./3.)*pl_sigma_0,
B=2./3.*pl_H))
md.add_macro("ksi",
-
"(1-1/max(1,mu*pow(J,-5/3)*Norm(devlogbetr)/Y0))/(2+{B}/(mu*pow(J,-5/3)))"
+ "(1-1/max(1,mu/J*Norm(devlogbetr)/Y0))/(2+{B}/(mu/J))"
.format(B=2./3.*pl_H))
md.add_macro("gamma", "gamma0+ksi*Norm(devlogbetr)")
md.add_macro("devlogbe", "(1-2*ksi)*devlogbetr")
-md.add_macro("tauD2d", "mu*pow(J,-2/3)*[1,0,0;0,1,0]*devlogbe*[1,0;0,1;0,0]")
-md.add_macro("tauD33", "mu*pow(J,-2/3)*devlogbe(3,3)")
+md.add_macro("tauD2d", "mu*[1,0,0;0,1,0]*devlogbe*[1,0;0,1;0,0]")
+md.add_macro("tauD33", "mu*devlogbe(3,3)")
md.add_nonlinear_generic_assembly_brick\
(mim,
"2*pi*X(2)*(((tauH*Id(2)+tauD2d)*Inv(F')):Grad_Test_u+(tauH+tauD33)/(X(2)+u(2))*Test_u(2))")
@@ -190,22 +189,21 @@ md.add_macro("invCp",
"(Inv(F3d)*Expm(-ksi*devlogbetr)*(F3d))*invCp0"\
md.add_filtered_fem_variable("dirmult", mfmult, XP_RG)
md.add_nonlinear_generic_assembly_brick(mim, "2*pi*X(2)*(disp-u(1))*dirmult",
XP_RG)
-print("Model dofs: %i" % md.nbdof())
-print("Displacement fem dofs: %i" % mfu.nbdof())
+print(f"Model dofs: {md.nbdof()}\nDisplacement fem dofs: {mfu.nbdof()}")
print("Dirichlet mult dofs: %i" % md.mesh_fem_of_variable("dirmult").nbdof())
shutil.copyfile(os.path.abspath(sys.argv[0]),resultspath+"/"+sys.argv[0])
starttime_overall = time.process_time()
-with open("%s/tension_axisymmetric.dat" % resultspath, "w") as f1:
+with open(f"{resultspath}/tension_axisymmetric.dat", "w") as f1:
for step in range(steps_t+1):
md.set_variable("disp", disp*step/float(steps_t))
print('STEP %i: Solving with disp = %g' % (step, md.variable("disp")))
starttime = time.process_time()
- md.solve('noisy', 'max_iter', 25, 'max_res', 1e-10,
- 'lsearch', 'simplest', 'alpha max ratio', 100, 'alpha min',
0.2, 'alpha mult', 0.6,
- 'alpha threshold res', 1e9)
- print('STEP %i COMPLETED IN %f SEC' % (step,
time.process_time()-starttime))
+ md.solve("noisy", "max_iter", 25, "max_res", 1e-10,
+ "lsearch", "simplest", "alpha max ratio", 100, "alpha min",
0.2, "alpha mult", 0.6,
+ "alpha threshold res", 1e9)
+ print("STEP %i COMPLETED IN %f SEC" % (step,
time.process_time()-starttime))
F = gf.asm_generic(mim, 0, "2*pi*X(2)*dirmult", XP_RG, md)
print("Displacement %g, total force %g" % (md.variable("disp"), F))
@@ -226,7 +224,7 @@ with open("%s/tension_axisymmetric.dat" % resultspath, "w")
as f1:
mfu, md.variable("u"), "Displacements",
mfout2, md.interpolation("dirmult", mfout2, XP_RG), "Nominal
reaction traction",
mfout2, md.local_projection(mim, "gamma", mfout2), "plastic
strain")
- mfout2.export_to_vtu("%s/tension_axisymmetric_%i.vtu" % (resultspath,
step), *output)
+ mfout2.export_to_vtu(f"{resultspath}/tension_axisymmetric_{step}.vtu",
*output)
md.set_variable("gamma0", md.interpolation("gamma", mimd1, -1))
md.set_variable("invCp0vec",
@@ -234,5 +232,5 @@ with open("%s/tension_axisymmetric.dat" % resultspath, "w")
as f1:
" [[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))
-print('OVERALL SOLUTION TIME IS %f SEC' %
(time.process_time()-starttime_overall))
+print("OVERALL SOLUTION TIME IS %f SEC" %
(time.process_time()-starttime_overall))
diff --git
a/contrib/continuum_mechanics/plasticity_fin_strain_lin_hardening_plane_strain.py
b/contrib/continuum_mechanics/plasticity_fin_strain_lin_hardening_plane_strain.py
index 3b3644c2..20bc6c79 100644
---
a/contrib/continuum_mechanics/plasticity_fin_strain_lin_hardening_plane_strain.py
+++
b/contrib/continuum_mechanics/plasticity_fin_strain_lin_hardening_plane_strain.py
@@ -1,8 +1,8 @@
#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
+# -*- coding: UTF8 -*-
# Python GetFEM interface
#
-# Copyright (C) 2020-2020 Konstantinos Poulios.
+# Copyright (C) 2020-2023 Konstantinos Poulios.
#
# This file is a part of GetFEM
#
@@ -62,14 +62,12 @@ mult_fem_order = 2 # dirichlet multipliers finite element
order
#integration_degree = 3 # 4 gauss points per quad
integration_degree = 5 # 9 gauss points per quad
-export_results = False;
-
#------------------------------------
resultspath = "./results"
-if (export_results and not os.path.exists(resultspath)):
+if not os.path.exists(resultspath):
os.makedirs(resultspath)
-tee = subprocess.Popen(["tee", "%s/tension_plane_strain.log" % resultspath],
+tee = subprocess.Popen(["tee", f"{resultspath}/tension_plane_strain.log"],
stdin=subprocess.PIPE)
sys.stdout.flush()
os.dup2(tee.stdin.fileno(), sys.stdout.fileno())
@@ -87,8 +85,7 @@ ymin = 0.
dx = L/2
dy = H/2
mesh = gf.Mesh("import", "structured",
- "GT='%s';ORG=[%f,%f];SIZES=[%f,%f];NSUBDIV=[%i,%i]"
- % (geotrans, xmin, ymin, dx, dy, N_L, N_H))
+
f"GT='{geotrans}';ORG=[{xmin},{ymin}];SIZES=[{dx},{dy}];NSUBDIV=[{N_L},{N_H}]")
N = mesh.dim()
outer_faces = mesh.outer_faces()
@@ -119,9 +116,8 @@ if dH > 0:
pts[0,i] = x
pts[1,i] -= (y*dH)/(2*H) * (1 + np.cos(2.*np.pi*x/L))
mesh.set_pts(pts)
-
-if (export_results):
- mesh.export_to_vtu("%s/mesh.vtu" % resultspath)
+
+mesh.export_to_vtu(f"{resultspath}/mesh.vtu")
# FEM
mfN = gf.MeshFem(mesh, N)
@@ -173,17 +169,17 @@ md.add_macro("invCp0", "[[[1,0,0],[0,0,0],[0,0,0]],"+\
md.add_macro("devlogbetr", "Deviator(Logm(F3d*invCp0*F3d'))")
md.add_macro("Y0","{A}+{B}*gamma0".format(A=np.sqrt(2./3.)*pl_sigma_0,
B=2./3.*pl_H))
md.add_macro("ksi",
-
"(1-1/max(1,mu*pow(J,-5/3)*Norm(devlogbetr)/Y0))/(2+{B}/(mu*pow(J,-5/3)))"
+ "(1-1/max(1,mu/J*Norm(devlogbetr)/Y0))/(2+{B}/(mu/J))"
.format(B=2./3.*pl_H))
md.add_macro("gamma", "gamma0+ksi*Norm(devlogbetr)")
md.add_macro("devlogbe", "(1-2*ksi)*devlogbetr")
-md.add_macro("tauD2d", "mu*pow(J,-2/3)*[1,0,0;0,1,0]*devlogbe*[1,0;0,1;0,0]")
+md.add_macro("tauD2d", "mu*[1,0,0;0,1,0]*devlogbe*[1,0;0,1;0,0]")
md.add_nonlinear_generic_assembly_brick\
(mim, "((tauH*Id(2)+tauD2d)*Inv(F')):Grad_Test_u")
-md.add_macro("sigmaD", "(mu*pow(J,-5/3)*devlogbe)")
-md.add_macro("sigma", "tauH/J*Id(3)+mu*pow(J,-5/3)*devlogbe")
+md.add_macro("sigmaD", "(mu/J*devlogbe)")
+md.add_macro("sigma", "tauH/J*Id(3)+mu/J*devlogbe")
md.add_macro("invCp", "(Inv(F3d)*Expm(-ksi*devlogbetr)*(F3d))*invCp0"\
"*(Inv(F3d)*Expm(-ksi*devlogbetr)*(F3d))'")
@@ -191,28 +187,21 @@ md.add_macro("invCp",
"(Inv(F3d)*Expm(-ksi*devlogbetr)*(F3d))*invCp0"\
md.add_filtered_fem_variable("dirmult", mfmult, XP_RG)
md.add_nonlinear_generic_assembly_brick(mim, "(disp-u(1))*dirmult", XP_RG)
-print("Model dofs: %i" % md.nbdof())
-print("Displacement fem dofs: %i" % mfu.nbdof())
+print(f"Model dofs: {md.nbdof()}\nDisplacement fem dofs: {mfu.nbdof()}")
print("Dirichlet mult dofs: %i" % md.mesh_fem_of_variable("dirmult").nbdof())
-if (export_results):
- shutil.copyfile(os.path.abspath(sys.argv[0]),resultspath+"/"+sys.argv[0])
-
+shutil.copyfile(os.path.abspath(sys.argv[0]),resultspath+"/"+sys.argv[0])
starttime_overall = time.process_time()
-
-if (export_results):
- f1 = open("%s/tension_plane_strain.dat" % resultspath, "w")
-
-try:
+with open(f"{resultspath}/tension_plane_strain.dat", "w") as f1:
for step in range(steps_t+1):
md.set_variable("disp", disp*step/float(steps_t))
print('STEP %i: Solving with disp = %g' % (step, md.variable("disp")))
starttime = time.process_time()
- md.solve('noisy', 'max_iter', 25, 'max_res', 1e-10,
- 'lsearch', 'simplest', 'alpha max ratio', 100, 'alpha min', 1,
'alpha mult', 0.6,
- 'alpha threshold res', 1e9)
- print('STEP %i COMPLETED IN %f SEC' % (step,
time.process_time()-starttime))
+ md.solve("noisy", "max_iter", 25, "max_res", 1e-10,
+ "lsearch", "simplest", "alpha max ratio", 100, "alpha min", 1,
"alpha mult", 0.6,
+ "alpha threshold res", 1e9)
+ print("STEP %i COMPLETED IN %f SEC" % (step,
time.process_time()-starttime))
F = gf.asm_generic(mim, 0, "dirmult", XP_RG, md)
print("Displacement %g, total force %g" % (md.variable("disp"), F))
@@ -220,12 +209,11 @@ try:
V = gf.asm_generic(mim, 0, "1", -1, md)
sigma11 = gf.asm_generic(mim, 0, "sigma(1,1)", -1, md)/V
gamma = gf.asm_generic(mim, 0, "gamma", -1, md)/V
- if (export_results):
- f1.write("%.10g %.10g %.10g %.10g %10g %10g\n"
- % (md.variable("disp"), F, A, F/A, sigma11, gamma))
- f1.flush()
+ f1.write("%.10g %.10g %.10g %.10g %10g %10g\n"
+ % (md.variable("disp"), F, A, F/A, sigma11, gamma))
+ f1.flush()
- output = (mfout1, md.local_projection(mim, "sqrt(1.5)*Norm(sigmaD)",
mfout1), "Von Mises Stress",
+ output = (mfout1, md.local_projection(mim, "sqrt(1.5)*Norm(sigmaD)",
mfout1), "Von Mises Stress",
mfout1, md.local_projection(mim, "J", mfout1), "J",
mfout1, md.local_projection(mim, "sigma(1,1)", mfout1),
"Cauchy stress 11",
mfout1, md.local_projection(mim, "sigma(2,2)", mfout1),
"Cauchy stress 22",
@@ -234,17 +222,13 @@ try:
mfu, md.variable("u"), "Displacements",
mfout2, md.interpolation("dirmult", mfout2, XP_RG), "Nominal
reaction traction",
mfout2, md.local_projection(mim, "gamma", mfout2), "plastic
strain")
-
- mfout2.export_to_vtu("%s/tension_plane_strain_%i.vtu" % (resultspath,
step), *output)
+ mfout2.export_to_vtu(f"{resultspath}/tension_plane_strain_{step}.vtu",
*output)
md.set_variable("gamma0", md.interpolation("gamma", 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))
-finally:
- if (export_results) :
- f1.close()
-print('OVERALL SOLUTION TIME IS %f SEC' %
(time.process_time()-starttime_overall))
+print("OVERALL SOLUTION TIME IS %f SEC" %
(time.process_time()-starttime_overall))
diff --git
a/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_3D.py
b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_3D.py
index 431131e8..4c9ec016 100644
---
a/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_3D.py
+++
b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_3D.py
@@ -2,7 +2,7 @@
# -*- coding: UTF8 -*-
# Python GetFEM interface
#
-# Copyright (C) 2020-2020 Konstantinos Poulios.
+# Copyright (C) 2020-2023 Konstantinos Poulios.
#
# This file is a part of GetFEM
#
@@ -67,7 +67,7 @@ resultspath = "./results"
if not os.path.exists(resultspath):
os.makedirs(resultspath)
-tee = subprocess.Popen(["tee", "%s/tension_3D.log" % resultspath],
+tee = subprocess.Popen(["tee", f"{resultspath}/tension_3D.log"],
stdin=subprocess.PIPE)
sys.stdout.flush()
os.dup2(tee.stdin.fileno(), sys.stdout.fileno())
@@ -82,8 +82,7 @@ ZM_RG = 4
# Mesh
mesh2d = gf.Mesh("import", "structured_ball",
- "GT='%s';ORG=[0,0];SIZES=[%f];NSUBDIV=[%i,%i];SYMMETRIES=2"
- % (geotrans2d, H/2., N_R1, N_R2))
+
f"GT='{geotrans2d}';ORG=[0,0];SIZES=[{H/2.}];NSUBDIV=[{N_R1},{N_R2}];SYMMETRIES=2")
mesh = gf.Mesh("prismatic", mesh2d, N_L, disp_fem_order)
trsf = np.zeros([3,3])
@@ -128,7 +127,7 @@ if dH > 0:
pts[2,i] -= (z*dH)/(2*H) * (1 + np.cos(2.*np.pi*x/L))
mesh.set_pts(pts)
-mesh.export_to_vtu("%s/mesh.vtu" % resultspath)
+mesh.export_to_vtu(f"{resultspath}/mesh.vtu")
# FEM
mfN = gf.MeshFem(mesh, N)
@@ -182,17 +181,17 @@ md.add_macro("invCp0", "[[[1,0,0],[0,0,0],[0,0,0]],"+\
md.add_macro("devlogbetr", "Deviator(Logm(F*invCp0*F'))")
md.add_macro("Y0","{A}+{B}*gamma0".format(A=np.sqrt(2./3.)*pl_sigma_0,
B=2./3.*pl_H))
md.add_macro("ksi",
-
"(1-1/max(1,mu*pow(J,-5/3)*Norm(devlogbetr)/Y0))/(2+{B}/(mu*pow(J,-5/3)))"
+ "(1-1/max(1,mu/J*Norm(devlogbetr)/Y0))/(2+{B}/(mu/J))"
.format(B=2./3.*pl_H))
md.add_macro("gamma", "gamma0+ksi*Norm(devlogbetr)")
md.add_macro("devlogbe", "(1-2*ksi)*devlogbetr")
-md.add_macro("tauD", "mu*pow(J,-2/3)*devlogbe")
+md.add_macro("tauD", "mu*devlogbe")
md.add_nonlinear_generic_assembly_brick\
(mim, "((tauH*Id(3)+tauD)*Inv(F')):Grad_Test_u")
-md.add_macro("sigmaD", "(mu*pow(J,-5/3)*devlogbe)")
-md.add_macro("sigma", "tauH/J*Id(3)+mu*pow(J,-5/3)*devlogbe")
+md.add_macro("sigmaD", "tauD/J")
+md.add_macro("sigma", "tauH/J*Id(3)+mu/J*devlogbe")
md.add_macro("invCp", "(Inv(F)*Expm(-ksi*devlogbetr)*(F))*invCp0"\
"*(Inv(F)*Expm(-ksi*devlogbetr)*(F))'")
@@ -200,22 +199,21 @@ md.add_macro("invCp",
"(Inv(F)*Expm(-ksi*devlogbetr)*(F))*invCp0"\
md.add_filtered_fem_variable("dirmult", mfmult, XP_RG)
md.add_nonlinear_generic_assembly_brick(mim, "(disp-u(1))*dirmult", XP_RG)
-print("Model dofs: %i" % md.nbdof())
-print("Displacement fem dofs: %i" % mfu.nbdof())
+print(f"Model dofs: {md.nbdof()}\nDisplacement fem dofs: {mfu.nbdof()}")
print("Dirichlet mult dofs: %i" % md.mesh_fem_of_variable("dirmult").nbdof())
shutil.copyfile(os.path.abspath(sys.argv[0]),resultspath+"/"+sys.argv[0])
starttime_overall = time.process_time()
-with open("%s/tension_3D.dat" % resultspath, "w") as f1:
+with open(f"{resultspath}/tension_3D.dat", "w") as f1:
for step in range(steps_t+1):
md.set_variable("disp", disp*step/float(steps_t))
print('STEP %i: Solving with disp = %g' % (step, md.variable("disp")))
starttime = time.process_time()
- md.solve('noisy', 'max_iter', 25, 'max_res', 1e-10,
- 'lsearch', 'simplest', 'alpha max ratio', 100, 'alpha min',
0.2, 'alpha mult', 0.6,
- 'alpha threshold res', 1e9)
- print('STEP %i COMPLETED IN %f SEC' % (step,
time.process_time()-starttime))
+ md.solve("noisy", "max_iter", 25, "max_res", 1e-10,
+ "lsearch", "simplest", "alpha max ratio", 100, "alpha min",
0.2, "alpha mult", 0.6,
+ "alpha threshold res", 1e9)
+ print("STEP %i COMPLETED IN %f SEC" % (step,
time.process_time()-starttime))
F = gf.asm_generic(mim, 0, "dirmult", XP_RG, md)
print("Displacement %g, total force %g" % (md.variable("disp"), F))
@@ -237,7 +235,7 @@ with open("%s/tension_3D.dat" % resultspath, "w") as f1:
mfu, md.variable("u"), "Displacements",
mfout2, md.interpolation("dirmult", mfout2, XP_RG), "Nominal
reaction traction",
mfout2, md.local_projection(mim, "gamma", mfout2), "plastic
strain")
- mfout2.export_to_vtu("%s/tension_3D_%i.vtu" % (resultspath, step),
*output)
+ mfout2.export_to_vtu(f"{resultspath}/tension_3D_{step}.vtu", *output)
md.set_variable("gamma0", md.interpolation("gamma", mimd1, -1))
md.set_variable("invCp0vec",
@@ -245,5 +243,5 @@ with open("%s/tension_3D.dat" % resultspath, "w") as f1:
" [[0,0,0,0.5,0,0],[0,1,0,0,0,0]
,[0,0,0,0,0.5,0]],"+\
"
[[0,0,0,0,0,0.5],[0,0,0,0,0.5,0],[0,0,1,0,0,0]]]:invCp", mimd6, -1))
-print('OVERALL SOLUTION TIME IS %f SEC' %
(time.process_time()-starttime_overall))
+print("OVERALL SOLUTION TIME IS %f SEC" %
(time.process_time()-starttime_overall))
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] [getfem-commits] branch master updated: Remove erroneous J^-2/3 factor from finite strain plasticity examples,
Konstantinos Poulios <=