Skip to content

Instantly share code, notes, and snippets.

@djachk
Created Feb 16, 2022
Embed
What would you like to do?
MPM versus FEM simulation of beam bending
// 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;
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment