Resources for Learning GetFEM++ with Python and GWFL

General
1

Hello everyone,

First of all, I apologize if this topic has already been discussed here.

I am new to GetFEM++ and successfully installed it using a Conda version available on GitHub. After that, I started exploring the official documentation, but I noticed that it contains only a few examples.

I would like to know if there are any more comprehensive resources for learning GetFEM++ with Python, such as tutorials or projects that I could directly use.

Additionally, I am particularly interested in the GWFL (Generic Weak Form Language), but I have found very little documentation on this feature, apart from a C++ version. Are there any specific documents or examples available for using GWFL with Python?

Thank you in advance for your help!

Best regards,

Djoumessi

1 Like
2

If you have not seen it, the place to start is our paper on GWFL
https://dl.acm.org/doi/10.1145/3412849

and this tutorial

For continuum mechanics I recommend the relatively up-to-date examples found in the folder

You can find a lot of useful tips in the snippet collection here
https://gitlab.gbar.dtu.dk/explore/snippets

3

Thank you verry much for your valuable ressources.
Best.

4

Dear Konstantinos,

I am currently reviewing in detail the information you sent me. One of my goals is to reproduce this result as an image (contact mechanics). Could you please provide me with some advice on this, if possible?

Thank you in advance for your help.

internalselfcontact
internalselfcontact1200Ă—1266 56.8 KB

5

if you want I can give you the code for that example, or I can give you the code for this similar example:

image
image1117Ă—732 77.4 KB

#!/usr/bin/env python
# -*- coding: UTF8 -*-
#
############################################################################
import numpy as np
import getfem as gf

gf.util_trace_level(1)
gf.util_warning_level(1)
PI = np.pi

# Input data
R = 50.     # ring radius
r = 5.      # cross section radius

NC = 26     # Number of elements along the ring circumference
Nr1 = 2     # Number of elements in the cross section core
Nr2 = 2     # Number of elements in the cross section peel

E = 1e5      # Young's modulus
nu = 0.45    # Poisson's ratio

mult_fem_order = 1  # contact multipliers finite element order
integration_degree = 5         # 27 gauss points per hexa
integration_contact_degree = 7 # 16 gauss points per quad

c_aug = 1e-2*E # Augmentation parameter
release_dist = r

# auxiliary constants
EXTERIOR_RG = 1
HANDLE1_RG = 2
HANDLE2_RG = 3


# Mesh generation:
# - create a 2D mesh for a circular disc
# - extrude the 2D mesh into a cylinder
# - wrap the cylinder into a closed ring
# - clone the ring mesh to get a second ring

m_ = gf.Mesh("import", "structured_ball",
              "GT='GT_QK(2,2)';ORG=[0,0];SIZES=[%g];NSUBDIV=[%i,%i]"
               % (r, Nr1, Nr2))
m0 = gf.Mesh("prismatic", m_, NC, 2) # second order

pts = m0.pts()
for i in range(pts.shape[1]):
  x = pts[0,i]
  z = pts[2,i]
  pts[0,i] = (R+x)*np.cos(2*PI*z)
  pts[2,i] = (R+x)*np.sin(2*PI*z)
m0.set_pts(pts)

m1 = gf.Mesh("empty", 3) 
m1.merge(m0, tol=1e-3) # this trick is used just for gluing the open ends of the ring in m0

m1.set_region(EXTERIOR_RG, m1.outer_faces())
m1.set_region(HANDLE1_RG, m1.outer_faces_in_box([-R-1.1*r,-1.1*r,-0.3*R],
                                                [-R+1.1*r, 1.1*r, 0.3*R]))
m1.set_region(HANDLE2_RG, m1.outer_faces_in_box([R-1.1*r,-1.1*r,-0.3*R],
                                                [R+1.1*r, 1.1*r, 0.3*R]))
m1.optimize_structure()
m1.export_to_vtu("mesh1.vtu")

m2 = gf.Mesh("clone", m1)
m2.translate([-R,2.5*r,0])
m2.export_to_vtu("mesh2.vtu")


# Finite element spaces and integration methods

mf1 = gf.MeshFem(m1, 3)
mf1.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(3)"))

mf2 = gf.MeshFem(m2, 3)
mf2.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(3)"))

mfmult1 = gf.MeshFem(m1, 3)
mfmult1.set_classical_fem(mult_fem_order)

mfmult2 = gf.MeshFem(m2, 3)
mfmult2.set_classical_fem(mult_fem_order)

