Another in place scatter method.

examples/cooks_membrane/script.jl

@@ -165,8 +165,10 @@ FEC.update_dofs!(
     (p_u.periodic_bcs, p_p.periodic_bcs)
 )
 
-pp_u = PostProcessor(mesh_u, "u.exo", u)
-pp_p = PostProcessor(mesh_p, "p.exo", p)
+# pp_u = PostProcessor(mesh_u, "u.exo", u)
+# pp_p = PostProcessor(mesh_p, "p.exo", p)
+
+
 # solver = NewtonSolver(DirectLinearSolver(asm))
 # integrator = QuasiStaticIntegrator(solver)
 
@@ -186,39 +188,49 @@ Uu = create_unknowns(asm)
 # 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)
+function test_func(asm, u, p)
+    FiniteElementContainers.update_time!(p)
+    FiniteElementContainers.update_bc_values!(p, asm)
+    assemble_vector!(asm, residual, u, p)
+    return nothing
+end
+
+@time test_func(asm, Uu, params)
+@time test_func(asm, Uu, params)
+
+# 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
+#         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)
+#         assemble_stiffness!(asm, stiffness, Uu, params)
+#         K = stiffness(asm)
 
-        ΔU = K \ R
-        Uu .-= ΔU
+#         ΔU = K \ R
+#         Uu .-= ΔU
 
-        println("iter = $iter, |R| = $(rnorm / r0), |ΔU| = $(norm(ΔU))")
+#         println("iter = $iter, |R| = $(rnorm / r0), |ΔU| = $(norm(ΔU))")
 
-    end
+#     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)
+#     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

examples/poisson/script.jl

@@ -35,7 +35,7 @@ physics = Poisson(f)
 props = create_properties(physics)
 u = ScalarFunction(V, "u")
 # asm = SparseMatrixAssembler(u; use_condensed = true)
-asm = SparseMatrixAssembler(u)
+asm = SparseMatrixAssembler(u; use_inplace_methods = true)
 
 
 # dbcs = nothing

src/Formulations.jl

@@ -191,6 +191,99 @@ function scatter_with_gradients_and_gradients!(
   return nothing
 end
 
+"""
+General stiffness like matrix scatter method
+"""
+function scatter_with_gradients_and_gradients!(
+  storage::AbstractVector,
+  form::AbstractElementFormulation{ND, NF},
+  e,
+  conns,
+  ∇N_X,
+  K::AbstractMatrix{T},
+) where {ND, NF, T <: Number}
+
+  @assert size(K, 1) == NF^2
+  @assert size(K, 2) == NF^2
+  @assert ND == size(∇N_X, 2)
+
+  N    = size(∇N_X, 1)
+  NDPE = N * NF
+
+  start_id = (e - 1) * NDPE^2 + 1
+  ids      = start_id:(start_id + NDPE^2 - 1)
+
+  inc = 1
+
+  for n2 in 1:N
+    for d2 in 1:NF
+      for n1 in 1:N
+        for d1 in 1:NF
+          contrib = zero(T)
+
+          for J in 1:ND
+            a = (J - 1) * ND + d1
+
+            for L in 1:ND
+              b = (L - 1) * ND + d2
+              contrib += ∇N_X[n1, J] * K[a, b] * ∇N_X[n2, L]
+            end
+          end
+
+          storage[ids[inc]] += contrib
+          inc += 1
+        end
+      end
+    end
+  end
+
+  return nothing
+end
+
+function scatter_gradient_gradient!(
+  storage::AbstractField,
+  form::AbstractElementFormulation{ND, NF},
+  e,
+  conns,
+  ∇N_X,
+  A::AbstractMatrix{T},
+  v_el::AbstractVector{T},
+) where {ND, NF, T <: Number}
+
+  @assert ND == size(∇N_X, 2)
+  @assert size(A) == (ND * NF, ND * NF)
+
+  N = size(∇N_X, 1)
+  NDPE = N * NF
+
+  @assert length(v_el) == NDPE
+
+  for n1 in 1:N
+    for d1 in 1:NF
+      global_id = NF * (conns[n1] - 1) + d1
+      contrib = zero(T)
+
+      for n2 in 1:N
+        for d2 in 1:NF
+          local_id = NF * (n2 - 1) + d2
+          for j1 in 1:ND
+            a = (j1 - 1) * NF + d1
+
+            for j2 in 1:ND
+                b = (j2 - 1) * NF + d2
+                contrib += ∇N_X[n1, j1] * A[a, b] * ∇N_X[n2, j2] * v_el[local_id]
+            end
+          end
+        end
+      end
+
+      fec_atomic_add!(storage, global_id, contrib)
+    end
+  end
+
+  return nothing
+end
+
 # implement for those that have it
 """
 Scalar equation specialization

test/poisson/TestPoissonCommon.jl

@@ -1,6 +1,8 @@
 # since this type of physics has no state
 # just return the old state which will be empty
 
+using ReferenceFiniteElements
+
 struct Poisson{F <: Function} <: AbstractPhysics{1, 0, 0}
   func::F
 end
@@ -91,7 +93,6 @@ end
   interps = map_interpolants(interps, x_el)
   (; X_q, N, ∇N_X, JxW) = interps
   ∇u_q = interpolate_field_gradients(physics, interps, u_el)
-  # R_q = ∇u_q * ∇N_X' - N' * physics.func(X_q, 0.0)
   form = GeneralFormulation{size(X_q, 1), num_fields(physics)}()
   scatter_with_gradients!(storage, form, e, conn, ∇N_X, JxW * ∇u_q)
   scatter_with_values!(storage, form, e, conn, N, -JxW * physics.func(X_q, 0.0))