Cleanup

Project.toml

@@ -11,10 +11,7 @@ Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
 GridapSolvers = "6d3209ee-5e3c-4db7-a716-942eb12ed534"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
-JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
@@ -35,13 +32,10 @@ Gridap = "0.19.8"
 GridapDistributed = "0.4"
 GridapGmsh = "0.7"
 GridapPETSc = "0.5"
-GridapSolvers = "0.6"
+GridapSolvers = "0.6, 0.7"
 IterativeSolvers = "0.9"
-JSON = "1.4"
 LinearAlgebra = "1.10"
-PartitionedArrays = "0.3"
-Roots = "2.2"
-SpecialFunctions = "2.7"
+PartitionedArrays = "≥0.3"
 StaticArrays = "1.9"
 UnPack = "1.0"
 WriteVTK = "1.21"

src/PhysicalModels/PhysicalModels.jl

@@ -6,11 +6,12 @@ using Gridap.Helpers
 using UnPack
 using ForwardDiff
 using LinearAlgebra
+using StaticArrays
+
 using ..TensorAlgebra
 using ..TensorAlgebra: _∂H∂F_2D
 using ..TensorAlgebra: trAA
-using StaticArrays
-using SpecialFunctions  # implements erf
+using ..TensorAlgebra: erf
 
 import Base: +
 import Gridap: update_state!

src/TensorAlgebra/Functions.jl

@@ -22,6 +22,28 @@ function logreg(J; Threshold=0.01)
 end
 
 
+"""
+Fast and dependency-free implementation of erf function, up to 1e-6 precision.
+"""
+@inline function erf(x::Float64)
+  p  = 0.3275911
+  a1 = 0.254829592
+  a2 = -0.284496736
+  a3 = 1.421413741
+  a4 = -1.453152027
+  a5 = 1.061405429
+
+  ax = abs(x)
+  t = 1.0 / (1.0 + p*ax)
+
+  y = (((((a5*t + a4)*t + a3)*t + a2)*t + a1)*t)
+
+  r = 1.0 - y*exp(-ax*ax)
+
+  return copysign(r, x)
+end
+
+
 function _∂H∂F_2D()
   TensorValue(0.0, 0.0, 0.0, 1.0, 0.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0)
 end

test/Project.toml

@@ -3,6 +3,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 GridapSolvers = "6d3209ee-5e3c-4db7-a716-942eb12ed534"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"

test/TestTensorAlgebra/TensorAlgebraTests.jl

@@ -1,6 +1,7 @@
 using Gridap.TensorValues
 using Gridap.Arrays
 using HyperFEM.TensorAlgebra
+using SpecialFunctions
 using Test
 
 
@@ -114,3 +115,10 @@ end
   cofA = det(A) * inv(A')
   @test isapprox(cof(A), cofA)
 end
+
+
+@testset "erf" begin
+  for x ∈ 0:.17234:2
+    @test isapprox(TensorAlgebra.erf(x), SpecialFunctions.erf(x), atol=1e-6)
+  end
+end

test/TestWeakForms/WeakForms.jl

@@ -129,76 +129,3 @@ end
 
 
 end
-
-# @testset "Assembly Jacobian Mechanics" begin
-
-#     pname = "test"
-#     meshfile = "test_cubicsphere.msh"
-#     simdir = datadir("sims", pname)
-#     setupfolder(simdir)
-
-#     geomodel = GmshDiscreteModel(datadir("models", meshfile))
-
-#     # Constitutive models
-#     ∇umacro = TensorValue(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0) * 1e-2
-
-#     Kmodel = EvolutiveKinematics(Mechano; F=(t) -> ((∇u) -> ∇u + one(∇u) + t * ∇umacro))
-#     physmodel_matrix = MooneyRivlin3D(λ=10.0, μ1=1.0, μ2=1.0, Kinematic=Kmodel)
-#     physmodel_inclussion = MooneyRivlin3D(λ=50.0, μ1=5.0, μ2=5.0, Kinematic=Kmodel)
-
-#     # Setup integration
-#     order = 1
-#     degree = 2 * order
-#     Ω = Triangulation(geomodel)
-#     dΩ = Measure(Ω, degree)
-#     Ω₁ = Interior(geomodel, tags="Phase0")
-#     dΩ₁ = Measure(Ω₁, degree)
-#     Ω₂ = Interior(geomodel, tags="Inclusions")
-#     dΩ₂ = Measure(Ω₂, degree)
-
-#     # Dirichlet conditions 
-#     evolu(Λ) = 1.0
-#     dir_u_tags = ["Corners"]
-#     dir_u_values = [[0.0, 0.0, 0.0]]
-#     dir_u_timesteps = [evolu]
-#     D_bc = DirichletBC(dir_u_tags, dir_u_values, dir_u_timesteps)
-
-#     #  FE spaces
-#     reffeu = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
-#     V = TestFESpace(geomodel, reffeu, D_bc, :H1)
-#     U = TrialFESpace(V, D_bc, 1.0)
-
-#     #  residual and jacobian function of load factor
-#     res(Λ) = (u, v) -> residual(physmodel_matrix, u, v, dΩ₁, Λ) + residual(physmodel_inclussion, u, v, dΩ₂, Λ)
-#     jac(Λ) = (u, du, v) -> jacobian(physmodel_matrix, u, du, v, dΩ₁, Λ) + jacobian(physmodel_inclussion, u, du, v, dΩ₂, Λ)
-
-#     uh = FEFunction(V, zeros(Float64, num_free_dofs(V)))
-
-#     function jach(uh, φh)
-#         jac((du, dφ), (v, vφ)) = jacobian(modelelectro, (uh, φh), (du, dφ), (v, vφ), dΩ)
-#     end
-#     function jac_mech(uh, φh)
-#         jac(du, v) = jacobian(modelelectro, Mechano, (uh, φh), du, v, dΩ)
-#     end
-#     function jac_elech(uh, φh)
-#         jac(dφ, vφ) = jacobian(modelelectro, Electro, (uh, φh), dφ, vφ, dΩ)
-#     end
-
-#     jac_ = assemble_matrix(jach(uh, φh), V, V)
-#     jac_m = assemble_matrix(jac_mech(uh, φh), Vu, Vu)
-#     jac_e = assemble_matrix(jac_elech(uh, φh), Vφ, Vφ)
-#     ≈
-#     @test norm(jac_) ≈ 18.934585248125135
-#     @test jac_[1] ≈ 0.7777777777777775
-#     @test jac_[end] ≈ -1.3333333333333326
-
-#     @test norm(jac_m) ≈ 16.420402846143244
-#     @test jac_m[1] ≈ 0.7777777777777775
-#     @test jac_m[end] ≈ 2.1111111111111094
-
-#     @test norm(jac_e) ≈ 4.216370213557837
-#     @test jac_e[1] ≈ -1.3333333333333321
-#     @test jac_e[end] ≈ -1.3333333333333326
-
-
-# end