Merge branch 'dev' into jc/add_subcell_ops

.github/workflows/Invalidations.yml

@@ -10,7 +10,7 @@ concurrency:
   cancel-in-progress: true
 
 jobs:
-  no_additional_invalidations:
+  evaluate:
     # Only run on PRs to the default branch.
     # In the PR trigger above branches can be specified only explicitly whereas this check should work for master, main, or any other default branch
     if: github.base_ref == github.event.repository.default_branch
@@ -30,11 +30,11 @@ jobs:
     - uses: julia-actions/julia-buildpkg@v1
     - uses: julia-actions/julia-invalidations@v1
       id: invs_default
-    
+
     - name: Report invalidation counts
       run: |
         echo "Invalidations on default branch: ${{ steps.invs_default.outputs.total }} (${{ steps.invs_default.outputs.deps }} via deps)" >> $GITHUB_STEP_SUMMARY
         echo "This branch: ${{ steps.invs_pr.outputs.total }} (${{ steps.invs_pr.outputs.deps }} via deps)" >> $GITHUB_STEP_SUMMARY
     - name: Check if the PR does increase number of invalidations
       if: steps.invs_pr.outputs.total > steps.invs_default.outputs.total
-      run: exit 1
+      run: exit 1

docs/src/RefElemData.md

@@ -6,6 +6,7 @@
 * `Nfaces`: number of faces on a given type of reference element.
 * `fv`: list of vertices defining faces, e.g., `[1,2], [2,3], [3,1]` for a triangle
 * `Fmask`: indices of interpolation nodes which lie on the faces
+* `VDM`: the generalized Vandermonde matrix, a square matrix whose columns are ``V_{ij} = \phi_{j}(x_i}``, where ``\phi_j`` are orthonormal basis functions and ``x_i`` are interpolation points.
 * `rst::NTuple{Dim, ...}`: tuple of vectors of length `N_p`, each of which contains coordinates of degree ``N`` optimized polynomial interpolation points.
 * `rstq::NTuple{Dim, ...}`,`wq`, `Vq`: tuple of volume quadrature points, vector of weights, and quadrature interpolation matrix. Each element of `rstq` and `wq` are vectors of length ``N_q``, and `Vq` is a matrix of size ``N_q \times N_p``.
 * `N_{\rm plot}`: the degree which determines the number of plotting points ``N_{p,{\rm plot}}``.

src/RefElemData.jl

@@ -207,21 +207,32 @@ Polynomial{T}() where {T} = Polynomial(T())
     MultidimensionalQuadrature
 
 A type parameter for `Polynomial` indicating that the quadrature 
-has no specific structure. 
+has no specific structure. Example usage: 
+```julia
+# these are both equivalent
+approximation_type = Polynomial{MultidimensionalQuadrature}() 
+approximation_type = Polynomial(MultidimensionalQuadrature())
+```
 """
 struct MultidimensionalQuadrature end
 
 """
     TensorProductQuadrature{T}
 
 A type parameter to `Polynomial` indicating that the quadrature has a tensor 
-product structure. 
+product structure. Example usage: 
+```julia
+# these are both equivalent
+approximation_type = Polynomial{TensorProductQuadrature}(gauss_quad(0, 0, 1)) 
+approximation_type = Polynomial(TensorProductQuadrature(gauss_quad(0, 0, 1)))
+```
 """
 struct TensorProductQuadrature{T}
     quad_rule_1D::T  # 1D quadrature nodes and weights (rq, wq)
 end
 
 TensorProductQuadrature(args...) = TensorProductQuadrature(args)
+Polynomial{TensorProductQuadrature}(args) = Polynomial(TensorProductQuadrature(args))
 
 """
     TensorProductGaussCollocation

src/RefElemData_polynomial.jl

@@ -1,4 +1,6 @@
-# The following functions determine what default quadrature type to use 
+# The following functions determine what default quadrature type to use for each element.
+# For tensor product elements, we default to TensorProductQuadrature. 
+# For simplices, wedges, and pyramids, we default to MultidimensionalQuadrature
 
 # simplices and pyramids default to multidimensional quadrature
 RefElemData(elem::Union{Line, Tri, Tet, Wedge, Pyr}, 
@@ -15,6 +17,13 @@ RefElemData(elem::Line, approx_type::Polynomial{<:TensorProductQuadrature}, N; k
     RefElemData(elem, Polynomial{MultidimensionalQuadrature}(), N; 
                 quad_rule_vol=approx_type.data.quad_rule_1D, kwargs...)
 
+function RefElemData(elem::Union{Tri, Tet, Wedge, Pyr}, 
+            approx_type::Polynomial{<:TensorProductQuadrature}, 
+            N; kwargs...)
+    error("Tensor product quadrature constructors not yet implemented " * 
+          "for Tri, Tet, Wedge, Pyr elements.")
+end
+
 """
     RefElemData(elem::Line, approximation_type, N;
                 quad_rule_vol = quad_nodes(elem, N+1))

src/explicit_timestep_utils.jl

@@ -6,8 +6,10 @@ Runge Kutta method. Coefficients evolve the residual, solution, and local time,
 
 # Example
 ```julia
-res = rk4a[i]*res + dt*rhs # i = RK stage
-@. u += rk4b[i]*res
+for i in eachindex(rk4a, rk4b)
+    @. res = rk4a[i] * res + dt * rhs # i = RK stage
+    @. u  += rk4b[i] * res
+end
 ```
 """
 function ck45()