Bump fec and cm

Project.toml

@@ -26,10 +26,10 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 [compat]
 Aqua = "0.8"
 Arpack = "0.5"
-ConstitutiveModels = "0.3"
+ConstitutiveModels = "0.4"
 Distributions = "0.25"
 DocStringExtensions = "0.9"
-FiniteElementContainers = "0.14"
+FiniteElementContainers = "0.15"
 ForwardDiff = "1"
 KernelAbstractions = "0.9"
 Krylov = "0.10"

examples/mechanics/hole_array/script_quasistatic.jl

@@ -8,23 +8,35 @@ mesh_file = Base.source_dir() * "/mesh/hole_array.exo"
 output_file = splitext(mesh_file)[1] * "-output.exo"
 
 # Times
-times = TimeStepper(0., 1., 20)
+t_end = 100.0
+times = TimeStepper(0., 2. * t_end, 40)
 
 # Physics
 physics = (;
-    Block1 = SolidMechanics(PlaneStrain(), NeoHookean())
+    # Block1 = SolidMechanics(PlaneStrain(), Hyperelastic(NeoHookean()))
+    Block1 = SolidMechanics(PlaneStrain(), FeFv())
+    # Block1 = SolidMechanics(PlaneStrain(), Hyperelastic(Hencky()))
 )
 props = (;
     Block1 = Dict{String, Any}(
         "density"         => 1.,
         "Young's modulus" => 1.,
         "Poisson's ratio" => 0.48,
         # "Jm"              => 3.
+        "G neq"           => 50.0,
+        "relaxation time" => 50.0
     )
 )
 
 # Boundary Conditions
-func_1(x, t) = -7.5 * t
+# func_1(x, t) = -9. * (t / t_end)
+function func_1(x, t)
+    if t >= t_end
+        return -9.
+    else
+        return -9. * (t / t_end)
+    end
+end
 func_2(x, t) = 0.0
 
 dirichlet_bcs = [

src/Physics.jl

@@ -25,11 +25,11 @@ end
     # constitutive
     θ = 0.0 # TODO
     ψ_q = ConstitutiveModels.helmholtz_free_energy(
-        physics.constitutive_model, props_el, dt, ∇u_q, θ, state_old_q, state_new_q
+        physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ
     )  
 
     σ_q = ConstitutiveModels.cauchy_stress(
-        physics.constitutive_model, props_el, dt, ∇u_q, θ, state_old_q, state_new_q
+        physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ
     )
     σ_q = σ_q |> symmetric
     return StandardMaterialOutput(ψ_q, σ_q, ∇u_q)
@@ -157,28 +157,43 @@ function setup_function(fspace::FunctionSpace, ::SolidMechanics)
     return VectorFunction(fspace, "displ")
 end
 
-function characteristic_element_length(
-    physics::SolidMechanics, interps, x_el,
-    t, dt, u_el, u_el_old,
-    state_old_q, state_new_q, props_el,
-)
+@inline function characteristic_element_length(x_el::SMatrix{2, N, T, N2}, u_el::SMatrix{2, N, T, N2}) where {N, T, N2}
+    x_cur = x_el + u_el
+    ndim = 2
+    nnpe = length(x_cur) ÷ ndim
+    cx = cy = zero(T)
+    for i in 1:nnpe
+        cx += x_cur[(i - 1) * ndim + 1]
+        cy += x_cur[(i - 1) * ndim + 2]
+    end
+    cx /= nnpe; cy /= nnpe
+    total = zero(T)
+    for i in 1:nnpe
+        dx = x_cur[(i - 1) * ndim + 1] - cx
+        dy = x_cur[(i - 1) * ndim + 2] - cy
+        total += sqrt(dx * dx + dy * dy)
+    end
+    return 2 * total / nnpe
+end
+
+@inline function characteristic_element_length(x_el::SMatrix{3, N, T, N3}, u_el::SMatrix{3, N, T, N3}) where {N, T, N3}
     x_cur = x_el + u_el
     ndim = 3
     nnpe = length(x_cur) ÷ ndim
-    T = eltype(x_cur)
+    # T = eltype(x_cur)
     cx = cy = cz = zero(T)
     for i in 1:nnpe
-        cx += x_cur[(i-1)*ndim + 1]
-        cy += x_cur[(i-1)*ndim + 2]
-        cz += x_cur[(i-1)*ndim + 3]
+        cx += x_cur[(i - 1) * ndim + 1]
+        cy += x_cur[(i - 1) * ndim + 2]
+        cz += x_cur[(i - 1) * ndim + 3]
     end
     cx /= nnpe; cy /= nnpe; cz /= nnpe
     total = zero(T)
     for i in 1:nnpe
-        dx = x_cur[(i-1)*ndim + 1] - cx
-        dy = x_cur[(i-1)*ndim + 2] - cy
-        dz = x_cur[(i-1)*ndim + 3] - cz
-        total += sqrt(dx*dx + dy*dy + dz*dz)
+        dx = x_cur[(i - 1) * ndim + 1] - cx
+        dy = x_cur[(i - 1) * ndim + 2] - cy
+        dz = x_cur[(i - 1) * ndim + 3] - cz
+        total += sqrt(dx * dx + dy * dy + dz * dz)
     end
     return 2 * total / nnpe
 end
@@ -193,12 +208,12 @@ end
     # kinematics
     ∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
 
-    mat_props = @views props_el[2:end]
+    # mat_props = @views props_el[2:end]
 
     # constitutive
     θ = 0.0 # TODO
     ψ_q = ConstitutiveModels.helmholtz_free_energy(
-        physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q
+        physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ
     )  
     return JxW * ψ_q
 end
@@ -215,12 +230,10 @@ end
     # kinematics
     ∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
 
-    mat_props = @views props_el[2:end]
-
     # constitutive
     θ = 0.0 # TODO
     ψ_q = func(
-        physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q
+        physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ
     )  
     return ψ_q
 end
@@ -236,12 +249,10 @@ end
     # kinematics
     ∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
 
-    mat_props = @views props_el[2:end]
-
     # constitutive
     θ = 0.0 # TODO
     ψ_q = func(
-        physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q
+        physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ
     )  
     return JxW * ψ_q
 end
@@ -289,19 +300,9 @@ end
 )
     interps = map_interpolants(interps, x_el)
     (; X_q, N, ∇N_X, JxW) = interps
