diff --git a/docs/Manifest.toml b/docs/Manifest.toml index 3e13ffc023..16dd7f0774 100644 --- a/docs/Manifest.toml +++ b/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" diff --git a/docs/Project.toml b/docs/Project.toml index 876724cc7f..9a6e50861c 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/src/literate-gallery/topology_optimization.jl b/docs/src/literate-gallery/topology_optimization.jl index a03917019f..cfe321a8aa 100644 --- a/docs/src/literate-gallery/topology_optimization.jl +++ b/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.