Merge pull request #67 from Cthonios/formulation-allocation-fix

fixing issues in large allocations for elements with nodes with large…
Cthonios · Nov 17, 2024 · 6ae2aff · 6ae2aff
2 parents 3c200cb + 5f943eb
commit 6ae2aff
Show file tree

Hide file tree

Showing 4 changed files with 75 additions and 65 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "FiniteElementContainers"
 uuid = "d08262e4-672f-4e7f-a976-f2cea5767631"
 authors = ["Craig M. Hamel <cmhamel32@gmail.com> and contributors"]
-version = "0.5.8"
+version = "0.5.9"
 
 [deps]
 Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"

diff --git a/src/formulations/IncompressiblePlaneStress.jl b/src/formulations/IncompressiblePlaneStress.jl
@@ -10,7 +10,8 @@ $(TYPEDSIGNATURES)
 """
 function discrete_gradient(::IncompressiblePlaneStress, ∇N_X)
   N   = size(∇N_X, 1)
-  tup = ntuple(i -> 0.0, Val(4 * 2 * N))
+  # tup = ntuple(i -> 0.0, Val(4 * 2 * N))
+  tup = zeros(SVector{4 * 2 * N, eltype(∇N_X)})
 
   for n in 1:N
     k = 2 * (n - 1) 
@@ -30,15 +31,16 @@ function discrete_gradient(::IncompressiblePlaneStress, ∇N_X)
     tup = setindex(tup, ∇N_X[n, 2], k + 2)
   end
 
-  return SMatrix{2 * N, 4, eltype(∇N_X), 2 * N * 4}(tup)
+  return SMatrix{2 * N, 4, eltype(∇N_X), 2 * N * 4}(tup.data)
 end
 
 """
 $(TYPEDSIGNATURES)
 """
 function discrete_symmetric_gradient(::IncompressiblePlaneStress, ∇N_X)
   N   = size(∇N_X, 1)
-  tup = ntuple(i -> 0.0, Val(3 * 2 * N))
+  # tup = ntuple(i -> 0.0, Val(3 * 2 * N))
+  tup = zeros(SVector{3 * 2 * N, eltype(∇N_X)})
 
   for n in 1:N
     k = 2 * (n - 1) 
@@ -54,15 +56,16 @@ function discrete_symmetric_gradient(::IncompressiblePlaneStress, ∇N_X)
     tup = setindex(tup, ∇N_X[n, 1], k + 2)
   end
 
-  return SMatrix{2 * N, 3, eltype(∇N_X), 2 * N * 3}(tup)
+  return SMatrix{2 * N, 3, eltype(∇N_X), 2 * N * 3}(tup.data)
 end
 
 """
 $(TYPEDSIGNATURES)
 """
 function discrete_values(::IncompressiblePlaneStress, N)
   N_nodes = size(N, 1)
-  tup = ntuple(i -> 0.0, Val(2 * N_nodes))
+  # tup = ntuple(i -> 0.0, Val(2 * N_nodes))
+  tup = zeros(SVector{2 * N_nodes, eltype(N)})
 
   for n in 1:N_nodes
     tup = setindex(tup, N[n], n)
@@ -72,7 +75,8 @@ function discrete_values(::IncompressiblePlaneStress, N)
     tup = setindex(tup, N[n], n + N_nodes)
   end
 
-  return SVector{2 * N_nodes, eltype(N)}(tup)
+  # return SVector{2 * N_nodes, eltype(N)}(tup)
+  return tup
 end
 
 """

diff --git a/src/formulations/PlaneStrain.jl b/src/formulations/PlaneStrain.jl
@@ -10,7 +10,7 @@ $(TYPEDSIGNATURES)
 """
 function discrete_gradient(::PlaneStrain, ∇N_X)
   N   = size(∇N_X, 1)
-  tup = ntuple(i -> 0.0, Val(4 * 2 * N))
+  tup = zeros(SVector{4 * 2 * N, eltype(∇N_X)})
 
   for n in 1:N
     k = 2 * (n - 1) 
@@ -29,16 +29,16 @@ function discrete_gradient(::PlaneStrain, ∇N_X)
     tup = setindex(tup, 0.0,        k + 1)
     tup = setindex(tup, ∇N_X[n, 2], k + 2)
   end
-
-  return SMatrix{2 * N, 4, eltype(∇N_X), 2 * N * 4}(tup)
+  return SMatrix{2 * N, 4, eltype(∇N_X), 2 * N * 4}(tup.data)
 end
 
 """
 $(TYPEDSIGNATURES)
 """
 function discrete_values(::PlaneStrain, N)
   N_nodes = size(N, 1)
-  tup = ntuple(i -> 0.0, Val(2 * N_nodes))
+  # tup = ntuple(i -> 0.0, Val(2 * N_nodes))
+  tup = zeros(SVector{2 * N_nodes, eltype(N)})
 
   for n in 1:N_nodes
     tup = setindex(tup, N[n], n)
@@ -48,15 +48,17 @@ function discrete_values(::PlaneStrain, N)
     tup = setindex(tup, N[n], n + N_nodes)
   end
 
-  return SVector{2 * N_nodes, eltype(N)}(tup)
+  # return SVector{2 * N_nodes, eltype(N)}(tup)
+  return tup
 end
 
 """
 $(TYPEDSIGNATURES)
 """
 function discrete_symmetric_gradient(::PlaneStrain, ∇N_X)
   N   = size(∇N_X, 1)
-  tup = ntuple(i -> 0.0, Val(3 * 2 * N))
+  # tup = ntuple(i -> 0.0, Val(3 * 2 * N))
+  tup = zeros(SVector{3 * 2 * N, eltype(∇N_X)})
 
   for n in 1:N
     k = 2 * (n - 1) 
@@ -72,7 +74,7 @@ function discrete_symmetric_gradient(::PlaneStrain, ∇N_X)
     tup = setindex(tup, ∇N_X[n, 1], k + 2)
   end
 
-  return SMatrix{2 * N, 3, eltype(∇N_X), 2 * N * 3}(tup)
+  return SMatrix{2 * N, 3, eltype(∇N_X), 2 * N * 3}(tup.data)
 end
 
 """