mim1 = gf.MeshIm(m1, integration_degree)
mim2 = gf.MeshIm(m2, integration_degree)
mim_contact1 = gf.MeshIm(m1, integration_contact_degree)
mim_contact2 = gf.MeshIm(m2, integration_contact_degree)


# Model object

md = gf.Model("real")

md.add_fem_variable("u1", mf1)
md.add_fem_variable("u2", mf2)
md.add_filtered_fem_variable("contact_mult1", mfmult1, EXTERIOR_RG)
md.add_filtered_fem_variable("contact_mult2", mfmult2, EXTERIOR_RG)
#mf_contact_mult = md.mesh_fem_of_variable("contact_mult")

# Elasticity
md.add_initialized_data("K", E/(3*(1-2*nu)))
md.add_initialized_data("G", E/(2*(1+nu)))
md.add_macro("F(u)", "Id(3)+Grad_u")
md.add_macro("tau(u)", "K*log(Det(F(u)))*Id(3)+G*pow(Det(F(u)),-2/3)*Deviator(F(u)*F(u)')")
md.add_nonlinear_term(mim1, "(tau(u1)*Inv(F(u1)')):Grad_Test_u1")
md.add_nonlinear_term(mim2, "(tau(u2)*Inv(F(u2)')):Grad_Test_u2")

# Lagrange multipliers for handles
md.add_variable("q11",3) # ring 1 handle 1
md.add_variable("q12",3) # ring 1 handle 2
md.add_variable("q21",3) # ring 2 handle 1
md.add_variable("q22",3) # ring 2 handle 2
md.add_initialized_data("p11", [0,0,0]) # we will keep this fixed
md.add_initialized_data("p12", [0,0,0]) # we will keep this fixed only at the beginning
md.add_initialized_data("p21", [0,0,0]) # we will drive this only at the beginning
md.add_initialized_data("p22", [0,0,0]) # we will drive this till the end

md.add_linear_term(mim1, "(u1-p11).q11", HANDLE1_RG)
md.add_linear_term(mim1, "(u1-p12).q12", HANDLE2_RG)
md.add_linear_term(mim2, "(u2-p21).q21", HANDLE1_RG)
md.add_linear_term(mim2, "(u2-p22).q22", HANDLE2_RG)

md.add_initialized_data("penal", 0.)
md.add_linear_term(mim1, "1e4*(u1-p11).Test_u1", HANDLE1_RG) # add some elastic support to handle 1 of ring 1 to avoid twisting
md.add_linear_term(mim1, "penal*(u1-p12).Test_u1", HANDLE2_RG) # add some elastic support to handle 2 of ring 1 to avoid twisting
md.add_linear_term(mim2, "penal*(u2-p21).Test_u2", HANDLE1_RG) # add some elastic support to handle 1 of ring 2 to avoid twisting
md.add_linear_term(mim2, "1e4*(u2-p22).Test_u2", HANDLE2_RG) # add some elastic support to handle 2 of ring 2 to avoid twisting


# Add contact terms to the model

md.add_initialized_data("c", c_aug)
if False: # use large sliding contact brick
  md.add_initialized_data("alpha", 1.)
  md.add_initialized_data("f", 0.)
  ibc = md.add_integral_large_sliding_contact_brick_raytracing("c", release_dist, "f", "alpha", 1) # 1 stands for the more symmetric version
  md.add_master_slave_contact_boundary_to_large_sliding_contact_brick(ibc, mim_contact1, EXTERIOR_RG, "u1", "contact_mult1", "")
  md.add_master_slave_contact_boundary_to_large_sliding_contact_brick(ibc, mim_contact2, EXTERIOR_RG, "u2", "contact_mult2", "")
