Update topology example

docs/Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "8722b6a67abf1b81e2c7808ddfcfb1d737cd0040"
+project_hash = "6f75568b4b6adf310d9c6b4556b4e69e4f996f31"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "3a6511b6e54550bcbc986c560921a8cd7761fcd8"
@@ -103,6 +103,18 @@ weakdeps = ["SparseArrays"]
     [deps.ArrayLayouts.extensions]
     ArrayLayoutsSparseArraysExt = "SparseArrays"
 
+[[deps.ArraysOfArrays]]
+deps = ["Statistics"]
+git-tree-sha1 = "ad745f4bc2cfc5fe58484e53e49c3e650103e4bd"
+uuid = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
+version = "0.6.5"
+weakdeps = ["Adapt", "ChainRulesCore", "StaticArraysCore"]
+
+    [deps.ArraysOfArrays.extensions]
+    ArraysOfArraysAdaptExt = "Adapt"
+    ArraysOfArraysChainRulesCoreExt = "ChainRulesCore"
+    ArraysOfArraysStaticArraysCoreExt = "StaticArraysCore"
+
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 

docs/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"

docs/src/literate-gallery/topology_optimization.jl

@@ -15,7 +15,7 @@
 # ## Introduction
 #
 # Topology optimization is the task of finding structures that are mechanically ideal.
-# In this example we cover the bending beam, where we specify a load, boundary conditions and the total mass. Then, our
+# In this example we cover the bending beam, where we specify a load, boundary conditions, and the total mass. Then, our
 # objective is to find the most suitable geometry within the design space minimizing the compliance (i.e. the inverse stiffness) of the structure.
 # We shortly introduce our simplified model for regular meshes. A detailed derivation of the method and advanced techniques
 # can be found in [JanHacJun2019regularizedthermotopopt](@cite) and
@@ -47,28 +47,28 @@
 # the value of the density field $\chi$. The advantage of this "split" approach is the very high computation speed.
 # The evolution equation consists of the driving force, the damping parameter $\eta$, the regularization parameter $\beta$ times the Laplacian,
 # which is necessary to avoid numerical issues like mesh dependence or checkerboarding, and the constraint parameters $\lambda$, to keep the mass constant,
-# and $\gamma$, to avoid leaving the set $[\chi_{\text{min}}, 1]$. By including gradient regularization, it becomes necessary to calculate the Laplacian.
+#, and $\gamma$, to avoid leaving the set $[\chi_{\text{min}}, 1]$. By including gradient regularization, it becomes necessary to calculate the Laplacian.
 # The Finite Difference Method for square meshes with the edge length $\Delta h$ approximates the Laplacian as follows:
 # ```math
 # \nabla^2 \chi_p = \frac{1}{(\Delta h)^2} (\chi_n + \chi_s + \chi_w + \chi_e - 4 \chi_p)
 # ```
-# Here, the indices refer to the different cardinal directions. Boundary element do not have neighbors in each direction. However, we can calculate
-# the central difference to fulfill Neumann boundary condition. For example, if the element is on the left boundary, we have to fulfill
+# Here, the indices refer to the different cardinal directions. Boundary elements do not have neighbors in each direction. However, we can calculate
+# the central difference to fulfill the Neumann boundary condition. For example, if the element is on the left boundary, we have to fulfill
 # ```math
 # \nabla \chi_p \cdot \textbf{n} = \frac{1}{\Delta h} (\chi_w - \chi_e) = 0
 # ```
-# from which follows $\chi_w = \chi_e$. Thus for boundary elements we can replace the value for the missing neighbor by the value of the opposite neighbor.
-# In order to find the corresponding neighbor elements, we will make use of Ferrites grid topology funcionalities.
+# from which follows $\chi_w = \chi_e$. Thus for boundary elements, we can replace the value for the missing neighbor with the value of the opposite neighbor.
+# In order to find the corresponding neighbor elements, we will make use of Ferrites grid topology functionalities.
 #
 # ## Commented Program
 # We now solve the problem in Ferrite. What follows is a program spliced with comments.
 #md # The full program, without comments, can be found in the next [section](@ref topology_optimization-plain-program).
 #
 # First we load all necessary packages.
-using Ferrite, SparseArrays, LinearAlgebra, Tensors, Printf
+using Ferrite, SparseArrays, LinearAlgebra, Tensors, Printf, ArraysOfArrays
 # Next, we create a simple square grid of the size 2x1. We apply a fixed Dirichlet boundary condition
 # to the left face set, called `clamped`. On the right face, we create a small set `traction`, where we
