Test cleanup

Project.toml

@@ -50,7 +50,6 @@ SparseArrays = "1"
 StaticArrays = "1"
 Tensors = "1"
 Test = "1"
-TestSetExtensions = "2"
 TimerOutputs = "0.5"
 julia = "1"
 
@@ -62,7 +61,6 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Exodus = "f57ae99e-f805-4780-bdca-96e224be1e5a"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
 
 [targets]
-test = ["AMDGPU", "Aqua", "CUDA", "Exodus", "MPI", "Test", "TestSetExtensions"]
+test = ["AMDGPU", "Aqua", "CUDA", "Exodus", "MPI", "Test"]

examples/electromechanics/script.jl

@@ -0,0 +1,119 @@
+using Exodus
+using FiniteElementContainers
+using ForwardDiff
+using StaticArrays
+using Tensors
+
+struct ElectroMechanics <: FiniteElementContainers.AbstractPhysics{4, 4, 0}
+end
+
+function FiniteElementContainers.create_properties(physics, props_dict)
+    return SVector{4, Float64}(
+        props_dict["shear modulus"],
+        props_dict["bulk modulus"],
+        props_dict["Jm"],
+        props_dict["permitivity"]
+    )
+end
+
+@inline function FiniteElementContainers.energy(
+    physics::ElectroMechanics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, props_el
+  )
+    G, K, Jm, ϵ = props_el
+
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    U_q, ∇U_q = interpolate_field_values_and_gradients(physics, interps, u_el)
+
+    ∇u_q = Tensor{2, 3, eltype(u_el), 9}(∇U_q[1:3, :])
+    F_q = ∇u_q + one(∇u_q)
+    J_q = det(F_q)
+    C_q = tdot(F_q)
+    C_inv_q = inv(C_q)
+    E_R_q = Vec{3, eltype(u_el)}(∇U_q[4, :])
+    D_R_q = ϵ * J_q * dot(C_inv_q, E_R_q)
+
+    # kinematics
+    I_1_bar = tr(J_q^(-2. / 3.) * C_q)
+
+    # constitutive
+    ψ_vol  = 0.5 * K * (0.5 * (J_q^2 - 1) - log(J_q))
+    ψ_dev  = -G * Jm / 2. * log(1. - (I_1_bar - 3.) / Jm)
+    ψ_mech = ψ_vol + ψ_dev
+    ψ_elec = (1. / (2 * ϵ * J_q)) * dot(D_R_q, C_q, D_R_q)
+    ψ_q    = ψ_mech + ψ_elec
+    state_new_q = copy(state_old_q)
+    return JxW * ψ_q, state_new_q
+end
+
+@inline function FiniteElementContainers.residual(
+    physics::ElectroMechanics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, props_el
+)
+    state_new_q = copy(state_old_q)
+    return ForwardDiff.gradient(
+        z -> energy(physics, interps, x_el, t, dt, z, u_el_old, state_old_q, props_el)[1], u_el
+    ), state_new_q
+end
+
+@inline function FiniteElementContainers.stiffness(
+    physics::ElectroMechanics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, props_el
+)
+    state_new_q = copy(state_old_q)
+    return ForwardDiff.hessian(
+        z -> energy(physics, interps, x_el, t, dt, z, u_el_old, state_old_q, props_el)[1], u_el
+    ), state_new_q
+end
+
+function run_electromechanics()
+    mesh_file = Base.source_dir() * "/cube.g"
+    output_file = Base.source_dir() * "/output.e"
+
+    zero_func(_, _) = 0.
+    potential_func(_, t) = 0.8 * t
+
+    mesh = UnstructuredMesh(mesh_file)
+    V = FunctionSpace(mesh, H1Field, Lagrange)
+    displ = VectorFunction(V, :displ)
+    φ = ScalarFunction(V, :phi)
+    u = FiniteElementContainers.GeneralFunction(displ, φ)
+    asm = SparseMatrixAssembler(u; use_condensed=true)
+    times = TimeStepper(0., 1., 11)
+
+    dbcs = DirichletBC[
+        DirichletBC("displ_x", "ssx-", zero_func)
+        DirichletBC("displ_y", "ssy-", zero_func)
+        DirichletBC("displ_z", "ssz-", zero_func)
+        DirichletBC("phi", "ssz-", zero_func)
+        DirichletBC("phi", "ssz+", potential_func)
+    ]
+
+    physics = ElectroMechanics()
+    props = Dict(
+        "shear modulus" => 1.0,
+        "bulk modulus"  => 100.0,
+        "Jm"            => 7.,
+        "permitivity"   => 1.
+    )
+    props = create_properties(physics, props)
+    p = create_parameters(mesh, asm, physics, props; dirichlet_bcs=dbcs, times=times)
+
+    # setup solver and integrator
+    # linear_solver = IterativeLinearSolver(asm, :CgSolver)
+    linear_solver = DirectLinearSolver(asm)
+    solver = NewtonSolver(linear_solver)
+    # solver = nsolver(lsolver(asm))
+    integrator = QuasiStaticIntegrator(solver)
+
+    pp = PostProcessor(mesh, output_file, u)
+    write_times(pp, 1, 0.0)
+
+    for n in 1:11
+        evolve!(integrator, p)
+        write_times(pp, n + 1, FiniteElementContainers.current_time(p.times))
+        write_field(pp, n + 1, ("displ_x", "displ_y", "displ_z", "phi"), p.h1_field)
+    end
+
+    close(pp)
+end
+
+run_electromechanics()

