Some work towards mixed space assembly. Its really hacky right now an…

examples/cooks_membrane/generate_meshes.sh

@@ -0,0 +1,2 @@
+gmsh geometry.geo -2 -setnumber element_order 1
+gmsh geometry.geo -2 -setnumber element_order 2

examples/cooks_membrane/geometry.geo

@@ -0,0 +1,53 @@
+////////////////////////////////////////////////////////////
+// Cook's membrane (Q1 quadrilateral mesh)
+//
+// Geometry:
+// (0,44) ----------- (48,60)
+//   |                   |
+//   |                   |
+// (0,0) ------------ (48,44)
+////////////////////////////////////////////////////////////
+
+SetFactory("OpenCASCADE");
+
+// Default if nothing is supplied on the command line
+DefineConstant[
+    element_order = 1
+];
+
+// Geometry
+Point(1) = {0,  0, 0};
+Point(2) = {48,44, 0};
+Point(3) = {48,60, 0};
+Point(4) = {0,44, 0};
+
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,1};
+
+Curve Loop(1) = {1,2,3,4};
+Plane Surface(1) = {1};
+
+// Structured mesh
+nx = 32;
+ny = 32;
+
+Transfinite Curve{1,3} = nx + 1;
+Transfinite Curve{2,4} = ny + 1;
+
+Transfinite Surface{1};
+Recombine Surface{1};
+
+// Physical groups
+Physical Surface("Domain") = {1};
+
+Physical Curve("Left")   = {4};
+Physical Curve("Right")  = {2};
+Physical Curve("Bottom") = {1};
+Physical Curve("Top")    = {3};
+
+Mesh.ElementOrder = element_order;
+Mesh.SecondOrderIncomplete = 0;
+
+Mesh 2;

examples/cooks_membrane/script.jl

