Defining a piecewise parameter with Heaviside in GWFL

[Djoumessi]

Hello everyone,

I’m trying to define a spatially varying shear modulus \mu(r) in a hollow cylinder using the GetFEM Weak Form Language (GWFL). The cylinder has an inner radius R_i​=2.0 and thickness 0.50.
I would like \mu to vary in the radial direction according to the following piecewise definition:

To achieve this, I used the following GWFL expression:
first approach

mu_expr = (
    "2.0 * (1 - Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.15))) + "
    "0.5 * (Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.15)) - "
    "Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.3))) + "
    "2.5 * Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.3))"
)

md1.add_macro("mu", mu_expr)

second approach

mu_expr = (
    "2.0 * (1 - Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.15))) + "
    "0.5 * (Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.15)) - "
    "Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.3))) + "
    "2.5 * Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.3))"
)
mf_mu = gf.MeshFem(m)
mf_mu.set_classical_fem(0)
md1.add_initialized_fem_data("mu", mf_mu, md1.interpolation(mu_expr, mf_mu)) 

third approach

mu_expr = (
    "2.0 * (1 - Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.15))) + "
    "1.0 * (Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.15)) - "
    "Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.3))) + "
    "2.5 * Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.3))"
)

mfp = gf.MeshFem(m, 1)  
mfp.set_fem(gf.Fem("FEM_PK(3,1)")) 

md1.add_initialized_fem_data("mu", mfp, md1.interpolation(mu_expr, mfp))

These somehow works, but the third approach perfom better. I’m not completely sure whether this is the most appropriate or numerically robust way to define such a piecewise-constant field in GetFEM.

Could you please tell me if this approach using Heaviside is correct and recommended for defining spatially varying material parameters?

Thank you very much for your help.

Best regards,

[Djoumessi]

I tested the third approach and performed an inflation test. I applied the same method with both GetFEM and FEniCS. The two images below show the distribution of \mu and the pressure, respectively.

For the \mu distribution:

• The GetFEM result is shown in the top figure
• The FEniCS result is shown in the bottom figure

Both results exhibit the same distribution.

For the pressure:

• The GetFEM result is shown in the top figure
• The FEniCS result is shown in the bottom figure

The pressure values match between both codes; however, the GetFEM pressure field appears asymmetric.

Could you please help me understand what might cause this asymmetry in GetFEM? Where could I have made a mistake?

Thank you in advance for your help.

[Konstantinos.Poulios]

hello, this is a good example, but please show us your mesh. I assume that your mesh is conformal with the boundaries R_i, R_i+0.15 and R_i+0.3.

Method 1 could work, but it is not needed to carry the entire expression in the macro every time you need mu.

Method 2 is correct, you can use that one. The expression md1.interpolation(mu_expr, mf_mu) will evaluate your shear modulus in the center of the 0-order element and assign this value to the entire element.

Method 3 is in principle not correct, because you are representing a discontinuous field on a continuous fem.

Method 4. If you use only 1 integration method (which will normally be the case), the best solution is to compute and store mu on the Gauss points with

mimd = gf.MeshImData(mim, -1, 1) # -1 is the region, 1 is for storing a scalar field
md.add_im_data("mu", mimd)
md.set_variable("mu", md.interpolation(mu_expr, mimd, -1))

Regarding the resulting pressure p what do you mean with this?

@Djoumessi

The pressure values match between both codes;

In paraview they do not seem to match. Please share your code and I will have a look at it.

Extra tip: For the expression sqrt(X(1)*X(1) + X(2)*X(2)) you could also use something like Norm(X.*[1,1,0]).

[Djoumessi]

Hi Konstantinos,

Here is the mesh I’m using. The Gmsh code is also available.
As for the GetFEM code, I’m currently extracting a minimal working example from my original project. I’ll share it as soon as I’m done.

the gmsh code

lc = 0.35; //characteristic length for the cell
Point(1) = {0,0,0,lc}; // center point (x,y,z,lc)
//outer ellipses
a = 2.5;
b = 2.5;
Point(2) = {a,0,0,lc};
Point(3) = {0,b,0,lc};
Point(4) = {-a,0,0,lc};
Point(5) = {0,-b,0,lc};
Ellipse(10) = {2,1,3,3};
Ellipse(11) = {3,1,4,4};
Ellipse(12) = {4,1,5,5};
Ellipse(13) = {5,1,2,2};

