MPM versus FEM simulation of beam bending

FEM_beam_bending.edp

 // Macro
 //Macros for the gradient of a vector field (u1, u2)
 macro grad11(u1, u2) (dx(u1)) //
 macro grad21(u1, u2) (dy(u1)) //
 macro grad12(u1, u2) (dx(u2)) //
 macro grad22(u1, u2) (dy(u2)) //
  
 //Macros for the deformation gradient
 macro F11(u1, u2) (1.0 + grad11(u1, u2)) //
 macro F12(u1, u2) (0.0 + grad12(u1, u2)) //
 macro F21(u1, u2) (0.0 + grad21(u1, u2)) //
 macro F22(u1, u2) (1.0 + grad22(u1, u2)) //
  
 //Macros for the incremental deformation gradient
 macro dF11(varu1, varu2) (grad11(varu1, varu2)) //
 macro dF12(varu1, varu2) (grad12(varu1, varu2)) //
 macro dF21(varu1, varu2) (grad21(varu1, varu2)) //
 macro dF22(varu1, varu2) (grad22(varu1, varu2)) //
 
 //Macro for the determinant of the deformation gradient
 macro J(u1, u2) (
       F11(u1, u2)*F22(u1, u2)
     - F12(u1, u2)*F21(u1, u2)
 ) //
 
 //Macros for the inverse of the deformation gradient
 macro Finv11 (u1, u2) (
       F22(u1, u2) / J(u1, u2)
 ) //
 macro Finv22 (u1, u2) (
       F11(u1, u2) / J(u1, u2)
 ) //
 macro Finv12 (u1, u2) (
     - F12(u1, u2) / J(u1, u2)
 ) //
 macro Finv21 (u1, u2) (
     - F21(u1, u2) / J(u1, u2)
 ) //
 
 //Macros for the square of the inverse of the deformation gradient
 macro FFinv11 (u1, u2) (
      Finv11(u1, u2)^2
     + Finv12(u1, u2)*Finv21(u1, u2)
 ) //
 
 macro FFinv12 (u1, u2) (
       Finv12(u1, u2)*(Finv11(u1, u2)
     + Finv22(u1, u2))
 ) //
 
 macro FFinv21 (u1, u2) (
       Finv21(u1, u2)*(Finv11(u1, u2)
     + Finv22(u1, u2))
 ) //
 
 macro FFinv22 (u1, u2) (
       Finv12(u1, u2)*Finv21(u1, u2)
     + Finv22(u1, u2)^2
 ) //
 
 //Macros for the inverse of the transpose of the deformation gradient
 macro FinvT11(u1, u2) (Finv11(u1, u2)) //
 macro FinvT12(u1, u2) (Finv21(u1, u2)) //
 macro FinvT21(u1, u2) (Finv12(u1, u2)) //
 macro FinvT22(u1, u2) (Finv22(u1, u2)) //
 
 //The left Cauchy-Green strain tensor
 macro B11(u1, u2) (
       F11(u1, u2)^2 + F12(u1, u2)^2
 )//
 
 macro B12(u1, u2) (
       F11(u1, u2)*F21(u1, u2)
     + F12(u1, u2)*F22(u1, u2)
 )//
 
 macro B21(u1, u2) (
       F11(u1, u2)*F21(u1, u2)
     + F12(u1, u2)*F22(u1, u2)
 )//
 
 macro B22(u1, u2)(
       F21(u1, u2)^2 + F22(u1, u2)^2
 )//
 
 //The macros for the auxiliary tensors (D0, D1, D2, ...): Begin
 //The tensor quantity D0 = F{n} (delta F{n+1})^T
 macro d0Aux11 (u1, u2, varu1, varu2) (
       dF11(varu1, varu2) * F11(u1, u2)
     + dF12(varu1, varu2) * F12(u1, u2)
 )//
 
 macro d0Aux12 (u1, u2, varu1, varu2) (
       dF21(varu1, varu2) * F11(u1, u2)
     + dF22(varu1, varu2) * F12(u1, u2)
 )//
 
 macro d0Aux21 (u1, u2, varu1, varu2) (
       dF11(varu1, varu2) * F21(u1, u2)
     + dF12(varu1, varu2) * F22(u1, u2)
)//
 
