diff --git a/Project.toml b/Project.toml index 6fe569e..6e11671 100644 --- a/Project.toml +++ b/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" diff --git a/src/PhysicalModels/PhysicalModels.jl b/src/PhysicalModels/PhysicalModels.jl index 6dd557a..b10756b 100644 --- a/src/PhysicalModels/PhysicalModels.jl +++ b/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! diff --git a/src/TensorAlgebra/Functions.jl b/src/TensorAlgebra/Functions.jl index 36e3f5b..fb8f2f6 100644 --- a/src/TensorAlgebra/Functions.jl +++ b/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 diff --git a/test/Project.toml b/test/Project.toml index d898967..1e56ebb 100644 --- a/test/Project.toml +++ b/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" diff --git a/test/TestTensorAlgebra/TensorAlgebraTests.jl b/test/TestTensorAlgebra/TensorAlgebraTests.jl index c84ecff..b0b191d 100644 --- a/test/TestTensorAlgebra/TensorAlgebraTests.jl +++ b/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 diff --git a/test/TestWeakForms/WeakForms.jl b/test/TestWeakForms/WeakForms.jl index 79c16ca..7709cca 100644 --- a/test/TestWeakForms/WeakForms.jl +++ b/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 \ No newline at end of file