Merge pull request #14 from MultiSimOLab/latest-gridap

Latest Gridap.jl

Author: jmartfrut

src/PhysicalModels/ViscousModels.jl

@@ -157,7 +157,7 @@ function JacobianReturnMapping(γα, Ce, Se, Se_trial, ∂Se∂Ce, F, λα)
     ∂res1_∂λα = -(1-γα) * Ge
     ∂res2_∂Ce = Ge
     res = [get_array(res1)[:]; res2]
-    ∂res = MMatrix{10,10}(zeros(10, 10))
+    ∂res = MMatrix{10,10}(zeros(10, 10))  # TODO: It'd be nice to use hvcat: ∂res = [∂res1_Ce ∂res1_∂λα; ∂res2_∂Ce 0.0]
     ∂res[1:9, 1:9] = get_array(∂res1_∂Ce)
     ∂res[1:9, 10] = get_array(∂res1_∂λα)[:]
     ∂res[10, 1:9] = (get_array(∂res2_∂Ce)[:])'
@@ -309,10 +309,9 @@ end
 
 function Energy(obj::ViscousIncompressible, Δt::Float64,
                 Ψe::Function, Se_::Function, ∂Se∂Ce_::Function,
-                F::TensorValue, Fn::TensorValue, stateVars::VectorValue)
-  state_vars = get_array(stateVars)
-  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))  # TODO: Update tensor slicing until next Gridap version has been released
-  λαn = state_vars[10]
+                F::TensorValue, Fn::TensorValue, A::VectorValue)
+  Uvn = TensorValue{3,3}(A[1:9]...)
+  λαn = A[10]
   #------------------------------------------
   # Get kinematics
   #------------------------------------------
@@ -343,17 +342,16 @@ end
 - `∂Se∂Ce_`: 2nd Piola Derivatives (function of C)
 - `F`: Current deformation gradient
 - `Fn`: Previous deformation gradient
-- `stateVars`: State variables (Uvα and λα)
+- `A`: State variables (Uvα and λα)
 
 # Return
 - `Pα::Gridap.TensorValues.TensorValue`
 """
 function Piola(obj::ViscousIncompressible, Δt::Float64,
                 Se_::Function, ∂Se∂Ce_::Function,
-                F::TensorValue, Fn::TensorValue, stateVars::VectorValue)
-  state_vars = get_array(stateVars)
-  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))  # TODO: Update tensor slicing until next Gridap version has been released
-  λαn = state_vars[10]
+                F::TensorValue, Fn::TensorValue, A::VectorValue)
+  Uvn = TensorValue{3,3}(A[1:9]...)
+  λαn = A[10]
   #------------------------------------------
   # Get kinematics
   #------------------------------------------
@@ -386,17 +384,16 @@ Visco-Elastic model for incompressible case
 - `∂Se∂Ce_`: 2nd Piola Derivatives (function of C)
 - `∇u_`: Current deformation gradient
 - `∇un_`: Previous deformation gradient
-- `stateVars`: State variables (Uvα and λα)
+- `A`: State variables (Uvα and λα)
 
 # Return
 - `Cα::Gridap.TensorValues.TensorValue`
 """
 function Tangent(obj::ViscousIncompressible, Δt::Float64,
                  Se_::Function, ∂Se∂Ce_::Function,
-                 F::TensorValue, Fn::TensorValue, stateVars::VectorValue)
-  state_vars = get_array(stateVars)
-  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))  # TODO: Update tensor slicing until next Gridap version has been released
-  λαn = state_vars[10]
+                 F::TensorValue, Fn::TensorValue, A::VectorValue)
+  Uvn = TensorValue{3,3}(A[1:9]...)
+  λαn = A[10]
   #------------------------------------------
   # Get kinematics
   #------------------------------------------
@@ -432,18 +429,17 @@ end
     - `∂Se∂Ce_::Function`: Piola Derivatives (function of C)
     - `∇u_::TensorValue`
     - `∇un_::TensorValue`
-    - `stateVars::VectorValue`: State variables (10-component vector gathering Uvα and λα)
+    - `A::VectorValue`: State variables (10-component vector gathering Uvα and λα)
 
     # Return
     - `::bool`: indicates whether the state variables should be updated
     - `::VectorValue`: State variables at new time
 """
 function ReturnMapping(obj::ViscousIncompressible, Δt::Float64,
                        Se_::Function, ∂Se∂Ce_::Function,
-                       F::TensorValue, Fn::TensorValue, stateVars::VectorValue)
-  state_vars = get_array(stateVars)
-  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))  # TODO: Update tensor slicing until next Gridap version has been released
-  λαn = state_vars[10]
+                       F::TensorValue, Fn::TensorValue, A::VectorValue)
+  Uvn = TensorValue{3,3}(A[1:9]...)
+  λαn = A[10]
   #------------------------------------------
   # Get kinematics
   #------------------------------------------
@@ -468,10 +464,9 @@ end
 
 function ViscousDissipation(obj::ViscousIncompressible, Δt::Float64,
                        Se_::Function, ∂Se∂Ce_::Function,
-                       F::TensorValue, Fn::TensorValue, stateVars::VectorValue)
-  state_vars = get_array(stateVars)
-  Uvn = TensorValue{3,3}(Tuple(state_vars[1:9]))  # TODO: Update tensor slicing until next Gridap version has been released
-  λαn = state_vars[10]
+                       F::TensorValue, Fn::TensorValue, A::VectorValue)
+  Uvn = TensorValue{3,3}(A[1:9]...)
+  λαn = A[10]
   #------------------------------------------
   # Get kinematics
   #------------------------------------------

src/TensorAlgebra/TensorAlgebra.jl

@@ -53,7 +53,7 @@ export Ellipsoid
 Compute the square root of a 3x3 matrix by means of eigen decomposition.
 """
 function sqrt(A::TensorValue{3})
