diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml new file mode 100644 index 0000000000..5b6e58d775 --- /dev/null +++ b/.github/workflows/lint.yml @@ -0,0 +1,13 @@ +name: Lint + +on: + pull_request: + push: + branches: [master] + +jobs: + pre-commit: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: pre-commit/action@2c7b3805fd2a0fd8c1884dcaebf91fc102a13ecd # v3.0.1 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000000..6c219d29c1 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,11 @@ +repos: +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v3.2.0 + hooks: + - id: check-added-large-files + - id: check-case-conflict + - id: check-toml + - id: check-yaml + - id: end-of-file-fixer + - id: mixed-line-ending + - id: trailing-whitespace diff --git a/CHANGELOG.md b/CHANGELOG.md index 1a93d2caa9..2717a6d6be 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -190,8 +190,8 @@ more discussion). + add!(dh, :u, Lagrange{RefTriangle, 1}()) ``` -- **VTK export**: Ferrite no longer extends methods from `WriteVTK.jl`, instead the new types - `VTKFile` and `VTKFileCollection` should be used instead. New methods exists for writing to +- **VTK export**: Ferrite no longer extends methods from `WriteVTK.jl`, instead the new types + `VTKFile` and `VTKFileCollection` should be used instead. New methods exists for writing to a `VTKFile`, e.g. `write_solution`, `write_cell_data`, `write_node_data`, and `write_projection`. See [#692][github-692]. @@ -267,14 +267,14 @@ more discussion). - VTK export now work with `QuadraticHexahedron` elements. ([#714][github-714]) -- The function `bounding_box(::AbstractGrid)` has been added. It computes the bounding box for - a given grid (based on its node coordinates), and returns the minimum and maximum vertices +- The function `bounding_box(::AbstractGrid)` has been added. It computes the bounding box for + a given grid (based on its node coordinates), and returns the minimum and maximum vertices of the bounding box. ([#880][github-880]) ### Changed - `create_sparsity_pattern` now supports cross-element dof coupling by passing kwarg - `topology` along with an optional `cross_coupling` matrix that behaves similar to + `topology` along with an optional `cross_coupling` matrix that behaves similar to the `coupling` kwarg. ([#710][github-#710]) - The `AbstractCell` interface has been reworked. This change should not affect user code, diff --git a/benchmark/benchmarks-assembly.jl b/benchmark/benchmarks-assembly.jl index 16293db85f..dddcd11530 100644 --- a/benchmark/benchmarks-assembly.jl +++ b/benchmark/benchmarks-assembly.jl @@ -67,7 +67,7 @@ for spatial_dim ∈ 1:3 LAGRANGE_SUITE["ritz-galerkin"]["face-flux"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_local_matrix($grid, $fsv, shape_gradient, shape_value, *) LAGRANGE_SUITE["petrov-galerkin"]["face-flux"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_local_matrix($grid, $fsv, shape_gradient, $fsv2, shape_value, *) - + ip = DiscontinuousLagrange{ref_type, order}() isv = InterfaceValues(qr_facet, ip, ip_geo); isv2 = InterfaceValues(qr_facet, ip, ip_geo); @@ -77,10 +77,10 @@ for spatial_dim ∈ 1:3 LAGRANGE_SUITE["ritz-galerkin"]["interface-{grad}⋅[[val]]"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_interfaces($dh, $isv, shape_gradient_average, shape_value_jump, ⋅) LAGRANGE_SUITE["petrov-galerkin"]["interface-{grad}⋅[[val]]"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_interfaces($dh, $isv, shape_gradient_average, $isv2, shape_value_jump, ⋅) - + LAGRANGE_SUITE["ritz-galerkin"]["interface-interior-penalty"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_interfaces($dh, $isv, shape_value_jump, shape_value_jump, ⋅) LAGRANGE_SUITE["petrov-galerkin"]["interface-interior-penalty"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_interfaces($dh, $isv, shape_value_jump, $isv2, shape_value_jump, ⋅) - + end end end diff --git a/docs/old_examples/README.md b/docs/old_examples/README.md index fec8fe738d..2dd9293b2b 100644 --- a/docs/old_examples/README.md +++ b/docs/old_examples/README.md @@ -1,2 +1,2 @@ The examples in this folder might be out of date. -For up to date examples, see the documentation. \ No newline at end of file +For up to date examples, see the documentation. diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index efbb3fd1c5..a8333e9665 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -64,7 +64,7 @@ @phdthesis{Cenanovic2017 year={2017}, } @misc{Kirby2017, - title={A general approach to transforming finite elements}, + title={A general approach to transforming finite elements}, author={Robert C. Kirby}, year={2017}, eprint={1706.09017}, diff --git a/docs/src/devdocs/FEValues.md b/docs/src/devdocs/FEValues.md index faacad82fa..8dc820187c 100644 --- a/docs/src/devdocs/FEValues.md +++ b/docs/src/devdocs/FEValues.md @@ -8,7 +8,7 @@ * [`FacetValues`](@ref) * [`BCValues`](@ref) * [`PointValues`](@ref) - + ## Internal types ```@docs @@ -20,7 +20,7 @@ Ferrite.BCValues ## Custom FEValues Custom FEValues, `fe_v::AbstractValues`, should normally implement the [`reinit!`](@ref) method. Subtypes of `AbstractValues` have default implementations for some functions, but require some lower-level access functions, specifically - + * [`function_value`](@ref), requires * [`shape_value`](@ref) * [`getnquadpoints`](@ref) @@ -33,7 +33,7 @@ Custom FEValues, `fe_v::AbstractValues`, should normally implement the [`reinit! * [`geometric_value`](@ref) * `getngeobasefunctions` * [`getnquadpoints`](@ref) - + ### Array bounds * Asking for the `n`th quadrature point must be inside array bounds if `1 <= n <= getnquadpoints(fe_v)`. (`checkquadpoint` can, alternatively, be dispatched to check that `n` is inbounds.) diff --git a/docs/src/devdocs/interpolations.md b/docs/src/devdocs/interpolations.md index c282d753a0..f8a58ad1d3 100644 --- a/docs/src/devdocs/interpolations.md +++ b/docs/src/devdocs/interpolations.md @@ -43,7 +43,7 @@ Ferrite.mapping_type ``` for all entities which exist on that reference element. The dof functions default to having no -dofs defined on a specific entity. Hence, not overloading of the dof functions will result in an -element with zero dofs. Also, it should always be double checked that everything is consistent as +dofs defined on a specific entity. Hence, not overloading of the dof functions will result in an +element with zero dofs. Also, it should always be double checked that everything is consistent as specified in the docstring of the corresponding function, as inconsistent implementations can lead to bugs which are really difficult to track down. diff --git a/docs/src/literate-gallery/helmholtz.jl b/docs/src/literate-gallery/helmholtz.jl index 8e01c44e70..17afe674cd 100644 --- a/docs/src/literate-gallery/helmholtz.jl +++ b/docs/src/literate-gallery/helmholtz.jl @@ -15,7 +15,7 @@ # ```math # n \cdot \nabla u = g_2 \quad x \in \Gamma_2 # ``` -# +# # Here Γ₁ is the union of the top and the right boundary of the square, # while Γ₂ is the union of the bottom and the left boundary. # @@ -150,7 +150,7 @@ function doassemble(cellvalues::CellValues, facetvalues::FacetValues, end end end - + celldofs!(global_dofs, cell) assemble!(assembler, global_dofs, fe, Ke) end diff --git a/docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl b/docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl index 40e699a3cf..6183808c1c 100644 --- a/docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl +++ b/docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl @@ -334,7 +334,7 @@ function solve(interpolation_u, interpolation_p) end; ## Save the solution fields - addstep!(pvd, t) do io + addstep!(pvd, t) do io write_solution(io, dh, w) end end; diff --git a/docs/src/literate-gallery/topology_optimization.jl b/docs/src/literate-gallery/topology_optimization.jl index b9e12f825e..894d8f067d 100644 --- a/docs/src/literate-gallery/topology_optimization.jl +++ b/docs/src/literate-gallery/topology_optimization.jl @@ -14,10 +14,10 @@ # # ## Introduction # -# Topology optimization is the task of finding structures that are mechanically ideal. +# 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 # 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 +# 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 # [BlaJanJun2022taylorwlsthermotopopt](@cite). # @@ -26,17 +26,17 @@ # means bulk material. Then, we use a SIMP ansatz (solid isotropic material with penalization) for the stiffness tensor # $C(\chi) = \chi^p C_0$, where $C_0$ is the stiffness of the bulk material. The SIMP exponent $p>1$ ensures that the # model prefers the density values void and bulk before the intermediate values. The variational formulation then yields -# the modified Gibbs energy +# the modified Gibbs energy # ```math # G = \int_{\Omega} \frac{1}{2} \chi^p \varepsilon : C : \varepsilon \; \text{d}V - \int_{\Omega} \boldsymbol{f} \cdot \boldsymbol{u} \; \text{d}V - \int_{\partial\Omega} \boldsymbol{t} \cdot \boldsymbol{u} \; \text{d}A. # ``` # Furthermore, we receive the evolution equation of the density # and the additional Neumann boundary condition in the strong form # ```math -# p_\chi + \eta \dot{\chi} + \lambda + \gamma - \beta \nabla^2 \chi \ni 0 \quad \forall \textbf{x} \in \Omega, +# p_\chi + \eta \dot{\chi} + \lambda + \gamma - \beta \nabla^2 \chi \ni 0 \quad \forall \textbf{x} \in \Omega, # ``` # ```math -# \beta \nabla \chi \cdot \textbf{n} = 0 \quad \forall \textbf{x} \in \partial \Omega, +# \beta \nabla \chi \cdot \textbf{n} = 0 \quad \forall \textbf{x} \in \partial \Omega, # ``` # with the thermodynamic driving force # ```math @@ -57,8 +57,8 @@ # ```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 by the value of the opposite neighbor. +# In order to find the corresponding neighbor elements, we will make use of Ferrites grid topology funcionalities. # # ## Commented Program # We now solve the problem in Ferrite. What follows is a program spliced with comments. @@ -70,16 +70,16 @@ using Ferrite, SparseArrays, LinearAlgebra, Tensors, Printf # 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. -function create_grid(n) +function create_grid(n) corners = [Vec{2}((0.0, 0.0)), Vec{2}((2.0, 0.0)), Vec{2}((2.0, 1.0)), Vec{2}((0.0, 1.0))] grid = generate_grid(Quadrilateral, (2*n, n), corners); - + ## node-/facesets for boundary conditions addnodeset!(grid, "clamped", x -> x[1] ≈ 0.0) - addfacetset!(grid, "traction", x -> x[1] ≈ 2.0 && norm(x[2]-0.5) <= 0.05); + addfacetset!(grid, "traction", x -> x[1] ≈ 2.0 && norm(x[2]-0.5) <= 0.05); return grid end #md nothing # hide @@ -95,7 +95,7 @@ function create_values() ip = Lagrange{RefQuadrilateral,1}()^2 cellvalues = CellValues(qr, ip) facetvalues = FacetValues(facet_qr, ip) - + return cellvalues, facetvalues end @@ -106,7 +106,7 @@ function create_dofhandler(grid) return dh end -function create_bc(dh) +function create_bc(dh) dbc = ConstraintHandler(dh) add!(dbc, Dirichlet(:u, getnodeset(dh.grid, "clamped"), (x,t) -> zero(Vec{2}), [1,2])) close!(dbc) @@ -122,14 +122,14 @@ end struct MaterialParameters{T, S <: SymmetricTensor{4, 2, T}} C::S - χ_min::T + χ_min::T p::T β::T η::T end #md nothing # hide -function MaterialParameters(E, ν, χ_min, p, β, η) +function MaterialParameters(E, ν, χ_min, p, β, η) δ(i,j) = i == j ? 1.0 : 0.0 # helper function G = E / 2(1 + ν) # =μ @@ -144,7 +144,7 @@ end # `MaterialState`. We add a constructor to initialize the struct. The function `update_material_states!` # updates the density values once we calculated the new values. -mutable struct MaterialState{T, S <: AbstractArray{SymmetricTensor{2, 2, T}, 1}} +mutable struct MaterialState{T, S <: AbstractArray{SymmetricTensor{2, 2, T}, 1}} χ::T # density ε::S # strain in each quadrature point end @@ -163,7 +163,7 @@ end # Next, we define a function to calculate the driving forces for all elements. # For this purpose, we iterate through all elements and calculate the average strain in each # element. Then, we compute the driving force from the formula introduced at the beginning. -# We create a second function to collect the density in each element. +# We create a second function to collect the density in each element. function compute_driving_forces(states, mp, dh, χn) pΨ = zeros(length(states)) @@ -185,18 +185,18 @@ function compute_densities(states, dh) end #md nothing # hide -# Now we calculate the Laplacian. For this purpose, we will later create the grid topology of +# Now we calculate the Laplacian. For this purpose, we will later create the grid topology of # the grid by using the function `ExclusiveTopology`. Then 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 -# face, as discussed in the introduction. Then, the approximation of the Laplacian reduces to the sum below. +# face, as discussed in the introduction. Then, the approximation of the Laplacian reduces to the sum below. function approximate_laplacian(dh, topology, χn, Δh) ∇²χ = zeros(getncells(dh.grid)) _nfacets = nfacets(dh.grid.cells[1]) opp = Dict(1=>3, 2=>4, 3=>1, 4=>2) nbg = zeros(Int,_nfacets) - + for element in CellIterator(dh) i = cellid(element) for j in 1:_nfacets @@ -207,7 +207,7 @@ function approximate_laplacian(dh, topology, χn, Δh) nbg[j] = first(getcells(getneighborhood(topology, dh.grid, FacetIndex(i,opp[j])))) end end - + ∇²χ[i] = (χn[nbg[1]]+χn[nbg[2]]+χn[nbg[3]]+χn[nbg[4]]-4*χn[i])/(Δh^2) end @@ -218,30 +218,30 @@ end # For the iterative computation of the solution, a function is needed to update the densities in each element. # To ensure that the mass is kept constant, we have to calculate the constraint # parameter $\lambda$, which we do via the bisection method. We repeat the calculation -# until the difference between the average density (calculated from the element-wise trial densities) and the target density nearly vanishes. +# until the difference between the average density (calculated from the element-wise trial densities) and the target density nearly vanishes. # By using the extremal values of $\Delta \chi$ as the starting interval, we guarantee that the method converges eventually. -function compute_χn1(χn, Δχ, ρ, ηs, χ_min) +function compute_χn1(χn, Δχ, ρ, ηs, χ_min) n_el = length(χn) - + χ_trial = zeros(n_el) ρ_trial = 0.0 - + λ_lower = minimum(Δχ) - ηs λ_upper = maximum(Δχ) + ηs λ_trial = 0.0 - + while(abs(ρ-ρ_trial)>1e-7) for i in 1:n_el Δχt = 1/ηs * (Δχ[i] - λ_trial) χ_trial[i] = maximum([χ_min, minimum([1.0, χn[i]+Δχt])]) end - + ρ_trial = 0.0 for i in 1:n_el - ρ_trial += χ_trial[i]/n_el + ρ_trial += χ_trial[i]/n_el end - + if(ρ_trial > ρ) λ_lower = λ_trial elseif(ρ_trial < ρ) @@ -256,65 +256,65 @@ end # Lastly, we use the following helper function to compute the average driving force, which is later # used to normalize the driving forces. This makes the used material parameters and numerical parameters independent -# of the problem. +# of the problem. function compute_average_driving_force(mp, pΨ, χn) n = length(pΨ) w = zeros(n) - + for i in 1:n w[i] = (χn[i]-mp.χ_min)*(1-χn[i]) end - + p_Ω = sum(w.*pΨ)/sum(w) # average driving force - + return p_Ω end #md nothing # hide -# Finally, we put everything together to update the density. The loop ensures the stability of the +# Finally, we put everything together to update the density. The loop ensures the stability of the # updated solution. function update_density(dh, states, mp, ρ, topology, Δh) n_j = Int(ceil(6*mp.β/(mp.η*Δh^2))) # iterations needed for stability - χn = compute_densities(states, dh) # old density field + χn = compute_densities(states, dh) # old density field χn1 = zeros(length(χn)) - + for j in 1:n_j ∇²χ = approximate_laplacian(dh, topology, χn, Δh) # Laplacian pΨ = compute_driving_forces(states, mp, dh, χn) # driving forces p_Ω = compute_average_driving_force(mp, pΨ, χn) # average driving force - - Δχ = pΨ/p_Ω + mp.β*∇²χ - χn1 = compute_χn1(χn, Δχ, ρ, mp.η, mp.χ_min) + Δχ = pΨ/p_Ω + mp.β*∇²χ + + χn1 = compute_χn1(χn, Δχ, ρ, mp.η, mp.χ_min) if(j 10 error("Reached maximum Newton iterations, aborting") break end - + ## current guess u .= un .+ Δu K, r = doassemble!(cellvalues, facetvalues, K, grid, dh, mp, u, states); - norm_r = norm(r[Ferrite.free_dofs(dbc)]) + norm_r = norm(r[Ferrite.free_dofs(dbc)]) if (norm_r) < NEWTON_TOL break - end + end apply_zero!(K, r, dbc) ΔΔu = Symmetric(K) \ r - + apply_zero!(ΔΔu, dbc) Δu .+= ΔΔu - end # of loop while NR-Iteration + end # of loop while NR-Iteration ## calculate compliance compliance = 1/2 * u' * K * u - + if(it==1) compliance_0 = compliance end - + ## check convergence criterium (twice!) if(abs(compliance-compliance_n)/compliance < tol) if(conv) @@ -477,16 +477,16 @@ function topopt(ra,ρ,n,filename; output=:false) else conv = :false end - + ## update density χ = update_density(dh, states, mp, ρ, topology, Δh) - + ## update old displacement, density and compliance un .= u Δu .= 0.0 update_material_states!(χ, states, dh) compliance_n = compliance - + ## output during calculation if(output) i = @sprintf("%3.3i", it) @@ -505,13 +505,13 @@ function topopt(ra,ρ,n,filename; output=:false) end end @printf "Rel. stiffness: %.4f \n" compliance^(-1)/compliance_0^(-1) - + return end #md nothing # hide # Lastly, we call our main function and compare the results. To create the -# complete output with all iteration steps, it is possible to set the output +# complete output with all iteration steps, it is possible to set the output # parameter to `true`. topopt(0.02, 0.5, 60, "small_radius"; output=:false); diff --git a/docs/src/literate-tutorials/dg_heat_equation.jl b/docs/src/literate-tutorials/dg_heat_equation.jl index 2e9d39f692..f5c2495837 100644 --- a/docs/src/literate-tutorials/dg_heat_equation.jl +++ b/docs/src/literate-tutorials/dg_heat_equation.jl @@ -71,7 +71,7 @@ # ```math # \int_\Gamma q \boldsymbol{\phi} \cdot \boldsymbol{n} \,\mathrm{d}\Gamma = \int_\Gamma \llbracket q\rrbracket \cdot \{\boldsymbol{\phi}\} \,\mathrm{d}\Gamma + \int_{\Gamma^0} \{q\} \llbracket \boldsymbol{\phi}\rrbracket \,\mathrm{d}\Gamma^0, # ``` -# where $\Gamma^0 : \Gamma \setminus \partial \Omega$, and the jump of the vector-valued field $\boldsymbol{\phi}$ is defined as +# where $\Gamma^0 : \Gamma \setminus \partial \Omega$, and the jump of the vector-valued field $\boldsymbol{\phi}$ is defined as # ```math # \llbracket \boldsymbol{\phi}\rrbracket = \boldsymbol{\phi}^+ \cdot \boldsymbol{n}^+ + \boldsymbol{\phi}^- \cdot \boldsymbol{n}^-\\ # ``` @@ -99,7 +99,7 @@ # ```math # \int_\Omega [\boldsymbol{\nabla} (u)] \cdot [\boldsymbol{\nabla} (\delta u)] \,\mathrm{d}\Omega + \int_\Gamma \llbracket \hat{u} - u\rrbracket \cdot \{\boldsymbol{\nabla} (\delta u)\} \,\mathrm{d}\Gamma + \int_{\Gamma^0} \{\hat{u} - u\} \llbracket \boldsymbol{\nabla} (\delta u)\rrbracket \,\mathrm{d}\Gamma^0 - \int_\Gamma \llbracket \delta u\rrbracket \cdot \{\hat{\boldsymbol{\sigma}}\} \,\mathrm{d}\Gamma - \int_{\Gamma^0} \{\delta u\} \llbracket \hat{\boldsymbol{\sigma}}\rrbracket \,\mathrm{d}\Gamma^0 = \int_\Omega \delta u \,\mathrm{d}\Omega,\\ # ``` -# The numerical fluxes chosen for the interior penalty method are $\boldsymbol{\hat{\sigma}} = \{\boldsymbol{\nabla} (u)\} - \alpha(\llbracket u\rrbracket)$ on $\Gamma$, $\hat{u} = \{u\}$ on the interfaces between elements $\Gamma^0 : \Gamma \setminus \partial \Omega$, +# The numerical fluxes chosen for the interior penalty method are $\boldsymbol{\hat{\sigma}} = \{\boldsymbol{\nabla} (u)\} - \alpha(\llbracket u\rrbracket)$ on $\Gamma$, $\hat{u} = \{u\}$ on the interfaces between elements $\Gamma^0 : \Gamma \setminus \partial \Omega$, # and $\hat{u} = 0$ on $\partial \Omega$. Such choice results in $\{\hat{\boldsymbol{\sigma}}\} = \{\boldsymbol{\nabla} (u)\} - \alpha(\llbracket u\rrbracket)$, $\llbracket \hat{u}\rrbracket = 0$, $\{\hat{u}\} = \{u\}$, $\llbracket \hat{\boldsymbol{\sigma}}\rrbracket = 0$ and the equation becomes # ```math # \int_\Omega [\boldsymbol{\nabla} (u)] \cdot [\boldsymbol{\nabla} (\delta u)] \,\mathrm{d}\Omega - \int_\Gamma \llbracket u\rrbracket \cdot \{\boldsymbol{\nabla} (\delta u)\} \,\mathrm{d}\Gamma - \int_\Gamma \llbracket \delta u\rrbracket \cdot \{\boldsymbol{\nabla} (u)\} - \llbracket \delta u\rrbracket \cdot \alpha(\llbracket u\rrbracket) \,\mathrm{d}\Gamma = \int_\Omega \delta u \,\mathrm{d}\Omega,\\ @@ -170,7 +170,7 @@ close!(dh); K = create_sparsity_pattern(dh, topology = topology, cross_coupling = trues(1,1)); # ### Boundary conditions -# The Dirichlet boundary conditions are treated +# The Dirichlet boundary conditions are treated # as usual by a `ConstraintHandler`. ch = ConstraintHandler(dh) add!(ch, Dirichlet(:u, getfacetset(grid, "right"), (x, t) -> 1.0)) diff --git a/docs/src/literate-tutorials/linear_shell.jl b/docs/src/literate-tutorials/linear_shell.jl index 54ef96bbb9..cbddd7b9b1 100644 --- a/docs/src/literate-tutorials/linear_shell.jl +++ b/docs/src/literate-tutorials/linear_shell.jl @@ -4,8 +4,8 @@ #- # ## Introduction # -# In this example we show how shell elements can be analyzed in Ferrite.jl. The shell implemented here comes from the book -# "The finite element method - Linear static and dynamic finite element analysis" by Hughes (1987), and a brief description of it is +# In this example we show how shell elements can be analyzed in Ferrite.jl. The shell implemented here comes from the book +# "The finite element method - Linear static and dynamic finite element analysis" by Hughes (1987), and a brief description of it is # given at the end of this tutorial. The first part of the tutorial explains how to set up the problem. # ## Setting up the problem @@ -15,7 +15,7 @@ using ForwardDiff function main() #wrap everything in a function... # First we generate a flat rectangular mesh. There is currently no built-in function for generating -# shell meshes in Ferrite, so we have to create our own simple mesh generator (see the +# shell meshes in Ferrite, so we have to create our own simple mesh generator (see the # function `generate_shell_grid` further down in this file). #+ nels = (10,10) @@ -23,8 +23,8 @@ size = (10.0, 10.0) grid = generate_shell_grid(nels, size) # Here we define the bi-linear interpolation used for the geometrical description of the shell. -# We also create two quadrature rules for the in-plane and out-of-plane directions. Note that we use -# under integration for the inplane integration, to avoid shear locking. +# We also create two quadrature rules for the in-plane and out-of-plane directions. Note that we use +# under integration for the inplane integration, to avoid shear locking. #+ ip = Lagrange{RefQuadrilateral,1}() qr_inplane = QuadratureRule{RefQuadrilateral}(1) @@ -38,8 +38,8 @@ add!(dh, :u, ip^3) add!(dh, :θ, ip^2) close!(dh) -# In order to apply our boundary conditions, we first need to create some facet- and vertex-sets. This -# is done with `addfacetset!` and `addvertexset!` +# In order to apply our boundary conditions, we first need to create some facet- and vertex-sets. This +# is done with `addfacetset!` and `addvertexset!` #+ addfacetset!(grid, "left", (x) -> x[1] ≈ 0.0) addfacetset!(grid, "right", (x) -> x[1] ≈ size[1]) @@ -63,7 +63,7 @@ add!(ch, Dirichlet(:θ, getvertexset(grid, "corner"), (x, t) -> (0.0), [2]) ) close!(ch) update!(ch, 0.0) -# Next we define relevant data for the shell, such as shear correction factor and stiffness matrix for the material. +# Next we define relevant data for the shell, such as shear correction factor and stiffness matrix for the material. # In this linear shell, plane stress is assumed, ie $\\sigma_{zz} = 0$. Therefor, the stiffness matrix is 5x5 (opposed to the normal 6x6). #+ κ = 5/6 # Shear correction factor @@ -99,7 +99,7 @@ for cell in CellIterator(grid) reinit!(cv, cell) celldofs!(celldofs, dh, cellid(cell)) getcoordinates!(cellcoords, grid, cellid(cell)) - + #Call the element routine integrate_shell!(ke, cv, qr_ooplane, cellcoords, data) @@ -120,7 +120,7 @@ end end; #end main functions # Below is the function that creates the shell mesh. It simply generates a 2d-quadrature mesh, and appends -# a third coordinate (z-direction) to the node-positions. +# a third coordinate (z-direction) to the node-positions. function generate_shell_grid(nels, size) _grid = generate_grid(Quadrilateral, nels, Vec((0.0,0.0)), Vec(size)) nodes = [(n.x[1], n.x[2], 0.0) |> Vec{3} |> Node for n in _grid.nodes] @@ -139,7 +139,7 @@ end; #md # !!! note #md # This element might experience various locking phenomenas, and should only be seen as a proof of concept. -# ##### Fiber coordinate system +# ##### Fiber coordinate system # The element uses two coordinate systems. The first coordianate system, called the fiber system, is created for each # element node, and is used as a reference frame for the rotations. The function below implements an algorithm that return the # fiber directions, $\boldsymbol{e}^{f}_{a1}$, $\boldsymbol{e}^{f}_{a2}$ and $\boldsymbol{e}^{f}_{a3}$, at each node $a$. @@ -153,7 +153,7 @@ function fiber_coordsys(Ps::Vector{Vec{3,Float64}}) j = 1 if a[1] > a[3]; a[3] = a[1]; j = 2; end if a[2] > a[3]; j = 3; end - + e3 = P e2 = Tensors.cross(P, basevec(Vec{3}, j)) e2 /= norm(e2) @@ -167,7 +167,7 @@ function fiber_coordsys(Ps::Vector{Vec{3,Float64}}) end; -# ##### Lamina coordinate system +# ##### Lamina coordinate system # The second coordinate system is the so called Lamina Coordinate system. It is # created for each integration point, and is defined to be tangent to the # mid-surface. It is in this system that we enforce that plane stress assumption, @@ -205,7 +205,7 @@ end; # A material point in the shell is defined as # ```math # \boldsymbol x(\xi, \eta, \zeta) = \sum_{a=1}^{N_{\text{nodes}}} N_a(\xi, \eta) \boldsymbol{\bar{x}}_{a} + ζ \frac{h}{2} \boldsymbol{\bar{p}_a} -# ``` +# ``` # where $\boldsymbol{\bar{x}}_{a}$ are nodal positions on the mid-surface, and $\boldsymbol{\bar{p}_a}$ is an vector that defines the fiber direction # on the reference surface. $N_a$ arethe shape functions. # @@ -235,13 +235,13 @@ end; # ``` # The displacement field is calculated as: # ```math -# \boldsymbol u = \sum_{a=1}^{N_{\text{nodes}}} N_a \bar{\boldsymbol u}_{a} + +# \boldsymbol u = \sum_{a=1}^{N_{\text{nodes}}} N_a \bar{\boldsymbol u}_{a} + # N_a ζ\frac{h}{2}(\theta_{a2} \boldsymbol e^{f}_{a1} - \theta_{a1} \boldsymbol e^{f}_{a2}) # # ``` # The gradient of the displacement (in the lamina coordinate system), then becomes: # ```math -# \frac{\partial u_{i}}{\partial x_j} = \sum_{m=1}^3 q_{im} \sum_{a=1}^{N_{\text{nodes}}} \frac{\partial N_a}{\partial x_j} \bar{u}_{am} + +# \frac{\partial u_{i}}{\partial x_j} = \sum_{m=1}^3 q_{im} \sum_{a=1}^{N_{\text{nodes}}} \frac{\partial N_a}{\partial x_j} \bar{u}_{am} + # \frac{\partial(N_a ζ)}{\partial x_j} \frac{h}{2} (\theta_{a2} e^{f}_{am1} - \theta_{a1} e^{f}_{am2}) # ``` function strain(dofvec::Vector{T}, N, dNdx, ζ, dζdx, q, ef1, ef2, h) where T @@ -262,7 +262,7 @@ function strain(dofvec::Vector{T}, N, dNdx, ζ, dζdx, q, ef1, ef2, h) where T end; # ##### Main element routine -# Below is the main routine that calculates the stiffness matrix of the shell element. +# Below is the main routine that calculates the stiffness matrix of the shell element. # Since it is a so called degenerate shell element, the code is similar to that for an standard continuum element. shape_reference_gradient(cv::CellValues, q_point, i) = cv.fun_values.dNdξ[i, q_point] diff --git a/docs/src/literate-tutorials/porous_media.jl b/docs/src/literate-tutorials/porous_media.jl index 6c8f6b78c1..81a807afbc 100644 --- a/docs/src/literate-tutorials/porous_media.jl +++ b/docs/src/literate-tutorials/porous_media.jl @@ -1,12 +1,12 @@ -# # Porous media +# # Porous media # Porous media is a two-phase material, consisting of solid parts and a liquid occupying -# the pores inbetween. -# Using the porous media theory, we can model such a material without explicitly -# resolving the microstructure, but by considering the interactions between the -# solid and liquid. In this example, we will additionally consider larger linear -# elastic solid aggregates that are impermeable. Hence, there is no liquids in -# these particles and the only unknown variable is the displacement field `:u`. +# the pores inbetween. +# Using the porous media theory, we can model such a material without explicitly +# resolving the microstructure, but by considering the interactions between the +# solid and liquid. In this example, we will additionally consider larger linear +# elastic solid aggregates that are impermeable. Hence, there is no liquids in +# these particles and the only unknown variable is the displacement field `:u`. # In the porous media, denoted the matrix, we have both the displacement field, # `:u`, as well as the liquid pressure, `:p`, as unknown. The simulation result # is shown below @@ -21,62 +21,62 @@ # \dot{\Phi}(\boldsymbol{\epsilon}, p) + \boldsymbol{w}(p) \cdot \boldsymbol{\nabla} &= 0 # \end{aligned} # ``` -# where -# ``\boldsymbol{\epsilon} = \left[\boldsymbol{u}\otimes\boldsymbol{\nabla}\right]^\mathrm{sym}`` -# The constitutive relationships are +# where +# ``\boldsymbol{\epsilon} = \left[\boldsymbol{u}\otimes\boldsymbol{\nabla}\right]^\mathrm{sym}`` +# The constitutive relationships are # ```math # \begin{aligned} # \boldsymbol{\sigma} &= \boldsymbol{\mathsf{C}}:\boldsymbol{\epsilon} - \alpha p \boldsymbol{I} \\ # \boldsymbol{w} &= - k \boldsymbol{\nabla} p \\ # \Phi &= \phi + \alpha \mathrm{tr}(\boldsymbol{\epsilon}) + \beta p # \end{aligned} -# ``` -# with +# ``` +# with # ``\boldsymbol{\mathsf{C}}=2G \boldsymbol{\mathsf{I}}^\mathrm{dev} + 3K \boldsymbol{I}\otimes\boldsymbol{I}``. -# The material parameters are then the -# shear modulus, ``G``, -# bulk modulus, ``K``, -# permeability, ``k``, +# The material parameters are then the +# shear modulus, ``G``, +# bulk modulus, ``K``, +# permeability, ``k``, # Biot's coefficient, ``\alpha``, and # liquid compressibility, ``\beta``. -# The porosity, ``\phi``, doesn't enter into the equations +# The porosity, ``\phi``, doesn't enter into the equations # (A different porosity leads to different skeleton stiffness and permeability). # -# +# # The variational (weak) form can then be derived for the variations ``\boldsymbol{\delta u}`` # and ``\delta p`` as # ```math # \begin{aligned} # \int_\Omega \left[\left[\boldsymbol{\delta u}\otimes\boldsymbol{\nabla}\right]^\mathrm{sym}: -# \boldsymbol{\mathsf{C}}:\boldsymbol{\epsilon} - \boldsymbol{\delta u} \cdot \boldsymbol{\nabla} \alpha p\right] \mathrm{d}\Omega +# \boldsymbol{\mathsf{C}}:\boldsymbol{\epsilon} - \boldsymbol{\delta u} \cdot \boldsymbol{\nabla} \alpha p\right] \mathrm{d}\Omega # &= \int_\Gamma \boldsymbol{\delta u} \cdot \boldsymbol{t} \mathrm{d} \Gamma \\ -# \int_\Omega \left[\delta p \left[\alpha \dot{\boldsymbol{u}} \cdot \boldsymbol{\nabla} + \beta \dot{p}\right] + -# \boldsymbol{\nabla}(\delta p) \cdot [k \boldsymbol{\nabla}]\right] \mathrm{d}\Omega -# &= \int_\Gamma \delta p w_\mathrm{n} \mathrm{d} \Gamma +# \int_\Omega \left[\delta p \left[\alpha \dot{\boldsymbol{u}} \cdot \boldsymbol{\nabla} + \beta \dot{p}\right] + +# \boldsymbol{\nabla}(\delta p) \cdot [k \boldsymbol{\nabla}]\right] \mathrm{d}\Omega +# &= \int_\Gamma \delta p w_\mathrm{n} \mathrm{d} \Gamma # \end{aligned} # ``` -# where ``\boldsymbol{t}=\boldsymbol{n}\cdot\boldsymbol{\sigma}`` is the traction and -# ``w_\mathrm{n} = \boldsymbol{n}\cdot\boldsymbol{w}`` is the normal flux. -# +# where ``\boldsymbol{t}=\boldsymbol{n}\cdot\boldsymbol{\sigma}`` is the traction and +# ``w_\mathrm{n} = \boldsymbol{n}\cdot\boldsymbol{w}`` is the normal flux. +# # ### Finite element form -# Discretizing in space using finite elements, we obtain the vector equation -# ``r_i = f_i^\mathrm{int} - f_{i}^\mathrm{ext}`` where ``f^\mathrm{ext}`` are the external -# "forces", and ``f_i^\mathrm{int}`` are the internal "forces". We split this into the -# displacement part ``r_i^\mathrm{u} = f_i^\mathrm{int,u} - f_{i}^\mathrm{ext,u}`` and +# Discretizing in space using finite elements, we obtain the vector equation +# ``r_i = f_i^\mathrm{int} - f_{i}^\mathrm{ext}`` where ``f^\mathrm{ext}`` are the external +# "forces", and ``f_i^\mathrm{int}`` are the internal "forces". We split this into the +# displacement part ``r_i^\mathrm{u} = f_i^\mathrm{int,u} - f_{i}^\mathrm{ext,u}`` and # pressure part ``r_i^\mathrm{p} = f_i^\mathrm{int,p} - f_{i}^\mathrm{ext,p}`` -# to obtain the discretized equation system +# to obtain the discretized equation system # ```math # \begin{aligned} -# f_i^\mathrm{int,u} &= \int_\Omega [\boldsymbol{\delta N}^\mathrm{u}_i\otimes\boldsymbol{\nabla}]^\mathrm{sym} : \boldsymbol{\mathsf{C}} : [\boldsymbol{u}\otimes\boldsymbol{\nabla}]^\mathrm{sym} \ -# - [\boldsymbol{\delta N}^\mathrm{u}_i \cdot \boldsymbol{\nabla}] \alpha p \mathrm{d}\Omega +# f_i^\mathrm{int,u} &= \int_\Omega [\boldsymbol{\delta N}^\mathrm{u}_i\otimes\boldsymbol{\nabla}]^\mathrm{sym} : \boldsymbol{\mathsf{C}} : [\boldsymbol{u}\otimes\boldsymbol{\nabla}]^\mathrm{sym} \ +# - [\boldsymbol{\delta N}^\mathrm{u}_i \cdot \boldsymbol{\nabla}] \alpha p \mathrm{d}\Omega # &= \int_\Gamma \boldsymbol{\delta N}^\mathrm{u}_i \cdot \boldsymbol{t} \mathrm{d} \Gamma \\ -# f_i^\mathrm{int,p} &= \int_\Omega \delta N_i^\mathrm{p} [\alpha [\dot{\boldsymbol{u}}\cdot\boldsymbol{\nabla}] + \beta\dot{p}] + \boldsymbol{\nabla}(\delta N_i^\mathrm{p}) \cdot [k \boldsymbol{\nabla}(p)] \mathrm{d}\Omega +# f_i^\mathrm{int,p} &= \int_\Omega \delta N_i^\mathrm{p} [\alpha [\dot{\boldsymbol{u}}\cdot\boldsymbol{\nabla}] + \beta\dot{p}] + \boldsymbol{\nabla}(\delta N_i^\mathrm{p}) \cdot [k \boldsymbol{\nabla}(p)] \mathrm{d}\Omega # &= \int_\Gamma \delta N_i^\mathrm{p} w_\mathrm{n} \mathrm{d} \Gamma # \end{aligned} # ``` # Approximating the time-derivatives, ``\dot{\boldsymbol{u}}\approx \left[\boldsymbol{u}-{}^n\boldsymbol{u}\right]/\Delta t`` -# and ``\dot{p}\approx \left[p-{}^np\right]/\Delta t``, we can implement the finite element equations in the residual form -# ``r_i(\boldsymbol{a}(t), t) = 0`` where the vector ``\boldsymbol{a}`` contains all unknown displacements ``u_i`` and pressures ``p_i``. +# and ``\dot{p}\approx \left[p-{}^np\right]/\Delta t``, we can implement the finite element equations in the residual form +# ``r_i(\boldsymbol{a}(t), t) = 0`` where the vector ``\boldsymbol{a}`` contains all unknown displacements ``u_i`` and pressures ``p_i``. # # The jacobian, ``K_{ij} = \partial r_i/\partial a_j``, is then split into four parts, # ```math @@ -87,14 +87,14 @@ # K_{ij}^\mathrm{pp} &= \frac{\partial r_i^\mathrm{p}}{\partial p_j} = \int_\Omega \delta N_i^\mathrm{p} \frac{N_j^\mathrm{p}}{\Delta t} + \boldsymbol{\nabla}(\delta N_i^\mathrm{p}) \cdot [k \boldsymbol{\nabla}(N_j^\mathrm{p})] \mathrm{d}\Omega # \end{aligned} # ``` -# -# We could assemble one stiffness matrix and one mass matrix, which would be constant, but for simplicity we only consider a single -# system matrix that depends on the time step, and assemble this for each step. The equations are still linear, so no iterations are required. -# +# +# We could assemble one stiffness matrix and one mass matrix, which would be constant, but for simplicity we only consider a single +# system matrix that depends on the time step, and assemble this for each step. The equations are still linear, so no iterations are required. +# # ## Implementation -# We now solve the problem step by step. The full program with fewer comments is found in +# We now solve the problem step by step. The full program with fewer comments is found in #md # the final [section](@ref porous-media-plain-program) -# +# # Required packages using Ferrite, FerriteMeshParser, Tensors @@ -130,8 +130,8 @@ function Elastic(;E=20.e3, ν=0.3) return Elastic(2G*I4dev + K*I4vol) end; -# Next, we define the element routine for the solid aggregates, where we dispatch on the -# `Elastic` material struct. Note that the unused inputs here are used for the porous matrix below. +# Next, we define the element routine for the solid aggregates, where we dispatch on the +# `Elastic` material struct. Note that the unused inputs here are used for the porous matrix below. function element_routine!(Ke, re, material::Elastic, cv, cell, a, args...) reinit!(cv, cell) n_basefuncs = getnbasefunctions(cv) @@ -143,7 +143,7 @@ function element_routine!(Ke, re, material::Elastic, cv, cell, a, args...) for i in 1:n_basefuncs δ∇N = shape_symmetric_gradient(cv, q_point, i) re[i] += (δ∇N ⊡ σ)*dΩ - for j in 1:n_basefuncs + for j in 1:n_basefuncs ∇N = shape_symmetric_gradient(cv, q_point, j) Ke[i, j] += (δ∇N ⊡ material.C ⊡ ∇N) * dΩ end @@ -152,7 +152,7 @@ function element_routine!(Ke, re, material::Elastic, cv, cell, a, args...) end; # ### PoroElasticity -# To define the poroelastic material, we re-use the elastic part from above for +# To define the poroelastic material, we re-use the elastic part from above for # the skeleton, and add the additional required material parameters. struct PoroElastic{T} elastic::Elastic{T} ## Skeleton stiffness @@ -163,8 +163,8 @@ struct PoroElastic{T} end PoroElastic(;elastic, k, ϕ, α, β) = PoroElastic(elastic, k, ϕ, α, β); -# The element routine requires a few more inputs since we have two fields, as well -# as the dependence on the rates of the displacements and pressure. +# The element routine requires a few more inputs since we have two fields, as well +# as the dependence on the rates of the displacements and pressure. # Again, we dispatch on the material type. function element_routine!(Ke, re, m::PoroElastic, cvs::Tuple, cell, a, a_old, Δt, sdh) ## Setup cellvalues and give easier names @@ -173,14 +173,14 @@ function element_routine!(Ke, re, m::PoroElastic, cvs::Tuple, cell, a, a_old, Δ dr_u = dof_range(sdh, :u) dr_p = dof_range(sdh, :p) - C = m.elastic.C ## Elastic stiffness + C = m.elastic.C ## Elastic stiffness ## Assemble stiffness and force vectors - for q_point in 1:getnquadpoints(cv_u) + for q_point in 1:getnquadpoints(cv_u) dΩ = getdetJdV(cv_u, q_point) p = function_value(cv_p, q_point, a, dr_p) p_old = function_value(cv_p, q_point, a_old, dr_p) - pdot = (p - p_old)/Δt + pdot = (p - p_old)/Δt ∇p = function_gradient(cv_p, q_point, a, dr_p) ϵ = function_symmetric_gradient(cv_u, q_point, a, dr_u) tr_ϵ_old = function_divergence(cv_u, q_point, a_old, dr_u) @@ -213,7 +213,7 @@ function element_routine!(Ke, re, m::PoroElastic, cvs::Tuple, cell, a, a_old, Δ ∇Np = shape_gradient(cv_p, q_point, jₚ) Np = shape_value(cv_p, q_point, jₚ) Ke[Iₚ,Jₚ] += (δNp * m.β*Np/Δt + m.k * (∇δNp ⋅ ∇Np) ) * dΩ - end + end end end end; @@ -226,8 +226,8 @@ struct FEDomain{M,CV,SDH<:SubDofHandler} sdh::SDH end; -# And then we can loop over a vector of such domains, allowing us to -# loop over each domain, to assemble the contributions from each +# And then we can loop over a vector of such domains, allowing us to +# loop over each domain, to assemble the contributions from each # cell in that domain (given by the `SubDofHandler`'s cellset) function doassemble!(K, r, domains::Vector{<:FEDomain}, a, a_old, Δt) assembler = start_assemble(K, r) @@ -240,9 +240,9 @@ end; # we can then loop over all cells in its cellset. Doing this # in a separate function (instead of a nested loop), ensures # that the calls to the `element_routine` are type stable, -# which can be important for good performance. +# which can be important for good performance. function doassemble!(assembler, domain::FEDomain, a, a_old, Δt) - material = domain.material + material = domain.material cv = domain.cellvalues sdh = domain.sdh n = ndofs_per_cell(sdh) @@ -253,7 +253,7 @@ function doassemble!(assembler, domain::FEDomain, a, a_old, Δt) for cell in CellIterator(sdh) ## copy values from a to ae map!(i->a[i], ae, celldofs(cell)) - map!(i->a_old[i], ae_old, celldofs(cell)) + map!(i->a_old[i], ae_old, celldofs(cell)) fill!(Ke, 0) fill!(re, 0) element_routine!(Ke, re, material, cv, cell, ae, ae_old, Δt, sdh) @@ -262,7 +262,7 @@ function doassemble!(assembler, domain::FEDomain, a, a_old, Δt) end; # ### Mesh import -# In this example, we import the mesh from the Abaqus input file, [`porous_media_0p25.inp`](porous_media_0p25.inp) using `FerriteMeshParser`'s +# In this example, we import the mesh from the Abaqus input file, [`porous_media_0p25.inp`](porous_media_0p25.inp) using `FerriteMeshParser`'s # `get_ferrite_grid` function. We then create one cellset for each phase (solid and porous) # for each element type. These 4 sets will later be used in their own `SubDofHandler` function get_grid() @@ -278,12 +278,12 @@ function get_grid() end; # ### Problem setup -# Define the finite element interpolation, integration, and boundary conditions. +# Define the finite element interpolation, integration, and boundary conditions. function setup_problem(;t_rise=0.1, u_max=-0.1) grid = get_grid() - ## Define materials + ## Define materials m_solid = Elastic(;E=20.e3, ν=0.3) m_porous = PoroElastic(;elastic=Elastic(;E=10e3, ν=0.3), β=1/15e3, α=0.9, k=5.0e-3, ϕ=0.8) @@ -292,8 +292,8 @@ function setup_problem(;t_rise=0.1, u_max=-0.1) ipu_tri = Lagrange{RefTriangle,2}()^2 ipp_quad = Lagrange{RefQuadrilateral,1}() ipp_tri = Lagrange{RefTriangle,1}() - - ## Quadrature rules + + ## Quadrature rules qr_quad = QuadratureRule{RefQuadrilateral}(2) qr_tri = QuadratureRule{RefTriangle}(2) @@ -311,11 +311,11 @@ function setup_problem(;t_rise=0.1, u_max=-0.1) ## Solid triangles sdh_solid_tri = SubDofHandler(dh, getcellset(grid,"solid3")) add!(sdh_solid_tri, :u, ipu_tri) - ## Porous quads + ## Porous quads sdh_porous_quad = SubDofHandler(dh, getcellset(grid, "porous4")) add!(sdh_porous_quad, :u, ipu_quad) add!(sdh_porous_quad, :p, ipp_quad) - ## Porous triangles + ## Porous triangles sdh_porous_tri = SubDofHandler(dh, getcellset(grid, "porous3")) add!(sdh_porous_tri, :u, ipu_tri) add!(sdh_porous_tri, :p, ipp_tri) @@ -330,8 +330,8 @@ function setup_problem(;t_rise=0.1, u_max=-0.1) ] ## Boundary conditions - ## Sliding for u, except top which is compressed - ## Sealed for p, except top with prescribed zero pressure + ## Sliding for u, except top which is compressed + ## Sealed for p, except top with prescribed zero pressure addfacetset!(dh.grid, "sides", x -> x[1] < 1e-6 || x[1] ≈ 5.0) addfacetset!(dh.grid, "top", x -> x[2]≈10.0) ch = ConstraintHandler(dh); @@ -345,7 +345,7 @@ function setup_problem(;t_rise=0.1, u_max=-0.1) end; # ### Solving -# Given the `DofHandler`, `ConstraintHandler`, and `CellValues`, +# Given the `DofHandler`, `ConstraintHandler`, and `CellValues`, # we can solve the problem by stepping through the time history function solve(dh, ch, domains; Δt=0.025, t_total=1.0) K = create_sparsity_pattern(dh); @@ -356,7 +356,7 @@ function solve(dh, ch, domains; Δt=0.025, t_total=1.0) for t in 0:Δt:t_total if t>0 update!(ch, t) - apply!(a, ch) + apply!(a, ch) doassemble!(K, r, domains, a, a_old, Δt) apply_zero!(K, r, ch) Δa = -K\r @@ -382,4 +382,4 @@ solve(dh, ch, domains); #md # #md # ```julia #md # @__CODE__ -#md # ``` \ No newline at end of file +#md # ``` diff --git a/docs/src/literate-tutorials/transient_heat_equation.jl b/docs/src/literate-tutorials/transient_heat_equation.jl index 2a865bed08..df64dc7913 100644 --- a/docs/src/literate-tutorials/transient_heat_equation.jl +++ b/docs/src/literate-tutorials/transient_heat_equation.jl @@ -182,7 +182,7 @@ A = (Δt .* K) + M; # by `get_rhs_data`. The function returns a `RHSData` struct, which contains all needed information to apply # the boundary conditions solely on the right-hand-side vector of the problem. rhsdata = get_rhs_data(ch, A); -# We set the values at initial time step, denoted by uₙ, to a bubble-shape described by +# We set the values at initial time step, denoted by uₙ, to a bubble-shape described by # $(x_1^2-1)(x_2^2-1)$, such that it is zero at the boundaries and the maximum temperature in the center. uₙ = zeros(length(f)); apply_analytical!(uₙ, dh, :u, x -> (x[1]^2 - 1) * (x[2]^2 - 1) * max_temp); diff --git a/docs/src/reference/fevalues.md b/docs/src/reference/fevalues.md index 26f050898d..2ef3f8204b 100644 --- a/docs/src/reference/fevalues.md +++ b/docs/src/reference/fevalues.md @@ -6,8 +6,8 @@ DocTestSetup = :(using Ferrite) # FEValues ## Main types -[`CellValues`](@ref) and [`FacetValues`](@ref) are the most common -subtypes of `Ferrite.AbstractValues`. For more details about how +[`CellValues`](@ref) and [`FacetValues`](@ref) are the most common +subtypes of `Ferrite.AbstractValues`. For more details about how these work, please see the related [topic guide](@ref fevalues_topicguide). ```@docs diff --git a/docs/src/topics/FEValues.md b/docs/src/topics/FEValues.md index d40d40bef3..df54430cbc 100644 --- a/docs/src/topics/FEValues.md +++ b/docs/src/topics/FEValues.md @@ -2,15 +2,15 @@ A key type of object in `Ferrite.jl` is the so-called `FEValues`, where the most common ones are `CellValues` and `FacetValues`. These objects are used inside the element routines and are used to query the integration weights, shape function values and gradients, and much more; see [`CellValues`](@ref) and [`FacetValues`](@ref). For these values to be correct, it is necessary to reinitialize these for the current cell by using the [`reinit!`](@ref) function. This function maps the values from the reference cell to the actual cell, a process described in detail below, see [Mapping of finite elements](@ref mapping_theory). After that, we show an implementation of a [`SimpleCellValues`](@ref SimpleCellValues) type to illustrate how `CellValues` work for the most standard case, excluding the generalizations and optimization that complicates the actual code. ## [Mapping of finite elements](@id mapping_theory) -The shape functions and gradients stored in an `FEValues` object, are reinitialized for each cell by calling the `reinit!` function. -The main part of this calculation, considers how to map the values and derivatives of the shape functions, +The shape functions and gradients stored in an `FEValues` object, are reinitialized for each cell by calling the `reinit!` function. +The main part of this calculation, considers how to map the values and derivatives of the shape functions, defined on the reference cell, to the actual cell. -The geometric mapping of a finite element from the reference coordinates to the real coordinates is shown in the following illustration. +The geometric mapping of a finite element from the reference coordinates to the real coordinates is shown in the following illustration. ![mapping_figure](https://raw.githubusercontent.com/Ferrite-FEM/Ferrite.jl/gh-pages/assets/fe_mapping.svg) -This mapping is given by the geometric shape functions, $\hat{N}_i^g(\boldsymbol{\xi})$, such that +This mapping is given by the geometric shape functions, $\hat{N}_i^g(\boldsymbol{\xi})$, such that ```math \begin{align*} \boldsymbol{x}(\boldsymbol{\xi}) =& \sum_{\alpha=1}^N \hat{\boldsymbol{x}}_\alpha \hat{N}_\alpha^g(\boldsymbol{\xi}) \\ @@ -24,8 +24,8 @@ where the defined $\boldsymbol{J}$ is the jacobian of the mapping, and in some c We require that the mapping from reference coordinates to real coordinates is [diffeomorphic](https://en.wikipedia.org/wiki/Diffeomorphism), meaning that we can express $\boldsymbol{x} = \boldsymbol{x}(\boldsymbol{\xi}(\boldsymbol{x}))$, such that ```math \begin{align*} - \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{x}} = \boldsymbol{I} &= \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{\xi}} \cdot \frac{\mathrm{d}\boldsymbol{\xi}}{\mathrm{d}\boldsymbol{x}} - \quad\Rightarrow\quad + \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{x}} = \boldsymbol{I} &= \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{\xi}} \cdot \frac{\mathrm{d}\boldsymbol{\xi}}{\mathrm{d}\boldsymbol{x}} + \quad\Rightarrow\quad \frac{\mathrm{d}\boldsymbol{\xi}}{\mathrm{d}\boldsymbol{x}} = \left[\frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{\xi}}\right]^{-1} = \boldsymbol{J}^{-1} \end{align*} ``` @@ -34,7 +34,7 @@ Depending on the function interpolation, we may want different types of mappings ### Identity mapping `Ferrite.IdentityMapping` -For scalar fields, we always use scalar base functions. For tensorial fields (non-scalar, e.g. vector-fields), the base functions can be constructed from scalar base functions, by using e.g. `VectorizedInterpolation`. From the perspective of the mapping, however, each component is mapped as an individual scalar base function. And for scalar base functions, we only require that the value of the base function is invariant to the element shape (real coordinate), and only depends on the reference coordinate, i.e. +For scalar fields, we always use scalar base functions. For tensorial fields (non-scalar, e.g. vector-fields), the base functions can be constructed from scalar base functions, by using e.g. `VectorizedInterpolation`. From the perspective of the mapping, however, each component is mapped as an individual scalar base function. And for scalar base functions, we only require that the value of the base function is invariant to the element shape (real coordinate), and only depends on the reference coordinate, i.e. ```math \begin{align*} N(\boldsymbol{x}) &= \hat{N}(\boldsymbol{\xi}(\boldsymbol{x}))\nonumber \\ @@ -45,7 +45,7 @@ For scalar fields, we always use scalar base functions. For tensorial fields (no ### Covariant Piola mapping, H(curl) `Ferrite.CovariantPiolaMapping` -The covariant Piola mapping of a vectorial base function preserves the tangential components. For the value, the mapping is defined as +The covariant Piola mapping of a vectorial base function preserves the tangential components. For the value, the mapping is defined as ```math \begin{align*} \boldsymbol{N}(\boldsymbol{x}) = \boldsymbol{J}^{-\mathrm{T}} \cdot \hat{\boldsymbol{N}}(\boldsymbol{\xi}(\boldsymbol{x})) @@ -66,7 +66,7 @@ which yields the gradient, \end{align*} ``` - Except for a few elements, $\boldsymbol{J}$ varies as a function of $\boldsymbol{x}$. The derivative can be calculated as + Except for a few elements, $\boldsymbol{J}$ varies as a function of $\boldsymbol{x}$. The derivative can be calculated as ```math \begin{align*} \frac{\mathrm{d} J^{-\mathrm{T}}_{ik}}{\mathrm{d} x_j} &= \frac{\mathrm{d} J^{-\mathrm{T}}_{ik}}{\mathrm{d} J_{mn}} \frac{\mathrm{d} J_{mn}}{\mathrm{d} x_j} = - J_{km}^{-1} J_{in}^{-T} \frac{\mathrm{d} J_{mn}}{\mathrm{d} x_j} \nonumber \\ @@ -77,7 +77,7 @@ which yields the gradient, ### Contravariant Piola mapping, H(div) `Ferrite.ContravariantPiolaMapping` -The covariant Piola mapping of a vectorial base function preserves the normal components. For the value, the mapping is defined as +The covariant Piola mapping of a vectorial base function preserves the normal components. For the value, the mapping is defined as ```math \begin{align*} \boldsymbol{N}(\boldsymbol{x}) = \frac{\boldsymbol{J}}{\det(\boldsymbol{J})} \cdot \hat{\boldsymbol{N}}(\boldsymbol{\xi}(\boldsymbol{x})) @@ -97,19 +97,19 @@ This gives the gradient ```math \begin{align*} \frac{\mathrm{d} N_i}{\mathrm{d} x_j} &= \frac{\mathrm{d}}{\mathrm{d} x_j} \left[\frac{J_{ik}}{\det(\boldsymbol{J})} \hat{N}_k\right] =\nonumber\\ - &= \frac{\mathrm{d} J_{ik}}{\mathrm{d} x_j} \frac{\hat{N}_k}{\det(\boldsymbol{J})} + &= \frac{\mathrm{d} J_{ik}}{\mathrm{d} x_j} \frac{\hat{N}_k}{\det(\boldsymbol{J})} - \frac{\mathrm{d} \det(\boldsymbol{J})}{\mathrm{d} x_j} \frac{J_{ik} \hat{N}_k}{\det(\boldsymbol{J})^2} + \frac{J_{ik}}{\det(\boldsymbol{J})} \frac{\mathrm{d} \hat{N}_k}{\mathrm{d} \xi_l} J_{lj}^{-1} \\ - &= \mathcal{H}_{ikl} J^{-1}_{lj} \frac{\hat{N}_k}{\det(\boldsymbol{J})} + &= \mathcal{H}_{ikl} J^{-1}_{lj} \frac{\hat{N}_k}{\det(\boldsymbol{J})} - J^{-T}_{mn} \mathcal{H}_{mnl} J^{-1}_{lj} \frac{J_{ik} \hat{N}_k}{\det(\boldsymbol{J})} + \frac{J_{ik}}{\det(\boldsymbol{J})} \frac{\mathrm{d} \hat{N}_k}{\mathrm{d} \xi_l} J_{lj}^{-1} \end{align*} ``` ## [Walkthrough: Creating `SimpleCellValues`](@id SimpleCellValues) -In the following, we walk through how to create a `SimpleCellValues` type which -works similar to `Ferrite.jl`'s `CellValues`, but is not performance optimized and not as general. The main purpose is to explain how the `CellValues` works for the standard case of `IdentityMapping` described above. -Please note that several internal functions are used, and these may change without a major version increment. Please see the [Developer documentation](@ref) for their documentation. +In the following, we walk through how to create a `SimpleCellValues` type which +works similar to `Ferrite.jl`'s `CellValues`, but is not performance optimized and not as general. The main purpose is to explain how the `CellValues` works for the standard case of `IdentityMapping` described above. +Please note that several internal functions are used, and these may change without a major version increment. Please see the [Developer documentation](@ref) for their documentation. ```@eval # Include the example here, but modify the Literate output to suit being embedded diff --git a/docs/src/topics/SimpleCellValues_literate.jl b/docs/src/topics/SimpleCellValues_literate.jl index 06ba44b779..e9fd664e55 100644 --- a/docs/src/topics/SimpleCellValues_literate.jl +++ b/docs/src/topics/SimpleCellValues_literate.jl @@ -2,8 +2,8 @@ using Ferrite, Test # Then, we define a simple version of the cell values object, which only supports -# * Scalar interpolations -# * Identity mapping from reference to physical cell. +# * Scalar interpolations +# * Identity mapping from reference to physical cell. # * The cell shape has the same dimension as the physical space (excludes so-called embedded cells). struct SimpleCellValues{T, dim} <: Ferrite.AbstractCellValues @@ -11,7 +11,7 @@ struct SimpleCellValues{T, dim} <: Ferrite.AbstractCellValues ## shape function number and q_point the integration point dNdξ::Matrix{Vec{dim,T}} # Precalculated shape gradients in the reference domain, dNdξ[i, q_point] dNdx::Matrix{Vec{dim,T}} # Cache for shape gradients in the physical domain, dNdx[i, q_point] - M::Matrix{T} # Precalculated geometric shape values, M[j, q_point] where j is the + M::Matrix{T} # Precalculated geometric shape values, M[j, q_point] where j is the ## geometric shape function number dMdξ::Matrix{Vec{dim,T}} # Precalculated geometric shape gradients, dMdξ[j, q_point] weights::Vector{T} # Given quadrature weights in the reference domain, weights[q_point] @@ -27,7 +27,7 @@ function SimpleCellValues(qr::QuadratureRule, ip_fun::Interpolation, ip_geo::Int weights = Ferrite.getweights(qr) n_qpoints = length(weights) T = eltype(weights) - + ## Function interpolation n_func_basefuncs = getnbasefunctions(ip_fun) N = zeros(T, n_func_basefuncs, n_qpoints) @@ -53,15 +53,15 @@ function SimpleCellValues(qr::QuadratureRule, ip_fun::Interpolation, ip_geo::Int SimpleCellValues(N, dNdξ, dNdx, M, dMdξ, weights, detJdV) end; -# To make our `SimpleCellValues` work in standard Ferrite code, +# To make our `SimpleCellValues` work in standard Ferrite code, # we need to dispatch some access functions: Ferrite.getnbasefunctions(cv::SimpleCellValues) = size(cv.N, 1) Ferrite.getnquadpoints(cv::SimpleCellValues) = size(cv.N, 2) Ferrite.shape_value(cv::SimpleCellValues, q_point::Int, i::Int) = cv.N[i, q_point] Ferrite.shape_gradient(cv::SimpleCellValues, q_point::Int, i::Int) = cv.dNdx[i, q_point]; -# The last step is then to dispatch `reinit!` for our `SimpleCellValues` to calculate -# the cached values `dNdx` and `detJdV` for the current cell according to the +# The last step is then to dispatch `reinit!` for our `SimpleCellValues` to calculate +# the cached values `dNdx` and `detJdV` for the current cell according to the # theory for `IdentityMapping` above. function Ferrite.reinit!(cv::SimpleCellValues, x::Vector{Vec{dim,T}}) where {dim,T} for (q_point, w) in pairs(cv.weights) # Loop over each quadrature point @@ -72,7 +72,7 @@ function Ferrite.reinit!(cv::SimpleCellValues, x::Vector{Vec{dim,T}}) where {dim end ## Calculate the correct integration weight for the current q_point cv.detJdV[q_point] = det(J)*w - ## map the shape gradients to the current geometry + ## map the shape gradients to the current geometry Jinv = inv(J) for i in 1:getnbasefunctions(cv) cv.dNdx[i, q_point] = cv.dNdξ[i, q_point] ⋅ Jinv @@ -98,4 +98,3 @@ ue = rand(getnbasefunctions(simple_cv)) q_point = 2 @test function_value(cv, q_point, ue) ≈ function_value(simple_cv, q_point, ue) @test function_gradient(cv, q_point, ue) ≈ function_gradient(simple_cv, q_point, ue) - diff --git a/docs/src/topics/boundary_conditions.md b/docs/src/topics/boundary_conditions.md index 84fa368d99..b87b61cf57 100644 --- a/docs/src/topics/boundary_conditions.md +++ b/docs/src/topics/boundary_conditions.md @@ -297,9 +297,9 @@ pdbc = PeriodicDirichlet( ## Initial Conditions -When solving time-dependent problems, initial conditions, different from zero, may be required. -For finite element formulations of ODE-type, -i.e. ``\boldsymbol{u}'(t) = \boldsymbol{f}(\boldsymbol{u}(t),t)``, +When solving time-dependent problems, initial conditions, different from zero, may be required. +For finite element formulations of ODE-type, +i.e. ``\boldsymbol{u}'(t) = \boldsymbol{f}(\boldsymbol{u}(t),t)``, where ``\boldsymbol{u}(t)`` are the degrees of freedom, initial conditions can be specified by the [`apply_analytical!`](@ref) function. For example, specify the initial pressure as a function of the y-coordinate @@ -319,4 +319,3 @@ See also [Transient heat equation](@ref tutorial-transient-heat-equation) for on equations (DAEs) need extra care during initialization. We refer to the paper ["Consistent Initial Condition Calculation for Differential-Algebraic Systems" by Brown et al.](https://dx.doi.org/10.1137/S1064827595289996) for more details on this matter. - diff --git a/docs/src/topics/constraints.md b/docs/src/topics/constraints.md index 8d232ae90a..ac2d39a19f 100644 --- a/docs/src/topics/constraints.md +++ b/docs/src/topics/constraints.md @@ -4,7 +4,7 @@ DocTestSetup = :(using Ferrite) # Constraints -PDEs can in general be subjected to a number of constraints, +PDEs can in general be subjected to a number of constraints, ```math g_I(\underline{a}) = 0, \quad I = 1 \text{ to } n_c @@ -12,7 +12,7 @@ g_I(\underline{a}) = 0, \quad I = 1 \text{ to } n_c where $g$ are (non-linear) constraint equations, $\underline{a}$ is a vector of the degrees of freedom, and $n_c$ is the number of constraints. There are many ways to -enforce these constraints, e.g. penalty methods and Lagrange multiplier methods. +enforce these constraints, e.g. penalty methods and Lagrange multiplier methods. ## Affine constraints @@ -35,7 +35,7 @@ add!(ch, lc1) add!(ch, lc2) ``` -Affine constraints will affect the sparsity pattern of the stiffness matrix, and as such, it is important to also include +Affine constraints will affect the sparsity pattern of the stiffness matrix, and as such, it is important to also include the `ConstraintHandler` as an argument when creating the sparsity pattern: ```julia @@ -43,9 +43,9 @@ K = create_sparsity_pattern(dh, ch) ``` ### Solving linear problems -To solve the system ``\underline{\underline{K}}\underline{a}=\underline{f}``, account for affine constraints the same way as for +To solve the system ``\underline{\underline{K}}\underline{a}=\underline{f}``, account for affine constraints the same way as for `Dirichlet` boundary conditions; first call `apply!(K, f, ch)`. This will condense `K` and `f` inplace (i.e -no new matrix will be created). Note however that we must also call `apply!` on the solution vector after +no new matrix will be created). Note however that we must also call `apply!` on the solution vector after solving the system to enforce the affine constraints: ```julia @@ -59,17 +59,17 @@ apply!(a, ch) # enforces affine constraints ``` ### Solving nonlinear problems -It is important to check the residual **after** applying boundary conditions when -solving nonlinear problems with affine constraints. -`apply_zero!(K, r, ch)` modifies the residual entries for dofs that are involved -in constraints to account for constraint forces. +It is important to check the residual **after** applying boundary conditions when +solving nonlinear problems with affine constraints. +`apply_zero!(K, r, ch)` modifies the residual entries for dofs that are involved +in constraints to account for constraint forces. The following pseudo-code shows a typical pattern for solving a non-linear problem with Newton's method: ```julia a = initial_guess(...) # Make any initial guess for a here, e.g. `a=zeros(ndofs(dh))` apply!(a, ch) # Make the guess fulfill all constraints in `ch` for iter in 1:maxiter doassemble!(K, r, ...) # Assemble the residual, r, and stiffness, K=∂r/∂a. - apply_zero!(K, r, ch) # Modify `K` and `r` to account for the constraints. + apply_zero!(K, r, ch) # Modify `K` and `r` to account for the constraints. check_convergence(r, ...) && break # Only check convergence after `apply_zero!(K, r, ch)` Δa = K \ r # Calculate the (negative) update apply_zero!(Δa, ch) # Change the constrained values in `Δa` such that `a-Δa` diff --git a/docs/src/topics/export.md b/docs/src/topics/export.md index 702561f522..e279d172c5 100644 --- a/docs/src/topics/export.md +++ b/docs/src/topics/export.md @@ -10,7 +10,7 @@ u = rand(ndofs(dh)); σ = rand(getncells(grid)) When the problem is solved, and the solution vector `u` is known we typically want to visualize it. The simplest way to do this is to write the solution to a VTK-file, which can be viewed in e.g. [`Paraview`](https://www.paraview.org/). -To write VTK-files, Ferrite comes with an export interface with a +To write VTK-files, Ferrite comes with an export interface with a [`WriteVTK.jl`](https://github.com/jipolanco/WriteVTK.jl) backend to simplify the exporting. @@ -20,7 +20,7 @@ VTKFile("my_solution", grid) do vtk write_solution(vtk, dh, u) end; ``` -where `write_solution` is just one example of the following functions that can be used +where `write_solution` is just one example of the following functions that can be used * [`write_solution`](@ref) * [`write_cell_data`](@ref) @@ -39,14 +39,14 @@ write_solution(vtk, dh, u) close(vtk); ``` -The data written by `write_solution`, `write_cell_data`, `write_node_data`, and `write_projection` may be either scalar (`Vector{<:Number}`) or tensor (`Vector{<:AbstractTensor}`) data. +The data written by `write_solution`, `write_cell_data`, `write_node_data`, and `write_projection` may be either scalar (`Vector{<:Number}`) or tensor (`Vector{<:AbstractTensor}`) data. -For simulations with multiple time steps, typically one `VTK` (`.vtu`) file is written +For simulations with multiple time steps, typically one `VTK` (`.vtu`) file is written for each time step. In order to connect the actual time with each of these files, a `VTKFileCollection` can be used, which will write one paraview datafile (`.pvd`) -file and one `VTKFile` (`.vtu`) for each time step. +file and one `VTKFile` (`.vtu`) for each time step. -```@example export +```@example export pvd = VTKFileCollection("my_results", grid) for t in range(0, 1, 5) # Do calculations to update u @@ -56,4 +56,4 @@ for t in range(0, 1, 5) end close(pvd); ``` -See [Transient heat equation](@ref tutorial-transient-heat-equation) for an example \ No newline at end of file +See [Transient heat equation](@ref tutorial-transient-heat-equation) for an example diff --git a/docs/src/topics/grid.md b/docs/src/topics/grid.md index f49f117835..3b69389979 100644 --- a/docs/src/topics/grid.md +++ b/docs/src/topics/grid.md @@ -6,8 +6,8 @@ DocTestSetup = :(using Ferrite) ## Mesh Reading -A Ferrite `Grid` can be generated with the [`generate_grid`](@ref) function. -More advanced meshes can be imported with the +A Ferrite `Grid` can be generated with the [`generate_grid`](@ref) function. +More advanced meshes can be imported with the [`FerriteMeshParser.jl`](https://github.com/Ferrite-FEM/FerriteMeshParser.jl) (from Abaqus input files), or even created and translated with the [`Gmsh.jl`](https://github.com/JuliaFEM/Gmsh.jl) and [`FerriteGmsh.jl`](https://github.com/Ferrite-FEM/FerriteGmsh.jl) package, respectively. @@ -43,7 +43,7 @@ FerriteMeshParser.get_ferrite_grid ``` If you are missing the translation of an Abaqus element that is equivalent to a `Ferrite.Cell`, -consider to open an [issue](https://github.com/Ferrite-FEM/FerriteMeshParser.jl/issues/new) or a pull request. +consider to open an [issue](https://github.com/Ferrite-FEM/FerriteMeshParser.jl/issues/new) or a pull request. ## `Grid` Datastructure @@ -71,12 +71,12 @@ julia> cells = [ ``` where each Quadrilateral, which is a subtype of `AbstractCell` saves in the field `nodes` the tuple of node IDs. -Additionally, the data structure `Grid` can hold node-, face- and cellsets. -All of these three sets are defined by a dictionary that maps a string key to a `Set`. -For the special case of node- and cellsets the dictionary's value is of type `Set{Int}`, i.e. a keyword is mapped to a node or cell ID, respectively. +Additionally, the data structure `Grid` can hold node-, face- and cellsets. +All of these three sets are defined by a dictionary that maps a string key to a `Set`. +For the special case of node- and cellsets the dictionary's value is of type `Set{Int}`, i.e. a keyword is mapped to a node or cell ID, respectively. Facesets are a more elaborate construction. They map a `String` key to a `Set{FaceIndex}`, where each `FaceIndex` consists of `(global_cell_id, local_face_id)`. -In order to understand the `local_face_id` properly, one has to consider the reference space of the element, which typically is spanned by a product of the interval ``[-1, 1]`` and in this particular example ``[-1, 1] \times [-1, 1]``. +In order to understand the `local_face_id` properly, one has to consider the reference space of the element, which typically is spanned by a product of the interval ``[-1, 1]`` and in this particular example ``[-1, 1] \times [-1, 1]``. In this space a local numbering of nodes and faces exists, i.e. @@ -113,13 +113,13 @@ julia> function compute_faceset(cells, global_faces, ip::Interpolation{dim}) whe d[ntuple(i-> cell.nodes[face[i]], nodes_per_face)] = FaceIndex(c, f) end end - + faces = Vector{FaceIndex}() for face in global_faces # lookup the element, local face combination for this face push!(faces, d[face]) end - + return faces end @@ -173,8 +173,8 @@ Ferrite.get_coordinate_type(::SmallGrid{dim}) where dim = Vec{dim,Float64} Ferrite.nnodes_per_cell(grid::SmallGrid, i::Int=1) = Ferrite.nnodes(grid.cells_test[i]) ``` -These definitions make many of `Ferrite`s functions work out of the box, e.g. you can now call -`getcoordinates(grid, cellid)` on the `SmallGrid`. +These definitions make many of `Ferrite`s functions work out of the box, e.g. you can now call +`getcoordinates(grid, cellid)` on the `SmallGrid`. Now, you would be able to assemble the heat equation example over the new custom `SmallGrid` type. Note that this particular subtype isn't able to handle boundary entity sets and so, you can't describe boundaries with it. diff --git a/docs/src/tutorials/index.md b/docs/src/tutorials/index.md index 80f3015721..cd135b7c79 100644 --- a/docs/src/tutorials/index.md +++ b/docs/src/tutorials/index.md @@ -126,9 +126,9 @@ Gmsh. #### [Tutorial 9: Porous media (SubDofHandler)](porous_media.md) -This tutorial introduces how to solve a complex linear problem, where there are different +This tutorial introduces how to solve a complex linear problem, where there are different fields on different subdomains, and different cell types in the grid. This requires using -the `SubDofHandler` interface. +the `SubDofHandler` interface. **Keywords**: Mixed grids, multiple fields, porous media, `SubDofHandler` diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 0cf3541974..72cc2f7bca 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -58,8 +58,8 @@ const DofCoefficients{T} = Vector{Pair{Int,T}} """ AffineConstraint(constrained_dof::Int, entries::Vector{Pair{Int,T}}, b::T) where T -Define an affine/linear constraint to constrain one degree of freedom, `u[i]`, -such that `u[i] = ∑(u[j] * a[j]) + b`, +Define an affine/linear constraint to constrain one degree of freedom, `u[i]`, +such that `u[i] = ∑(u[j] * a[j]) + b`, where `i=constrained_dof` and each element in `entries` are `j => a[j]` """ struct AffineConstraint{T} @@ -177,7 +177,7 @@ isclosed(ch::ConstraintHandler) = ch.closed[] free_dofs(ch::ConstraintHandler) = ch.free_dofs prescribed_dofs(ch::ConstraintHandler) = ch.prescribed_dofs -# Equivalent to `copy!(out, setdiff(1:n_entries, diff))`, but requires that +# Equivalent to `copy!(out, setdiff(1:n_entries, diff))`, but requires that # `issorted(diff)` and that `1 ≤ diff[1] ≤ diff[end] ≤ n_entries` function _sorted_setdiff!(out::Vector{Int}, n_entries::Int, diff::Vector{Int}) n_diff = length(diff) @@ -528,12 +528,12 @@ apply! apply_zero!(K::SparseMatrixCSC, rhs::AbstractVector, ch::ConstraintHandler) Adjust the matrix `K` and the right hand side `rhs` to account for prescribed Dirichlet -boundary conditions and affine constraints such that `du = K \\ rhs` gives the expected +boundary conditions and affine constraints such that `du = K \\ rhs` gives the expected result (e.g. `du` zero for all prescribed degrees of freedom). apply_zero!(v::AbstractVector, ch::ConstraintHandler) -Zero-out values in `v` corresponding to prescribed degrees of freedom and update values +Zero-out values in `v` corresponding to prescribed degrees of freedom and update values prescribed by affine constraints, such that if `a` fulfills the constraints, `a ± v` also will. @@ -553,9 +553,9 @@ apply_zero!(ΔΔu, ch) # Make sure values are exactly zero ``` !!! note - The last call to `apply_zero!` is only strictly necessary for affine constraints. - However, even if the Dirichlet boundary conditions should be fulfilled after - `apply!(K, g, ch)`, solvers of linear systems are not exact. + The last call to `apply_zero!` is only strictly necessary for affine constraints. + However, even if the Dirichlet boundary conditions should be fulfilled after + `apply!(K, g, ch)`, solvers of linear systems are not exact. `apply!(ΔΔu, ch)` can be used to make sure the values for the prescribed degrees of freedom are fulfilled exactly. """ @@ -1263,7 +1263,7 @@ system. For other types of periodicities the `transform` function can be used. T `transform` function is applied on the coordinates of the image facet, and is expected to transform the coordinates to the matching locations in the mirror set. -The keyword `tol` specifies the tolerance (i.e. distance and deviation in facet-normals) +The keyword `tol` specifies the tolerance (i.e. distance and deviation in facet-normals) between a image-facet and mirror-facet, for them to be considered matched. See also: [`collect_periodic_facets!`](@ref), [`PeriodicDirichlet`](@ref). diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index 9c6eabf65c..a115147d8c 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -6,7 +6,7 @@ abstract type AbstractDofHandler end Access some grid representation for the dof handler. !!! note - This API function is currently not well-defined. It acts as the interface between + This API function is currently not well-defined. It acts as the interface between distributed assembly and assembly on a single process, because most parts of the functionality can be handled by only acting on the locally owned cell set. """ @@ -28,17 +28,17 @@ end """ SubDofHandler(dh::AbstractDofHandler, cellset::AbstractVecOrSet{Int}) -Create an `sdh::SubDofHandler` from the parent `dh`, pertaining to the -cells in `cellset`. This allows you to add fields to parts of the domain, or using -different interpolations or cell types (e.g. `Triangles` and `Quadrilaterals`). All +Create an `sdh::SubDofHandler` from the parent `dh`, pertaining to the +cells in `cellset`. This allows you to add fields to parts of the domain, or using +different interpolations or cell types (e.g. `Triangles` and `Quadrilaterals`). All fields and cell types must be the same in one `SubDofHandler`. After construction any number of discrete fields can be added to the SubDofHandler using [`add!`](@ref). Construction is finalized by calling [`close!`](@ref) on the parent `dh`. # Examples -We assume we have a `grid` containing "Triangle" and "Quadrilateral" cells, -including the cellsets "triangles" and "quadilaterals" for to these cells. +We assume we have a `grid` containing "Triangle" and "Quadrilateral" cells, +including the cellsets "triangles" and "quadilaterals" for to these cells. ```julia dh = DofHandler(grid) @@ -50,7 +50,7 @@ sdh_quad = SubDofHandler(dh, getcellset(grid, "quadilaterals")) ip_quad = Lagrange{RefQuadrilateral, 2}()^2 # vector interpolation for a field u add!(sdh_quad, :u, ip_quad) -close!(dh) # Finalize by closing the parent +close!(dh) # Finalize by closing the parent ``` """ function SubDofHandler(dh::DH, cellset::AbstractVecOrSet{Int}) where {DH <: AbstractDofHandler} @@ -270,7 +270,7 @@ function add!(sdh::SubDofHandler, name::Symbol, ip::Interpolation) # TODO: warn if interpolation type is not the same? end end - + # Check that interpolation is compatible with cells it it added to refshape_sdh = getrefshape(getcells(sdh.dh.grid, first(sdh.cellset))) if refshape_sdh !== getrefshape(ip) @@ -333,7 +333,7 @@ For the `DofHandler` each `SubDofHandler` is visited in the order they were adde For each field in the `SubDofHandler` create dofs for the cell. This means that dofs on a particular cell will be numbered in groups for each field, so first the dofs for field 1 are distributed, then field 2, etc. -For each cell dofs are first distributed on its vertices, then on the interior of edges (if applicable), then on the +For each cell dofs are first distributed on its vertices, then on the interior of edges (if applicable), then on the interior of faces (if applicable), and finally on the cell interior. The entity ordering follows the geometrical ordering found in [`vertices`](@ref), [`faces`](@ref) and [`edges`](@ref). """ @@ -360,7 +360,7 @@ function __close!(dh::DofHandler{dim}) where {dim} edgedicts = [Dict{NTuple{2, Int}, Int}() for _ in 1:numfields] # `facedict` keeps track of the visited faces. We only need to store the first dof we - # add to the face since currently more dofs per face isn't supported. + # add to the face since currently more dofs per face isn't supported. # A face is uniquely determined by 3 vertex nodes, see sortface facedicts = [Dict{NTuple{3, Int}, Int}() for _ in 1:numfields] @@ -487,14 +487,14 @@ function _distribute_dofs_for_cell!(dh::DofHandler{sdim}, cell::AbstractCell, ip ip_info.nvertexdofs, nextdof, ip_info.n_copies, ) - # Distribute dofs for edges + # Distribute dofs for edges nextdof = add_edge_dofs( dh.cell_dofs, cell, edgedict, ip_info.nedgedofs, nextdof, ip_info.adjust_during_distribution, ip_info.n_copies, ) - # Distribute dofs for faces. + # Distribute dofs for faces. nextdof = add_face_dofs( dh.cell_dofs, cell, facedict, ip_info.nfacedofs, nextdof, @@ -674,7 +674,7 @@ Here the unique representation is the sorted node index tuple. Note that in 3D we only need indices to uniquely identify a face, so the unique representation is always a tuple length 3. """ -function sortface end +function sortface end """ sortface_fast(face::Tuple{Int}) @@ -694,9 +694,9 @@ function sortface_fast end For more details we refer to [1] as we follow the methodology described therein. -[1] Scroggs, M. W., Dokken, J. S., Richardson, C. N., & Wells, G. N. (2022). - Construction of arbitrary order finite element degree-of-freedom maps on - polygonal and polyhedral cell meshes. ACM Transactions on Mathematical +[1] Scroggs, M. W., Dokken, J. S., Richardson, C. N., & Wells, G. N. (2022). + Construction of arbitrary order finite element degree-of-freedom maps on + polygonal and polyhedral cell meshes. ACM Transactions on Mathematical Software (TOMS), 48(2), 1-23. !!!TODO citation via software. @@ -762,16 +762,16 @@ end """ sortfacet_fast(facet::NTuple{N, Int}) -Returns the unique representation of the `facet` by sorting its node indices. +Returns the unique representation of the `facet` by sorting its node indices. Dispatches on `sortedges_fast` or `sortfaces_fast` depending on `N` """ -function sortfacet_fast end +function sortfacet_fast end # Vertex sortfacet_fast(facet::Tuple{Int}) = facet -# Edge +# Edge sortfacet_fast(facet::NTuple{2, Int}) = sortedge_fast(facet) -# Face +# Face sortfacet_fast(facet::NTuple{3, Int}) = sortface_fast(facet) sortfacet_fast(facet::NTuple{4, Int}) = sortface_fast(facet) @@ -816,7 +816,7 @@ end """ _find_field(sdh::SubDofHandler, field_name::Symbol)::Int -Return the index of the field with name `field_name` in the `SubDofHandler` `sdh`. Return +Return the index of the field with name `field_name` in the `SubDofHandler` `sdh`. Return `nothing` if the field is not found. See also: [`find_field(dh::DofHandler, field_name::Symbol)`](@ref), [`find_field(sdh::SubDofHandler, field_name::Symbol)`](@ref). diff --git a/src/Dofs/apply_analytical.jl b/src/Dofs/apply_analytical.jl index 847fd58f8c..ed355401f6 100644 --- a/src/Dofs/apply_analytical.jl +++ b/src/Dofs/apply_analytical.jl @@ -6,7 +6,7 @@ end """ apply_analytical!( - a::AbstractVector, dh::AbstractDofHandler, fieldname::Symbol, + a::AbstractVector, dh::AbstractDofHandler, fieldname::Symbol, f::Function, cellset=1:getncells(get_grid(dh))) Apply a solution `f(x)` by modifying the values in the degree of freedom vector `a` @@ -18,7 +18,7 @@ and for vector fields with dimension `dim`, `f(x)::Vec{dim}`. This function can be used to apply initial conditions for time dependent problems. !!! note - + This function only works for standard nodal finite element interpolations when the function value at the (algebraic) node is equal to the corresponding degree of freedom value. diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl index 5df796035f..bdca14c191 100644 --- a/src/Export/VTK.jl +++ b/src/Export/VTK.jl @@ -3,8 +3,8 @@ VTKFile(filename::AbstractString, grid::AbstractGrid; kwargs...) VTKFile(filename::AbstractString, dh::DofHandler; kwargs...) -Create a `VTKFile` that contains an unstructured VTK grid. -The keyword arguments are forwarded to `WriteVTK.vtk_grid`, see +Create a `VTKFile` that contains an unstructured VTK grid. +The keyword arguments are forwarded to `WriteVTK.vtk_grid`, see [Data Formatting Options](https://juliavtk.github.io/WriteVTK.jl/stable/grids/syntax/#Data-formatting-options) This file handler can be used to to write data with @@ -17,7 +17,7 @@ This file handler can be used to to write data with * [`Ferrite.write_nodeset`](@ref) * [`Ferrite.write_constraints`](@ref) -It is necessary to call `close(::VTKFile)` to save the data after writing +It is necessary to call `close(::VTKFile)` to save the data after writing to the file handler. Using the supported `do`-block does this automatically: ```julia VTKFile(filename, grid) do vtk @@ -58,8 +58,8 @@ end VTKFileCollection(name::String, grid::AbstractGrid; kwargs...) VTKFileCollection(name::String, dh::DofHandler; kwargs...) -Create a paraview data file (.pvd) that can be used to -save multiple vtk file along with a time stamp. The keyword arguments +Create a paraview data file (.pvd) that can be used to +save multiple vtk file along with a time stamp. The keyword arguments are forwarded to `WriteVTK.paraview_collection`. See [`addstep!`](@ref) for examples for how to use `VTKFileCollection`. @@ -81,7 +81,7 @@ Base.close(pvd::VTKFileCollection) = WriteVTK.vtk_save(pvd.pvd) """ addstep!(f::Function, pvd::VTKFileCollection, t::Real, [grid_or_dh]; kwargs...) -Add a step at time `t` by writing a `VTKFile` to `pvd`. +Add a step at time `t` by writing a `VTKFile` to `pvd`. The keyword arguments are forwarded to `WriteVTK.vtk_grid`. If required, a new grid can be used by supplying the grid or dofhandler as the last argument. Should be used in a do-block, e.g. @@ -90,7 +90,7 @@ filename = "myoutput" pvd = VTKFileCollection(filename, grid) for (n, t) in pairs(timevector) # Calculate, e.g., the solution `u` and the stress `σeff` - addstep!(pvd, t) do io + addstep!(pvd, t) do io write_cell_data(io, σeff, "Effective stress") write_solution(io, dh, u) end @@ -109,8 +109,8 @@ end """ addstep!(pvd::VTKFileCollection, vtk::VTKFile, t) -As an alternative to using the `addstep!(pvd, t) do` block, it is -also possible to add a specific `vtk` at time `t` to `pvd`. +As an alternative to using the `addstep!(pvd, t) do` block, it is +also possible to add a specific `vtk` at time `t` to `pvd`. Note that this will close the `vtk`. Example ```julia filename = "myoutput" @@ -239,13 +239,13 @@ end write_solution(vtk::VTKFile, dh::AbstractDofHandler, u::Vector, suffix="") Save the values at the nodes in the degree of freedom vector `u` to `vtk`. -Each field in `dh` will be saved separately, and `suffix` can be used to append +Each field in `dh` will be saved separately, and `suffix` can be used to append to the fieldname. -`u` can also contain tensorial values, but each entry in `u` must correspond to a -degree of freedom in `dh`, see [`write_node_data`](@ref write_node_data) for details. -Use `write_node_data` directly when exporting values that are already -sorted by the nodes in the grid. +`u` can also contain tensorial values, but each entry in `u` must correspond to a +degree of freedom in `dh`, see [`write_node_data`](@ref write_node_data) for details. +Use `write_node_data` directly when exporting values that are already +sorted by the nodes in the grid. """ function write_solution(vtk::VTKFile, dh::AbstractDofHandler, u::Vector, suffix="") fieldnames = Ferrite.getfieldnames(dh) # all primary fields @@ -280,13 +280,13 @@ end """ write_node_data(vtk::VTKFile, nodedata::Vector{Real}, name) write_node_data(vtk::VTKFile, nodedata::Vector{<:AbstractTensor}, name) - + Write the `nodedata` that is ordered by the nodes in the grid to `vtk`. -When `nodedata` contains `Tensors.Vec`s, each component is exported. +When `nodedata` contains `Tensors.Vec`s, each component is exported. Two-dimensional vectors are padded with zeros. -When `nodedata` contains second order tensors, the index order, +When `nodedata` contains second order tensors, the index order, `[11, 22, 33, 23, 13, 12, 32, 31, 21]`, follows the default Voigt order in Tensors.jl. """ function write_node_data(vtk::VTKFile, nodedata, name) @@ -311,9 +311,9 @@ end write_cellset(vtk, grid::AbstractGrid, cellset::String) write_cellset(vtk, grid::AbstractGrid, cellsets::Union{AbstractVector{String},AbstractSet{String}) -Write all cell sets in the grid with name according to their keys and -celldata 1 if the cell is in the set, and 0 otherwise. It is also possible to -only export a single `cellset`, or multiple `cellsets`. +Write all cell sets in the grid with name according to their keys and +celldata 1 if the cell is in the set, and 0 otherwise. It is also possible to +only export a single `cellset`, or multiple `cellsets`. """ function write_cellset(vtk, grid::AbstractGrid, cellsets=keys(getcellsets(getgrid(vtk)))) z = zeros(getncells(grid)) @@ -332,7 +332,7 @@ write_cellset(vtk, grid::AbstractGrid, cellset::String) = write_cellset(vtk, gri Saves the dirichlet boundary conditions to a vtkfile. Values will have a 1 where bcs are active and 0 otherwise """ -function write_constraints(vtk, ch::ConstraintHandler) +function write_constraints(vtk, ch::ConstraintHandler) unique_fields = [] for dbc in ch.dbcs push!(unique_fields, dbc.field_name) diff --git a/src/FEValues/CellValues.jl b/src/FEValues/CellValues.jl index 5383381ed7..e550062d36 100644 --- a/src/FEValues/CellValues.jl +++ b/src/FEValues/CellValues.jl @@ -41,9 +41,9 @@ struct CellValues{FV, GM, QR, detT} <: AbstractCellValues qr::QR # QuadratureRule detJdV::detT # AbstractVector{<:Number} or Nothing end -function CellValues(::Type{T}, qr::QuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation; - update_gradients::Union{Bool,Nothing} = nothing, update_detJdV::Union{Bool,Nothing} = nothing) where T - +function CellValues(::Type{T}, qr::QuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation; + update_gradients::Union{Bool,Nothing} = nothing, update_detJdV::Union{Bool,Nothing} = nothing) where T + _update_detJdV = update_detJdV === nothing ? true : update_detJdV FunDiffOrder = update_gradients === nothing ? 1 : convert(Int, update_gradients) # Logic must change when supporting update_hessian kwargs GeoDiffOrder = max(required_geo_diff_order(mapping_type(ip_fun), FunDiffOrder), _update_detJdV) @@ -70,7 +70,7 @@ geometric_interpolation(cv::CellValues) = geometric_interpolation(cv.geo_mapping getdetJdV(cv::CellValues, q_point::Int) = cv.detJdV[q_point] getdetJdV(::CellValues{<:Any, <:Any, <:Any, Nothing}, ::Int) = throw(ArgumentError("detJdV is not saved in CellValues")) -# Accessors for function values +# Accessors for function values getnbasefunctions(cv::CellValues) = getnbasefunctions(cv.fun_values) function_interpolation(cv::CellValues) = function_interpolation(cv.fun_values) function_difforder(cv::CellValues) = function_difforder(cv.fun_values) @@ -81,7 +81,7 @@ shape_gradient_type(cv::CellValues) = shape_gradient_type(cv.fun_values) @propagate_inbounds shape_gradient(cv::CellValues, q_point::Int, i::Int) = shape_gradient(cv.fun_values, q_point, i) @propagate_inbounds shape_symmetric_gradient(cv::CellValues, q_point::Int, i::Int) = shape_symmetric_gradient(cv.fun_values, q_point, i) -# Access quadrature rule values +# Access quadrature rule values getnquadpoints(cv::CellValues) = getnquadpoints(cv.qr) @inline function _update_detJdV!(detJvec::AbstractVector, q_point::Int, w, mapping) @@ -99,7 +99,7 @@ function reinit!(cv::CellValues, cell::Union{AbstractCell, Nothing}, x::Abstract geo_mapping = cv.geo_mapping fun_values = cv.fun_values n_geom_basefuncs = getngeobasefunctions(geo_mapping) - + check_reinit_sdim_consistency(:CellValues, shape_gradient_type(cv), eltype(x)) if cell === nothing && !isa(mapping_type(fun_values), IdentityMapping) throw(ArgumentError("The cell::AbstractCell input is required to reinit! non-identity function mappings")) @@ -126,6 +126,6 @@ function Base.show(io::IO, d::MIME"text/plain", cv::CellValues) print(io, "CellValues(", vstr, ", rdim=$rdim, and sdim=$sdim): ") print(io, getnquadpoints(cv), " quadrature points") print(io, "\n Function interpolation: "); show(io, d, ip_fun) - print(io, "\nGeometric interpolation: "); + print(io, "\nGeometric interpolation: "); sdim === nothing ? show(io, d, ip_geo) : show(io, d, ip_geo^sdim) -end \ No newline at end of file +end diff --git a/src/FEValues/FacetValues.jl b/src/FEValues/FacetValues.jl index 4c8c81a866..e87ade1962 100644 --- a/src/FEValues/FacetValues.jl +++ b/src/FEValues/FacetValues.jl @@ -40,12 +40,12 @@ struct FacetValues{FV, GM, FQR, detT, nT, V_FV<:AbstractVector{FV}, V_GM<:Abstra current_facet::ScalarWrapper{Int} end -function FacetValues(::Type{T}, fqr::FacetQuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation{sdim} = default_geometric_interpolation(ip_fun); - update_gradients::Union{Bool,Nothing} = nothing) where {T,sdim} - +function FacetValues(::Type{T}, fqr::FacetQuadratureRule, ip_fun::Interpolation, ip_geo::VectorizedInterpolation{sdim} = default_geometric_interpolation(ip_fun); + update_gradients::Union{Bool,Nothing} = nothing) where {T,sdim} + FunDiffOrder = update_gradients === nothing ? 1 : convert(Int, update_gradients) # Logic must change when supporting update_hessian kwargs GeoDiffOrder = max(required_geo_diff_order(mapping_type(ip_fun), FunDiffOrder), 1) - # Not sure why the type-asserts are required here but not for CellValues, + # Not sure why the type-asserts are required here but not for CellValues, # but they solve the type-instability for FacetValues geo_mapping = [GeometryMapping{GeoDiffOrder}(T, ip_geo.ip, qr) for qr in fqr.face_rules]::Vector{<:GeometryMapping{GeoDiffOrder}} fun_values = [FunctionValues{FunDiffOrder}(T, ip_fun, qr, ip_geo) for qr in fqr.face_rules]::Vector{<:FunctionValues{FunDiffOrder}} @@ -104,7 +104,7 @@ getnormal(fv::FacetValues, qp::Int) = fv.normals[qp] nfacets(fv::FacetValues) = length(fv.geo_mapping) function set_current_facet!(fv::FacetValues, face_nr::Int) - # Checking face_nr before setting current_facet allows us to use @inbounds + # Checking face_nr before setting current_facet allows us to use @inbounds # when indexing by getcurrentfacet(fv) in other places! checkbounds(Bool, 1:nfacets(fv), face_nr) || throw(ArgumentError("Face index out of range.")) fv.current_facet[] = face_nr @@ -121,7 +121,7 @@ function reinit!(fv::FacetValues, cell::Union{AbstractCell, Nothing}, x::Abstrac if !checkbounds(Bool, x, 1:n_geom_basefuncs) || length(x) != n_geom_basefuncs throw_incompatible_coord_length(length(x), n_geom_basefuncs) end - + geo_mapping = get_geo_mapping(fv) fun_values = get_fun_values(fv) @@ -137,7 +137,7 @@ function reinit!(fv::FacetValues, cell::Union{AbstractCell, Nothing}, x::Abstrac detJ = norm(weight_norm) detJ > 0.0 || throw_detJ_not_pos(detJ) @inbounds fv.detJdV[q_point] = detJ * w - @inbounds fv.normals[q_point] = weight_norm / norm(weight_norm) + @inbounds fv.normals[q_point] = weight_norm / norm(weight_norm) apply_mapping!(fun_values, q_point, mapping, cell) end end diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index 6a21ea5d75..0e664dacea 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -28,8 +28,8 @@ typeof_dNdξ(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{s """ FunctionValues{DiffOrder}(::Type{T}, ip_fun, qr::QuadratureRule, ip_geo::VectorizedInterpolation) -Create a `FunctionValues` object containing the shape values and gradients (up to order `DiffOrder`) -for both the reference cell (precalculated) and the real cell (updated in `reinit!`). +Create a `FunctionValues` object containing the shape values and gradients (up to order `DiffOrder`) +for both the reference cell (precalculated) and the real cell (updated in `reinit!`). """ FunctionValues @@ -49,7 +49,7 @@ end function FunctionValues{DiffOrder}(::Type{T}, ip::Interpolation, qr::QuadratureRule, ip_geo::VectorizedInterpolation) where {DiffOrder, T} n_shape = getnbasefunctions(ip) n_qpoints = getnquadpoints(qr) - + Nξ = zeros(typeof_N(T, ip, ip_geo), n_shape, n_qpoints) Nx = isa(mapping_type(ip), IdentityMapping) ? Nξ : similar(Nξ) @@ -110,21 +110,21 @@ function check_reinit_sdim_consistency(valname, ::Val{sdim_val}, ::Val{sdim_x}) throw(ArgumentError("The $valname (sdim=$sdim_val) and coordinates (sdim=$sdim_x) have different spatial dimensions.")) end -# Mapping types -struct IdentityMapping end +# Mapping types +struct IdentityMapping end # Not yet implemented: # struct CovariantPiolaMapping end # PR798 # struct ContravariantPiolaMapping end # PR798 -# struct DoubleCovariantPiolaMapping end -# struct DoubleContravariantPiolaMapping end +# struct DoubleCovariantPiolaMapping end +# struct DoubleContravariantPiolaMapping end mapping_type(fv::FunctionValues) = mapping_type(fv.ip) """ required_geo_diff_order(fun_mapping, fun_diff_order::Int) -Return the required order of geometric derivatives to map -the function values and gradients from the reference cell +Return the required order of geometric derivatives to map +the function values and gradients from the reference cell to the physical cell geometry. """ required_geo_diff_order(::IdentityMapping, fun_diff_order::Int) = fun_diff_order @@ -167,6 +167,6 @@ end return nothing end -# TODO in PR798, apply_mapping! for +# TODO in PR798, apply_mapping! for # * CovariantPiolaMapping # * ContravariantPiolaMapping diff --git a/src/FEValues/GeometryMapping.jl b/src/FEValues/GeometryMapping.jl index b6346b7550..e017cae2ce 100644 --- a/src/FEValues/GeometryMapping.jl +++ b/src/FEValues/GeometryMapping.jl @@ -1,9 +1,9 @@ """ MappingValues(J, H) -The mapping values are calculated based on a +The mapping values are calculated based on a `geometric_mapping::GeometryMapping` along with the cell coordinates, -and the stored jacobian, `J`, and potentially hessian, `H`, are +and the stored jacobian, `J`, and potentially hessian, `H`, are used when mapping the `FunctionValues` to the current cell during `reinit!`. """ MappingValues @@ -13,14 +13,14 @@ struct MappingValues{JT, HT} H::HT # dJ/dξ # Hessian end -@inline getjacobian(mv::MappingValues{<:Union{AbstractTensor, SMatrix}}) = mv.J +@inline getjacobian(mv::MappingValues{<:Union{AbstractTensor, SMatrix}}) = mv.J # @inline gethessian(mv::MappingValues{<:Any,<:AbstractTensor}) = mv.H # PR798 """ GeometryMapping{DiffOrder}(::Type{T}, ip_geo, qr::QuadratureRule) -Create a `GeometryMapping` object which contains the geometric +Create a `GeometryMapping` object which contains the geometric * shape values * gradient values (if DiffOrder ≥ 1) @@ -31,7 +31,7 @@ Create a `GeometryMapping` object which contains the geometric GeometryMapping struct GeometryMapping{DiffOrder, IP, M_t, dMdξ_t, d2Mdξ2_t} - ip::IP # ::Interpolation Geometric interpolation + ip::IP # ::Interpolation Geometric interpolation M::M_t # ::AbstractMatrix{<:Number} Values of geometric shape functions dMdξ::dMdξ_t # ::AbstractMatrix{<:Vec} Gradients of geometric shape functions in ref-domain d2Mdξ2::d2Mdξ2_t # ::AbstractMatrix{<:Tensor{2}} Hessians of geometric shape functions in ref-domain @@ -47,8 +47,8 @@ struct GeometryMapping{DiffOrder, IP, M_t, dMdξ_t, d2Mdξ2_t} return new{1, IP, M_t, dMdξ_t, Nothing}(ip, M, dMdξ, nothing) end #= function GeometryMapping( - ip::IP, M::M_t, dMdξ::dMdξ_t, d2Mdξ2::d2Mdξ2_t) where - {IP <: ScalarInterpolation, M_t<:AbstractMatrix{<:Number}, + ip::IP, M::M_t, dMdξ::dMdξ_t, d2Mdξ2::d2Mdξ2_t) where + {IP <: ScalarInterpolation, M_t<:AbstractMatrix{<:Number}, dMdξ_t <: AbstractMatrix{<:Vec}, d2Mdξ2_t <: AbstractMatrix{<:Tensor{2}}} return new{2, IP, M_t, dMdξ_t, d2Mdξ2_t}(ip, M, dMdξ, d2Mdξ2) end =# # PR798 @@ -63,7 +63,7 @@ end function GeometryMapping{1}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRule) where T n_shape = getnbasefunctions(ip) n_qpoints = getnquadpoints(qr) - + M = zeros(T, n_shape, n_qpoints) dMdξ = zeros(Vec{getdim(ip),T}, n_shape, n_qpoints) @@ -74,7 +74,7 @@ end #= function GeometryMapping{2}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRule) where T n_shape = getnbasefunctions(ip) n_qpoints = getnquadpoints(qr) - + M = zeros(T, n_shape, n_qpoints) dMdξ = zeros(Vec{getdim(ip),T}, n_shape, n_qpoints) d2Mdξ2 = zeros(Tensor{2,getdim(ip),T}, n_shape, n_qpoints) diff --git a/src/FEValues/InterfaceValues.jl b/src/FEValues/InterfaceValues.jl index e2a9b307da..661a25aa36 100644 --- a/src/FEValues/InterfaceValues.jl +++ b/src/FEValues/InterfaceValues.jl @@ -134,11 +134,11 @@ function reinit!( transform_interface_points!(quad_points_b, quad_points_a, interface_transformation) # TODO: This is the bottleneck, cache it? @assert length(quad_points_a) <= length(quad_points_b) - + # Re-evaluate shape functions in the transformed quadrature points precompute_values!(get_fun_values(iv.there), quad_points_b) precompute_values!(get_geo_mapping(iv.there), quad_points_b) - + # reinit! the "there" side reinit!(iv.there, cell_there, coords_there, facet_there) return iv @@ -485,7 +485,7 @@ end end @inline function _get_transformation_matrix(::NTuple{N,Int}, ::InterfaceOrientationInfo) where N - throw(ArgumentError("transformation is not implemented")) + throw(ArgumentError("transformation is not implemented")) end @doc raw""" diff --git a/src/FEValues/common_values.jl b/src/FEValues/common_values.jl index 93704d83a9..c7a694a748 100644 --- a/src/FEValues/common_values.jl +++ b/src/FEValues/common_values.jl @@ -33,14 +33,14 @@ end Update the `CellValues`/`FacetValues` object for a cell or face with coordinates `x`. The derivatives of the shape functions, and the new integration weights are computed. -For interpolations with non-identity mappings, the current `cell` is also required. +For interpolations with non-identity mappings, the current `cell` is also required. """ reinit! """ getnquadpoints(fe_v::AbstractValues) -Return the number of quadrature points. For `FacetValues`, +Return the number of quadrature points. For `FacetValues`, this is the number for the current face. """ function getnquadpoints end @@ -71,7 +71,7 @@ shape_value(fe_v::AbstractValues, q_point::Int, base_function::Int) """ geometric_value(fe_v::AbstractValues, q_point, base_function::Int) -Return the value of the geometric shape function `base_function` evaluated in +Return the value of the geometric shape function `base_function` evaluated in quadrature point `q_point`. """ geometric_value(fe_v::AbstractValues, q_point::Int, base_function::Int) @@ -292,7 +292,7 @@ function spatial_coordinate(fe_v::AbstractValues, q_point::Int, x::AbstractVecto end -# Utility functions used by GeometryMapping, FunctionValues +# Utility functions used by GeometryMapping, FunctionValues _copy_or_nothing(x) = copy(x) _copy_or_nothing(::Nothing) = nothing @@ -315,4 +315,3 @@ function shape_hessians_gradients_and_values!(hessians::AbstractMatrix, gradient end end =# - diff --git a/src/Grid/coloring.jl b/src/Grid/coloring.jl index 9b4e18ffd8..7c41700c82 100644 --- a/src/Grid/coloring.jl +++ b/src/Grid/coloring.jl @@ -68,7 +68,7 @@ end # See Appendix A in https://www.math.colostate.edu/%7Ebangerth/publications/2013-pattern.pdf function workstream_coloring(incidence_matrix, cellset) - + if length(cellset) == 0 return Vector{Int}[] elseif length(cellset) == 1 diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index c40bf37140..251f081c77 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -16,7 +16,7 @@ Node(x::NTuple{dim,T}) where {dim,T} = Node(Vec{dim,T}(x)) """ get_node_coordinate(::Node) - + Get the value of the node coordinate. """ get_node_coordinate(n::Node) = n.x @@ -58,7 +58,7 @@ nnodes( c::AbstractCell) = length(get_node_ids(c)) Ferrite.vertices(::AbstractCell) Returns a tuple with the node indices (of the nodes in a grid) for each vertex in a given cell. -This function induces the [`VertexIndex`](@ref), where the second index +This function induces the [`VertexIndex`](@ref), where the second index corresponds to the local index into this tuple. """ vertices(::AbstractCell) @@ -67,7 +67,7 @@ vertices(::AbstractCell) Ferrite.edges(::AbstractCell) Returns a tuple of 2-tuples containing the ordered node indices (of the nodes in a grid) corresponding to -the vertices that define an *oriented edge*. This function induces the +the vertices that define an *oriented edge*. This function induces the [`EdgeIndex`](@ref), where the second index corresponds to the local index into this tuple. Note that the vertices are sufficient to define an edge uniquely. @@ -91,7 +91,7 @@ reference_faces(::AbstractRefShape) Ferrite.faces(::AbstractCell) Returns a tuple of n-tuples containing the ordered node indices (of the nodes in a grid) corresponding to -the vertices that define an *oriented face*. This function induces the +the vertices that define an *oriented face*. This function induces the [`FaceIndex`](@ref), where the second index corresponds to the local index into this tuple. An *oriented face* is a face with the first node having the local index and the other @@ -105,7 +105,7 @@ faces(::AbstractCell) Ferrite.facets(::AbstractCell) Returns a tuple of n-tuples containing the ordered node indices (of the nodes in a grid) corresponding to -the vertices that define an oriented facet. This function induces the +the vertices that define an oriented facet. This function induces the [`FacetIndex`](@ref), where the second index corresponds to the local index into this tuple. See also [`vertices`](@ref), [`edges`](@ref), and [`faces`](@ref) @@ -120,13 +120,13 @@ facets(::AbstractCell) Ferrite.reference_facets(::Type{<:AbstractRefShape}) Returns a tuple of n-tuples containing the ordered local node indices corresponding to -the vertices that define an oriented facet. +the vertices that define an oriented facet. See also [`reference_vertices`](@ref), [`reference_edges`](@ref), and [`reference_faces`](@ref) """ reference_facets(::Type{<:AbstractRefShape}) -@inline reference_facets(refshape::Type{<:AbstractRefShape{1}}) = map(i -> (i,), reference_vertices(refshape)) +@inline reference_facets(refshape::Type{<:AbstractRefShape{1}}) = map(i -> (i,), reference_vertices(refshape)) @inline reference_facets(refshape::Type{<:AbstractRefShape{2}}) = reference_edges(refshape) @inline reference_facets(refshape::Type{<:AbstractRefShape{3}}) = reference_faces(refshape) @@ -146,8 +146,8 @@ Default implementation: `c.nodes`. """ get_node_ids(c::AbstractCell) = c.nodes -# Default implementations of = vertices/edges/faces that work as long as get_node_ids -# and `reference_` are correctly implemented for the cell / reference shape. +# Default implementations of = vertices/edges/faces that work as long as get_node_ids +# and `reference_` are correctly implemented for the cell / reference shape. function vertices(c::Ferrite.AbstractCell{RefShape}) where RefShape ns = Ferrite.get_node_ids(c) @@ -169,9 +169,9 @@ function faces(c::Ferrite.AbstractCell{RefShape}) where RefShape end # RefLine (refdim = 1) -reference_vertices(::Type{RefLine}) = (1, 2) +reference_vertices(::Type{RefLine}) = (1, 2) reference_edges(::Type{RefLine}) = ((1, 2),) # e1 -reference_faces(::Type{RefLine}) = () # - +reference_faces(::Type{RefLine}) = () # - # RefTriangle (refdim = 2) reference_vertices(::Type{RefTriangle}) = (1, 2, 3) @@ -195,7 +195,7 @@ function reference_edges(::Type{RefHexahedron}) (7, 8), (8, 5), (1, 5), (2, 6), (3, 7), (4, 8)) # e7 ... e12 end function reference_faces(::Type{RefHexahedron}) - return ((1, 4, 3, 2), (1, 2, 6, 5), (2, 3, 7, 6), # f1, f2, f3 + return ((1, 4, 3, 2), (1, 2, 6, 5), (2, 3, 7, 6), # f1, f2, f3 (3, 4, 8, 7), (1, 5, 8, 4), (5, 6, 7, 8)) # f4, f5, f6 end @@ -206,7 +206,7 @@ function reference_edges(::Type{RefPrism}) (3, 6), (4, 5), (4, 6), (6, 5)) # e6, e7, e8, e9 end function reference_faces(::Type{RefPrism}) - return ((1, 3, 2), (1, 2, 5, 4), (3, 1, 4, 6), # f1, f2, f3 + return ((1, 3, 2), (1, 2, 5, 4), (3, 1, 4, 6), # f1, f2, f3 (2, 3, 6, 5), (4, 5, 6)) # f4, f5 end @@ -315,7 +315,7 @@ function Grid(cells::Vector{C}, facesets=nothing, # deprecated vertexsets::Dict{String, <:AbstractVecOrSet{VertexIndex}}=Dict{String,OrderedSet{VertexIndex}}(), boundary_matrix = nothing) where {dim,C,T} - if facesets !== nothing + if facesets !== nothing if isempty(facetsets) @warn "facesets in Grid is deprecated, use facetsets instead" maxlog=1 for (key, set) in facesets @@ -362,13 +362,13 @@ getsdim(::AbstractGrid{sdim}) where sdim = sdim """ get_reference_dimension(grid::AbstractGrid) -> Union{Int, Symbol} -Get information about the reference dimensions of the cells in the grid. -If all cells have the same reference dimension, `rdim::Int` is returned. +Get information about the reference dimensions of the cells in the grid. +If all cells have the same reference dimension, `rdim::Int` is returned. For grids with mixed reference dimensions, `:mixed` is returned. Used internally to dispatch facet-calls to the correct entity when `rdim isa Int`. """ get_reference_dimension(g::AbstractGrid) = _get_reference_dimension(getcells(g)) -_get_reference_dimension(::AbstractVector{C}) where C <: AbstractCell{<:AbstractRefShape{rdim}} where rdim = rdim # Fast path for single rdim inferable from eltype +_get_reference_dimension(::AbstractVector{C}) where C <: AbstractCell{<:AbstractRefShape{rdim}} where rdim = rdim # Fast path for single rdim inferable from eltype function _get_reference_dimension(cells::AbstractVector{<:AbstractCell}) # Could make fast-path for eltype being union of cells with different rdims, but @KristofferC recommends against that, # https://discourse.julialang.org/t/iterating-through-types-of-a-union-in-a-type-stable-manner/58285/3 @@ -413,7 +413,7 @@ to a Node. "Returns the number of nodes in the grid." @inline getnnodes(grid::AbstractGrid) = length(grid.nodes) "Returns the number of nodes of the `i`-th cell." -function nnodes_per_cell(grid::AbstractGrid) +function nnodes_per_cell(grid::AbstractGrid) if !isconcretetype(getcelltype(grid)) error("There are different celltypes in the `grid`. Use `nnodes_per_cell(grid, cellid::Int)` instead") end @@ -491,7 +491,7 @@ end """ getcoordinates(grid::AbstractGrid, idx::Union{Int,CellIndex}) getcoordinates(cache::CellCache) - + Get a vector with the coordinates of the cell corresponding to `idx` or `cache` """ @inline function getcoordinates(grid::AbstractGrid, idx::Int) @@ -506,7 +506,7 @@ end """ getcoordinates!(x::Vector{<:Vec}, grid::AbstractGrid, idx::Union{Int,CellIndex}) getcoordinates!(x::Vector{<:Vec}, grid::AbstractGrid, cell::AbstractCell) - + Mutate `x` to the coordinates of the cell corresponding to `idx` or `cell`. """ @inline function getcoordinates!(x::Vector{Vec{dim,T}}, grid::Ferrite.AbstractGrid, cell::Ferrite.AbstractCell) where {dim,T} @@ -516,7 +516,7 @@ Mutate `x` to the coordinates of the cell corresponding to `idx` or `cell`. end return x end -@inline function getcoordinates!(x::Vector{Vec{dim,T}}, grid::Ferrite.AbstractGrid, cellid::Int) where {dim,T} +@inline function getcoordinates!(x::Vector{Vec{dim,T}}, grid::Ferrite.AbstractGrid, cellid::Int) where {dim,T} cell = getcells(grid, cellid) getcoordinates!(x, grid, cell) end @@ -524,7 +524,7 @@ end """ get_node_coordinate(grid::AbstractGrid, n::Int) - + Return the coordinate of the `n`th node in `grid` """ get_node_coordinate(grid, n) = get_node_coordinate(getnodes(grid, n)) @@ -565,12 +565,12 @@ boundaryfunction(::Type{VertexIndex}) = vertices boundaryfunction(::Type{FacetIndex}) = facets for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex, :FacetIndex) - @eval begin + @eval begin #Constructor ($INDEX)(a::Int, b::Int) = ($INDEX)((a,b)) Base.getindex(I::($INDEX), i::Int) = I.idx[i] - + #To be able to do a,b = faceidx Base.iterate(I::($INDEX), state::Int=1) = (state==3) ? nothing : (I[state], state+1) @@ -610,8 +610,8 @@ end """ SurfaceOrientationInfo -Orientation information for 2D entities. Such an entity can be -possibly flipped (i.e. the defining vertex order is reverse to the +Orientation information for 2D entities. Such an entity can be +possibly flipped (i.e. the defining vertex order is reverse to the spanning vertex order) and the vertices can be rotated against each other. Take for example the faces ``` @@ -625,7 +625,7 @@ which are rotated against each other by 90° (shift index is 1) or the faces | A | | B | 4---3 3---4 ``` -which are flipped against each other. Any combination of these can happen. +which are flipped against each other. Any combination of these can happen. The combination to map this local face to the defining face is encoded with this data structure via ``rotate \\circ flip`` where the rotation is indiced by the shift index. @@ -652,7 +652,7 @@ is called *regular*, indicated by `flipped=false`, while the oriented path ``` is called *inverted*, indicated by `flipped=true`. -2D entities can be flipped (i.e. the defining vertex order is reverse to the +2D entities can be flipped (i.e. the defining vertex order is reverse to the spanning vertex order) and the vertices can be rotated against each other. The reference entity is a one with it's first node is the lowest index vertex diff --git a/src/Grid/grid_generators.jl b/src/Grid/grid_generators.jl index b1bb88c932..42eca02dc1 100644 --- a/src/Grid/grid_generators.jl +++ b/src/Grid/grid_generators.jl @@ -217,7 +217,7 @@ function generate_grid(::Type{Hexahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Ve foreach(s -> sort!(s, by = x -> x.idx), values(facetsets)) return Grid(cells, nodes, facetsets=facetsets) -end +end # Wedge function generate_grid(::Type{Wedge}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T} @@ -238,9 +238,9 @@ function generate_grid(::Type{Wedge}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}( node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y, n_nodes_z)) cells = Wedge[] for k in 1:nel_z, j in 1:nel_y, i in 1:nel_x - push!(cells, Wedge((node_array[i,j,k], node_array[i+1,j,k], node_array[i,j+1,k], + push!(cells, Wedge((node_array[i,j,k], node_array[i+1,j,k], node_array[i,j+1,k], node_array[i,j,k+1], node_array[i+1,j,k+1], node_array[i,j+1,k+1]))) # ◺ - push!(cells, Wedge((node_array[i+1,j,k], node_array[i+1,j+1,k], node_array[i,j+1,k], + push!(cells, Wedge((node_array[i+1,j,k], node_array[i+1,j+1,k], node_array[i,j+1,k], node_array[i+1,j,k+1], node_array[i+1,j+1,k+1], node_array[i,j+1,k+1]))) # ◹ end @@ -249,7 +249,7 @@ function generate_grid(::Type{Wedge}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}( @views le = map(x -> FacetIndex(x,3), c_nxyz[1, 1, :, :][:]) @views ri = map(x -> FacetIndex(x,2), c_nxyz[2, end, :, :][:]) - @views fr = map(x -> FacetIndex(x,2), c_nxyz[1, :, 1, :][:]) + @views fr = map(x -> FacetIndex(x,2), c_nxyz[1, :, 1, :][:]) @views ba = map(x -> FacetIndex(x,4), c_nxyz[2, :, end, :][:]) @views bo = [map(x -> FacetIndex(x,1), c_nxyz[1, :, :, 1][:]) ; map(x -> FacetIndex(x,1), c_nxyz[2, :, :, 1][:])] @views to = [map(x -> FacetIndex(x,5), c_nxyz[1, :, :, end][:]) ; map(x -> FacetIndex(x,5), c_nxyz[2, :, :, end][:])] @@ -265,7 +265,7 @@ function generate_grid(::Type{Wedge}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}( foreach(s -> sort!(s, by = x -> x.idx), values(facetsets)) return Grid(cells, nodes, facetsets=facetsets) -end +end #Pyramid function generate_grid(::Type{Pyramid}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T} @@ -288,20 +288,20 @@ function generate_grid(::Type{Pyramid}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3 midy = 0.5(coords_y[j+1] + coords_y[j]) midz = 0.5(coords_z[k+1] + coords_z[k]) push!(nodes, Node((midx, midy, midz))) - end + end # Generate cells node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y, n_nodes_z)) cells = Pyramid[] - midnodecounter = n_nodes_x*n_nodes_y*n_nodes_z + midnodecounter = n_nodes_x*n_nodes_y*n_nodes_z for k in 1:nel_z, j in 1:nel_y, i in 1:nel_x midnodecounter += 1 pyramid1 = Pyramid((node_array[i,j,k], node_array[i+1,j,k], node_array[i,j+1,k], node_array[i+1,j+1,k], midnodecounter )) # bottom - pyramid2 = Pyramid((node_array[i,j,k], node_array[i,j,k+1], node_array[i+1,j,k], node_array[i+1,j,k+1], midnodecounter )) # front + pyramid2 = Pyramid((node_array[i,j,k], node_array[i,j,k+1], node_array[i+1,j,k], node_array[i+1,j,k+1], midnodecounter )) # front pyramid3 = Pyramid((node_array[i+1,j,k], node_array[i+1,j,k+1], node_array[i+1,j+1,k], node_array[i+1,j+1,k+1], midnodecounter )) # right pyramid4 = Pyramid((node_array[i,j+1,k], node_array[i+1,j+1,k], node_array[i,j+1,k+1], node_array[i+1,j+1,k+1], midnodecounter )) # back - pyramid5 = Pyramid((node_array[i,j,k], node_array[i,j+1,k], node_array[i,j,k+1], node_array[i,j+1,k+1], midnodecounter )) # left - pyramid6 = Pyramid((node_array[i,j,k+1], node_array[i,j+1,k+1], node_array[i+1,j,k+1], node_array[i+1,j+1,k+1], midnodecounter )) # top + pyramid5 = Pyramid((node_array[i,j,k], node_array[i,j+1,k], node_array[i,j,k+1], node_array[i,j+1,k+1], midnodecounter )) # left + pyramid6 = Pyramid((node_array[i,j,k+1], node_array[i,j+1,k+1], node_array[i+1,j,k+1], node_array[i+1,j+1,k+1], midnodecounter )) # top push!(cells, pyramid1, pyramid2, pyramid3, pyramid4, pyramid5, pyramid6) end @@ -311,9 +311,9 @@ function generate_grid(::Type{Pyramid}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3 @views le = map(x -> FacetIndex(x,1), c_nxyz[5, 1, :, :][:]) @views ri = map(x -> FacetIndex(x,1), c_nxyz[3, end, :, :][:]) - @views fr = map(x -> FacetIndex(x,1), c_nxyz[2, :, 1, :][:]) + @views fr = map(x -> FacetIndex(x,1), c_nxyz[2, :, 1, :][:]) @views ba = map(x -> FacetIndex(x,1), c_nxyz[4, :, end, :][:]) - @views bo = map(x -> FacetIndex(x,1), c_nxyz[1, :, :, 1][:]) + @views bo = map(x -> FacetIndex(x,1), c_nxyz[1, :, :, 1][:]) @views to = map(x -> FacetIndex(x,1), c_nxyz[6, :, :, end][:]) facetsets = Dict( @@ -327,7 +327,7 @@ function generate_grid(::Type{Pyramid}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3 foreach(s -> sort!(s, by = x -> x.idx), values(facetsets)) return Grid(cells, nodes, facetsets=facetsets) -end +end function Ferrite.generate_grid(::Type{SerendipityQuadraticHexahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T} nel_x = nel[1]; nel_y = nel[2]; nel_z = nel[3]; nel_tot = nel_x*nel_y*nel_z @@ -353,9 +353,9 @@ function Ferrite.generate_grid(::Type{SerendipityQuadraticHexahedron}, nel::NTup # Generate cells cells = SerendipityQuadraticHexahedron[] - for k in 1:2:2nel_z, j in 1:2:2nel_y, i in 1:2:2nel_x + for k in 1:2:2nel_z, j in 1:2:2nel_y, i in 1:2:2nel_x push!(cells, SerendipityQuadraticHexahedron(( - node_array[i,j,k], node_array[i+2,j,k], node_array[i+2,j+2,k], node_array[i,j+2,k], # vertices bot + node_array[i,j,k], node_array[i+2,j,k], node_array[i+2,j+2,k], node_array[i,j+2,k], # vertices bot node_array[i,j,k+2], node_array[i+2,j,k+2], node_array[i+2,j+2,k+2], node_array[i,j+2,k+2], # vertices top node_array[i+1,j,k], node_array[i+2,j+1,k], node_array[i+1,j+2,k], node_array[i,j+1,k], # edges horizontal bottom node_array[i+1,j,k+2], node_array[i+2,j+1,k+2], node_array[i+1,j+2,k+2], node_array[i,j+1,k+2], # edges horizontal top @@ -383,7 +383,7 @@ function Ferrite.generate_grid(::Type{SerendipityQuadraticHexahedron}, nel::NTup facetsets["top"] = OrderedSet{FacetIndex}(boundary[(1:length(cell_array[:,:,end][:])) .+ offset]); offset += length(cell_array[:,:,end][:]) foreach(s -> sort!(s, by = x -> x.idx), values(facetsets)) - return Grid(cells, nodes, facetsets=facetsets) + return Grid(cells, nodes, facetsets=facetsets) end # Triangle diff --git a/src/Grid/topology.jl b/src/Grid/topology.jl index 31560453f1..c6fbe9f831 100644 --- a/src/Grid/topology.jl +++ b/src/Grid/topology.jl @@ -53,7 +53,7 @@ The struct saves the highest dimensional neighborhood, i.e. if something is conn - `face_neighbor::Matrix{EntityNeighborhood,Int}`: `face_neighbor[cellid,local_face_id]` -> neighboring face - `vertex_neighbor::Matrix{EntityNeighborhood,Int}`: `vertex_neighbor[cellid,local_vertex_id]` -> neighboring vertex - `edge_neighbor::Matrix{EntityNeighborhood,Int}`: `edge_neighbor[cellid_local_vertex_id]` -> neighboring edge -- `face_skeleton::Union{Vector{FaceIndex}, Nothing}`: +- `face_skeleton::Union{Vector{FaceIndex}, Nothing}`: !!! note Currently mixed-dimensional queries do not work at the moment. They will be added back later. """ @@ -346,11 +346,11 @@ end """ _create_skeleton(neighborhood::Matrix{EntityNeighborhood{BI}}) where BI <: Union{FaceIndex, EdgeIndex, VertexIndex} -Materializes the skeleton from the `neighborhood` information by returning a `Vector{BI}` with `BI`s describing -the unique entities in the grid. +Materializes the skeleton from the `neighborhood` information by returning a `Vector{BI}` with `BI`s describing +the unique entities in the grid. -*Example:* With `BI=EdgeIndex`, and an edge between cells and 1 and 2, with vertices 2 and 5, could be described by either -`EdgeIndex(1, 2)` or `EdgeIndex(2, 4)`, but only one of these will be in the vector returned by this function. +*Example:* With `BI=EdgeIndex`, and an edge between cells and 1 and 2, with vertices 2 and 5, could be described by either +`EdgeIndex(1, 2)` or `EdgeIndex(2, 4)`, but only one of these will be in the vector returned by this function. """ function _create_skeleton(neighborhood::Matrix{EntityNeighborhood{BI}}) where BI <: Union{FaceIndex, EdgeIndex, VertexIndex} i = 1 @@ -367,9 +367,9 @@ end """ vertexskeleton(top::ExclusiveTopology, ::AbstractGrid) -> Vector{VertexIndex} -Materializes the skeleton from the `neighborhood` information by returning a `Vector{VertexIndex}` -describing the unique vertices in the grid. (One unique vertex may have multiple `VertexIndex`, but only -one is included in the returned `Vector`) +Materializes the skeleton from the `neighborhood` information by returning a `Vector{VertexIndex}` +describing the unique vertices in the grid. (One unique vertex may have multiple `VertexIndex`, but only +one is included in the returned `Vector`) """ function vertexskeleton(top::ExclusiveTopology, ::Union{AbstractGrid,Nothing}=nothing) if top.vertex_skeleton === nothing @@ -381,9 +381,9 @@ end """ edgeskeleton(top::ExclusiveTopology, ::AbstractGrid) -> Vector{EdgeIndex} -Materializes the skeleton from the `neighborhood` information by returning a `Vector{EdgeIndex}` -describing the unique edge in the grid. (One unique edge may have multiple `EdgeIndex`, but only -one is included in the returned `Vector`) +Materializes the skeleton from the `neighborhood` information by returning a `Vector{EdgeIndex}` +describing the unique edge in the grid. (One unique edge may have multiple `EdgeIndex`, but only +one is included in the returned `Vector`) """ function edgeskeleton(top::ExclusiveTopology, ::Union{AbstractGrid,Nothing}=nothing) if top.edge_skeleton === nothing @@ -395,9 +395,9 @@ end """ faceskeleton(top::ExclusiveTopology, ::AbstractGrid) -> Vector{FaceIndex} -Materializes the skeleton from the `neighborhood` information by returning a `Vector{FaceIndex}` -describing the unique faces in the grid. (One unique face may have multiple `FaceIndex`, but only -one is included in the returned `Vector`) +Materializes the skeleton from the `neighborhood` information by returning a `Vector{FaceIndex}` +describing the unique faces in the grid. (One unique face may have multiple `FaceIndex`, but only +one is included in the returned `Vector`) """ function faceskeleton(top::ExclusiveTopology, ::Union{AbstractGrid,Nothing}=nothing) if top.face_skeleton === nothing @@ -409,12 +409,12 @@ end """ facetskeleton(top::ExclusiveTopology, grid::AbstractGrid) -Materializes the skeleton from the `neighborhood` information by returning a `Vector{BI}` where -`BI <: Union{VertexIndex, EdgeIndex, FaceIndex}`. -It describes the unique facets in the grid, and allows for dimension-independent code in the case -that all cells have the same reference dimension. For cells with different reference dimensions, +Materializes the skeleton from the `neighborhood` information by returning a `Vector{BI}` where +`BI <: Union{VertexIndex, EdgeIndex, FaceIndex}`. +It describes the unique facets in the grid, and allows for dimension-independent code in the case +that all cells have the same reference dimension. For cells with different reference dimensions, [`Ferrite.vertexskeleton`](@ref), [`Ferrite.edgeskeleton`](@ref), or [`Ferrite.faceskeleton`](@ref) -must be used explicitly. +must be used explicitly. """ function facetskeleton(top::ExclusiveTopology, grid::AbstractGrid) rdim = get_reference_dimension(grid) diff --git a/src/Grid/utils.jl b/src/Grid/utils.jl index 4fc83a6df7..fe96df897c 100644 --- a/src/Grid/utils.jl +++ b/src/Grid/utils.jl @@ -11,7 +11,7 @@ Adds a cellset to the grid with key `name`. Cellsets are typically used to define subdomains of the problem, e.g. two materials in the computational domain. The `DofHandler` can construct different fields which live not on the whole domain, but rather on a cellset. `all=true` implies that `f(x)` must return `true` for all nodal coordinates `x` in the cell if the cell -should be added to the set, otherwise it suffices that `f(x)` returns `true` for one node. +should be added to the set, otherwise it suffices that `f(x)` returns `true` for one node. ```julia addcellset!(grid, "left", Set((1,3))) #add cells with id 1 and 3 to cellset left @@ -28,41 +28,41 @@ end """ addnodeset!(grid::AbstractGrid, name::String, nodeid::AbstractVecOrSet{Int}) - addnodeset!(grid::AbstractGrid, name::String, f::Function) + addnodeset!(grid::AbstractGrid, name::String, f::Function) -Adds a `nodeset::OrderedSet{Int}` to the `grid`'s `nodesets` with key `name`. Has the same interface as `addcellset`. +Adds a `nodeset::OrderedSet{Int}` to the `grid`'s `nodesets` with key `name`. Has the same interface as `addcellset`. However, instead of mapping a cell id to the `String` key, a set of node ids is returned. """ -addnodeset!(grid::AbstractGrid, name::String, nodeid::AbstractVecOrSet{Int}) = +addnodeset!(grid::AbstractGrid, name::String, nodeid::AbstractVecOrSet{Int}) = _addset!(grid, name, nodeid, getnodesets(grid)) -addnodeset!(grid::AbstractGrid, name::String, f::Function) = +addnodeset!(grid::AbstractGrid, name::String, f::Function) = _addset!(grid, name, create_nodeset(grid, f), getnodesets(grid)) """ addfacetset!(grid::AbstractGrid, name::String, faceid::AbstractVecOrSet{FacetIndex}) - addfacetset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) + addfacetset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) Adds a facetset to the grid with key `name`. A facetset maps a `String` key to a `OrderedSet` of tuples corresponding to `(global_cell_id, local_facet_id)`. Facetsets can be used to initialize `Dirichlet` boundary conditions for the `ConstraintHandler`. `all=true` implies that `f(x)` must return `true` for all nodal coordinates `x` on the facet if the facet -should be added to the set, otherwise it suffices that `f(x)` returns `true` for one node. +should be added to the set, otherwise it suffices that `f(x)` returns `true` for one node. ```julia addfacetset!(grid, "right", Set((FacetIndex(2,2), FacetIndex(4,2)))) #see grid manual example for reference addfacetset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0) #see incompressible elasticity example for reference ``` """ -addfacetset!(grid::AbstractGrid, name::String, set::AbstractVecOrSet{FacetIndex}) = +addfacetset!(grid::AbstractGrid, name::String, set::AbstractVecOrSet{FacetIndex}) = _addset!(grid, name, set, getfacetsets(grid)) -addfacetset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) = +addfacetset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) = _addset!(grid, name, create_facetset(grid, f; all=all), getfacetsets(grid)) """ addvertexset!(grid::AbstractGrid, name::String, faceid::AbstractVecOrSet{FaceIndex}) - addvertexset!(grid::AbstractGrid, name::String, f::Function) + addvertexset!(grid::AbstractGrid, name::String, f::Function) Adds a vertexset to the grid with key `name`. A vertexset maps a `String` key to a `OrderedSet` of tuples corresponding to `(global_cell_id, local_vertex_id)`. @@ -73,10 +73,10 @@ addvertexset!(grid, "right", Set((VertexIndex(2,2), VertexIndex(4,2)))) addvertexset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0) ``` """ -addvertexset!(grid::AbstractGrid, name::String, set::AbstractVecOrSet{VertexIndex}) = +addvertexset!(grid::AbstractGrid, name::String, set::AbstractVecOrSet{VertexIndex}) = _addset!(grid, name, set, getvertexsets(grid)) -addvertexset!(grid::AbstractGrid, name::String, f::Function) = +addvertexset!(grid::AbstractGrid, name::String, f::Function) = _addset!(grid, name, create_vertexset(grid, f; all=true), getvertexsets(grid)) function _addset!(grid::AbstractGrid, name::String, _set::AbstractVecOrSet, dict::Dict) @@ -133,7 +133,7 @@ function _create_set(f::Function, grid::AbstractGrid, ::Type{BI}; all=true) wher return set end -# Given a boundary index, for example EdgeIndex(2, 1), add this to `set`, as well as any other `EdgeIndex` in the grid +# Given a boundary index, for example EdgeIndex(2, 1), add this to `set`, as well as any other `EdgeIndex` in the grid # pointing to the same edge (i.e. indices belong to neighboring cells) function push_entity_instances!(set::OrderedSet{BI}, grid::AbstractGrid, top::ExclusiveTopology, entity::BI) where {BI <: BoundaryIndex} push!(set, entity) # Add the given entity @@ -150,11 +150,11 @@ function push_entity_instances!(set::OrderedSet{BI}, grid::AbstractGrid, top::Ex return set end -# Create a `OrderedSet{BI}` whose entities are a subset of facets which do not have neighbors, i.e. that are on the boundary. -# Note that this may include entities from cells incident to the facet, e.g. +# Create a `OrderedSet{BI}` whose entities are a subset of facets which do not have neighbors, i.e. that are on the boundary. +# Note that this may include entities from cells incident to the facet, e.g. # ____ consider the case of a vertex boundary set, with f(x) only true on the right side. Then also the VertexIndex -# |\ | belong to the lower left cell, in its lower right corner, is on the boundary, so this should be added too. -# |_\| That is done by the `push_entity_instances!` function. +# |\ | belong to the lower left cell, in its lower right corner, is on the boundary, so this should be added too. +# |_\| That is done by the `push_entity_instances!` function. function _create_boundaryset(f::Function, grid::AbstractGrid, top::ExclusiveTopology, ::Type{BI}; all = true) where {BI <: BoundaryIndex} # Function barrier as get_facet_facet_neighborhood is not always type stable function _makeset(ff_nh) @@ -193,14 +193,14 @@ function create_cellset(grid::AbstractGrid, f::Function; all::Bool=true) end pass && push!(cells, i) end - return cells + return cells end function create_nodeset(grid::AbstractGrid, f::Function) nodes = OrderedSet{Int}() for (i, n) in pairs(getnodes(grid)) f(get_node_coordinate(n)) && push!(nodes, i) end - return nodes + return nodes end create_vertexset(grid::AbstractGrid, f::Function; kwargs...) = _create_set(f, grid, VertexIndex; kwargs...) create_edgeset( grid::AbstractGrid, f::Function; kwargs...) = _create_set(f, grid, EdgeIndex; kwargs...) @@ -215,7 +215,7 @@ create_boundaryfacetset( grid::AbstractGrid, top::ExclusiveTopology, f::Function """ bounding_box(grid::AbstractGrid) -Computes the axis-aligned bounding box for a given grid, based on its node coordinates. +Computes the axis-aligned bounding box for a given grid, based on its node coordinates. Returns the minimum and maximum vertex coordinates of the bounding box. """ function bounding_box(grid::AbstractGrid{dim}) where {dim} diff --git a/src/PointEvalHandler.jl b/src/PointEvalHandler.jl index 7b7e2e14c3..304a425843 100644 --- a/src/PointEvalHandler.jl +++ b/src/PointEvalHandler.jl @@ -55,7 +55,7 @@ function _get_cellcoords(points::AbstractVector{Vec{dim,T}}, grid::AbstractGrid, # set up tree structure for finding nearest nodes to points kdtree = KDTree(reinterpret(Vec{dim,T}, getnodes(grid))) - nearest_nodes, _ = knn(kdtree, points, search_nneighbors, true) + nearest_nodes, _ = knn(kdtree, points, search_nneighbors, true) cells = Vector{Union{Nothing, Int}}(nothing, length(points)) local_coords = Vector{Union{Nothing, Vec{dim, T}}}(nothing, length(points)) diff --git a/src/Quadrature/gaussquad_prism_table.jl b/src/Quadrature/gaussquad_prism_table.jl index 36ec8b7662..ce626dcd35 100644 --- a/src/Quadrature/gaussquad_prism_table.jl +++ b/src/Quadrature/gaussquad_prism_table.jl @@ -1,6 +1,6 @@ # Symmetric quadrature rules takes from -# Witherden, Freddie D., and Peter E. Vincent. "On the identification of -# symmetric quadrature rules for finite element methods." Computers & +# Witherden, Freddie D., and Peter E. Vincent. "On the identification of +# symmetric quadrature rules for finite element methods." Computers & # Mathematics with Applications 69.10 (2015): 1232-1241. # Note that the original rule is defined on [-1,1]^3 while our reference prism is defined on [0,1]^3, hence we transform in the end. function _get_gauss_prismdata_polyquad(n::Int) diff --git a/src/Quadrature/gaussquad_pyramid_table.jl b/src/Quadrature/gaussquad_pyramid_table.jl index 043e03b87d..77f773d2c0 100644 --- a/src/Quadrature/gaussquad_pyramid_table.jl +++ b/src/Quadrature/gaussquad_pyramid_table.jl @@ -1,8 +1,8 @@ # Symmetric quadrature rules takes from -# Witherden, Freddie D., and Peter E. Vincent. "On the identification of -# symmetric quadrature rules for finite element methods." Computers & +# Witherden, Freddie D., and Peter E. Vincent. "On the identification of +# symmetric quadrature rules for finite element methods." Computers & # Mathematics with Applications 69.10 (2015): 1232-1241. -# TODO: Implement quadrature data from: +# TODO: Implement quadrature data from: # https://www.sciencedirect.com/science/article/pii/S0168874X1200203X?via%3Dihub#s0065 function _get_gauss_pyramiddata_polyquad(n::Int) if n == 1 @@ -85,7 +85,7 @@ function _get_gauss_pyramiddata_polyquad(n::Int) else throw(ArgumentError("unsupported order for prism polyquad integration")) end - # + # # The above quadrature rule is defined for a pyramid spanning [-1,1] × [-1,1] × [-1,1], with volume 8/3 and with 5th node in center. # The reference pyramid in ferrite spans [0,1] × [0,1] × [0,1], with volume 1/3 and with 5th node in corner. # Here we map thequadrature points to the pyramid defined in Ferrite. @@ -104,4 +104,4 @@ function _get_gauss_pyramiddata_polyquad(n::Int) end return xw -end \ No newline at end of file +end diff --git a/src/Quadrature/gaussquad_tri_table.jl b/src/Quadrature/gaussquad_tri_table.jl index d6af16fed4..a3ee7ad7d2 100644 --- a/src/Quadrature/gaussquad_tri_table.jl +++ b/src/Quadrature/gaussquad_tri_table.jl @@ -1,13 +1,13 @@ -# Order 1 to 8 heights points / wheights have been suggested in +# Order 1 to 8 heights points / wheights have been suggested in # Dunavant, D. A. (1985), High degree efficient symmetrical Gaussian quadrature # rules for the triangle. Int. J. Numer. Meth. Engng., 21: 1129–1148. doi: # 10.1002/nme.1620210612 # # Quadrature rules for orders 9 to 20 have been obtained using the -# basix.make_quadrature(basix.CellType.triangle, n) calls of the +# basix.make_quadrature(basix.CellType.triangle, n) calls of the # FEniCS / basix python package # -# see +# see # https://docs.fenicsproject.org/basix/main/python/_autosummary/basix.html?highlight=quadraturetype#basix.make_quadrature diff --git a/src/assembler.jl b/src/assembler.jl index 573f4437dc..dd4caa4b84 100644 --- a/src/assembler.jl +++ b/src/assembler.jl @@ -24,7 +24,7 @@ Create an `Assembler` object which can be used to assemble element contributions global sparse matrix. Use [`assemble!`](@ref) for each element, and [`finish_assemble`](@ref), to finalize the assembly and return the sparse matrix. -Note that giving a sparse matrix as input can be more efficient. See below and +Note that giving a sparse matrix as input can be more efficient. See below and as described in the [manual](@ref man-assembly). !!! note diff --git a/src/deprecations.jl b/src/deprecations.jl index 9fa590abad..fab94c791b 100644 --- a/src/deprecations.jl +++ b/src/deprecations.jl @@ -59,14 +59,14 @@ Base.@deprecate_binding Line3D Line Base.@deprecate_binding Quadrilateral3D Quadrilateral export Line2D, Line3D, Quadrilateral3D -using WriteVTK: vtk_grid +using WriteVTK: vtk_grid export vtk_grid # To give better error function WriteVTK.vtk_grid(::String, ::Union{AbstractGrid,AbstractDofHandler}; kwargs...) error(join(("The vtk interface has been updated in Ferrite v1.0.", "See https://github.com/Ferrite-FEM/Ferrite.jl/pull/679.", "Use VTKFile to open a vtk file, and the functions", - "write_solution, write_cell_data, and write_projection to save data."), + "write_solution, write_cell_data, and write_projection to save data."), "\n")) end @@ -354,7 +354,7 @@ end @deprecate transform! transform_coordinates! export addfaceset! # deprecated, export for backwards compatibility. -# Use warn to show for standard users. +# Use warn to show for standard users. function addfaceset!(grid::AbstractGrid, name, set::Union{Set{FaceIndex}, Vector{FaceIndex}}) @warn "addfaceset! is deprecated, use addfacetset! instead. Interpreting FaceIndex as FacetIndex" new_set = Set(FacetIndex(idx[1], idx[2]) for idx in set) diff --git a/src/exports.jl b/src/exports.jl index ebca1fdb5d..0e834f88a6 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -165,7 +165,7 @@ export write_node_data, VTKFileCollection, addstep!, - + # L2 Projection project, L2Projector, diff --git a/src/interpolations.jl b/src/interpolations.jl index 0ff442deb3..089211e38c 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -237,7 +237,7 @@ end shape_hessian_gradient_and_value(ip::Interpolation, ξ::Vec, i::Int) Optimized version combining the evaluation [`Ferrite.shape_value(::Interpolation)`](@ref), -[`Ferrite.shape_gradient(::Interpolation)`](@ref), and the gradient of the latter. +[`Ferrite.shape_gradient(::Interpolation)`](@ref), and the gradient of the latter. """ function shape_hessian_gradient_and_value(ip::Interpolation, ξ::Vec, i::Int) return hessian(x -> shape_value(ip, x, i), ξ, :all) @@ -252,7 +252,7 @@ and indices corresponding to the indices of a dof in [`vertices`](@ref), [`faces [`edges`](@ref). Only required for nodal interpolations. - + TODO: Separate nodal and non-nodal interpolations. """ reference_coordinates(::Interpolation) @@ -354,7 +354,7 @@ edge dofs are included here. The dofs appearing in the tuple must be continuous and increasing! The first dof must be the computed via "last edge interior dof index + 1", if face dofs exist. """ -facedof_interior_indices(::Interpolation) +facedof_interior_indices(::Interpolation) """ volumedof_interior_indices(ip::Interpolation) @@ -389,7 +389,7 @@ facetdof_indices(ip::InterpolationByDim{2}) = Ferrite.edgedof_indices(ip) facetdof_indices(ip::InterpolationByDim{1}) = Ferrite.vertexdof_indices(ip) facetdof_interior_indices(ip::InterpolationByDim{3}) = Ferrite.facedof_interior_indices(ip) facetdof_interior_indices(ip::InterpolationByDim{2}) = Ferrite.edgedof_interior_indices(ip) -facetdof_interior_indices(ip::InterpolationByDim{1}) = ntuple(_ -> (), nvertices(ip)) +facetdof_interior_indices(ip::InterpolationByDim{1}) = ntuple(_ -> (), nvertices(ip)) dirichlet_facetdof_indices(ip::InterpolationByDim{3}) = dirichlet_facedof_indices(ip) dirichlet_facetdof_indices(ip::InterpolationByDim{2}) = dirichlet_edgedof_indices(ip) dirichlet_facetdof_indices(ip::InterpolationByDim{1}) = dirichlet_vertexdof_indices(ip) @@ -1060,19 +1060,19 @@ end getnbasefunctions(::Lagrange{RefPrism,2}) = 18 facedof_indices(::Lagrange{RefPrism,2}) = ( - #Vertices| Edges | Face + #Vertices| Edges | Face (1,3,2 , 8,10,7 ), - (1,2,5,4, 7,11,13,9, 16), + (1,2,5,4, 7,11,13,9, 16), (3,1,4,6, 8,9,14,12, 17), (2,3,6,5, 10,12,15,11, 18), (4,5,6 , 13,15,14 ), ) facedof_interior_indices(::Lagrange{RefPrism,2}) = ( - #Vertices| Edges | Face - (), - (16,), - (17,), - (18,), + #Vertices| Edges | Face + (), + (16,), + (17,), + (18,), (), ) edgedof_indices(::Lagrange{RefPrism,2}) = ( @@ -1154,7 +1154,7 @@ end getnbasefunctions(::Lagrange{RefPyramid,1}) = 5 facedof_indices(::Lagrange{RefPyramid,1}) = ((1,3,4,2), (1,2,5), (1,5,3), (2,4,5), (3,5,4), ) edgedof_indices(::Lagrange{RefPyramid,1}) = ((1,2), (1,3), (1,5), (2,4), (2,5), (4,3), (3,5), (4,5)) - + function reference_coordinates(::Lagrange{RefPyramid,1}) return [Vec{3, Float64}((0.0, 0.0, 0.0)), Vec{3, Float64}((1.0, 0.0, 0.0)), @@ -1180,28 +1180,28 @@ end getnbasefunctions(::Lagrange{RefPyramid,2}) = 14 facedof_indices(::Lagrange{RefPyramid,2}) = ( - #Vertices | Edges | Face - (1,3,4,2, 7,11,9,6, 14), - (1,2,5 , 6,10,8 ), - (1,5,3 , 7,12,8 ), - (2,4,5 , 9,13,10 ), - (3,5,4 , 12,13,11 ), + #Vertices | Edges | Face + (1,3,4,2, 7,11,9,6, 14), + (1,2,5 , 6,10,8 ), + (1,5,3 , 7,12,8 ), + (2,4,5 , 9,13,10 ), + (3,5,4 , 12,13,11 ), ) facedof_interior_indices(::Lagrange{RefPyramid,2}) = ( - (14,), - (), - (), - (), + (14,), + (), + (), + (), (), ) edgedof_indices(::Lagrange{RefPyramid,2}) = ( - (1,2,6), - (1,3,7), - (1,5,8), - (2,4,9), - (2,5,10), - (4,3,11), - (3,5,12), + (1,2,6), + (1,3,7), + (1,5,8), + (2,4,9), + (2,5,10), + (4,3,11), + (3,5,12), (4,5,13) ) edgedof_interior_indices(::Lagrange{RefPyramid,2}) = ( @@ -1306,7 +1306,7 @@ struct Serendipity{shape, order, unused} <: ScalarInterpolation{shape,order} end end -# Note that the edgedofs for high order serendipity elements are defined in terms of integral moments, +# Note that the edgedofs for high order serendipity elements are defined in terms of integral moments, # so no permutation exists in general. See e.g. Scroggs et al. [2022] for an example. # adjust_dofs_during_distribution(::Serendipity{refshape, order}) where {refshape, order} = false adjust_dofs_during_distribution(::Serendipity{<:Any, 2}) = false @@ -1444,8 +1444,8 @@ end Classical non-conforming Crouzeix–Raviart element. For details we refer to the original paper: -M. Crouzeix and P. Raviart. "Conforming and nonconforming finite element -methods for solving the stationary Stokes equations I." ESAIM: Mathematical Modelling +M. Crouzeix and P. Raviart. "Conforming and nonconforming finite element +methods for solving the stationary Stokes equations I." ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique 7.R3 (1973): 33-75. """ struct CrouzeixRaviart{shape, order, unused} <: ScalarInterpolation{shape, order} @@ -1467,7 +1467,7 @@ function reference_coordinates(::CrouzeixRaviart) Vec{2, Float64}((0.5, 0.0))] end -function shape_value(ip::CrouzeixRaviart{RefTriangle, 1}, ξ::Vec{2}, i::Int) +function shape_value(ip::CrouzeixRaviart{RefTriangle, 1}, ξ::Vec{2}, i::Int) ξ_x = ξ[1] ξ_y = ξ[2] i == 1 && return 2*ξ_x + 2*ξ_y - 1 @@ -1550,7 +1550,7 @@ is_discontinuous(::Type{<:VectorizedInterpolation{<:Any, <:Any, <:Any, ip}}) whe Get the type of mapping from the reference cell to the real cell for an interpolation `ip`. Subtypes of `ScalarInterpolation` and `VectorizedInterpolation` return `IdentityMapping()`, but other non-scalar interpolations may request different -mapping types. +mapping types. """ function mapping_type end diff --git a/src/utils.jl b/src/utils.jl index c77469b357..b6350598f2 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -46,7 +46,7 @@ Base.copy(s::ScalarWrapper{T}) where {T} = ScalarWrapper{T}(copy(s.x)) convert_to_orderedset(set::AbstractVector{T}) where T = OrderedSet{T}(set) convert_to_orderedset(set::AbstractSet{T}) where T = convert(OrderedSet{T}, set) -function convert_to_orderedsets(namedsets::Dict{String, <: AbstractVecOrSet{T}}) where T +function convert_to_orderedsets(namedsets::Dict{String, <: AbstractVecOrSet{T}}) where T return Dict{String,OrderedSet{T}}(k => convert_to_orderedset(v) for (k,v) in namedsets) end convert_to_orderedsets(namedsets::Dict{String, <: OrderedSet}) = namedsets diff --git a/test/integration/test_simple_scalar_convergence.jl b/test/integration/test_simple_scalar_convergence.jl index 37d5f594fa..8b0ae5c265 100644 --- a/test/integration/test_simple_scalar_convergence.jl +++ b/test/integration/test_simple_scalar_convergence.jl @@ -102,7 +102,7 @@ function check_and_compute_convergence_norms(dh, u, cellvalues, testatol) ∇uₐₙₐ = gradient(x-> prod(cos, x*π/2), x) ∇uₐₚₚᵣₒₓ = function_gradient(cellvalues, q_point, uₑ) ∇L2norm += norm(∇uₐₙₐ-∇uₐₚₚᵣₒₓ)^2*dΩ - + # Pointwise convergence @test uₐₙₐ ≈ uₐₚₚᵣₒₓ atol=testatol end @@ -198,7 +198,7 @@ end dh, ch, cellvalues = ConvergenceTestHelper.setup_poisson_problem(grid, interpolation, interpolation_geo, qr) u = ConvergenceTestHelper.solve(dh, ch, cellvalues) L2₁, H1₁, _ = ConvergenceTestHelper.check_and_compute_convergence_norms(dh, u, cellvalues, 1e-2) - + # "Fine case" N₂ = 2*N₁ grid = generate_grid(geometry, ntuple(x->N₂, getdim(geometry))); diff --git a/test/runtests.jl b/test/runtests.jl index 39b13a6293..ba74624411 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,7 +21,7 @@ end const RUN_JET_TESTS = VERSION >= v"1.9" && isempty(VERSION.prerelease) if RUN_JET_TESTS - using Pkg: Pkg + using Pkg: Pkg Pkg.add("JET") using JET: @test_call else @@ -57,10 +57,10 @@ HAS_EXTENSIONS && include("blockarrays.jl") include("test_examples.jl") @test all(x -> isdefined(Ferrite, x), names(Ferrite)) # Test that all exported symbols are defined -#= See which is not defined if fails +#= See which is not defined if fails for name in names(Ferrite) isdefined(Ferrite, name) || @warn "Ferrite.$name is not defined but $name is exported" -end +end =# # Integration tests diff --git a/test/test_abstractgrid.jl b/test/test_abstractgrid.jl index 2c368a184e..f0896adc77 100644 --- a/test/test_abstractgrid.jl +++ b/test/test_abstractgrid.jl @@ -27,7 +27,7 @@ ip = Lagrange{RefQuadrilateral, 1}() qr = QuadratureRule{RefQuadrilateral}(2) cellvalues = CellValues(qr, ip); - + dhs = [DofHandler(grid) for grid in (subtype_grid, reference_grid)] u1 = Vector{Float64}(undef, 9) u2 = Vector{Float64}(undef, 9) diff --git a/test/test_apply_analytical.jl b/test/test_apply_analytical.jl index 45140a4109..16d6470200 100644 --- a/test/test_apply_analytical.jl +++ b/test/test_apply_analytical.jl @@ -109,8 +109,8 @@ apply_analytical!(a, dh, :u, f) @test sum(a)/length(a) ≈ num_udofs/(num_udofs+num_pdofs) - # If not super/subparametric, compare with ConstraintHandler and node set - if ip_order_u==ip_order_p==getcellorder(CT) + # If not super/subparametric, compare with ConstraintHandler and node set + if ip_order_u==ip_order_p==getcellorder(CT) fill!(a, 0) a_ch = copy(a) fp(x) = norm(x)^2 @@ -126,8 +126,8 @@ apply_analytical!(a, dh, :u, fu) apply_analytical!(a, dh, :p, fp) - @test a ≈ a_ch - end + @test a ≈ a_ch + end end end end diff --git a/test/test_apply_rhs.jl b/test/test_apply_rhs.jl index f56dbce4e9..bcd2517d22 100644 --- a/test/test_apply_rhs.jl +++ b/test/test_apply_rhs.jl @@ -3,15 +3,15 @@ function test_apply_rhs() ip = Lagrange{RefQuadrilateral,1}() qr = QuadratureRule{RefQuadrilateral}(2) cellvalues = CellValues(qr, ip) - + dh = DofHandler(grid) add!(dh, :u, ip) close!(dh) - + K = create_sparsity_pattern(dh) - + ch = ConstraintHandler(dh) - + ∂Ω = union(getfacetset.((grid,), ["left", "right"])...) dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0) add!(ch, dbc); @@ -19,32 +19,32 @@ function test_apply_rhs() ∂Ω = union(getfacetset.((grid,), ["top", "bottom"])...) dbc = Dirichlet(:u, ∂Ω, (x, t) -> 2) add!(ch, dbc); - + close!(ch) update!(ch, 0.0); - + function doassemble!( cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler, ) - + n_basefuncs = getnbasefunctions(cellvalues) Ke = zeros(n_basefuncs, n_basefuncs) fe = zeros(n_basefuncs) - + f = zeros(ndofs(dh)) assembler = start_assemble(K, f) - + @inbounds for cell in CellIterator(dh) fill!(Ke, 0) fill!(fe, 0) - + reinit!(cellvalues, cell) - + for q_point = 1:getnquadpoints(cellvalues) dΩ = getdetJdV(cellvalues, q_point) - + for i = 1:n_basefuncs v = shape_value(cellvalues, q_point, i) ∇v = shape_gradient(cellvalues, q_point, i) @@ -55,17 +55,17 @@ function test_apply_rhs() end end end - + assemble!(assembler, celldofs(cell), fe, Ke) end return K, f end - + K, f = doassemble!(cellvalues, K, dh) A = create_sparsity_pattern(dh) A, g = doassemble!(cellvalues, A, dh) rhsdata = get_rhs_data(ch, A) - + apply!(K, f, ch) apply!(A, ch) # need to apply bcs to A once apply_rhs!(rhsdata, g, ch) @@ -73,6 +73,6 @@ function test_apply_rhs() u₂ = A \ g return u₁, u₂ end - + u1, u2 = test_apply_rhs() @test u1 == u2 diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index dc3fca82d3..f4e06030b6 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -104,7 +104,7 @@ @test v == vc end end - # Test that qr and detJdV is copied as expected. + # Test that qr and detJdV is copied as expected. # Note that qr remain aliased, as defined by `copy(qr)=qr`, see quadrature.jl. for fname in (:qr, :detJdV) v = getfield(cv, fname) @@ -147,13 +147,13 @@ end fsv = FacetValues(qr_f, ip) fvv = FacetValues(qr_f, VectorizedInterpolation(ip)) fsv_embedded = FacetValues(qr_f, ip, ip^3) - + x, n = valid_coordinates_and_normals(ip) reinit!(csv, x) reinit!(cvv, x) reinit!(fsv, x, 1) reinit!(fvv, x, 1) - + # Wrong number of coordinates xx = [x; x] @test_throws ArgumentError reinit!(csv, xx) @@ -166,7 +166,7 @@ end @test_throws ArgumentError spatial_coordinate(fsv, qp, xx) @test_throws ArgumentError spatial_coordinate(fvv, qp, xx) - # Wrong dimension of coordinates + # Wrong dimension of coordinates @test_throws ArgumentError reinit!(csv_embedded, x) @test_throws ArgumentError reinit!(fsv_embedded, x, 1) @@ -306,7 +306,7 @@ end @testset "CellValues constructor entry points" begin qr = QuadratureRule{RefTriangle}(1) - + for fun_ip in (Lagrange{RefTriangle, 1}(), Lagrange{RefTriangle, 2}()^2) value_type(T) = fun_ip isa ScalarInterpolation ? T : Vec{2, T} grad_type(T) = fun_ip isa ScalarInterpolation ? Vec{2, T} : Tensor{2, 2, T, 4} @@ -356,13 +356,13 @@ end end @testset "CustomCellValues" begin - + @testset "SimpleCellValues" begin include(joinpath(@__DIR__, "../docs/src/topics/SimpleCellValues_literate.jl")) end - + @testset "TestCustomCellValues" begin - + struct TestCustomCellValues{CV<:CellValues} <: Ferrite.AbstractValues cv::CV end @@ -382,7 +382,7 @@ end ae = rand(getnbasefunctions(cv)) q_point = rand(1:getnquadpoints(cv)) cv_custom = TestCustomCellValues(cv) - for fun in (function_value, function_gradient, + for fun in (function_value, function_gradient, function_divergence, function_symmetric_gradient, function_curl) @test fun(cv_custom, q_point, ae) == fun(cv, q_point, ae) end diff --git a/test/test_constraints.jl b/test/test_constraints.jl index 94161ea981..8bf4e19fee 100644 --- a/test/test_constraints.jl +++ b/test/test_constraints.jl @@ -17,7 +17,7 @@ @test_throws ErrorException("components not sorted: [2, 1]") Dirichlet(:u, Γ, (x, t) -> 0, Int[2, 1]) @test_throws ErrorException("components not unique: [2, 2]") Dirichlet(:u, Γ, (x, t) -> 0, Int[2, 2]) @test_throws ErrorException("No dof prescribed for order 0 interpolations") add!(ch, Dirichlet(:z, Γ, (x, t) -> 0)) - for (s, v) in [(:s, :v), (:sd, :vd)] + for (s, v) in [(:s, :v), (:sd, :vd)] ## Scalar dbc = Dirichlet(s, Γ, (x, t) -> 0) add!(ch, dbc) @@ -169,7 +169,7 @@ end #Shell mesh edge bcs - nodes = [Node{3,Float64}(Vec(0.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,0.0,0.0)), + nodes = [Node{3,Float64}(Vec(0.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,1.0,0.0)), Node{3,Float64}(Vec(0.0,1.0,0.0)), Node{3,Float64}(Vec(2.0,0.0,0.0)), Node{3,Float64}(Vec(2.0,2.0,0.0))] @@ -217,7 +217,7 @@ end # This can be merged with the continuous test or removed. # Shell mesh edge bcs - nodes = [Node{3,Float64}(Vec(0.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,0.0,0.0)), + nodes = [Node{3,Float64}(Vec(0.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,1.0,0.0)), Node{3,Float64}(Vec(0.0,1.0,0.0)), Node{3,Float64}(Vec(2.0,0.0,0.0)), Node{3,Float64}(Vec(2.0,2.0,0.0))] @@ -384,7 +384,7 @@ end @test a ≈ aa ≈ al ≈ a_rhs1 ≈ a_rhs2 end - # Test nonlinear solution procedure (on linear problem) with affine constraints + # Test nonlinear solution procedure (on linear problem) with affine constraints # using standard assembly (i.e. not local condensation) @testset "nonlinear" begin params = (k=1.0, f=1.0, a=1.0, b=0.2, tol=1e-10, maxiter=2) @@ -392,7 +392,7 @@ end dh = DofHandler(grid); add!(dh, :u, Lagrange{RefLine,1}()); close!(dh) function doassemble!(K, r, dh, a, params) - # Spring elements + # Spring elements k = params.k Ke = [k -k; -k k] # Quick and dirty assem @@ -415,12 +415,12 @@ end a = zeros(ndofs(dh)) # Nonlinear solution - apply!(a, ch) + apply!(a, ch) for niter = 0:params.maxiter doassemble!(K, r, dh, a, params) apply_zero!(K, r, ch) norm(r) < params.tol && break - Δa = -K\r + Δa = -K\r apply_zero!(Δa, ch) a .+= Δa end diff --git a/test/test_deprecations.jl b/test/test_deprecations.jl index 3062e6b747..c7567e811b 100644 --- a/test/test_deprecations.jl +++ b/test/test_deprecations.jl @@ -97,7 +97,7 @@ addfaceset!(grid, "right_face_explicit", Set(Ferrite.FaceIndex(fi[1], fi[2]) for end @testset "vtk_grid" begin - # Ensure no MethodError on pre v1. + # Ensure no MethodError on pre v1. @test_throws ErrorException vtk_grid("old", generate_grid(Line, (1,))) end diff --git a/test/test_dofs.jl b/test/test_dofs.jl index 2f90bd1ef8..348b30e2b9 100644 --- a/test/test_dofs.jl +++ b/test/test_dofs.jl @@ -87,7 +87,7 @@ end @testset "Dofs for quad in 3d (shell)" begin -nodes = [Node{3,Float64}(Vec(0.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,0.0,0.0)), +nodes = [Node{3,Float64}(Vec(0.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,1.0,0.0)), Node{3,Float64}(Vec(0.0,1.0,0.0)), Node{3,Float64}(Vec(2.0,0.0,0.0)), Node{3,Float64}(Vec(2.0,2.0,0.0))] @@ -543,38 +543,38 @@ end # reshape.(Iterators.product(fill([true, false], 9)...) |> collect |> vec .|> collect, Ref((3,3))), [ true true true - true true true - true true true + true true true + true true true ], [ true false false - false true false - false false true + false true false + false false true ], [ true true false - true true true - false true true + true true true + false true true ], # Component coupling [ true true true true - true true true true true true true true - true true true true + true true true true + true true true true ], [ true false false false - false true false false + false true false false false false true false - false false false true + false false false true ], [ true true true false - true true true true true true true true - false true true true + true true true true + false true true true ], ] function is_stored(A, i, j) @@ -589,7 +589,7 @@ end i_dofs = dof_range(sdh, field1_idx) ip1 = sdh.field_interpolations[field1_idx] vdim[1] = typeof(ip1) <: VectorizedInterpolation && size(coupling)[1] == 4 ? Ferrite.get_n_copies(ip1) : 1 - for dim1 in 1:vdim[1] + for dim1 in 1:vdim[1] for cell2_idx in neighbors sdh2 = dh.subdofhandlers[dh.cell_to_subdofhandler[cell2_idx]] coupling_idx[2] = 1 @@ -638,7 +638,7 @@ end close!(dh) for coupling in couplings, cross_coupling in couplings K = create_sparsity_pattern(dh; coupling=coupling, topology = topology, cross_coupling = cross_coupling) - all(coupling) && @test K == create_sparsity_pattern(dh, topology = topology, cross_coupling = cross_coupling) + all(coupling) && @test K == create_sparsity_pattern(dh, topology = topology, cross_coupling = cross_coupling) check_coupling(dh, topology, K, coupling, cross_coupling) end @@ -646,7 +646,7 @@ end @test_throws ErrorException("coupling not square") create_sparsity_pattern(dh; coupling=[true true]) @test_throws ErrorException("coupling not symmetric") create_symmetric_sparsity_pattern(dh; coupling=[true true; false true]) @test_throws ErrorException("could not create coupling") create_symmetric_sparsity_pattern(dh; coupling=falses(100, 100)) - + # Test coupling with subdomains # Note: `check_coupling` works for this case only because the second domain has dofs from the first domain in order. Otherwise tests like in continuous ip are required. grid = generate_grid(Quadrilateral, (2, 1)) @@ -685,8 +685,8 @@ end # Node numbering: # 3 ____ 4 4 - # | | | - # | | | (Beam attached to facet) + # | | | + # | | | (Beam attached to facet) # 1 ____ 2 2 dim = 2 @@ -712,15 +712,15 @@ end celldofs!(dofsbeam, dh, 2) @test dofsbeam == [2, 3, 6] - # Node numbering: + # Node numbering: # 5--------7 # / /| # / / | # 6--------8 | # | | 3 <-- Shell attached on face (4, 3, 7, 8) - # | | / - # | |/ - # 2--------4 + # | | / + # | |/ + # 2--------4 dim = 2 grid = generate_grid(Hexahedron, (1,1,1)) @@ -753,4 +753,4 @@ end @test dofsshell[5:8] == [11,20,15,19] #Shared face dof @test dofsshell[9] == 24 -end \ No newline at end of file +end diff --git a/test/test_facevalues.jl b/test/test_facevalues.jl index 5ee67c3d78..6ed9eea8f3 100644 --- a/test/test_facevalues.jl +++ b/test/test_facevalues.jl @@ -114,7 +114,7 @@ for (scalar_interpol, quad_rule) in ( end end end - # Test that qr, detJdV, normals, and current_facet are copied as expected. + # Test that qr, detJdV, normals, and current_facet are copied as expected. # Note that qr remain aliased, as defined by `copy(qr)=qr`, see quadrature.jl. # Make it easy to test scalar wrapper equality _mock_isequal(a, b) = a == b @@ -132,7 +132,7 @@ for (scalar_interpol, quad_rule) in ( end @testset "show" begin - # Just smoke test to make sure show doesn't error. + # Just smoke test to make sure show doesn't error. fv = FacetValues(FacetQuadratureRule{RefQuadrilateral}(2), Lagrange{RefQuadrilateral,2}()) showstring = sprint(show, MIME"text/plain"(), fv) @test startswith(showstring, "FacetValues(scalar, rdim=2, sdim=2): 2 quadrature points per face") diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl index 3ba3e18487..168e629a2c 100644 --- a/test/test_grid_dofhandler_vtk.jl +++ b/test/test_grid_dofhandler_vtk.jl @@ -94,7 +94,7 @@ end minv, maxv = Ferrite.bounding_box(grid) @test minv ≈ 2left @test maxv ≈ 2right - + end end # of testset @@ -241,7 +241,7 @@ end grid = Grid(cells, nodes) ip1 = DiscontinuousLagrange{RefQuadrilateral, 1}() ip2 = DiscontinuousLagrange{RefTriangle, 1}() - dh = DofHandler(grid); + dh = DofHandler(grid); sdh1 = SubDofHandler(dh, Set([1])); add!(sdh1, :u, ip1); sdh2 = SubDofHandler(dh, Set([2])); add!(sdh2, :u, ip2); close!(dh) @@ -304,7 +304,7 @@ end quadlinefaceskeleton = Ferrite.facetskeleton(quadlinetopo, quadlinegrid) # Test faceskeleton @test Set(linefaceskeleton) == Set(quadlinefaceskeleton) == Set([ - VertexIndex(1,1), VertexIndex(1,2), VertexIndex(2,2),VertexIndex(3,2), + VertexIndex(1,1), VertexIndex(1,2), VertexIndex(2,2),VertexIndex(3,2), ]) # (11) @@ -360,9 +360,9 @@ end @test Set(faceskeleton) == Set(quadfaceskeleton) == Set([ EdgeIndex(1,1), EdgeIndex(1,2), EdgeIndex(1,3), EdgeIndex(1,4), EdgeIndex(2,1), EdgeIndex(2,2), EdgeIndex(2,3), - EdgeIndex(3,2), EdgeIndex(3,3), EdgeIndex(3,4), + EdgeIndex(3,2), EdgeIndex(3,3), EdgeIndex(3,4), EdgeIndex(4,2), EdgeIndex(4,3), - EdgeIndex(5,2), EdgeIndex(5,3), EdgeIndex(5,4), + EdgeIndex(5,2), EdgeIndex(5,3), EdgeIndex(5,4), EdgeIndex(6,2), EdgeIndex(6,3), ]) @@ -412,8 +412,8 @@ end @test Set(faceskeleton) == Set(sfaceskeleton) == Set([ FaceIndex(1,1), FaceIndex(1,2), FaceIndex(1,3), FaceIndex(1,4), FaceIndex(1,5), FaceIndex(1,6), FaceIndex(2,1), FaceIndex(2,2), FaceIndex(2,3), FaceIndex(2,4), FaceIndex(2,6), - FaceIndex(3,1), FaceIndex(3,3), FaceIndex(3,4), FaceIndex(3,5), FaceIndex(3,6), - FaceIndex(4,1), FaceIndex(4,3), FaceIndex(4,4), FaceIndex(4,6), + FaceIndex(3,1), FaceIndex(3,3), FaceIndex(3,4), FaceIndex(3,5), FaceIndex(3,6), + FaceIndex(4,1), FaceIndex(4,3), FaceIndex(4,4), FaceIndex(4,6), ]) # +-----+-----+ @@ -538,8 +538,8 @@ end quadratic_quadgrid = generate_grid(QuadraticQuadrilateral,(3,3)) quadgrid_topology = ExclusiveTopology(quadratic_quadgrid) #quadface_skeleton = Ferrite.faceskeleton(topology, quadgrid) - #@test quadface_skeleton == face_skeleton - + #@test quadface_skeleton == face_skeleton + # add more regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/518 @test all(quadgrid_topology.edge_edge_neighbor .== topology.edge_edge_neighbor) @test all(quadgrid_topology.vertex_vertex_neighbor .== topology.vertex_vertex_neighbor) @@ -665,7 +665,7 @@ end @testset "1d" begin grid = generate_grid(Line, (2,)) - + Ferrite.vertexdof_indices(::VectorLagrangeTest{RefLine,1,2}) = ((1,2),(3,4)) dh1 = DofHandler(grid) add!(dh1, :u, VectorLagrangeTest{RefLine,1,2}()) @@ -710,7 +710,7 @@ end end @testset "VTKFileCollection" begin - @testset "equivalence of addstep methods" begin + @testset "equivalence of addstep methods" begin grid = generate_grid(Triangle, (2,2)) celldata = rand(getncells(grid)) fname = "addstep" @@ -727,7 +727,7 @@ end @test !(isopen(vtk.vtk)) end close.((pvd1, pvd2)) - @test pvd1.step == pvd2.step # Same nr of steps added + @test pvd1.step == pvd2.step # Same nr of steps added for (n, t) in pairs(timesteps) fname1 = string(fname, "_", n, ".vtu") fname2 = string(fname, "2_", n, ".vtu") @@ -738,7 +738,7 @@ end end rm(string(fname, ".pvd")) # Solving https://github.com/Ferrite-FEM/Ferrite.jl/issues/397 - # would allow checking the pvd files as well. + # would allow checking the pvd files as well. end @testset "kwargs forwarding" begin grid = generate_grid(Quadrilateral, (10,10)) diff --git a/test/test_interfacevalues.jl b/test/test_interfacevalues.jl index 2a1955c9d2..a53986b3fe 100644 --- a/test/test_interfacevalues.jl +++ b/test/test_interfacevalues.jl @@ -116,7 +116,7 @@ for (cell_shape, scalar_interpol, quad_rule) in ( #TODO: update interfaces for lines (Line, DiscontinuousLagrange{RefLine, 1}(), FacetQuadratureRule{RefLine}(2)), - (QuadraticLine, DiscontinuousLagrange{RefLine, 2}(), FacetQuadratureRule{RefLine}(2)), + (QuadraticLine, DiscontinuousLagrange{RefLine, 2}(), FacetQuadratureRule{RefLine}(2)), (Quadrilateral, DiscontinuousLagrange{RefQuadrilateral, 1}(), FacetQuadratureRule{RefQuadrilateral}(2)), (QuadraticQuadrilateral, DiscontinuousLagrange{RefQuadrilateral, 2}(), FacetQuadratureRule{RefQuadrilateral}(2)), (Triangle, DiscontinuousLagrange{RefTriangle, 1}(), FacetQuadratureRule{RefTriangle}(2)), @@ -148,7 +148,7 @@ end func_interpol = scalar_interpol for func_interpol in (scalar_interpol, VectorizedInterpolation(scalar_interpol)) - iv = cell_shape ∈ (QuadraticLine, QuadraticQuadrilateral, QuadraticTriangle, QuadraticTetrahedron) ? + iv = cell_shape ∈ (QuadraticLine, QuadraticQuadrilateral, QuadraticTriangle, QuadraticTetrahedron) ? InterfaceValues(quad_rule, func_interpol, ip) : InterfaceValues(quad_rule, func_interpol) test_interfacevalues(grid, iv) end @@ -199,10 +199,10 @@ # end @testset "Unordered nodes 3D" begin dim = 2 - nodes = [Node((-1.0, 0.0, 0.0)), Node((0.0, 0.0, 0.0)), Node((1.0, 0.0, 0.0)), - Node((-1.0, 1.0, 0.0)), Node((0.0, 1.0, 0.0)), Node((1.0, 1.0, 0.0)), - Node((-1.0, 0.0, 1.0)), Node((0.0, 0.0, 1.0)), Node((1.0, 0.0, 1.0)), - Node((-1.0, 1.0, 1.0)), Node((0.0, 1.0, 1.0)), Node((1.0, 1.0, 1.0)), + nodes = [Node((-1.0, 0.0, 0.0)), Node((0.0, 0.0, 0.0)), Node((1.0, 0.0, 0.0)), + Node((-1.0, 1.0, 0.0)), Node((0.0, 1.0, 0.0)), Node((1.0, 1.0, 0.0)), + Node((-1.0, 0.0, 1.0)), Node((0.0, 0.0, 1.0)), Node((1.0, 0.0, 1.0)), + Node((-1.0, 1.0, 1.0)), Node((0.0, 1.0, 1.0)), Node((1.0, 1.0, 1.0)), ] cells = [ Hexahedron((1,2,5,4,7,8,11,10)), diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 367c6268e1..2c883dbb79 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -104,7 +104,7 @@ # The total number of dofs must match the number of base functions totaldofs = sum(length.(Ferrite.vertexdof_indices(interpolation));init=0) totaldofs += sum(length.(Ferrite.facedof_interior_indices(interpolation));init=0) - totaldofs += sum(length.(Ferrite.edgedof_interior_indices(interpolation));init=0) + totaldofs += sum(length.(Ferrite.edgedof_interior_indices(interpolation));init=0) totaldofs += length(Ferrite.volumedof_interior_indices(interpolation)) @test totaldofs == n_basefuncs @@ -164,7 +164,7 @@ for (facenodes, normal) in zip(_bfunc, reference_normals(interpolation)) @test __outward_normal(coords, facenodes) ≈ normal end - end + end # regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/520 #=interpolation_type = typeof(interpolation).name.wrapper diff --git a/test/test_l2_projection.jl b/test/test_l2_projection.jl index 7f046306cb..980cae91c1 100644 --- a/test/test_l2_projection.jl +++ b/test/test_l2_projection.jl @@ -133,7 +133,7 @@ function test_projection_mixedgrid() # use a SymmetricTensor here for testing the symmetric version of project f(x) = SymmetricTensor{2,2,Float64}((1 + x[1]^2, 2x[2]^2, x[1]*x[2])) xe = getcoordinates(mesh, 1) - + # analytical values qp_values = [[f(spatial_coordinate(cv, qp, xe)) for qp in 1:getnquadpoints(cv)]] qp_values_matrix = reduce(hcat, qp_values) diff --git a/test/test_mixeddofhandler.jl b/test/test_mixeddofhandler.jl index eb369444c6..334eed8043 100644 --- a/test/test_mixeddofhandler.jl +++ b/test/test_mixeddofhandler.jl @@ -110,7 +110,7 @@ function test_2d_mixed_1_el() @test ndofs(dh) == 12 @test ndofs_per_cell(dh, 1) == 12 @test celldofs(dh, 1) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] - + @test Set(Ferrite.getfieldnames(dh)) == Set(Ferrite.getfieldnames(dh.subdofhandlers[1])) end @@ -430,12 +430,12 @@ end function test_evaluate_at_grid_nodes() # 5_______6 - # |\ | + # |\ | # | \ | # 3______\4 # | | # | | - # 1_______2 + # 1_______2 nodes = [Node((0.0, 0.0)), Node((1.0, 0.0)), @@ -481,16 +481,16 @@ function test_mixed_grid_show() @test occursin("2 Quadrilateral/Triangle cells", str) end -# regression tests for https://github.com/KristofferC/JuAFEM.jl/issues/315 +# regression tests for https://github.com/KristofferC/JuAFEM.jl/issues/315 function test_subparametric_quad() #linear geometry grid = generate_grid(Quadrilateral, (1,1)) ip = Lagrange{RefQuadrilateral,2}() - + dh = DofHandler(grid) add!(dh, :u, ip^2) close!(dh) - + ch = ConstraintHandler(dh) dbc1 = Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 0.0, 2) add!(ch, dbc1) @@ -505,11 +505,11 @@ function test_subparametric_triangle() grid = generate_grid(Triangle, (1,1)) ip = Lagrange{RefTriangle,2}() - + dh = DofHandler(grid) add!(dh, :u, ip^2) close!(dh) - + ch = ConstraintHandler(dh) dbc1 = Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 0.0, 2) add!(ch, dbc1) @@ -537,12 +537,12 @@ end function test_separate_fields_on_separate_domains() # 5_______6 - # |\ | + # |\ | # | \ | # 3______\4 # | | # | | - # 1_______2 + # 1_______2 # Given: a vector field :q defined on the quad and a scalar field :t defined on the triangles nodes = [Node((0.0, 0.0)), Node((1.0, 0.0)), @@ -640,7 +640,7 @@ function test_vtk_export() end sha = bytes2hex(open(SHA.sha1, filename*".vtu")) @test sha == "339ab8a8a613c2f38af684cccd695ae816671607" - rm(filename*".vtu") # clean up + rm(filename*".vtu") # clean up end @testset "DofHandler" begin diff --git a/test/test_notebooks.jl b/test/test_notebooks.jl index 9982369afc..6ae90aae2a 100644 --- a/test/test_notebooks.jl +++ b/test/test_notebooks.jl @@ -55,5 +55,3 @@ module Cook end end end - - diff --git a/test/test_pointevaluation.jl b/test/test_pointevaluation.jl index 607b5d00b6..7ccf57cf88 100644 --- a/test/test_pointevaluation.jl +++ b/test/test_pointevaluation.jl @@ -137,16 +137,16 @@ function dofhandler2(;three_dimensional=true) ip_f = Lagrange{RefHexahedron,2}() ip_f_v = ip_f^3 qr = QuadratureRule{RefHexahedron}(3) - else + else mesh = generate_grid(Quadrilateral, (20, 20)) f_s = x -> 1.0 + x[1] + x[2] + x[1] * x[2] f_v = x -> Vec{2}((1.0 + x[1] + x[2] + x[1] * x[2], 2.0 - x[1] - x[2] - x[1] * x[2])) points = [Vec((x, x, )) for x in range(0; stop=1, length=100)] ip_f = Lagrange{RefQuadrilateral,2}() ip_f_v = ip_f^2 - qr = QuadratureRule{RefQuadrilateral}(3) + qr = QuadratureRule{RefQuadrilateral}(3) end - + csv = CellValues(qr, ip_f) cvv = CellValues(qr, ip_f_v) dh = DofHandler(mesh); @@ -160,7 +160,7 @@ function dofhandler2(;three_dimensional=true) fe = zeros(ndofs_per_cell(dh)) s_dofs = dof_range(dh, :s) v_dofs = dof_range(dh, :v) - + for cell in CellIterator(dh) fill!(me, 0) fill!(fe, 0) @@ -219,15 +219,15 @@ function dofhandler2(;three_dimensional=true) end function mixed_grid() - ## Mixed grid where not all cells have the same fields + ## Mixed grid where not all cells have the same fields # 5_______6 - # |\ | + # |\ | # | \ | # 3______\4 # | | # | | - # 1_______2 + # 1_______2 nodes = [Node((0.0, 0.0)), Node((1.0, 0.0)), @@ -261,7 +261,7 @@ function mixed_grid() end end - # construct projector + # construct projector projector = L2Projector(ip_quad, mesh; set=getcellset(mesh, "quads")) points = [Vec((x, 2x)) for x in range(0.0; stop=1.0, length=10)] @@ -330,7 +330,7 @@ function first_point_missing() mesh = generate_grid(Quadrilateral, (1, 1)) points = [Vec(2.0, 0.0), Vec(0.0, 0.0)] ph = PointEvalHandler(mesh, points; warn=false) - + @test isnothing(ph.local_coords[1]) @test ph.local_coords[2] == Vec(0.0, 0.0) end diff --git a/test/test_utils.jl b/test/test_utils.jl index 0ba523e99b..c0fb2e9e83 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -12,12 +12,12 @@ reference_face_area(fs::VectorizedInterpolation, f::Int) = reference_face_area(f reference_face_area(fs::Interpolation{Ferrite.RefHypercube{dim}}, face::Int) where {dim} = 2^(dim-1) reference_face_area(fs::Interpolation{RefTriangle}, face::Int) = face == 1 ? sqrt(2) : 1.0 reference_face_area(fs::Interpolation{RefTetrahedron}, face::Int) = face == 3 ? sqrt(2 * 1.5) / 2.0 : 0.5 -function reference_face_area(fs::Interpolation{RefPrism}, face::Int) +function reference_face_area(fs::Interpolation{RefPrism}, face::Int) face == 4 && return √2 face ∈ [1,5] && return 0.5 face ∈ [2,3] && return 1.0 end -function reference_face_area(fs::Interpolation{RefPyramid}, face::Int) +function reference_face_area(fs::Interpolation{RefPyramid}, face::Int) face == 1 && return 1.0 face ∈ [2,3] && return 0.5 face ∈ [4,5] && return sqrt(2)/2 @@ -183,7 +183,7 @@ function calculate_volume(::Lagrange{RefHexahedron, 1}, x::Vector{Vec{3, T}}) wh end function calculate_volume(::Lagrange{RefPrism, order}, x::Vector{Vec{3, T}}) where {T, order} - vol = norm((x[4] - x[1]) ⋅ ((x[2] - x[1]) × (x[3] - x[1]))) / 2.0 + vol = norm((x[4] - x[1]) ⋅ ((x[2] - x[1]) × (x[3] - x[1]))) / 2.0 return vol end @@ -254,7 +254,7 @@ check_equal_or_nan(a::Union{Tensor, Array}, b::Union{Tensor, Array}) = all(check # Hypercube is simply ⨂ᵈⁱᵐ Line :) sample_random_point(::Type{Ferrite.RefHypercube{ref_dim}}) where {ref_dim} = Vec{ref_dim}(2.0 .* rand(Vec{ref_dim}) .- 1.0) # Dirichlet type sampling -function sample_random_point(::Type{Ferrite.RefSimplex{ref_dim}}) where {ref_dim} +function sample_random_point(::Type{Ferrite.RefSimplex{ref_dim}}) where {ref_dim} ξ = rand(ref_dim+1) ξₜ = -log.(ξ) return Vec{ref_dim}(ntuple(i->ξₜ[i], ref_dim) ./ sum(ξₜ)) @@ -290,7 +290,7 @@ module DummyRefShapes struct RefDodecahedron <: Ferrite.AbstractRefShape{3} end function Ferrite.reference_faces(::Type{RefDodecahedron}) return ( - (1, 5, 4, 3, 2), + (1, 5, 4, 3, 2), ) end end diff --git a/test/test_vtk_export.jl b/test/test_vtk_export.jl index 97c88fe20f..f6a3dc467f 100644 --- a/test/test_vtk_export.jl +++ b/test/test_vtk_export.jl @@ -1,6 +1,6 @@ -@testset "VTKFile" begin #TODO: Move all vtk tests here +@testset "VTKFile" begin #TODO: Move all vtk tests here @testset "show(::VTKFile)" begin - mktempdir() do tmp + mktempdir() do tmp grid = generate_grid(Quadrilateral, (2,2)) vtk = VTKFile(joinpath(tmp, "showfile"), grid) showstring_open = sprint(show, MIME"text/plain"(), vtk)