-# will later apply a force in negative y-direction.
+# will later apply a force in the negative y-direction.
 
 function create_grid(n)
     corners = [Vec{2}((0.0, 0.0)),
@@ -84,7 +84,7 @@ function create_grid(n)
 end
 #md nothing # hide
 
-# Next, we create the FE values, the DofHandler and the Dirichlet boundary condition.
+# Next, we create the FE values, the DofHandler, and the Dirichlet boundary condition.
 
 function create_values()
     ## quadrature rules
@@ -106,9 +106,9 @@ function create_dofhandler(grid)
     return dh
 end
 
-function create_bc(dh)
+function create_bc(dh, grid)
     dbc = ConstraintHandler(dh)
-    add!(dbc, Dirichlet(:u, getnodeset(dh.grid, "clamped"), (x,t) -> zero(Vec{2}), [1,2]))
+    add!(dbc, Dirichlet(:u, getnodeset(grid, "clamped"), (x,t) -> zero(Vec{2}), [1,2]))
     close!(dbc)
     t = 0.0
     update!(dbc, t)
@@ -142,7 +142,7 @@ end
 
 # To store the density and the strain required to calculate the driving forces, we create the struct
 # `MaterialState`. We add a constructor to initialize the struct. The function `update_material_states!`
-# updates the density values once we calculated the new values.
+# updates the density values once we calculate the new values.
 
 mutable struct MaterialState{T, S <: AbstractArray{SymmetricTensor{2, 2, T, 3}, 1}}
     χ::T # density
@@ -185,37 +185,36 @@ function compute_densities(states, dh)
 end
 #md nothing # hide
 
-# For the Laplacian we need some neighboorhood information which is constant throughout the analysis so we compute it once and cache it.
+# For the Laplacian we need some neighborhood information which is constant throughout the analysis so we compute it once and cache it.
 # We iterate through each face of each element,
 # obtaining the neighboring element by using the `getneighborhood` function. For boundary faces,
-# the function call will return an empty object. In that case we use the dictionary to instead find the opposite
+# the function call will return an empty object. In that case, we use the `opp` tuple to find instead the opposite
 # face, as discussed in the introduction.
 
-function cache_neighborhood(dh, topology)
-    nbgs = Vector{Vector{Int}}(undef, getncells(dh.grid))
-    _nfacets = nfacets(dh.grid.cells[1])
-    opp = Dict(1=>3, 2=>4, 3=>1, 4=>2)
-
+function compute_neighborhood(dh, grid, topology)
+    opp = (3, 4, 1, 2)
+    neighborhood = Ferrite.get_facet_facet_neighborhood(topology, grid)
+    nbgs = VectorOfArrays{Int, 1}()
+    nbg = Int[]
     for element in CellIterator(dh)
-        nbg = zeros(Int,_nfacets)
+        _nfacets = nfacets(element)
+        resize!(nbg, _nfacets)
         i = cellid(element)
         for j in 1:_nfacets
-            nbg_cellid = getneighborhood(topology, dh.grid, FacetIndex(i,j))
-            if(!isempty(nbg_cellid))
-                nbg[j] = first(nbg_cellid)[1] # assuming only one face neighbor per cell
+            neighbor_facets = neighborhood[i, j]
+            if(!isempty(neighbor_facets))
+                nbg[j] = neighbor_facets[1][1] # assuming only one facet neighbor per cell
             else # boundary face
-                nbg[j] = first(getneighborhood(topology, dh.grid, FacetIndex(i,opp[j])))[1]
+                nbg[j] = neighborhood[i, opp[j]][1][1]
             end
         end
-
-        nbgs[i] = nbg
+        push!(nbgs, nbg)
     end
-
     return nbgs
 end
 #md nothing # hide
 
-# Now we calculate the Laplacian using the previously cached neighboorhood information.
+# Now we calculate the Laplacian using the previously cached neighborhood information.
 function approximate_laplacian(nbgs, χn, Δh)
     ∇²χ = zeros(length(nbgs))
     for i in 1:length(nbgs)
@@ -287,13 +286,13 @@ end
 # Finally, we put everything together to update the density. The loop ensures the stability of the
 # updated solution.
 
-function update_density(dh, states, mp, ρ,  neighboorhoods, Δh)
+function update_density(dh, states, mp, ρ,  neighborhoods, Δh)
     n_j = Int(ceil(6*mp.β/(mp.η*Δh^2))) # iterations needed for stability
     χn = compute_densities(states, dh) # old density field
     χn1 = zeros(length(χn))
 
     for j in 1:n_j
-        ∇²χ = approximate_laplacian(neighboorhoods, χn, Δh) # Laplacian
+        ∇²χ = approximate_laplacian(neighborhoods, χn, Δh) # Laplacian
         pΨ = compute_driving_forces(states, mp, dh, χn) # driving forces
         p_Ω = compute_average_driving_force(mp, pΨ, χn) # average driving force
 
@@ -387,18 +386,18 @@ end
 #md nothing # hide
 
 # We put everything together in the main function. Here the user may choose the radius parameter, which
-# is related to the regularization parameter as $\beta = ra^2$, the starting density, the number of elements in vertical direction and finally the
+# is related to the regularization parameter as $\beta = ra^2$, the starting density, the number of elements in the vertical direction, and finally the
 # name of the output. Additionally, the user may choose whether only the final design (default)
 # or every iteration step is saved.
 #
-# First, we compute the material parameters and create the grid, DofHandler, boundary condition and FE values.
+# First, we compute the material parameters and create the grid, DofHandler, boundary condition, and FE values.
 # Then we prepare the iterative Newton-Raphson method by pre-allocating all important vectors. Furthermore,
 # we create material states for each element and construct the topology of the grid.
 #
 # During each iteration step, first we solve our FE problem in the Newton-Raphson loop. With the solution of the
-# elastomechanic problem, we check for convergence of our topology design. The criteria has to be fulfilled twice in
+# elastomechanic problem, we check for convergence of our topology design. The criteria have to be fulfilled twice in
 # a row to avoid oscillations. If no convergence is reached yet, we update our design and prepare the next iteration step.
-# Finally, we output the results in paraview and calculate the relative stiffness of the final design, i.e. how much how
+# Finally, we output the results in paraview and calculate the relative stiffness of the final design, i.e. how much
 # the stiffness increased compared to the starting point.
 
 function topopt(ra,ρ,n,filename; output=:false)
@@ -409,7 +408,7 @@ function topopt(ra,ρ,n,filename; output=:false)
     grid = create_grid(n)
     dh = create_dofhandler(grid)
     Δh = 1/n # element edge length
-    dbc = create_bc(dh)
+    dbc = create_bc(dh, grid)
 
     ## cellvalues
     cellvalues, facetvalues = create_values()
@@ -423,9 +422,9 @@ function topopt(ra,ρ,n,filename; output=:false)
     ΔΔu = zeros(n_dofs) # new displacement correction
 
     ## create material states
-    states = [MaterialState(ρ, getnquadpoints(cellvalues)) for _ in 1:getncells(dh.grid)]
+    states = [MaterialState(ρ, getnquadpoints(cellvalues)) for _ in 1:getncells(grid)]
 
-    χ = zeros(getncells(dh.grid))
+    χ = zeros(getncells(grid))
 
     r = zeros(n_dofs) # residual
     K = allocate_matrix(dh) # stiffness matrix
@@ -438,7 +437,7 @@ function topopt(ra,ρ,n,filename; output=:false)
     conv = :false
 
     topology = ExclusiveTopology(grid)
-    neighboorhoods = cache_neighborhood(dh, topology)
+    neighborhoods = compute_neighborhood(dh, grid, topology)
 
     ## Newton-Raphson loop
     NEWTON_TOL = 1e-8
@@ -491,7 +490,7 @@ function topopt(ra,ρ,n,filename; output=:false)
         end
 
         ## update density
-        χ = update_density(dh, states, mp, ρ, neighboorhoods, Δh)
+        χ = update_density(dh, states, mp, ρ, neighborhoods, Δh)
 
         ## update old displacement, density and compliance
         un .= u
@@ -535,7 +534,7 @@ end
 #
 # ![](bending.png)
 #
-# *Figure 2*: Optimization results of the bending beam for smaller (left) and larger (right) value of the regularization parameter $\beta$.
+# *Figure 2*: Optimization results of the bending beam for smaller (left) and larger (right) values of the regularization parameter $\beta$.
 #
 # To prove mesh independence, the user could vary the mesh resolution and compare the results.