-  λ, Q = eigen(get_array(A))  # TODO: the get_array must be removed as long as it is supported after https://github.com/gridap/Gridap.jl/pull/1157
+  λ, Q = eigen(A)
   λ = sqrt.(λ)
   TensorValue{3}(λ[1]*Q[1:3]*Q[1:3]' + λ[2]*Q[4:6]*Q[4:6]' + λ[3]*Q[7:9]*Q[7:9]')
 end

test/TestConstitutiveModels/ViscousModelsBenchmarks.jl

@@ -0,0 +1,28 @@
+using Gridap.TensorValues
+using HyperFEM.PhysicalModels
+using BenchmarkTools
+
+
+function benchmark_viscous_model(model)
+  Ψ, ∂Ψu, ∂Ψuu = model(Δt = 1e-2)
+  F = TensorValue(1.:9...) * 1e-3 + I3
+  Fn = TensorValue(1.:9...) * 5e-4 + I3
+  Uvn = TensorValue(1.,2.,3.,2.,4.,5.,3.,5.,6.) * 2e-4 + I3
+  J = det(F)
+  Uvn *= J^(-1/3)
+  λvn = 1e-3
+  Avn = VectorValue(Uvn.data..., λvn)
+  print("Ψ(F, Fn, Avn)    |")
+  @btime $Ψ($F, $Fn, $Avn)
+  print("∂Ψu(F, Fn, Avn)  |")
+  @btime $∂Ψu($F, $Fn, $Avn)
+  print("∂Ψuu(F, Fn, Avn) |")
+  @btime $∂Ψuu($F, $Fn, $Avn)
+end
+
+
+elasto = NeoHookean3D(λ=λ, μ=μ)
+visco = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μ1), τ1)
+visco_elastic = GeneralizedMaxwell(elasto, visco)
+benchmark_viscous_model(visco);
+benchmark_viscous_model(visco_elastic);

test/TestConstitutiveModels/ViscousModelsTests.jl

@@ -5,8 +5,6 @@ using HyperFEM.TensorAlgebra
 using HyperFEM.PhysicalModels
 using StaticArrays
 using Test
-import Base:≈
-≈(a::MultiValue, b::MultiValue; kwargs...) = ≈(get_array(a), get_array(b); kwargs...)
 
 
 μ      = 1.367e4  # Pa

test/TestTensorAlgebra/TensorAlgebra.jl

@@ -45,7 +45,6 @@ end
 # @code_warntype (A ⊗₁₃²⁴ B)
 # @code_warntype (D × A)
 
- 
 
 @testset "cross" begin
   A = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
@@ -60,11 +59,7 @@ end
   @test norm(B × C) == 1.104276381618298e-6
   @test norm(D × A) == 0.00012378691368638284
   @test norm(get_array(A) × get_array(B))== 6.246230863488799e-5
-
-
-  
-
- end
+end
 
 
 @testset "inner" begin
@@ -124,3 +119,10 @@ end
   @test contraction_IP_PJKL(A,H) == TensorValue(7.0, 10.0, 15.0, 22.0, 23.0, 34.0, 31.0, 46.0, 39.0, 58.0, 47.0, 70.0, 55.0, 82.0, 63.0, 94.0)
   @test contraction_IP_JPKL(A,H) == TensorValue(10.0, 14.0, 14.0, 20.0, 26.0, 38.0, 30.0, 44.0, 42.0, 62.0, 46.0, 68.0, 58.0, 86.0, 62.0, 92.0)
 end
+
+
+@testset "sqrt" begin
+  A = TensorValue(1.:9...)
+  A = A'*A + I3
+  @test isapprox(sqrt(A), TensorValue(sqrt(get_array(A))), rtol=1e-14)
+end

test/runtests.jl

@@ -9,10 +9,6 @@ using Gridap.TensorValues
 using HyperFEM: jacobian, IterativeSolver, solve!
 using BenchmarkTools
 
-import Base.isapprox
-
-isapprox(A::MultiValue, B::MultiValue; kwargs...) = isapprox(get_array(A), get_array(B); kwargs...)
-
 
 @testset "HyperFEMTests" verbose = true begin