Large sliding contact brick raytracing with surface mesh as Master surface

Questions and issues
1

Hello everyone!

I’m quite new to getfem and the field of FEM simulation in general, so sorry in advance if my question is obvious.

I have a problem, where I want to inflate an object inside of another object, say a ball inside a hollow box. The ball shall not be able to penetrate the box, I want to implement a contact problem. I’ve put together this little test script:

import getfem as gf
import numpy as np

VOL = True

mo_sphere = gf.MesherObject('ball', [0.0, 0.0, 1.0], 1.0)
mo_bigbox = gf.MesherObject('rectangle', [-2.0, -2.0, -0.5], [2.0, 2.0, 2.5])
mo_smallbox = gf.MesherObject('rectangle', [-1.0, -1.0, 0.0], [1.0, 1.0, 2.0])
mo_box = gf.MesherObject('set minus', mo_bigbox, mo_smallbox)

m_slave = gf.Mesh('generate', mo_sphere, 0.5)
if VOL:
    m_master = gf.Mesh('generate', mo_box, 0.5)
else:
    m_master = gf.Mesh('empty', 3)
    pts = np.array([
        [-2.0, -2.0, -0.5], # 0: Links-Vorne-Unten
        [ 2.0, -2.0, -0.5], # 1: Rechts-Vorne-Unten
        [ 2.0,  2.0, -0.5], # 2: Rechts-Hinten-Unten
        [-2.0,  2.0, -0.5], # 3: Links-Hinten-Unten
        [-2.0, -2.0,  2.5], # 4: Links-Vorne-Oben
        [ 2.0, -2.0,  2.5], # 5: Rechts-Vorne-Oben
        [ 2.0,  2.0,  2.5], # 6: Rechts-Hinten-Oben
        [-2.0,  2.0,  2.5]  # 7: Links-Hinten-Oben
    ]).T
    triangles = np.array([
        [0, 2, 1], [0, 3, 2], # Boden
        [4, 5, 6], [4, 6, 7], # Deckel
        [0, 1, 5], [0, 5, 4], # Vorne
        [1, 2, 6], [1, 6, 5], # Rechts
        [2, 3, 7], [2, 7, 6], # Hinten
        [3, 0, 4], [3, 4, 7] # Links
    ], dtype=np.int32)

    m_master.add_point(pts)
    for t in triangles:
        m_master.add_convex(gf.GeoTrans("GT_PK(2,1)"), pts[:, t])


if VOL:
    m_master.set_region(MASTER_CONTACT := 1, m_master.outer_faces_in_box([-1.1, -1.1, -0.1], [1.1, 1.1, 2.1]))
else:
    m_master.set_region(MASTER_CONTACT := 1, np.vstack([m_master.cvid(), np.full_like(m_master.cvid(), -1)]))
m_slave.set_region(SLAVE_CONTACT := 2, m_slave.outer_faces())
m_slave.set_region(SLAVE_DIRICHLET := 3, m_slave.outer_faces_with_direction([-0.1, -0.1, -0.1], 100))

print(f"MASTER_REGION:   {m_master.region(MASTER_CONTACT).shape}")
print(f"SLAVE_REGION:    {m_slave.region(SLAVE_CONTACT).shape}")
print(f"SLAVE_DIRICHLET: {m_slave.region(SLAVE_DIRICHLET).shape}")

mim_slave = gf.MeshIm(m_slave, 4)
mim_master = gf.MeshIm(m_master, 4)
mf_slave = gf.MeshFem(m_slave, 3)
mf_master = gf.MeshFem(m_master, 3)
mf_slave.set_classical_fem(1)
mf_master.set_classical_fem(1)

mf_debug = gf.MeshFem(m_master, 1) 
mf_debug.set_classical_fem(1) 
region_mask = np.zeros(mf_debug.nb_basic_dof())
dofs_on_box = mf_debug.basic_dof_on_region(MASTER_CONTACT)
region_mask[dofs_on_box] = 1.0
mf_master.export_to_vtk("master_mesh.vtk", mf_debug, region_mask, "RegionMarker")

mf_debug = gf.MeshFem(m_slave, 1)
mf_debug.set_classical_fem(1)
region_mask = np.zeros(mf_debug.nb_basic_dof())
dofs_on_box = mf_debug.basic_dof_on_region(SLAVE_CONTACT)
region_mask[dofs_on_box] = 1.0
mf_debug2 = gf.MeshFem(m_slave, 1)
mf_debug2.set_classical_fem(1)
region_mask2 = np.zeros(mf_debug2.nb_basic_dof())
dofs_on_box2 = mf_debug2.basic_dof_on_region(SLAVE_DIRICHLET)
region_mask2[dofs_on_box2] = 1.0
mf_slave.export_to_vtk("slave_mesh.vtk", mf_debug, region_mask, "RegionMarker", mf_debug2, region_mask2, "DirichletMarker")