macro d0Aux22 (u1, u2, varu1, varu2) (
      dF21(varu1, varu2) * F21(u1, u2)
    + dF22(varu1, varu2) * F22(u1, u2)
)//

///The tensor quantity D1 = D0 + D0^T
macro d1Aux11 (u1, u2, varu1, varu2) (
      2.0 * d0Aux11 (u1, u2, varu1, varu2)
)//

macro d1Aux12 (u1, u2, varu1, varu2) (
      d0Aux12 (u1, u2, varu1, varu2)
    + d0Aux21 (u1, u2, varu1, varu2)
)//

macro d1Aux21 (u1, u2, varu1, varu2) (
      d1Aux12 (u1, u2, varu1, varu2)
)//

macro d1Aux22 (u1, u2, varu1, varu2)(
      2.0 * d0Aux22 (u1, u2, varu1, varu2)
)//

///The tensor quantity D2 = F^{-T}_{n} dF_{n+1}
macro d2Aux11 (u1, u2, varu1, varu2) (
      dF11(varu1, varu2) * FinvT11(u1, u2)
    + dF21(varu1, varu2) * FinvT12(u1, u2)
)//

macro d2Aux12 (u1, u2, varu1, varu2) (
      dF12(varu1, varu2) * FinvT11(u1, u2)
    + dF22(varu1, varu2) * FinvT12(u1, u2)
)//

macro d2Aux21 (u1, u2, varu1, varu2) (
      dF11(varu1, varu2) * FinvT21(u1, u2)
    + dF21(varu1, varu2) * FinvT22(u1, u2)
)//

macro d2Aux22 (u1, u2, varu1, varu2) (
      dF12(varu1, varu2) * FinvT21(u1, u2)
    + dF22(varu1, varu2) * FinvT22(u1, u2)
)//
 
///The tensor quantity D3 = F^{-2}_{n} dF_{n+1}
macro d3Aux11 (u1, u2, varu1, varu2, w1, w2) (
      dF11(varu1, varu2) * FFinv11(u1, u2) * grad11(w1, w2)
    + dF21(varu1, varu2) * FFinv12(u1, u2) * grad11(w1, w2)
    + dF11(varu1, varu2) * FFinv21(u1, u2) * grad12(w1, w2)
    + dF21(varu1, varu2) * FFinv22(u1, u2) * grad12(w1, w2)
)//
 
macro d3Aux12 (u1, u2, varu1, varu2, w1, w2) (
      dF12(varu1, varu2) * FFinv11(u1, u2) * grad11(w1, w2)
    + dF22(varu1, varu2) * FFinv12(u1, u2) * grad11(w1, w2)
    + dF12(varu1, varu2) * FFinv21(u1, u2) * grad12(w1, w2)
    + dF22(varu1, varu2) * FFinv22(u1, u2) * grad12(w1, w2)
)//

macro d3Aux21 (u1, u2, varu1, varu2, w1, w2) (
      dF11(varu1, varu2) * FFinv11(u1, u2) * grad21(w1, w2)
    + dF21(varu1, varu2) * FFinv12(u1, u2) * grad21(w1, w2)
    + dF11(varu1, varu2) * FFinv21(u1, u2) * grad22(w1, w2)
     + dF21(varu1, varu2) * FFinv22(u1, u2) * grad22(w1, w2)
)//

macro d3Aux22 (u1, u2, varu1, varu2, w1, w2) (
      dF12(varu1, varu2) * FFinv11(u1, u2) * grad21(w1, w2)
     + dF22(varu1, varu2) * FFinv12(u1, u2) * grad21(w1, w2)
    + dF12(varu1, varu2) * FFinv21(u1, u2) * grad22(w1, w2)
    + dF22(varu1, varu2) * FFinv22(u1, u2) * grad22(w1, w2)
)//

///The tensor quantity D4 = (grad w) * Finv
macro d4Aux11 (w1, w2, u1, u2) (
      Finv11(u1, u2)*grad11(w1, w2)
    + Finv21(u1, u2)*grad12(w1, w2)
)//

