Add tensor_view function and show usage in 2D-1 benchmark example.

docs/make.jl

@@ -61,7 +61,7 @@ function make_all(; with_examples::Bool = true, run_examples::Bool = true, run_n
 	makedocs(
 		modules = [ExtendableFEM],
 		sitename = "ExtendableFEM.jl",
-		authors = "Christian Merdon",
+		authors = "Christian Merdon, Jan Philipp Thiele",
 		repo = "github.com/chmerdon/ExtendableFEM.jl",
 		clean = false,
 		checkdocs = :all,
@@ -71,6 +71,7 @@ function make_all(; with_examples::Bool = true, run_examples::Bool = true, run_n
 			"Index" => "package_index.md",
 			"Problem Description" => [
 				"problemdescription.md",
+                "tensordescription.md",
 				"nonlinearoperator.md",
 				"bilinearoperator.md",
 				"linearoperator.md",

docs/src/nonlinearoperator.md

@@ -7,6 +7,14 @@ constructed with special constructors for BilinearOperator or LinearOperator.
 
 ## Constructor
 
+To describe a NonlinearOperator we have to specify a kernel function. 
+These functions are 'flat' in the sense that the input and output vector 
+contain the components of the test-function values and derivatives
+as specified by `oa_test` and `oa_args` respectively.
+The assembly of the local matrix will be done internally 
+by multiplying the subvectors of result with its test-function counterparts.
+For a more detailed explanation of this see the following
+
 ```@autodocs
 Modules = [ExtendableFEM]
 Pages = ["common_operators/nonlinear_operator.jl"]
@@ -15,15 +23,76 @@ Order   = [:type, :function]
 
 ## Example - NSE convection operator
 
-For the 2D Navier--Stokes equations, a kernel function for the convection operator could look like
+
+For the Navier--Stokes equations, we need a kernel function for the nonlinear
+convection term
+```math
+\begin{equation}
+(v,u\cdot\nabla u) = (v,\nabla u^T u)
+\end{equation}
+```
+In 2D the input (as specified below) will contain the two
+components of ``u=(u_1,u_2)'`` and the four components of the gradient 
+``\nabla u = \begin{pmatrix} u_{11} & u_{12} \\ u_{21} & u_{22}\end{pmatrix}``
+in order, i.e. ``(u_1,u_2,u_{11},u_{12},u_{21},u_{22})``.
+As the convection term is tested with ``v``, 
+the ouptut vector ``o`` only has to contain what should be tested with each component
+of ``v``, i.e.
+```math
+\begin{equation}
+    A_\text{local} = (v_1,v_2)^T(o_1,o_2) = 
+        \begin{pmatrix}
+            v_1o_1 & v_1o_2\\
+            v_2o_1 & v_2o_2
+        \end{pmatrix}.
+\end{equation}
+```
+To construct the kernel there are two options, 
+component-wise and based on `tensor_view`.
+For the first we have to write the convection term as individual components
+```math
+\begin{equation}
+o = 
+    \begin{pmatrix}
+        u_1\cdot u_{11}+u_2\cdot u_{12}\\
+        u_1\cdot u_{21}+u_2\cdot u_{22}\\
+    \end{pmatrix}
+= 
+\begin{pmatrix}
+    u\cdot (u_11,u_12)^T\\
+    u\cdot (u_21,u_22)^T
+\end{pmatrix}.
+\end{equation}
+```
+To make our lives a bit easier we will extract the subcompontents of 
+`input` as views, such that `∇u[3]` actually accesses `input[5]`,
+which corresponds to the third entry ``u_{21}`` of ``\nabla u``. 
 ```julia
 function kernel!(result, input, qpinfo)
     u, ∇u = view(input, 1:2), view(input,3:6)
     result[1] = dot(u, view(∇u,1:2))
     result[2] = dot(u, view(∇u,3:4))
+    return nothing
 end
 ```
-and the coressponding NonlinearOperator constructor call reads
+To improve readability of the kernels and to make them easier to understand,
+we provide the function `tensor_view` which constructs a view and reshapes 
+it into an object matching the given `TensorDescription`.
+See the [table](@ref "Which tensor for which unknown?") 
+to see which tensor size is needed for which derivative of a scalar, vector 
+or matrix-valued variable.
+```julia
+function kernel!(result, input, qpinfo)
+    u = tensor_view(input,1,TDVector(2))
+    v = tensor_view(result,1,TDVector(2))
+    ∇u = tensor_view(input,3,TDMatrix(2))
+    tmul!(v,∇u,u)
+    return nothing
+end
+```
+
+The coressponding NonlinearOperator constructor call is the same in both cases 
+and reads
 ```julia
 u = Unknown("u"; name = "velocity")
 NonlinearOperator(kernel!, [id(u)], [id(u),grad(u)])
@@ -55,6 +124,21 @@ triggers that the ```result``` vector of the kernel is multiplied with the Ident
 
     Kernels are allowed to depend on region numbers, space and time coordinates via the qpinfo argument.
 
+!!! note "Dimension independent kernels"
+
+    If done correctly, the operator-based approach allows us to write a kernel 
+    that is 'independent' of the spatial dimension, 
+    i.e. one instead of up to three kernels.
+    Assuming `dim` is a known variable we can re-write the kernel from above as
+    ```julia
+    function kernel!(result, input, qpinfo)
+        u = tensor_view(input,1,TDVector(dim))
+        v = tensor_view(result,1,TDVector(dim))
+        ∇u = tensor_view(input,1+dim,TDMatrix(dim))
+        tmul!(v,∇u,u)
+        return nothing
+    end
+    ```
 
 ## Newton by local jacobians of kernel
 

docs/src/tensordescription.md

@@ -0,0 +1,68 @@
+
+# Tensor Description
+
+To be able to construct [reshaped views](@ref "Reshaped views") 
+of the test functions and their derivates, we can describe the 
+shape of the view through a [`TensorDescription{R,D}`](@ref ExtendableFEM.TensorDescription{R,D}) 
+where `R` is the *rank* of the tensor and `D` is the dimension 
+or extent of the tensor in each of the `R` directions. 
+That means a real valued `R`-tensor is an element of 
+``\underbrace{\mathbb{R}^D\times\cdots\times\mathbb{R}^D}_{R \text{ times}}``. 
+Specifically, we can identify the following mathematical objects with 
+tensors of different ranks:
+
+| math. object                                 | `R`-Tensor | 
+| :------------------------------------------- | :--------- |
+| scalar ``\in\mathbb{R}``                     | 0-Tensor   | 
+| vector ``\in\mathbb{R}^D``                   | 1-Tensor   | 
+| matrix ``\in\mathbb{R}^D\times\mathbb{R}^D`` | 2-Tensor   | 
+
+For finite elements, `D` usually matches the spatial dimension of 
+the problem we want to solve, i.e. `D=2` for 2D and `D=3` for 3D.
+
+## Tensor Types
+
+```@docs
+ExtendableFEM.TensorDescription
+ExtendableFEM.TDScalar
+ExtendableFEM.TDVector
+ExtendableFEM.TDMatrix
+ExtendableFEM.TDRank3
+ExtendableFEM.TDRank4
+```
+
+
+## Reshaped views
+
+```@autodocs
+Modules = [ExtendableFEM]
+Pages = ["tensors.jl"]
+Order   = [:function]
+```
+
+## Which tensor for which unknown?
+For an unknown variable `u` of tensor rank `r` 
+a derivative of order `n` has rank `r+n`,
+i.e. the hessian (n=2) of a scalar unknown (rank 0)
+and the gradient (n=1) of a vector valued (rank 1) 
+variable are both matrices (rank 2).
+
+For a more comprehensive list see the following table
+
+| derivative order   | scalar-valued    | vector-valued            | matrix-valued            |
+| :----------------- | :--------------- | :----------------------- | :----------------------- |
+| 0 (value/`id`)     | `TDScalar(D)`    | `TDVector(D)`            | `TDMatrix(D)`            |
+| 1 (`grad`)         | `TDVector(D)`    | `TDMatrix(D)`            | `TDRank3(D)`             |
+| 2 (`hessian`)      | `TDMatrix(D)`    | `TDRank3(D)`             | `TDRank4(D)`             |
+| 3                  | `TDRank3(D)`     | `TDRank4(D)`             | `TensorDescription(5,D)` |
+| 4                  | `TDRank4(D)`     | `TensorDescription(5,D)` | `TensorDescription(6,D)` |
+| ``\vdots``         | ``\vdots``       | ``\vdots``               | ``\vdots``               | 
+
+
+
+## Helpers
+
+
+```@docs
+tmul!
+```

examples/Example245_NSEFlowAroundCylinder.jl

@@ -43,19 +43,47 @@ function inflow!(result, qpinfo)
 	result[2] = 0.0
 end
 
+## hand constructed identity matrix for kernel to avoid allocations
+const II = [1 0; 0 1]
+
+
+## Example of a kernel using tensor_view() function to allow for an operator 
+## based style of writing the semilinear form.
+## For comparison we also provide the kernel_nonlinear_flat! function below 
+## that uses a component-wise style of writing the semilinear form.
+
+## 
+## the scalar product ``(\nabla v, \mu \nabla u - p)`` will be evaluated
+## so in general `a = b` corresponds to ``(a,b)``. 
+
+## Note that the order of vector entries between the kernel and the call to 
+## NonlinearOperator have to match.
 function kernel_nonlinear!(result, u_ops, qpinfo)
-	u, ∇u, p = view(u_ops, 1:2), view(u_ops, 3:6), view(u_ops, 7)
-	μ = qpinfo.params[1]
-	result[1] = dot(u, view(∇u, 1:2))
-	result[2] = dot(u, view(∇u, 3:4))
-	result[3] = μ * ∇u[1] - p[1]
-	result[4] = μ * ∇u[2]
-	result[5] = μ * ∇u[3]
-	result[6] = μ * ∇u[4] - p[1]
-	result[7] = -(∇u[1] + ∇u[4])
-	return nothing
+    # Shape values of vectorial u are starting at index 1
+    # view as 1-tensor(vector) of length dim=2 in 2D
+    u = tensor_view(u_ops, 1, TDVector(2))
+    v = tensor_view(result, 1, TDVector(2))
+    # gradients of vectorial u are starting at index 3
+    # view as 2-tensor of size 2x2 in 2D
+    ∇u = tensor_view(u_ops, 3, TDMatrix(2))
+    ∇v = tensor_view(result, 3, TDMatrix(2))
+    # values of scalar p are starting at index 7
+    # view as 0-tensor (single value) 
+    p = tensor_view(u_ops, 7, TDScalar())
+    q = tensor_view(result, 7, TDScalar())
+    # get viscosity at current quadrature point
+    μ = qpinfo.params[1]
+    # Note that all operators should be element-wise to avoid allocations
+    # `(v,u⋅∇u) = (v,∇u^T⋅u)`
+    tmul!(v,∇u,u)
+    # `(∇v,μ∇u-p)`
+    ∇v .= μ .* ∇u .- p[1] .* II
+    # `(q,-∇⋅u)`
+    q[1] = -dot(∇u, II)
+    return nothing
 end
 
+
 ## everything is wrapped in a main function
 function main(; Plotter = nothing, μ = 1e-3, maxvol = 1e-3, reconstruct = true, kwargs...)
 

src/ExtendableFEM.jl

@@ -73,6 +73,17 @@ export assign_reduction!
 
 include("helper_functions.jl")
 export get_periodic_coupling_info
+export tmul!
+
+include("tensors.jl")
+export TensorDescription
+export TDScalar
+export TDVector
+export TDMatrix
+export TDRank3
+export TDRank4
+export tensor_view
+
 
 include("solver_config.jl")
 export SolverConfiguration

src/helper_functions.jl

@@ -210,3 +210,68 @@ function get_dofmap(FES, xgrid, AT)
 	return FES.dofgrid === xgrid ? FES[DM] : FES[ParentDofmap4Dofmap(DM)]
 end
 
+
+# """
+# ````
+# function tensor_view(input, i, rank, dim)
+# ````
+
+# Returns a view of input[i] and following entries 
+# reshaped as a tensor of rank `rank`.
+# The parameter `dim` specifies the size of a tensor in each direction, 
+# e.g. a 1-Tensor (Vector) of length(dim) or a dim x dim 2-Tensor (matrix).
+# As an example `tensor_view(v,5,2,5)` returns a view of `v(5:29)`  
+# as a 5x5 matrix.
+
+# """
+# function tensor_view(input, i::Int, rank::Int, dim::Int)
+#     if rank == 0
+#         return view(input, i:i)
+#     elseif rank == 1
+#         return view(input, i:i+dim-1)
+#     elseif rank == 2
+#         return reshape(view(input, i:(i+(dim*dim)-1)), (dim,dim))
+#     elseif rank == 3
+#         return reshape(view(input, i:(i+(dim*dim*dim)-1)), (dim, dim,dim))
+#     else
+#         @warn "tensor_view for rank > 3 is a general implementation that needs allocations!"
+#         return reshape(view(input, i:(i+(dim^rank)-1)),ntuple(i->dim,rank))
+#     end
+# end
+
+
+"""
+````
+function tmul!(y,A,x,α=1.0,β=0.0)
+````
+
+Combined inplace  matrix-vector multiply-add ``A^T x α + y β``.
+The result is stored in `y` by overwriting it.  Note that `y` must not be
+aliased with either `A` or `x`.
+
+"""
+function tmul!(y, A, x, α = 1.0, β = 0.0)
+    for i in eachindex(y)
+        y[i] *= β
+        for j in eachindex(x)
+            y[i] += α * A[j, i] * x[j]
+        end
+    end
+end
+
+"""
+````
+function tmul!(y::AbstractVector{T}, A::AbstractMatrix{T}, x::AbstractVector{T}, α=1.0, β=0.0) where {T<:AbstractFloat}
+````
+
+Overload of the generic function for types supported by 
+`LinearAlgebra.BLAS.gemv!` to avoid slow run times for large inputs.
+"""
+function tmul!(
+    y::AbstractVector{T},
+    A::AbstractMatrix{T},
+    x::AbstractVector{T},
+    α=1.0,
+    β=0.0) where {T<:AbstractFloat}
+    LinearAlgebra.BLAS.gemv!('T', α, A, x, β, y)
+end