md = gf.Model("real")
md.add_initialized_data("pressure", 3.0)
md.add_initialized_data("lambda", [1])
md.add_initialized_data("mu", [1])
md.add_initialized_data("r", 1)
md.add_initialized_data("alpha", 1)
md.add_initialized_data("f", 0)

md.add_fem_variable("ub", mf_slave)
md.add_filtered_fem_variable("uq", mf_master, MASTER_CONTACT)
md.add_filtered_fem_variable("egkdsfhgldfskjghlambda", mf_slave, SLAVE_CONTACT)

md.add_isotropic_linearized_elasticity_brick(mim_slave, "ub", "lambda", "mu")
md.add_source_term_brick(mim_slave, "ub", f"pressure * Normal", SLAVE_CONTACT)

md.add_Dirichlet_condition_with_multipliers(mim_slave, "ub", 1, SLAVE_DIRICHLET)
md.add_Dirichlet_condition_with_multipliers(mim_master, "uq", 1, MASTER_CONTACT)

ind = md.add_integral_large_sliding_contact_brick_raytracing("r", 5, "f", "alpha", 0)
md.add_master_contact_boundary_to_large_sliding_contact_brick(ind, mim_master, MASTER_CONTACT, "uq")
md.add_slave_contact_boundary_to_large_sliding_contact_brick(ind, mim_slave, SLAVE_CONTACT, "ub", "egkdsfhgldfskjghlambda")


md.solve("noisy", "max_iter", 500, "max_res", 1E-8)
ub = md.variable("ub")
mf_slave.export_to_vtk("displacement.vtk", ub, "ub")

This works, as long as VOL is set to True. With VOL set to true, the box containing the ball is volumetric, and I use only the inner box surface for the contact problem, found by outer_faces_in_box(). However, in the original problem, my outer boundary is a surface mesh. It is inconvenient to volumetrize the surface, just to able to find an inner surface with outer_faces_in_box(). I figured, it thought be possible to do the contact problem surface mesh directly. However this fails:
Trace 2 in getfem_models.cc, line 4399: Mass term assembly for Dirichlet condition
Trace 2 in getfem_models.cc, line 4399: Mass term assembly for Dirichlet condition
Trace 2 in getfem_models.cc, line 6051: Linearized isotropic elasticity: generic matrix assembly
Trace 2 in getfem_models.cc, line 3308: Generic source term assembly
Trace 2 in getfem_models.cc, line 3319: Source term: generic source term assembly
Trace 2 in getfem_models.cc, line 4399: Mass term assembly for Dirichlet condition
Trace 2 in getfem_models.cc, line 4399: Mass term assembly for Dirichlet condition
Trace 2 in getfem_models.cc, line 2652: Global generic assembly RHS
[1] 16857 segmentation fault (core dumped) python ./problem_demo.py
My question is now: Am I doing something wrong? If yes, what? Or is the contact brick not build for this application, is it not even possible?

Thanks in advance, I appreciate every help!
P.S.:

I’m aware of the fact that the Lame constants for my material don’t make any sense, this was simply not my focus yet.

P.P.S:

The mesh generation for the volumetric mesh sometimes fails and sometimes works, with different error messages, either:


QH6019 qhull input error (qh_scalelast): can not scale last coordinate to [   0,    0].  Input is cocircular or cospherical.   Use option ‘Qz’ to add a point at infinity.

While executing:  | qhull QJ d Qbb Pp T0
Options selected for Qhull 2020.2.r 2020/08/31:
run-id 2054137973  delaunay  Qbbound-last  Pprecision-ignore  Pgood  _run 1
QJoggle 6e-11  _joggle-seed 16807  _maxoutside  0
logic_error exception caught
Traceback (most recent call last):
File “/home/..././problem_demo.py”, line 13, in 
m_master = gf.Mesh(‘generate’, mo_box, 0.5)
File “/home/.../getfem-5.4.4/.venv/lib/python3.13/site-packages/getfem/getfem.py”, line 1164, in init
generic_constructor(self,‘mesh’,*args)
~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^
File “/home/.../getfem-5.4.4/.venv/lib/python3.13/site-packages/getfem/getfem.py”, line 65, in generic_constructor
self.id = getfem_from_constructor(clname,*args)
~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^
RuntimeError: (Getfem::InterfaceError) – Error in gmm/gmm_dense_qr.h, line 786 void gmm::symmetric_qr_algorithm(const MAT1&, const VECT&, const MAT2&, typename number_traits<typename linalg_traits::value_type>::magnitude_type, bool) [with MAT1 = dense_matrix; VECT = std::vector; MAT2 = dense_matrix; typename number_traits<typename linalg_traits::value_type>::magnitude_type = double; typename linalg_traits::value_type = double]:
QR algorithm failed.

