diff --git a/docs/Manifest.toml b/docs/Manifest.toml index 69af303a8e..a6d1f7372b 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -251,6 +251,12 @@ git-tree-sha1 = "05ba0d07cd4fd8b7a39541e31a7b0254704ea581" uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9" version = "0.1.13" +[[deps.CodeTracking]] +deps = ["InteractiveUtils", "REPL", "UUIDs"] +git-tree-sha1 = "cfb7a2e89e245a9d5016b70323db412b3a7438d5" +uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" +version = "3.0.2" + [[deps.CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] git-tree-sha1 = "962834c22b66e32aa10f7611c08c8ca4e20749a9" @@ -315,6 +321,11 @@ weakdeps = ["Dates", "LinearAlgebra"] [deps.Compat.extensions] CompatLinearAlgebraExt = "LinearAlgebra" +[[deps.Compiler]] +git-tree-sha1 = "382d79bfe72a406294faca39ef0c3cef6e6ce1f1" +uuid = "807dbc54-b67e-4c79-8afb-eafe4df6f2e1" +version = "0.1.1" + [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" @@ -1009,6 +1020,12 @@ git-tree-sha1 = "c0c9b76f3520863909825cbecdef58cd63de705a" uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8" version = "3.1.5+0" +[[deps.JuliaInterpreter]] +deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"] +git-tree-sha1 = "58927c485919bf17ea308d9d82156de1adf4b006" +uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a" +version = "0.10.12" + [[deps.JuliaSyntaxHighlighting]] deps = ["StyledStrings"] uuid = "ac6e5ff7-fb65-4e79-a425-ec3bc9c03011" @@ -1274,6 +1291,12 @@ git-tree-sha1 = "f00544d95982ea270145636c181ceda21c4e2575" uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" version = "1.2.0" +[[deps.LoweredCodeUtils]] +deps = ["CodeTracking", "Compiler", "JuliaInterpreter"] +git-tree-sha1 = "5d4278f755440f70648d80cc6225f51e78e94094" +uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b" +version = "3.5.1" + [[deps.METIS_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] git-tree-sha1 = "2eefa8baa858871ae7770c98c3c2a7e46daba5b4" @@ -1815,6 +1838,16 @@ git-tree-sha1 = "62389eeff14780bfe55195b7204c0d8738436d64" uuid = "ae029012-a4dd-5104-9daa-d747884805df" version = "1.3.1" +[[deps.Revise]] +deps = ["CodeTracking", "FileWatching", "InteractiveUtils", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Preferences", "REPL", "UUIDs"] +git-tree-sha1 = "5f4f629c085b87e71125eec6773f5f872c74a47a" +uuid = "295af30f-e4ad-537b-8983-00126c2a3abe" +version = "3.14.2" +weakdeps = ["Distributed"] + + [deps.Revise.extensions] + DistributedExt = "Distributed" + [[deps.RuntimeGeneratedFunctions]] deps = ["ExprTools", "SHA", "Serialization"] git-tree-sha1 = "cfcdc949c4660544ab0fdeed169561cb22f835f4" diff --git a/docs/Project.toml b/docs/Project.toml index 3a40d906bf..5477d9bd47 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -22,6 +22,7 @@ Optim = "429524aa-4258-5aef-a3af-852621145aeb" OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TaskLocalValues = "ed4db957-447d-4319-bfb6-7fa9ae7ecf34" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" diff --git a/docs/src/literate-tutorials/computational_homogenization.jl b/docs/src/literate-tutorials/computational_homogenization.jl index 5165d5b3e1..6039070fc0 100644 --- a/docs/src/literate-tutorials/computational_homogenization.jl +++ b/docs/src/literate-tutorials/computational_homogenization.jl @@ -323,14 +323,14 @@ function doassemble!(cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler, for q_point in 1:getnquadpoints(cellvalues) dΩ = getdetJdV(cellvalues, q_point) for i in 1:n_basefuncs - δεi = shape_symmetric_gradient(cellvalues, q_point, i) + δNεᵢ = shape_symmetric_gradient(cellvalues, q_point, i) for j in 1:n_basefuncs - δεj = shape_symmetric_gradient(cellvalues, q_point, j) - Ke[i, j] += (δεi ⊡ E ⊡ δεj) * dΩ + Nεⱼ = shape_symmetric_gradient(cellvalues, q_point, j) + Ke[i, j] += (δNεᵢ ⊡ E ⊡ Nεⱼ) * dΩ end for (rhs, ε) in enumerate(εᴹ) σᴹ = E ⊡ ε - fe[i, rhs] += (- δεi ⊡ σᴹ) * dΩ + fe[i, rhs] += (- δNεᵢ ⊡ σᴹ) * dΩ end end end diff --git a/docs/src/literate-tutorials/dg_heat_equation.jl b/docs/src/literate-tutorials/dg_heat_equation.jl index 4b5ef9fae5..c3b9a3ec89 100644 --- a/docs/src/literate-tutorials/dg_heat_equation.jl +++ b/docs/src/literate-tutorials/dg_heat_equation.jl @@ -213,15 +213,15 @@ function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues) dΩ = getdetJdV(cellvalues, q_point) ## Loop over test shape functions for i in 1:n_basefuncs - δu = shape_value(cellvalues, q_point, i) - ∇δu = shape_gradient(cellvalues, q_point, i) + δNᵢ = shape_value(cellvalues, q_point, i) + ∇δNᵢ = shape_gradient(cellvalues, q_point, i) ## Add contribution to fe - fe[i] += δu * dΩ + fe[i] += δNᵢ * dΩ ## Loop over trial shape functions for j in 1:n_basefuncs - ∇u = shape_gradient(cellvalues, q_point, j) + ∇Nⱼ = shape_gradient(cellvalues, q_point, j) ## Add contribution to Ke - Ke[i, j] += (∇δu ⋅ ∇u) * dΩ + Ke[i, j] += (∇δNᵢ ⋅ ∇Nⱼ) * dΩ end end end @@ -240,15 +240,15 @@ function assemble_interface!(Ki::Matrix, iv::InterfaceValues, μ::Float64) ## Loop over test shape functions for i in 1:getnbasefunctions(iv) ## Multiply the jump by the negative normal to get the definition from the theory section. - δu_jump = shape_value_jump(iv, q_point, i) * (-normal) - ∇δu_avg = shape_gradient_average(iv, q_point, i) + δNᵢ_jump = shape_value_jump(iv, q_point, i) * (-normal) + ∇δNᵢ_avg = shape_gradient_average(iv, q_point, i) ## Loop over trial shape functions for j in 1:getnbasefunctions(iv) ## Multiply the jump by the negative normal to get the definition from the theory section. - u_jump = shape_value_jump(iv, q_point, j) * (-normal) - ∇u_avg = shape_gradient_average(iv, q_point, j) + Nⱼ = shape_value_jump(iv, q_point, j) * (-normal) + ∇Nⱼ_avg = shape_gradient_average(iv, q_point, j) ## Add contribution to Ki - Ki[i, j] += -(δu_jump ⋅ ∇u_avg + ∇δu_avg ⋅ u_jump) * dΓ + μ * (δu_jump ⋅ u_jump) * dΓ + Ki[i, j] += -(δNᵢ_jump ⋅ ∇Nⱼ_avg + ∇δNᵢ_avg ⋅ Nⱼ) * dΓ + μ * (δNᵢ_jump ⋅ Nⱼ) * dΓ end end end @@ -266,9 +266,9 @@ function assemble_boundary!(fe::Vector, fv::FacetValues) ∂Ω = getdetJdV(fv, q_point) ## Loop over test shape functions for i in 1:getnbasefunctions(fv) - δu = shape_value(fv, q_point, i) + δNᵢ = shape_value(fv, q_point, i) boundary_flux = normal[2] - fe[i] = boundary_flux * δu * ∂Ω + fe[i] = boundary_flux * δNᵢ * ∂Ω end end return fe diff --git a/docs/src/literate-tutorials/heat_equation.jl b/docs/src/literate-tutorials/heat_equation.jl index d132f4ea72..9c77095e6b 100644 --- a/docs/src/literate-tutorials/heat_equation.jl +++ b/docs/src/literate-tutorials/heat_equation.jl @@ -132,11 +132,9 @@ close!(ch) # # !!! note "Notation" # Comparing with the brief finite element introduction in [Introduction to FEM](@ref), -# the variables `δu`, `∇δu` and `∇u` are actually $\phi_i(\textbf{x}_q)$, $\nabla +# the variables `δNᵢ`, `∇δNᵢ` and `∇Nⱼ` are actually $\phi_i(\textbf{x}_q)$, $\nabla # \phi_i(\textbf{x}_q)$ and $\nabla \phi_j(\textbf{x}_q)$, i.e. the evaluation of the -# trial and test functions in the quadrature point ``\textbf{x}_q``. However, to -# underline the strong parallel between the weak form and the implementation, this -# example uses the symbols appearing in the weak form. +# trial and test functions in the quadrature point ``\textbf{x}_q``. function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues) n_basefuncs = getnbasefunctions(cellvalues) @@ -149,15 +147,15 @@ function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues) dΩ = getdetJdV(cellvalues, q_point) ## Loop over test shape functions for i in 1:n_basefuncs - δu = shape_value(cellvalues, q_point, i) - ∇δu = shape_gradient(cellvalues, q_point, i) + δNᵢ = shape_value(cellvalues, q_point, i) + ∇δNᵢ = shape_gradient(cellvalues, q_point, i) ## Add contribution to fe - fe[i] += δu * dΩ + fe[i] += δNᵢ * dΩ ## Loop over trial shape functions for j in 1:n_basefuncs - ∇u = shape_gradient(cellvalues, q_point, j) + ∇Nⱼ = shape_gradient(cellvalues, q_point, j) ## Add contribution to Ke - Ke[i, j] += (∇δu ⋅ ∇u) * dΩ + Ke[i, j] += (∇δNᵢ ⋅ ∇Nⱼ) * dΩ end end end diff --git a/docs/src/literate-tutorials/hyperelasticity.jl b/docs/src/literate-tutorials/hyperelasticity.jl index 6ab05b3138..55caa02e92 100644 --- a/docs/src/literate-tutorials/hyperelasticity.jl +++ b/docs/src/literate-tutorials/hyperelasticity.jl @@ -262,16 +262,16 @@ function assemble_element!(ke, ge, cell, cv, fv, mp, ue, ΓN) ## Loop over test functions for i in 1:ndofs ## Test function and gradient - δui = shape_value(cv, qp, i) - ∇δui = shape_gradient(cv, qp, i) + δNᵢ = shape_value(cv, qp, i) + ∇δNᵢ = shape_gradient(cv, qp, i) ## Add contribution to the residual from this test function - ge[i] += (∇δui ⊡ P - δui ⋅ b) * dΩ + ge[i] += (∇δNᵢ ⊡ P - δNᵢ ⋅ b) * dΩ - ∇δui∂P∂F = ∇δui ⊡ ∂P∂F # Hoisted computation + ∇δNᵢ∂P∂F = ∇δNᵢ ⊡ ∂P∂F # Hoisted computation for j in 1:ndofs - ∇δuj = shape_gradient(cv, qp, j) + ∇δNⱼ = shape_gradient(cv, qp, j) ## Add contribution to the tangent - ke[i, j] += (∇δui∂P∂F ⊡ ∇δuj) * dΩ + ke[i, j] += (∇δNᵢ∂P∂F ⊡ ∇δNⱼ) * dΩ end end end @@ -284,8 +284,8 @@ function assemble_element!(ke, ge, cell, cv, fv, mp, ue, ΓN) t = tn * getnormal(fv, q_point) dΓ = getdetJdV(fv, q_point) for i in 1:ndofs - δui = shape_value(fv, q_point, i) - ge[i] -= (δui ⋅ t) * dΓ + δNᵢ = shape_value(fv, q_point, i) + ge[i] -= (δNᵢ ⋅ t) * dΓ end end end diff --git a/docs/src/literate-tutorials/incompressible_elasticity.jl b/docs/src/literate-tutorials/incompressible_elasticity.jl index 58b0e8b178..31eff9e140 100644 --- a/docs/src/literate-tutorials/incompressible_elasticity.jl +++ b/docs/src/literate-tutorials/incompressible_elasticity.jl @@ -137,14 +137,14 @@ function assemble_up!(Ke, fe, cell, cellvalues, facetvalues_u, grid, mp, t, dofr end for (iₚ, Iₚ) in pairs(dofrange_p) - δp = shape_value(cellvalues.p, q_point, iₚ) + δNpᵢ = shape_value(cellvalues.p, q_point, iₚ) for (jᵤ, Jᵤ) in pairs(dofrange_u) - divδu = shape_divergence(cellvalues.u, q_point, jᵤ) - Ke[Iₚ, Jᵤ] += -δp * divδu * dΩ + divδNuⱼ = shape_divergence(cellvalues.u, q_point, jᵤ) + Ke[Iₚ, Jᵤ] += -δNpᵢ * divδNuⱼ * dΩ end for (jₚ, Jₚ) in pairs(dofrange_p[1:iₚ]) - p = shape_value(cellvalues.p, q_point, jₚ) - Ke[Iₚ, Jₚ] += - 1 / mp.K * δp * p * dΩ + Npⱼ = shape_value(cellvalues.p, q_point, jₚ) + Ke[Iₚ, Jₚ] += - 1 / mp.K * δNpᵢ * Npⱼ * dΩ end end @@ -161,8 +161,8 @@ function assemble_up!(Ke, fe, cell, cellvalues, facetvalues_u, grid, mp, t, dofr for q_point in 1:getnquadpoints(facetvalues_u) dΓ = getdetJdV(facetvalues_u, q_point) for (iᵤ, Iᵤ) in pairs(dofrange_u) - δu = shape_value(facetvalues_u, q_point, iᵤ) - fe[Iᵤ] += (δu ⋅ t) * dΓ + δNuᵢ = shape_value(facetvalues_u, q_point, iᵤ) + fe[Iᵤ] += (δNuᵢ ⋅ t) * dΓ end end end diff --git a/docs/src/literate-tutorials/linear_elasticity.jl b/docs/src/literate-tutorials/linear_elasticity.jl index a5574f9bd4..7d9af1ed78 100644 --- a/docs/src/literate-tutorials/linear_elasticity.jl +++ b/docs/src/literate-tutorials/linear_elasticity.jl @@ -192,8 +192,8 @@ function assemble_external_forces!(f_ext, dh, facetset, facetvalues, prescribed_ ## Get the integration weight for the current quadrature point. dΓ = getdetJdV(facetvalues, qp) for i in 1:getnbasefunctions(facetvalues) - Nᵢ = shape_value(facetvalues, qp, i) - fe_ext[i] += tₚ ⋅ Nᵢ * dΓ + δNᵢ = shape_value(facetvalues, qp, i) + fe_ext[i] += tₚ ⋅ δNᵢ * dΓ end end ## Add the local contributions to the correct indices in the global external force vector @@ -267,11 +267,11 @@ function assemble_cell!(ke, cellvalues, C) dΩ = getdetJdV(cellvalues, q_point) for i in 1:getnbasefunctions(cellvalues) ## Gradient of the test function - ∇Nᵢ = shape_gradient(cellvalues, q_point, i) + ∇δNᵢ = shape_gradient(cellvalues, q_point, i) for j in 1:getnbasefunctions(cellvalues) ## Symmetric gradient of the trial function ∇ˢʸᵐNⱼ = shape_symmetric_gradient(cellvalues, q_point, j) - ke[i, j] += (∇Nᵢ ⊡ C ⊡ ∇ˢʸᵐNⱼ) * dΩ + ke[i, j] += (∇δNᵢ ⊡ C ⊡ ∇ˢʸᵐNⱼ) * dΩ end end end diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index 72dc235aae..77b9df8674 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -265,10 +265,10 @@ function assemble_mass_matrix(cv::MultiFieldCellValues, M::SparseMatrixCSC, dh:: ## Remember that we assemble a vector mass term, hence the dot product. ## There is only one time derivative on the left hand side, so only one mass block is non-zero. for i in 1:n_basefuncs_v - φᵢ = shape_value(cv.v, q_point, i) + Nvᵢ = shape_value(cv.v, q_point, i) for j in 1:n_basefuncs_v - φⱼ = shape_value(cv.v, q_point, j) - Mₑ[BlockIndex((v▄, v▄), (i, j))] += φᵢ ⋅ φⱼ * dΩ + δNvⱼ = shape_value(cv.v, q_point, j) + Mₑ[BlockIndex((v▄, v▄), (i, j))] += Nvᵢ ⋅ δNvⱼ * dΩ end end end @@ -309,20 +309,20 @@ function assemble_stokes_matrix(cv::MultiFieldCellValues, ν, K::SparseMatrixCSC # Assemble local viscosity block of $A$ #+ for i in 1:n_basefuncs_v - ∇φᵢ = shape_gradient(cv.v, q_point, i) + ∇Nvᵢ = shape_gradient(cv.v, q_point, i) for j in 1:n_basefuncs_v - ∇φⱼ = shape_gradient(cv.v, q_point, j) - Kₑ[BlockIndex((v▄, v▄), (i, j))] -= ν * ∇φᵢ ⊡ ∇φⱼ * dΩ + ∇δNvⱼ = shape_gradient(cv.v, q_point, j) + Kₑ[BlockIndex((v▄, v▄), (i, j))] -= ν * ∇Nvᵢ ⊡ ∇δNvⱼ * dΩ end end # Assemble local pressure and incompressibility blocks of $B^{\textrm{T}}$ and $B$. #+ for j in 1:n_basefuncs_p - ψ = shape_value(cv.p, q_point, j) + δNpⱼ = shape_value(cv.p, q_point, j) for i in 1:n_basefuncs_v - divφ = shape_divergence(cv.v, q_point, i) - Kₑ[BlockIndex((v▄, p▄), (i, j))] += (divφ * ψ) * dΩ - Kₑ[BlockIndex((p▄, v▄), (j, i))] += (ψ * divφ) * dΩ + div_δNvᵢ = shape_divergence(cv.v, q_point, i) + Kₑ[BlockIndex((v▄, p▄), (i, j))] += (div_δNvᵢ * δNpⱼ) * dΩ + Kₑ[BlockIndex((p▄, v▄), (j, i))] += (δNpⱼ * div_δNvᵢ) * dΩ end end end @@ -411,7 +411,7 @@ function navierstokes_rhs_element!(dvₑ, vₑ, cv) ∇v = function_gradient(cv.v, q_point, vₑ) v = function_value(cv.v, q_point, vₑ) for j in 1:n_basefuncs - φⱼ = shape_value(cv.v, q_point, j) + δNvⱼ = shape_value(cv.v, q_point, j) # Note that in Tensors.jl the definition $\textrm{grad} v = \nabla v$ holds. # With this information it can be quickly shown in index notation that # ```math @@ -419,7 +419,7 @@ function navierstokes_rhs_element!(dvₑ, vₑ, cv) # ``` # where we should pay attentation to the transpose of the gradient. #+ - dvₑ[j] -= v ⋅ ∇v' ⋅ φⱼ * dΩ + dvₑ[j] -= v ⋅ ∇v' ⋅ δNvⱼ * dΩ end end return @@ -468,7 +468,7 @@ function navierstokes_jac_element!(Jₑ, vₑ, cv) ∇v = function_gradient(cv.v, q_point, vₑ) v = function_value(cv.v, q_point, vₑ) for j in 1:n_basefuncs - φⱼ = shape_value(cv.v, q_point, j) + δNvⱼ = shape_value(cv.v, q_point, j) # Note that in Tensors.jl the definition $\textrm{grad} v = \nabla v$ holds. # With this information it can be quickly shown in index notation that # ```math @@ -477,9 +477,9 @@ function navierstokes_jac_element!(Jₑ, vₑ, cv) # where we should pay attentation to the transpose of the gradient. #+ for i in 1:n_basefuncs - φᵢ = shape_value(cv.v, q_point, i) - ∇φᵢ = shape_gradient(cv.v, q_point, i) - Jₑ[j, i] -= (φᵢ ⋅ ∇v' + v ⋅ ∇φᵢ') ⋅ φⱼ * dΩ + Nvᵢ = shape_value(cv.v, q_point, i) + ∇Nvᵢ = shape_gradient(cv.v, q_point, i) + Jₑ[j, i] -= (Nvᵢ ⋅ ∇v' + v ⋅ ∇Nvᵢ') ⋅ δNvⱼ * dΩ end end end diff --git a/docs/src/literate-tutorials/plasticity.jl b/docs/src/literate-tutorials/plasticity.jl index c99a199cc5..092e66a298 100644 --- a/docs/src/literate-tutorials/plasticity.jl +++ b/docs/src/literate-tutorials/plasticity.jl @@ -211,11 +211,11 @@ function assemble_cell!(Ke, re, cell, cellvalues, material, ue, state, state_old dΩ = getdetJdV(cellvalues, q_point) for i in 1:n_basefuncs - δϵ = shape_symmetric_gradient(cellvalues, q_point, i) - re[i] += (δϵ ⊡ σ) * dΩ # add internal force to residual + δN_ϵ = shape_symmetric_gradient(cellvalues, q_point, i) + re[i] += (δN_ϵ ⊡ σ) * dΩ # add internal force to residual for j in 1:i # loop only over lower half - Δϵ = shape_symmetric_gradient(cellvalues, q_point, j) - Ke[i, j] += δϵ ⊡ D ⊡ Δϵ * dΩ + ΔN_ϵ = shape_symmetric_gradient(cellvalues, q_point, j) + Ke[i, j] += δN_ϵ ⊡ D ⊡ ΔN_ϵ * dΩ end end end @@ -243,8 +243,8 @@ function doassemble_neumann!(r, dh, facetset, facetvalues, t) for q_point in 1:getnquadpoints(facetvalues) dΓ = getdetJdV(facetvalues, q_point) for i in 1:n_basefuncs - δu = shape_value(facetvalues, q_point, i) - re[i] -= (δu ⋅ t) * dΓ + δNuᵢ = shape_value(facetvalues, q_point, i) + re[i] -= (δNuᵢ ⋅ t) * dΓ end end assemble!(r, celldofs(fc), re) diff --git a/docs/src/literate-tutorials/porous_media.jl b/docs/src/literate-tutorials/porous_media.jl index 73219232d2..c5704058c0 100644 --- a/docs/src/literate-tutorials/porous_media.jl +++ b/docs/src/literate-tutorials/porous_media.jl @@ -122,11 +122,11 @@ function element_routine!(Ke, re, material::Elastic, cv::CellValues, a, args...) ϵ = function_symmetric_gradient(cv, q_point, a) σ = material.C ⊡ ϵ for i in 1:n_basefuncs - δ∇N = shape_symmetric_gradient(cv, q_point, i) - re[i] += (δ∇N ⊡ σ) * dΩ + ∇δNᵢ = shape_symmetric_gradient(cv, q_point, i) + re[i] += (∇δNᵢ ⊡ σ) * dΩ for j in 1:n_basefuncs - ∇N = shape_symmetric_gradient(cv, q_point, j) - Ke[i, j] += (δ∇N ⊡ material.C ⊡ ∇N) * dΩ + ∇Nⱼ = shape_symmetric_gradient(cv, q_point, j) + Ke[i, j] += (∇δNᵢ ⊡ material.C ⊡ ∇Nⱼ) * dΩ end end end @@ -167,31 +167,31 @@ function element_routine!(Ke, re, m::PoroElastic, cv::MultiFieldCellValues, a, a σ_eff = C ⊡ ϵ ## Variation of u_i for (iᵤ, Iᵤ) in pairs(dr_u) - ∇δNu = shape_symmetric_gradient(cv.u, q_point, iᵤ) - div_δNu = shape_divergence(cv.u, q_point, iᵤ) - re[Iᵤ] += (∇δNu ⊡ σ_eff - div_δNu * p * m.α) * dΩ + ∇δNuᵢ = shape_symmetric_gradient(cv.u, q_point, iᵤ) + div_δNuᵢ = shape_divergence(cv.u, q_point, iᵤ) + re[Iᵤ] += (∇δNuᵢ ⊡ σ_eff - div_δNuᵢ * p * m.α) * dΩ for (jᵤ, Jᵤ) in pairs(dr_u) - ∇Nu = shape_symmetric_gradient(cv.u, q_point, jᵤ) - Ke[Iᵤ, Jᵤ] += (∇δNu ⊡ C ⊡ ∇Nu) * dΩ + ∇Nuⱼ = shape_symmetric_gradient(cv.u, q_point, jᵤ) + Ke[Iᵤ, Jᵤ] += (∇δNuᵢ ⊡ C ⊡ ∇Nuⱼ) * dΩ end for (jₚ, Jₚ) in pairs(dr_p) - Np = shape_value(cv.p, q_point, jₚ) - Ke[Iᵤ, Jₚ] -= (div_δNu * m.α * Np) * dΩ + Npⱼ = shape_value(cv.p, q_point, jₚ) + Ke[Iᵤ, Jₚ] -= (div_δNuᵢ * m.α * Npⱼ) * dΩ end end ## Variation of p_i for (iₚ, Iₚ) in pairs(dr_p) - δNp = shape_value(cv.p, q_point, iₚ) - ∇δNp = shape_gradient(cv.p, q_point, iₚ) - re[Iₚ] += (δNp * (m.α * tr_ϵ_dot + m.β * pdot) + m.k * (∇δNp ⋅ ∇p)) * dΩ + δNpᵢ = shape_value(cv.p, q_point, iₚ) + ∇δNpᵢ = shape_gradient(cv.p, q_point, iₚ) + re[Iₚ] += (δNpᵢ * (m.α * tr_ϵ_dot + m.β * pdot) + m.k * (∇δNpᵢ ⋅ ∇p)) * dΩ for (jᵤ, Jᵤ) in pairs(dr_u) - div_Nu = shape_divergence(cv.u, q_point, jᵤ) - Ke[Iₚ, Jᵤ] += δNp * (m.α / Δt) * div_Nu * dΩ + div_Nuⱼ = shape_divergence(cv.u, q_point, jᵤ) + Ke[Iₚ, Jᵤ] += δNpᵢ * (m.α / Δt) * div_Nuⱼ * dΩ end for (jₚ, Jₚ) in pairs(dr_p) - ∇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Ω + ∇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 diff --git a/docs/src/literate-tutorials/reactive_surface.jl b/docs/src/literate-tutorials/reactive_surface.jl index e02c219753..e2fe880dc2 100644 --- a/docs/src/literate-tutorials/reactive_surface.jl +++ b/docs/src/literate-tutorials/reactive_surface.jl @@ -105,13 +105,13 @@ function assemble_element_mass!(Me::Matrix, cellvalues::CellValues) dΩ = getdetJdV(cellvalues, q_point) ## Loop over test shape functions for i in 1:n_basefuncs - δuᵢ = shape_value(cellvalues, q_point, i) + δNᵢ = shape_value(cellvalues, q_point, i) ## Loop over trial shape functions for j in 1:n_basefuncs - δuⱼ = shape_value(cellvalues, q_point, j) + Nⱼ = shape_value(cellvalues, q_point, j) ## Add contribution to Ke - Me₁[i, j] += (δuᵢ * δuⱼ) * dΩ - Me₂[i, j] += (δuᵢ * δuⱼ) * dΩ + Me₁[i, j] += (δNᵢ * Nⱼ) * dΩ + Me₂[i, j] += (δNᵢ * Nⱼ) * dΩ end end end @@ -136,13 +136,13 @@ function assemble_element_diffusion!(De::Matrix, cellvalues::CellValues, materia dΩ = getdetJdV(cellvalues, q_point) ## Loop over test shape functions for i in 1:n_basefuncs - ∇δuᵢ = shape_gradient(cellvalues, q_point, i) + ∇δNᵢ = shape_gradient(cellvalues, q_point, i) ## Loop over trial shape functions for j in 1:n_basefuncs - ∇δuⱼ = shape_gradient(cellvalues, q_point, j) + ∇Nⱼ = shape_gradient(cellvalues, q_point, j) ## Add contribution to Ke - De₁[i, j] += D₁ * (∇δuᵢ ⋅ ∇δuⱼ) * dΩ - De₂[i, j] += D₂ * (∇δuᵢ ⋅ ∇δuⱼ) * dΩ + De₁[i, j] += D₁ * (∇δNᵢ ⋅ ∇Nⱼ) * dΩ + De₂[i, j] += D₂ * (∇δNᵢ ⋅ ∇Nⱼ) * dΩ end end end diff --git a/docs/src/literate-tutorials/stokes-flow.jl b/docs/src/literate-tutorials/stokes-flow.jl index 422cffea48..61c02b0c9b 100644 --- a/docs/src/literate-tutorials/stokes-flow.jl +++ b/docs/src/literate-tutorials/stokes-flow.jl @@ -389,10 +389,10 @@ function assemble_system!(K, f, dh, cvu, cvp) ndofs_u = length(range_u) range_p = dof_range(dh, :p) ndofs_p = length(range_p) - ϕᵤ = Vector{Vec{2, Float64}}(undef, ndofs_u) - ∇ϕᵤ = Vector{Tensor{2, 2, Float64, 4}}(undef, ndofs_u) - divϕᵤ = Vector{Float64}(undef, ndofs_u) - ϕₚ = Vector{Float64}(undef, ndofs_p) + Nᵤ = Vector{Vec{2, Float64}}(undef, ndofs_u) + ∇Nᵤ = Vector{Tensor{2, 2, Float64, 4}}(undef, ndofs_u) + divNᵤ = Vector{Float64}(undef, ndofs_u) + Nₚ = Vector{Float64}(undef, ndofs_p) for cell in CellIterator(dh) reinit!(cvu, cell) reinit!(cvp, cell) @@ -401,31 +401,31 @@ function assemble_system!(K, f, dh, cvu, cvp) for qp in 1:getnquadpoints(cvu) dΩ = getdetJdV(cvu, qp) for i in 1:ndofs_u - ϕᵤ[i] = shape_value(cvu, qp, i) - ∇ϕᵤ[i] = shape_gradient(cvu, qp, i) - divϕᵤ[i] = shape_divergence(cvu, qp, i) + Nᵤ[i] = shape_value(cvu, qp, i) + ∇Nᵤ[i] = shape_gradient(cvu, qp, i) + divNᵤ[i] = shape_divergence(cvu, qp, i) end for i in 1:ndofs_p - ϕₚ[i] = shape_value(cvp, qp, i) + Nₚ[i] = shape_value(cvp, qp, i) end ## u-u for (i, I) in pairs(range_u), (j, J) in pairs(range_u) - ke[I, J] += (∇ϕᵤ[i] ⊡ ∇ϕᵤ[j]) * dΩ + ke[I, J] += (∇Nᵤ[i] ⊡ ∇Nᵤ[j]) * dΩ end ## u-p for (i, I) in pairs(range_u), (j, J) in pairs(range_p) - ke[I, J] += (-divϕᵤ[i] * ϕₚ[j]) * dΩ + ke[I, J] += (-divNᵤ[i] * Nₚ[j]) * dΩ end ## p-u for (i, I) in pairs(range_p), (j, J) in pairs(range_u) - ke[I, J] += (-divϕᵤ[j] * ϕₚ[i]) * dΩ + ke[I, J] += (-divNᵤ[j] * Nₚ[i]) * dΩ end ## rhs for (i, I) in pairs(range_u) x = spatial_coordinate(cvu, qp, getcoordinates(cell)) b = exp(-100 * norm(x - Vec{2}((0.75, 0.1)))^2) bv = Vec{2}((b, 0.0)) - fe[I] += (ϕᵤ[i] ⋅ bv) * dΩ + fe[I] += (Nᵤ[i] ⋅ bv) * dΩ end end assemble!(assembler, celldofs(cell), ke, fe) diff --git a/docs/src/literate-tutorials/transient_heat_equation.jl b/docs/src/literate-tutorials/transient_heat_equation.jl index 9f5a8a8c7a..646884e232 100644 --- a/docs/src/literate-tutorials/transient_heat_equation.jl +++ b/docs/src/literate-tutorials/transient_heat_equation.jl @@ -35,27 +35,34 @@ # ``` # The semidiscrete weak form is given by # ```math -# \int_{\Omega}v \frac{\partial u}{\partial t} \ \mathrm{d}\Omega + \int_{\Omega} k \nabla v \cdot \nabla u \ \mathrm{d}\Omega = \int_{\Omega} f v \ \mathrm{d}\Omega, +# \int_{\Omega}\delta u \frac{\partial u}{\partial t} \ \mathrm{d}\Omega + \int_{\Omega} k \nabla \delta u \cdot \nabla u \ \mathrm{d}\Omega = \int_{\Omega} f \delta u \ \mathrm{d}\Omega, # ``` -# where $v$ is a suitable test function. Now, we still need to discretize the time derivative. An implicit Euler scheme is applied, +# where $\delta u$ is a suitable test function. Now, we still need to discretize the time derivative. An implicit Euler scheme is applied, # which yields: # ```math -# \int_{\Omega} v\, u_{n+1}\ \mathrm{d}\Omega + \Delta t\int_{\Omega} k \nabla v \cdot \nabla u_{n+1} \ \mathrm{d}\Omega = \Delta t\int_{\Omega} f v \ \mathrm{d}\Omega + \int_{\Omega} v \, u_{n} \ \mathrm{d}\Omega. +# \int_{\Omega} \delta u\, u_{n+1}\ \mathrm{d}\Omega + \Delta t\int_{\Omega} k \nabla \delta u \cdot \nabla u_{n+1} \ \mathrm{d}\Omega = \Delta t\int_{\Omega} f \delta u \ \mathrm{d}\Omega + \int_{\Omega} \delta u \, u_{n} \ \mathrm{d}\Omega. # ``` # If we assemble the discrete operators, we get the following algebraic system: # ```math # \mathbf{M} \mathbf{u}_{n+1} + Δt \mathbf{K} \mathbf{u}_{n+1} = Δt \mathbf{f} + \mathbf{M} \mathbf{u}_{n} # ``` +# Here, `M` refers to the so called mass matrix, which always occurs in time related terms, i.e. +# ```math +# M_{ij} = \int_{\Omega} \delta N_i \, N_j \ \mathrm{d}\Omega, +# ``` +# where $N_j$ and $\delta N_i$ are the shape functions of the trial and test functions, respectively. # In this example we apply the boundary conditions to the assembled discrete operators (mass matrix $\mathbf{M}$ and stiffnes matrix $\mathbf{K}$) # only once. We utilize the fact that in finite element computations Dirichlet conditions can be applied by # zero out rows and columns that correspond -# to a prescribed dof in the system matrix ($\mathbf{A} = Δt \mathbf{K} + \mathbf{M}$) and setting the value of the right-hand side vector to the value -# of the Dirichlet condition. Thus, we only need to apply in every time step the Dirichlet condition to the right-hand side of the problem. +# to a prescribed dof in the system matrix ($\mathbf{A} = Δt \mathbf{K} + \mathbf{M}$) and setting the value of the +# right-hand side vector to the value of the Dirichlet condition. Thus, we only need to apply in every time step the +# Dirichlet condition to the right-hand side of the problem. For more details +# on the derivation and discretisation, see [Introduction to FEM](@ref). #- # ## Commented program # # Now we solve the problem in Ferrite. What follows is a program spliced with comments. -#md # The full program, without comments, can be found in the next [section](@ref heat_equation-plain-program). +#md # The full program, without comments, can be found in the next [section](@ref transient_heat_equation-plain-program). # # First we load Ferrite, and some other packages we need. using Ferrite, SparseArrays, WriteVTK @@ -75,11 +82,6 @@ add!(dh, :u, ip) close!(dh); # By means of the `DofHandler` we can allocate the needed `SparseMatrixCSC`. -# `M` refers here to the so called mass matrix, which always occurs in time related terms, i.e. -# ```math -# M_{ij} = \int_{\Omega} v_i \, u_j \ \mathrm{d}\Omega, -# ``` -# where $u_i$ and $v_j$ are trial and test functions, respectively. K = allocate_matrix(dh); M = allocate_matrix(dh); # We also preallocate the right hand side @@ -107,76 +109,48 @@ close!(ch) update!(ch, 0.0); # ### Assembling the linear system -# As in the heat equation example we define a `doassemble!` function that assembles the diffusion parts of the equation: -function doassemble_K!(K::SparseMatrixCSC, f::Vector, cellvalues::CellValues, dh::DofHandler) - +# As in [Tutorial 1: Heat equation](heat_equation.md), we define a `doassemble!` function that assembles the +# diffusion and diffusive parts of the equation: +function doassemble!(K::SparseMatrixCSC, M::SparseMatrixCSC, f::Vector, cellvalues::CellValues, dh::DofHandler) n_basefuncs = getnbasefunctions(cellvalues) Ke = zeros(n_basefuncs, n_basefuncs) + Me = zeros(n_basefuncs, n_basefuncs) fe = zeros(n_basefuncs) - + #Initiate assembler for matrices and vector assembler = start_assemble(K, f) - + assembler_M = start_assemble(M) + #Iterate over cells for cell in CellIterator(dh) - fill!(Ke, 0) + fill!(Me, 0) fill!(fe, 0) - reinit!(cellvalues, cell) - + #Assemble local contributions for q_point in 1:getnquadpoints(cellvalues) dΩ = getdetJdV(cellvalues, q_point) - for i in 1:n_basefuncs - v = shape_value(cellvalues, q_point, i) - ∇v = shape_gradient(cellvalues, q_point, i) - fe[i] += 0.1 * v * dΩ + δNᵢ = shape_value(cellvalues, q_point, i) + ∇δNᵢ = shape_gradient(cellvalues, q_point, i) + fe[i] += 0.1 * δNᵢ * dΩ for j in 1:n_basefuncs - ∇u = shape_gradient(cellvalues, q_point, j) - Ke[i, j] += 1.0e-3 * (∇v ⋅ ∇u) * dΩ + Nⱼ = shape_value(cellvalues, q_point, j) + ∇Nⱼ = shape_gradient(cellvalues, q_point, j) + Ke[i, j] += 1.0e-3 * (∇δNᵢ ⋅ ∇Nⱼ) * dΩ + Me[i, j] += (δNᵢ * Nⱼ) * dΩ end end end - + #Update global matrices assemble!(assembler, celldofs(cell), Ke, fe) + assemble!(assembler_M, celldofs(cell), Me) end - return K, f + return K, M, f end #md nothing # hide -# In addition to the diffusive part, we also need a function that assembles the mass matrix `M`. -function doassemble_M!(M::SparseMatrixCSC, cellvalues::CellValues, dh::DofHandler) - - n_basefuncs = getnbasefunctions(cellvalues) - Me = zeros(n_basefuncs, n_basefuncs) - - assembler = start_assemble(M) - for cell in CellIterator(dh) - - fill!(Me, 0) - - reinit!(cellvalues, cell) - - for q_point in 1:getnquadpoints(cellvalues) - dΩ = getdetJdV(cellvalues, q_point) - - for i in 1:n_basefuncs - v = shape_value(cellvalues, q_point, i) - for j in 1:n_basefuncs - u = shape_value(cellvalues, q_point, j) - Me[i, j] += (v * u) * dΩ - end - end - end - - assemble!(assembler, celldofs(cell), Me) - end - return M -end -#md nothing # hide # ### Solution of the system # We first assemble all parts in the prior allocated `SparseMatrixCSC`. -K, f = doassemble_K!(K, f, cellvalues, dh) -M = doassemble_M!(M, cellvalues, dh) +K, M, f = doassemble!(K, M, f, cellvalues, dh) A = (Δt .* K) + M; # Now, we need to save all boundary condition related values of the unaltered system matrix `A`, which is done # by `get_rhs_data`. The function returns a `RHSData` struct, which contains all needed information to apply diff --git a/docs/src/topics/assembly.md b/docs/src/topics/assembly.md index e0f052d660..364e7fca14 100644 --- a/docs/src/topics/assembly.md +++ b/docs/src/topics/assembly.md @@ -111,8 +111,8 @@ using Ferrite ip = Lagrange{RefTriangle, 2}() # DofHandler -const N = 100 -grid = generate_grid(Triangle, (N, N)) +const N_grid = 100 +grid = generate_grid(Triangle, (N_grid, N_grid)) const dh = DofHandler(grid) add!(dh, :u, ip) close!(dh) @@ -253,7 +253,7 @@ and extract the dofs for that element. Finally, an assembler is created with ```@example assembly-perf dofs_per_cell = ndofs_per_cell(dh) const Ke = rand(dofs_per_cell, dofs_per_cell) -const dofs = celldofs(dh, N * N ÷ 2) +const dofs = celldofs(dh, N_grid * N_grid ÷ 2) const assembler = start_assemble(K) nothing # hide @@ -320,10 +320,10 @@ function assemble_system!(assembler_function::F, K, dh, cv) where {F} for qp in 1:getnquadpoints(cv) dΩ = getdetJdV(cv, qp) for i in 1:n - ∇ϕi = shape_gradient(cv, qp, i) + ∇δNᵢ = shape_gradient(cv, qp, i) for j in 1:n - ∇ϕj = shape_gradient(cv, qp, j) - ke[i, j] += ( ∇ϕi ⋅ ∇ϕj ) * dΩ + ∇Nⱼ = shape_gradient(cv, qp, j) + ke[i, j] += ( ∇δNᵢ ⋅ ∇Nⱼ ) * dΩ end end end diff --git a/docs/src/topics/boundary_conditions.md b/docs/src/topics/boundary_conditions.md index dbd03b7dec..d578481b59 100644 --- a/docs/src/topics/boundary_conditions.md +++ b/docs/src/topics/boundary_conditions.md @@ -135,8 +135,8 @@ for fc in FacetIterator(dh, getfacetset(grid, "right")) for q_point in 1:getnquadpoints(fv) dΓ = getdetJdV(fv, q_point) for i in 1:getnbasefunctions(fv) - δu = shape_value(fv, q_point, i) - fe[i] += δu * qn * dΓ + δNᵢ = shape_value(fv, q_point, i) + fe[i] += δNᵢ * qn * dΓ end end assemble!(f, celldofs(fc), fe) @@ -149,7 +149,7 @@ through the local `fe` vector and then using `assemble!`): # ... dofs = celldofs(fc) for i in 1:getnbasefunctions(fv) - f[dofs[i]] += δu * qn * dΓ + f[dofs[i]] += δN * qn * dΓ end ``` @@ -164,8 +164,8 @@ for facet in 1:nfacets(cell) for q_point in 1:getnquadpoints(facetvalues) dΓ = getdetJdV(facetvalues, q_point) for i in 1:getnbasefunctions(facetvalues) - δu = shape_value(facetvalues, q_point, i) - fe[i] += δu * qn * dΓ + δNᵢ = shape_value(facetvalues, q_point, i) + fe[i] += δNᵢ * qn * dΓ end end end diff --git a/docs/src/topics/fe_intro.md b/docs/src/topics/fe_intro.md index 0c44f944f6..3512bd05e7 100644 --- a/docs/src/topics/fe_intro.md +++ b/docs/src/topics/fe_intro.md @@ -47,16 +47,16 @@ simplicity we will consider only constant conductivity $k$. ## Weak form The solution to the equation above is usually calculated from the corresponding weak -form. By multiplying the equation with an arbitrary *test function* $\delta u$, integrating +form. By multiplying the equation with an arbitrary *test function* $\varphi$, integrating over the domain and using partial integration we obtain the *weak form*. Now our problem can be stated as: Find $u \in \mathbb{U}$ s.t. ```math -\int_\Omega \nabla \delta u \cdot (k \nabla u) \, \mathrm{d}\Omega = -\int_{\Gamma_\mathrm{N}} \delta u \, q^\mathrm{p} \, \mathrm{d}\Gamma + -\int_\Omega \delta u \, f \, \mathrm{d}\Omega \quad \forall \, \delta u \in \mathbb{T} +\int_\Omega \nabla \varphi \cdot (k \nabla u) \, \mathrm{d}\Omega = +\int_{\Gamma_\mathrm{N}} \varphi \, q^\mathrm{p} \, \mathrm{d}\Gamma + +\int_\Omega \varphi \, f \, \mathrm{d}\Omega \quad \forall \, \varphi \in \mathbb{T} ``` where $\mathbb{U}, \mathbb{T}$ are suitable function spaces with sufficiently regular @@ -80,14 +80,14 @@ Ferrite supports different approximations of the finite element fields and the g such basic numbering. For more details, see the [Ferrite numbering rules](@ref "Global-DoF-indices"). Note that *shape functions* are sometimes referred to as *basis functions* or *trial functions*, -and instead of $\phi_i$ they are sometimes denoted $N_i$. In this example we choose to approximate +and instead of $\phi_i$ they are sometimes denoted by $N_i$. In this example we choose to approximate the test function in the same way. This approach is known as the *Bubnov-Galerkin finite element method*. Formally we write the evaluation of our approximations at a specific point $\mathbf{x}$ in our domain $\Omega$ as: ```math u_\mathrm{h}(\mathbf{x}) = \sum_{i=1}^{\mathrm{N}} \phi_i(\mathbf{x}) \, \hat{u}_i,\qquad -\delta u_\mathrm{h}(\mathbf{x}) = \sum_{i=1}^{\mathrm{N}} \phi_i(\mathbf{x}) \, \delta \hat{u}_i \, . +\varphi_\mathrm{h}(\mathbf{x}) = \sum_{i=1}^{\mathrm{N}} \phi_i(\mathbf{x}) \, \delta \hat{u}_i \, . ``` Since test and trial functions are usually chosen in such a way, that they build the basis of @@ -101,7 +101,7 @@ We may now insert these approximations in the weak form, which results in \int_{\Omega_\mathrm{h}} \phi_i \, f \, \mathrm{d}\Omega \right) \, . ``` -Since this equation must hold for arbitrary $\delta u_\mathrm{h}$, the equation must especially +Since this equation must hold for arbitrary $\varphi_\mathrm{h}$, the equation must especially hold for the specific choice that only one of the nodal values $\delta \hat{u}_i$ is fixed to 1 while an all other coefficients are fixed to 0. Repeating this argument for all $i$ from 1 to N we obtain N linear equations. This way the discrete problem can be written as a system of linear equations @@ -169,6 +169,11 @@ being the chosen approximation when changing from the integral to the finite sum For an example of the implementation to solve a heat problem with Ferrite check out [this thoroughly commented example](@ref tutorial-heat-equation). +### Notation (tutorials) +For a function, `u` such that $u \approx \sum N_i a_i$, the corresponding test function is represented as `δu` such that $\delta u \approx \sum \delta N_i c_i$. In code, `N` / `δN` denote the `shape_value` for the function interpolation. +In coupled problems, the function's name is appended to the shape functions, i.e., `Nu` / `δNu`. +For assembly, these are subscripted with `i` or `j` based on the loop's variable, as in `Nuᵢ` / `δNuᵢ`. + ## More details We finally want to note that this quick introduction barely scratches the surface of the finite element