src/FiniteElementContainers.jl

@@ -182,6 +182,7 @@ include("PostProcessors.jl")
 #
 include("TimeSteppers.jl")
 include("Parameters.jl")
+include("integrals/Integrals.jl")
 include("solvers/Solvers.jl")
 include("integrators/Integrators.jl")
 

src/Functions.jl

@@ -178,6 +178,28 @@ function Base.show(io::IO, ::SymmetricTensorFunction{S, F}) where {S, F}
   println(io, "  names: $S")
 end
 
+struct GeneralFunction{S, F} <: AbstractFunction{S, F}
+  fspace::F
+end
+
+function GeneralFunction(args...)
+  fspace = args[1].fspace
+  for arg in args
+    @assert typeof(arg.fspace) == typeof(fspace)
+  end
+  # syms = mapreduce(names, vcat, args)
+  syms = ()
+  for arg in args
+    syms = (syms..., names(arg)...)
+  end
+  # @show syms
+  return GeneralFunction{syms, typeof(fspace)}(fspace)
+end
+
+function Base.show(io::IO, ::GeneralFunction{S, F}) where {S, F}
+  println(io, "GeneralFunction:")
+  println(io, "  names: $S")
+end
 
 # struct StateFunction{S, F, NS, NQ} <: AbstractFunction{S, F}
 #   fspace::F

src/assemblers/SparseMatrixAssembler.jl

@@ -8,15 +8,13 @@ struct SparseMatrixAssembler{
   Condensed,
   NumArrDims,
   NumFields,
-  # BlockPattern      <: NamedTuple,
   IV                <: AbstractArray{Int, 1},
   RV                <: AbstractArray{Float64, 1},
   Var               <: AbstractFunction,
   FieldStorage      <: AbstractField{Float64, NumArrDims, RV, NumFields}, # should we make 2 not hardcoded? E.g. for 
   QuadratureStorage <: NamedTuple
 } <: AbstractAssembler{DofManager{Condensed, Int, IV, Var}}
   dof::DofManager{Condensed, Int, IV, Var}
-  # matrix_pattern::SparseMatrixPattern{IV, BlockPattern, RV}
   matrix_pattern::SparseMatrixPattern{IV, RV}
   vector_pattern::SparseVectorPattern{IV}
   constraint_storage::RV

src/assemblers/SparseMatrixAssemblerNew.jl

@@ -0,0 +1,11 @@
+struct SparseMatrixAssembler{
+    Condensed,
+    IV  <: AbstractArray{Int, 1},
+    RV  <: AbstractArray{Float64, 1},
+    Var <: AbstractFunction
+} <: AbstractAssembler{DofManager{Condensed, Int, IV, Var}}
+    constraint_storage::RV
+    dof::DofManager{Condensed, Int, IV, Var}
+    matrix_pattern::SparseMatrixPattern{IV, BlockPattern, RV}
+    vector_pattern::SparseVectorPattern{IV}
+end

src/integrals/Integrals.jl

@@ -0,0 +1,59 @@
+abstract type AbstractIntegral{
+    A         <: AbstractAssembler,
+    Integrand <: Function
+} end
+
+function integrate end
+
+struct ScalarIntegral{A, I, C <: NamedTuple} <: AbstractIntegral{A, I}
+    assembler::A
+    cache::C
+    integrand::I
+end
+
+function ScalarIntegral(asm, integrand)
+    fspace = function_space(asm.dof)
+    cache = Matrix{Float64}[]
+    for (key, val) in pairs(fspace.ref_fes)
+        NQ = ReferenceFiniteElements.num_quadrature_points(val)
+        NE = size(getfield(fspace.elem_conns, key), 2)
+        field = L2ElementField(zeros(Float64, NQ, NE))
+        push!(cache, field)
+    end
+    cache = NamedTuple{keys(fspace.ref_fes)}(tuple(scalar_quadarature_storage...))
+    return ScalarIntegral(asm, cache, integrand)
+end
+
+# function gradient(integral::ScalarIntegral)
+#     # func(physics, interps, x, t, dt, u, u_n, state_old, props) = ForwardDiff.gradient(
+#     #     z -> integral.integrand(physics, interps, x, t, dt, z, u_n, state_old, props)[1],
+#     # )
+#     function integrand_grad(physics, interps, x, t, dt, u, u_n, state_old, props)
+#         return ForwardDiff.gradient()
+#     # return VectorIntegral
+# end
+
+function integrate(integral::ScalarIntegral, U, p)
+    cache, dof = integral.cache, integral.assembler.dof
+    func = integral.integrand
+    assemble_quadrature_quantity!(cache, dof, func, U, p)
+    return mapreduce(sum, sum, values(cache))
+end
+
+struct VectorIntegral{A, I, C <: AbstractField} <: AbstractIntegral{A, I}
+    assembler::A
+    cache::C
+    integrand::I
+end
+
+function VectorIntegral(asm, integrand)
+    cache = create_field(asm)
+    return VectorIntegral(asm, cache, integrand)
+end
+
+function integrate(integral::VectorIntegral, U, p)
+    cache, dof = integral.cache, integral.assembler.dof
+    func = integral.integrand
+    assemble_vector!(cache, dof, func, U, p)
+    return cache
+end

