Some work towards block matrix assembly with a very basic but still v…

.gitignore

@@ -15,3 +15,6 @@ Manifest.toml
 *.e.*.*
 *.out
 build/
+*.exo
+*.exo.*
+*.msh

examples/cooks_membrane/generate_meshes.sh

@@ -0,0 +1,2 @@
+gmsh geometry.geo -2 -setnumber element_order 1 -o mesh_p1.msh
+gmsh geometry.geo -2 -setnumber element_order 2 -o mesh_p2.msh

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,224 @@
+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_el, ∇u_q, p_q[1])
+    P_vol = pk1_stress_vol(props_el, ∇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_el, ∇u_q, p_q[1])
+    A_vol = material_tangent_vol(props_el, ∇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[1, 2] .* interps_u.∇N_X[:, 2]
+    )
+    G_pu_y = J_q .* (
+        F_inv_T_q[2, 1] .* 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() * "/mesh_p2.msh")
+mesh_p = UnstructuredMesh(Base.source_dir() * "/mesh_p1.msh")
+
+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.01 * 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, 10)
+
+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)
+)
+
+pp_u = PostProcessor(mesh_u, "u.exo", u)
+pp_p = PostProcessor(mesh_p, "p.exo", p)
+# solver = NewtonSolver(DirectLinearSolver(asm))
+# integrator = QuasiStaticIntegrator(solver)
+
+# for n in 1:20
+#     evolve!(integrator, params)
+# end
+
+Uu = create_unknowns(asm)
+# # U  = create_field(asm)
+
+# assemble_stiffness!(asm, stiffness, Uu, params)
+# K = stiffness(asm)
+
+# K_up = K.blocks[1, 2]
+# K_pu = K.blocks[2, 1]
+
+# temp = K_up .- K_pu'
+# display(K_up .- K_pu')
+
+for n in 1:10
+    FiniteElementContainers.update_time!(params)
+    FiniteElementContainers.update_bc_values!(params, asm)
+    r0 = -1e6
+    for iter in 1:10
+        assemble_vector!(asm, residual, Uu, params)
+        R = residual(asm)
+
+        rnorm = norm(R)
+        if iter == 1
+            r0 = rnorm
+        end
+
+        if rnorm / r0 < 1e-8
+            break
+        end
+
+        assemble_stiffness!(asm, stiffness, Uu, params)
+        K = stiffness(asm)
+
+        ΔU = K \ R
+        Uu .-= ΔU
+
+        println("iter = $iter, |R| = $(rnorm / r0), |ΔU| = $(norm(ΔU))")
+
+    end
+
+    write_times(pp_u, n + 1, params[1].times.time_current)
+    write_times(pp_p, n + 1, params[2].times.time_current)
+    write_field(pp_u, n + 1, ("displ_x", "displ_y"), params[1].field)
+    write_field(pp_p, n + 1, ("pressure",), params[2].field)
+end
+close(pp_u)
+close(pp_p)
+
+# Δu = K \ R

src/Parameters.jl

@@ -427,6 +427,21 @@ function update_bc_values!(p::AbstractParameters, assembler)
   return nothing
 end
 
+function update_bc_values!(params::Tuple, assembler)
+  for n in 1:length(params)
+    p = params[n]
+    X = coordinates(p)
+    t = current_time(p)
+    update_bc_values!(p.dirichlet_bcs, X, t)
+    # update_bc_values!(p.neumann_bcs, assembler, X, t)
+    # update_bc_values!(p.periodic_bcs, X, t)
+    # # update_bc_values!(p.robin_bcs, assembler, X, t, p.field)
+    # update_source_values!(p.sources, assembler, X, t)
+
+    # update_bc_values!(p[n], assembler)
+  end
+end
+
 function update_bc_values!(p::TypeStableParameters, assembler)
   X = coordinates(p)
   t = current_time(p)
@@ -454,6 +469,18 @@ function update_dofs!(asm::AbstractAssembler, p::Parameters)
   return nothing
 end
 
+function _update_field!(p::AbstractParameters)
+  p.field_old.data .= p.field.data
+  return nothing
+end
+
+function _update_field!(p::Tuple)
+  for n in 1:length(p)
+    _update_field!(p[n])
+  end
+  return nothing
+end
+
 function _update_for_assembly!(p::AbstractParameters, dof::DofManager, Uu)
   update_field_dirichlet_bcs!(p.field, p.dirichlet_bcs)
   update_field_unknowns!(p.field, dof, Uu)
@@ -477,6 +504,12 @@ function _update_for_assembly!(p::AbstractParameters, dof::DofManager, Uu, Vu)
   return nothing
 end
 
+function _update_for_assembly!(p, dof::Tuple, Uu)
+  for n in 1:length(p)
+    _update_for_assembly!(p[n], dof[n], view(Uu, BlockArrays.Block(n)))
+  end
+end
+
 # Full-DOF flavor: caller is responsible for assembling the merged
 # vectors U_full = [Uu; U_BC] and v_full = [v_free; v_BC] themselves.
 # Unlike the free-DOF flavors above, we do NOT call
@@ -501,3 +534,10 @@ function update_time!(p::AbstractParameters)
   p.times.time_current = current_time(p.times) + time_step(p.times)
   return nothing
 end
+
+function update_time!(p::Tuple)
+  for n in 1:length(p)
+    update_time!(p[n])
+  end
+  return nothing
+end

src/Solvers.jl

@@ -49,6 +49,17 @@ struct DirectLinearSolver{
     # what's the best way to do this with general solvers?
     ΔUu::U
 
+    function DirectLinearSolver(assembler::BlockSparseMatrixAssembler)
+        preconditioner = I
+        ΔUu = similar(assembler.residual_unknowns)
+        fill!(ΔUu, zero(eltype(ΔUu)))
+        new{typeof(assembler), typeof(preconditioner), typeof(ΔUu)}(
+            assembler, preconditioner, 
+            DirectLinearSolverSettings(), 
+            TimerOutput(), ΔUu
+        )
+    end 
+
     function DirectLinearSolver(assembler::SparseMatrixAssembler)
         preconditioner = I
         ΔUu = similar(assembler.residual_unknowns)