else: # implement the contact equations directly
  md.add_raytracing_transformation("raytrace", release_dist)
  md.define_variable_group("u", "u1", "u2")
  md.add_master_contact_boundary_to_raytracing_transformation("raytrace", m1, "u", EXTERIOR_RG)
  md.add_master_contact_boundary_to_raytracing_transformation("raytrace", m2, "u", EXTERIOR_RG)

  md.add_nonlinear_term(mim_contact1, "-contact_mult1.Test_u1", EXTERIOR_RG)
  md.add_nonlinear_term(mim_contact2, "-contact_mult2.Test_u2", EXTERIOR_RG)
  md.add_nonlinear_term(mim_contact1, "-(1/c)*contact_mult1.Test_contact_mult1", EXTERIOR_RG)
  md.add_nonlinear_term(mim_contact2, "-(1/c)*contact_mult2.Test_contact_mult2", EXTERIOR_RG)
  md.add_macro("normal1", "Normalized(Inv(F(u1)').Normal)")
  md.add_macro("normal2", "Normalized(Inv(F(u2)').Normal)")
  md.add_macro("gap1", "(Interpolate(X,raytrace)+Interpolate(u,raytrace)-X-u1).normal1")
  md.add_macro("gap2", "(Interpolate(X,raytrace)+Interpolate(u,raytrace)-X-u2).normal2")
  md.add_nonlinear_term(mim_contact1,
    "Interpolate_filter(raytrace,(-1/c)*neg_part(contact_mult1.normal1+c*gap1)*normal1.Test_contact_mult1, 1)", EXTERIOR_RG)
  md.add_nonlinear_term(mim_contact2,
    "Interpolate_filter(raytrace,(-1/c)*neg_part(contact_mult2.normal2+c*gap2)*normal2.Test_contact_mult2, 1)", EXTERIOR_RG)
  md.add_nonlinear_term(mim_contact1,
                        "Interpolate_filter(raytrace, -neg_part(contact_mult1.normal1+c*gap1)*normal1.Interpolate(Test_u,raytrace), 1)",
                        EXTERIOR_RG)
  md.add_nonlinear_term(mim_contact2,
                        "Interpolate_filter(raytrace, -neg_part(contact_mult2.normal2+c*gap2)*normal2.Interpolate(Test_u,raytrace), 1)",
                        EXTERIOR_RG) 


# Perform the simulation

print("Model dofs: %i" % md.nbdof())
print("Displacement fem dofs: %i" % (mf1.nbdof()+mf2.nbdof()))
#print("Contact mult dofs: %i" % mf_contact_mult.nbdof())

for step in range(90):
  if step <= 30: # initial phase
    t = step/30
    md.set_variable("p21", [1.8*R*(max(t,min(4*t,0.2))),
                            -5*r*min(t/0.4,1), #-5*r*(np.cos(PI*t)-1)/2,
                            0])
    md.set_variable("p22", [1.8*R*max(0,(t-0.4)/0.6)**2,
                            -5*r*min(t/0.4,1),
                            0])
  elif step <= 38: # second phase, releasing handles 21 and 22
    if step == 31:
      md.disable_variable("q12")
      md.disable_variable("q21")
      md.set_variable("q12", [0,0,0])
      md.set_variable("q21", [0,0,0])
    md.set_variable("penal", (1e2,3e1,1e1,3e0,1e0,1e-1,1e-3,0.)[step-31])
  else: # final phase
    t = 1+(step-38)/100
    if t > 1.25: # slow down
      t = 1.25 + 0.2*(t-1.25)
    md.set_variable("p22", [1.8*R*max(0,(t-0.4)/0.6)**2,
                            -5*r,
                            0])
  print("Step", step, "t =", t)
  md.solve("noisy", "max_iter", 20, "max_res", 1e-8,
           "lsearch", "simplest", "alpha max ratio", 3., "alpha min", 0.02, "alpha mult", 0.6)
  mf1.export_to_vtu(f"ring1_{step}.vtu", mf1, md.variable("u1"), "Displacements")
  mf2.export_to_vtu(f"ring2_{step}.vtu", mf2, md.variable("u2"), "Displacements")
6

Dear Konstantinos

Thank you so much for the code you shared—it’s truly impressive! I really appreciate the level of detail and expertise in it. I would take my time to understand, it is really an avance case. Looking at the code, I see that if I’m not wrong, the implemented method is the Lagrange multiplier method with friction.

I was wondering if this code could be adapted to my self-contact case? And if you happen to have a specific code for that example, I would be very interested, as it would help me understand and progress much faster.

Thanks again for your valuable help!

7

Dear Konstantinos,
I took the time to read and understand the code you provided. My question about self-contact is now resolved. This code snippet effectively takes self-contact into account. Please correct me if I’m wrong. I will reach out again if I need more information.

Thank you in advance!

md.add_master_contact_boundary_to_raytracing_transformation("raytrace", m1, "u", EXTERIOR_RG)
md.add_master_contact_boundary_to_raytracing_transformation("raytrace", m2, "u", EXTERIOR_RG)
8

