# Defining a piecewise parameter with Heaviside in GWFL

**URL:** <https://getfem.discourse.group/t/defining-a-piecewise-parameter-with-heaviside-in-gwfl/104>\
**Category:** General\
**Created:** [November 1, 2025, 2:34pm UTC](https://getfem.discourse.group/t/defining-a-piecewise-parameter-with-heaviside-in-gwfl/104 "2025-11-01T14:34:58Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Djoumessi](https://avatars.discourse-cdn.com/v4/letter/d/da6949/32.png) [@Djoumessi](https://getfem.discourse.group/u/Djoumessi)\
**Post date:** [November 1, 2025, 2:34pm UTC](https://getfem.discourse.group/t/defining-a-piecewise-parameter-with-heaviside-in-gwfl/104/1 "2025-11-01T14:34:58Z")

</div>

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:

![capture](https://global.discourse-cdn.com/free1/uploads/getfem/original/1X/660fd8f3da01884a0ab2b9c23580539463efdce0.png)

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

```auto
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

```auto
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

```auto
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,

---

<div class="post-metadata">

**Author:** ![Djoumessi](https://avatars.discourse-cdn.com/v4/letter/d/da6949/32.png) [@Djoumessi](https://getfem.discourse.group/u/Djoumessi)\
**Post date:** [November 1, 2025, 3:47pm UTC](https://getfem.discourse.group/t/defining-a-piecewise-parameter-with-heaviside-in-gwfl/104/2 "2025-11-01T15:47:25Z")

</div>

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:

 ![same](https://global.discourse-cdn.com/free1/uploads/getfem/original/1X/6d310a6b62ac042ced7bca089eb3eb3e852e06a5.png)

- 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:

 ![prendre](https://global.discourse-cdn.com/free1/uploads/getfem/original/1X/92ceccec8fad5f008dadf192fe2e451f264187b1.png)

- 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.

---

<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:** [November 3, 2025, 8:51am UTC](https://getfem.discourse.group/t/defining-a-piecewise-parameter-with-heaviside-in-gwfl/104/3 "2025-11-03T08:51:14Z")

</div>

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

```python
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])`.

---

<div class="post-metadata">

**Author:** ![Djoumessi](https://avatars.discourse-cdn.com/v4/letter/d/da6949/32.png) [@Djoumessi](https://getfem.discourse.group/u/Djoumessi)\
**Post date:** [November 3, 2025, 4:17pm UTC](https://getfem.discourse.group/t/defining-a-piecewise-parameter-with-heaviside-in-gwfl/104/4 "2025-11-03T16:17:00Z")

</div>

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.

 ![Screenshot from 2025-11-03 16-05-23](https://global.discourse-cdn.com/free1/uploads/getfem/original/1X/98f87ac8d85e99016ecff69248a82c6c3b9e367d.png)  
 ![Screenshot from 2025-11-03 17-09-10](https://global.discourse-cdn.com/free1/uploads/getfem/original/1X/43de13a78bc134fe6768f71af434b262e9450f7e.png)  
the gmsh code

```auto
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

---

<div class="post-metadata">

**Author:** ![Djoumessi](https://avatars.discourse-cdn.com/v4/letter/d/da6949/32.png) [@Djoumessi](https://getfem.discourse.group/u/Djoumessi)\
**Post date:** [January 9, 2026, 11:36am UTC](https://getfem.discourse.group/t/defining-a-piecewise-parameter-with-heaviside-in-gwfl/104/5 "2026-01-09T11:36:18Z")

</div>

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.

```auto
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")

```

 ![maze](https://global.discourse-cdn.com/free1/uploads/getfem/original/1X/415b3386cfde47036748b4217dfee311ea54f1b0.webp)