macro d4Aux12 (w1, w2, u1, u2) (
      Finv12(u1, u2)*grad11(w1, w2)
    + Finv22(u1, u2)*grad12(w1, w2)
)//

macro d4Aux21 (w1, w2, u1, u2) (
      Finv11(u1, u2)*grad21(w1, w2)
    + Finv21(u1, u2)*grad22(w1, w2)
)//

macro d4Aux22 (w1, w2, u1, u2) (
      Finv12(u1, u2)*grad21(w1, w2)
    + Finv22(u1, u2)*grad22(w1, w2)
)//
//The macros for the auxiliary tensors (D0, D1, D2, ...): End

//The macros for the various stress measures: BEGIN
//The Kirchhoff stress tensor
macro StressK11(u1, u2) (                   
      mu * (B11(u1, u2) - 1.0) + 2 * la * (J(u1,u2)- 1)*J(u1,u2)
)//

//The Kirchhoff stress tensor
macro StressK12(u1, u2) (
      mu * B12(u1, u2)
)//

//The Kirchhoff stress tensor
macro StressK21(u1, u2) (
       mu * B21(u1, u2)
)//

//The Kirchhoff stress tensor
macro StressK22(u1, u2) (
      mu * (B22(u1, u2) - 1.0) + 2 * la * (J(u1,u2)- 1)*J(u1,u2)
)//

//The tangent Kirchhoff stress tensor
macro TanK11(u1, u2, varu1, varu2) (
      mu * d1Aux11(u1, u2, varu1, varu2)  + 2 * la * (2*J(u1,u2) - 1) * J(u1,u2) * d2Aux11(u1,u2,varu1,varu2)
)//

macro TanK12(u1, u2, varu1, varu2) (
      mu * d1Aux12(u1, u2, varu1, varu2)  + 2 * la * (2*J(u1,u2) - 1) * J(u1,u2) * d2Aux12(u1,u2,varu1,varu2)
)//

macro TanK21(u1, u2, varu1, varu2) (
      mu * d1Aux21(u1, u2, varu1, varu2)  + 2 * la * (2*J(u1,u2) - 1) * J(u1,u2) * d2Aux21(u1,u2,varu1,varu2)
)//
macro TanK22(u1, u2, varu1, varu2) (
      mu * d1Aux22(u1, u2, varu1, varu2)  + 2 * la * (2*J(u1,u2) - 1) * J(u1,u2) * d2Aux22(u1,u2,varu1,varu2)
)//
//The macros for the stress tensor components: END

// Parameters
//real mu = 5.e2; //Elastic coefficients (kg/cm^2)
real D = 1.e3; //1.e3; //(1 / compressibility)
real Pa = 0.0; //Stress loads   -3.e2

 real E = 500;   // rajouté...
 real nu = 0.35; //0.4 -0.2;   //rajouté...
 real mu = E/(2*(1+nu));
 real la = E * nu / ((1 + nu) * (1 - 2 * nu)) ;
//la=0.0;

cout << "mu= "<<mu << " lambda= "<< la <<endl;
 
real tol = 1.e-4; //Tolerance (L^2)

int m = 40, n = 20;
real f = -15;    

mesh Th = square(80, 40, [0.6*x, 0.1*y + 0.8]);
cout << "premier plot"<<endl;
plot(Th, wait=true, cmm="TOUT PREMIER DESSIN");

// Fespace
fespace Wh(Th, P1dc);  //P1dc
fespace Vh(Th, [P1, P1]);
Vh [w1, w2], [u1n,u2n], [varu1, varu2];
Vh [ehat1x, ehat1y], [ehat2x, ehat2y];
Vh [auxVec1, auxVec2]; //The individual elements of the total 1st Piola-Kirchoff stress
Vh [ef1, ef2];

fespace Sh(Th, P1);
Sh p, ppp;
Sh StrK11, StrK12, StrK21, StrK22;
Sh u1, u2;

// Problem
varf vfMass1D(p, q) = int2d(Th)(p*q);
matrix Mass1D = vfMass1D(Sh, Sh, solver=CG);
 
p[] = 1; //1
ppp[] = Mass1D * p[];

real DomainMass = ppp[].sum;
cout << "DomainMass = " << DomainMass << endl;