// inner ellipses
move = 0.00; //"move" center
Point(101) = {move,0,0,lc};
c = 2.0;
d = 2.0;
Point(6) = {c+move,0,0,lc};
Point(7) = {move,d,0,lc};
Point(8) = {-c+move,0,0,lc};
Point(9) = {move,-d,0,lc};
Ellipse(14) = {6,101,7,7};
Ellipse(15) = {7,101,8,8};
Ellipse(16) = {8,101,9,9};
Ellipse(17) = {9,101,6,6};
//Line(40) = {7,3};
//Line(41) = {6,2};
//Line(42) = {9,5};
//Line(43) = {8,4};

// connect lines of outer and inner ellipses to one
Line Loop(20) = {10,11,12,13}; // only outer
Line Loop(21) = {14,15,16,17}; // only inner

 

Plane Surface(32) = {20,21};

length = 3.5;
csf[] = Extrude(0,0,length){Surface{32};};

Physical Surface (92) = {74};
Physical Surface (93) = {73, 69, 65, 61}; // inner
Physical Surface (91) = {57, 53, 49, 45};// outer
Physical Volume ("vol") = {1};
//+
Physical Surface(95) = {32};
//Mesh.Algorithm = 8;
//Mesh 2;
//Mesh 3;

Thanks

[Djoumessi]

Dear Konstantinons,

Apologies for the delay. While reviewing our previous conversation, I realized that I had only sent the Gmsh code.

Please find attached a minimal working example (MWE) that I have prepared. I also took the opportunity to add an additional part, namely the construction of mechanical parameters from an image and their mapping onto the computational domain .

For this purpose, I relied on an example previously published in the discussion group, which I attempted to adapt to my specific case. The corresponding implementation is included directly in the code.

I apologize if some comments in the code are written in French.

import numpy as np
import getfem as gf
import os
import time
from PIL import Image, ImageOps
from scipy.ndimage import map_coordinates

start_time = time.time()

gf.util_trace_level(1)
gf.util_warning_level(1)

m = gf.Mesh('import', 'gmsh', 'GI.msh')
mf = gf.MeshFem(m, 1)
mf.set_fem(gf.Fem("FEM_PK(3,2)"))
mim = gf.MeshIm(m, 4) 


integration_degree = 3         # 27 gauss points per hexa

#mesh coordinate

coords = m.pts()

print(coords.shape)
print(coords)


#----------------------- from gray scale image --------------------
# --- Charger image ---
img = Image.open("maze.jpg").convert("L")
data = np.array(img)/255.0
NY, NX = data.shape
######################
md1 = gf.Model("real")
mf1 = gf.MeshFem(m, 1)
mf1.set_fem(gf.Fem("FEM_PK(3,2)"))
mfu = gf.MeshFem(m, m.dim())   # displacement
mfu.set_fem(gf.Fem("FEM_PK(3,2)")) # displacement
mfp = gf.MeshFem(m, 1)  # pressure
mfp.set_fem(gf.Fem("FEM_PK(3,1)")) # presssure

mim2 = gf.MeshIm(m, integration_degree)

md1.add_fem_variable("u", mfu) # displacement
md1.add_fem_variable('p', mfp) # pressure variable
# mechanical parameters
E = 10.0
nu = 0.3
mu = E / (2 * (1 + nu))
lmbda = 100*E * nu / ((1 + nu) * (1 - 2 * nu))

#Expression to project
mu_expr = (
    "2.0 * (1 - Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.15))) + "
    "0.5 * (Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.15)) - "
    "Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.3))) + "
    "2.5 * Heaviside(sqrt(X(1)*X(1) + X(2)*X(2)) - (2.0 + 0.3))"
)

#md1.add_macro("mu", mu_expr) # method 1
mfE = gf.MeshFem(m)
mfE.set_classical_fem(1)
md1.add_initialized_fem_data("mu", mfE, md1.interpolation(mu_expr, mfE)) # method 2
#md1.add_initialized_fem_data("mu", mfp, md1.interpolation(mu_expr, mfp)) # method 3

#mimd = gf.MeshImData(mim2, -1, 1) # -1 is the region, 1 is for storing a scalar field
#md1.add_im_data("mu", mimd)
#md1.set_variable("mu", md1.interpolation(mu_expr, mimd, -1)) # method 4
#md1.add_initialized_data("mu",2.5)
md1.add_initialized_data("pEpi", 0.0)

