Implement 3D Axisymmetric kernels on a 1D mesh

Project.toml

@@ -11,6 +11,7 @@ IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"

src/DDMFaultSlip.jl

@@ -2,6 +2,7 @@ module DDMFaultSlip
 
 using StaticArrays
 using SparseArrays
+using SpecialFunctions
 using LinearAlgebra
 using Statistics: mean
 using InteractiveUtils
@@ -13,7 +14,7 @@ using WriteVTK
 using DelimitedFiles
 
 include("mesh.jl")
-export Point2D, Point3D, DDMesh1D, DDMesh2D
+export Point2D, Point3D, DDMesh, DDMesh1D, DDMesh2D
 
 include("collocation.jl")
 export DD3DShearElasticMatrix

src/assembly.jl

@@ -153,7 +153,7 @@ function update!(problem::NormalDDProblem{T}, solver::DDSolver{R,T}) where {R,T<
     problem.σ.value = problem.σ.value_old + mul!(zeros(T, size(solver.solution)), solver.E, solver.solution; global_index=true, threads=false)
     # Fluid coupling
     if hasFluidCoupling(problem)
-        problem.σ.value -= problem.fluid_coupling.p + problem.fluid_coupling.p_old
+        problem.σ.value -= problem.fluid_coupling.p - problem.fluid_coupling.p_old
     end
 
     return nothing

src/collocation.jl

@@ -1,8 +1,8 @@
-function F(e::Point3D{T}, r::Point3D{T}) where {T<:Real}
+function F(e::Point3D{T}, r::Point3D{T})::T where {T<:Real}
     return (r[1] * e[1] + r[2] * e[2]) / norm(r)
 end
 
-function integralI003(ei::DDTriangleElem{T}, ej::DDTriangleElem{T}) where {T<:Real}
+function integralI003(ei::DDTriangleElem{T}, ej::DDTriangleElem{T})::T where {T<:Real}
     I = 0.0
     for p in 1:3
         ap = ej.nodes[p]
@@ -11,16 +11,12 @@
         else
             bp = ej.nodes[1]
         end
-        # ra = ei.X - ap
-        # rb = ei.X - bp
-        # ep = bp - ap
-        # I += 1.0 / (ra[2] * ep[1] - ra[1] * ep[2]) * (F(ep, rb) - F(ep, ra))
         I += 1.0 / ((ei.X - ap)[2] * (bp - ap)[1] - (ei.X - ap)[1] * (bp - ap)[2]) * (F(bp - ap, ei.X - bp) - F(bp - ap, ei.X - ap))
     end
     return I
 end
 
-function integralI205(ei::DDTriangleElem{T}, ej::DDTriangleElem{T}) where {T<:Real}
+function integralI205(ei::DDTriangleElem{T}, ej::DDTriangleElem{T})::T where {T<:Real}
     I = 0.0
     for p in 1:3
         ap = ej.nodes[p]
@@ -29,16 +25,12 @@
         else
             bp = ej.nodes[1]
         end
-        # ra = ei.X - ap
-        # rb = ei.X - bp
-        # ep = bp - ap
-        # I += (ra[2] * ep[1] - 2.0 * ra[1] * ep[2]) / (3 * (ra[2] * ep[1] - ra[1] * ep[2])^2) * (F(ep, rb) - F(ep, ra)) + ep[1] * ep[2] / (3 * (ra[2] * ep[1] - ra[1] * ep[2])^2) * (F(ra, rb) - F(ra, ra))
         I += ((ei.X - ap)[2] * (bp - ap)[1] - 2.0 * (ei.X - ap)[1] * (bp - ap)[2]) / (3 * ((ei.X - ap)[2] * (bp - ap)[1] - (ei.X - ap)[1] * (bp - ap)[2])^2) * (F(bp - ap, ei.X - bp) - F(bp - ap, ei.X - ap)) + (bp - ap)[1] * (bp - ap)[2] / (3 * ((ei.X - ap)[2] * (bp - ap)[1] - (ei.X - ap)[1] * (bp - ap)[2])^2) * (F(ei.X - ap, ei.X - bp) - F(ei.X - ap, ei.X - ap))
     end
     return I
 end
 
-function integralI025(ei::DDTriangleElem{T}, ej::DDTriangleElem{T}) where {T<:Real}
+function integralI025(ei::DDTriangleElem{T}, ej::DDTriangleElem{T})::T where {T<:Real}
     I = 0.0
     for p in 1:3
         ap = ej.nodes[p]
@@ -47,16 +39,12 @@
         else
             bp = ej.nodes[1]
         end
-        # ra = ei.X - ap
-        # rb = ei.X - bp
-        # ep = bp - ap
-        # I += (2.0 * ra[2] * ep[1] - ra[1] * ep[2]) / (3 * (ra[2] * ep[1] - ra[1] * ep[2])^2) * (F(ep, rb) - F(ep, ra)) - ep[1] * ep[2] / (3 * (ra[2] * ep[1] - ra[1] * ep[2])^2) * (F(ra, rb) - F(ra, ra))
         I += (2.0 * (ei.X - ap)[2] * (bp - ap)[1] - (ei.X - ap)[1] * (bp - ap)[2]) / (3 * ((ei.X - ap)[2] * (bp - ap)[1] - (ei.X - ap)[1] * (bp - ap)[2])^2) * (F(bp - ap, ei.X - bp) - F(bp - ap, ei.X - ap)) - (bp - ap)[1] * (bp - ap)[2] / (3 * ((ei.X - ap)[2] * (bp - ap)[1] - (ei.X - ap)[1] * (bp - ap)[2])^2) * (F(ei.X - ap, ei.X - bp) - F(ei.X - ap, ei.X - ap))
     end
     return I
 end
 
-function integralI115(ei::DDTriangleElem{T}, ej::DDTriangleElem{T}) where {T<:Real}
+function integralI115(ei::DDTriangleElem{T}, ej::DDTriangleElem{T})::T where {T<:Real}
     I = 0.0
     for p in 1:3
         ap = ej.nodes[p]
@@ -65,14 +53,15 @@
         else
             bp = ej.nodes[1]
         end
-        ra = ei.X - ap
-        rb = ei.X - bp
-        ep = bp - ap
-        I += (ra[1] * ep[1] - ra[2] * ep[2]) / (6 * (ra[2] * ep[1] - ra[1] * ep[2])^2) * (F(ep, rb) - F(ep, ra)) + (ep[2]^2 - ep[1]^2) / (6 * (ra[2] * ep[1] - ra[1] * ep[2])^2) * (F(ra, rb) - F(ra, ra))
+        I += ((ei.X - ap)[1] * (bp - ap)[1] - (ei.X - ap)[2] * (bp - ap)[2]) / (6 * ((ei.X - ap)[2] * (bp - ap)[1] - (ei.X - ap)[1] * (bp - ap)[2])^2) * (F(bp - ap, ei.X - bp) - F(bp - ap, ei.X - ap)) + ((bp - ap)[2]^2 - (bp - ap)[1]^2) / (6 * ((ei.X - ap)[2] * (bp - ap)[1] - (ei.X - ap)[1] * (bp - ap)[2])^2) * (F(ei.X - ap, ei.X - bp) - F(ei.X - ap, ei.X - ap))
     end
     return I
 end
 
+function m(ri::Point2D{T}, rj::Point2D{T}) where {T<:Real}
+    return 4.0 * norm(ri) * norm(rj) / (norm(ri) + norm(rj))^2
+end
+
 abstract type ElasticKernelMatrix{T<:Real} <: AbstractMatrix{T} end
 
 """
@@ -97,6 +86,28 @@
     return K.n, K.n
 end
 
+"""
+Piecewise constant (PWC) three-dimensional axisymmetric elastic kernel matrix
+"""
+struct DD3DAxisymmetricElasticMatrix{T<:Real} <: ElasticKernelMatrix{T}
+    e::Vector{DDEdgeElem{T}}
+    μ::T
+    n::Int
+
+    # Constructor
+    function DD3DAxisymmetricElasticMatrix(mesh::DDMesh1D{T}, μ::T) where {T<:Real}
+        return new{T}(mesh.elems, μ, length(mesh.elems))
+    end
+end
+
+function Base.getindex(K::DD3DAxisymmetricElasticMatrix, i::Int, j::Int)
+    return K.μ / (2 * π) * (ellipe(m(K.e[i].X, K.e[j].nodes[1])) / (norm(K.e[j].nodes[1]) - norm(K.e[i].X)) - ellipe(m(K.e[i].X, K.e[j].nodes[2])) / (norm(K.e[j].nodes[2]) - norm(K.e[i].X)) + ellipk(m(K.e[i].X, K.e[j].nodes[1])) / (norm(K.e[j].nodes[1]) + norm(K.e[i].X)) - ellipk(m(K.e[i].X, K.e[j].nodes[2])) / (norm(K.e[j].nodes[2]) + norm(K.e[i].X)))
+end
+
+function Base.size(K::DD3DAxisymmetricElasticMatrix)
+    return K.n, K.n
+end
+
 """
 Piecewise constant (PWC) three-dimensional normal elastic kernel matrix
 """
@@ -158,22 +169,22 @@
 """
 Piecewise constant (PWC) three-dimensional shear axis symmetric elastic kernel matrix
 """
-struct DD3DShearAxisSymmetricElasticMatrix{T<:Real} <: ElasticKernelMatrix{T}
+struct DD3DShearAxisymmetricElasticMatrix{T<:Real} <: ElasticKernelMatrix{T}
     e::Vector{DDTriangleElem{T}}
     μ::T
     n::Int
 
     # Constructor
-    function DD3DShearAxisSymmetricElasticMatrix(mesh::DDMesh2D{T}, μ::T) where {T<:Real}
+    function DD3DShearAxisymmetricElasticMatrix(mesh::DDMesh2D{T}, μ::T) where {T<:Real}
         return new{T}(mesh.elems, μ, length(mesh.elems))
     end
 end
 
-function Base.getindex(K::DD3DShearAxisSymmetricElasticMatrix, i::Int, j::Int)
+function Base.getindex(K::DD3DShearAxisymmetricElasticMatrix, i::Int, j::Int)
     return K.μ / (4 * π) * integralI003(K.e[i], K.e[j])
 end
 
-function Base.size(K::DD3DShearAxisSymmetricElasticMatrix)
+function Base.size(K::DD3DShearAxisymmetricElasticMatrix)
     return K.n, K.n
 end
 
@@ -204,6 +215,33 @@
     return K.n, K.n
 end
 
+"""
+Piecewise constant (PWC) three-dimensional axisymmetric jacobian matrix
+"""
+struct DD3DAxisymmetricJacobianMatrix{T<:Real} <: ElasticKernelMatrix{T}
+    e::Vector{DDEdgeElem{T}}
+    mat_loc::Matrix{Vector{T}}
+    μ::T
+    n::Int
+
+    # Constructor
+    function DD3DAxisymmetricJacobianMatrix(mesh::DDMesh1D{T}, mat_loc::Matrix{Vector{T}}, μ::T) where {T<:Real}
+        return new{T}(mesh.elems, mat_loc, μ, length(mesh.elems))
+    end
+end
+
+function Base.getindex(K::DD3DAxisymmetricJacobianMatrix, i::Int, j::Int)
+    if (i == j)
+        return K.μ / (2 * π) * (ellipe(m(K.e[i].X, K.e[j].nodes[1])) / (norm(K.e[j].nodes[1]) - norm(K.e[i].X)) - ellipe(m(K.e[i].X, K.e[j].nodes[2])) / (norm(K.e[j].nodes[2]) - norm(K.e[i].X)) + ellipk(m(K.e[i].X, K.e[j].nodes[1])) / (norm(K.e[j].nodes[1]) + norm(K.e[i].X)) - ellipk(m(K.e[i].X, K.e[j].nodes[2])) / (norm(K.e[j].nodes[2]) + norm(K.e[i].X))) + K.mat_loc[1,1][i]
+    else
+        return K.μ / (2 * π) * (ellipe(m(K.e[i].X, K.e[j].nodes[1])) / (norm(K.e[j].nodes[1]) - norm(K.e[i].X)) - ellipe(m(K.e[i].X, K.e[j].nodes[2])) / (norm(K.e[j].nodes[2]) - norm(K.e[i].X)) + ellipk(m(K.e[i].X, K.e[j].nodes[1])) / (norm(K.e[j].nodes[1]) + norm(K.e[i].X)) - ellipk(m(K.e[i].X, K.e[j].nodes[2])) / (norm(K.e[j].nodes[2]) + norm(K.e[i].X)))
+    end
+end
+
+function Base.size(K::DD3DAxisymmetricJacobianMatrix)
+    return K.n, K.n
+end
+
 """
 Piecewise constant (PWC) three-dimensional normal jacobian matrix
 """
@@ -233,29 +271,29 @@
 end
 
 """
-Piecewise constant (PWC) three-dimensional shear axis symmetric jacobian matrix
+Piecewise constant (PWC) three-dimensional shear axisymmetric jacobian matrix
 """
-struct DD3DShearAxisSymmetricJacobianMatrix{T<:Real} <: ElasticKernelMatrix{T}
+struct DD3DShearAxisymmetricJacobianMatrix{T<:Real} <: ElasticKernelMatrix{T}
     e::Vector{DDTriangleElem{T}}
     mat_loc::Matrix{Vector{T}}
     μ::T
     n::Int
 
     # Constructor
-    function DD3DShearAxisSymmetricJacobianMatrix(mesh::DDMesh2D{T}, mat_loc::Matrix{Vector{T}}, μ::T) where {T<:Real}
+    function DD3DShearAxisymmetricJacobianMatrix(mesh::DDMesh2D{T}, mat_loc::Matrix{Vector{T}}, μ::T) where {T<:Real}
         return new{T}(mesh.elems, mat_loc, μ, length(mesh.elems))
     end
 end
 
-function Base.getindex(K::DD3DShearAxisSymmetricJacobianMatrix, i::Int, j::Int)
+function Base.getindex(K::DD3DShearAxisymmetricJacobianMatrix, i::Int, j::Int)
     if (i == j)
         return K.μ / (4 * π) * integralI003(K.e[i], K.e[j]) + K.mat_loc[1,1][i]
     else
         return K.μ / (4 * π) * integralI003(K.e[i], K.e[j])
     end
 end
 
-function Base.size(K::DD3DShearAxisSymmetricJacobianMatrix)
+function Base.size(K::DD3DShearAxisymmetricJacobianMatrix)
     return K.n, K.n
 end
 

src/constraints/friction.jl

@@ -75,7 +75,7 @@
             return Δp
         end
     end
-    throw(ErrorException("Plastic update failed after $(iter) iterations!"))
+    throw(ErrorException("Plastic update failed after $(cst.max_iter) iterations!"))
 end
 
 # Default traction calculations - 2D

src/outputs.jl

@@ -69,7 +69,7 @@ struct VTKDomainOutput{T<:Real} <: AbstractOutput
 end
 
 function createVTK(out::VTKDomainOutput{T}, problem::AbstractDDProblem{T}, time::T, dt::T, time_step::Int) where {T<:Real}
-    vtk = vtk_grid(string(out.filename_base, "_", string(time_step)), out.vtk_points, out.vtk_cells; append=true, ascii=false, compress=true)
+    vtk = vtk_grid(string(out.filename_base, "_", string(time_step)), out.vtk_points, out.vtk_cells; append=false, ascii=false, compress=true)
     if hasproperty(problem, :w)
         vtk["w", VTKCellData()] = problem.w.value
         if (dt != 0.0)

src/problem.jl

@@ -19,6 +19,9 @@
     " A boolean to specify if the problem is transient"
     transient::Bool
 
+    " Axisymmetric geometry"
+    axisymmetric::Bool
+
     " The normal DD variable"
     w::Variable{T}
 
@@ -37,9 +40,16 @@
     " Pressure coupling"
     fluid_coupling::AbstractFluidCoupling
 
+    " Coordinate system"
+    coord_type::String
+
     " Constructor"
-    function NormalDDProblem(mesh::DDMesh{T}; transient::Bool=false, μ::T=1.0, ν::T=0.0) where {T<:Real}
-        return new{T}(mesh, μ, ν, length(mesh.elems), length(mesh.elems), transient,
+    function NormalDDProblem(mesh::DDMesh{T}; transient::Bool=false, axisymmetric::Bool=false, μ::T=1.0, ν::T=0.0) where {T<:Real}
+        # Check Poisson ratio and axisymmetric geometry
+        if ((ν != 0.0) && (axisymmetric))
+            warn("The Poisson's ratio is set ν ≠ 0 with an axisymmetric geometry. The axisymmetric geometry will be ignored!")
+        end
+        return new{T}(mesh, μ, ν, length(mesh.elems), length(mesh.elems), transient, axisymmetric,
             Variable(T, :w, length(mesh.elems)),
             AuxVariable(T, :σ, length(mesh.elems)),
             DefaultConstraint(),
@@ -69,6 +79,9 @@
     " A boolean to specify if the problem is transient"
     transient::Bool
 
+    " Axisymmetric geometry"
+    axisymmetric::Bool
+
     " The shear DD variables"
     δ::Variable{T}
 
@@ -91,9 +104,13 @@
     fluid_coupling::AbstractFluidCoupling
 
     " Constructor"
-    function ShearDDProblem(mesh::DDMesh{T}; transient::Bool=false, μ::T=1.0, ν::T=0.0) where {T<:Real}
+    function ShearDDProblem(mesh::DDMesh{T}; transient::Bool=false, axisymmetric::Bool=false, μ::T=1.0, ν::T=0.0) where {T<:Real}   
         if isa(mesh, DDMesh1D)
-            return new{T}(mesh, μ, ν,  length(mesh.elems), length(mesh.elems), transient,
+            # Check Poisson ratio and axisymmetric geometry
+            if (ν != 0.0)
+                warn("The Poisson's ratio is set ν ≠ 0 on a 1D mesh. Its value will be ignored!")
+            end
+            return new{T}(mesh, μ, ν,  length(mesh.elems), length(mesh.elems), transient, axisymmetric,
                 Variable(T, :δ, length(mesh.elems)),
                 AuxVariable(T, :σ, length(mesh.elems)),
                 AuxVariable(T, :τ, length(mesh.elems)),
@@ -104,7 +121,7 @@
             )
         else
             if (ν == 0.0)
-                return new{T}(mesh, μ, ν,  length(mesh.elems), length(mesh.elems), transient,
+                return new{T}(mesh, μ, ν,  length(mesh.elems), length(mesh.elems), transient, axisymmetric,
                     Variable(T, :δ, length(mesh.elems)),
                     AuxVariable(T, :σ, length(mesh.elems)),
                     AuxVariable(T, :τ, length(mesh.elems)),
@@ -114,7 +131,10 @@
                     DefaultFluidCoupling(),
                 )
             else
-                return new{T}(mesh, μ, ν,  length(mesh.elems), 2*length(mesh.elems), transient,
+                if (axisymmetric)
+                    warn("The Poisson's ratio is set ν ≠ 0 with an axisymmetric geometry. The axisymmetric geometry will be ignored!")
+                end
+                return new{T}(mesh, μ, ν,  length(mesh.elems), 2*length(mesh.elems), transient, axisymmetric,
                     Variable(T, :δ, 2*length(mesh.elems)),
                     AuxVariable(T, :σ, length(mesh.elems)),
                     AuxVariable(T, :τ, 2*length(mesh.elems)),

src/solver.jl

@@ -82,8 +82,13 @@ mutable struct DDSolver{R,T<:Real}
 
         # Assemble collocation matrix
         if isa(problem.mesh, DDMesh1D)
-            K = DD2DElasticMatrix(problem.mesh, problem.μ)
-            Kj = DD2DJacobianMatrix(problem.mesh, mat_loc, problem.μ)
+            if (problem.axisymmetric)
+                K = DD3DAxisymmetricElasticMatrix(problem.mesh, problem.μ)
+                Kj = DD3DAxisymmetricJacobianMatrix(problem.mesh, mat_loc, problem.μ)
+            else
+                K = DD2DElasticMatrix(problem.mesh, problem.μ)
+                Kj = DD2DJacobianMatrix(problem.mesh, mat_loc, problem.μ)
+            end
         elseif isa(problem.mesh, DDMesh2D)
             K = DD3DNormalElasticMatrix(problem.mesh, problem.μ, problem.ν)
             Kj = DD3DNormalJacobianMatrix(problem.mesh, mat_loc, problem.μ, problem.ν)
@@ -95,7 +100,7 @@ mutable struct DDSolver{R,T<:Real}
         # Compatibility
         comp = PartialACA(; atol=hmat_atol)
         # Assemble H-matrix
-        E = assemble_hmatrix(K, Xclt, Xclt; adm, comp, global_index=true, threads=true, distributed=false)
+        E = assemble_hmatrix(K, Xclt, Xclt; adm, comp, global_index=true, threads=false, distributed=false)
         mat = E
 
         # Preconditioner
@@ -134,14 +139,19 @@ mutable struct DDSolver{R,T<:Real}
 
         # Assemble collocation matrix
         if isa(problem.mesh, DDMesh1D)
-            K = DD2DElasticMatrix(problem.mesh, problem.μ)
-            Kj = DD2DJacobianMatrix(problem.mesh, mat_loc, problem.μ)
+            if (problem.axisymmetric)
+                K = DD3DAxisymmetricElasticMatrix(problem.mesh, problem.μ)
+                Kj = DD3DAxisymmetricJacobianMatrix(problem.mesh, mat_loc, problem.μ)
+            else
+                K = DD2DElasticMatrix(problem.mesh, problem.μ)
+                Kj = DD2DJacobianMatrix(problem.mesh, mat_loc, problem.μ)
+            end
             # Cluster tree
             Xclt = ClusterTree([problem.mesh.elems[i].X for i in eachindex(problem.mesh.elems)])
         else
             if (problem.ν == 0.0)
-                K = DD3DShearAxisSymmetricElasticMatrix(problem.mesh, problem.μ)
-                Kj = DD3DShearAxisSymmetricJacobianMatrix(problem.mesh, mat_loc, problem.μ)
+                K = DD3DShearAxisymmetricElasticMatrix(problem.mesh, problem.μ)
+                Kj = DD3DShearAxisymmetricJacobianMatrix(problem.mesh, mat_loc, problem.μ)
                 Xclt = ClusterTree([problem.mesh.elems[i].X for i in eachindex(problem.mesh.elems)])
             else
                 K = DD3DShearElasticMatrix(problem.mesh, problem.μ, problem.ν)