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

Project.toml

@@ -1,12 +1,13 @@
 name = "FiniteElementContainers"
 uuid = "d08262e4-672f-4e7f-a976-f2cea5767631"
-version = "0.14.2"
+version = "0.14.3"
 authors = ["Craig M. Hamel <cmhamel32@gmail.com> and contributors"]
 
 [deps]
 AcceleratedKernels = "6a4ca0a5-0e36-4168-a932-d9be78d558f1"
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 DynamicExpressions = "a40a106e-89c9-4ca8-8020-a735e8728b6b"
 Exodus = "f57ae99e-f805-4780-bdca-96e224be1e5a"
@@ -26,7 +27,6 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 [weakdeps]
 AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
 AbaqusReader = "bc6b9049-e460-56d6-94b4-a597b2c0390d"
-BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
 Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
@@ -35,7 +35,6 @@ PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 [extensions]
 AMDGPUExt = "AMDGPU"
 AbaqusReaderExt = "AbaqusReader"
-BlockArraysExt = "BlockArrays"
 CUDAExt = "CUDA"
 GmshExt = "Gmsh"
 MetisExt = "Metis"

examples/cooks_membrane/geometry_q1.geo

@@ -0,0 +1,48 @@
+////////////////////////////////////////////////////////////
+// Cook's membrane (Q1 quadrilateral mesh)
+//
+// Geometry:
+// (0,44) ----------- (48,60)
+//   |                   |
+//   |                   |
+// (0,0) ------------ (48,44)
+////////////////////////////////////////////////////////////
+
+SetFactory("OpenCASCADE");
+
+// 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 = 1;
+Mesh.SecondOrderIncomplete = 0;
+
+Mesh 2;

examples/cooks_membrane/geometry_q2.geo

@@ -0,0 +1,44 @@
+////////////////////////////////////////////////////////////
+// Cook's membrane (Q2 quadrilateral mesh)
+////////////////////////////////////////////////////////////
+
+SetFactory("OpenCASCADE");
+
+// 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};
+
+// Second-order quadrilateral elements (9-node)
+Mesh.ElementOrder = 2;
+Mesh.SecondOrderIncomplete = 0;
+Mesh.HighOrderOptimize = 1;
+
+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)