StartUpDG 1.0 (#160)

* Make triangle and tet node orderings consistent (#153)

* simplifying test

* fixing vertex orderings

* splitting wedge-pyr MeshData tests out

* fix triangulate tests

* fix wedge mesh ordering

* updating hardcoded values in VTK tests

* update tests

* bump compat for NodesAndModes to 1

* Bump Julia and RecursiveArrayTools compat, remove NamedArrayPartition  (#159)

* remove NamedArrayPartition files

d

* bump compat

* add sparsearrays compat entry

tests seem to be failing on CI without it?

* bump Julia compat

* bump julia CI version

* set lower compat bound for SummationByPartsOperators

* bump doc Julia version

* reexporting NamedArrayPartition

* Refactor `RefElemData` (#157)

* add comments

* comments and slight refactor

* comments and renaming Gauss -> TensorProductGaussCollocation

* update docstrings

* fix comment grammar

* remove unnecessary specialization

* removing outdated comments

* minor reorganization

* removing unnecessary functions

* reorganizing helper functions

* minor comment cleanup

* comments for Kubatko SBP

* comments

* introduce MultidimensionalQuadrature dispatch type

* fix hybrid mesh RefElemData

* fix tests

* simplify SBP RefElemData

* update invalidations

* improve comments

* more descriptive error message

for RefElemData constructors combining Tri/Tet/Wedge/Pyr with TensorProductQuadrature. This should be easy to implement in the future (Stroud)

* Add new cut cell `MeshData`  (#165)

* add some temporary test files

* fix docstrings

* remove cruft

* update plot

made nicer plot of Caratheodory pruning for Christina's proposal

* update file name

Makes it mroe clear this file is meant for plotting

* fixing cut cell demo for v1.0+

* adding Caratheodory pruning

* cleaning up cutcell demo

* improve efficiency slightly

* refactoring functions

* improve comments

* more refactoring

* generalize map_to_interval

* add new version of "generate_sampling_points"

* add dispatch to preserve old version of MeshData

* remove cruft

* add new routines to create a cut cell MeshData with positive weights

* refactoring

* fix construction of cut cell face node indices

previously assumed same number of nodes on all faces

* fix precomputation of operators

* clean up test of SBP property

* add test of SBP property using new MeshData

* add MomentFitting dispatch for old cut cell MeshData

d

* remove outdated todo

* changing wJf to wf internally

* add face node index array

* make face centroid computation more compact

* specialize connect_mesh

* add new MeshData based on subtriangulations

* add Subtriangulation() as a quadrature type

d

* add quadrature type to CutCellMesh meshtype

* test both Subtriangulation and MomentFitting

* add some docs

* allow specifying the target cut cell quadrature degree via keyword arg

* committing scratch testing files before deleting

d

* removing some scratch test files

* add tests for the weak SBP property

* fix an error when num_cartesian_cells = 0

* remove scratchpad files

* Adding hybridized SBP operators (#166)

* move old MomentFitting cut cell code into separate file

* fix default target_degree for volume quadrature

* formatting and modifying default parameters for cut interpolation nodes

* setting default boundary quadrature + renaming variables for clarity

* code cleanup

* storing cut_face_node_indices_by_elem_per_face

* adding hybridized SBP operators on cut cells

* adding cut cell hybridized SBP tests

* update test tolerances

* update hybridized SBP test tolerance

* Improve efficiency of cut-cell `MeshData` (#167)

* diagm -> Diagonal

* use views

* preallocate sorting permutation vector p

* improve efficiency of caratheodory

* remove type instability for PhysicalFrame basis

* bump compat for PathIntersections - should be more type stable

* fix compat for PathIntersections

* add docstring

* release <= 3.4 compat bound on RecursiveArrayTools

* add some NEWS.md updates

* remove SimpleUnpack and usages of @unpack

* Remove Requires.jl (#169)

* remove Requires statements in SummationByPartsOperators.jl ext

* remove requires statements in StartUpDG.jl

* add Plots extension

* remove unnecessary [extras] section (since we have test/Project.toml)

* remove unnecessary compat restriction on SummationByPartsOperators.jl

* add TriangulatePlotExt

* remove Requires.jl

* add to NEWS.md

.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

.github/workflows/ci.yml

@@ -13,7 +13,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.7' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
+          - '1.10' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
           - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
           - 'nightly'
         os:

.github/workflows/documentation.yml

@@ -16,7 +16,7 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1
         with:
-          version: '1.7'
+          version: '1.10'
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy

NEWS.md

@@ -2,6 +2,25 @@
 
 StartUpDG.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1) used in the Julia ecosystem. Recent changes will be documented in this file for human readability.
 
+## Changes when updating to v1.0.0
+
+Most of the major changes are tracked in this [PR](https://github.com/jlchan/StartUpDG.jl/pull/160). Some descriptions of other changes are listed below. 
+
+#### Added
+
+* Generation of cut-cell meshes using `Subtriangulation` quadrature by default, which ensures positive quadrature weights. The old behavior is retained by specifying a `MomentFitting` quadrature type. 
+* Added `subcell_limiting_operators`, which constructs multi-dimensional versions of subcell operators used in [Lin, Chan (2024)](https://doi.org/10.1016/j.jcp.2023.112677). These subcell operators are constructed from sparse operators returned by `sparse_low_order_SBP_operators(rd::RefElemData)`.
+
+#### Changed
+
+* `NamedArrayPartition` was moved to RecursiveArrayTools.jl. 
+* The required Julia version was increased to v1.10. This was to make StartUpDG.jl compatibility with RecursiveArrayTools.jl v3.4+ (see above).
+* Removed SimpleUnpack.jl as a dependency. Loading StartUpDG.jl will no longer reexport `@unpack`, since destructuring via `(; propertyname) = x` is supported natively in Julia 1.7 and up.
+* Updated to NodesAndModes v1.0+, which changed the ordering of triangle nodes to make them consistent with tet node ordering. 
+* We introduced a `MultidimensionalQuadrature` type. All `Polynomial` approximation types now utilize either `MultidimensionalQuadrature` or `TensorProductQuadrature` as a type parameter. The previous type parameter `DefaultPolynomialType` is now simply used to determine the default quadrature type parameter. Note that this is internal behavior and should not impact standard usage of StartUpDG.jl.
+* Removed Requires.jl in favor of [package extensions](https://pkgdocs.julialang.org/v1/creating-packages/#Conditional-loading-of-code-in-packages-(Extensions)) for Plots.jl and SummationByPartsOperators.jl dependencies. 
+
+
 ## Changes when updating to v0.17
 
 #### Added

Project.toml

@@ -1,7 +1,7 @@
 name = "StartUpDG"
 uuid = "472ebc20-7c99-4d4b-9470-8fde4e9faa0f"
 authors = ["Jesse Chan", "Yimin Lin"]
-version = "0.17.7"
+version = "1.0.0"
 
 [deps]
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
@@ -15,38 +15,34 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
-Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
-SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [weakdeps]
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240"
 
 [extensions]
 StartUpDGSummationByPartsOperatorsExt = "SummationByPartsOperators"
+TriangulatePlotsExt = "Plots"
 
 [compat]
 ConstructionBase = "1"
 FillArrays = "0.13, 1"
 HDF5 = "0.16, 0.17"
 Kronecker = "0.5"
-NodesAndModes = "0.9"
-PathIntersections = "0.1"
+NodesAndModes = "1"
+PathIntersections = "0.1, 0.2"
 RecipesBase = "1"
-RecursiveArrayTools = "2"
+RecursiveArrayTools = "3"
 Reexport = "1"
-Requires = "1"
 Setfield = "1"
-SimpleUnPack = "1"
+SparseArrays = "1"
 StaticArrays = "1"
 SummationByPartsOperators = "0.5"
 Triangulate = "2"
 WriteVTK = "1"
-julia = "1.7"
-
-[extras]
-SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240"
+julia = "1.10"

docs/src/more_meshes.md

@@ -72,10 +72,15 @@ circle = PresetGeometries.Circle(R=0.33, x0=0, y0=0)
 cells_per_dimension_x, cells_per_dimension_y = 4, 4
 
 rd = RefElemData(Quad(), N=3)
-md = MeshData(rd, (circle, ), cells_per_dimension_x, cells_per_dimension_y; precompute_operators=true)
+md = MeshData(rd, (circle, ), cells_per_dimension_x, cells_per_dimension_y, Subtriangulation(); precompute_operators=true)
 ```
-The interpolation points on cut cells `md.x.cut` are determined from sampled points and a pivoted QR decomposition. The quadrature points on cut cells `md.xq.cut` are determined similarly. However, the cut-cell quadrature weights `md.wJq.cut` are not currently positive. The optional keyword argument `precompute_operators` specifies 
-whether to precompute differentiation, face interpolation, mass, and lifting matrices for each cut cell. If 
+Here, the final argument `quadrature_type = Subtriangulation()` determines how the quadrature on cut cells is determined. For `Subtriangulation()`, the quadrature on cut cells is constructed from a curved isoparametric subtriangulation of the cut cell. The number of quadrature points on a cut cell is then reduced (while preserving positivity) using Caratheodory pruning. If not specified, the `quadrature_type` argument defaults to `Subtriangulation()`. 
+
+Quadrature rules can also be constructed by specifying `quadrature_type = MomentFitting()`. The quadrature points on cut cells `md.xq.cut` are determined from sampling and a pivoted QR decomposition. This is not recommended, as it can be both slower, and the cut-cell quadrature weights `md.wJq.cut` are not guaranteed to be positive. 
+
+The interpolation points on cut cells `md.x.cut` are determined from sampled points and a pivoted QR decomposition. 
+
+The optional keyword argument `precompute_operators` specifies whether to precompute differentiation, face interpolation, mass, and lifting matrices for each cut cell. If 
 `precompute_operators=true`, these are stored in `md.mesh_type.cut_cell_operators`. 
 
 As with hybrid meshes, the nodal coordinates `md.x`, `md.y` are `NamedArrayPartition`s with `cartesian` and `cut` fields. For example, `md.x.cartesian` and `md.x.cut` are the x-coordinates of the Cartesian and cut cells, respectively. Likewise, `md.mapP` indexes linearly into the array of face coordinates and specifies exterior node indices. For example, we can interpolate a function to face nodes and compute exterior values via the following code:

ext/StartUpDGSummationByPartsOperatorsExt.jl

@@ -6,26 +6,14 @@ using SparseArrays: sparse, droptol!, spzeros
 using StartUpDG
 
 # Required for visualization code
-if isdefined(Base, :get_extension)
-    using SummationByPartsOperators:
-        SummationByPartsOperators,
-        DerivativeOperator,
-        grid,
-        AbstractDerivativeOperator,
-        AbstractNonperiodicDerivativeOperator,
-        PeriodicDerivativeOperator,
-        AbstractPeriodicDerivativeOperator
-else
-    # Until Julia v1.9 is the minimum required version for Trixi.jl, we still support Requires.jl
-    using ..SummationByPartsOperators
-    using ..SummationByPartsOperators:
-        AbstractDerivativeOperator,
-        AbstractPeriodicDerivativeOperator,
-        AbstractNonperiodicDerivativeOperator,
-        DerivativeOperator,
-        PeriodicDerivativeOperator,
-        grid
-end
+using SummationByPartsOperators:
+    SummationByPartsOperators,
+    DerivativeOperator,
+    grid,
+    AbstractDerivativeOperator,
+    AbstractNonperiodicDerivativeOperator,
+    PeriodicDerivativeOperator,
+    AbstractPeriodicDerivativeOperator
 
 function construct_1d_operators(D::AbstractDerivativeOperator, tol)
     nodes_1d = collect(grid(D))

src/TriangulatePlots.jl → ext/TriangulatePlotsExt.jl

@@ -1,8 +1,8 @@
-module TriangulatePlots
+module TriangulatePlotsExt
 
 using StartUpDG: BoundaryTagPlotter, RecipesBase
 
-using ..Plots: Plots
+using Plots: Plots
 
 RecipesBase.@recipe function f(m::BoundaryTagPlotter)
   triout = m.triout

src/RefElemData.jl

@@ -76,7 +76,7 @@ function Base.propertynames(x::RefElemData{3}, private::Bool = false)
 end
 
 # convenience unpacking routines
-function Base.getproperty(x::RefElemData{Dim, ElementType, ApproxType}, s::Symbol) where {Dim, ElementType, ApproxType}
+function Base.getproperty(x::RefElemData, s::Symbol) 
     if s==:r
         return getfield(x, :rst)[1]
     elseif s==:s
@@ -133,16 +133,18 @@ function Base.getproperty(x::RefElemData{Dim, ElementType, ApproxType}, s::Symbo
 end
 
 """
-    function RefElemData(elem; N, kwargs...)
-    function RefElemData(elem, approx_type; N, kwargs...)
+    RefElemData(elem; N, kwargs...)
+    RefElemData(elem, approx_type; N, kwargs...)
 
 Keyword argument constructor for RefElemData (to "label" `N` via `rd = RefElemData(Line(), N=3)`)
 """
 RefElemData(elem; N, kwargs...) = RefElemData(elem, N; kwargs...)
-RefElemData(elem, approx_type; N, kwargs...) = RefElemData(elem, approx_type, N; kwargs...)
+RefElemData(elem, approx_type; N, kwargs...) = 
+    RefElemData(elem, approx_type, N; kwargs...)
 
 # default to Polynomial-type RefElemData
-RefElemData(elem, N::Int; kwargs...) = RefElemData(elem, Polynomial(), N; kwargs...)
+RefElemData(elem, N::Int; kwargs...) = 
+    RefElemData(elem, Polynomial(), N; kwargs...)
 
 
 @inline Base.ndims(::Line) = 1
@@ -170,12 +172,12 @@ RefElemData(elem, N::Int; kwargs...) = RefElemData(elem, Polynomial(), N; kwargs
 
 # Wedges have different types of faces depending on the face. 
 # We define the first three faces to be quadrilaterals and the 
-# last two faces are triangles.
+# last two faces to be triangles.
 @inline face_type(::Wedge, id) = (id <= 3) ? Quad() : Tri()
 
 # Pyramids have different types of faces depending on the face. 
 # We define the first four faces to be triangles and the 
-# last face to be a quadrilateral. 
+# last face to be the quadrilateral face. 
 @inline face_type(::Pyr, id) = (id <= 4) ? Tri() : Quad()
 
 # ====================================================
@@ -197,20 +199,49 @@ end
 struct DefaultPolynomialType end
 Polynomial() = Polynomial{DefaultPolynomialType}(DefaultPolynomialType())
 
+# this constructor enables us to construct a `Polynomial` type via
+# Polynomial{TensorProductQuadrature}() or Polynomial{MultidimensionalQuadrature}().  
+Polynomial{T}() where {T} = Polynomial(T())
+
+"""
+    MultidimensionalQuadrature
+
+A type parameter for `Polynomial` indicating that the quadrature 
+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 
+A type parameter to `Polynomial` indicating that the quadrature has a tensor 
+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(r1D, w1D) = TensorProductQuadrature((r1D, w1D))
+TensorProductQuadrature(args...) = TensorProductQuadrature(args)
+Polynomial{TensorProductQuadrature}(args) = Polynomial(TensorProductQuadrature(args))
 
-# Polynomial{Gauss} type indicates (N+1)-point Gauss quadrature on tensor product elements
-struct Gauss end 
-Polynomial{Gauss}() = Polynomial(Gauss())
+"""
+    TensorProductGaussCollocation
+
+Polynomial{TensorProductGaussCollocation} type indicates a tensor product 
+# (N+1)-point Gauss quadrature on tensor product elements. 
+"""
+struct TensorProductGaussCollocation end 
+const Gauss = TensorProductGaussCollocation
 
 # ========= SBP approximation types ============
 
@@ -223,7 +254,7 @@ struct TensorProductLobatto end
 struct Hicken end 
 struct Kubatko{FaceNodeType} end
 
-# face node types for Kubatko
+# face node types for Kubatko nodes
 struct LegendreFaceNodes end
 struct LobattoFaceNodes end
 
@@ -268,6 +299,7 @@ _short_typeof(x) = typeof(x)
 _short_typeof(approx_type::Wedge) = "Wedge"
 _short_typeof(approx_type::Pyr) = "Pyr"
 
-_short_typeof(approx_type::Polynomial{<:DefaultPolynomialType}) = "Polynomial"
-_short_typeof(approx_type::Polynomial{<:Gauss}) = "Polynomial{Gauss}"
+# _short_typeof(approx_type::Polynomial{<:DefaultPolynomialType}) = "Polynomial"
+_short_typeof(approx_type::Polynomial{<:MultidimensionalQuadrature}) = "Polynomial"
+_short_typeof(approx_type::Polynomial{<:TensorProductGaussCollocation}) = "Polynomial{Gauss}"
 _short_typeof(approx_type::Polynomial{<:TensorProductQuadrature}) = "Polynomial{TensorProductQuadrature}"