yes, it does. The unusual thing with this model is that it defines two contact multiplier fields one on each body, which is redundant in the case of ring1-ring2 contact, and theoretically makes the problem not well defined, because there are infinite ways to split the contact traction between the two multiplier fields. But it works quite well in practice, nevertheless.

9

Dear Konstantinos,

Thanks to the code you provided, I was able to adapt it for this example shown in the image. My goal was to define the inner surface of the upper cylinder as a self-contact surface. However, when I set it up this way, the code fails, and I am not sure why. Do you have any advice?

Thank you in advance!

Additionally, I have an unrelated question. Could you provide me with a resource on time-dependent problems?. I took the time to read what is suggested in the documentation but I need more infos.

selfcontact
selfcontact476Ă—425 233 KB

10

well, in what sense does it “fail”? To begin with, what kind of elements are you using? Is it linear tetras or quadratic tetras? Try to use quadratic hexas, either 20 node (.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(3)"))) or 27-node (.set_classical_fem(2)). If you insist on using tetras, please at least avoid using linear ones.

11

Thank you for you reply I would try your suggestions.
here is actually the code. I really need tetras.

#!/usr/bin/env python
# -*- coding: UTF8 -*-
#
############################################################################
import numpy as np
import getfem as gf

gf.util_trace_level(1)
gf.util_warning_level(1)
PI = np.pi

ro = 24    # external radius. 
ri = 23.2   # internal radius
r = 0.8    # thickness
Length = 25  # length
h = 1.

Nr = 2           # Subdivisions in radial 
Ntheta = 2       # Subdivisions in angular
Nz = 20           # Subdivisions along z (axe z)


E1 = 1e5      # Young's modulus of solid 1
nu1 = 0.3    # Poisson's ratio of solid 1

E2 = 1e4      # Young's modulus of solid 1
nu2 = 0.3    # Poisson's ratio of solid 1

mult_fem_order = 1  # contact multipliers finite element order
integration_degree = 5         # 27 gauss points per hexa
integration_contact_degree = 7 # 16 gauss points per quad

c_aug = 1e-2*E1 # Augmentation parameter
release_dist = r



# create mesh
# rotation.
axis = np.array([0, 1, 0])
angle = PI/2
rotationMat = np.array([
    [np.cos(angle), 0.0, np.sin(angle)],
    [0, 1, 0],
    [-np.sin(angle), 0.0, np.cos(angle)]
])
# Import gmsh mesh.

m1 = gf.Mesh('import', 'gmsh', 'tube.msh')
m2 = gf.Mesh('import', 'gmsh', 'tube1.msh')
m2.translate([-50,ro*2+5,-50])
m2.transform(rotationMat)

print(m1.regions())

m1.optimize_structure()
m1.export_to_vtu("tube.vtu")
m2.optimize_structure()
m2.export_to_vtu("tube1.vtu")

regionM1, regionM2 = m1.regions(), m2.regions()
'''region_id = regionM1[0]
print(region_id)'''
# Regions 
EXTERIOR_Tu1 = regionM1[0]
INTERIOR_Tu2 = regionM2[3]
EXTERIOR_Tu2 = regionM2[0]
HANDLE1_Tu2 = regionM2[1]
HANDLE1_Tu1 = regionM1[1]
HANDLE2_Tu1 = regionM1[4]

# Finite element spaces and integration methods

mf1 = gf.MeshFem(m1, 3)
mf1.set_fem(gf.Fem("FEM_PK(3,1)"))

mf2 = gf.MeshFem(m2, 3)
mf2.set_fem(gf.Fem("FEM_PK(3,1)"))

mfmult1 = gf.MeshFem(m1, 3)
mfmult1.set_classical_fem(mult_fem_order)

mfmult2 = gf.MeshFem(m2, 3)
mfmult2.set_classical_fem(mult_fem_order)

#mfmult3 = gf.MeshFem(m2, 3)
#mfmult3.set_classical_fem(mult_fem_order) # for self contact

mim1 = gf.MeshIm(m1, integration_degree)
mim2 = gf.MeshIm(m2, integration_degree)
mim_contact1 = gf.MeshIm(m1, integration_contact_degree)
mim_contact2 = gf.MeshIm(m2, integration_contact_degree)
#mim_contact3 = gf.MeshIm(m2, integration_contact_degree) # for self contact


# Model object

md = gf.Model("real")

md.add_fem_variable("u1", mf1)
md.add_fem_variable("u2", mf2)
md.add_filtered_fem_variable("contact_mult1", mfmult1, 91)
md.add_filtered_fem_variable("contact_mult2", mfmult2, 911)
#md.add_filtered_fem_variable("contact_mult3", mfmult3, 933) # for self contact

# Elasticity
md.add_initialized_data("K1", E1/(3*(1-2*nu1)))
md.add_initialized_data("G1", E1/(2*(1+nu1)))
md.add_initialized_data("K2", E2/(3*(1-2*nu2)))
md.add_initialized_data("G2", E2/(2*(1+nu2)))
md.add_macro("F(u)", "Id(3)+Grad_u")
md.add_macro("tau(u)", "K1*log(Det(F(u)))*Id(3)+G1*pow(Det(F(u)),-2/3)*Deviator(F(u)*F(u)')")
md.add_macro("tau2(u)", "K2*log(Det(F(u)))*Id(3)+G2*pow(Det(F(u)),-2/3)*Deviator(F(u)*F(u)')")
md.add_nonlinear_term(mim1, "(tau(u1)*Inv(F(u1)')):Grad_Test_u1")
md.add_nonlinear_term(mim2, "(tau2(u2)*Inv(F(u2)')):Grad_Test_u2")

# Lagrange multipliers for handles
md.add_variable("q11",3) # tube 1 handle 1
md.add_variable("q12",3) # tube 1 handle 2
md.add_variable("q22",3) # tube 2 handle 2
md.add_initialized_data("p11", [0,0,0]) # we will keep this fixed
md.add_initialized_data("p12", [0,0,0]) # we will keep this fixed only at the beginning
md.add_initialized_data("p22", [0,0,0]) # we will drive this till the end

md.add_linear_term(mim1, "(u1-p11).q11", 95)
md.add_linear_term(mim1, "(u1-p12).q12", 92)
md.add_linear_term(mim2, "(u2-p22).q22", 912)

md.add_initialized_data("penal", 0.)
md.add_linear_term(mim1, "penal*(u1-p11).Test_u1", 95) # add some elastic support to handle 1 of tube 1 to avoid twisting
md.add_linear_term(mim1, "penal*(u1-p12).Test_u1", 92) # add some elastic support to handle 2 of tube 1 to avoid twisting
md.add_linear_term(mim2, "1e4*(u2-p22).Test_u2", 912) # add some elastic support to handle 2 of tube 2 to avoid twisting

# Add contact terms to the model

md.add_initialized_data("c", c_aug)

md.add_initialized_data("alpha", 1.)
md.add_initialized_data("f", 0.)
ibc = md.add_integral_large_sliding_contact_brick_raytracing("c", release_dist, "f", "alpha", 1) # 1 stands for the more symmetric version
md.add_master_slave_contact_boundary_to_large_sliding_contact_brick(ibc, mim_contact1,91, "u1", "contact_mult1", "")
md.add_master_slave_contact_boundary_to_large_sliding_contact_brick(ibc, mim_contact2, 911, "u2", "contact_mult2", "")
#md.add_master_slave_contact_boundary_to_large_sliding_contact_brick(ibc, mim_contact2, 933, "u2", "contact_mult2", "") for self contact

# Perform the simulation

print("Model dofs: %i" % md.nbdof())
print("Displacement fem dofs: %i" % (mf1.nbdof()+mf2.nbdof()))
#print("Contact mult dofs: %i" % mf_contact_mult.nbdof())

for step in range(90):
    t = step/30
    md.set_variable("p22", [0, -10*t, 0])
    print("Step", step, "t =", t)
    md.solve("noisy", "max_iter", 20, "max_res", 1e-8,
            "lsearch", "simplest", "alpha max ratio", 3., "alpha min", 0.02, "alpha mult", 0.6)
    mf1.export_to_vtu(f"Tube1/tube1_{step}.vtu", mf1, md.variable("u1"), "Displacements")
    mf2.export_to_vtu(f"Tube1/tube2_{step}.vtu", mf2, md.variable("u2"), "Displacements")

unfortunately I cannot upload the mesh file.

Thanks.

12

Linear tetras should not be used for any FE-modelling. They are just too bad. You need to rethink your strategy before continuing.

13

Thanks for your advised, I was able to solve the problem, I think my problem was on the definition of contact surfaces.
Thank you very much.

1000266173
10002661731220Ă—751 55.4 KB

1 Like