Package hygiene (#82)

* package hygiene

* add `Mass` and `Stiffness` types

Author: dkarrasch

Project.toml

@@ -1,6 +1,6 @@
 name = "CoherentStructures"
 uuid = "0c1513b4-3a13-56f1-9cd2-8312eaec88c4"
-version = "0.4.4"
+version = "0.4.5"
 
 [deps]
 ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"

docs/Manifest.toml

@@ -48,7 +48,7 @@ ctx = regularTriangularGrid((100,100), [0.0,0.0],[2π,2π])
 pred  = (x,y) -> ((x[1] - y[1]) % 2π) < 1e-9 && ((x[2] - y[2]) % 2π) < 1e-9
 bdata = BoundaryData(ctx,pred)
 
-id2 = one(Tensors.Tensor{2,2}) # 2D identity tensor
+id2 = one(Tensor{2,2}) # 2D identity tensor
 cgfun = x -> 0.5*(id2 +  dott(inv(DstandardMap(x))))
 
 K = assembleStiffnessMatrix(ctx,cgfun,bdata=bdata)

docs/makeimages.jl

@@ -162,7 +162,8 @@ u[20] = 2.0; u[38] = 3.0; u[56] = 4.0
 plot_u(ctx, u, 200, 200, bdata=bdata, colorbar=:none)
 ```
 
-To apply boundary conditions to a stiffness/mass matrix, use the `applyBCS` function. Note that `assembleStiffnessMatrix` and `assembleMassMatrix` take a `bdata` argument that does this internally.
+To apply boundary conditions to a stiffness/mass matrix, use the `applyBCS` function. Note
+that both `assemble` functions take a `bdata` argument that does this internally.
 
 ## Plotting and Videos
 
@@ -172,8 +173,9 @@ The simplest way to plot is using the [`plot_u`](@ref) function. Plots and video
 
 ## Parallelisation
 
-Many of the plotting functions support parallelism internally.
-Tensor fields can be constructed in parallel, and then passed to [`assembleStiffnessMatrix`](@ref). For an example that does this, see
+Many of the plotting functions support parallelism internally. Tensor fields can be
+constructed in parallel, and then passed to [`assemble(Stiffness(), _)`](@ref). For an
+example that does this, see
 TODO: Add this example
 
 ## FEM-API
@@ -185,8 +187,7 @@ CurrentModule = CoherentStructures
 ### Stiffness and Mass Matrices
 
 ```@docs
-assembleStiffnessMatrix
-assembleMassMatrix
+assemble
 adaptiveTOCollocationStiffnessMatrix
 ```
 

docs/src/fem.md

@@ -2,42 +2,41 @@ module CoherentStructures
 
 # use standard libraries
 using LinearAlgebra
-using ProgressMeter
+using LinearAlgebra: checksquare
+import ProgressMeter
 using SparseArrays
+using SparseArrays: nzvalview, dropzeros!
 using Distributed
 using SharedArrays: SharedArray
 using Statistics: mean
 
 # import data type packages
-import StaticArrays
-using StaticArrays: SVector, @SVector, SArray, SMatrix, @SMatrix
-import Tensors
-using Tensors: Vec, Tensor, SymmetricTensor
-import AxisArrays
+using StaticArrays: SVector, @SVector, SMatrix, @SMatrix
+using Tensors: Vec, Tensor, SymmetricTensor, basevec, dott, tdot, otimes, symmetric
 using AxisArrays: AxisArray, ClosedInterval, axisvalues
 
-import DiffEqBase
-import OrdinaryDiffEq
-const ODE = OrdinaryDiffEq
-import Distances
-const Dists = Distances
-import NearestNeighbors
-const NN = NearestNeighbors
-import Interpolations
-const ITP = Interpolations
+using DiffEqBase: DiffEqBase, initialize!, isconstant, update_coefficients!, @..
+using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, ODEFunction, ContinuousCallback,
+    terminate!, solve, Tsit5, BS5, OrdinaryDiffEqNewtonAlgorithm, DEFAULT_LINSOLVE,
+    alg_order, OrdinaryDiffEqMutableCache, alg_cache, @muladd, perform_step!, @unpack,
+    unwrap_alg, is_mass_matrix_alg
+
+using Distances: Distances, PreMetric, SemiMetric, Metric, Euclidean, PeriodicEuclidean,
+    pairwise, pairwise!, colwise, colwise!, result_type
+
+using NearestNeighbors: BallTree, KDTree, inrange, knn, MinkowskiMetric
+
+using Interpolations: Interpolations, LinearInterpolation, CubicSplineInterpolation,
+    interpolate, scale, BSpline, Linear, Cubic, Natural, OnGrid, Free
 
 # import linear algebra related packages
-import LinearMaps
-const LMs = LinearMaps
-import IterativeSolvers
-import ArnoldiMethod
+using LinearMaps: LinearMap
+using IterativeSolvers: cg
+using ArnoldiMethod: partialschur, partialeigen
 
 # import geometry related packages
-import GeometricalPredicates
-const GP = GeometricalPredicates
-import VoronoiDelaunay
-const VD = VoronoiDelaunay
-
+using GeometricalPredicates: GeometricalPredicates, Point, AbstractPoint2D, Point2D, getx, gety
+using VoronoiDelaunay: DelaunayTessellation2D, findindex, isexternal, max_coord, min_coord
 
 # Ferrite
 import Ferrite

src/CoherentStructures.jl

@@ -1,40 +1,45 @@
 # Strongly inspired by an example provided on Ferrite's github page, modified and
 # extended by Nathanael Schilling
 
+struct Stiffness end
+struct Mass end
+
 #Works in n=2 and n=3
-tensorIdentity(x::Vec{dim}, _, p) where {dim} = one(SymmetricTensor{2,dim,Float64,3(dim-1)})
+tensorIdentity(x::Vec{dim,T}, _, p) where {dim,T} = one(SymmetricTensor{2,dim,T})
 
 """
-    assembleStiffnessMatrix(ctx, A, p=nothing; bdata=BoundaryData())
+    assemble(Stiffness(), ctx; A=Id, p=nothing, bdata=BoundaryData())
 
 Assemble the stiffness-matrix for a symmetric bilinear form
 ```math
 a(u,v) = \\int \\nabla u(x)\\cdot A(x)\\nabla v(x)f(x) dx.
 ```
 The integral is approximated using numerical quadrature. `A` is a function that returns a
 `SymmetricTensor{2,dim}` object and must have one of the following signatures:
-   * `A(x::Vector{Float64})`;
-   * `A(x::Vec{dim})`;
-   * `A(x::Vec{dim}, index::Int, p)`. Here, `x` is equal to `ctx.quadrature_points[index]`,
-     and `p` is the one passed to `assembleStiffnessMatrix`.
+
+* `A(x::Vector{Float64})`;
+* `A(x::Vec{dim})`;
+* `A(x::Vec{dim}, index::Int, p)`. Here, `x` is equal to `ctx.quadrature_points[index]`,
+  and `p` is some parameter, think of some precomputed object that is indexed via `index`.
 
 The ordering of the result is in dof order, except that boundary conditions from `bdata` are
 applied. The default is natural (homogeneous Neumann) boundary conditions.
 """
-function assembleStiffnessMatrix(ctx::GridContext, A=tensorIdentity, p=nothing; bdata=BoundaryData())
+function assemble(::Stiffness, ctx::GridContext; A=tensorIdentity, p=nothing, bdata=BoundaryData())
     if A === tensorIdentity
-        return _assembleStiffnessMatrix(ctx, A, p, bdata=bdata)
+        return _assembleStiffnessMatrix(ctx, A, p, bdata = bdata)
     elseif !isempty(methods(A, (Vec,)))
         return _assembleStiffnessMatrix(ctx, (qp, i, p) -> A(qp), p, bdata=bdata)
     elseif !isempty(methods(A, (Vec, Int, Any)))
         return _assembleStiffnessMatrix(ctx, (qp, i, p) -> A(qp, i, p), p, bdata=bdata)
     elseif !isempty(methods(A, (Vector{Float64},)))
-        return _assembleStiffnessMatrix(ctx, (qp, i, p) -> A(Vector{Float64}(qp)), p, bdata=bdata)
+        return _assembleStiffnessMatrix(ctx, (qp, i, p) -> A(convert(Vector{Float64}, qp)), p, bdata=bdata)
     end
-    error("Function parameter A does not accept types supported by assembleStiffnessMatrix")
+    error("function argument `A` does not admit any of the accepted signatures")
 end
+@deprecate assembleStiffnessMatrix(ctx::GridContext, A=tensorIdentity, p=nothing; bdata=BoundaryData()) assemble(Stiffness(), ctx; A=A, p=p, bdata=bdata)
 
-function _assembleStiffnessMatrix(ctx::GridContext, A, p; bdata=BoundaryData())
+function _assembleStiffnessMatrix(ctx, A, p; bdata=BoundaryData())
     cv = FEM.CellScalarValues(ctx.qr, ctx.ip, ctx.ip_geom)
     dh = ctx.dh
     K = FEM.create_sparsity_pattern(dh)
@@ -70,7 +75,7 @@ end
 
 
 """
-    assembleMassMatrix(ctx; bdata=BoundaryData(), lumped=false)
+    assemble(Mass(), ctx; bdata=BoundaryData(), lumped=false)
 
 Assemble the mass matrix
 ```math
@@ -86,10 +91,10 @@ Returns a lumped mass matrix if `lumped=true`.
 # Example
 ```
 ctx.mass_weights = map(f, ctx.quadrature_points)
-M = assembleMassMatrix(ctx)
+M = assemble(Mass(), ctx)
 ```
 """
-function assembleMassMatrix(ctx::GridContext; bdata=BoundaryData(), lumped=false)
+function assemble(::Mass, ctx::GridContext; bdata=BoundaryData(), lumped=false)
     cv = FEM.CellScalarValues(ctx.qr, ctx.ip, ctx.ip_geom)
     dh = ctx.dh
     M = FEM.create_sparsity_pattern(dh)
@@ -123,12 +128,12 @@ function assembleMassMatrix(ctx::GridContext; bdata=BoundaryData(), lumped=false
     M = applyBCS(ctx, M, bdata)
     return lumped ? spdiagm(0 => dropdims(reduce(+, M; dims=1); dims=1)) : M
 end
+@deprecate assembleMassMatrix(ctx::GridContext; bdata=BoundaryData(), lumped=false) assemble(Mass(), ctx; bdata=bdata, lumped=lumped)
 
 """
-    getQuadPointsPoints(ctx)
+    getQuadPoints(ctx)
 
-Compute the coordinates of all quadrature points on a grid.
-Helper function.
+Compute the coordinates of all quadrature points on a grid. Helper function.
 """
 function getQuadPoints(ctx::GridContext{dim}) where {dim}
     cv = FEM.CellScalarValues(ctx.qr, ctx.ip_geom)

src/FEMassembly.jl

@@ -50,11 +50,11 @@ function nonAdaptiveTOCollocation(
 
     #Calculate the integral of shape function in the domain (dof order)
     shape_function_weights_domain_original = undoBCS(ctx_domain,
-            vec(sum(assembleMassMatrix(ctx_domain, bdata=bdata_domain), dims=1)),
+            vec(sum(assemble(Mass(), ctx_domain, bdata=bdata_domain), dims=1)),
             bdata_domain)
     #Calculate the integral of shape functions in the codomain (dof order)
     #Do not include boundary conditions here, as we end up summing over this later
-    shape_function_weights_codomain = vec(sum(assembleMassMatrix(ctx_codomain),dims=1))
+    shape_function_weights_codomain = vec(sum(assemble(Mass(), ctx_codomain), dims=1))
 
     #This variable will hold approximate pullback (in the measure sense)
     #of shape_function_weights_codomain to domain via finv
@@ -102,23 +102,28 @@ function nonAdaptiveTOCollocation(
 end
 
 """
-    adaptiveTOCollocationStiffnessMatrix(ctx,flow_maps,times=nothing; [quadrature_order, on_torus,on_cylinder, LL, UR, bdata, volume_preserving=true,flow_map_mode=0] )
+    adaptiveTOCollocationStiffnessMatrix(ctx, flow_maps, times=nothing; [quadrature_order, on_torus,on_cylinder, LL, UR, bdata, volume_preserving=true, flow_map_mode=0] )
 
-Calculate the matrix-representation of the bilinear form ``a(u,v) = 1/N \\sum_n^N a_1(I_hT_nu,I_hT_nv)`` where
-``I_h`` is pointwise interpolation of the grid obtained by doing Delaunay triangulation on images of grid points from ctx
-and ``T_n`` is the Transfer-operator for ``x \\mapsto flow_maps(x,times)[n]`` and ``a_1`` is the weak form of the Laplacian on the codomain. Moreover,
-``N`` in the equation above is equal to `length(times)` and ``t_n`` ranges over the elements of `times`.
+Calculate the matrix-representation of the bilinear form ``a(u,v) = 1/N \\sum_n^N a_1(I_hT_nu,I_hT_nv)``
+where ``I_h`` is pointwise interpolation of the grid obtained by doing Delaunay
+triangulation on images of grid points from `ctx` and ``T_n`` is the transfer operator for
+``x \\mapsto flow_maps(x, times)[n]`` and ``a_1`` is the weak form of the Laplacian on the
+codomain. Moreover, ``N`` in the equation above is equal to `length(times)` and ``t_n``
+ranges over the elements of `times`.
 
-If `times==nothing`, take ``N=1`` above and use the map ``x \\mapsto flow_maps(x)` instead of the version with `t_n`.
+If `times==nothing`, take ``N=1`` above and use the map ``x \\mapsto flow_maps(x)`` instead
+of the version with `t_n`.
 
-If `on_torus` is true, the Delaunay Triangulation is done on the torus.
-If `on_cylinder` is true, then triangulation is done on cylinder (periodic) x. In both of these cases we require `bdata` for boundary information
-on the original domain as well as `LL` and `UR` as lower-left and upper-right corners of the image.
+If `on_torus` is true, the Delaunay triangulation is done on the torus.
+If `on_cylinder` is true, then triangulation is done on an x-periodic cylinder. In both of
+these cases we require `bdata` for boundary information on the original domain as well as
+`LL` and `UR` as lower-left and upper-right corners of the image.
 
 If `volume_preserving == false`, add a volume_correction term to ``a_1`` (See paper by Froyland & Junge).
 
-If `flow_map_mode==0`, apply flow map to nodal basis function coordinates.
-If `flow_map_mode==1`, apply flow map to nodal basis function index number (allows for precomputed trajectories).
+If `flow_map_mode=0`, apply flow map to nodal basis function coordinates.
+If `flow_map_mode=1`, apply flow map to nodal basis function index number (allows for
+precomputed trajectories).
 """
 function adaptiveTOCollocationStiffnessMatrix(ctx::GridContext{2}, flow_maps, times=nothing;
                                         quadrature_order=default_quadrature_order,
@@ -183,58 +188,52 @@ function adaptiveTOCollocationStiffnessMatrix(ctx::GridContext{2}, flow_maps, ti
             #the old grid
             new_density_nodevals = new_density_bcdofvals
             while length(new_density_nodevals) != new_ctx.n
-                push!(new_density_nodevals,0.0)
+                push!(new_density_nodevals, 0.0)
             end
 
-            new_ctx.mass_weights = [evaluate_function_from_node_or_cellvals(new_ctx,new_density_nodevals,q)
-                       for q in new_ctx.quadrature_points ]
+            new_ctx.mass_weights = [evaluate_function_from_node_or_cellvals(new_ctx, new_density_nodevals, q)
+                for q in new_ctx.quadrature_points]
         end
-        I, J, V = findnz(assembleStiffnessMatrix(new_ctx,bdata=new_bdata))
+        I, J, V = findnz(assemble(Stiffness(), new_ctx, bdata=new_bdata))
         I .= translation_table[I]
         J .= translation_table[J]
         n = ctx.n - length(bdata.periodic_dofs_from)
-        push!(As,sparse(I,J,V,n,n))
+        push!(As, sparse(I, J, V, n, n))
     end
     return mean(As)
 end
 
 
 #Reordering in the periodic case is slightly more tricky
-function bcdofNewToBcdofOld(
-        old_ctx::GridContext{dim},bdata::BoundaryData,
-        new_ctx::GridContext{dim},new_bdata::BoundaryData,K
-        ) where {dim}
-
-        I, J ,V = findnz(K)
-
-        #Here I,J are in pdof order for new_ctx
-
-        bcdof_to_node_new = bcdof_to_node(new_ctx,new_bdata)
-        node_to_bcdof_old = node_to_bcdof(old_ctx,bdata)
-        I .= node_to_bcdof_old[bcdof_to_node_new[I]]
-        J .= node_to_bcdof_old[bcdof_to_node_new[J]]
-
-        old_pdof_n = old_ctx.n - length(bdata.periodic_dofs_from)
-
-        return sparse(I,J,V,old_pdof_n,old_pdof_n)
+function bcdofNewToBcdofOld(old_ctx::GridContext{dim}, bdata::BoundaryData,
+                            new_ctx::GridContext{dim}, new_bdata::BoundaryData, K) where {dim}
+    I, J ,V = findnz(K)
+
+    # Here I,J are in pdof order for new_ctx
+    bcdof_to_node_new = bcdof_to_node(new_ctx, new_bdata)
+    node_to_bcdof_old = node_to_bcdof(old_ctx, bdata)
+    I .= node_to_bcdof_old[bcdof_to_node_new[I]]
+    J .= node_to_bcdof_old[bcdof_to_node_new[J]]
+    old_pdof_n = old_ctx.n - length(bdata.periodic_dofs_from)
+    return sparse(I, J, V, old_pdof_n, old_pdof_n)
 end
 
-function node_to_bcdof(ctx::GridContext{dim},bdata::BoundaryData) where {dim}
-        n_nodes = ctx.n - length(bdata.periodic_dofs_from)
-        bdata_table = BCTable(ctx,bdata)
-        return bdata_table[ctx.node_to_dof]
+function node_to_bcdof(ctx::GridContext{dim}, bdata::BoundaryData) where {dim}
+    n_nodes = ctx.n - length(bdata.periodic_dofs_from)
+    bdata_table = BCTable(ctx, bdata)
+    return bdata_table[ctx.node_to_dof]
 end
 
 function bcdof_to_node(ctx::GridContext{dim},bdata::BoundaryData) where {dim}
-        n_nodes = ctx.n - length(bdata.periodic_dofs_from)
-        bdata_table = BCTable(ctx,bdata)
-        result = zeros(Int,n_nodes)
-        for i in 1:ctx.n
-            if result[bdata_table[ctx.node_to_dof[i]]] == 0
-                result[bdata_table[ctx.node_to_dof[i]]] = i
-            end
+    n_nodes = ctx.n - length(bdata.periodic_dofs_from)
+    bdata_table = BCTable(ctx,bdata)
+    result = zeros(Int,n_nodes)
+    for i in 1:ctx.n
+        if iszero(result[bdata_table[ctx.node_to_dof[i]]])
+            result[bdata_table[ctx.node_to_dof[i]]] = i
         end
-        return result
+    end
+    return result
 end
 
 
@@ -274,8 +273,8 @@ function adaptiveTOFutureGrid(ctx::GridContext{dim}, flow_map;
     #Do volume corrections
     #All values are in node order for new_ctx, which is the same as bcdof order for ctx2
     n_nodes = length(new_nodes_in_bcdof_order)
-    vols_new = sum(assembleMassMatrix(new_ctx, bdata=new_bdata), dims=1)[1, node_to_bcdof(new_ctx, new_bdata)[1:n_nodes]]
-    vols_old = sum(assembleMassMatrix(ctx, bdata=bdata), dims=1)[1, 1:n_nodes]
+    vols_new = sum(assemble(Mass(), new_ctx, bdata=new_bdata), dims=1)[1, node_to_bcdof(new_ctx, new_bdata)[1:n_nodes]]
+    vols_old = sum(assemble(Mass(), ctx, bdata=bdata), dims=1)[1, 1:n_nodes]
 
     return new_ctx, new_bdata, vols_new ./ vols_old
 end
@@ -319,12 +318,12 @@ function adaptiveTOCollocation(ctx::GridContext{dim}, flow_map;
             throw(AssertionError("Invalid projection_method"))
         end
         if volume_preserving
-            L = sparse(I, npoints, npoints)[node_to_bcdof(ctx,bdata)[bcdof_to_node(ctx_new,bdata_new)],:]
+            L = sparse(I, npoints, npoints)[node_to_bcdof(ctx, bdata)[bcdof_to_node(ctx_new,bdata_new)],:]
             result = ALPHA_bc*L
         else
             volume_change_pdof = volume_change[bcdof_to_node(ctx,bdata)]
-            K = assembleStiffnessMatrix(ctx,bdata=bdata)
-            M = assembleMassMatrix(ctx,bdata=bdata)
+            K = assemble(Stiffness(), ctx, bdata = bdata)
+            M = assemble(Mass(), ctx, bdata = bdata)
             volume_change_pdof = (M - 1e-2K)\(M*volume_change_pdof)
             volume_change = volume_change_pdof[node_to_bcdof(ctx,bdata)]