@@ -0,0 +1,172 @@
+import FiniteElementContainers as FEC
+using FiniteElementContainers
+using Gmsh
+using StaticArrays
+using Tensors
+
+struct TwoFieldSolidMechanics{NF, NP, NS} <: AbstractPhysics{NF, NP, NS}
+end
+
+struct Displ <: AbstractPhysics{2, 3, 0}
+end
+
+struct Pressure <: AbstractPhysics{1, 3, 0}
+end
+
+function FiniteElementContainers.create_properties(::TwoFieldSolidMechanics)
+    ρ = 1e3
+    K = 1.e9
+    G = 1.e6
+    return SVector{3, Float64}(ρ, K, G)
+end
+
+function jacobian(∇u)
+    return det(∇u + one(∇u))
+end
+
+function pk1_stress_iso(props, ∇u, p)
+    κ, μ    = props[2], props[3]
+    F       = ∇u + one(∇u)
+    J       = det(F)
+    J_m_13  = 1. / cbrt(J)
+    J_m_23  = J_m_13 * J_m_13
+    I_1     = tr(tdot(F))
+    F_inv_T = inv(F)'
+    P_iso   = μ * J_m_23 * (F - (1. / 3.) * I_1 * F_inv_T)
+    return P_iso
+end
+
+function pk1_stress_vol(props, ∇u, p)
+    κ, μ    = props[2], props[3]
+    F       = ∇u + one(∇u)
+    J       = det(F)
+    F_inv_T = inv(F)'
+    # P_vol   = 0.5 * κ * (J * J - 1.) * F_inv_T
+    P_vol   = p * J * F_inv_T
+    return P_vol
+end
+
+material_tangent_iso(props, ∇u, p) = Tensors.gradient(z -> pk1_stress_iso(props, z, p), ∇u) 
+material_tangent_vol(props, ∇u, p) = Tensors.gradient(z -> pk1_stress_vol(props, z, p), ∇u) 
+
+@inline function FiniteElementContainers.residual(
+    physics::TwoFieldSolidMechanics, interps, x_el, t, dt,
+    u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    u_el, p_el = u_el
+    interps_u, interps_p = interps
+    x_el_u, x_el_p = x_el
+    interps_u = map_interpolants(interps_u, x_el_u)
+    interps_p = map_interpolants(interps_p, x_el_p)
+    JxW_u = interps_u.JxW
+    JxW_p = interps_p.JxW
+    ∇u_q = interpolate_field_gradients(Displ(), interps_u, u_el)
+    ∇u_q = modify_field_gradients(PlaneStrain(), ∇u_q)
+    p_q = interpolate_field_values(Pressure(), interps_p, p_el)
+
+    # constitutive
+    P_iso = pk1_stress_iso(props, ∇u_q, p_q[1])
+    P_vol = pk1_stress_vol(props, ∇u_q, p_q[1])
+    J = jacobian(∇u_q)
+
+    P_q = extract_stress(PlaneStrain(), P_iso + P_vol)
+    G_q = discrete_gradient(PlaneStrain(), interps_u.∇N_X)
+    R_u = JxW_u * G_q * P_q
+    R_p = JxW_p * (J - one(J)) * interps_p.N
+    return R_u, R_p
+end
+
+@inline function FiniteElementContainers.stiffness(
+    physics::TwoFieldSolidMechanics, interps, x_el, t, dt,
+    u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    u_el, p_el = u_el
+    interps_u, interps_p = interps
+    x_el_u, x_el_p = x_el
+    interps_u = map_interpolants(interps_u, x_el_u)
+    interps_p = map_interpolants(interps_p, x_el_p)
+    JxW_u = interps_u.JxW
+    JxW_p = interps_p.JxW
+    ∇u_q = interpolate_field_gradients(Displ(), interps_u, u_el)
+    ∇u_q = modify_field_gradients(PlaneStrain(), ∇u_q)
+    p_q = interpolate_field_values(Pressure(), interps_p, p_el)
+    J_q = jacobian(∇u_q)
+    F_q = ∇u_q + one(∇u_q)
+    F_inv_T_q = inv(F_q)'
+    dPdp_q = extract_stress(PlaneStrain(), J_q * F_inv_T_q)
+    A_iso = material_tangent_iso(props, ∇u_q, p_q[1])
+    A_vol = material_tangent_vol(props, ∇u_q, p_q[1])
+    G_q = discrete_gradient(PlaneStrain(), interps_u.∇N_X)
+    G_pu_x = J_q .* (
+        F_inv_T_q[1, 1] * interps_u.∇N_X[: ,1] +
+        F_inv_T_q[2, 1] * interps_u.∇N_X[: ,2]
+    )
+
+    G_pu_y = J_q * (
+        F_inv_T_q[1, 2] * interps_u.∇N_X[:, 1] +
+        F_inv_T_q[2, 2] * interps_u.∇N_X[:, 2]
+    )
+    Nd = length(G_pu_x)
+
+    tup = MVector{2 * Nd, eltype(G_pu_x)}(undef)
+
+    for i in 1:Nd
+        tup[2i-1] = G_pu_x[i]
+        tup[2i]   = G_pu_y[i]
+    end
+
+    G_pu_q = SVector{2 * Nd, eltype(G_pu_x)}(tup)
+    K_uu = JxW_u * G_q * extract_stiffness(PlaneStrain(), A_iso + A_vol) * G_q'
+    K_up = JxW_u * G_q * dPdp_q * interps_p.N'
+    K_pu = JxW_p * interps_p.N * G_pu_q'
+    K_pp = zero(SMatrix{length(p_el), length(p_el), Float64, length(p_el)^2})
+    return (
+        (K_uu, K_up),
+        (K_pu, K_pp)
+    )
+end
+
+mesh_u = UnstructuredMesh(Base.source_dir() * "/geometry_q2.geo")
+mesh_p = UnstructuredMesh(Base.source_dir() * "/geometry_q1.geo")
+
+V_u = FunctionSpace(mesh_u, H1Field, Lagrange)
+V_p = FunctionSpace(mesh_p, H1Field, Lagrange)
+# V_J = FunctionSpace(mesh_p, H1Field, Lagrange)
+
+u = VectorFunction(V_u, "displ")
+p = ScalarFunction(V_p, "pressure")
+# J = ScalarFunction(V_J, "jacobian")
+
+zero_func(_, _) = 0.0
+displ_func(_, t) = 0.1 * t
+dbcs_u = DirichletBC[
+    DirichletBC("displ_x", zero_func; nodeset_name = "Left")
+    DirichletBC("displ_y", zero_func; nodeset_name = "Left")
+    DirichletBC("displ_x", zero_func; nodeset_name = "Right")
+    DirichletBC("displ_y", displ_func; nodeset_name = "Right")
+]
+physics = TwoFieldSolidMechanics{3, 0, 0}()
+props = create_properties(physics)
+times = TimeStepper(0.0, 1.0, 20)
+
+dof_u, dof_p = DofManager(u), DofManager(p)
+dof = (dof_u, dof_p)
+# dof_u, dof_p, dof_J = DofManager(u), DofManager(p), DofManager(J)
+# dof = (dof_u, dof_p, dof_J)
+asm = FEC.BlockSparseMatrixAssembler(dof)
+
+p_u    = create_parameters(mesh_u, SparseMatrixAssembler(dof_u), physics, props; dirichlet_bcs = dbcs_u, times = times)
+p_p    = create_parameters(mesh_p, SparseMatrixAssembler(dof_p), physics, props; times = times)
+params = (p_u, p_p)
+
+FEC.update_dofs!(
+    asm,
+    (p_u.dirichlet_bcs, p_p.dirichlet_bcs),
+    (p_u.periodic_bcs, p_p.periodic_bcs)
+)
+
+Uu = create_unknowns(asm)
+U  = create_field(asm)
+
+assemble_vector!(asm, residual, Uu, params)
+assemble_matrix!(asm, stiffness, Uu, params)

src/Enzyme.jl

@@ -144,9 +144,7 @@ function _assemble_vector_block_enzyme_safe!(
   for e in 1:conns_all.nelems[b]
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old)
-
     props_el = properties(props, e, b)
-  #   # val_el = _element_scratch(return_type, ref_fe, U)
 
     for q in 1:num_cell_quadrature_points(ref_fe)
       interps = _cell_interpolants(ref_fe, q)

src/FiniteElementContainers.jl

@@ -222,6 +222,7 @@ include("Constraints.jl")
 include("InitialConditions.jl")
 
 include("Formulations.jl")
+include("Gather.jl")
 include("Physics.jl")
 include("assemblers/Assemblers.jl")
 #

src/FunctionSpaces.jl

@@ -42,6 +42,8 @@ const _el_name_to_juliac_safe_id = Dict{String, Int}(
   "TETRA4"  => 6,
   "TETRA10" => 7
 )
+# Lagrange elements
+# defaulting to fully integrated elements for now
 const _juliac_safe_ref_fes = (
   ReferenceFE(Hex{Lagrange, 1}(), GaussLobattoLegendre(2, 2)),  # HEX8 for Lagrange
   ReferenceFE(Quad{Lagrange, 1}(), GaussLobattoLegendre(2, 2)), # QUAD4 for Lagrange
@@ -137,7 +139,6 @@ struct FunctionSpace{
   block_to_ref_fe_id::BTRE
   coords::H1Field{RT, RV, ND}
   elem_conns::Connectivity{IT, IV}
-  # need to remove this mabye?
   elem_id_maps::Vector{Vector{IT}} # TODO create new type for ID map similar to connectivity
   elem_to_facets::Union{Nothing, Connectivity{IT, IV}}
   facet_orientation::Union{Nothing, Connectivity{IT, IV}}
@@ -146,7 +147,7 @@ struct FunctionSpace{
   ref_fes::RefFEs
 
   function FunctionSpace{is_juliac_safe, FT}(
-    block_names, block_to_ref_fe_id, coords, elem_conns,
+    block_names, block_to_ref_fe_id, coords, elem_conns, 
     elem_id_maps, elem_to_facets, facet_orientation, node_id_map, ref_fes
   ) where {is_juliac_safe, FT}
     new{
@@ -272,7 +273,7 @@ function block_entity_size(fspace::FunctionSpace, b::Int)
   return (num_entities_per_element(fspace, b), num_elements(fspace, b))
 end
 
-function block_reference_element(fspace::FunctionSpace{false, FT, I, V, BTRE, C, R}, block_id::Int) where {I, FT, V, BTRE, C, R}
+function block_reference_element(fspace::FunctionSpace{false, FT, I, V, BTRE, C, R}, block_id::Int) where {FT, I, V, BTRE, C, R}
   return fspace.ref_fes[block_id]
 end
 
@@ -360,7 +361,6 @@ function num_entities(fspace::FunctionSpace)
 end
 
 function num_entities_per_element(fspace::FunctionSpace, b::Int)
-  return num_entities_per_element(fspace.elem_conns, b)
   if _field_type(fspace) == L2Field
     return block_quadrature_size(fspace, b)[1]
   else