Here we discuss about the development of GetFEM and possible new features.
1 Like
Here, as discussed, GetFEM shell elements are not adapted to 3D geometry and are currently only compatible with plates.
https://lists.nongnu.org/archive/html/getfem-users/2018-10/msg00003.html
In order to be able to define shell elements in any orientation, it is necessary to be able to define the local coordinates of each convexity and rotate by direction cosine matrix.
What is the best GetFEM development policy to achieve this?
not a shell model, but we can use this pinched hemisphere model with solid elements and a mixed displacement-pressure formulation as a reference
#!/usr/bin/env python
# -*- coding: UTF8 -*-
#
# Copyright (C) 2023-2023 Konstantinos Poulios.
#
############################################################################
import getfem as gf
import numpy as np
import os, sys, subprocess
import time
import shutil
gf.util_trace_level(1)
# Input data
T = 0.04 # sphere thickness
R = 10. # sphere midsurface radius
N_C = 20 # number of elements per 45 deg of the circumference
N_T = 1 # number of elements through the thickness
E = 6.825e7 # Young's modulus (initial)
nu = 0.3 # Poisson's ratio (initial)
steps = 1
force = 1e0
#------------------------------------
geotrans = "GT_QK(3,2)" # geometric transformation
disp_fem_order = 2 # displacements finite element order
press_fem_order = 1 # pressure finite element order (if positive or zero, then a mixed u-p formulation is used)
integration_degree = 5 # 9 Gauss points per quad
#------------------------------------
mixed_up = (press_fem_order >= 0)
resultspath = "./sphere_solid"
if not os.path.exists(resultspath): os.makedirs(resultspath)
tee = subprocess.Popen(["tee", f"{resultspath}/sphere_solid.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())
# auxiliary constants
XM_RG = 3
YM_RG = 4
ZM_RG = 5
INT_RG = 7
EXT_RG = 8
mtmp = gf.Mesh("import", "structured",
f"GT='{geotrans}';ORG=[0,0,0];SIZES=[1,1,1];NSUBDIV=[{N_C},{N_C},{N_T}]")
pts = mtmp.pts()
for i in range(pts.shape[1]):
[x,y,z] = pts[:,i]
x = np.tan(x*np.pi/4)
y = np.tan(y*np.pi/4)
aux = (R - T/2 + T*z)/np.sqrt(x*x+y*y+1)
pts[:,i] = aux * np.array([1,x,y])
mtmp.set_pts(pts)
mesh = gf.Mesh("clone", mtmp)
mtmp.transform([[0,1,0],[0,0,1],[1,0,0]])
mesh.merge(mtmp)
mtmp.transform([[0,1,0],[0,0,1],[1,0,0]])
mesh.merge(mtmp)
mesh.optimize_structure()
N = mesh.dim()
mesh.set_region(XM_RG, mesh.outer_faces_in_box([-1e-5,-1e-5,-1e-5], [1e-5,R+T/2+1e-5,R+T/2+1e-5]))
mesh.set_region(YM_RG, mesh.outer_faces_in_box([-1e-5,-1e-5,-1e-5], [R+T/2+1e-5,1e-5,R+T/2+1e-5]))
mesh.set_region(ZM_RG, mesh.outer_faces_in_box([-1e-5,-1e-5,-1e-5], [R+T/2+1e-5,R+T/2+1e-5,1e-5]))
mesh.set_region(EXT_RG, mesh.outer_faces_with_direction([1,1,1], np.pi/2))
mesh.set_region(INT_RG, mesh.outer_faces())
mesh.region_subtract(INT_RG, XM_RG)
mesh.region_subtract(INT_RG, YM_RG)
mesh.region_subtract(INT_RG, ZM_RG)
mesh.region_subtract(INT_RG, EXT_RG)
mesh.export_to_vtu(f"{resultspath}/mesh.vtu")
# FEM
mfN = gf.MeshFem(mesh, N)
mfN.set_classical_fem(disp_fem_order)
xloaded_dof = list(set.intersection(set(mfN.basic_dof_on_region(YM_RG)[0::N]),
set(mfN.basic_dof_on_region(ZM_RG)[0::N]),
set(mfN.basic_dof_on_region(EXT_RG)[0::N])))
yloaded_dof = list(set.intersection(set(mfN.basic_dof_on_region(XM_RG)[1::N]),
set(mfN.basic_dof_on_region(ZM_RG)[1::N]),
set(mfN.basic_dof_on_region(EXT_RG)[1::N])))
zfixed_dof = list(set.intersection(set(mfN.basic_dof_on_region(XM_RG)[2::N]),
set(mfN.basic_dof_on_region(YM_RG)[2::N]),
set(mfN.basic_dof_on_region(EXT_RG)[2::N])))
kept_dofs = list(set(range(mfN.nbdof()))-set(mfN.basic_dof_on_region(XM_RG)[0::N]) # remove all x-dofs on XM_RG
-set(mfN.basic_dof_on_region(YM_RG)[1::N]) # remove all y-dofs on YM_RG
# -set(mfN.basic_dof_on_region(ZM_RG)[2::N]) # remove all z-dofs on ZM_RG
-set(zfixed_dof)) # remove z-dof on the pole
mfu = gf.MeshFem("partial", mfN, kept_dofs)
if mixed_up:
mfp = gf.MeshFem(mesh, 1)
mfp.set_classical_fem(press_fem_order)
mfout = gf.MeshFem(mesh)
mfout.set_classical_discontinuous_fem(disp_fem_order)
mim = gf.MeshIm(mesh, integration_degree)
# Model
md = gf.Model("real")
md.add_initialized_data("K", E/(3*(1-2*nu)))
md.add_initialized_data("G", E/(2*(1+nu)))
md.add_fem_variable("u", mfu) # displacements
md.add_macro("F", "(Id(3)+Grad_u)")
md.add_macro("J", "Det(F)")
if mixed_up:
# md.add_fem_variable("p_o_J", mfp)
# md.add_macro("p", "-p_o_J*J")
# md.add_nonlinear_term(mim, "(p_o_J/K+log(J)/J)*Test_p_o_J")
md.add_fem_variable("p", mfp)
md.add_nonlinear_term(mim, "(p-K*log(J))*Test_p")
else:
md.add_macro("p", "K*log(J)")
md.add_macro("devtau", "G*pow(J,-2/3)*Deviator(F*F')")
md.add_macro("tau", "p*Id(3)+devtau")
md.add_nonlinear_term(mim, "(tau*Inv(F)'):Grad(Test_u)")
ib_load = md.add_explicit_rhs("u", np.zeros(mfu.nbdof()))
LOAD = np.zeros(mfu.nb_basic_dof())
LOAD[xloaded_dof] = force
LOAD[yloaded_dof] = -force
LOAD = mfu.reduce_vector(LOAD)
print(f"Displacement dofs: {mfu.nbdof()}")
print(f"Total model dofs: {md.nbdof()}")
mass_mat = gf.asm_mass_matrix(mim, mfout)
shutil.copyfile(os.path.abspath(sys.argv[0]),resultspath+"/"+sys.argv[0])
starttime_overall = time.process_time()
with open(f"{resultspath}/sphere_solid.dat", "w") as f:
for step in range(steps+1):
print("SOLVING STEP %i" % step)
starttime = time.process_time()
md.set_private_rhs(ib_load, LOAD*(step/steps))
iters = 40
nit = \
md.solve("noisy", "lsolver", "mumps", "max_iter", iters, "max_res", 1e-7, #)[0]
"lsearch", "simplest", "alpha max ratio", 1., "alpha min", 0.5, "alpha mult", 0.6)[0]
print("STEP %i COMPLETED IN %f SEC" % (step, time.process_time()-starttime))
args = ("t", True, mfout, np.zeros(mfout.nbdof()), "select_output", "t")
VM = gf.asm_generic(mim, 1, "sqrt(1.5)*Norm(devtau)/J*Test_t", -1, md, *args)
VM = np.transpose(gf.linsolve_mumps(mass_mat, VM))
output = (mfout, VM, "Von Mises Stress")
for i,j in [[1,1],[2,2],[3,3],[1,2],[2,1]]:
TAU = gf.asm_generic(mim, 1, f"tau({i},{j})*Test_t", -1, md, *args)
SIGMA = gf.asm_generic(mim, 1, f"tau({i},{j})/J*Test_t", -1, md, *args)
TAU = np.transpose(gf.linsolve_mumps(mass_mat, TAU))
SIGMA = np.transpose(gf.linsolve_mumps(mass_mat, SIGMA))
output += (mfout, TAU, f"Kirchhoff stress matrix {i}{j}")
output += (mfout, SIGMA, f"Cauchy Stress matrix {i}{j}")
for i,j in [[1,1],[2,2],[3,3],[1,2],[2,1],[2,3],[3,2],[3,1],[1,3],]:
B = gf.asm_generic(mim, 1, f"(F*F')({i},{j})*Test_t", -1, md, *args)
B = np.transpose(gf.linsolve_mumps(mass_mat, B))
output += (mfout, B, f"Left Cauchy Green matrix {i}{j}")
if mixed_up:
# output = (mfp, md.variable("p_o_J"), "Pressure over J") + output
output = (mfp, md.variable("p"), "Pressure") + output
output = (mfu, md.variable("u"), "Displacements") + output
mfout.export_to_vtu(f"{resultspath}/sphere_solid_elements_{step}.vtu", *output)
U = mfu.extend_vector(md.variable("u"))
print(U[xloaded_dof])
f.write(f"step={step} force={force*step/steps} disp={U[xloaded_dof]}\n")
f.flush()
print("OVERALL SOLUTION TIME IS %f SEC" % (time.process_time()-starttime_overall))
1 Like
If I understand correctly, does this script indicate that the mass matrix can be used as a direction cosine matrix? (I will do an assembly later to check my understanding.)
https://getfem.org/userdoc/model_mass.html
no the “mass_mat” matrix appearing in the script is used only for post-processing, in order to project results on the mfout mesh_fem (projection of fields makes more sense than interpolation). But it is a minor detail, it is not part of the model.
1 Like
Sorry for my lack of understanding. I’ll read the script a little more carefully.