-    # ∇u_q = interpolate_field_gradients(physics, interps, u_el)
-
-    # kinematics
-    # ∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
 
-    # mat_props = @views props_el[2:end]
     ρ = props_el[1]
 
-    # constitutive
-    # θ = 0. # TODO
-    # A_q = ConstitutiveModels.material_tangent(
-    #     physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q
-    # )
     scatter_with_values_and_values!(
       storage, physics.formulation, e, conn, N, JxW * ρ
     )
@@ -346,16 +347,13 @@ end
     # kinematics
     ∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
 
-    mat_props = @views props_el[2:end]
-
     # constitutive
     θ = 0.0 # TODO
     P_q = ConstitutiveModels.pk1_stress(
-        physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q
+        physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ
     )    
     # turn into voigt notation
     P_q = extract_stress(physics.formulation, P_q)
-    # P_q = tovoigt(SVector, P_q)
     G_q = discrete_gradient(physics.formulation, ∇N_X)
     f_q = G_q * P_q
     return JxW * f_q[:]
@@ -374,17 +372,30 @@ end
     # kinematics
     ∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
 
-    mat_props = @views props_el[2:end]
-
     # constitutive
     θ = 0.0 # TODO
     P_q = ConstitutiveModels.pk1_stress(
-        physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q
+        physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ
     )
     scatter_with_gradients!(storage, physics.formulation, e, conn, ∇N_X, JxW * P_q)
     return nothing
 end
 
+@inline function stable_time_step(
+    physics::SolidMechanics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)   
+    N, D = size(interps.∇N_ξ)
+    x_el = SMatrix{D, N, eltype(x_el), N * D}(x_el.data)
+    u_el = SMatrix{D, N, eltype(u_el), N * D}(u_el.data)
+    l = characteristic_element_length(x_el, u_el)
+    ρ = props_el[1]
+    M = p_wave_modulus(physics.constitutive_model, props_el)
+    c_p = sqrt(M / ρ)
+    CFL = 0.1
+    stable_dt = CFL * l / c_p
+    return stable_dt
+end
+
 @inline function FiniteElementContainers.stiffness(  
     physics::SolidMechanics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
 )
@@ -395,12 +406,10 @@ end
     # kinematics
     ∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
 
-    mat_props = @views props_el[2:end]
-
     # constitutive
     θ = 0. # TODO
     A_q = ConstitutiveModels.material_tangent(
-        physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q
+        physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ
     )
     # turn into voigt notation
     K_q = extract_stiffness(physics.formulation, A_q)
@@ -421,12 +430,10 @@ end
     # kinematics
     ∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
 
