Minor changes in transient_heat_equation.jl and change in notation of test functions

docs/src/literate-tutorials/heat_equation.jl

@@ -30,9 +30,9 @@
 # where $\partial \Omega$ denotes the boundary of $\Omega$.
 # The resulting weak form is given as follows: Find ``u \in \mathbb{U}`` such that
 # ```math
-# \int_{\Omega} \nabla \delta u \cdot \nabla u \ d\Omega = \int_{\Omega} \delta u \ d\Omega \quad \forall \delta u \in \mathbb{T},
+# \int_{\Omega} \nabla \varphi \cdot \nabla u \ d\Omega = \int_{\Omega} \varphi \ d\Omega \quad \forall \varphi \in \mathbb{T},
 # ```
-# where $\delta u$ is a test function, and where $\mathbb{U}$ and $\mathbb{T}$ are suitable
+# where $\varphi$ is a test function, and where $\mathbb{U}$ and $\mathbb{T}$ are suitable
 # trial and test function sets, respectively.
 #-
 # ## Commented program
@@ -127,12 +127,12 @@ close!(ch)
 # then reused for all elements) so we first need to make sure that they are all zeroes at
 # the start of the function by using `fill!`. Then we loop over all the quadrature points,
 # and for each quadrature point we loop over all the (local) shape functions. We need the
-# value and gradient of the test function, `δu` and also the gradient of the trial function
+# value and gradient of the test function, `ϕ` and also the gradient of the trial function
 # `u`. We get all of these from `cellvalues`.
 #
 # !!! 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 `ϕ`, `∇ϕ` and `∇u` 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
@@ -149,15 +149,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)
+            ϕ = shape_value(cellvalues, q_point, i)
+            ∇ϕ = shape_gradient(cellvalues, q_point, i)
             ## Add contribution to fe
-            fe[i] += δu * dΩ
+            fe[i] += ϕ * dΩ
             ## Loop over trial shape functions
             for j in 1:n_basefuncs
                 ∇u = shape_gradient(cellvalues, q_point, j)
                 ## Add contribution to Ke
-                Ke[i, j] += (∇δu ⋅ ∇u) * dΩ
+                Ke[i, j] += (∇ϕ ⋅ ∇u) * dΩ
             end
         end
     end

docs/src/literate-tutorials/transient_heat_equation.jl

@@ -35,22 +35,29 @@
 # ```
 # 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} \varphi_i \, u_j \ \mathrm{d}\Omega,
+# ```
+# where $u_j$ and $\varphi_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 fe-intro).
 #-
 # ## Commented program
 #
@@ -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 the [heat equation example](@ref heat_equation.jl) 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Ω
+                ϕ = shape_value(cellvalues, q_point, i)
+                ∇ϕ = shape_gradient(cellvalues, q_point, i)
+                fe[i] += 0.1 * ϕ * dΩ
                 for j in 1:n_basefuncs
+                    u = shape_value(cellvalues, q_point, j)
                     ∇u = shape_gradient(cellvalues, q_point, j)
-                    Ke[i, j] += 1.0e-3 * (∇v ⋅ ∇u) * dΩ
+                    Ke[i, j] += 1.0e-3 * (∇ϕ ⋅ ∇u) * dΩ
+                    Me[i, j] += (ϕ * u) * 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

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Γ
+            ϕ = shape_value(fv, q_point, i)
+            fe[i] += ϕ * 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]] += ϕ * 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Γ
+                ϕ = shape_value(facetvalues, q_point, i)
+                fe[i] += ϕ * qn * dΓ
             end
         end
     end

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
@@ -87,7 +87,7 @@ 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