or:
python: bgeot_kdtree.cc:115: std::unique_ptr<bgeot::kdtree_elt_base> bgeot::build_tree_(ITER, ITER, unsigned int): Assertion `(itmedian).n[dir] >= median && median >= ((itmedian-1)).n[dir]’ failed.
[1] 18907 IOT instruction (core dumped) python ./problem_demo.py

I don’t know if I have done something wrong or whether this is an issue. But as it succeeds sometimes, I can test, and I don’t have to use the mesh generator later, so it’s not a concern for me now.

2

hello and welcome.

Do you want to just model contact between a ball and a rigid box?

If the box is supposed to be rigid, you can just use 6 times the add_rigid_obstacle_to_large_sliding_contact_brick function. Once per each face of the box.

Some off-topic hints,

  1. avoid linear tetrahedrals at all cost
  2. replace vtk with vtu
  3. avoid the pre-implemented “bricks”, nowadays if you know your math an the GWFL syntax you only need the “add_linear_term” or “add_nonlinear_term” functions
  4. try to addgf.util_trace_level(1) and gf.util_warning_level(1) in the beginning of your script
  5. for general mesh generation prefer gmsh over getfem. Only the structures meshes are quite nice in GetFEM, we have structured meshes e.g. for cuboid, disk, sphere, and cylinders are easy to create as well.
  6. get inspiration from my snippet collection available at kopo · GitLab
3

Hello, thanks for the kind welcome and the quick response!

The problem is indeed, that I want it for arbitrary shapes. Some context:
I want to do FEM simulation on organ models extracted from CT images, so no, it shall not be a box in the end. The ball box is simply a dummy I’ve put together to illustrate my problem. As far as I understand, add_rigid_obstacle_to_large_sliding_contact_brick does only work for analytically defined objects, where I can have an explicit distance function. So it won’t help me with my problem, where I have an arbitrary hollow mesh, used as a cage for an expanding inner model, right? (The cage, in this case box, shall indeed be rigid)
A side note on gmsh: My current workflow (if it deserves that name) is exporting my model from 3dSlicer, the medical imaging software in use, and volumetrizing it with gmsh, means putting a box around it and meshing the space between the box and the model. However, 3dSlicer and gmsh don’t work together very well, gmsh is not usable from within 3dSlicer, Slicer crashes, most likely because both softwares want to take control of the same resources. I created a workaround by exporting the gmsh commands into a temporary script and running it in a subprocess, but this probably doesn’t deserve to be called a solution. Long story short: I want to avoid gmsh, and I want to avoid the volumetrization in total. As I only use the inner surface of the box as the master surface for the contact, the rest of the model is useless, I thought it must be possible to avoid this heckle.

On your hints:
Thanks for all of them! I appreciate them and I will surely follow them!
On 3., would you say I should also avoid the contact brick and do it via add_raytracing_transformationinstead?

4

ok, it helps to understand what your objective is.

Regarding use of add_integral_large_sliding_contact_brick_raytracing or add_raytracing_transformation+GWFL, I would say in this case it is in general justifiable to use the brick, since it actually simplifies things significantly. However, when you have a non-standard use case like yours, there are things you can only achieve using the add_raytracing_transformation+GWFL solution.

Now coming to your case, in GetFEM you cannot avoid the volumetrization of the boundary. This is because the “Normal” in GWFL is defined always on the edges of a surface and tangent to the surface if your domain is a surface even in 3D. To get “Normal” to be normal on the surface, this surface needs to be a boundary of a volume domain. Sorry for this inconvenience, but it is not a big deal. I suggest you start with your surface mesh and you give it some thickness outwards. So if you start with triangle mesh, you will end up with a mesh of prisms/wedges (extruded triangles).

If the master side is always rigid, it will not matter that you basically double the nodes of the surface mesh. Then the next step is to define the gap function using add_raytracing_transformation. Conceptually your code should look something like this for a scalar Lagrange multiplier field (with penalty it would be even simpler):

md.add_fem_variable("u", mf_slave)
md.add_fem_data("u_rigid", mf_master) # u_rigid is just data (zero), does not require any equation
md.add_filtered_fem_variable("contact_mult", mfmult, SLAVE_RG) # mfmult is scalar mesh_fem

...
md.add_raytracing_transformation("raytrace", release_dist)
md.add_master_contact_boundary_to_raytracing_transformation\
  ("raytrace", m_master, "u_rigid", MASTER_RG)
md.add_slave_contact_boundary_to_raytracing_transformation\
  ("raytrace", m_slave, "u", SLAVE_RG)

md.add_macro("normal", "Normalized(Inv(F(u)').Normal)")
md.add_macro("gap", "(Interpolate(X,raytrace)-X-u).normal")
md.add_nonlinear_term(mim_contact, # everywhere
  "-1/c*contact_mult*Test_contact_mult", SLAVE_RG) # c is the augmentation parameter
md.add_nonlinear_term(mim_contact, # wherever contact is detected
  "Interpolate_filter(raytrace,"
                     "-1/c*neg_part(contact_mult+c*gap)*Test_contact_mult, 1)",
  SLAVE_RG)

From our 2014 paper:

image
image874×217 28.3 KB