# General development and new features ideas

**URL:** <https://getfem.discourse.group/t/general-development-and-new-features-ideas/20>\
**Category:** Development\
**Tags:** development\
**Created:** [March 28, 2023, 11:00am UTC](https://getfem.discourse.group/t/general-development-and-new-features-ideas/20 "2023-03-28T11:00:02Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [March 28, 2023, 11:00am UTC](https://getfem.discourse.group/t/general-development-and-new-features-ideas/20/1 "2023-03-28T11:00:02Z")

</div>

Here we discuss about the development of GetFEM and possible new features.

---

<div class="post-metadata">

**Author:** ![tkoyama010](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/tkoyama010/32/29_2.png) [@tkoyama010](https://getfem.discourse.group/u/tkoyama010)\
**Post date:** [March 27, 2023, 2:57am UTC](https://getfem.discourse.group/t/general-development-and-new-features-ideas/20/2 "2023-03-27T02:57:52Z")

</div>

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](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?

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [March 29, 2023, 7:37am UTC](https://getfem.discourse.group/t/general-development-and-new-features-ideas/20/3 "2023-03-29T07:37:55Z")

</div>

not a shell model, but we can use this pinched hemisphere model with solid elements and a mixed displacement-pressure formulation as a reference

```python
#!/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))

```

---

<div class="post-metadata">

**Author:** ![tkoyama010](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/tkoyama010/32/29_2.png) [@tkoyama010](https://getfem.discourse.group/u/tkoyama010)\
**Post date:** [March 29, 2023, 11:04am UTC](https://getfem.discourse.group/t/general-development-and-new-features-ideas/20/4 "2023-03-29T11:04:07Z")

</div>

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](https://getfem.org/userdoc/model_mass.html)

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [March 29, 2023, 11:19am UTC](https://getfem.discourse.group/t/general-development-and-new-features-ideas/20/5 "2023-03-29T11:19:31Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![tkoyama010](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/tkoyama010/32/29_2.png) [@tkoyama010](https://getfem.discourse.group/u/tkoyama010)\
**Post date:** [March 29, 2023, 11:47am UTC](https://getfem.discourse.group/t/general-development-and-new-features-ideas/20/6 "2023-03-29T11:47:23Z")

</div>

Sorry for my lack of understanding. I’ll read the script a little more carefully.