-
-
Save djachk/28007c9d37ef6e448719263a9b37ee14 to your computer and use it in GitHub Desktop.
MPM versus FEM simulation of beam bending
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| // 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; |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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