varf vmass ([u1, u2], [v1, v2], solver=CG)
    = int2d(Th)( (u1*v1 + u2*v2) / DomainMass );

matrix Mass = vmass(Vh, Vh);

matrix Id = vmass(Vh, Vh);

//Define the standard Euclidean basis functions
[ehat1x, ehat1y] = [1.0, 0.0];
[ehat2x, ehat2y] = [0.0, 1.0];

real ContParam, dContParam;
 
problem neoHookeanInc ([varu1, varu2], [w1, w2], solver=LU)
    = int2d(Th, qforder=1)(
        - (
              StressK11 (u1n, u2n) * d3Aux11(u1n, u2n, varu1, varu2, w1, w2)
            + StressK12 (u1n, u2n) * d3Aux12(u1n, u2n, varu1, varu2, w1, w2)
            + StressK21 (u1n, u2n) * d3Aux21(u1n, u2n, varu1, varu2, w1, w2)
            + StressK22 (u1n, u2n) * d3Aux22(u1n, u2n, varu1, varu2, w1, w2)
        )
        + TanK11 (u1n, u2n, varu1, varu2) * d4Aux11(w1, w2, u1n, u2n)
        + TanK12 (u1n, u2n, varu1, varu2) * d4Aux12(w1, w2, u1n, u2n)
         + TanK21 (u1n, u2n, varu1, varu2) * d4Aux21(w1, w2, u1n, u2n)
        + TanK22 (u1n, u2n, varu1, varu2) * d4Aux22(w1, w2, u1n, u2n)
    )
    + int2d(Th, qforder=1)(
          StressK11 (u1n, u2n) * d4Aux11(w1, w2, u1n, u2n)
        + StressK12 (u1n, u2n) * d4Aux12(w1, w2, u1n, u2n)
        + StressK21 (u1n, u2n) * d4Aux21(w1, w2, u1n, u2n)
        + StressK22 (u1n, u2n) * d4Aux22(w1, w2, u1n, u2n)
     )
     //Choose one of the following two boundary conditions involving Pa:
     // Load vectors normal to the boundary:
    - int1d(Th, 1)(              //1
          Pa * (w1*N.x + w2*N.y)   
    )
        - int2d(Th)(   //rajouté... gravité
         f*w2
    )
    //Load vectors tangential to the boundary:
    //- int1d(Th, 1)(
     //    Pa * (w1*N.y - w2*N.x)
     //)
     + on(4, varu1=0, varu2=0)  //2
    ;

//Auxiliary variables
matrix auxMat;

// Newton's method
 ContParam = 0.;
dContParam = 0.01;

//Initialization:
[varu1, varu2] = [0., 0.];
[u1n, u2n] = [0., 0.];
real res = 2. * tol;
real eforceres;
int loopcount = 0;
int loopmax = 400;//45;

