diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index f8e1788c..afc95110 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-05T12:11:12","documenter_version":"1.5.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-05T13:31:18","documenter_version":"1.5.0"}} \ No newline at end of file diff --git a/dev/examples/dam_break/index.html b/dev/examples/dam_break/index.html index 4d036a45..5b73799c 100644 --- a/dev/examples/dam_break/index.html +++ b/dev/examples/dam_break/index.html @@ -355,4 +355,4 @@ # Extract the activated degrees of freedom submatrix(A, dofmap) -end

This page was generated using Literate.jl.

+end

This page was generated using Literate.jl.

diff --git a/dev/examples/elastic_impact/index.html b/dev/examples/elastic_impact/index.html index 743a54ec..017fea29 100644 --- a/dev/examples/elastic_impact/index.html +++ b/dev/examples/elastic_impact/index.html @@ -189,4 +189,4 @@ end end end -end

This page was generated using Literate.jl.

+end

This page was generated using Literate.jl.

diff --git a/dev/examples/implicit_jacobian_based/index.html b/dev/examples/implicit_jacobian_based/index.html index f80ac814..59e49bfe 100644 --- a/dev/examples/implicit_jacobian_based/index.html +++ b/dev/examples/implicit_jacobian_based/index.html @@ -72,7 +72,7 @@ end # Sparse matrix - A = create_sparse_matrix(Vec{3}, it, grid.X) + A = create_sparse_matrix(it, grid.X) # Outputs outdir = mkpath(joinpath("output", "implicit_jacobian_based")) @@ -174,4 +174,4 @@ end submatrix(A, dofmap) -end

This page was generated using Literate.jl.

+end

This page was generated using Literate.jl.

diff --git a/dev/examples/implicit_jacobian_free/index.html b/dev/examples/implicit_jacobian_free/index.html index 4e81a5cc..7b540743 100644 --- a/dev/examples/implicit_jacobian_free/index.html +++ b/dev/examples/implicit_jacobian_free/index.html @@ -182,4 +182,4 @@ @. JδU = δU - β*Δt^2 * $dofmap(grid.f) * $dofmap(grid.m⁻¹) end -end

This page was generated using Literate.jl.

+end

This page was generated using Literate.jl.

diff --git a/dev/examples/rigid_body_contact/index.html b/dev/examples/rigid_body_contact/index.html index a5c1260c..98de005d 100644 --- a/dev/examples/rigid_body_contact/index.html +++ b/dev/examples/rigid_body_contact/index.html @@ -186,4 +186,4 @@ σ = dev(σ) end σ -end

This page was generated using Literate.jl.

+end

This page was generated using Literate.jl.

diff --git a/dev/examples/tlmpm_vortex/index.html b/dev/examples/tlmpm_vortex/index.html index 12c54b08..2664ca72 100644 --- a/dev/examples/tlmpm_vortex/index.html +++ b/dev/examples/tlmpm_vortex/index.html @@ -158,4 +158,4 @@ end end end -end

This page was generated using Literate.jl.

+end

This page was generated using Literate.jl.