test/TestAssemblers.jl

@@ -1,16 +1,11 @@
-f(X, _) = 2. * π^2 * sin(π * X[1]) * sin(π * X[2])
-bc_func(_, _) = 0.
-# TODO this creates a warning during testing due to
-# reinclude
-include("poisson/TestPoissonCommon.jl") 
-
-
-@testset ExtendedTestSet "Sparse matrix assembler" begin
+function test_sparse_matrix_assembler()
   # create very simple poisson problem
   mesh_file = "./poisson/poisson.g"
   mesh = UnstructuredMesh(mesh_file)
   V = FunctionSpace(mesh, H1Field, Lagrange) 
-  physics = Poisson()
+  f(X, _) = 2. * π^2 * sin(π * X[1]) * sin(π * X[2])
+  bc_func(_, _) = 0.
+  physics = Poisson(f)
   props = SVector{0, Float64}()
   u = ScalarFunction(V, :u)
   @show asm = SparseMatrixAssembler(u)
@@ -54,3 +49,7 @@ include("poisson/TestPoissonCommon.jl")
   # mainly just to test the show method for parameters
   @show p
 end
+
+@testset "Sparse matrix assembler" begin
+  test_sparse_matrix_assembler()
+end

test/TestBCs.jl

@@ -43,7 +43,7 @@ function test_neumann_bc_container_init()
   @show bc
 end
 
-@testset ExtendedTestSet "BoundaryConditions" begin
+@testset "BoundaryConditions" begin
   test_dirichlet_bc_input()
   test_dirichlet_bc_container_init()
   test_neumann_bc_input()

test/TestConnectivities.jl

@@ -1,4 +1,4 @@
-@testset ExtendedTestSet "Connectivities" begin
+function test_connectivities()
   conn_in = [
     1 5 9;
     2 6 10;
@@ -42,3 +42,7 @@
   # conn_temp = connectivity(conn, 1)
   # @test SVector{4, Int64}(conn_temp[:, 1]) ≈ conn_in[:, 1]
 end
+
+@testset "Connectivities" begin
+  test_connectivities()
+end

test/TestDofManagers.jl

@@ -61,7 +61,7 @@ end
 #   @test FiniteElementContainers.num_bcs(dof) == 5
 # end
 
-@testset ExtendedTestSet "DofManager" begin
+@testset "DofManager" begin
   # mesh = FileMesh(FiniteElementContainers.ExodusMesh, "./poisson/poisson.g")
   mesh = UnstructuredMesh("./poisson/poisson.g")
   fspace = FunctionSpace(mesh, H1Field, Lagrange)

test/TestFields.jl

@@ -1,4 +1,4 @@
-@testset ExtendedTestSet "L2 Element Field" begin
+function test_l2_element_field()
   data = rand(2, 20)
   field = L2ElementField(data)
 
@@ -31,7 +31,7 @@
   end
 end
 
-@testset ExtendedTestSet "H1 Field" begin
+function test_h1_field()
   data = rand(2, 20)
   field = H1Field(data)
 
@@ -63,7 +63,7 @@ end
   end
 end
 
-@testset ExtendedTestSet "Quadrature Field" begin
+function test_l2_quadrature_field()
   data = rand(Float64, 2, 3, 100)
   field = L2QuadratureField(data)
 
@@ -94,3 +94,9 @@ end
   @test size(field) == (0, 3, 0)
   @test axes(field) == (Base.OneTo(0), Base.OneTo(3), Base.OneTo(0))
 end
+
+@testset "Fields" begin
+  test_h1_field()
+  test_l2_element_field()
+  test_l2_quadrature_field()
+end

test/TestFormulations.jl

@@ -369,7 +369,7 @@ function test_three_dimensional(∇N_X, ∇u_q, A_q)
 end
 
 # TODO test formulations on various element types
-@testset ExtendedTestSet "Formulations" begin
+function test_formulations()
   X = [
     0.0 0.0;
     1.0 0.0;
@@ -397,3 +397,7 @@ end
   ∇u_q = SMatrix{3, 3, Float64, 9}((1., 2., 3., 4., 5., 6., 7., 8., 9.))
   test_three_dimensional(∇N_X, ∇u_q, A_q)
 end
+
+@testset "Formulations" begin
+  test_formulations()
+end