// Iterations
while (loopcount <= loopmax && res >= tol){
//while (res >= tol){      
    loopcount ++;
    cout << "Loop " << loopcount << endl;

    // Solve
    neoHookeanInc;

    // Update
    u1 = varu1;
    u2 = varu2;

    // Residual
    w1[] = Mass*varu1[];
    res = sqrt(w1[]' * varu1[]); //L^2 norm of [varu1, varu2]
    cout << " L^2 residual = " << res << endl;

    // Newton
     u1n[] += varu1[];

    // Plot
    plot([u1n,u2n], cmm="displacement");
}

// Plot
plot(Th, wait=true, cmm="PREMIER PLOT");

// Movemesh
mesh Th1 = movemesh(Th, [x+u1n, y+u2n]);

// Plot
plot(Th1, wait=true, cmm="DEUXIEME PLOT");
//plot([u1n,u2n], cmm="TROISIEME PLOT");

// write mesh vertices to file 
ofstream ff("fichier_plot_beam_FEM_neohookean.txt");
int NbVertices = Th1.nv;
cout << "Number of vertices = " << NbVertices << endl;
for (int i = 0; i < NbVertices; i++)
      ff << Th1(i).x << " " << Th1(i).y << endl;

MPM-FEM-beam.ipynb

MPM_beam_bending.py

import numpy as np

import taichi as ti

ti.init(arch=ti.cpu)  # Try to run on GPU

quality = 1  # Use a larger value for higher-res simulations
n_particles, n_grid = 100000 * quality**2, 128 * quality
dx, inv_dx = 1 / n_grid, float(n_grid)
dt = 0.1e-4 / quality   #1e-4
elapsed_time =ti.field(dtype=float, shape=())
n_substeps =ti.field(dtype=int, shape=())
indice_temoin =ti.field(dtype=int, shape=())
p_vol, p_rho = (dx * 0.5)**2, 1
p_mass = p_vol * p_rho
E, nu = 500, 0.35  # Young's modulus and Poisson's ratio  0.1e5, 0.2

mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / (
    (1 + nu) * (1 - 2 * nu))  # Lame parameters
x = ti.Vector.field(2, dtype=float, shape=n_particles)  # position
v = ti.Vector.field(2, dtype=float, shape=n_particles)  # velocity
C = ti.Matrix.field(2, 2, dtype=float,
                    shape=n_particles)  # affine velocity field
F = ti.Matrix.field(2, 2, dtype=float,
                    shape=n_particles)  # deformation gradient                 
material = ti.field(dtype=int, shape=n_particles)  # material id
Jp = ti.field(dtype=float, shape=n_particles)  # plastic deformation
grid_v = ti.Vector.field(2, dtype=float,
                         shape=(n_grid, n_grid))  # grid node momentum/velocity

grid_m = ti.field(dtype=float, shape=(n_grid, n_grid))  # grid node mass

# geometrical dimensions of beam 
length_beam = 0.6; 
width_beam = 0.10
point_bas = 0.8  
indice_bas = int(point_bas * n_grid)
indice_haut = int(indice_bas + width_beam * n_grid)
gamma = 0.9999  # damping coefficient
g = 15 

@ti.kernel
def substep():
    for i, j in grid_m:
        grid_v[i, j] = [0, 0]
        grid_m[i, j] = 0
    for p in x:  # Particle state update and scatter to grid (P2G)
        base = (x[p] * inv_dx - 0.5).cast(int)
        fx = x[p] * inv_dx - base.cast(float)
        # Quadratic kernels  [http://mpm.graphics   Eqn. 123, with x=fx, fx-1,fx-2]
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]
        F[p] = (ti.Matrix.identity(float, 2) +
                dt * C[p]) @ F[p]  # deformation gradient update            
        h = ti.exp( 10 * (1.0 - Jp[p]))  # Hardening coefficient: snow gets harder when compressed
        if material[p] == 1:  # jelly, make it softer
            h = 0.3
        h = 1.0 # no change
        mu, la = mu_0 * h, lambda_0 * h       
        if material[p] == 0:  # liquid
            mu = 0.0
        U, sig, V = ti.svd(F[p])
        J = 1.0
        for d in ti.static(range(2)):
            new_sig = sig[d, d]
            if material[p] == 2:  # Snow
                new_sig = min(max(sig[d, d], 1 - 2.5e-2),
                              1 + 4.5e-3)  # Plasticity
            Jp[p] *= sig[d, d] / new_sig
            sig[d, d] = new_sig
            J *= new_sig
        if material[p] == 0:  # Reset deformation gradient to avoid numerical instability
            F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J)
        elif material[p] == 2:
            F[p] = U @ sig @ V.transpose(
            )  # Reconstruct elastic deformation gradient after plasticity
        stress = 2 * mu * (F[p] - U @ V.transpose()) @ F[p].transpose(
        ) + ti.Matrix.identity(float, 2) * la * J * (J - 1)                     
        stress = (-dt * p_vol * 4 * inv_dx * inv_dx) * stress
        affine = stress + p_mass * C[p]                                         
        for i, j in ti.static(ti.ndrange(
                3, 3)):  # Loop over 3x3 grid node neighborhood
            offset = ti.Vector([i, j])
            dpos = (offset.cast(float) - fx) * dx
            weight = w[i][0] * w[j][1]
            grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos)   
            grid_m[base + offset] += weight * p_mass

    for i, j in grid_m:        
        if grid_m[i, j] > 0:  # No need for epsilon here
            grid_v[i,
                   j] = (1 / grid_m[i, j]) * grid_v[i,
                                                    j]  # Momentum to velocity
            grid_v[i, j][1] -= dt * g  # gravity 50
            if i < 3 and grid_v[i, j][0] < 0:
                grid_v[i, j][0] = 0  # Boundary conditions
            if i > n_grid - 3 and grid_v[i, j][0] > 0: grid_v[i, j][0] = 0
            if j < 3 and grid_v[i, j][1] < 0: grid_v[i, j][1] = 0
            if j > n_grid - 3 and grid_v[i, j][1] > 0: grid_v[i, j][1] = 0
            
            # definition of beam
            if i < 3 and (j >= indice_bas and j <= indice_haut):
                grid_v[i, j][0] = 0
                grid_v[i, j][1] = 0
    for p in x:  # grid to particle (G2P)
        base = (x[p] * inv_dx - 0.5).cast(int)
        fx = x[p] * inv_dx - base.cast(float)
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1.0)**2, 0.5 * (fx - 0.5)**2]
        new_v = ti.Vector.zero(float, 2)
        new_C = ti.Matrix.zero(float, 2, 2)
        for i, j in ti.static(ti.ndrange(
                3, 3)):  # loop over 3x3 grid node neighborhood
            dpos = ti.Vector([i, j]).cast(float) - fx
            g_v = grid_v[base + ti.Vector([i, j])]
            weight = w[i][0] * w[j][1]
            new_v += weight * g_v * gamma
            new_C += 4 * inv_dx * weight * g_v.outer_product(dpos)
        v[p], C[p] = new_v, new_C
        x[p] += dt * v[p]  # advection
    old_elapsed_time = elapsed_time[None]
    elapsed_time[None] += dt
    n_substeps[None] += 1
    if int(elapsed_time[None]/0.1) > int(old_elapsed_time/0.1) :
        print('elapsed time = ',elapsed_time[None])
    

