Finite-element discretizations of the dynamic Laplacian

assemble(Stiffness(), ctx; A=Id, p=nothing, bdata=BoundaryData())

Assemble the stiffness-matrix for a symmetric bilinear form

\[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 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.

assemble(Mass(), ctx; bdata=BoundaryData(), lumped=false)

Assemble the mass matrix

\[M_{i,j} = \int \varphi_j(x) \varphi_i(x) f(x)d\lambda^d\]

The integral is approximated using numerical quadrature. The values of f(x) are taken from ctx.mass_weights, and should be ordered in the same way as ctx.quadrature_points.

The result is ordered to be usable with a stiffness matrix with boundary data bdata.

Returns a lumped mass matrix if lumped=true.

Example

ctx.mass_weights = map(f, ctx.quadrature_points)
M = assemble(Mass(), ctx)

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.

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 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 volumecorrection term to ``a1`` (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).

regular1dP1Grid(numnodes, left=0.0, right=1.0; [quadrature_order])

Create a regular grid with numnodes nodes on the interval [left, right] in 1d, with P1-Lagrange basis functions.

regular1dP2Grid(numnodes, left=0.0, right=1.0; [quadrature_order])

Create a regular grid with numnodes non-interior nodes on the interval [left, right], with P2-Lagrange elements.

regular2dGrid(gridType, numnodes, LL=(0.0,0.0), UR=(1.0,1.0); quadrature_order=default_quadrature_order)

Constructs a regular grid. gridType should be one from regular2dGridTypes.

regularTriangularGrid(numnodes, LL, UR; [quadrature_order, PC=false])

Create a regular triangular grid on a rectangle; does not use Delaunay triangulation internally. If

regularDelaunayGrid(numnodes=(25,25), LL=(0.0,0.0), UR=(1.0,1.0); [quadrature_order, on_torus=false, on_cylinder=false, nudge_epsilon=1e-5, PC=false])

Create a regular grid on a square with lower left corner LL and upper-right corner UR. Uses Delaunay triangulation internally. If on_torus==true, uses a periodic Delaunay triangulation in both directions. If on_cylinder==true uses a periodic Delaunay triangulatin in x direction only. To avoid degenerate special cases, all nodes are given a random nudge, the strength of which depends on numnodes and nudge_epsilon. If PC==true, returns a piecewise constant grid. Else returns a P1-Lagrange grid.

irregularDelaunayGrid(nodes_in; [on_torus=false, on_cylinder=false, LL, UR, PC=false, ...])

Triangulate the nodes nodes_in and return a GridContext and bdata for them. If on_torus==true, the triangulation is done on a torus. If PC==true, return a mesh with piecewise constant shape-functions, else P1 Lagrange.

regularP2TriangularGrid(numnodes=(25,25), LL=(0.0,0.0), UR=(1.0,1.0), quadrature_order=default_quadrature_order)

Create a regular P2 triangular grid on a rectangle. Does not use Delaunay triangulation internally.

regularP2DelaunayGrid(numnodes=(25,25), LL=(0.0,0.0), UR=(1.0,1.0), quadrature_order=default_quadrature_order)

Create a regular P2 triangular grid with numnodes being the number of (non-interior) nodes in each direction.

regularP2QuadrilateralGrid(numnodes, LL, UR; [quadrature_order, PC=false])

Create a regular P1 quadrilateral grid on a rectangle. If PC==true, use piecewise constant shape functions, otherwise use P1 Lagrange.

regularP2QuadrilateralGrid(numnodes=(25,25), LL=(0.0,0.0), UR=(1.0,1.0), quadrature_order=default_quadrature_order)

Create a regular P2 quadrilateral grid on a rectangle.

regularTetrahedralGrid(numnodes=(10,10,10), LL=(0.0,0.0,0.0), UR=(1.0,1.0,1.0), quadrature_order=default_quadrature_order3D)

Create a regular P1 tetrahedral grid on a Cuboid in 3D. Does not use Delaunay triangulation internally.

regularP2TetrahedralGrid(numnodes=(10,10,10), LL=(0.0,0.0,0.0), UR=(1.0,1.0,1.0), quadrature_order=default_quadrature_order3D)

Create a regular P2 tetrahedral grid on a Cuboid in 3D. Does not use Delaunay triangulation internally.

mutable struct GridContext

Stores everything needed as "context" to be able to work on a FEM grid based on the Ferrite package. Adds a point-locator API which facilitates plotting functions defined on the grid within Julia.

Fields

• grid::FEM.Grid, ip::FEM.Interpolation, ip_geom::FEM.Interpolation, qr::FEM.QuadratureRule; see the Ferrite package.

• loc::PointLocator: object used for point location on the grid.

• node_to_dof::Vector{Int}: lookup table for dof index of a node (for Lagrange elements)

• dof_to_node::Vector{Int}: inverse of nodetodof

• cell_to_dof::Vector{Int}: lookup table for dof index of a cell (for piecewise constant elements)

• dof_to_cell::Vector{Int}: inverse of celltodof

• num_nodes: number of nodes on the grid

• num_cells: number of elements (e.g. triangles,quadrilaterals, ...) on the grid

• n: number of degrees of freedom (== num_nodes for Lagrange Elements, and == num_cells for piecewise constant elements)

• quadrature_points::Vector{Vec{dim,Float64}}: lists of all quadrature points on the grid, in a specific order.

• mass_weights::Vector{Float64}: weighting for stiffness/mass matrices.

• spatialBounds: if available, the corners of a bounding box of a domain. For regular grids, the bounds are tight.

• numberOfPointsInEachDirection: for regular grids, how many (non-interior) nodes make up the regular grid.

• gridType: a string describing what kind of grid this is (e.g. "regular triangular grid")

• no_precompute=false: whether to precompute objects like quadrature points. Only enable this if you know what you are doing.

mutable struct BoundaryData

Represent (a combination of) homogeneous Dirichlet and periodic boundary conditions.

Fields

• dbc_dofs: list of dofs that should have homogeneous Dirichlet boundary conditions. Must be sorted.
• periodic_dofs_from and periodic_dofs_to are both Vector{Int}. The former must be strictly increasing, both must have the same length. periodic_dofs_from[i] is identified with periodic_dofs_to[i]. periodic_dofs_from[i] must be strictly larger than periodic_dofs_to[i]. Multiple dofs can be identified with the same dof. If some dof is identified with another dof and one of them is in dbc_dofs, both points must be in dbc_dofs.

getHomDBCS(ctx, which="all")

Return BoundaryData object corresponding to homogeneous Dirichlet boundary conditions for a set of facesets. which="all" is shorthand for ["left", "right", "top", "bottom"].

undoBCS(ctx, u, bdata)

Given a vector u in dof order with boundary conditions applied, return the corresponding u in dof order without the boundary conditions.

applyBCS(ctx_row, K, bdata_row; [ctx_col, bdata_col, bdata_row, add_vals=true])

Apply the boundary conditions from bdata_row and bdata_col to the sparse matrix K. Only applies boundary conditions accross columns (rows) if bdata_row==nothing (bdata_col==nothing). If add_vals==true (the default), then values in rows that should be cominbed are added. Otherwise, one of the rows is discarded and the values of the other are used.

dof2node(ctx,u)

Interprets u as an array of coefficients ordered in dof order, and reorders them to be in node order.

getDofCoordinates(ctx,dofindex)

Return the coordinates of the node corresponding to the dof with index dofindex.

evaluate_function_from_dofvals(ctx, dofvals, x_in; outside_value=NaN,project_in=false)

Evaluate the function at point x_in with coefficients of dofs given by dofvals (in dof-order). Return outside_value if point is out of bounds. Project the point into the domain if project_in==true. For evaluation at many points, or for many dofvals, the function evaluate_function_from_dofvals_multiple is more efficient.

evaluate_function_from_node_or_cellvals(ctx, vals, x_in; outside_value=0, project_in=false)

Like evaluate_function_from_dofvals, but the coefficients from vals are assumed to be in node order. This is more efficient than evaluate_function_from_dofvals.

evaluate_function_from_node_or_cellvals_multiple(ctx, vals, xin; is_diag=false, kwargs...)

Like evaluate_function_from_dofvals_multiple but uses node- (or cell- if piecewise constant interpolation) ordering for vals, which makes it slightly more efficient. If vals is a diagonal matrix, set is_diag to true for much faster evaluation.

nodal_interpolation(ctx,f)

Perform nodal interpolation of a function. Returns a vector of coefficients in dof order.

getH(ctx)

Return the mesh width of a regular grid.

plot_u(ctx, dof_vals, nx=100, ny=100, LL, UR; bdata=nothing, kwargs...)

Plot the function with coefficients (in dof order, possible boundary conditions in bdata) given by dof_vals on the grid ctx. The domain to be plotted on is given by ctx.spatialBounds. The function is evaluated on a regular nx by ny grid, the resulting plot is a heatmap. Keyword arguments are passed down to plot_u_eulerian, which this function calls internally.

plot_u_eulerian(ctx, dof_vals, inverse_flow_map, LL, UR[, nx, ny];
    euler_to_lagrange_points=nothing,
    only_get_lagrange_points=false,
    z=nothing,
    postprocessor=nothing,
    bdata=nothing,
    throw_errors=false,
    kwargs...)

Plot a heatmap of a function in Eulerian coordinates, i.e., the pushforward of $f$. This is given by $f \circ \Phi^{-1}$, $f$ is a function defined on the grid ctx, represented by coefficients given by dof_vals (with possible boundary conditions given in bdata).

The argument inverse_flow_map is $\Phi^{-1}$.

The resulting plot is on a regular nx x ny grid with lower left corner LL and upper right corner UR. Points that fall outside of the domain represented by ctx are plotted as NaN, which results in transparency.

One can pass values to be plotted directly by providing them in an array in the argument z. postprocessor can modify the values being plotted, return_scalar_field results in these values being returned. See the source code for further details. Additional keyword arguments are passed to Plots.heatmap.

Inverse flow maps are computed in parallel if there are multiple workers.

eulerian_videos(ctx, us, inverse_flow_map_t, t0, tf, nx, ny, nt, LL,UR, num_videos=1;
    extra_kwargs_fun=nothing, ...)

Create num_videos::Int videos in eulerian coordinates, i.e., where the time $t$ is varied, plot $f_i \circ \Phi_t^0$ for $f_1, \dots$.

Arguments

• us(i,t) is a vector of dofs to be plotted at time t for the ith video.
• inverse_flow_map_t(t,x) is $\Phi_t^0(x)$.
• t0, tf are initial and final time.
• LL and UR are the coordinates of the domain's corners.
• nx, ny, and nt determine the number of points in each direction.
• extra_kwargs_fun(i,t) can be used to provide additional keyword arguments to Plots.heatmap().

Additional kwargs are passed on to plot_eulerian_video.

As much as possible is done in parallel.

Returns a Vector of iterables result. Call Plots.animate(result[i]) to get an animation.

eulerian_video(ctx, u, inverse_flow_map_t, t0, tf, nx, ny, nt, LL, UR; extra_kwargs_fun=nothing, ...)

Like eulerian_videos, but u(t) is a vector of dofs, and extra_kwargs_fun(t) gives extra keyword arguments. Returns only one result, on which Plots.animate() can be applied for an animation.

plot_ftle(odefun, p, tspan, LL, UR, nx, ny;
    δ=1e-9, tolerance=1e-4, solver=OrdinaryDiffEq.Tsit5(),
    existing_plot=nothing, flip_y=false, check_inbounds=always_true, pass_on_errors=false)

Make a heatmap of a FTLE field using finite differences. If existing_plot is given a value, plot using heatmap! on top of it. If flip_y is true, then flip the y-coordinate (needed sometimes due to a bug in Plots). Points where check_inbounds(x[1], x[2], p) == false are set to NaN, i.e., plotted transparently. Unless pass_on_errors is set to true, errors from calculating FTLE values are caught and ignored.