# for the projection
mf_mu = gf.MeshFem(m, 1)
mf_mu.set_fem(gf.Fem("FEM_PK(3,1)")) # mu
Muu = md1.interpolation(mu_expr, mfp) # use to visualize the Mu from expression



# ------------------- partant de l'image ------------------
# --- MeshFem pour rho ---
mfrho = gf.MeshFem(m, 1)
mfrho.set_classical_fem(0)
#mfrho.set_fem(gf.Fem("FEM_PK(3,1)")) 
md1.add_fem_data("rho", mfrho)

# --- Coordonnées ---
Xp = md1.interpolation("X(1)", mfrho)
Yp = md1.interpolation("X(2)", mfrho)
Zp = md1.interpolation("X(3)", mfrho)

# --- Coord. cylindriques ---
theta = np.arctan2(Yp, Xp)
theta[theta < 0] += 2*np.pi
R = np.sqrt(Xp**2 + Yp**2)
L = Zp.max() - Zp.min()

# --- Mapping ---
i = (theta / (2*np.pi) * (NX-1)).astype(int)
j = ((Zp - Zp.min()) / L * (NY-1)).astype(int)
i = np.clip(i, 0, NX-1)
j = np.clip(j, 0, NY-1)

rho = data[NY-1-j, i]

pts = mfrho.basic_dof_nodes()
#pts = m.pts()
Xp, Yp, Zp = pts[0, :], pts[1, :], pts[2, :]

# --- Convert ---
theta = np.arctan2(Yp, Xp)
theta[theta < 0] += 2 * np.pi
R = np.sqrt(Xp**2 + Yp**2)
L = Zp.max() - Zp.min()

# --- Normalisation to image index ---
xi = (theta / (2 * np.pi)) * (NX - 1)
yi = ((Zp - Zp.min()) / L) * (NY - 1)

# ---  map_coordinates  ---

rho = map_coordinates(data, [NY - 1 - yi, xi], order=1, mode='nearest')


md1.set_variable("rho", rho) # set variable

# Kinematics

md1.add_macro("F", "Id(3)+Grad_u")
md1.add_macro("J", "Det(F)") 
md1.add_macro("C", "F'*F")
md1.add_macro("B", "F*F'")
md1.add_macro("I1", "Trace(C)")

# Stress tensor
#md1.add_macro("sigma", "2*0.000035*rho*B - p*Id(3)")

md1.add_macro("sigma", "2*mu*B - p*Id(3)")

md1.add_initialized_data("stab", 100)
md1.add_macro("Piola", "J*sigma*Inv(F)'")
md1.add_nonlinear_term(mim2, "Piola:Grad_Test_u")
md1.add_nonlinear_term(mim2, "Test_p*(J-1)")
md1.add_nonlinear_term(mim2, "stab*(Grad_p.Grad_Test_p)")
md1.add_macro("normal2", "Normalized(Inv(F').Normal)")
md1.add_macro("normal3", "Inv(F').Normal")
md1.add_nonlinear_term(mim2, "-pEpi*Det(F)*normal2.Test_u", 93)

md1.add_initialized_data("end92", [0.0,0.0,0.0])
md1.add_initialized_data("end95", [0.0,0.0,0.0])
md1.add_Dirichlet_condition_with_multipliers(mim2, 'u', mfu, 92, 'end92')
md1.add_Dirichlet_condition_with_multipliers(mim2, 'u', mfu, 95, 'end95')

# Dossier de résultats
if not os.path.exists('outFW'):
    os.makedirs('outFW')
tf = 2.0
to = 0.0
n = 0
dt = 0.1
while to <= tf:
    print("time=: %i" % to)
    #md.set_variable("p22", [0, -0*t, 0])
    md1.set_variable("pEpi", -0.3*to)
    md1.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"internalS/flat_tube_{step}.vtu", mf1, md.variable("u1"), "Displacements")
    output = (mfu, md1.variable("u"), "u",
                  mfp, md1.variable("p"), "p",
                  mf_mu, Muu, "Mu",
                  mfE, md1.variable("mu"), "mu",
                  mfrho, md1.variable("rho"), "rho")

    mfu.export_to_vtu(f"outFW/EM_{n:04d}.vtu", *output)
    to += dt    
    n += 1
end_time = time.time()  # fin du chronométrage

print("Temps d'exécution :", end_time - start_time, "secondes")