@ti.kernel
def initialize():
    for i in range(n_particles):
        x[i]=[ti.random()*length_beam, ti.random()*width_beam + point_bas]
        if x[i][0] > 0.98 * length_beam and indice_temoin[None] == 0: # for printing evolution of one particle close to the end of beam
            indice_temoin[None] = i
        material[i] = 1  # 0: fluid 1: jelly 2: snow        
        v[i] = ti.Matrix([0, 0])
        F[i] = ti.Matrix([[1, 0], [0, 1]])
        Jp[i] = 1
    elapsed_time[None] = 0.0
    n_substeps[None] = 0


initialize()
print("info: oscillations will be small enough after elapsed time > 0.5; you can stop before by typing escape on beam image")
gui = ti.GUI("Taichi MLS-MPM-99", res=1024, background_color=0x112F41)
limit_range = 30 #int(max(30,2e-3 // dt))

while (not gui.get_event(ti.GUI.ESCAPE, ti.GUI.EXIT)) and  (elapsed_time[None] < 0.6):    
    for s in range(limit_range):   #range(int(2e-3 // dt))
        substep()
        if n_substeps[None]%400 == 0 :
            print("nb_steps = ", n_substeps[None], ", particule_temoin X = ", x[indice_temoin[None]][0], ", particule_temoin Y = ",x[indice_temoin[None]][1] )
    gui.circles(x.to_numpy(), radius=1.5, color=0x068587)
    gui.show(
    )  # Change to gui.show(f'{frame:06d}.png') to write images to disk

# opening file to write contour of beam
fich = open("fichier_plot_beam_MPM.txt", "w" )
scatter_plot_MPM = []
liste_p = []

n_pas = 1024
pas = 1/n_pas
for p in range(n_particles) :
    indice_i = int(x[p][0]/pas)
    liste_p.append([indice_i, x[p][1]])

# select contour particles
for i in range(n_pas):
    ligne_v = [el[1] for el in liste_p if el[0]==i ]
    if len(ligne_v) > 0 :
        min_y = min(ligne_v )
        max_y = max(ligne_v )
        scatter_plot_MPM.extend([[i*pas, min_y], [i*pas, max_y]]) 

# write contour of beam to file
for p in scatter_plot_MPM:
    fich.write(str(p[0]) + " " + str(p[1]) + "\n")
  
fich.close()