diff --git a/src/formulations/ThreeDimensional.jl b/src/formulations/ThreeDimensional.jl
@@ -10,108 +10,111 @@ $(TYPEDSIGNATURES)
 """
 function discrete_gradient(::ThreeDimensional, ∇N_X)
   N   = size(∇N_X, 1)
-  tup = ntuple(i -> 0.0, Val(9 * 3 * N))
+  # tup = ntuple(i -> 0.0, Val(9 * 3 * N))
+  tup = zeros(SVector{9 * 3 * N, eltype(∇N_X)})
 
   for n in 1:N
     k = 3 * (n - 1) 
-    tup = set(tup, ∇N_X[n, 1], k + 1)
-    tup = set(tup, 0.0,        k + 2)
-    tup = set(tup, 0.0,        k + 3)
+    tup = setindex(tup, ∇N_X[n, 1], k + 1)
+    tup = setindex(tup, 0.0,        k + 2)
+    tup = setindex(tup, 0.0,        k + 3)
 
     k = 3 * (n - 1) + 3 * N
-    tup = set(tup, 0.0,        k + 1)
-    tup = set(tup, ∇N_X[n, 1], k + 2)
-    tup = set(tup, 0.0,        k + 3)
+    tup = setindex(tup, 0.0,        k + 1)
+    tup = setindex(tup, ∇N_X[n, 1], k + 2)
+    tup = setindex(tup, 0.0,        k + 3)
 
     k = 3 * (n - 1) + 2 * 3 * N
-    tup = set(tup, 0.0,        k + 1)
-    tup = set(tup, 0.0,        k + 2)
-    tup = set(tup, ∇N_X[n, 1], k + 3)
+    tup = setindex(tup, 0.0,        k + 1)
+    tup = setindex(tup, 0.0,        k + 2)
+    tup = setindex(tup, ∇N_X[n, 1], k + 3)
 
     #
     k = 3 * (n - 1) + 3 * 3 * N
-    tup = set(tup, ∇N_X[n, 2], k + 1)
-    tup = set(tup, 0.0,        k + 2)
-    tup = set(tup, 0.0,        k + 3)
+    tup = setindex(tup, ∇N_X[n, 2], k + 1)
+    tup = setindex(tup, 0.0,        k + 2)
+    tup = setindex(tup, 0.0,        k + 3)
 
     k = 3 * (n - 1) + 4 * 3 * N
-    tup = set(tup, 0.0,        k + 1)
-    tup = set(tup, ∇N_X[n, 2], k + 2)
-    tup = set(tup, 0.0,        k + 3)
+    tup = setindex(tup, 0.0,        k + 1)
+    tup = setindex(tup, ∇N_X[n, 2], k + 2)
+    tup = setindex(tup, 0.0,        k + 3)
 
     k = 3 * (n - 1) + 5 * 3 * N
-    tup = set(tup, 0.0,        k + 1)
-    tup = set(tup, 0.0,        k + 2)
-    tup = set(tup, ∇N_X[n, 2], k + 3)
+    tup = setindex(tup, 0.0,        k + 1)
+    tup = setindex(tup, 0.0,        k + 2)
+    tup = setindex(tup, ∇N_X[n, 2], k + 3)
 
     #
     k = 3 * (n - 1) + 6 * 3 * N
-    tup = set(tup, ∇N_X[n, 3], k + 1)
-    tup = set(tup, 0.0,        k + 2)
-    tup = set(tup, 0.0,        k + 3)
+    tup = setindex(tup, ∇N_X[n, 3], k + 1)
+    tup = setindex(tup, 0.0,        k + 2)
+    tup = setindex(tup, 0.0,        k + 3)
 
     k = 3 * (n - 1) + 7 * 3 * N
-    tup = set(tup, 0.0,        k + 1)
-    tup = set(tup, ∇N_X[n, 3], k + 2)
-    tup = set(tup, 0.0,        k + 3)
+    tup = setindex(tup, 0.0,        k + 1)
+    tup = setindex(tup, ∇N_X[n, 3], k + 2)
+    tup = setindex(tup, 0.0,        k + 3)
 
     k = 3 * (n - 1) + 8 * 3 * N
-    tup = set(tup, 0.0,        k + 1)
-    tup = set(tup, 0.0,        k + 2)
-    tup = set(tup, ∇N_X[n, 3], k + 3)
+    tup = setindex(tup, 0.0,        k + 1)
+    tup = setindex(tup, 0.0,        k + 2)
+    tup = setindex(tup, ∇N_X[n, 3], k + 3)
   end
 
-  return SMatrix{3 * N, 9, eltype(∇N_X), 3 * N * 9}(tup)
+  return SMatrix{3 * N, 9, eltype(∇N_X), 3 * N * 9}(tup.data)
 end
 
 """
 $(TYPEDSIGNATURES)
 """
 function discrete_symmetric_gradient(::ThreeDimensional, ∇N_X)
   N   = size(∇N_X, 1)
-  tup = ntuple(i -> 0.0, Val(6 * 3 * N))
+  # tup = ntuple(i -> 0.0, Val(6 * 3 * N))
+  tup = zeros(SVector{6 * 3 * N, eltype(∇N_X)})
 
   for n in 1:N
     k = 3 * (n - 1)
-    tup = set(tup, ∇N_X[n, 1], k + 1)
-    tup = set(tup, 0.0,        k + 2)
-    tup = set(tup, 0.0,        k + 3)
+    tup = setindex(tup, ∇N_X[n, 1], k + 1)
+    tup = setindex(tup, 0.0,        k + 2)
+    tup = setindex(tup, 0.0,        k + 3)
 
     k = 3 * (n - 1) + 3 * N
-    tup = set(tup, 0.0,        k + 1)
-    tup = set(tup, ∇N_X[n, 2], k + 2)
-    tup = set(tup, 0.0,        k + 3)
+    tup = setindex(tup, 0.0,        k + 1)
+    tup = setindex(tup, ∇N_X[n, 2], k + 2)
+    tup = setindex(tup, 0.0,        k + 3)
 
     k = 3 * (n - 1) + 2 * 3 * N
-    tup = set(tup, 0.0,        k + 1)
-    tup = set(tup, 0.0,        k + 2)
-    tup = set(tup, ∇N_X[n, 3], k + 3)
+    tup = setindex(tup, 0.0,        k + 1)
+    tup = setindex(tup, 0.0,        k + 2)
+    tup = setindex(tup, ∇N_X[n, 3], k + 3)
 
     k = 3 * (n - 1) + 3 * 3 * N
-    tup = set(tup, ∇N_X[n, 2], k + 1)
-    tup = set(tup, ∇N_X[n, 1], k + 2)
-    tup = set(tup, 0.0,        k + 3)
+    tup = setindex(tup, ∇N_X[n, 2], k + 1)
+    tup = setindex(tup, ∇N_X[n, 1], k + 2)
+    tup = setindex(tup, 0.0,        k + 3)
 
     k = 3 * (n - 1) + 4 * 3 * N
-    tup = set(tup, 0.0,        k + 1)
-    tup = set(tup, ∇N_X[n, 3], k + 2)
-    tup = set(tup, ∇N_X[n, 2], k + 3)
+    tup = setindex(tup, 0.0,        k + 1)
+    tup = setindex(tup, ∇N_X[n, 3], k + 2)
+    tup = setindex(tup, ∇N_X[n, 2], k + 3)
 
     k = 3 * (n - 1) + 5 * 3 * N
-    tup = set(tup, ∇N_X[n, 3], k + 1)
-    tup = set(tup, 0.0,        k + 2)
-    tup = set(tup, ∇N_X[n, 1], k + 3)
+    tup = setindex(tup, ∇N_X[n, 3], k + 1)
+    tup = setindex(tup, 0.0,        k + 2)
+    tup = setindex(tup, ∇N_X[n, 1], k + 3)
 
   end
-  return SMatrix{3 * N, 6, eltype(∇N_X), 3 * N * 6}(tup)
+  return SMatrix{3 * N, 6, eltype(∇N_X), 3 * N * 6}(tup.data)
 end
 
 """
 $(TYPEDSIGNATURES)
 """
 function discrete_values(::ThreeDimensional, N)
   N_nodes = size(N, 1)
-  tup = ntuple(i -> 0.0, Val(3 * N_nodes))
+  # tup = ntuple(i -> 0.0, Val(3 * N_nodes))
+  tup = zeros(SVector{3 * N_nodes, eltype(N)})
 
   for n in 1:N_nodes
     tup = setindex(tup, N[n], n)
@@ -125,7 +128,8 @@ function discrete_values(::ThreeDimensional, N)
     tup = setindex(tup, N[n], n + 2 * N_nodes)
   end
 
-  return SVector{3 * N_nodes, eltype(N)}(tup)
+  # return SVector{3 * N_nodes, eltype(N)}(tup)
+  return tup
 end
 
 """