-    mat_props = @views props_el[2:end]
-
     # constitutive
     θ = 0. # TODO
     A_q = ConstitutiveModels.material_tangent(
-        physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q
+        physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ
     )
     scatter_with_gradients_and_gradients!(
       storage, physics.formulation, e, conn, ∇N_X, JxW * A_q
@@ -442,12 +449,10 @@ end
     ∇u_q = interpolate_field_gradients(physics, interps, u_el)
     ∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
 
-    mat_props = @views props_el[2:end]
-
     # constitutive
     θ = 0. # TODO
     A_q = ConstitutiveModels.material_tangent(
-        physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q
+        physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ
     )
     # turn into voigt notation
     K_q = extract_stiffness(physics.formulation, A_q)
@@ -466,12 +471,10 @@ end
     ∇u_q = interpolate_field_gradients(physics, interps, u_el)
     ∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
 
-    mat_props = @views props_el[2:end]
-
     # constitutive
     θ = 0. # TODO
     A_q = ConstitutiveModels.material_tangent(
-        physics.constitutive_model, mat_props, dt, ∇u_q, θ, state_old_q, state_new_q
+        physics.constitutive_model, props_el, state_old_q, state_new_q, dt, ∇u_q, θ
     )
     scatter_with_gradients_and_gradients!(storage, physics.formulation, e, conn, ∇N_X, JxW * A_q, v_el)
     return nothing

src/objectives/ExplicitDynamicsObjective.jl

@@ -171,41 +171,17 @@ end
 # copied from Carina.jl
 function _compute_stable_dt(asm, p, CFL)
     fspace = FEC.function_space(asm.dof)
-
-    # Pre-allocate per-block storage for element char lengths (nq × nelem)
-    # char_len_storage = Array{Float64, 3}[]
-    # for (b, ref_fe) in enumerate(fspace.ref_fes)
-    # for b in 1:FEC.num_blocks(fspace)
-    #     ref_fe = FEC.block_reference_element(fspace, b)
-    #     nquad = num_cell_quadrature_points(ref_fe)
-    #     nelem = num_elements(fspace, b)
-    #     push!(char_len_storage, zeros(Float64, 1, nquad, nelem))
-    # end
-    # char_len_storage = NamedTuple{keys(fspace.ref_fes)}(char_len_storage)
-    char_len_storage = L2Field(undef, Float64, 1, FEC.block_quadrature_sizes(fspace))
-    fill!(char_len_storage, 0.0)
+    stable_dts = L2Field(undef, Float64, 1, FEC.block_quadrature_sizes(fspace))
+    fill!(stable_dts, 0.0)
 
     # Assemble per-element char lengths on device
     U_zeros = zeros(Float64, length(asm.dof.unknown_dofs))
     FEC.assemble_quadrature_quantity!(
-        char_len_storage, nothing, asm.dof,
-        characteristic_element_length,
+        stable_dts, nothing, asm.dof,
+        stable_time_step,
         U_zeros, p
     )
-
-    # Min-reduction over all blocks
-    stable_dt = Inf
-    for (b, (block_physics, props)) in enumerate(zip(
-        values(p.physics), values(p.properties),
-    ))
-        block_storage = FEC.block_view(char_len_storage, b)
-        ρ   = props[1]
-        M   = p_wave_modulus(block_physics.constitutive_model, props)
-        c_p = sqrt(M / ρ)
-        h_min = minimum(block_storage)   # GPU-native reduction if on device
-        block_dt = CFL * h_min / c_p
-        stable_dt = min(stable_dt, block_dt)
-    end
+    stable_dt = reduce(min, stable_dts)
     return stable_dt
 end
 

test/TestExplicitDynamics.jl

@@ -8,12 +8,12 @@ function test_explicit_dynamics()
 
     # Physics
     physics = (;
-        cube = SolidMechanics(
-            ThreeDimensional(), NeoHookean()
+        sphere = SolidMechanics(
+            ThreeDimensional(), Hyperelastic(NeoHookean())
         )
     )
     props = (;
-        cube = Dict{String, Any}(
+        sphere = Dict{String, Any}(
             "density"         => 1.e3,
             "Young's modulus" => 1.e4,
             "Poisson's ratio" => 0.33

test/TestHeatConduction.jl

@@ -5,10 +5,10 @@ function test_heat_conduction()
     times = TimeStepper(0.0, 75.0, 100)
 
     physics = (;
-        unnamed_block_1 = Cthonios.HeatConduction()
+        var"" = Cthonios.HeatConduction()
     )
     props = (;
-        unnamed_block_1 = Dict{String, Any}(
+        var"" = Dict{String, Any}(
             "density"              => 1.0,
             "heat capacity"        => 20.0,
             "thermal conductivity" => 1.0

test/TestImplicitDynamics.jl

@@ -9,12 +9,12 @@ function test_implicit_dynamics()
 
     # Physics
     physics = (;
-        cube = SolidMechanics(
-            ThreeDimensional(), NeoHookean()
+        sphere = SolidMechanics(
+            ThreeDimensional(), Hyperelastic(NeoHookean())
         )
     )
     props = (;
-        cube = Dict{String, Any}(
+        sphere = Dict{String, Any}(
             "density"         => 1.e3,
             "Young's modulus" => 1.e4,
             "Poisson's ratio" => 0.33

test/TestQuasiStatic.jl

@@ -8,7 +8,7 @@ function test_quasistatic()
 
     # Physics
     physics = (;
-        Block1 = SolidMechanics(PlaneStrain(), NeoHookean())
+        Block1 = SolidMechanics(PlaneStrain(), Hyperelastic(NeoHookean()))
     )
     props = (;
         Block1 = Dict{String, Any}(