diff --git a/dev/getting_started/index.html b/dev/getting_started/index.html index 43ce99b5..38f00b58 100644 --- a/dev/getting_started/index.html +++ b/dev/getting_started/index.html @@ -100,4 +100,4 @@ aspect_ratio = :equal, legend = false, ) -end every 100Example block output +end every 100Example block output diff --git a/dev/index.html b/dev/index.html index 162c1bf5..5de686e7 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,2 +1,2 @@ -Home · Tesserae.jl
GitHub
+Home · Tesserae.jl
GitHub
diff --git a/dev/search_index.js b/dev/search_index.js index b34cd2b4..40355ce4 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"EditURL = \"../../literate/examples/tlmpm_vortex.jl\"","category":"page"},{"location":"examples/tlmpm_vortex/#Total-Lagrangian-MPM","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"","category":"section"},{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"This example demonstrates the total lagrangian material point method[1]. The implementation solves generalized vortex problem[1] using a linear kernel.","category":"page"},{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"note: Note\nCurrently, the Bernstein function used in the paper[1] has not been implemented.","category":"page"},{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"[1]: de Vaucorbeil, A., Nguyen, V.P. and Hutchinson, C.R., 2020. A Total-Lagrangian Material Point Method for solid mechanics problems involving large deformations. Computer Methods in Applied Mechanics and Engineering, 360, p.112783.","category":"page"},{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"using Tesserae\n\nfunction main()\n\n # Simulation parameters\n h = 0.02 # Grid spacing\n T = 1.0 # Time span\n CFL = 0.1 # Courant number\n α = 0.99 # PIC-FLIP parameter\n\n # Material constants\n E = 1e6 # Young's modulus\n ν = 0.3 # Poisson's ratio\n λ = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter\n μ = E / 2(1 + ν) # Shear modulus\n ρ⁰ = 1e3 # Initial density\n\n # Geometry\n Rᵢ = 0.75\n Rₒ = 1.25\n\n # Equations for vortex\n G = π\n R̄ = (Rᵢ + Rₒ) / 2\n function calc_b_Rθ(R, t)\n local h′′, h′, h = hessian(R -> 1-8((R-R̄)/(Rᵢ-Rₒ))^2+16((R-R̄)/(Rᵢ-Rₒ))^4, R, :all)\n local g′′, g′, g = hessian(t -> G*sin(π*t/T), t, :all)\n β = g * h\n b_R = ( μ/ρ⁰*(3g*h′+R*g*h′′) - R*g′′*h)*sin(β) + (μ/ρ⁰*R*(g*h′)^2 - R*(g′*h)^2)*cos(β)\n b_θ = (-μ/ρ⁰*(3g*h′+R*g*h′′) + R*g′′*h)*cos(β) + (μ/ρ⁰*R*(g*h′)^2 + R*(g′*h)^2)*sin(β)\n Vec(b_R, b_θ)\n end\n isinside(x::Vec) = Rᵢ^2 < x⋅x < Rₒ^2\n\n GridProp = @NamedTuple begin\n X :: Vec{2, Float64}\n m :: Float64\n m⁻¹ :: Float64\n mv :: Vec{2, Float64}\n fint :: Vec{2, Float64}\n fext :: Vec{2, Float64}\n b :: Vec{2, Float64}\n v :: Vec{2, Float64}\n vⁿ :: Vec{2, Float64}\n end\n ParticleProp = @NamedTuple begin\n x :: Vec{2, Float64}\n X :: Vec{2, Float64}\n m :: Float64\n V⁰ :: Float64\n v :: Vec{2, Float64}\n ṽ :: Vec{2, Float64}\n ã :: Vec{2, Float64}\n P :: SecondOrderTensor{2, Float64, 4}\n F :: SecondOrderTensor{2, Float64, 4}\n end\n\n # Background grid\n grid = generate_grid(GridProp, CartesianMesh(h, (-1.5,1.5), (-1.5,1.5)))\n outside_gridinds = findall(!isinside, grid.X)\n\n # Particles\n particles = generate_particles(ParticleProp, grid.X; alg=GridSampling(), spacing=1)\n particles.V⁰ .= volume(grid.X) / length(particles)\n\n filter!(pt->isinside(pt.x), particles)\n\n @. particles.X = particles.x\n @. particles.m = ρ⁰ * particles.V⁰\n @. particles.F = one(particles.F)\n @show length(particles)\n\n # Precompute linear kernel values\n mpvalues = generate_mpvalues(LinearBSpline(), grid.X, length(particles))\n for p in eachindex(particles, mpvalues)\n update!(mpvalues[p], particles.x[p], grid.X)\n end\n\n # Outputs\n outdir = mkpath(joinpath(\"output\", \"tlmpm_vortex\"))\n pvdfile = joinpath(outdir, \"paraview\")\n closepvd(openpvd(pvdfile)) # create file\n\n t = 0.0\n step = 0\n fps = 60\n savepoints = collect(LinRange(t, T, round(Int, T*fps)+1))\n\n Tesserae.@showprogress while t < T\n\n # Calculate timestep based on the wave speed\n vmax = maximum(@. sqrt((λ+2μ) / (particles.m/(particles.V⁰ * det(particles.F)))) + norm(particles.v))\n Δt = CFL * spacing(grid) / vmax\n\n # Compute grid body forces\n for i in eachindex(grid)\n if isinside(grid.X[i])\n (x, y) = grid.X[i]\n R = sqrt(x^2 + y^2)\n θ = atan(y, x)\n grid.b[i] = rotmat(θ) ⋅ calc_b_Rθ(R, t)\n end\n end\n\n # Particle-to-grid transfer\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * v[p]\n fint[i] = @∑ -V⁰[p] * P[p] ⋅ ∇w[ip]\n end\n\n # Update grid velocity\n @. grid.m⁻¹ = inv(grid.m) * !iszero(grid.m)\n @. grid.fext = grid.m * grid.b\n @. grid.vⁿ = grid.mv * grid.m⁻¹\n @. grid.v = grid.vⁿ + Δt * (grid.fint + grid.fext) * grid.m⁻¹\n grid.v[outside_gridinds] .= zero(eltype(grid.v))\n\n # Update particle velocity and position\n @G2P grid=>i particles=>p mpvalues=>ip begin\n ṽ[p] = @∑ w[ip] * v[i]\n ã[p] = @∑ w[ip] * (v[i] - vⁿ[i]) / Δt\n v[p] = (1-α)*ṽ[p] + α*(v[p] + Δt*ã[p])\n x[p] += Δt * ṽ[p]\n end\n\n # Remap updated velocity to grid (MUSL)\n @P2G grid=>i particles=>p mpvalues=>ip begin\n mv[i] = @∑ w[ip] * m[p] * v[p]\n v[i] = mv[i] * m⁻¹[i]\n end\n grid.v[outside_gridinds] .= zero(eltype(grid.v))\n\n # Update stress\n @G2P grid=>i particles=>p mpvalues=>ip begin\n F[p] += @∑ Δt * v[i] ⊗ ∇w[ip]\n P[p] = μ * (F[p] - inv(F[p])') + λ * log(det(F[p])) * inv(F[p])'\n end\n\n t += Δt\n step += 1\n\n if t > first(savepoints)\n popfirst!(savepoints)\n openpvd(pvdfile; append=true) do pvd\n openvtm(string(pvdfile, step)) do vtm\n angle(x) = atan(x[2], x[1])\n openvtk(vtm, particles.x) do vtk\n vtk[\"velocity\"] = particles.v\n vtk[\"initial angle\"] = angle.(particles.X)\n end\n openvtk(vtm, grid.X) do vtk\n vtk[\"external force\"] = grid.fext\n end\n pvd[t] = vtm\n end\n end\n end\n end\nend","category":"page"},{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"","category":"page"},{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"This page was generated using Literate.jl.","category":"page"},{"location":"examples/implicit_jacobian_based/","page":"Jacobian-based implicit method","title":"Jacobian-based implicit method","text":"EditURL = \"../../literate/examples/implicit_jacobian_based.jl\"","category":"page"},{"location":"examples/implicit_jacobian_based/#Jacobian-based-implicit-method","page":"Jacobian-based implicit method","title":"Jacobian-based implicit method","text":"","category":"section"},{"location":"examples/implicit_jacobian_based/","page":"Jacobian-based implicit method","title":"Jacobian-based implicit method","text":"using Tesserae\n\nfunction main()\n\n # Simulation parameters\n h = 0.05 # Grid spacing\n T = 1.0 # Time span\n g = 20.0 # Gravity acceleration\n Δt = 0.05 # Timestep\n\n # Material constants\n E = 1e6 # Young's modulus\n ν = 0.3 # Poisson's ratio\n λ = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter\n μ = E / 2(1 + ν) # Shear modulus\n ρ⁰ = 500.0 # Density\n\n # Newmark-beta integration\n β = 1/4\n γ = 1/2\n\n GridProp = @NamedTuple begin\n X :: Vec{3, Float64}\n m :: Float64\n m⁻¹ :: Float64\n v :: Vec{3, Float64}\n vⁿ :: Vec{3, Float64}\n mv :: Vec{3, Float64}\n a :: Vec{3, Float64}\n aⁿ :: Vec{3, Float64}\n ma :: Vec{3, Float64}\n u :: Vec{3, Float64}\n f :: Vec{3, Float64}\n end\n ParticleProp = @NamedTuple begin\n x :: Vec{3, Float64}\n m :: Float64\n V⁰ :: Float64\n v :: Vec{3, Float64}\n a :: Vec{3, Float64}\n b :: Vec{3, Float64}\n ∇u :: SecondOrderTensor{3, Float64, 9}\n F :: SecondOrderTensor{3, Float64, 9}\n ΔF⁻¹ :: SecondOrderTensor{3, Float64, 9}\n τ :: SecondOrderTensor{3, Float64, 9}\n c :: FourthOrderTensor{3, Float64, 81}\n end\n\n # Background grid\n grid = generate_grid(GridProp, CartesianMesh(h, (0.0,1.2), (0.0,2.0), (-0.2,0.2)))\n\n # Particles\n beam = Tesserae.Box((0,1), (0.85,1.15), (-0.15,0.15))\n particles = generate_particles(ParticleProp, grid.X; domain=beam, alg=GridSampling())\n particles.V⁰ .= volume(beam) / length(particles)\n @. particles.m = ρ⁰ * particles.V⁰\n @. particles.F = one(particles.F)\n @. particles.b = Vec(0,-g,0)\n @show length(particles)\n\n # Interpolation\n # Use the kernel correction to properly handle the boundary conditions\n it = KernelCorrection(QuadraticBSpline())\n mpvalues = generate_mpvalues(it, grid.X, length(particles))\n\n # Neo-Hookean model\n function kirchhoff_stress(F)\n J = det(F)\n b = symmetric(F ⋅ F')\n μ*(b-I) + λ*log(J)*I\n end\n\n # Sparse matrix\n A = create_sparse_matrix(Vec{3}, it, grid.X)\n\n # Outputs\n outdir = mkpath(joinpath(\"output\", \"implicit_jacobian_based\"))\n pvdfile = joinpath(outdir, \"paraview\")\n closepvd(openpvd(pvdfile)) # Create file\n\n t = 0.0\n step = 0\n\n Tesserae.@showprogress while t < T\n\n for p in eachindex(particles, mpvalues)\n update!(mpvalues[p], particles.x[p], grid.X)\n end\n\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * v[p]\n ma[i] = @∑ w[ip] * m[p] * a[p]\n end\n\n # Compute the grid velocity and acceleration at t = tⁿ\n @. grid.m⁻¹ = inv(grid.m) * !iszero(grid.m)\n @. grid.vⁿ = grid.mv * grid.m⁻¹\n @. grid.aⁿ = grid.ma * grid.m⁻¹\n\n # Update a dof map\n dofmask = trues(3, size(grid)...)\n for i in eachindex(grid)\n dofmask[:,i] .= !iszero(grid.m[i])\n end\n for i in eachindex(grid)[1,:,:]\n dofmask[:,i] .= false\n grid.vⁿ[i] = zero(Vec{3})\n grid.aⁿ[i] = zero(Vec{3})\n end\n dofmap = DofMap(dofmask)\n\n # Solve the nonlinear equation\n state = (; grid, particles, mpvalues, kirchhoff_stress, β, γ, A, dofmap, Δt)\n @. grid.u = zero(grid.u) # Set zero dispacement for the first guess of the solution\n U = copy(dofmap(grid.u)) # Convert grid data to plain vector data\n compute_residual(U) = residual(U, state)\n compute_jacobian(U) = jacobian(U, state)\n Tesserae.newton!(U, compute_residual, compute_jacobian)\n\n # Grid dispacement, velocity and acceleration have been updated during Newton's iterations\n @G2P grid=>i particles=>p mpvalues=>ip begin\n ∇u[p] = @∑ u[i] ⊗ ∇w[ip]\n a[p] = @∑ w[ip] * a[i]\n v[p] += @∑ Δt * w[ip] * ((1-γ)*a[p] + γ*a[i])\n x[p] = @∑ w[ip] * (X[i] + u[i])\n F[p] = (I + ∇u[p]) ⋅ F[p]\n end\n\n t += Δt\n step += 1\n\n openpvd(pvdfile; append=true) do pvd\n openvtk(string(pvdfile, step), particles.x) do vtk\n vtk[\"velocity\"] = particles.v\n vtk[\"von Mises\"] = @. vonmises(particles.τ / det(particles.F))\n pvd[t] = vtk\n end\n end\n end\nend\n\nfunction residual(U::AbstractVector, state)\n (; grid, particles, mpvalues, kirchhoff_stress, β, γ, dofmap, Δt) = state\n\n dofmap(grid.u) .= U\n @. grid.a = (1/(β*Δt^2))*grid.u - (1/(β*Δt))*grid.vⁿ - (1/2β-1)*grid.aⁿ\n @. grid.v = grid.vⁿ + Δt*((1-γ)*grid.aⁿ + γ*grid.a)\n\n transposing_tensor(σ) = @einsum (i,j,k,l) -> σ[i,l] * one(σ)[j,k]\n @G2P grid=>i particles=>p mpvalues=>ip begin\n # In addition to updating the stress tensor, the stiffness tensor,\n # which is utilized in the Jacobian-vector product, is also updated.\n ∇u[p] = @∑ u[i] ⊗ ∇w[ip]\n ΔF⁻¹[p] = inv(I + ∇u[p])\n c[p], τ[p] = gradient(∇u -> kirchhoff_stress((I + ∇u) ⋅ F[p]), ∇u[p], :all)\n c[p] = c[p] - transposing_tensor(τ[p] ⋅ ΔF⁻¹[p]')\n end\n @P2G grid=>i particles=>p mpvalues=>ip begin\n f[i] = @∑ -V⁰[p] * τ[p] ⋅ (∇w[ip] ⋅ ΔF⁻¹[p]) + w[ip] * m[p] * b[p]\n end\n\n @. $dofmap(grid.m) * $dofmap(grid.a) - $dofmap(grid.f)\nend\n\nfunction jacobian(U::AbstractVector, state)\n (; grid, particles, mpvalues, β, A, dofmap, Δt) = state\n\n I(i,j) = ifelse(i===j, one(Mat{3,3}), zero(Mat{3,3}))\n dotdot(a,C,b) = @einsum (i,j) -> a[k] * C[i,k,j,l] * b[l]\n @P2G_Matrix grid=>(i,j) particles=>p mpvalues=>(ip,jp) begin\n A[i,j] = @∑ (dotdot(∇w[ip] ⋅ ΔF⁻¹[p], c[p], ∇w[jp])) * V⁰[p] + 1/(β*Δt^2) * I(i,j) * m[p] * w[jp]\n end\n\n submatrix(A, dofmap)\nend","category":"page"},{"location":"examples/implicit_jacobian_based/","page":"Jacobian-based implicit method","title":"Jacobian-based implicit method","text":"","category":"page"},{"location":"examples/implicit_jacobian_based/","page":"Jacobian-based implicit method","title":"Jacobian-based implicit method","text":"This page was generated using Literate.jl.","category":"page"},{"location":"getting_started/#Getting-Started","page":"Getting Started","title":"Getting Started","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"using Tesserae\nimport Plots\n\n# Material constants\nE = 500 # Young's modulus\nν = 0.3 # Poisson's ratio\nλ = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter\nμ = E / 2(1 + ν) # Shear modulus\nρ = 1000 # Density\nr = 0.2 # Radius of disk\n\n# Properties for grid and particles\nstruct GridProp\n x :: Vec{2, Float64}\n m :: Float64\n mv :: Vec{2, Float64}\n f :: Vec{2, Float64}\n v :: Vec{2, Float64}\n vⁿ :: Vec{2, Float64}\nend\nstruct ParticleProp\n x :: Vec{2, Float64}\n m :: Float64\n V :: Float64\n v :: Vec{2, Float64}\n ∇v :: SecondOrderTensor{2, Float64, 4}\n σ :: SymmetricSecondOrderTensor{2, Float64, 3}\nend\n\n# Mesh\nmesh = CartesianMesh(0.05, (0,1), (0,1))\n\n# Background grid\ngrid = generate_grid(GridProp, mesh)\n\n# Particles\nparticles = let\n pts = generate_particles(ParticleProp, mesh; alg=GridSampling())\n pts.V .= volume(mesh) / length(pts)\n\n # Left disk\n lhs = findall(pts.x) do (x,y)\n (x-r)^2 + (y-r)^2 < r^2\n end\n\n # Right disk\n s = 1-r\n rhs = findall(pts.x) do (x,y)\n (x-s)^2 + (y-s)^2 < r^2\n end\n\n pts.v[lhs] .= Vec( 0.1, 0.1)\n pts.v[rhs] .= Vec(-0.1,-0.1)\n \n pts[[lhs; rhs]]\nend\n@. particles.m = ρ * particles.V\n\n# Interpolation\nmpvalues = [MPValue(LinearBSpline(), mesh) for _ in 1:length(particles)]\n\n# Plot results by `Plots.@gif`\nΔt = 0.001\nPlots.@gif for t in range(0, 4-Δt, step=Δt)\n\n # Update interpolation values\n for p in 1:length(particles)\n update!(mpvalues[p], particles[p], mesh)\n end\n\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * v[p]\n f[i] = @∑ -V[p] * σ[p] ⋅ ∇w[ip]\n vⁿ[i] = mv[i] / m[i]\n v[i] = vⁿ[i] + Δt * (f[i] / m[i])\n end\n\n @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] += @∑ w[ip] * (v[i] - vⁿ[i])\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n x[p] += @∑ Δt * (w[ip] * v[i])\n end\n\n for p in 1:length(particles)\n Δϵₚ = Δt * symmetric(particles.∇v[p])\n Δσₚ = λ*tr(Δϵₚ)*I + 2μ*Δϵₚ\n particles.V[p] *= 1 + tr(Δϵₚ)\n particles.σ[p] += Δσₚ\n end\n\n # plot results\n Plots.scatter(\n reinterpret(Tuple{Float64,Float64}, particles.x),\n lims = (0,1),\n ticks = 0:0.2:1,\n minorgrid = true,\n minorticks = 4,\n aspect_ratio = :equal,\n legend = false,\n )\nend every 100","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"EditURL = \"../../literate/examples/elastic_impact.jl\"","category":"page"},{"location":"examples/elastic_impact/#Transfer-schemes","page":"Transfer schemes","title":"Transfer schemes","text":"","category":"section"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"In this example, the following transfer schemes are demonstrated:","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"PIC–FLIP mixed transfer[1]\nAffine PIC (APIC) transfer[2]\nTaylor PIC (TPIC) transfer[3]","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"The problem evolves the elastic impact between two rings, which is consistent with previous studies[4].","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"[1]: Zhu, Y. and Bridson, R., 2005. Animating sand as a fluid. ACM Transactions on Graphics (TOG), 24(3), pp.965-972.","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"[2]: Jiang, C., Schroeder, C., Selle, A., Teran, J. and Stomakhin, A., 2015. The affine particle-in-cell method. ACM Transactions on Graphics (TOG), 34(4), pp.1-10.","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"[3]: Nakamura, K., Matsumura, S. and Mizutani, T., 2023. Taylor particle-in-cell transfer and kernel correction for material point method. Computer Methods in Applied Mechanics and Engineering, 403, p.115720.","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"[4]: Li, X., Fang, Y., Li, M. and Jiang, C., 2022. BFEMP: Interpenetration-free MPM–FEM coupling with barrier contact. Computer Methods in Applied Mechanics and Engineering, 390, p.114350.","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"using Tesserae\n\nabstract type Transfer end\nstruct FLIP <: Transfer α::Float64 end\nstruct APIC <: Transfer end\nstruct TPIC <: Transfer end\n\nfunction main(transfer::Transfer = FLIP(1.0))\n\n # Simulation parameters\n h = 0.1 # Grid spacing\n T = 0.6 # Time span\n CFL = 0.8 # Courant number\n\n # Material constants\n E = 100e6 # Young's modulus\n ν = 0.2 # Poisson's ratio\n λ = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter\n μ = E / 2(1 + ν) # Shear modulus\n ρ⁰ = 1e3 # Initial density\n\n # Geometry\n L = 20.0 # Length of domain\n W = 15.0 # Width of domain\n rᵢ = 3.0 # Inner radius of rings\n rₒ = 4.0 # Outer radius of rings\n\n GridProp = @NamedTuple begin\n x :: Vec{2, Float64}\n m :: Float64\n m⁻¹ :: Float64\n mv :: Vec{2, Float64}\n f :: Vec{2, Float64}\n v :: Vec{2, Float64}\n vⁿ :: Vec{2, Float64}\n end\n ParticleProp = @NamedTuple begin\n x :: Vec{2, Float64}\n m :: Float64\n V⁰ :: Float64\n V :: Float64\n v :: Vec{2, Float64}\n ∇v :: SecondOrderTensor{2, Float64, 4}\n σ :: SymmetricSecondOrderTensor{2, Float64, 3}\n F :: SecondOrderTensor{2, Float64, 4}\n B :: SecondOrderTensor{2, Float64, 4} # for APIC\n end\n\n # Background grid\n grid = generate_grid(GridProp, CartesianMesh(h, (-L,L), (-W/2,W/2)))\n\n # Particles\n particles = let\n pts = generate_particles(ParticleProp, grid.x)\n pts.V .= pts.V⁰ .= volume(grid.x) / length(pts)\n\n lhs = findall(pts.x) do (x, y)\n rᵢ^2 < (x+L/4)^2+y^2 < rₒ^2\n end\n rhs = findall(pts.x) do (x, y)\n rᵢ^2 < (x-L/4)^2+y^2 < rₒ^2\n end\n\n # Set initial velocities\n @. pts.v[lhs] = Vec(30, 0)\n @. pts.v[rhs] = -Vec(30, 0)\n\n pts[[lhs; rhs]]\n end\n @. particles.m = ρ⁰ * particles.V⁰\n @. particles.F = one(particles.F)\n @show length(particles)\n\n # Interpolation\n mpvalues = generate_mpvalues(QuadraticBSpline(), grid.x, length(particles))\n\n # Material model (neo-Hookean)\n function caucy_stress(F)\n b = F ⋅ F'\n J = det(F)\n (μ*(b-I) + λ*log(J)*I) / J\n end\n\n # Outputs\n outdir = mkpath(joinpath(\"output\", \"elastic_impact\"))\n pvdfile = joinpath(outdir, \"paraview\")\n closepvd(openpvd(pvdfile)) # create file\n\n t = 0.0\n step = 0\n fps = 120\n savepoints = collect(LinRange(t, T, round(Int, T*fps)+1))\n\n Tesserae.@showprogress while t < T\n\n # Calculate timestep based on the wave speed\n vmax = maximum(@. sqrt((λ+2μ) / (particles.m/particles.V)) + norm(particles.v))\n Δt = CFL * spacing(grid) / vmax\n\n # Update interpolation values\n for p in eachindex(particles, mpvalues)\n update!(mpvalues[p], particles.x[p], grid.x)\n end\n\n # Particle-to-grid transfer\n if transfer isa FLIP\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * v[p]\n f[i] = @∑ -V[p] * σ[p] ⋅ ∇w[ip]\n end\n elseif transfer isa APIC\n local Dₚ⁻¹ = inv(1/4 * h^2 * I)\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * (v[p] + B[p] ⋅ Dₚ⁻¹ ⋅ (x[i] - x[p]))\n f[i] = @∑ -V[p] * σ[p] ⋅ ∇w[ip]\n end\n elseif transfer isa TPIC\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * (v[p] + ∇v[p] ⋅ (x[i] - x[p]))\n f[i] = @∑ -V[p] * σ[p] ⋅ ∇w[ip]\n end\n end\n\n # Update grid velocity\n @. grid.m⁻¹ = inv(grid.m) * !iszero(grid.m)\n @. grid.vⁿ = grid.mv * grid.m⁻¹\n @. grid.v = grid.vⁿ + Δt * grid.f * grid.m⁻¹\n\n # Grid-to-particle transfer\n if transfer isa FLIP\n local α = transfer.α\n @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] = @∑ w[ip] * ((1-α)*v[i] + α*(v[p] + (v[i]-vⁿ[i])))\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n x[p] += @∑ Δt * (w[ip] * v[i])\n\n end\n elseif transfer isa APIC\n @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] = @∑ w[ip] * v[i]\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n B[p] = @∑ w[ip] * v[i] ⊗ (x[i]-x[p])\n x[p] += Δt * v[p]\n end\n elseif transfer isa TPIC\n @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] = @∑ w[ip] * v[i]\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n x[p] += Δt * v[p]\n end\n end\n\n # Update other particle properties\n for p in eachindex(particles)\n ∇uₚ = Δt * particles.∇v[p]\n Fₚ = (I + ∇uₚ) ⋅ particles.F[p]\n σₚ = caucy_stress(Fₚ)\n particles.σ[p] = σₚ\n particles.F[p] = Fₚ\n particles.V[p] = det(Fₚ) * particles.V⁰[p]\n end\n\n t += Δt\n step += 1\n\n if t > first(savepoints)\n popfirst!(savepoints)\n openpvd(pvdfile; append=true) do pvd\n openvtm(string(pvdfile, step)) do vtm\n function stress3x3(F)\n z = zero(Mat{2,1})\n F3x3 = [F z\n z' 1]\n caucy_stress(F3x3)\n end\n openvtk(vtm, particles.x) do vtk\n vtk[\"velocity\"] = particles.v\n vtk[\"von Mises\"] = @. vonmises(stress3x3(particles.F))\n end\n openvtk(vtm, grid.x) do vtk\n vtk[\"velocity\"] = grid.v\n end\n pvd[t] = vtm\n end\n end\n end\n end\nend","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"This page was generated using Literate.jl.","category":"page"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"EditURL = \"../../literate/examples/dam_break.jl\"","category":"page"},{"location":"examples/dam_break/#Stabilized-mixed-MPM-for-incompressible-fluid","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"","category":"section"},{"location":"examples/dam_break/#Main-function","page":"Stabilized mixed MPM for incompressible fluid","title":"Main function","text":"","category":"section"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"using Tesserae\nusing LinearAlgebra\n\nabstract type Transfer end\nstruct FLIP <: Transfer α::Float64 end\nstruct TPIC <: Transfer end\n\nfunction main(transfer::Transfer = FLIP(1.0))\n\n # Simulation parameters\n h = 0.02 # Grid spacing\n T = 7.0 # Time span\n g = 9.81 # Gravity acceleration\n Δt = 2.5e-3 # Timestep\n\n # Material constants\n ρ = 1.0e3 # Initial density\n μ = 1.01e-3 # Dynamic viscosity (Pa⋅s)\n\n # Newmark-beta method\n β = 0.5\n γ = 1.0\n\n # Utils\n cellnodes(cell) = cell:(cell+oneunit(cell))\n cellcenter(cell, mesh) = mean(mesh[cellnodes(cell)])\n\n # Properties for grid and particles\n GridProp = @NamedTuple begin\n X :: Vec{2, Float64}\n x :: Vec{2, Float64}\n m :: Float64\n m⁻¹ :: Float64\n v :: Vec{2, Float64}\n vⁿ :: Vec{2, Float64}\n mv :: Vec{2, Float64}\n a :: Vec{2, Float64}\n aⁿ :: Vec{2, Float64}\n ma :: Vec{2, Float64}\n u :: Vec{2, Float64}\n p :: Float64\n # δ-correction\n V :: Float64\n Ṽ :: Float64\n E :: Float64\n # Residuals\n u_p :: Vec{3, Float64}\n R_mom :: Vec{2, Float64}\n R_mas :: Float64\n end\n ParticleProp = @NamedTuple begin\n x :: Vec{2, Float64}\n m :: Float64\n V :: Float64\n v :: Vec{2, Float64}\n ∇v :: SecondOrderTensor{2, Float64, 4}\n a :: Vec{2, Float64}\n ∇a :: SecondOrderTensor{2, Float64, 4}\n p :: Float64\n ∇p :: Vec{2, Float64}\n s :: SymmetricSecondOrderTensor{2, Float64, 3}\n b :: Vec{2, Float64}\n # δ-correction\n ∇E² :: Vec{2, Float64}\n # Stabilization\n τ₁ :: Float64\n τ₂ :: Float64\n end\n\n # Background grid\n grid = generate_grid(GridProp, CartesianMesh(h, (0,3.22), (0,2.5)))\n for cell in CartesianIndices(size(grid).-1)\n for i in cellnodes(cell)\n grid.V[i] += (spacing(grid)/2)^2\n end\n end\n\n # Particles\n particles = generate_particles(ParticleProp, grid.X; spacing=1/3)\n particles.V .= volume(grid.X) / length(particles)\n filter!(pt -> pt.x[1]<1.2 && pt.x[2]<0.6, particles)\n @. particles.m = ρ * particles.V\n @. particles.b = Vec(0,-g)\n @show length(particles)\n\n # Interpolation\n it = KernelCorrection(QuadraticBSpline())\n mpvalues = map(p -> MPValue(it, grid.X), eachindex(particles))\n mpvalues_cell = map(CartesianIndices(size(grid).-1)) do cell\n mp = MPValue(it, grid.X)\n xc = cellcenter(cell, grid.X)\n update!(mp, xc, grid.X)\n end\n\n # BlockSpace for threaded computation\n blockspace = BlockSpace(grid.X)\n\n # Sparse matrix\n A = create_sparse_matrix(Vec{3}, it, grid.X)\n\n # Output\n outdir = mkpath(joinpath(\"output\", \"dam_break\"))\n pvdfile = joinpath(outdir, \"paraview\")\n closepvd(openpvd(pvdfile)) # Create file\n\n t = 0.0\n step = 0\n fps = 30\n savepoints = collect(LinRange(t, T, round(Int, T*fps)+1))\n\n Tesserae.@showprogress while t < T\n\n # Update interpolation values based on the nodes of active cells\n # where the particles are located\n activenodes = falses(size(grid))\n for p in eachindex(particles)\n cell = whichcell(particles.x[p], grid.X)\n activenodes[cellnodes(cell)] .= true\n end\n for p in eachindex(mpvalues, particles)\n update!(mpvalues[p], particles.x[p], grid.X, activenodes)\n end\n for cell in CartesianIndices(size(grid) .- 1)\n xc = cellcenter(cell, grid.X)\n update!(mpvalues_cell[cell], xc, grid.X, activenodes)\n end\n\n update!(blockspace, particles.x)\n\n if transfer isa FLIP\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * v[p]\n ma[i] = @∑ w[ip] * m[p] * a[p]\n end\n elseif transfer isa TPIC\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * (v[p] + ∇v[p] ⋅ (X[i] - x[p]))\n ma[i] = @∑ w[ip] * m[p] * (a[p] + ∇a[p] ⋅ (X[i] - x[p]))\n end\n end\n\n # Remove negative mass occasionally generated due to the kernel correction\n # by forcing negative values to zero. While not an ideal solution, this sufficiently\n # reduces spurious pressure oscillations without sacrificing simulation accuracy.\n @. grid.m = max(grid.m, 0)\n\n @. grid.m⁻¹ = inv(grid.m) * !iszero(grid.m)\n @. grid.vⁿ = grid.mv * grid.m⁻¹\n @. grid.aⁿ = grid.ma * grid.m⁻¹\n\n # Update a dof map\n dofmask = trues(3, size(grid)...)\n for i in eachindex(grid)\n dofmask[:,i] .= !iszero(grid.m[i])\n end\n for i in @view eachindex(grid)[[begin,end],:] # Walls\n grid.vⁿ[i] = grid.vⁿ[i] .* (false,true)\n grid.aⁿ[i] = grid.aⁿ[i] .* (false,true)\n dofmask[1,i] = false\n end\n for i in @view eachindex(grid)[:,begin] # Floor\n grid.vⁿ[i] = grid.vⁿ[i] .* (true,false)\n grid.aⁿ[i] = grid.aⁿ[i] .* (true,false)\n dofmask[2,i] = false\n end\n dofmap = DofMap(dofmask)\n\n # Solve grid position, dispacement, velocity, acceleration and pressure by VMS method\n state = (; grid, particles, mpvalues, mpvalues_cell, blockspace, ρ, μ, β, γ, A, dofmap, Δt)\n Δt′ = variational_multiscale_method(state)\n\n if transfer isa FLIP\n local α = transfer.α\n @threaded @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] = @∑ ((1-α)*v[i] + α*(v[p] + Δt*((1-γ)*a[p] + γ*a[i]))) * w[ip]\n a[p] = @∑ a[i] * w[ip]\n x[p] = @∑ x[i] * w[ip]\n end\n elseif transfer isa TPIC\n @threaded @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] = @∑ v[i] * w[ip]\n a[p] = @∑ a[i] * w[ip]\n x[p] = @∑ x[i] * w[ip]\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n ∇a[p] = @∑ a[i] ⊗ ∇w[ip]\n end\n end\n\n # Particle shifting based on the δ-correction\n particle_shifting(state)\n\n # Remove particles that accidentally move outside of the mesh\n outside = findall(x->!isinside(x, grid.X), particles.x)\n deleteat!(particles, outside)\n deleteat!(mpvalues, outside)\n\n t += Δt′\n step += 1\n\n # Write results\n if t > first(savepoints)\n popfirst!(savepoints)\n openpvd(pvdfile; append=true) do pvd\n openvtm(string(pvdfile, step)) do vtm\n vorticity(∇v) = ∇v[2,1] - ∇v[1,2]\n openvtk(vtm, particles.x) do vtk\n vtk[\"pressure\"] = particles.p\n vtk[\"velocity\"] = particles.v\n vtk[\"vorticity\"] = vorticity.(particles.∇v)\n end\n openvtk(vtm, grid.X) do vtk\n end\n pvd[t] = vtm\n end\n end\n end\n end\nend","category":"page"},{"location":"examples/dam_break/#Variational-multiscale-method","page":"Stabilized mixed MPM for incompressible fluid","title":"Variational multiscale method","text":"","category":"section"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"function variational_multiscale_method(state)\n\n # The simulation might fail occasionally due to regions with very small masses,\n # such as splashes[^1]. Therefore, in this script, if Newton's method doesn't converge,\n # a half time step is applied.\n\n (; grid, dofmap, Δt) = state\n @. grid.u_p = zero(grid.u_p)\n\n solved = false\n while !solved\n # Reconstruct state using the current time step\n state = merge(state, (; Δt))\n\n # Compute VMS stabilization coefficients using current grid velocity,\n # which is used for Jacobian matrix\n grid.v .= grid.vⁿ\n compute_VMS_stabilization_coefficients(state)\n\n # Try computing the Jacobian matrix and performing its LU decomposition.\n # In this formulation[^1], the initial Jacobian matrix is used in all Newton's iterations.\n # If the computation fails, use a smaller time step.\n J = lu(jacobian(state); check=false)\n issuccess(J) || (Δt /= 2; continue)\n\n # Solve nonlinear system\n U = zeros(ndofs(dofmap)) # Initialize nodal dispacement and pressure with zero\n solved = Tesserae.newton!(U, U->residual(U,state), U->J;\n linsolve=(x,A,b)->ldiv!(x,A,b), atol=1e-10, rtol=1e-10)\n\n # If the simulation fails to solve, retry with a smaller time step\n solved || (Δt /= 2)\n end\n\n # Update the positions of grid nodes\n @. grid.x = grid.X + grid.u\n\n # Return the acutally applied time step\n Δt\nend","category":"page"},{"location":"examples/dam_break/#VMS-stabilization-coefficients","page":"Stabilized mixed MPM for incompressible fluid","title":"VMS stabilization coefficients","text":"","category":"section"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"function compute_VMS_stabilization_coefficients(state)\n (; grid, particles, mpvalues_cell, ρ, μ, Δt) = state\n\n c₁ = 4.0\n c₂ = 2.0\n τdyn = 1.0\n h = sqrt(4*spacing(grid)^2/π)\n\n # In the following computation, `@G2P` is unavailable\n # due to the use of `mpvalues_cell`\n for p in eachindex(particles)\n v̄ₚ = zero(eltype(particles.v))\n mp = mpvalues_cell[whichcell(particles.x[p], grid.X)]\n gridindices = neighboringnodes(mp, grid)\n for ip in eachindex(gridindices)\n i = gridindices[ip]\n v̄ₚ += mp.w[ip] * grid.v[i]\n end\n τ₁ = inv(ρ*τdyn/Δt + c₂*ρ*norm(v̄ₚ)/h + c₁*μ/h^2)\n τ₂ = h^2 / (c₁*τ₁)\n particles.τ₁[p] = τ₁\n particles.τ₂[p] = τ₂\n end\nend","category":"page"},{"location":"examples/dam_break/#δ-correction","page":"Stabilized mixed MPM for incompressible fluid","title":"δ-correction","text":"","category":"section"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"function particle_shifting(state)\n (; grid, particles, mpvalues) = state\n\n @P2G grid=>i particles=>p mpvalues=>ip begin\n Ṽ[i] = @∑ V[p] * w[ip]\n E[i] = max(0, -V[i] + Ṽ[i])\n end\n\n E² = sum(E->E^2, grid.E)\n @G2P grid=>i particles=>p mpvalues=>ip begin\n ∇E²[p] = @∑ 2V[p] * E[i] * ∇w[ip]\n end\n\n b₀ = E² / sum(∇E²->∇E²⋅∇E², particles.∇E²)\n\n for p in eachindex(particles)\n xₚ = particles.x[p] - b₀ * particles.∇E²[p]\n if isinside(xₚ, grid.X)\n particles.x[p] = xₚ\n end\n end\nend","category":"page"},{"location":"examples/dam_break/#Residual-vector","page":"Stabilized mixed MPM for incompressible fluid","title":"Residual vector","text":"","category":"section"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"function residual(U, state)\n (; grid, particles, mpvalues, μ, β, γ, dofmap, Δt) = state\n\n # Map `U` to grid dispacement and pressure\n dofmap(grid.u_p) .= U\n grid.u .= map(x->@Tensor(x[1:2]), grid.u_p)\n grid.p .= map(x->x[3], grid.u_p)\n\n # Recompute nodal velocity and acceleration based on the Newmark-beta method\n @. grid.v = γ/(β*Δt)*grid.u - (γ/β-1)*grid.vⁿ - Δt/2*(γ/β-2)*grid.aⁿ\n @. grid.a = 1/(β*Δt^2)*grid.u - 1/(β*Δt)*grid.vⁿ - (1/2β-1)*grid.aⁿ\n\n # Recompute particle properties for residual vector\n @G2P grid=>i particles=>p mpvalues=>ip begin\n a[p] = @∑ a[i] * w[ip]\n p[p] = @∑ p[i] * w[ip]\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n ∇p[p] = @∑ p[i] * ∇w[ip]\n s[p] = 2μ * symmetric(∇v[p])\n end\n\n # Compute VMS stabilization coefficients based on the current nodal velocity\n compute_VMS_stabilization_coefficients(state)\n\n # Compute residual values\n @P2G grid=>i particles=>p mpvalues=>ip begin\n R_mom[i] = @∑ V[p]*s[p]⋅∇w[ip] - m[p]*b[p]*w[ip] - V[p]*p[p]*∇w[ip] + τ₂[p]*V[p]*tr(∇v[p])*∇w[ip]\n R_mas[i] = @∑ V[p]*tr(∇v[p])*w[ip] + τ₁[p]*m[p]*(a[p]-b[p])⋅∇w[ip] + τ₁[p]*V[p]*∇p[p]⋅∇w[ip]\n R_mom[i] += m[i]*a[i]\n end\n\n # Map grid values to vector `R`\n dofmap(map(vcat, grid.R_mom, grid.R_mas))\nend","category":"page"},{"location":"examples/dam_break/#Jacobian-matrix","page":"Stabilized mixed MPM for incompressible fluid","title":"Jacobian matrix","text":"","category":"section"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"function jacobian(state)\n (; grid, particles, mpvalues, blockspace, ρ, μ, β, γ, A, dofmap, Δt) = state\n\n # Construct the Jacobian matrix\n cₚ = 2μ * one(SymmetricFourthOrderTensor{2})\n I(i,j) = ifelse(i===j, one(Mat{2,2}), zero(Mat{2,2}))\n @threaded @P2G_Matrix grid=>(i,j) particles=>p mpvalues=>(ip,jp) blockspace begin\n A[i,j] = @∑ begin\n Kᵤᵤ = (γ/(β*Δt) * ∇w[ip] ⋅ cₚ ⋅ ∇w[jp]) * V[p] + 1/(β*Δt^2) * I(i,j) * m[p] * w[jp]\n Kᵤₚ = -∇w[ip] * w[jp] * V[p]\n Kₚᵤ = (γ/(β*Δt)) * w[ip] * ∇w[jp] * V[p]\n K̂ᵤᵤ = γ/(β*Δt) * τ₂[p] * ∇w[ip] ⊗ ∇w[jp] * V[p]\n K̂ₚᵤ = 1/(β*Δt^2) * τ₁[p] * ρ * ∇w[ip] * w[jp] * V[p]\n K̂ₚₚ = τ₁[p] * ∇w[ip] ⋅ ∇w[jp] * V[p]\n [Kᵤᵤ+K̂ᵤᵤ Kᵤₚ\n Mat{1,2}(Kₚᵤ+K̂ₚᵤ) K̂ₚₚ]\n end\n end\n\n # Extract the activated degrees of freedom\n submatrix(A, dofmap)\nend","category":"page"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"","category":"page"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"This page was generated using Literate.jl.","category":"page"},{"location":"#Tesserae","page":"Home","title":"Tesserae","text":"","category":"section"},{"location":"examples/implicit_jacobian_free/","page":"Jacobian-free Newton–Krylov method","title":"Jacobian-free Newton–Krylov method","text":"EditURL = \"../../literate/examples/implicit_jacobian_free.jl\"","category":"page"},{"location":"examples/implicit_jacobian_free/#Jacobian-free-Newton–Krylov-method","page":"Jacobian-free Newton–Krylov method","title":"Jacobian-free Newton–Krylov method","text":"","category":"section"},{"location":"examples/implicit_jacobian_free/","page":"Jacobian-free Newton–Krylov method","title":"Jacobian-free Newton–Krylov method","text":"using Tesserae\n\nusing IterativeSolvers: gmres!\nusing LinearMaps: LinearMap\n\nfunction main()\n\n # Simulation parameters\n h = 0.05 # Grid spacing\n T = 1.0 # Time span\n g = 20.0 # Gravity acceleration\n Δt = 0.02 # Timestep\n\n # Material constants\n E = 1e6 # Young's modulus\n ν = 0.3 # Poisson's ratio\n λ = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter\n μ = E / 2(1 + ν) # Shear modulus\n ρ⁰ = 500.0 # Density\n\n # Newmark-beta integration\n β = 1/4\n γ = 1/2\n\n GridProp = @NamedTuple begin\n X :: Vec{3, Float64}\n m :: Float64\n m⁻¹ :: Float64\n v :: Vec{3, Float64}\n vⁿ :: Vec{3, Float64}\n mv :: Vec{3, Float64}\n a :: Vec{3, Float64}\n aⁿ :: Vec{3, Float64}\n ma :: Vec{3, Float64}\n u :: Vec{3, Float64}\n f :: Vec{3, Float64}\n end\n ParticleProp = @NamedTuple begin\n x :: Vec{3, Float64}\n m :: Float64\n V⁰ :: Float64\n v :: Vec{3, Float64}\n a :: Vec{3, Float64}\n b :: Vec{3, Float64}\n ∇u :: SecondOrderTensor{3, Float64, 9}\n F :: SecondOrderTensor{3, Float64, 9}\n ΔF⁻¹ :: SecondOrderTensor{3, Float64, 9}\n τ :: SecondOrderTensor{3, Float64, 9}\n c :: FourthOrderTensor{3, Float64, 81}\n end\n\n # Background grid\n grid = generate_grid(GridProp, CartesianMesh(h, (0.0,1.2), (0.0,2.0), (-0.2,0.2)))\n\n # Particles\n beam = Tesserae.Box((0,1), (0.85,1.15), (-0.15,0.15))\n particles = generate_particles(ParticleProp, grid.X; domain=beam, alg=GridSampling())\n particles.V⁰ .= volume(beam) / length(particles)\n @. particles.m = ρ⁰ * particles.V⁰\n @. particles.F = one(particles.F)\n @. particles.b = Vec(0,-g,0)\n @show length(particles)\n\n # Interpolation\n # Use the kernel correction to properly handle the boundary conditions\n mpvalues = generate_mpvalues(KernelCorrection(QuadraticBSpline()), grid.X, length(particles))\n\n # Neo-Hookean model\n function kirchhoff_stress(F)\n J = det(F)\n b = symmetric(F ⋅ F')\n μ*(b-I) + λ*log(J)*I\n end\n\n # Outputs\n outdir = mkpath(joinpath(\"output\", \"implicit_jacobian_free\"))\n pvdfile = joinpath(outdir, \"paraview\")\n closepvd(openpvd(pvdfile)) # Create file\n\n t = 0.0\n step = 0\n\n Tesserae.@showprogress while t < T\n\n for p in eachindex(particles, mpvalues)\n update!(mpvalues[p], particles.x[p], grid.X)\n end\n\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * v[p]\n ma[i] = @∑ w[ip] * m[p] * a[p]\n end\n\n # Compute the grid velocity and acceleration at t = tⁿ\n @. grid.m⁻¹ = inv(grid.m) * !iszero(grid.m)\n @. grid.vⁿ = grid.mv * grid.m⁻¹\n @. grid.aⁿ = grid.ma * grid.m⁻¹\n\n # Update a dof map\n dofmask = trues(3, size(grid)...)\n for i in eachindex(grid)\n dofmask[:,i] .= !iszero(grid.m[i])\n end\n for i in eachindex(grid)[1,:,:]\n dofmask[:,i] .= false\n grid.vⁿ[i] = zero(Vec{3})\n grid.aⁿ[i] = zero(Vec{3})\n end\n dofmap = DofMap(dofmask)\n\n # Solve the nonlinear equation\n state = (; grid, particles, mpvalues, kirchhoff_stress, β, γ, dofmap, Δt)\n @. grid.u = zero(grid.u) # Set zero dispacement for the first guess of the solution\n U = copy(dofmap(grid.u)) # Convert grid data to plain vector data\n compute_residual(U) = residual(U, state)\n compute_jacobian(U) = jacobian(U, state)\n Tesserae.newton!(U, compute_residual, compute_jacobian; linsolve = (x,A,b)->gmres!(x,A,b))\n\n # Grid dispacement, velocity and acceleration have been updated during Newton's iterations\n @G2P grid=>i particles=>p mpvalues=>ip begin\n ∇u[p] = @∑ u[i] ⊗ ∇w[ip]\n a[p] = @∑ w[ip] * a[i]\n v[p] += @∑ Δt * w[ip] * ((1-γ)*a[p] + γ*a[i])\n x[p] = @∑ w[ip] * (X[i] + u[i])\n F[p] = (I + ∇u[p]) ⋅ F[p]\n end\n\n t += Δt\n step += 1\n\n openpvd(pvdfile; append=true) do pvd\n openvtk(string(pvdfile, step), particles.x) do vtk\n vtk[\"velocity\"] = particles.v\n vtk[\"von Mises\"] = @. vonmises(particles.τ / det(particles.F))\n pvd[t] = vtk\n end\n end\n end\nend\n\nfunction residual(U::AbstractVector, state)\n (; grid, particles, mpvalues, kirchhoff_stress, β, γ, dofmap, Δt) = state\n\n dofmap(grid.u) .= U\n @. grid.a = (1/(β*Δt^2))*grid.u - (1/(β*Δt))*grid.vⁿ - (1/2β-1)*grid.aⁿ\n @. grid.v = grid.vⁿ + Δt*((1-γ)*grid.aⁿ + γ*grid.a)\n\n transposing_tensor(σ) = @einsum (i,j,k,l) -> σ[i,l] * one(σ)[j,k]\n @G2P grid=>i particles=>p mpvalues=>ip begin\n # In addition to updating the stress tensor, the stiffness tensor,\n # which is utilized in the Jacobian-vector product, is also updated.\n ∇u[p] = @∑ u[i] ⊗ ∇w[ip]\n ΔF⁻¹[p] = inv(I + ∇u[p])\n c[p], τ[p] = gradient(∇u -> kirchhoff_stress((I + ∇u) ⋅ F[p]), ∇u[p], :all)\n c[p] = c[p] - transposing_tensor(τ[p] ⋅ ΔF⁻¹[p]')\n end\n @P2G grid=>i particles=>p mpvalues=>ip begin\n f[i] = @∑ -V⁰[p] * τ[p] ⋅ (∇w[ip] ⋅ ΔF⁻¹[p]) + w[ip] * m[p] * b[p]\n end\n\n @. β*Δt^2 * ($dofmap(grid.a) - $dofmap(grid.f) * $dofmap(grid.m⁻¹))\nend\n\nfunction jacobian(U::AbstractVector, state)\n (; grid, particles, mpvalues, β, dofmap, Δt) = state\n\n # Create a linear map to represent Jacobian-vector product J*δU.\n # `U` is acutally not used because the stiffness tensor is already calculated\n # when computing the residual vector.\n LinearMap(ndofs(dofmap)) do JδU, δU\n dofmap(grid.u) .= δU\n\n @G2P grid=>i particles=>p mpvalues=>ip begin\n ∇u[p] = @∑ u[i] ⊗ ∇w[ip]\n τ[p] = c[p] ⊡ ∇u[p]\n end\n @P2G grid=>i particles=>p mpvalues=>ip begin\n f[i] = @∑ -V⁰[p] * τ[p] ⋅ (∇w[ip] ⋅ ΔF⁻¹[p])\n end\n\n @. JδU = δU - β*Δt^2 * $dofmap(grid.f) * $dofmap(grid.m⁻¹)\n end\nend","category":"page"},{"location":"examples/implicit_jacobian_free/","page":"Jacobian-free Newton–Krylov method","title":"Jacobian-free Newton–Krylov method","text":"","category":"page"},{"location":"examples/implicit_jacobian_free/","page":"Jacobian-free Newton–Krylov method","title":"Jacobian-free Newton–Krylov method","text":"This page was generated using Literate.jl.","category":"page"},{"location":"examples/rigid_body_contact/","page":"Frictional contact with rigid body","title":"Frictional contact with rigid body","text":"EditURL = \"../../literate/examples/rigid_body_contact.jl\"","category":"page"},{"location":"examples/rigid_body_contact/#Frictional-contact-with-rigid-body","page":"Frictional contact with rigid body","title":"Frictional contact with rigid body","text":"","category":"section"},{"location":"examples/rigid_body_contact/","page":"Frictional contact with rigid body","title":"Frictional contact with rigid body","text":"using Tesserae\n\nmutable struct Disk{dim, T}\n x::Vec{dim, T}\n v::Vec{dim, T}\nend\n\nfunction main()\n\n # Simulation parameters\n h = 0.004 # Grid spacing\n T = 5.0 # Time span\n g = 9.81 # Gravity acceleration\n CFL = 1.0 # Courant number\n\n # Material constants\n E = 1000e3 # Young's modulus\n ν = 0.49 # Poisson's ratio\n λ = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter\n G = E / 2(1 + ν) # Shear modulus\n σy = 1.0e3 # Yield stress\n ρ⁰ = 1.0e3 # Initial density\n\n # Disk\n D = 0.04\n disk = Disk(Vec(0, 7.5D), Vec(0, -0.25D))\n\n # Contact parameters\n k = 1e6 # Penalty coefficient\n μ = 0.6 # Friction coefficient\n\n # Utils\n @inline function contact_force_normal(x, r, x_disk)\n d = x - x_disk\n k * max(D/2 - (norm(d)-r), 0) * normalize(d)\n end\n @inline function contact_force_tangential(fₙ, v, m, Δt)\n iszero(fₙ) && return zero(fₙ)\n n = normalize(fₙ)\n fₜ = -m * (v-(v⋅n)*n) / Δt # Sticking force\n min(1, μ*norm(fₙ)/norm(fₜ)) * fₜ\n end\n\n # Properties for grid and particles\n GridProp = @NamedTuple begin\n x :: Vec{2, Float64}\n m :: Float64\n m⁻¹ :: Float64\n mv :: Vec{2, Float64}\n fint :: Vec{2, Float64}\n fext :: Vec{2, Float64}\n v :: Vec{2, Float64}\n vⁿ :: Vec{2, Float64}\n end\n ParticleProp = @NamedTuple begin\n x :: Vec{2, Float64}\n m :: Float64\n V :: Float64\n r :: Float64\n v :: Vec{2, Float64}\n ∇v :: SecondOrderTensor{2, Float64, 4}\n F :: SecondOrderTensor{3, Float64, 9}\n σ :: SymmetricSecondOrderTensor{3, Float64, 6}\n ϵ :: SymmetricSecondOrderTensor{3, Float64, 6}\n b :: Vec{2, Float64}\n end\n\n # Background grid\n H = 7D # Ground height\n grid = generate_grid(GridProp, CartesianMesh(h, (0,5D), (0,H+D)))\n\n # Particles\n particles = generate_particles(ParticleProp, grid.x)\n disk_points = filter(particles.x) do (x,y) # Points representing a disk just for visualization\n x^2 + (y-7.5D)^2 < (D/2)^2\n end\n particles.V .= volume(grid.x) / length(particles)\n filter!(pt -> pt.x[2] < H, particles)\n @. particles.m = ρ⁰ * particles.V\n @. particles.r = particles.V^(1/2) / 2\n @. particles.b = Vec(0,-g)\n @. particles.F = one(particles.F)\n for p in eachindex(particles)\n x, y = particles.x[p]\n σ_y = -ρ⁰ * g * (H - y)\n σ_x = ν/(1-ν) * σ_y\n particles.σ[p] = [σ_x 0.0 0.0\n 0.0 σ_y 0.0\n 0.0 0.0 σ_x]\n end\n @show length(particles)\n\n # Interpolation\n mpvalues = generate_mpvalues(KernelCorrection(QuadraticBSpline()), grid.x, length(particles))\n\n # Output\n outdir = mkpath(joinpath(\"output\", \"rigid_body_contact\"))\n pvdfile = joinpath(outdir, \"paraview\")\n closepvd(openpvd(pvdfile)) # Create file\n\n t = 0.0\n step = 0\n fps = 20\n savepoints = collect(LinRange(t, T, round(Int, T*fps)+1))\n\n Tesserae.@showprogress while t < T\n\n vmax = maximum(@. sqrt((λ+2G) / (particles.m/particles.V)) + norm(particles.v))\n Δt = CFL * spacing(grid) / vmax\n\n for p in eachindex(particles, mpvalues)\n update!(mpvalues[p], particles.x[p], grid.x)\n end\n\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * (v[p] + ∇v[p] ⋅ (x[i] - x[p])) # Taylor transfer\n fint[i] = @∑ -V[p] * resizedim(σ[p], Val(2)) ⋅ ∇w[ip] + w[ip] * m[p] * b[p]\n fext[i] = @∑ w[ip] * contact_force_normal(x[p], r[p], disk.x)\n end\n\n @. grid.m⁻¹ = inv(grid.m) * !iszero(grid.m)\n @. grid.vⁿ = grid.mv * grid.m⁻¹\n @. grid.v = grid.vⁿ + Δt * grid.fint * grid.m⁻¹\n @. grid.fext += contact_force_tangential(grid.fext, grid.v-disk.v, grid.m, Δt)\n @. grid.v += Δt * grid.fext * grid.m⁻¹\n\n for i in eachindex(grid)[[begin,end],:]\n grid.v[i] = grid.v[i] .* (false,true)\n end\n for i in eachindex(grid)[:,[begin,end]]\n grid.v[i] = grid.v[i] .* (false,false)\n end\n\n @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] = @∑ w[ip] * v[i] # PIC transfer\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n x[p] += Δt * v[p]\n end\n\n for p in 1:length(particles)\n Δϵ = resizedim(Δt * symmetric(particles.∇v[p]), Val(3))\n particles.σ[p] = vonmises_model(particles.σ[p], Δϵ; λ, G, σy)\n particles.V[p] *= 1 + tr(Δϵ)\n particles.ϵ[p] += Δϵ\n end\n\n disk.x += disk.v * Δt\n @. disk_points += disk.v * Δt\n\n t += Δt\n step += 1\n\n if t > first(savepoints)\n popfirst!(savepoints)\n openpvd(pvdfile; append=true) do pvd\n openvtm(string(pvdfile, step)) do vtm\n deviatoric_strain(ϵ) = sqrt(2/3 * dev(ϵ) ⊡ dev(ϵ))\n openvtk(vtm, particles.x) do vtk\n vtk[\"vonmises\"] = @. vonmises(particles.σ)\n vtk[\"deviatoric strain\"] = @. deviatoric_strain(particles.ϵ)\n end\n openvtk(vtm, disk_points) do vtk\n end\n pvd[t] = vtm\n end\n end\n end\n end\nend\n\nfunction vonmises_model(σⁿ, Δϵ; λ, G, σy)\n δ = one(SymmetricSecondOrderTensor{3})\n I = one(SymmetricFourthOrderTensor{3})\n cᵉ = λ*δ⊗δ + 2G*I\n σ_trial = σⁿ + cᵉ ⊡ Δϵ\n dfdσ, f_trial = gradient(σ -> vonmises(σ) - σy, σ_trial, :all)\n if f_trial > 0\n dλ = f_trial / (dfdσ ⊡ cᵉ ⊡ dfdσ)\n σ = σ_trial - cᵉ ⊡ (dλ * dfdσ)\n else\n σ = σ_trial\n end\n if mean(σ) > 0 # simple tension cut-off\n σ = dev(σ)\n end\n σ\nend","category":"page"},{"location":"examples/rigid_body_contact/","page":"Frictional contact with rigid body","title":"Frictional contact with rigid body","text":"","category":"page"},{"location":"examples/rigid_body_contact/","page":"Frictional contact with rigid body","title":"Frictional contact with rigid body","text":"This page was generated using Literate.jl.","category":"page"}] +[{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"EditURL = \"../../literate/examples/tlmpm_vortex.jl\"","category":"page"},{"location":"examples/tlmpm_vortex/#Total-Lagrangian-MPM","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"","category":"section"},{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"This example demonstrates the total lagrangian material point method[1]. The implementation solves generalized vortex problem[1] using a linear kernel.","category":"page"},{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"note: Note\nCurrently, the Bernstein function used in the paper[1] has not been implemented.","category":"page"},{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"[1]: de Vaucorbeil, A., Nguyen, V.P. and Hutchinson, C.R., 2020. A Total-Lagrangian Material Point Method for solid mechanics problems involving large deformations. Computer Methods in Applied Mechanics and Engineering, 360, p.112783.","category":"page"},{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"using Tesserae\n\nfunction main()\n\n # Simulation parameters\n h = 0.02 # Grid spacing\n T = 1.0 # Time span\n CFL = 0.1 # Courant number\n α = 0.99 # PIC-FLIP parameter\n\n # Material constants\n E = 1e6 # Young's modulus\n ν = 0.3 # Poisson's ratio\n λ = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter\n μ = E / 2(1 + ν) # Shear modulus\n ρ⁰ = 1e3 # Initial density\n\n # Geometry\n Rᵢ = 0.75\n Rₒ = 1.25\n\n # Equations for vortex\n G = π\n R̄ = (Rᵢ + Rₒ) / 2\n function calc_b_Rθ(R, t)\n local h′′, h′, h = hessian(R -> 1-8((R-R̄)/(Rᵢ-Rₒ))^2+16((R-R̄)/(Rᵢ-Rₒ))^4, R, :all)\n local g′′, g′, g = hessian(t -> G*sin(π*t/T), t, :all)\n β = g * h\n b_R = ( μ/ρ⁰*(3g*h′+R*g*h′′) - R*g′′*h)*sin(β) + (μ/ρ⁰*R*(g*h′)^2 - R*(g′*h)^2)*cos(β)\n b_θ = (-μ/ρ⁰*(3g*h′+R*g*h′′) + R*g′′*h)*cos(β) + (μ/ρ⁰*R*(g*h′)^2 + R*(g′*h)^2)*sin(β)\n Vec(b_R, b_θ)\n end\n isinside(x::Vec) = Rᵢ^2 < x⋅x < Rₒ^2\n\n GridProp = @NamedTuple begin\n X :: Vec{2, Float64}\n m :: Float64\n m⁻¹ :: Float64\n mv :: Vec{2, Float64}\n fint :: Vec{2, Float64}\n fext :: Vec{2, Float64}\n b :: Vec{2, Float64}\n v :: Vec{2, Float64}\n vⁿ :: Vec{2, Float64}\n end\n ParticleProp = @NamedTuple begin\n x :: Vec{2, Float64}\n X :: Vec{2, Float64}\n m :: Float64\n V⁰ :: Float64\n v :: Vec{2, Float64}\n ṽ :: Vec{2, Float64}\n ã :: Vec{2, Float64}\n P :: SecondOrderTensor{2, Float64, 4}\n F :: SecondOrderTensor{2, Float64, 4}\n end\n\n # Background grid\n grid = generate_grid(GridProp, CartesianMesh(h, (-1.5,1.5), (-1.5,1.5)))\n outside_gridinds = findall(!isinside, grid.X)\n\n # Particles\n particles = generate_particles(ParticleProp, grid.X; alg=GridSampling(), spacing=1)\n particles.V⁰ .= volume(grid.X) / length(particles)\n\n filter!(pt->isinside(pt.x), particles)\n\n @. particles.X = particles.x\n @. particles.m = ρ⁰ * particles.V⁰\n @. particles.F = one(particles.F)\n @show length(particles)\n\n # Precompute linear kernel values\n mpvalues = generate_mpvalues(LinearBSpline(), grid.X, length(particles))\n for p in eachindex(particles, mpvalues)\n update!(mpvalues[p], particles.x[p], grid.X)\n end\n\n # Outputs\n outdir = mkpath(joinpath(\"output\", \"tlmpm_vortex\"))\n pvdfile = joinpath(outdir, \"paraview\")\n closepvd(openpvd(pvdfile)) # create file\n\n t = 0.0\n step = 0\n fps = 60\n savepoints = collect(LinRange(t, T, round(Int, T*fps)+1))\n\n Tesserae.@showprogress while t < T\n\n # Calculate timestep based on the wave speed\n vmax = maximum(@. sqrt((λ+2μ) / (particles.m/(particles.V⁰ * det(particles.F)))) + norm(particles.v))\n Δt = CFL * spacing(grid) / vmax\n\n # Compute grid body forces\n for i in eachindex(grid)\n if isinside(grid.X[i])\n (x, y) = grid.X[i]\n R = sqrt(x^2 + y^2)\n θ = atan(y, x)\n grid.b[i] = rotmat(θ) ⋅ calc_b_Rθ(R, t)\n end\n end\n\n # Particle-to-grid transfer\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * v[p]\n fint[i] = @∑ -V⁰[p] * P[p] ⋅ ∇w[ip]\n end\n\n # Update grid velocity\n @. grid.m⁻¹ = inv(grid.m) * !iszero(grid.m)\n @. grid.fext = grid.m * grid.b\n @. grid.vⁿ = grid.mv * grid.m⁻¹\n @. grid.v = grid.vⁿ + Δt * (grid.fint + grid.fext) * grid.m⁻¹\n grid.v[outside_gridinds] .= zero(eltype(grid.v))\n\n # Update particle velocity and position\n @G2P grid=>i particles=>p mpvalues=>ip begin\n ṽ[p] = @∑ w[ip] * v[i]\n ã[p] = @∑ w[ip] * (v[i] - vⁿ[i]) / Δt\n v[p] = (1-α)*ṽ[p] + α*(v[p] + Δt*ã[p])\n x[p] += Δt * ṽ[p]\n end\n\n # Remap updated velocity to grid (MUSL)\n @P2G grid=>i particles=>p mpvalues=>ip begin\n mv[i] = @∑ w[ip] * m[p] * v[p]\n v[i] = mv[i] * m⁻¹[i]\n end\n grid.v[outside_gridinds] .= zero(eltype(grid.v))\n\n # Update stress\n @G2P grid=>i particles=>p mpvalues=>ip begin\n F[p] += @∑ Δt * v[i] ⊗ ∇w[ip]\n P[p] = μ * (F[p] - inv(F[p])') + λ * log(det(F[p])) * inv(F[p])'\n end\n\n t += Δt\n step += 1\n\n if t > first(savepoints)\n popfirst!(savepoints)\n openpvd(pvdfile; append=true) do pvd\n openvtm(string(pvdfile, step)) do vtm\n angle(x) = atan(x[2], x[1])\n openvtk(vtm, particles.x) do vtk\n vtk[\"velocity\"] = particles.v\n vtk[\"initial angle\"] = angle.(particles.X)\n end\n openvtk(vtm, grid.X) do vtk\n vtk[\"external force\"] = grid.fext\n end\n pvd[t] = vtm\n end\n end\n end\n end\nend","category":"page"},{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"","category":"page"},{"location":"examples/tlmpm_vortex/","page":"Total Lagrangian MPM","title":"Total Lagrangian MPM","text":"This page was generated using Literate.jl.","category":"page"},{"location":"examples/implicit_jacobian_based/","page":"Jacobian-based implicit method","title":"Jacobian-based implicit method","text":"EditURL = \"../../literate/examples/implicit_jacobian_based.jl\"","category":"page"},{"location":"examples/implicit_jacobian_based/#Jacobian-based-implicit-method","page":"Jacobian-based implicit method","title":"Jacobian-based implicit method","text":"","category":"section"},{"location":"examples/implicit_jacobian_based/","page":"Jacobian-based implicit method","title":"Jacobian-based implicit method","text":"using Tesserae\n\nfunction main()\n\n # Simulation parameters\n h = 0.05 # Grid spacing\n T = 1.0 # Time span\n g = 20.0 # Gravity acceleration\n Δt = 0.05 # Timestep\n\n # Material constants\n E = 1e6 # Young's modulus\n ν = 0.3 # Poisson's ratio\n λ = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter\n μ = E / 2(1 + ν) # Shear modulus\n ρ⁰ = 500.0 # Density\n\n # Newmark-beta integration\n β = 1/4\n γ = 1/2\n\n GridProp = @NamedTuple begin\n X :: Vec{3, Float64}\n m :: Float64\n m⁻¹ :: Float64\n v :: Vec{3, Float64}\n vⁿ :: Vec{3, Float64}\n mv :: Vec{3, Float64}\n a :: Vec{3, Float64}\n aⁿ :: Vec{3, Float64}\n ma :: Vec{3, Float64}\n u :: Vec{3, Float64}\n f :: Vec{3, Float64}\n end\n ParticleProp = @NamedTuple begin\n x :: Vec{3, Float64}\n m :: Float64\n V⁰ :: Float64\n v :: Vec{3, Float64}\n a :: Vec{3, Float64}\n b :: Vec{3, Float64}\n ∇u :: SecondOrderTensor{3, Float64, 9}\n F :: SecondOrderTensor{3, Float64, 9}\n ΔF⁻¹ :: SecondOrderTensor{3, Float64, 9}\n τ :: SecondOrderTensor{3, Float64, 9}\n c :: FourthOrderTensor{3, Float64, 81}\n end\n\n # Background grid\n grid = generate_grid(GridProp, CartesianMesh(h, (0.0,1.2), (0.0,2.0), (-0.2,0.2)))\n\n # Particles\n beam = Tesserae.Box((0,1), (0.85,1.15), (-0.15,0.15))\n particles = generate_particles(ParticleProp, grid.X; domain=beam, alg=GridSampling())\n particles.V⁰ .= volume(beam) / length(particles)\n @. particles.m = ρ⁰ * particles.V⁰\n @. particles.F = one(particles.F)\n @. particles.b = Vec(0,-g,0)\n @show length(particles)\n\n # Interpolation\n # Use the kernel correction to properly handle the boundary conditions\n it = KernelCorrection(QuadraticBSpline())\n mpvalues = generate_mpvalues(it, grid.X, length(particles))\n\n # Neo-Hookean model\n function kirchhoff_stress(F)\n J = det(F)\n b = symmetric(F ⋅ F')\n μ*(b-I) + λ*log(J)*I\n end\n\n # Sparse matrix\n A = create_sparse_matrix(it, grid.X)\n\n # Outputs\n outdir = mkpath(joinpath(\"output\", \"implicit_jacobian_based\"))\n pvdfile = joinpath(outdir, \"paraview\")\n closepvd(openpvd(pvdfile)) # Create file\n\n t = 0.0\n step = 0\n\n Tesserae.@showprogress while t < T\n\n for p in eachindex(particles, mpvalues)\n update!(mpvalues[p], particles.x[p], grid.X)\n end\n\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * v[p]\n ma[i] = @∑ w[ip] * m[p] * a[p]\n end\n\n # Compute the grid velocity and acceleration at t = tⁿ\n @. grid.m⁻¹ = inv(grid.m) * !iszero(grid.m)\n @. grid.vⁿ = grid.mv * grid.m⁻¹\n @. grid.aⁿ = grid.ma * grid.m⁻¹\n\n # Update a dof map\n dofmask = trues(3, size(grid)...)\n for i in eachindex(grid)\n dofmask[:,i] .= !iszero(grid.m[i])\n end\n for i in eachindex(grid)[1,:,:]\n dofmask[:,i] .= false\n grid.vⁿ[i] = zero(Vec{3})\n grid.aⁿ[i] = zero(Vec{3})\n end\n dofmap = DofMap(dofmask)\n\n # Solve the nonlinear equation\n state = (; grid, particles, mpvalues, kirchhoff_stress, β, γ, A, dofmap, Δt)\n @. grid.u = zero(grid.u) # Set zero dispacement for the first guess of the solution\n U = copy(dofmap(grid.u)) # Convert grid data to plain vector data\n compute_residual(U) = residual(U, state)\n compute_jacobian(U) = jacobian(U, state)\n Tesserae.newton!(U, compute_residual, compute_jacobian)\n\n # Grid dispacement, velocity and acceleration have been updated during Newton's iterations\n @G2P grid=>i particles=>p mpvalues=>ip begin\n ∇u[p] = @∑ u[i] ⊗ ∇w[ip]\n a[p] = @∑ w[ip] * a[i]\n v[p] += @∑ Δt * w[ip] * ((1-γ)*a[p] + γ*a[i])\n x[p] = @∑ w[ip] * (X[i] + u[i])\n F[p] = (I + ∇u[p]) ⋅ F[p]\n end\n\n t += Δt\n step += 1\n\n openpvd(pvdfile; append=true) do pvd\n openvtk(string(pvdfile, step), particles.x) do vtk\n vtk[\"velocity\"] = particles.v\n vtk[\"von Mises\"] = @. vonmises(particles.τ / det(particles.F))\n pvd[t] = vtk\n end\n end\n end\nend\n\nfunction residual(U::AbstractVector, state)\n (; grid, particles, mpvalues, kirchhoff_stress, β, γ, dofmap, Δt) = state\n\n dofmap(grid.u) .= U\n @. grid.a = (1/(β*Δt^2))*grid.u - (1/(β*Δt))*grid.vⁿ - (1/2β-1)*grid.aⁿ\n @. grid.v = grid.vⁿ + Δt*((1-γ)*grid.aⁿ + γ*grid.a)\n\n transposing_tensor(σ) = @einsum (i,j,k,l) -> σ[i,l] * one(σ)[j,k]\n @G2P grid=>i particles=>p mpvalues=>ip begin\n # In addition to updating the stress tensor, the stiffness tensor,\n # which is utilized in the Jacobian-vector product, is also updated.\n ∇u[p] = @∑ u[i] ⊗ ∇w[ip]\n ΔF⁻¹[p] = inv(I + ∇u[p])\n c[p], τ[p] = gradient(∇u -> kirchhoff_stress((I + ∇u) ⋅ F[p]), ∇u[p], :all)\n c[p] = c[p] - transposing_tensor(τ[p] ⋅ ΔF⁻¹[p]')\n end\n @P2G grid=>i particles=>p mpvalues=>ip begin\n f[i] = @∑ -V⁰[p] * τ[p] ⋅ (∇w[ip] ⋅ ΔF⁻¹[p]) + w[ip] * m[p] * b[p]\n end\n\n @. $dofmap(grid.m) * $dofmap(grid.a) - $dofmap(grid.f)\nend\n\nfunction jacobian(U::AbstractVector, state)\n (; grid, particles, mpvalues, β, A, dofmap, Δt) = state\n\n I(i,j) = ifelse(i===j, one(Mat{3,3}), zero(Mat{3,3}))\n dotdot(a,C,b) = @einsum (i,j) -> a[k] * C[i,k,j,l] * b[l]\n @P2G_Matrix grid=>(i,j) particles=>p mpvalues=>(ip,jp) begin\n A[i,j] = @∑ (dotdot(∇w[ip] ⋅ ΔF⁻¹[p], c[p], ∇w[jp])) * V⁰[p] + 1/(β*Δt^2) * I(i,j) * m[p] * w[jp]\n end\n\n submatrix(A, dofmap)\nend","category":"page"},{"location":"examples/implicit_jacobian_based/","page":"Jacobian-based implicit method","title":"Jacobian-based implicit method","text":"","category":"page"},{"location":"examples/implicit_jacobian_based/","page":"Jacobian-based implicit method","title":"Jacobian-based implicit method","text":"This page was generated using Literate.jl.","category":"page"},{"location":"getting_started/#Getting-Started","page":"Getting Started","title":"Getting Started","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"using Tesserae\nimport Plots\n\n# Material constants\nE = 500 # Young's modulus\nν = 0.3 # Poisson's ratio\nλ = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter\nμ = E / 2(1 + ν) # Shear modulus\nρ = 1000 # Density\nr = 0.2 # Radius of disk\n\n# Properties for grid and particles\nstruct GridProp\n x :: Vec{2, Float64}\n m :: Float64\n mv :: Vec{2, Float64}\n f :: Vec{2, Float64}\n v :: Vec{2, Float64}\n vⁿ :: Vec{2, Float64}\nend\nstruct ParticleProp\n x :: Vec{2, Float64}\n m :: Float64\n V :: Float64\n v :: Vec{2, Float64}\n ∇v :: SecondOrderTensor{2, Float64, 4}\n σ :: SymmetricSecondOrderTensor{2, Float64, 3}\nend\n\n# Mesh\nmesh = CartesianMesh(0.05, (0,1), (0,1))\n\n# Background grid\ngrid = generate_grid(GridProp, mesh)\n\n# Particles\nparticles = let\n pts = generate_particles(ParticleProp, mesh; alg=GridSampling())\n pts.V .= volume(mesh) / length(pts)\n\n # Left disk\n lhs = findall(pts.x) do (x,y)\n (x-r)^2 + (y-r)^2 < r^2\n end\n\n # Right disk\n s = 1-r\n rhs = findall(pts.x) do (x,y)\n (x-s)^2 + (y-s)^2 < r^2\n end\n\n pts.v[lhs] .= Vec( 0.1, 0.1)\n pts.v[rhs] .= Vec(-0.1,-0.1)\n \n pts[[lhs; rhs]]\nend\n@. particles.m = ρ * particles.V\n\n# Interpolation\nmpvalues = [MPValue(LinearBSpline(), mesh) for _ in 1:length(particles)]\n\n# Plot results by `Plots.@gif`\nΔt = 0.001\nPlots.@gif for t in range(0, 4-Δt, step=Δt)\n\n # Update interpolation values\n for p in 1:length(particles)\n update!(mpvalues[p], particles[p], mesh)\n end\n\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * v[p]\n f[i] = @∑ -V[p] * σ[p] ⋅ ∇w[ip]\n vⁿ[i] = mv[i] / m[i]\n v[i] = vⁿ[i] + Δt * (f[i] / m[i])\n end\n\n @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] += @∑ w[ip] * (v[i] - vⁿ[i])\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n x[p] += @∑ Δt * (w[ip] * v[i])\n end\n\n for p in 1:length(particles)\n Δϵₚ = Δt * symmetric(particles.∇v[p])\n Δσₚ = λ*tr(Δϵₚ)*I + 2μ*Δϵₚ\n particles.V[p] *= 1 + tr(Δϵₚ)\n particles.σ[p] += Δσₚ\n end\n\n # plot results\n Plots.scatter(\n reinterpret(Tuple{Float64,Float64}, particles.x),\n lims = (0,1),\n ticks = 0:0.2:1,\n minorgrid = true,\n minorticks = 4,\n aspect_ratio = :equal,\n legend = false,\n )\nend every 100","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"EditURL = \"../../literate/examples/elastic_impact.jl\"","category":"page"},{"location":"examples/elastic_impact/#Transfer-schemes","page":"Transfer schemes","title":"Transfer schemes","text":"","category":"section"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"In this example, the following transfer schemes are demonstrated:","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"PIC–FLIP mixed transfer[1]\nAffine PIC (APIC) transfer[2]\nTaylor PIC (TPIC) transfer[3]","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"The problem evolves the elastic impact between two rings, which is consistent with previous studies[4].","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"[1]: Zhu, Y. and Bridson, R., 2005. Animating sand as a fluid. ACM Transactions on Graphics (TOG), 24(3), pp.965-972.","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"[2]: Jiang, C., Schroeder, C., Selle, A., Teran, J. and Stomakhin, A., 2015. The affine particle-in-cell method. ACM Transactions on Graphics (TOG), 34(4), pp.1-10.","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"[3]: Nakamura, K., Matsumura, S. and Mizutani, T., 2023. Taylor particle-in-cell transfer and kernel correction for material point method. Computer Methods in Applied Mechanics and Engineering, 403, p.115720.","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"[4]: Li, X., Fang, Y., Li, M. and Jiang, C., 2022. BFEMP: Interpenetration-free MPM–FEM coupling with barrier contact. Computer Methods in Applied Mechanics and Engineering, 390, p.114350.","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"using Tesserae\n\nabstract type Transfer end\nstruct FLIP <: Transfer α::Float64 end\nstruct APIC <: Transfer end\nstruct TPIC <: Transfer end\n\nfunction main(transfer::Transfer = FLIP(1.0))\n\n # Simulation parameters\n h = 0.1 # Grid spacing\n T = 0.6 # Time span\n CFL = 0.8 # Courant number\n\n # Material constants\n E = 100e6 # Young's modulus\n ν = 0.2 # Poisson's ratio\n λ = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter\n μ = E / 2(1 + ν) # Shear modulus\n ρ⁰ = 1e3 # Initial density\n\n # Geometry\n L = 20.0 # Length of domain\n W = 15.0 # Width of domain\n rᵢ = 3.0 # Inner radius of rings\n rₒ = 4.0 # Outer radius of rings\n\n GridProp = @NamedTuple begin\n x :: Vec{2, Float64}\n m :: Float64\n m⁻¹ :: Float64\n mv :: Vec{2, Float64}\n f :: Vec{2, Float64}\n v :: Vec{2, Float64}\n vⁿ :: Vec{2, Float64}\n end\n ParticleProp = @NamedTuple begin\n x :: Vec{2, Float64}\n m :: Float64\n V⁰ :: Float64\n V :: Float64\n v :: Vec{2, Float64}\n ∇v :: SecondOrderTensor{2, Float64, 4}\n σ :: SymmetricSecondOrderTensor{2, Float64, 3}\n F :: SecondOrderTensor{2, Float64, 4}\n B :: SecondOrderTensor{2, Float64, 4} # for APIC\n end\n\n # Background grid\n grid = generate_grid(GridProp, CartesianMesh(h, (-L,L), (-W/2,W/2)))\n\n # Particles\n particles = let\n pts = generate_particles(ParticleProp, grid.x)\n pts.V .= pts.V⁰ .= volume(grid.x) / length(pts)\n\n lhs = findall(pts.x) do (x, y)\n rᵢ^2 < (x+L/4)^2+y^2 < rₒ^2\n end\n rhs = findall(pts.x) do (x, y)\n rᵢ^2 < (x-L/4)^2+y^2 < rₒ^2\n end\n\n # Set initial velocities\n @. pts.v[lhs] = Vec(30, 0)\n @. pts.v[rhs] = -Vec(30, 0)\n\n pts[[lhs; rhs]]\n end\n @. particles.m = ρ⁰ * particles.V⁰\n @. particles.F = one(particles.F)\n @show length(particles)\n\n # Interpolation\n mpvalues = generate_mpvalues(QuadraticBSpline(), grid.x, length(particles))\n\n # Material model (neo-Hookean)\n function caucy_stress(F)\n b = F ⋅ F'\n J = det(F)\n (μ*(b-I) + λ*log(J)*I) / J\n end\n\n # Outputs\n outdir = mkpath(joinpath(\"output\", \"elastic_impact\"))\n pvdfile = joinpath(outdir, \"paraview\")\n closepvd(openpvd(pvdfile)) # create file\n\n t = 0.0\n step = 0\n fps = 120\n savepoints = collect(LinRange(t, T, round(Int, T*fps)+1))\n\n Tesserae.@showprogress while t < T\n\n # Calculate timestep based on the wave speed\n vmax = maximum(@. sqrt((λ+2μ) / (particles.m/particles.V)) + norm(particles.v))\n Δt = CFL * spacing(grid) / vmax\n\n # Update interpolation values\n for p in eachindex(particles, mpvalues)\n update!(mpvalues[p], particles.x[p], grid.x)\n end\n\n # Particle-to-grid transfer\n if transfer isa FLIP\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * v[p]\n f[i] = @∑ -V[p] * σ[p] ⋅ ∇w[ip]\n end\n elseif transfer isa APIC\n local Dₚ⁻¹ = inv(1/4 * h^2 * I)\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * (v[p] + B[p] ⋅ Dₚ⁻¹ ⋅ (x[i] - x[p]))\n f[i] = @∑ -V[p] * σ[p] ⋅ ∇w[ip]\n end\n elseif transfer isa TPIC\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * (v[p] + ∇v[p] ⋅ (x[i] - x[p]))\n f[i] = @∑ -V[p] * σ[p] ⋅ ∇w[ip]\n end\n end\n\n # Update grid velocity\n @. grid.m⁻¹ = inv(grid.m) * !iszero(grid.m)\n @. grid.vⁿ = grid.mv * grid.m⁻¹\n @. grid.v = grid.vⁿ + Δt * grid.f * grid.m⁻¹\n\n # Grid-to-particle transfer\n if transfer isa FLIP\n local α = transfer.α\n @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] = @∑ w[ip] * ((1-α)*v[i] + α*(v[p] + (v[i]-vⁿ[i])))\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n x[p] += @∑ Δt * (w[ip] * v[i])\n\n end\n elseif transfer isa APIC\n @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] = @∑ w[ip] * v[i]\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n B[p] = @∑ w[ip] * v[i] ⊗ (x[i]-x[p])\n x[p] += Δt * v[p]\n end\n elseif transfer isa TPIC\n @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] = @∑ w[ip] * v[i]\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n x[p] += Δt * v[p]\n end\n end\n\n # Update other particle properties\n for p in eachindex(particles)\n ∇uₚ = Δt * particles.∇v[p]\n Fₚ = (I + ∇uₚ) ⋅ particles.F[p]\n σₚ = caucy_stress(Fₚ)\n particles.σ[p] = σₚ\n particles.F[p] = Fₚ\n particles.V[p] = det(Fₚ) * particles.V⁰[p]\n end\n\n t += Δt\n step += 1\n\n if t > first(savepoints)\n popfirst!(savepoints)\n openpvd(pvdfile; append=true) do pvd\n openvtm(string(pvdfile, step)) do vtm\n function stress3x3(F)\n z = zero(Mat{2,1})\n F3x3 = [F z\n z' 1]\n caucy_stress(F3x3)\n end\n openvtk(vtm, particles.x) do vtk\n vtk[\"velocity\"] = particles.v\n vtk[\"von Mises\"] = @. vonmises(stress3x3(particles.F))\n end\n openvtk(vtm, grid.x) do vtk\n vtk[\"velocity\"] = grid.v\n end\n pvd[t] = vtm\n end\n end\n end\n end\nend","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"","category":"page"},{"location":"examples/elastic_impact/","page":"Transfer schemes","title":"Transfer schemes","text":"This page was generated using Literate.jl.","category":"page"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"EditURL = \"../../literate/examples/dam_break.jl\"","category":"page"},{"location":"examples/dam_break/#Stabilized-mixed-MPM-for-incompressible-fluid","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"","category":"section"},{"location":"examples/dam_break/#Main-function","page":"Stabilized mixed MPM for incompressible fluid","title":"Main function","text":"","category":"section"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"using Tesserae\nusing LinearAlgebra\n\nabstract type Transfer end\nstruct FLIP <: Transfer α::Float64 end\nstruct TPIC <: Transfer end\n\nfunction main(transfer::Transfer = FLIP(1.0))\n\n # Simulation parameters\n h = 0.02 # Grid spacing\n T = 7.0 # Time span\n g = 9.81 # Gravity acceleration\n Δt = 2.5e-3 # Timestep\n\n # Material constants\n ρ = 1.0e3 # Initial density\n μ = 1.01e-3 # Dynamic viscosity (Pa⋅s)\n\n # Newmark-beta method\n β = 0.5\n γ = 1.0\n\n # Utils\n cellnodes(cell) = cell:(cell+oneunit(cell))\n cellcenter(cell, mesh) = mean(mesh[cellnodes(cell)])\n\n # Properties for grid and particles\n GridProp = @NamedTuple begin\n X :: Vec{2, Float64}\n x :: Vec{2, Float64}\n m :: Float64\n m⁻¹ :: Float64\n v :: Vec{2, Float64}\n vⁿ :: Vec{2, Float64}\n mv :: Vec{2, Float64}\n a :: Vec{2, Float64}\n aⁿ :: Vec{2, Float64}\n ma :: Vec{2, Float64}\n u :: Vec{2, Float64}\n p :: Float64\n # δ-correction\n V :: Float64\n Ṽ :: Float64\n E :: Float64\n # Residuals\n u_p :: Vec{3, Float64}\n R_mom :: Vec{2, Float64}\n R_mas :: Float64\n end\n ParticleProp = @NamedTuple begin\n x :: Vec{2, Float64}\n m :: Float64\n V :: Float64\n v :: Vec{2, Float64}\n ∇v :: SecondOrderTensor{2, Float64, 4}\n a :: Vec{2, Float64}\n ∇a :: SecondOrderTensor{2, Float64, 4}\n p :: Float64\n ∇p :: Vec{2, Float64}\n s :: SymmetricSecondOrderTensor{2, Float64, 3}\n b :: Vec{2, Float64}\n # δ-correction\n ∇E² :: Vec{2, Float64}\n # Stabilization\n τ₁ :: Float64\n τ₂ :: Float64\n end\n\n # Background grid\n grid = generate_grid(GridProp, CartesianMesh(h, (0,3.22), (0,2.5)))\n for cell in CartesianIndices(size(grid).-1)\n for i in cellnodes(cell)\n grid.V[i] += (spacing(grid)/2)^2\n end\n end\n\n # Particles\n particles = generate_particles(ParticleProp, grid.X; spacing=1/3)\n particles.V .= volume(grid.X) / length(particles)\n filter!(pt -> pt.x[1]<1.2 && pt.x[2]<0.6, particles)\n @. particles.m = ρ * particles.V\n @. particles.b = Vec(0,-g)\n @show length(particles)\n\n # Interpolation\n it = KernelCorrection(QuadraticBSpline())\n mpvalues = map(p -> MPValue(it, grid.X), eachindex(particles))\n mpvalues_cell = map(CartesianIndices(size(grid).-1)) do cell\n mp = MPValue(it, grid.X)\n xc = cellcenter(cell, grid.X)\n update!(mp, xc, grid.X)\n end\n\n # BlockSpace for threaded computation\n blockspace = BlockSpace(grid.X)\n\n # Sparse matrix\n A = create_sparse_matrix(Vec{3}, it, grid.X)\n\n # Output\n outdir = mkpath(joinpath(\"output\", \"dam_break\"))\n pvdfile = joinpath(outdir, \"paraview\")\n closepvd(openpvd(pvdfile)) # Create file\n\n t = 0.0\n step = 0\n fps = 30\n savepoints = collect(LinRange(t, T, round(Int, T*fps)+1))\n\n Tesserae.@showprogress while t < T\n\n # Update interpolation values based on the nodes of active cells\n # where the particles are located\n activenodes = falses(size(grid))\n for p in eachindex(particles)\n cell = whichcell(particles.x[p], grid.X)\n activenodes[cellnodes(cell)] .= true\n end\n for p in eachindex(mpvalues, particles)\n update!(mpvalues[p], particles.x[p], grid.X, activenodes)\n end\n for cell in CartesianIndices(size(grid) .- 1)\n xc = cellcenter(cell, grid.X)\n update!(mpvalues_cell[cell], xc, grid.X, activenodes)\n end\n\n update!(blockspace, particles.x)\n\n if transfer isa FLIP\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * v[p]\n ma[i] = @∑ w[ip] * m[p] * a[p]\n end\n elseif transfer isa TPIC\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * (v[p] + ∇v[p] ⋅ (X[i] - x[p]))\n ma[i] = @∑ w[ip] * m[p] * (a[p] + ∇a[p] ⋅ (X[i] - x[p]))\n end\n end\n\n # Remove negative mass occasionally generated due to the kernel correction\n # by forcing negative values to zero. While not an ideal solution, this sufficiently\n # reduces spurious pressure oscillations without sacrificing simulation accuracy.\n @. grid.m = max(grid.m, 0)\n\n @. grid.m⁻¹ = inv(grid.m) * !iszero(grid.m)\n @. grid.vⁿ = grid.mv * grid.m⁻¹\n @. grid.aⁿ = grid.ma * grid.m⁻¹\n\n # Update a dof map\n dofmask = trues(3, size(grid)...)\n for i in eachindex(grid)\n dofmask[:,i] .= !iszero(grid.m[i])\n end\n for i in @view eachindex(grid)[[begin,end],:] # Walls\n grid.vⁿ[i] = grid.vⁿ[i] .* (false,true)\n grid.aⁿ[i] = grid.aⁿ[i] .* (false,true)\n dofmask[1,i] = false\n end\n for i in @view eachindex(grid)[:,begin] # Floor\n grid.vⁿ[i] = grid.vⁿ[i] .* (true,false)\n grid.aⁿ[i] = grid.aⁿ[i] .* (true,false)\n dofmask[2,i] = false\n end\n dofmap = DofMap(dofmask)\n\n # Solve grid position, dispacement, velocity, acceleration and pressure by VMS method\n state = (; grid, particles, mpvalues, mpvalues_cell, blockspace, ρ, μ, β, γ, A, dofmap, Δt)\n Δt′ = variational_multiscale_method(state)\n\n if transfer isa FLIP\n local α = transfer.α\n @threaded @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] = @∑ ((1-α)*v[i] + α*(v[p] + Δt*((1-γ)*a[p] + γ*a[i]))) * w[ip]\n a[p] = @∑ a[i] * w[ip]\n x[p] = @∑ x[i] * w[ip]\n end\n elseif transfer isa TPIC\n @threaded @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] = @∑ v[i] * w[ip]\n a[p] = @∑ a[i] * w[ip]\n x[p] = @∑ x[i] * w[ip]\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n ∇a[p] = @∑ a[i] ⊗ ∇w[ip]\n end\n end\n\n # Particle shifting based on the δ-correction\n particle_shifting(state)\n\n # Remove particles that accidentally move outside of the mesh\n outside = findall(x->!isinside(x, grid.X), particles.x)\n deleteat!(particles, outside)\n deleteat!(mpvalues, outside)\n\n t += Δt′\n step += 1\n\n # Write results\n if t > first(savepoints)\n popfirst!(savepoints)\n openpvd(pvdfile; append=true) do pvd\n openvtm(string(pvdfile, step)) do vtm\n vorticity(∇v) = ∇v[2,1] - ∇v[1,2]\n openvtk(vtm, particles.x) do vtk\n vtk[\"pressure\"] = particles.p\n vtk[\"velocity\"] = particles.v\n vtk[\"vorticity\"] = vorticity.(particles.∇v)\n end\n openvtk(vtm, grid.X) do vtk\n end\n pvd[t] = vtm\n end\n end\n end\n end\nend","category":"page"},{"location":"examples/dam_break/#Variational-multiscale-method","page":"Stabilized mixed MPM for incompressible fluid","title":"Variational multiscale method","text":"","category":"section"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"function variational_multiscale_method(state)\n\n # The simulation might fail occasionally due to regions with very small masses,\n # such as splashes[^1]. Therefore, in this script, if Newton's method doesn't converge,\n # a half time step is applied.\n\n (; grid, dofmap, Δt) = state\n @. grid.u_p = zero(grid.u_p)\n\n solved = false\n while !solved\n # Reconstruct state using the current time step\n state = merge(state, (; Δt))\n\n # Compute VMS stabilization coefficients using current grid velocity,\n # which is used for Jacobian matrix\n grid.v .= grid.vⁿ\n compute_VMS_stabilization_coefficients(state)\n\n # Try computing the Jacobian matrix and performing its LU decomposition.\n # In this formulation[^1], the initial Jacobian matrix is used in all Newton's iterations.\n # If the computation fails, use a smaller time step.\n J = lu(jacobian(state); check=false)\n issuccess(J) || (Δt /= 2; continue)\n\n # Solve nonlinear system\n U = zeros(ndofs(dofmap)) # Initialize nodal dispacement and pressure with zero\n solved = Tesserae.newton!(U, U->residual(U,state), U->J;\n linsolve=(x,A,b)->ldiv!(x,A,b), atol=1e-10, rtol=1e-10)\n\n # If the simulation fails to solve, retry with a smaller time step\n solved || (Δt /= 2)\n end\n\n # Update the positions of grid nodes\n @. grid.x = grid.X + grid.u\n\n # Return the acutally applied time step\n Δt\nend","category":"page"},{"location":"examples/dam_break/#VMS-stabilization-coefficients","page":"Stabilized mixed MPM for incompressible fluid","title":"VMS stabilization coefficients","text":"","category":"section"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"function compute_VMS_stabilization_coefficients(state)\n (; grid, particles, mpvalues_cell, ρ, μ, Δt) = state\n\n c₁ = 4.0\n c₂ = 2.0\n τdyn = 1.0\n h = sqrt(4*spacing(grid)^2/π)\n\n # In the following computation, `@G2P` is unavailable\n # due to the use of `mpvalues_cell`\n for p in eachindex(particles)\n v̄ₚ = zero(eltype(particles.v))\n mp = mpvalues_cell[whichcell(particles.x[p], grid.X)]\n gridindices = neighboringnodes(mp, grid)\n for ip in eachindex(gridindices)\n i = gridindices[ip]\n v̄ₚ += mp.w[ip] * grid.v[i]\n end\n τ₁ = inv(ρ*τdyn/Δt + c₂*ρ*norm(v̄ₚ)/h + c₁*μ/h^2)\n τ₂ = h^2 / (c₁*τ₁)\n particles.τ₁[p] = τ₁\n particles.τ₂[p] = τ₂\n end\nend","category":"page"},{"location":"examples/dam_break/#δ-correction","page":"Stabilized mixed MPM for incompressible fluid","title":"δ-correction","text":"","category":"section"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"function particle_shifting(state)\n (; grid, particles, mpvalues) = state\n\n @P2G grid=>i particles=>p mpvalues=>ip begin\n Ṽ[i] = @∑ V[p] * w[ip]\n E[i] = max(0, -V[i] + Ṽ[i])\n end\n\n E² = sum(E->E^2, grid.E)\n @G2P grid=>i particles=>p mpvalues=>ip begin\n ∇E²[p] = @∑ 2V[p] * E[i] * ∇w[ip]\n end\n\n b₀ = E² / sum(∇E²->∇E²⋅∇E², particles.∇E²)\n\n for p in eachindex(particles)\n xₚ = particles.x[p] - b₀ * particles.∇E²[p]\n if isinside(xₚ, grid.X)\n particles.x[p] = xₚ\n end\n end\nend","category":"page"},{"location":"examples/dam_break/#Residual-vector","page":"Stabilized mixed MPM for incompressible fluid","title":"Residual vector","text":"","category":"section"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"function residual(U, state)\n (; grid, particles, mpvalues, μ, β, γ, dofmap, Δt) = state\n\n # Map `U` to grid dispacement and pressure\n dofmap(grid.u_p) .= U\n grid.u .= map(x->@Tensor(x[1:2]), grid.u_p)\n grid.p .= map(x->x[3], grid.u_p)\n\n # Recompute nodal velocity and acceleration based on the Newmark-beta method\n @. grid.v = γ/(β*Δt)*grid.u - (γ/β-1)*grid.vⁿ - Δt/2*(γ/β-2)*grid.aⁿ\n @. grid.a = 1/(β*Δt^2)*grid.u - 1/(β*Δt)*grid.vⁿ - (1/2β-1)*grid.aⁿ\n\n # Recompute particle properties for residual vector\n @G2P grid=>i particles=>p mpvalues=>ip begin\n a[p] = @∑ a[i] * w[ip]\n p[p] = @∑ p[i] * w[ip]\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n ∇p[p] = @∑ p[i] * ∇w[ip]\n s[p] = 2μ * symmetric(∇v[p])\n end\n\n # Compute VMS stabilization coefficients based on the current nodal velocity\n compute_VMS_stabilization_coefficients(state)\n\n # Compute residual values\n @P2G grid=>i particles=>p mpvalues=>ip begin\n R_mom[i] = @∑ V[p]*s[p]⋅∇w[ip] - m[p]*b[p]*w[ip] - V[p]*p[p]*∇w[ip] + τ₂[p]*V[p]*tr(∇v[p])*∇w[ip]\n R_mas[i] = @∑ V[p]*tr(∇v[p])*w[ip] + τ₁[p]*m[p]*(a[p]-b[p])⋅∇w[ip] + τ₁[p]*V[p]*∇p[p]⋅∇w[ip]\n R_mom[i] += m[i]*a[i]\n end\n\n # Map grid values to vector `R`\n dofmap(map(vcat, grid.R_mom, grid.R_mas))\nend","category":"page"},{"location":"examples/dam_break/#Jacobian-matrix","page":"Stabilized mixed MPM for incompressible fluid","title":"Jacobian matrix","text":"","category":"section"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"function jacobian(state)\n (; grid, particles, mpvalues, blockspace, ρ, μ, β, γ, A, dofmap, Δt) = state\n\n # Construct the Jacobian matrix\n cₚ = 2μ * one(SymmetricFourthOrderTensor{2})\n I(i,j) = ifelse(i===j, one(Mat{2,2}), zero(Mat{2,2}))\n @threaded @P2G_Matrix grid=>(i,j) particles=>p mpvalues=>(ip,jp) blockspace begin\n A[i,j] = @∑ begin\n Kᵤᵤ = (γ/(β*Δt) * ∇w[ip] ⋅ cₚ ⋅ ∇w[jp]) * V[p] + 1/(β*Δt^2) * I(i,j) * m[p] * w[jp]\n Kᵤₚ = -∇w[ip] * w[jp] * V[p]\n Kₚᵤ = (γ/(β*Δt)) * w[ip] * ∇w[jp] * V[p]\n K̂ᵤᵤ = γ/(β*Δt) * τ₂[p] * ∇w[ip] ⊗ ∇w[jp] * V[p]\n K̂ₚᵤ = 1/(β*Δt^2) * τ₁[p] * ρ * ∇w[ip] * w[jp] * V[p]\n K̂ₚₚ = τ₁[p] * ∇w[ip] ⋅ ∇w[jp] * V[p]\n [Kᵤᵤ+K̂ᵤᵤ Kᵤₚ\n Mat{1,2}(Kₚᵤ+K̂ₚᵤ) K̂ₚₚ]\n end\n end\n\n # Extract the activated degrees of freedom\n submatrix(A, dofmap)\nend","category":"page"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"","category":"page"},{"location":"examples/dam_break/","page":"Stabilized mixed MPM for incompressible fluid","title":"Stabilized mixed MPM for incompressible fluid","text":"This page was generated using Literate.jl.","category":"page"},{"location":"#Tesserae","page":"Home","title":"Tesserae","text":"","category":"section"},{"location":"examples/implicit_jacobian_free/","page":"Jacobian-free Newton–Krylov method","title":"Jacobian-free Newton–Krylov method","text":"EditURL = \"../../literate/examples/implicit_jacobian_free.jl\"","category":"page"},{"location":"examples/implicit_jacobian_free/#Jacobian-free-Newton–Krylov-method","page":"Jacobian-free Newton–Krylov method","title":"Jacobian-free Newton–Krylov method","text":"","category":"section"},{"location":"examples/implicit_jacobian_free/","page":"Jacobian-free Newton–Krylov method","title":"Jacobian-free Newton–Krylov method","text":"using Tesserae\n\nusing IterativeSolvers: gmres!\nusing LinearMaps: LinearMap\n\nfunction main()\n\n # Simulation parameters\n h = 0.05 # Grid spacing\n T = 1.0 # Time span\n g = 20.0 # Gravity acceleration\n Δt = 0.02 # Timestep\n\n # Material constants\n E = 1e6 # Young's modulus\n ν = 0.3 # Poisson's ratio\n λ = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter\n μ = E / 2(1 + ν) # Shear modulus\n ρ⁰ = 500.0 # Density\n\n # Newmark-beta integration\n β = 1/4\n γ = 1/2\n\n GridProp = @NamedTuple begin\n X :: Vec{3, Float64}\n m :: Float64\n m⁻¹ :: Float64\n v :: Vec{3, Float64}\n vⁿ :: Vec{3, Float64}\n mv :: Vec{3, Float64}\n a :: Vec{3, Float64}\n aⁿ :: Vec{3, Float64}\n ma :: Vec{3, Float64}\n u :: Vec{3, Float64}\n f :: Vec{3, Float64}\n end\n ParticleProp = @NamedTuple begin\n x :: Vec{3, Float64}\n m :: Float64\n V⁰ :: Float64\n v :: Vec{3, Float64}\n a :: Vec{3, Float64}\n b :: Vec{3, Float64}\n ∇u :: SecondOrderTensor{3, Float64, 9}\n F :: SecondOrderTensor{3, Float64, 9}\n ΔF⁻¹ :: SecondOrderTensor{3, Float64, 9}\n τ :: SecondOrderTensor{3, Float64, 9}\n c :: FourthOrderTensor{3, Float64, 81}\n end\n\n # Background grid\n grid = generate_grid(GridProp, CartesianMesh(h, (0.0,1.2), (0.0,2.0), (-0.2,0.2)))\n\n # Particles\n beam = Tesserae.Box((0,1), (0.85,1.15), (-0.15,0.15))\n particles = generate_particles(ParticleProp, grid.X; domain=beam, alg=GridSampling())\n particles.V⁰ .= volume(beam) / length(particles)\n @. particles.m = ρ⁰ * particles.V⁰\n @. particles.F = one(particles.F)\n @. particles.b = Vec(0,-g,0)\n @show length(particles)\n\n # Interpolation\n # Use the kernel correction to properly handle the boundary conditions\n mpvalues = generate_mpvalues(KernelCorrection(QuadraticBSpline()), grid.X, length(particles))\n\n # Neo-Hookean model\n function kirchhoff_stress(F)\n J = det(F)\n b = symmetric(F ⋅ F')\n μ*(b-I) + λ*log(J)*I\n end\n\n # Outputs\n outdir = mkpath(joinpath(\"output\", \"implicit_jacobian_free\"))\n pvdfile = joinpath(outdir, \"paraview\")\n closepvd(openpvd(pvdfile)) # Create file\n\n t = 0.0\n step = 0\n\n Tesserae.@showprogress while t < T\n\n for p in eachindex(particles, mpvalues)\n update!(mpvalues[p], particles.x[p], grid.X)\n end\n\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * v[p]\n ma[i] = @∑ w[ip] * m[p] * a[p]\n end\n\n # Compute the grid velocity and acceleration at t = tⁿ\n @. grid.m⁻¹ = inv(grid.m) * !iszero(grid.m)\n @. grid.vⁿ = grid.mv * grid.m⁻¹\n @. grid.aⁿ = grid.ma * grid.m⁻¹\n\n # Update a dof map\n dofmask = trues(3, size(grid)...)\n for i in eachindex(grid)\n dofmask[:,i] .= !iszero(grid.m[i])\n end\n for i in eachindex(grid)[1,:,:]\n dofmask[:,i] .= false\n grid.vⁿ[i] = zero(Vec{3})\n grid.aⁿ[i] = zero(Vec{3})\n end\n dofmap = DofMap(dofmask)\n\n # Solve the nonlinear equation\n state = (; grid, particles, mpvalues, kirchhoff_stress, β, γ, dofmap, Δt)\n @. grid.u = zero(grid.u) # Set zero dispacement for the first guess of the solution\n U = copy(dofmap(grid.u)) # Convert grid data to plain vector data\n compute_residual(U) = residual(U, state)\n compute_jacobian(U) = jacobian(U, state)\n Tesserae.newton!(U, compute_residual, compute_jacobian; linsolve = (x,A,b)->gmres!(x,A,b))\n\n # Grid dispacement, velocity and acceleration have been updated during Newton's iterations\n @G2P grid=>i particles=>p mpvalues=>ip begin\n ∇u[p] = @∑ u[i] ⊗ ∇w[ip]\n a[p] = @∑ w[ip] * a[i]\n v[p] += @∑ Δt * w[ip] * ((1-γ)*a[p] + γ*a[i])\n x[p] = @∑ w[ip] * (X[i] + u[i])\n F[p] = (I + ∇u[p]) ⋅ F[p]\n end\n\n t += Δt\n step += 1\n\n openpvd(pvdfile; append=true) do pvd\n openvtk(string(pvdfile, step), particles.x) do vtk\n vtk[\"velocity\"] = particles.v\n vtk[\"von Mises\"] = @. vonmises(particles.τ / det(particles.F))\n pvd[t] = vtk\n end\n end\n end\nend\n\nfunction residual(U::AbstractVector, state)\n (; grid, particles, mpvalues, kirchhoff_stress, β, γ, dofmap, Δt) = state\n\n dofmap(grid.u) .= U\n @. grid.a = (1/(β*Δt^2))*grid.u - (1/(β*Δt))*grid.vⁿ - (1/2β-1)*grid.aⁿ\n @. grid.v = grid.vⁿ + Δt*((1-γ)*grid.aⁿ + γ*grid.a)\n\n transposing_tensor(σ) = @einsum (i,j,k,l) -> σ[i,l] * one(σ)[j,k]\n @G2P grid=>i particles=>p mpvalues=>ip begin\n # In addition to updating the stress tensor, the stiffness tensor,\n # which is utilized in the Jacobian-vector product, is also updated.\n ∇u[p] = @∑ u[i] ⊗ ∇w[ip]\n ΔF⁻¹[p] = inv(I + ∇u[p])\n c[p], τ[p] = gradient(∇u -> kirchhoff_stress((I + ∇u) ⋅ F[p]), ∇u[p], :all)\n c[p] = c[p] - transposing_tensor(τ[p] ⋅ ΔF⁻¹[p]')\n end\n @P2G grid=>i particles=>p mpvalues=>ip begin\n f[i] = @∑ -V⁰[p] * τ[p] ⋅ (∇w[ip] ⋅ ΔF⁻¹[p]) + w[ip] * m[p] * b[p]\n end\n\n @. β*Δt^2 * ($dofmap(grid.a) - $dofmap(grid.f) * $dofmap(grid.m⁻¹))\nend\n\nfunction jacobian(U::AbstractVector, state)\n (; grid, particles, mpvalues, β, dofmap, Δt) = state\n\n # Create a linear map to represent Jacobian-vector product J*δU.\n # `U` is acutally not used because the stiffness tensor is already calculated\n # when computing the residual vector.\n LinearMap(ndofs(dofmap)) do JδU, δU\n dofmap(grid.u) .= δU\n\n @G2P grid=>i particles=>p mpvalues=>ip begin\n ∇u[p] = @∑ u[i] ⊗ ∇w[ip]\n τ[p] = c[p] ⊡ ∇u[p]\n end\n @P2G grid=>i particles=>p mpvalues=>ip begin\n f[i] = @∑ -V⁰[p] * τ[p] ⋅ (∇w[ip] ⋅ ΔF⁻¹[p])\n end\n\n @. JδU = δU - β*Δt^2 * $dofmap(grid.f) * $dofmap(grid.m⁻¹)\n end\nend","category":"page"},{"location":"examples/implicit_jacobian_free/","page":"Jacobian-free Newton–Krylov method","title":"Jacobian-free Newton–Krylov method","text":"","category":"page"},{"location":"examples/implicit_jacobian_free/","page":"Jacobian-free Newton–Krylov method","title":"Jacobian-free Newton–Krylov method","text":"This page was generated using Literate.jl.","category":"page"},{"location":"examples/rigid_body_contact/","page":"Frictional contact with rigid body","title":"Frictional contact with rigid body","text":"EditURL = \"../../literate/examples/rigid_body_contact.jl\"","category":"page"},{"location":"examples/rigid_body_contact/#Frictional-contact-with-rigid-body","page":"Frictional contact with rigid body","title":"Frictional contact with rigid body","text":"","category":"section"},{"location":"examples/rigid_body_contact/","page":"Frictional contact with rigid body","title":"Frictional contact with rigid body","text":"using Tesserae\n\nmutable struct Disk{dim, T}\n x::Vec{dim, T}\n v::Vec{dim, T}\nend\n\nfunction main()\n\n # Simulation parameters\n h = 0.004 # Grid spacing\n T = 5.0 # Time span\n g = 9.81 # Gravity acceleration\n CFL = 1.0 # Courant number\n\n # Material constants\n E = 1000e3 # Young's modulus\n ν = 0.49 # Poisson's ratio\n λ = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter\n G = E / 2(1 + ν) # Shear modulus\n σy = 1.0e3 # Yield stress\n ρ⁰ = 1.0e3 # Initial density\n\n # Disk\n D = 0.04\n disk = Disk(Vec(0, 7.5D), Vec(0, -0.25D))\n\n # Contact parameters\n k = 1e6 # Penalty coefficient\n μ = 0.6 # Friction coefficient\n\n # Utils\n @inline function contact_force_normal(x, r, x_disk)\n d = x - x_disk\n k * max(D/2 - (norm(d)-r), 0) * normalize(d)\n end\n @inline function contact_force_tangential(fₙ, v, m, Δt)\n iszero(fₙ) && return zero(fₙ)\n n = normalize(fₙ)\n fₜ = -m * (v-(v⋅n)*n) / Δt # Sticking force\n min(1, μ*norm(fₙ)/norm(fₜ)) * fₜ\n end\n\n # Properties for grid and particles\n GridProp = @NamedTuple begin\n x :: Vec{2, Float64}\n m :: Float64\n m⁻¹ :: Float64\n mv :: Vec{2, Float64}\n fint :: Vec{2, Float64}\n fext :: Vec{2, Float64}\n v :: Vec{2, Float64}\n vⁿ :: Vec{2, Float64}\n end\n ParticleProp = @NamedTuple begin\n x :: Vec{2, Float64}\n m :: Float64\n V :: Float64\n r :: Float64\n v :: Vec{2, Float64}\n ∇v :: SecondOrderTensor{2, Float64, 4}\n F :: SecondOrderTensor{3, Float64, 9}\n σ :: SymmetricSecondOrderTensor{3, Float64, 6}\n ϵ :: SymmetricSecondOrderTensor{3, Float64, 6}\n b :: Vec{2, Float64}\n end\n\n # Background grid\n H = 7D # Ground height\n grid = generate_grid(GridProp, CartesianMesh(h, (0,5D), (0,H+D)))\n\n # Particles\n particles = generate_particles(ParticleProp, grid.x)\n disk_points = filter(particles.x) do (x,y) # Points representing a disk just for visualization\n x^2 + (y-7.5D)^2 < (D/2)^2\n end\n particles.V .= volume(grid.x) / length(particles)\n filter!(pt -> pt.x[2] < H, particles)\n @. particles.m = ρ⁰ * particles.V\n @. particles.r = particles.V^(1/2) / 2\n @. particles.b = Vec(0,-g)\n @. particles.F = one(particles.F)\n for p in eachindex(particles)\n x, y = particles.x[p]\n σ_y = -ρ⁰ * g * (H - y)\n σ_x = ν/(1-ν) * σ_y\n particles.σ[p] = [σ_x 0.0 0.0\n 0.0 σ_y 0.0\n 0.0 0.0 σ_x]\n end\n @show length(particles)\n\n # Interpolation\n mpvalues = generate_mpvalues(KernelCorrection(QuadraticBSpline()), grid.x, length(particles))\n\n # Output\n outdir = mkpath(joinpath(\"output\", \"rigid_body_contact\"))\n pvdfile = joinpath(outdir, \"paraview\")\n closepvd(openpvd(pvdfile)) # Create file\n\n t = 0.0\n step = 0\n fps = 20\n savepoints = collect(LinRange(t, T, round(Int, T*fps)+1))\n\n Tesserae.@showprogress while t < T\n\n vmax = maximum(@. sqrt((λ+2G) / (particles.m/particles.V)) + norm(particles.v))\n Δt = CFL * spacing(grid) / vmax\n\n for p in eachindex(particles, mpvalues)\n update!(mpvalues[p], particles.x[p], grid.x)\n end\n\n @P2G grid=>i particles=>p mpvalues=>ip begin\n m[i] = @∑ w[ip] * m[p]\n mv[i] = @∑ w[ip] * m[p] * (v[p] + ∇v[p] ⋅ (x[i] - x[p])) # Taylor transfer\n fint[i] = @∑ -V[p] * resizedim(σ[p], Val(2)) ⋅ ∇w[ip] + w[ip] * m[p] * b[p]\n fext[i] = @∑ w[ip] * contact_force_normal(x[p], r[p], disk.x)\n end\n\n @. grid.m⁻¹ = inv(grid.m) * !iszero(grid.m)\n @. grid.vⁿ = grid.mv * grid.m⁻¹\n @. grid.v = grid.vⁿ + Δt * grid.fint * grid.m⁻¹\n @. grid.fext += contact_force_tangential(grid.fext, grid.v-disk.v, grid.m, Δt)\n @. grid.v += Δt * grid.fext * grid.m⁻¹\n\n for i in eachindex(grid)[[begin,end],:]\n grid.v[i] = grid.v[i] .* (false,true)\n end\n for i in eachindex(grid)[:,[begin,end]]\n grid.v[i] = grid.v[i] .* (false,false)\n end\n\n @G2P grid=>i particles=>p mpvalues=>ip begin\n v[p] = @∑ w[ip] * v[i] # PIC transfer\n ∇v[p] = @∑ v[i] ⊗ ∇w[ip]\n x[p] += Δt * v[p]\n end\n\n for p in 1:length(particles)\n Δϵ = resizedim(Δt * symmetric(particles.∇v[p]), Val(3))\n particles.σ[p] = vonmises_model(particles.σ[p], Δϵ; λ, G, σy)\n particles.V[p] *= 1 + tr(Δϵ)\n particles.ϵ[p] += Δϵ\n end\n\n disk.x += disk.v * Δt\n @. disk_points += disk.v * Δt\n\n t += Δt\n step += 1\n\n if t > first(savepoints)\n popfirst!(savepoints)\n openpvd(pvdfile; append=true) do pvd\n openvtm(string(pvdfile, step)) do vtm\n deviatoric_strain(ϵ) = sqrt(2/3 * dev(ϵ) ⊡ dev(ϵ))\n openvtk(vtm, particles.x) do vtk\n vtk[\"vonmises\"] = @. vonmises(particles.σ)\n vtk[\"deviatoric strain\"] = @. deviatoric_strain(particles.ϵ)\n end\n openvtk(vtm, disk_points) do vtk\n end\n pvd[t] = vtm\n end\n end\n end\n end\nend\n\nfunction vonmises_model(σⁿ, Δϵ; λ, G, σy)\n δ = one(SymmetricSecondOrderTensor{3})\n I = one(SymmetricFourthOrderTensor{3})\n cᵉ = λ*δ⊗δ + 2G*I\n σ_trial = σⁿ + cᵉ ⊡ Δϵ\n dfdσ, f_trial = gradient(σ -> vonmises(σ) - σy, σ_trial, :all)\n if f_trial > 0\n dλ = f_trial / (dfdσ ⊡ cᵉ ⊡ dfdσ)\n σ = σ_trial - cᵉ ⊡ (dλ * dfdσ)\n else\n σ = σ_trial\n end\n if mean(σ) > 0 # simple tension cut-off\n σ = dev(σ)\n end\n σ\nend","category":"page"},{"location":"examples/rigid_body_contact/","page":"Frictional contact with rigid body","title":"Frictional contact with rigid body","text":"","category":"page"},{"location":"examples/rigid_body_contact/","page":"Frictional contact with rigid body","title":"Frictional contact with rigid body","text":"This page was generated using Literate.jl.","category":"page"}] }