diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml index ba81f83e..cfcaa3af 100644 --- a/.github/workflows/Invalidations.yml +++ b/.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 \ No newline at end of file diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 4690dd6e..0ce96dc6 100644 --- a/.github/workflows/ci.yml +++ b/.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: diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 8f552204..56242b26 100644 --- a/.github/workflows/documentation.yml +++ b/.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 diff --git a/NEWS.md b/NEWS.md index d6020d5f..44a2416f 100644 --- a/NEWS.md +++ b/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 diff --git a/Project.toml b/Project.toml index 89ae6440..b8e65a56 100644 --- a/Project.toml +++ b/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" \ No newline at end of file diff --git a/docs/src/more_meshes.md b/docs/src/more_meshes.md index 8190b26b..5f1e3bda 100644 --- a/docs/src/more_meshes.md +++ b/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: diff --git a/ext/StartUpDGSummationByPartsOperatorsExt.jl b/ext/StartUpDGSummationByPartsOperatorsExt.jl index e036cda8..f0db42c2 100644 --- a/ext/StartUpDGSummationByPartsOperatorsExt.jl +++ b/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)) diff --git a/src/TriangulatePlots.jl b/ext/TriangulatePlotsExt.jl similarity index 96% rename from src/TriangulatePlots.jl rename to ext/TriangulatePlotsExt.jl index 1fcc6669..7c44d4b9 100644 --- a/src/TriangulatePlots.jl +++ b/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 diff --git a/src/RefElemData.jl b/src/RefElemData.jl index ac14339b..ed1b9084 100644 --- a/src/RefElemData.jl +++ b/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}" diff --git a/src/RefElemData_SBP.jl b/src/RefElemData_SBP.jl index dbc780c6..14fba196 100644 --- a/src/RefElemData_SBP.jl +++ b/src/RefElemData_SBP.jl @@ -1,8 +1,8 @@ """ - function RefElemData(elementType::Line, approxType::SBP, N) - function RefElemData(elementType::Quad, approxType::SBP, N) - function RefElemData(elementType::Hex, approxType::SBP, N) - function RefElemData(elementType::Tri, approxType::SBP, N) + RefElemData(elementType::Line, approxType::SBP, N) + RefElemData(elementType::Quad, approxType::SBP, N) + RefElemData(elementType::Hex, approxType::SBP, N) + RefElemData(elementType::Tri, approxType::SBP, N) SBP reference element data for `Quad()`, `Hex()`, and `Tri()` elements. @@ -20,38 +20,11 @@ function RefElemData(elementType::Line, approxType::SBP{TensorProductLobatto}, N return _convert_RefElemData_fields_to_SBP(rd, approxType) end -function RefElemData(elementType::Quad, approxType::SBP{TensorProductLobatto}, N; tol = 100*eps(), kwargs...) +function RefElemData(elementType::Union{Quad, Hex}, approxType::SBP{TensorProductLobatto}, N; tol = 100*eps(), kwargs...) - # make 2D SBP nodes/weights - r1D, w1D = gauss_lobatto_quad(0, 0, N) - sq, rq = vec.(NodesAndModes.meshgrid(r1D)) # this is to match ordering of nrstJ - wr, ws = vec.(NodesAndModes.meshgrid(w1D)) - wq = wr .* ws - quad_rule_vol = (rq, sq, wq) - quad_rule_face = (r1D, w1D) - - rd = RefElemData(elementType, N; quad_rule_vol = quad_rule_vol, quad_rule_face = quad_rule_face, kwargs...) - - rd = @set rd.Vf = droptol!(sparse(rd.Vf), tol) - rd = @set rd.LIFT = Diagonal(rd.wq) \ (rd.Vf' * Diagonal(rd.wf)) # TODO: make this more efficient with LinearMaps? - - return _convert_RefElemData_fields_to_SBP(rd, approxType) -end - -function RefElemData(elementType::Hex, approxType::SBP{TensorProductLobatto}, N; tol = 100*eps(), kwargs...) - - # make 2D SBP nodes/weights - r1D, w1D = gauss_lobatto_quad(0, 0, N) - rf, sf = vec.(NodesAndModes.meshgrid(r1D, r1D)) - wr, ws = vec.(NodesAndModes.meshgrid(w1D, w1D)) - wf = wr .* ws - sq, rq, tq = vec.(NodesAndModes.meshgrid(r1D, r1D, r1D)) # this is to match ordering of nrstJ - wr, ws, wt = vec.(NodesAndModes.meshgrid(w1D, w1D, w1D)) - wq = wr .* ws .* wt - quad_rule_vol = (rq, sq, tq, wq) - quad_rule_face = (rf, sf, wf) - - rd = RefElemData(elementType, N; quad_rule_vol = quad_rule_vol, quad_rule_face = quad_rule_face, kwargs...) + rd = RefElemData(elementType, + Polynomial(TensorProductQuadrature(gauss_lobatto_quad(0, 0, N))), + N; kwargs...) rd = @set rd.Vf = droptol!(sparse(rd.Vf), tol) rd = @set rd.LIFT = Diagonal(rd.wq) \ (rd.Vf' * Diagonal(rd.wf)) # TODO: make this more efficient with LinearMaps? @@ -221,158 +194,5 @@ function hybridized_SBP_operators(rd) return Qrsth, VhP, Ph, Vh end -# default to doing nothing -map_nodes_to_symmetric_element(element_type, rst...) = rst - -# for triangles, map to an equilateral triangle -function map_nodes_to_symmetric_element(::Tri, r, s) - # biunit right triangular vertices - v1, v2, v3 = SVector{2}.(zip(nodes(Tri(), 1)...)) - - denom = (v2[2] - v3[2]) * (v1[1] - v3[1]) + (v3[1] - v2[1]) * (v1[2] - v3[2]) - L1 = @. ((v2[2] - v3[2]) * (r - v3[1]) + (v3[1] - v2[1]) * (s - v3[2])) / denom - L2 = @. ((v3[2] - v1[2]) * (r - v3[1]) + (v1[1] - v3[1]) * (s - v3[2])) / denom - L3 = @. 1 - L1 - L2 - - # equilateral vertices - v1 = SVector{2}(2 * [-.5, -sqrt(3) / 6]) - v2 = SVector{2}(2 * [.5, -sqrt(3)/6]) - v3 = SVector{2}(2 * [0, sqrt(3)/3]) - x = @. v1[1] * L1 + v2[1] * L2 + v3[1] * L3 - y = @. v1[2] * L1 + v2[2] * L2 + v3[2] * L3 - return x, y -end - - -""" - function sparse_low_order_SBP_operators(rd) - -Constructs sparse low order SBP operators given a `RefElemData`. -Returns operators `Qrst..., E ≈ Vf * Pq` that satisfy a generalized -summation-by-parts (GSBP) property: - - `Q_i + Q_i^T = E' * B_i * E` -""" -function sparse_low_order_SBP_operators(rd::RefElemData{NDIMS}) where {NDIMS} - (; Pq, Vf, rstq, wf, nrstJ) = rd - - # if element is a simplex, convert nodes to an equilateral triangle for symmetry - rstq = map_nodes_to_symmetric_element(rd.element_type, rstq...) - - # build distance and adjacency matrices - D = [norm(SVector{NDIMS}(getindex.(rstq, i)) - SVector{NDIMS}(getindex.(rstq, j))) - for i in eachindex(first(rstq)), j in eachindex(first(rstq))] - A = zeros(Int, length(first(rstq)), length(first(rstq))) - for i in axes(D, 1) - # heuristic cutoff - we take NDIMS + 1 neighbors, but the smallest distance = 0, - # so we need to access the first NDIMS + 2 sorted distance entries. - dist = sort(view(D, i, :))[NDIMS + 2] - for j in findall(view(D, i, :) .< 1.01 * dist) - A[i, j] = 1 - end - end - A = (A + A' .> 0) # symmetrize - L = Diagonal(vec(sum(A, dims=2))) - A - sorted_eigvals = sort(eigvals(L)) - - # check that there's only one zero null vector - @assert sorted_eigvals[2] > 100 * eps() - - E_dense = Vf * Pq - E = zeros(size(E_dense)) - for i in axes(E, 1) - # find all j such that E[i,j] ≥ 0.5, e.g., points which positively contribute to at least half of the - # interpolation. These seem to be associated with volume points "j" that are close to face point "i". - ids = findall(E_dense[i, :] .>= 0.5) - E[i, ids] .= E_dense[i, ids] ./ sum(E_dense[i, ids]) # normalize so sum(E, dims=2) = [1, 1, ...] still. - end - Brst = (nJ -> diagm(wf .* nJ)).(nrstJ) - - # compute potential - e = ones(size(L, 2)) - right_hand_sides = map(B -> 0.5 * sum(E' * B, dims=2), Brst) - psi_augmented = map(b -> [L e; e' 0] \ [b; 0], right_hand_sides) - psi = map(x -> x[1:end-1], psi_augmented) - - # construct sparse skew part - function construct_skew_matrix_from_potential(ψ) - S = zeros(length(ψ), length(ψ)) - for i in axes(S, 1), j in axes(S, 2) - if A[i,j] > 0 - S[i,j] = ψ[j] - ψ[i] - end - end - return S - end - - Srst = construct_skew_matrix_from_potential.(psi) - Qrst = map((S, B) -> S + 0.5 * E' * B * E, Srst, Brst) - return sparse.(Qrst), sparse(E) -end - -function sparse_low_order_SBP_1D_operators(rd::RefElemData) - E = zeros(2, rd.N+1) - E[1, 1] = 1 - E[2, end] = 1 - - # create volume operators - Q = diagm(1 => ones(rd.N), -1 => -ones(rd.N)) - Q[1,1] = -1 - Q[end,end] = 1 - Q = 0.5 * Q - - return (sparse(Q),), sparse(E) -end - -sparse_low_order_SBP_operators(rd::RefElemData{1, Line, <:Union{<:SBP, <:Polynomial{Gauss}}}) = - sparse_low_order_SBP_1D_operators(rd) - -function diagonal_1D_mass_matrix(N, ::SBP) - _, w1D = gauss_lobatto_quad(0, 0, N) - return Diagonal(w1D) -end - -function diagonal_1D_mass_matrix(N, ::Polynomial{Gauss}) - _, w1D = gauss_quad(0, 0, N) - return Diagonal(w1D) -end - -function sparse_low_order_SBP_operators(rd::RefElemData{2, Quad, <:Union{<:SBP, <:Polynomial{Gauss}}}) - (Q1D,), E1D = sparse_low_order_SBP_1D_operators(rd) - - # permute face node ordering for the first 2 faces - ids = reshape(1:(rd.N+1) * 2, :, 2) - Er = zeros((rd.N+1) * 2, rd.Np) - Er[vec(ids'), :] .= kron(I(rd.N+1), E1D) - Es = kron(E1D, I(rd.N+1)) - E = vcat(Er, Es) - - M1D = diagonal_1D_mass_matrix(rd.N, rd.approximation_type) - Qr = kron(M1D, Q1D) - Qs = kron(Q1D, M1D) - - return sparse.((Qr, Qs)), sparse(E) -end - -function sparse_low_order_SBP_operators(rd::RefElemData{3, Hex, <:Union{<:SBP, <:Polynomial{Gauss}}}) - (Q1D,), E1D = sparse_low_order_SBP_1D_operators(rd) - - # permute face nodes for the faces in the ±r and ±s directions - ids = reshape(1:(rd.N+1)^2 * 2, rd.N+1, :, 2) - Er, Es, Et = ntuple(_ -> zeros((rd.N+1)^2 * 2, rd.Np), 3) - Er[vec(permutedims(ids, [2, 3, 1])), :] .= kron(I(rd.N+1), E1D, I(rd.N+1)) - Es[vec(permutedims(ids, [3, 2, 1])), :] .= kron(I(rd.N+1), I(rd.N+1), E1D) - Et .= kron(E1D, I(rd.N+1), I(rd.N+1)) - - # create boundary extraction matrix - E = vcat(Er, Es, Et) - - M1D = diagonal_1D_mass_matrix(rd.N, rd.approximation_type) - Qr = kron(M1D, Q1D, M1D) - Qs = kron(M1D, M1D, Q1D) - Qt = kron(Q1D, M1D, M1D) - - return sparse.((Qr, Qs, Qt)), sparse(E) -end diff --git a/src/RefElemData_TensorProductWedge.jl b/src/RefElemData_TensorProductWedge.jl index 5fdbe4e1..ac39f3ac 100644 --- a/src/RefElemData_TensorProductWedge.jl +++ b/src/RefElemData_TensorProductWedge.jl @@ -95,7 +95,8 @@ function RefElemData(elem::Wedge, approximation_type::TensorProductWedge; kwargs wt, wrs = _wedge_tensor_product(line.wq, tri.wq) wq = wt .* wrs - # `line.Vq` is a `UniformScaling` type for `RefElemData` built from SummationByPartsOperators.jl in Trixi.jl + # `line.Vq` is a `UniformScaling` type for `RefElemData` built + # from SummationByPartsOperators.jl Vq = line.Vq isa UniformScaling ? kron(I(num_line_nodes), tri.Vq) : kron(line.Vq, tri.Vq) M = Vq' * diagm(wq) * Vq Pq = M \ (Vq' * diagm(wq)) @@ -104,7 +105,7 @@ function RefElemData(elem::Wedge, approximation_type::TensorProductWedge; kwargs # tensor product plotting nodes tp, rp = _wedge_tensor_product(line.rp, tri.rp) _, sp = _wedge_tensor_product(line.rp, tri.sp) - # `line.Vp` is a `UniformScaling` type for `RefElemData` from SummationByPartsOperators.jl in Trixi.jl + # `line.Vp` is a `UniformScaling` type for `RefElemData` from SummationByPartsOperators.jl Vp = line.Vp isa UniformScaling ? kron(I(num_line_nodes), tri.Vp) : kron(line.Vp, tri.Vp) # set the polynomial degree as the tuple of the line and triangle degree for now diff --git a/src/RefElemData_polynomial.jl b/src/RefElemData_polynomial.jl index f8c6fe70..45564982 100644 --- a/src/RefElemData_polynomial.jl +++ b/src/RefElemData_polynomial.jl @@ -1,94 +1,65 @@ -function init_face_data(elem::Tri; quad_rule_face = gauss_quad(0,0,N)) - r1D, w1D = quad_rule_face - e = ones(size(r1D)) - z = zeros(size(r1D)) - rf, sf = map_face_nodes(elem, r1D) - wf = vec(repeat(w1D, 3, 1)); - nrJ = [z; e; -e] - nsJ = [-e; e; z] - return rf,sf,wf,nrJ,nsJ +# 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}, + approx_type::Polynomial{DefaultPolynomialType}, + N; kwargs...) = + RefElemData(elem, Polynomial{MultidimensionalQuadrature}(), N; kwargs...) + +# on quad and hex elements, default to a tensor product quadrature +RefElemData(elem::Union{Quad, Hex}, approximation_type::Polynomial{DefaultPolynomialType}, N; kwargs...) = + RefElemData(elem, Polynomial(TensorProductQuadrature(gauss_quad(0, 0, N+1))), N; kwargs...) + +# special case: for lines, tensor product and multidimensional quadrature are the same +RefElemData(elem::Line, approx_type::Polynomial{<:TensorProductQuadrature}, N; kwargs...) = + 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 -function init_face_data(elem::Quad; quad_rule_face=gauss_quad(0, 0, N)) - Nfaces = 4 - r1D, w1D = quad_rule_face - e = ones(size(r1D)) - z = zeros(size(r1D)) - rf, sf = map_face_nodes(elem, r1D) - wf = vec(repeat(w1D, Nfaces, 1)) - nrJ = [-e; e; z; z] - nsJ = [z; z; -e; e] - - return rf, sf, wf, nrJ, nsJ -end - -function init_face_data(elem::Hex; quad_rule_face=quad_nodes(Quad(), N)) - rquad, squad, wquad = quad_rule_face - e = ones(size(rquad)) - zz = zeros(size(rquad)) - rf, sf, tf = map_face_nodes(elem, rquad, squad) - Nfaces = 6 - wf = vec(repeat(wquad, Nfaces, 1)); - nrJ = [-e; e; zz; zz; zz; zz] - nsJ = [zz; zz; -e; e; zz; zz] - ntJ = [zz; zz; zz; zz; -e; e] - return rf, sf, tf, wf, nrJ, nsJ, ntJ -end - -function init_face_data(elem::Tet; quad_rule_face=quad_nodes(Tri(), N)) - rquad, squad, wquad = quad_rule_face - e = ones(size(rquad)) - zz = zeros(size(rquad)) - rf, sf, tf = map_face_nodes(elem, rquad, squad) - Nfaces = 4 - wf = vec(repeat(wquad, Nfaces, 1)); - nrJ = [zz; e; -e; zz] - nsJ = [-e; e; zz; zz] - ntJ = [zz; e; zz; -e] - return rf, sf, tf, wf, nrJ, nsJ, ntJ -end - - """ - RefElemData(elem::Line, N; + RefElemData(elem::Line, approximation_type, N; quad_rule_vol = quad_nodes(elem, N+1)) - RefElemData(elem::Union{Tri, Quad}, N; - quad_rule_vol = quad_nodes(elem, N), - quad_rule_face = gauss_quad(0, 0, N)) - RefElemData(elem::Union{Hex, Tet}, N; - quad_rule_vol = quad_nodes(elem, N), - quad_rule_face = quad_nodes(Quad(), N)) - RefElemData(elem; N, kwargs...) # version with keyword args + RefElemData(elem, approximation_type, N; + quad_rule_vol = quad_nodes(elem, N), + quad_rule_face = quad_nodes(face_type(elem), N)) Constructor for `RefElemData` for different element types. """ -function RefElemData(elem::Line, approx_type::Polynomial{DefaultPolynomialType}, N; - quad_rule_vol=quad_nodes(elem, N+1), - Nplot=10) +function RefElemData(elem::Line, + approx_type::Polynomial{MultidimensionalQuadrature}, + N; quad_rule_vol=quad_nodes(elem, N+1), Nplot=10) fv = face_vertices(elem) - # Construct matrices on reference elements + # reference element nodes r = nodes(elem, N) Fmask = [1 N+1] + + # compute operators VDM = vandermonde(elem, N, r) Dr = grad_vandermonde(elem, N, r)/VDM + V1 = vandermonde(elem, 1, r) / vandermonde(elem, 1, nodes(elem, 1)) - V1 = vandermonde(elem, 1, r) / vandermonde(elem, 1, [-1; 1]) - + # quadrature operators rq, wq = quad_rule_vol + rf, wf = [-1.0; 1.0], [1.0; 1.0] + nrJ = [-1.0; 1.0] Vq = vandermonde(elem, N, rq) / VDM M = Vq' * diagm(wq) * Vq Pq = M \ (Vq' * diagm(wq)) - - rf = [-1.0; 1.0] - nrJ = [-1.0; 1.0] - wf = [1.0; 1.0] Vf = vandermonde(elem, N, rf) / VDM LIFT = M \ (Vf') # lift matrix # plotting nodes - rp = equi_nodes(elem, Nplot) + rp = equi_nodes(elem, Nplot) Vp = vandermonde(elem, N, rp) / VDM return RefElemData(elem, approx_type, N, fv, V1, @@ -99,8 +70,9 @@ function RefElemData(elem::Line, approx_type::Polynomial{DefaultPolynomialType}, M, Pq, tuple(Dr), LIFT) end -function RefElemData(elem::Union{Tri, Quad}, - approx_type::Polynomial{DefaultPolynomialType}, N; + +function RefElemData(elem::Union{Tri, Quad}, + approx_type::Polynomial{MultidimensionalQuadrature}, N; quad_rule_vol=quad_nodes(elem, N), quad_rule_face=quad_nodes(face_type(elem), N), Nplot=10) @@ -111,26 +83,27 @@ function RefElemData(elem::Union{Tri, Quad}, r, s = nodes(elem, N) Fmask = hcat(find_face_nodes(elem, r, s)...) - VDM, Vr, Vs = basis(elem, N, r, s) - Dr = Vr / VDM - Ds = Vs / VDM - # low order interpolation nodes - r1, s1 = nodes(elem, 1) + r1, s1 = nodes(elem, 1) V1 = vandermonde(elem, 1, r, s) / vandermonde(elem, 1, r1, s1) - rf, sf, wf, nrJ, nsJ = init_face_data(elem, quad_rule_face = quad_rule_face) + # differentiation operators + VDM, Vr, Vs = basis(elem, N, r, s) + Dr = Vr / VDM + Ds = Vs / VDM + # quadrature nodes rq, sq, wq = quad_rule_vol + rf, sf, wf, nrJ, nsJ = init_face_data(elem; quad_rule_face) + + # quadrature-based operators Vq = vandermonde(elem, N, rq, sq) / VDM + Vf = vandermonde(elem, N, rf, sf) / VDM # interpolates from nodes to face nodes M = Vq' * diagm(wq) * Vq Pq = M \ (Vq' * diagm(wq)) - - Vf = vandermonde(elem, N, rf, sf) / VDM # interpolates from nodes to face nodes LIFT = M \ (Vf' * diagm(wf)) # lift matrix used in rhs evaluation - # plotting nodes - rp, sp = equi_nodes(elem, Nplot) + rp, sp = equi_nodes(elem, Nplot) # plotting nodes Vp = vandermonde(elem, N, rp, sp) / VDM return RefElemData(elem, approx_type, N, fv, V1, @@ -141,13 +114,14 @@ function RefElemData(elem::Union{Tri, Quad}, M, Pq, (Dr, Ds), LIFT) end -function RefElemData(elem::Union{Hex, Tet}, approx_type::Polynomial{DefaultPolynomialType}, N; +function RefElemData(elem::Union{Tet, Hex}, + approx_type::Polynomial{MultidimensionalQuadrature}, N; quad_rule_vol=quad_nodes(elem, N), quad_rule_face=quad_nodes(face_type(elem), N), Nplot=10) if elem isa Hex && N > 4 - @warn "Since N > 4, we suggest using `RefElemData(Hex(), TensorProductQuadrature(gauss_quad(0, 0, $N+1)), $N)`, " * + @warn "Since N > 4, we suggest using `RefElemData(Hex(), Polynomial(TensorProductQuadrature(gauss_quad(0, 0, $N+1))), $N)`, " * "which is more efficient." end @@ -164,14 +138,13 @@ function RefElemData(elem::Union{Hex, Tet}, approx_type::Polynomial{DefaultPolyn V1 = vandermonde(elem, 1, r, s, t) / vandermonde(elem, 1, r1, s1, t1) # Nodes on faces, and face node coordinate - rf, sf, tf, wf, nrJ, nsJ, ntJ = init_face_data(elem, quad_rule_face = quad_rule_face) - - # quadrature nodes - build from 1D nodes. rq, sq, tq, wq = quad_rule_vol + rf, sf, tf, wf, nrJ, nsJ, ntJ = init_face_data(elem; quad_rule_face) + + # quadrature operators Vq = vandermonde(elem, N, rq, sq, tq) / VDM M = Vq' * diagm(wq) * Vq Pq = M \ (Vq' * diagm(wq)) - Vf = vandermonde(elem, N, rf, sf, tf) / VDM LIFT = M \ (Vf' * diagm(wf)) @@ -187,74 +160,6 @@ function RefElemData(elem::Union{Hex, Tet}, approx_type::Polynomial{DefaultPolyn M, Pq, (Dr, Ds, Dt), LIFT) end -""" - RefElemData(elem::Hex, ::TensorProductQuadrature, N; Nplot = 10) - RefElemData(elem::Hex, approximation_type::Polynomial{<:TensorProductQuadrature}, N; Nplot = 10) - -Constructor for hexahedral `RefElemData` where the quadrature is assumed to have tensor product structure. -""" -RefElemData(elem::Hex, approximation_parameter::TensorProductQuadrature, N; Nplot = 10) = - RefElemData(elem, Polynomial(approximation_parameter), N; Nplot) - -function RefElemData(elem::Hex, - approximation_type::Polynomial{<:TensorProductQuadrature}, N; - Nplot = 10) - - fv = face_vertices(elem) - - # Construct matrices on reference elements - r, s, t = nodes(elem, N) - Fmask = hcat(find_face_nodes(elem, r, s, t)...) - - # construct 1D operator for faster Kronecker solves - r1D = nodes(Line(), N) - rq1D, wq1D = approximation_type.data.quad_rule_1D - VDM_1D = vandermonde(Line(), N, r1D) - Vq1D = vandermonde(Line(), N, rq1D) / VDM_1D - invVDM_1D = inv(VDM_1D) - invM_1D = VDM_1D * VDM_1D' - M1D = Vq1D' * diagm(wq1D) * Vq1D - - # form kronecker products of multidimensional matrices to invert/multiply - VDM = kronecker(VDM_1D, VDM_1D, VDM_1D) - invVDM = kronecker(invVDM_1D, invVDM_1D, invVDM_1D) - invM = kronecker(invM_1D, invM_1D, invM_1D) - - M = kronecker(M1D, M1D, M1D) - - _, Vr, Vs, Vt = basis(elem, N, r, s, t) - Dr, Ds, Dt = (A -> A * invVDM).((Vr, Vs, Vt)) - - # low order interpolation nodes - r1, s1, t1 = nodes(elem, 1) - V1 = vandermonde(elem, 1, r, s, t) / vandermonde(elem, 1, r1, s1, t1) - - # Nodes on faces, and face node coordinate - quad_rule_face = tensor_product_quadrature(face_type(elem), approximation_type.data.quad_rule_1D...) - rf, sf, tf, wf, nrJ, nsJ, ntJ = init_face_data(elem, quad_rule_face=quad_rule_face) - - # quadrature nodes - build from 1D nodes. - rq, sq, tq, wq = tensor_product_quadrature(elem, approximation_type.data.quad_rule_1D...) - Vq = kronecker(Vq1D, Vq1D, Vq1D) # vandermonde(elem, N, rq, sq, tq) * invVDM - Pq = invM * (Vq' * diagm(wq)) - - Vf = vandermonde(elem, N, rf, sf, tf) * invVDM - LIFT = invM * (Vf' * diagm(wf)) - - # plotting nodes - rp1D = LinRange(-1, 1, Nplot + 1) - Vp1D = vandermonde(Line(), N, rp1D) / VDM_1D - Vp = kronecker(Vp1D, Vp1D, Vp1D) - rp, sp, tp = vec.(StartUpDG.NodesAndModes.meshgrid(rp1D, rp1D, rp1D)) - - return RefElemData(elem, approximation_type, N, fv, V1, - tuple(r, s, t), VDM, vec(Fmask), - tuple(rp, sp, tp), Vp, - tuple(rq, sq, tq), wq, Vq, - tuple(rf, sf, tf), wf, Vf, tuple(nrJ, nsJ, ntJ), - M, Pq, (Dr, Ds, Dt), LIFT) -end - """ RefElemData(elem::Wedge, approximation_type::Polynomial, N; quad_rule_vol=quad_nodes(elem, N), @@ -265,7 +170,8 @@ end Builds operators for prisms/wedges """ -function RefElemData(elem::Wedge, approximation_type::Polynomial, N; +function RefElemData(elem::Wedge, + approximation_type::Polynomial{MultidimensionalQuadrature}, N; quad_rule_vol=quad_nodes(elem, N), quad_rule_face_quad=quad_nodes(Quad(), N), quad_rule_face_tri=quad_nodes(Tri(), N), @@ -333,7 +239,8 @@ function RefElemData(elem::Wedge, approximation_type::Polynomial, N; end """ - RefElemData(elem::Pyr, approximation_type::Polynomial, N; + RefElemData(elem::Pyr, + approximation_type::Polynomial, N; quad_rule_vol=quad_nodes(elem, N), quad_rule_face_quad=quad_nodes(Quad(), N), quad_rule_face_tri=quad_nodes(Tri(), N), @@ -342,7 +249,8 @@ end Builds operators for pyramids. """ -function RefElemData(elem::Pyr, approximation_type::Polynomial, N; +function RefElemData(elem::Pyr, + approximation_type::Polynomial{MultidimensionalQuadrature}, N; quad_rule_vol=quad_nodes(elem, N), quad_rule_face_quad=quad_nodes(Quad(), N), quad_rule_face_tri=quad_nodes(Tri(), N), @@ -353,10 +261,8 @@ function RefElemData(elem::Pyr, approximation_type::Polynomial, N; fv = face_vertices(elem) #Get interpolation nodes of degree N - r, s, t = nodes(elem, N) - - VDM = vandermonde(elem, N, r, s, t) - + r, s, t = nodes(elem, N) + VDM = vandermonde(elem, N, r, s, t) Fmask = find_face_nodes(elem, r, s, t) # low order interpolation nodes @@ -392,7 +298,7 @@ function RefElemData(elem::Pyr, approximation_type::Polynomial, N; rq, sq, tq, wq = quad_rule_vol Vq, Vrq, Vsq, Vtq = map(A -> A / VDM, basis(elem, N, rq, sq, tq)) - M = Vq' * diagm(wq) * Vq + M = Vq' * diagm(wq) * Vq Pq = M \ (Vq' * diagm(wq)) # We define nodal differentiation matrices using quadrature instead of using @@ -402,12 +308,12 @@ function RefElemData(elem::Pyr, approximation_type::Polynomial, N; # (Du, v) = (du/dx, v) ∀v ∈ pyramid_space Drst = map(Vderiv -> M \ (Vq' * diagm(wq) * Vderiv), (Vrq, Vsq, Vtq)) + LIFT = M \ (Vf' * diagm(wf)) + # plotting nodes rp, sp, tp = equi_nodes(elem, Nplot) Vp = vandermonde(elem, N, rp, sp, tp) / VDM - LIFT = M \ (Vf' * diagm(wf)) - return RefElemData(Pyr(node_ids_by_face), approximation_type, N, fv, V1, tuple(r, s, t), VDM, Fmask, tuple(rp, sp, tp), Vp, @@ -416,19 +322,150 @@ function RefElemData(elem::Pyr, approximation_type::Polynomial, N; M, Pq, Drst, LIFT) end +""" + RefElemData(elem::Union{Quad, Hex}, approximation_type::TensorProductQuadrature, N) + RefElemData(elem::Union{Quad, Hex}, approximation_type::Polynomial{<:TensorProductQuadrature}, N) + +Constructor for quadrilateral and hexahedral `RefElemData` where the quadrature is assumed to have a +tensor product structure. +""" +function RefElemData(elem::Quad, + approximation_type::Polynomial{<:TensorProductQuadrature}, N; + quad_rule_face = approximation_type.data.quad_rule_1D, + Nplot = 10) + + fv = face_vertices(elem) + + # Construct matrices on reference elements + r, s = nodes(elem, N) + Fmask = hcat(find_face_nodes(elem, r, s)...) + + # construct 1D operator for faster Kronecker solves + r1D = nodes(Line(), N) + rq1D, wq1D = approximation_type.data.quad_rule_1D + VDM_1D = vandermonde(Line(), N, r1D) + Vq1D = vandermonde(Line(), N, rq1D) / VDM_1D + invVDM_1D = inv(VDM_1D) + invM_1D = VDM_1D * VDM_1D' + M1D = Vq1D' * diagm(wq1D) * Vq1D + + # form kronecker products of multidimensional matrices to invert/multiply + VDM = kronecker(VDM_1D, VDM_1D) + invVDM = kronecker(invVDM_1D, invVDM_1D) + invM = kronecker(invM_1D, invM_1D) + + M = kronecker(M1D, M1D) + + _, Vr, Vs = basis(elem, N, r, s) + Dr, Ds = (A -> A * invVDM).((Vr, Vs)) + + # low order interpolation nodes + r1, s1 = nodes(elem, 1) + V1 = vandermonde(elem, 1, r, s) / vandermonde(elem, 1, r1, s1) + + # Nodes on faces, and face node coordinate + rf, sf, wf, nrJ, nsJ = init_face_data(elem, quad_rule_face=quad_rule_face) + + # quadrature nodes - build from 1D nodes. + rq, sq, wq = tensor_product_quadrature(elem, approximation_type.data.quad_rule_1D...) + Vq = kronecker(Vq1D, Vq1D) # vandermonde(elem, N, rq, sq, tq) * invVDM + Pq = invM * (Vq' * diagm(wq)) + + Vf = vandermonde(elem, N, rf, sf) * invVDM + LIFT = invM * (Vf' * diagm(wf)) + + # plotting nodes + rp1D = LinRange(-1, 1, Nplot + 1) + Vp1D = vandermonde(Line(), N, rp1D) / VDM_1D + Vp = kronecker(Vp1D, Vp1D, Vp1D) + rp, sp = vec.(StartUpDG.NodesAndModes.meshgrid(rp1D, rp1D)) + + return RefElemData(elem, approximation_type, N, fv, V1, + tuple(r, s), VDM, vec(Fmask), + tuple(rp, sp), Vp, + tuple(rq, sq), wq, Vq, + tuple(rf, sf), wf, Vf, tuple(nrJ, nsJ), + M, Pq, (Dr, Ds), LIFT) +end + +function RefElemData(elem::Hex, + approximation_type::Polynomial{<:TensorProductQuadrature}, N; + quad_rule_face = + tensor_product_quadrature(face_type(elem), + approximation_type.data.quad_rule_1D...), + Nplot = 10) + + fv = face_vertices(elem) + + # Construct matrices on reference elements + r, s, t = nodes(elem, N) + Fmask = hcat(find_face_nodes(elem, r, s, t)...) + + # construct 1D operator for faster Kronecker solves + r1D = nodes(Line(), N) + rq1D, wq1D = approximation_type.data.quad_rule_1D + VDM_1D = vandermonde(Line(), N, r1D) + Vq1D = vandermonde(Line(), N, rq1D) / VDM_1D + invVDM_1D = inv(VDM_1D) + invM_1D = VDM_1D * VDM_1D' + M1D = Vq1D' * diagm(wq1D) * Vq1D + + # form kronecker products of multidimensional matrices to invert/multiply + VDM = kronecker(VDM_1D, VDM_1D, VDM_1D) + invVDM = kronecker(invVDM_1D, invVDM_1D, invVDM_1D) + invM = kronecker(invM_1D, invM_1D, invM_1D) + + M = kronecker(M1D, M1D, M1D) + + _, Vr, Vs, Vt = basis(elem, N, r, s, t) + Dr, Ds, Dt = (A -> A * invVDM).((Vr, Vs, Vt)) + + # low order interpolation nodes + r1, s1, t1 = nodes(elem, 1) + V1 = vandermonde(elem, 1, r, s, t) / vandermonde(elem, 1, r1, s1, t1) + + # Nodes on faces, and face node coordinate + rf, sf, tf, wf, nrJ, nsJ, ntJ = init_face_data(elem, quad_rule_face=quad_rule_face) + + # quadrature nodes - build from 1D nodes. + rq, sq, tq, wq = tensor_product_quadrature(elem, approximation_type.data.quad_rule_1D...) + Vq = kronecker(Vq1D, Vq1D, Vq1D) # vandermonde(elem, N, rq, sq, tq) * invVDM + Pq = invM * (Vq' * diagm(wq)) + + Vf = vandermonde(elem, N, rf, sf, tf) * invVDM + LIFT = invM * (Vf' * diagm(wf)) + + # plotting nodes + rp1D = LinRange(-1, 1, Nplot + 1) + Vp1D = vandermonde(Line(), N, rp1D) / VDM_1D + Vp = kronecker(Vp1D, Vp1D, Vp1D) + rp, sp, tp = vec.(StartUpDG.NodesAndModes.meshgrid(rp1D, rp1D, rp1D)) + + return RefElemData(elem, approximation_type, N, fv, V1, + tuple(r, s, t), VDM, vec(Fmask), + tuple(rp, sp, tp), Vp, + tuple(rq, sq, tq), wq, Vq, + tuple(rf, sf, tf), wf, Vf, tuple(nrJ, nsJ, ntJ), + M, Pq, (Dr, Ds, Dt), LIFT) +end + +RefElemData(elem::Hex, approximation_parameter::TensorProductQuadrature, N; Nplot = 10) = + RefElemData(elem, Polynomial(approximation_parameter), N; Nplot) + + function tensor_product_quadrature(::Line, r1D, w1D) return r1D, w1D end -function tensor_product_quadrature(::Quad, r1D, w1D) +function tensor_product_quadrature(::Union{Tri, Quad}, r1D, w1D) sq, rq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D)) ws, wr = vec.(StartUpDG.NodesAndModes.meshgrid(w1D)) wq = wr .* ws return rq, sq, wq end -function tensor_product_quadrature(::Hex, r1D, w1D) +function tensor_product_quadrature(::Union{Tet, Hex}, r1D, w1D) rq, sq, tq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D, r1D, r1D)) wr, ws, wt = vec.(StartUpDG.NodesAndModes.meshgrid(w1D, w1D, w1D)) wq = wr .* ws .* wt @@ -438,15 +475,16 @@ end """ RefElemData(elem::Union{Line, Quad, Hex}, approximation_type::Polynomial{Gauss}, N) -Builds `rd::RefElemData` with (N+1)-point Gauss quadrature in each dimension. +Builds a `rd::RefElemData` with (N+1)-point Gauss quadrature in each dimension. """ function RefElemData(element_type::Union{Line, Quad, Hex}, - approximation_type::Polynomial{<:Gauss}, N; kwargs...) - # explicitly specify Gauss quadrature rule with (N+1)^d points - quad_rule_vol = tensor_product_quadrature(element_type, StartUpDG.gauss_quad(0, 0, N)...) - rd = RefElemData(element_type, Polynomial(), N; quad_rule_vol) - return @set rd.approximation_type = approximation_type + approximation_type::Polynomial{<:TensorProductGaussCollocation}, + N; kwargs...) + + quadrature_type = TensorProductQuadrature(gauss_quad(0, 0, N)) + rd = RefElemData(element_type, Polynomial(quadrature_type), N; kwargs...) + return @set rd.approximation_type = approximation_type end -RefElemData(element_type::Union{Line, Quad, Hex}, ::Gauss, N; kwargs...) = - RefElemData(element_type, Polynomial{Gauss}(), N; kwargs...) +RefElemData(element_type::Union{Line, Quad, Hex}, ::TensorProductGaussCollocation, N; kwargs...) = + RefElemData(element_type, Polynomial{TensorProductGaussCollocation}(), N; kwargs...) diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index 41b098f0..59664a0b 100644 --- a/src/StartUpDG.jl +++ b/src/StartUpDG.jl @@ -7,28 +7,28 @@ using FillArrays: Fill using HDF5: h5open # used to read in SBP triangular node data using Kronecker: kronecker # for Hex element matrix manipulations using LinearAlgebra: - cond, diagm, eigvals, Diagonal, UniformScaling, I, mul!, norm, qr, ColumnNorm, Symmetric + cond, diagm, eigvals, Diagonal, UniformScaling, I, mul!, norm, qr, ColumnNorm, Symmetric, nullspace, pinv using NodesAndModes: meshgrid, find_face_nodes, face_vertices @reexport using NodesAndModes # for basis functions using PathIntersections: PathIntersections @reexport using PathIntersections: PresetGeometries using Printf: @sprintf using RecipesBase: RecipesBase +@reexport using RecursiveArrayTools: NamedArrayPartition using StaticArrays: SVector, SMatrix using Setfield: setproperties, @set # for "modifying" structs (setproperties) -@reexport using SimpleUnPack: @unpack -using SparseArrays: sparse, droptol!, blockdiag +using SparseArrays: sparse, droptol!, blockdiag, nnz using Triangulate: Triangulate, TriangulateIO, triangulate @reexport using WriteVTK @inline mean(x) = sum(x) / length(x) -# reference element utility functions include("RefElemData.jl") include("RefElemData_polynomial.jl") export RefElemData, Polynomial -export TensorProductQuadrature, Gauss +export MultidimensionalQuadrature, TensorProductQuadrature +export TensorProductGaussCollocation, Gauss include("RefElemData_TensorProductWedge.jl") export TensorProductWedge @@ -36,7 +36,11 @@ export TensorProductWedge include("RefElemData_SBP.jl") export SBP, DefaultSBPType, TensorProductLobatto, Hicken, Kubatko # types for SBP node dispatch export LobattoFaceNodes, LegendreFaceNodes # type parameters for SBP{Kubatko{...}} -export hybridized_SBP_operators, sparse_low_order_SBP_operators +export hybridized_SBP_operators + +include("low_order_sbp.jl") +export sparse_low_order_SBP_operators +export subcell_limiting_operators export inverse_trace_constant, face_type include("ref_elem_utils.jl") @@ -55,16 +59,17 @@ export make_periodic include("boundary_utils.jl") export boundary_face_centroids, tag_boundary_faces, tag_boundary_nodes -# helper array type for cut cell and hybrid meshes -include("named_array_partition.jl") -export NamedArrayPartition - include("hybrid_meshes.jl") export num_faces, num_vertices, HybridMeshExample include("physical_frame_basis.jl") include("cut_cell_meshes.jl") export PhysicalFrame, equi_nodes +export Subtriangulation + +# ! this will be deprecated in a future release +include("cut_cell_moment_fitting.jl") +export MomentFitting include("state_redistribution.jl") export StateRedistribution, apply! @@ -101,20 +106,4 @@ export CircularDomain, PartialCircularDomain include("explicit_timestep_utils.jl") export ck45 # LSERK 45 - -using Requires: @require - -function __init__() - @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin - include("TriangulatePlots.jl") - end - - # Until Julia v1.9 is the minimum required version for StartUpDG.jl, we still support Requires.jl - @static if !isdefined(Base, :get_extension) - @require SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" begin - include("../ext/StartUpDGSummationByPartsOperatorsExt.jl") - end - end -end - end # module diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 3a3a1bc4..c1420a33 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -20,24 +20,15 @@ face interpolation, and lifting operators). The field `cut_cell_data` contains additional data from PathIntersections. """ -struct CutCellMesh{T1, T2, T3, T4, T5} +struct CutCellMesh{T1, T2, T3, T4, T5, T6} physical_frame_elements::T1 cut_face_nodes::T2 objects::T3 cut_cell_operators::T4 cut_cell_data::T5 + quadrature_type::T6 end -# TODO: add isoparametric cut cell mesh with positive quadrature points -# # This mesh type has a polynomial representation of objects, so we don't store the curve info -# struct IsoparametricCutCellMesh{T1, T2, T3, T4} -# physical_frame_elements::T1 -# cut_face_nodes::T2 -# cut_cell_operators::T3 -# cut_cell_data::T4 -# end - - function Base.show(io::IO, ::MIME"text/plain", md::MeshData{DIM, <:CutCellMesh}) where {DIM} @nospecialize md print(io,"Cut-cell MeshData of dimension $DIM with $(md.num_elements) elements " * @@ -45,17 +36,7 @@ function Base.show(io::IO, ::MIME"text/plain", md::MeshData{DIM, <:CutCellMesh}) end # maps x ∈ [-1,1] to [a,b] -map_to_interval(x, a, b) = a + (b-a) * 0.5 * (1 + x) - -function count_cut_faces(cutcells) - num_cut_faces = zeros(Int, length(cutcells)) - for e in eachindex(cutcells) - curve = cutcells[e] - stop_points = curve.stop_pts - num_cut_faces[e] = length(stop_points) - 1 - end - return num_cut_faces -end +map_to_interval(x, a, b) = @. a + (b-a) * 0.5 * (1 + x) function neighbor_across_face(f, ex, ey) if f==1 @@ -71,39 +52,6 @@ function neighbor_across_face(f, ex, ey) end end -function compute_face_centroids(rd, xf, yf, cutcell_data) - - (; region_flags, cut_faces_per_cell, cut_face_offsets ) = cutcell_data - num_cut_cells = length(cut_faces_per_cell) - num_cartesian_cells = sum(region_flags .== 0) - num_cut_faces = sum(cut_faces_per_cell) - - num_points_per_face = length(rd.rf) ÷ num_faces(rd.element_type) - - face_centroids_x = NamedArrayPartition(cartesian=zeros(num_faces(rd.element_type), num_cartesian_cells), - cut=zeros(num_cut_faces)) - face_centroids_y = similar(face_centroids_x) - - for e in 1:num_cartesian_cells - xf_element = reshape(view(xf.cartesian, :, e), num_points_per_face, num_faces(rd.element_type)) - yf_element = reshape(view(yf.cartesian, :, e), num_points_per_face, num_faces(rd.element_type)) - view(face_centroids_x.cartesian, :, e) .= vec(sum(xf_element, dims=1) / num_points_per_face) - view(face_centroids_y.cartesian, :, e) .= vec(sum(yf_element, dims=1) / num_points_per_face) - end - - for e in 1:num_cut_cells - face_node_ids = (1:(num_points_per_face * cut_faces_per_cell[e])) .+ cut_face_offsets[e] * num_points_per_face - xf_element = reshape(view(xf.cut, face_node_ids), num_points_per_face, cut_faces_per_cell[e]) - yf_element = reshape(view(yf.cut, face_node_ids), num_points_per_face, cut_faces_per_cell[e]) - - face_ids = (1:cut_faces_per_cell[e]) .+ cut_face_offsets[e] - view(face_centroids_x.cut, face_ids) .= vec(sum(xf_element, dims=1) / num_points_per_face) - view(face_centroids_y.cut, face_ids) .= vec(sum(yf_element, dims=1) / num_points_per_face) - end - - return face_centroids_x, face_centroids_y -end - function compute_element_indices(region_flags) cells_per_dimension_x, cells_per_dimension_y = size(region_flags) @@ -142,12 +90,9 @@ function is_inside_domain(ex, ey, regions) end end -# generates at least Np_target sampling points within a cut cell defined by `curve` -# returns both x_sampled, y_sampled (physical points inside the cut cell), as well as -# r_sampled, y_sampled (reference points which correspond to x_sampled, y_sampled). -function generate_sampling_points(objects, elem, rd, Np_target; N_sampled = 4 * rd.N) +function generate_sampling_points(objects, elem, Np_target::Int, N_sampled::Int) - r_sampled, s_sampled = equi_nodes(rd.element_type, N_sampled) # oversampled nodes + r_sampled, s_sampled = equi_nodes(Quad(), N_sampled) # oversampled nodes # map sampled points to the background Cartesian cell x_sampled, y_sampled = map_nodes_to_background_cell(elem, r_sampled, s_sampled) @@ -160,7 +105,7 @@ function generate_sampling_points(objects, elem, rd, Np_target; N_sampled = 4 * while sum(is_in_domain) < Np_target N_sampled *= 2 # double degree of sampling - r_sampled, s_sampled = equi_nodes(rd.element_type, N_sampled) # oversampled nodes + r_sampled, s_sampled = equi_nodes(Quad(), N_sampled) # oversampled nodes x_sampled, y_sampled = map_nodes_to_background_cell(elem, r_sampled, s_sampled) # check if all the points are in all the objects @@ -175,240 +120,550 @@ function generate_sampling_points(objects, elem, rd, Np_target; N_sampled = 4 * return x_sampled[ids_in_element], y_sampled[ids_in_element] end -# returns points (xf, yf), scaled normals (nxJ, nyJ), and face Jacobian (Jf) -# for a curve returned from PathIntersections. -# `out` should hold `xf, yf, nxJ, nyJ, Jf = ntuple(_ -> similar(points, (length(points), num_faces)), 5)` -function map_points_to_cut_cell_faces(points, curve) - num_faces = length(curve.subcurves) - out = ntuple(_ -> similar(points, (length(points) * num_faces)), 5) - return map_points_to_cut_cell_faces!(out, points, curve) +is_Cartesian(flag) = flag==0 ? true : false +is_cut(flag) = flag > 0 + +# returns the 1D quadrature used to build a RefElemData surface quadrature +function get_1d_quadrature(rd::RefElemData{2, Quad}) + rf = reshape(rd.rf, :, num_faces(rd.element_type)) + wf = reshape(rd.wf, :, num_faces(rd.element_type)) + + # face ordering on a quad is -/+ x, -/+ y. face 3 = -y + return rf[:, 3], wf[:, 3] end -function map_points_to_cut_cell_faces!(out, points, curve) - xf, yf, nxJ, nyJ, Jf = out - stop_points = curve.stop_pts - r1D = points - for f in 1:length(stop_points)-1 - for i in eachindex(r1D) - fid = i + (f-1) * length(r1D) +function calculate_cutcells(vx, vy, objects, ds = 1e-3, arc_tol = 1e-10, corner_tol = 1e-10) - s = map_to_interval(r1D[i], stop_points[f], stop_points[f+1]) - - # compute tangent vector at a node, face Jacobian, and normals - x_node, y_node = curve(s) - xf[fid], yf[fid] = x_node, y_node - tangent_vector = PathIntersections.ForwardDiff.derivative(curve, s) - - # the face Jacobian involves scaling between mapped and reference face - # reference face = [-1, 1] - scaling = (stop_points[f+1] - stop_points[f]) / 2 - Jf[fid] = norm(tangent_vector) * scaling - - # we have to flip the sign to get the outward normal. - # note: we compute the scaled normal nxJ for consistency with other meshes. - normal_node = SVector{2}(tangent_vector[2], -tangent_vector[1]) - nxJ[fid], nyJ[fid] = (-normal_node / norm(normal_node)) .* Jf[fid] - end + stop_pts = PathIntersections.find_mesh_intersections((vx, vy), objects, ds, arc_tol, corner_tol, + closed_list=true, closure_tol=1e-12) + + # Calculate cutcells + region_flags, cutcell_indices, cutcells = + PathIntersections.define_regions((vx, vy), objects, stop_pts, binary_regions=false) + + cells_per_dimension_x = length(vx) - 1 + cells_per_dimension_y = length(vy) - 1 + + # sort the vector of cut cells so that they match the ordering when + # iterating through Cartesian mesh indices via (ex, ey). + cutcell_ordering = zeros(Int, length(cutcells)) + e = 1 + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + if is_cut(region_flags[ex, ey]) + cutcell_ordering[e] = cutcell_indices[ex, ey] + e += 1 + end end - return vec.((xf, yf, nxJ, nyJ, Jf)) + permute!(cutcells, cutcell_ordering) + + return region_flags, cutcells end +""" + subtriangulated_cutcell_quadrature(cutcell, rd_tri::RefElemData, + r1D = gauss_lobatto_quad(0,0,rd_tri.N)) -# Computes face geometric terms from a RefElemData, `quad_rule_face = (r1D, w1D)`, -# the vectors of the 1D vertex nodes `vx` and `vy`, and named tuple -# `cutcell_data is a NamedTuple containing `objects`, `region_flags`, `stop_pts``, `cutcells`. -function compute_geometric_data(rd::RefElemData{2, Quad}, quad_rule_face, - vx, vy, cutcell_data; tol=100 * eps()) +Constructs a quadrature from subtriangulations. The degree of both the quadrature +rule and isoparametric mapping are specified by `rd_tri`. - # domain size and reference face weights - cells_per_dimension_x, cells_per_dimension_y = length(vx) - 1, length(vy) - 1 - LX, LY = (x -> x[2] - x[1]).(extrema.((vx, vy))) +The optional argument `r1D` specifies the nodal points used for constructing a +curved mapping via interpolatory warp and blend. +""" +function subtriangulated_cutcell_quadrature(cutcell, rd_tri, + r1D = NodesAndModes.nodes(Line(), rd_tri.N)) + # vxy = matrix of the form [x_coordinates, y_coordinates] + vxy = hcat(getindex.(cutcell.(cutcell.stop_pts[1:end-1]), 1), + getindex.(cutcell.(cutcell.stop_pts[1:end-1]), 2)) + (VX, VY), EToV = StartUpDG.triangulate_points(permutedims(vxy)) + + # allocate arrays for face interp coordinates + tri_face_coords = (rd_tri.r[vec(rd_tri.Fmask)], rd_tri.s[vec(rd_tri.Fmask)]) + tri_face_coords_x = zeros(length(r1D), 3) + tri_face_coords_y = zeros(length(r1D), 3) + + # allocate output arrays + xq, yq, wJq = ntuple(_ -> zeros(length(rd_tri.wq), size(EToV, 1)), 3) + + # loop over each triangle, map 1D interpolation points to faces + for e in axes(EToV, 1) + ids = view(EToV, e, :) + for (face_index, face_vertices) in enumerate(SVector{2}.(rd_tri.fv)) + vertices_on_face = sort(ids[face_vertices]) + + # map face interpolation points to a physical element. + + # This assumes PathIntersections.jl uses a clockwise ordering of stop curve points. + # Since StartUpDG uses a CCW ordering, we reverse the order for + for i in eachindex(r1D) + # If the vertex indices are far apart, it's the last face/boundary curve + if (x->abs(x[2]-x[1]))(vertices_on_face) == length(VX) - 1 + s = map_to_interval(r1D[i], reverse(cutcell.stop_pts[end-1:end])...) + point = cutcell(s) + + # if vertex indices are consecutive, it's a boundary face + elseif (x->x[2]-x[1])(vertices_on_face) == 1 + + curve_id = minimum(ids[face_vertices]) + s = map_to_interval(r1D[i], reverse(cutcell.stop_pts[curve_id:curve_id+1])...) + point = cutcell(s) + + else # it's a non-boundary face, it's a straight line + point = SVector{2}.(map_to_interval(r1D[i], VX[ids[face_vertices]]...), + map_to_interval(r1D[i], VY[ids[face_vertices]]...)) + end + tri_face_coords_x[i, face_index] = point[1] + tri_face_coords_y[i, face_index] = point[2] + end + end + + # this operator performs a least squares fit, and is equivalent to isoparametric + # warp and blend. can be moved outside of the cutcell loop for efficiency. + warp_face_points_to_interp = + face_basis(Tri(), rd_tri.N, rd_tri.rst...) / + face_basis(Tri(), rd_tri.N, tri_face_coords...) + + # this performs a least squares fit interpolation by the face basis. It's + # equivalent to isoparametric warp and blend if the face node locations are + # representable by the face basis (e.g., polynomial). + tri_warped_coords_x = warp_face_points_to_interp * vec(tri_face_coords_x) + tri_warped_coords_y = warp_face_points_to_interp * vec(tri_face_coords_y) + + _, _, _, _, Jq_e = + StartUpDG.geometric_factors(tri_warped_coords_x, tri_warped_coords_y, + rd_tri.Vq * rd_tri.Dr, rd_tri.Vq * rd_tri.Ds) + + view(xq, :, e) .= rd_tri.Vq * tri_warped_coords_x + view(yq, :, e) .= rd_tri.Vq * tri_warped_coords_y + @. wJq[:,e] = rd_tri.wq * Jq_e + end - r1D, w1D = quad_rule_face + return xq, yq, wJq +end - (; objects, region_flags, cutcells, cut_faces_per_cell ) = cutcell_data +# quadrature type for cut cell meshes +struct Subtriangulation end - # count number of cells and cut face nodes - num_cartesian_cells = sum(region_flags .== 0) - num_cut_cells = sum(region_flags .> 0) - nodes_per_face = length(r1D) - num_cut_face_nodes = nodes_per_face * sum(cut_faces_per_cell) +""" + MeshData(rd::RefElemData, objects, + vx::AbstractVector, vy::AbstractVector, + quadrature_type=Subtriangulation(); + quad_rule_face=get_1d_quadrature(rd), + precompute_operators=false) - # compute face data - xf, yf, nxJ, nyJ, Jf = ntuple(_ -> NamedArrayPartition(cartesian=zeros(rd.Nfq, num_cartesian_cells), - cut=zeros(num_cut_face_nodes)), 5) +Constructor for MeshData utilizing moment fitting. Does not guarantee positive +quadrature weights, and is slower due to the use of adaptive sampling +to construct - # the face Jacobian involves scaling between mapped and reference domain - # this is precomputed here since it's needed to compute the normals - face_ids_left_right = 1:(length(rd.rf) ÷ 2) - face_ids_top_bottom = ((length(rd.rf) ÷ 2) + 1):length(rd.rf) - Jf.cartesian[face_ids_top_bottom, :] .= LX / (cells_per_dimension_x * sum(w1D)) - Jf.cartesian[face_ids_left_right, :] .= LY / (cells_per_dimension_y * sum(w1D)) +!!! Warning: this may be deprecated or removed in future versions. +""" +MeshData(rd::RefElemData, objects, cells_per_dimension, + quadrature_type=Subtriangulation(); kwargs...) = + MeshData(rd, objects, + cells_per_dimension, cells_per_dimension, + quadrature_type; kwargs...) - # compute face data - e = 1 - for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y - if is_Cartesian(region_flags[ex, ey]) - vx_element = SVector(vx[ex], vx[ex + 1], vx[ex], vx[ex + 1]) - vy_element = SVector(vy[ey], vy[ey], vy[ey + 1], vy[ey + 1]) - x_element, y_element = map(x -> rd.V1 * x, (vx_element, vy_element)) - mul!(view(xf.cartesian, :, e), rd.Vf, x_element) - mul!(view(yf.cartesian, :, e), rd.Vf, y_element) - view(nxJ.cartesian, :, e) .= rd.nrJ .* view(Jf.cartesian, :, e) - view(nyJ.cartesian, :, e) .= rd.nsJ .* view(Jf.cartesian, :, e) - e = e + 1 - end - end +function MeshData(rd::RefElemData, objects, + cells_per_dimension_x::Int, cells_per_dimension_y::Int, + quadrature_type=Subtriangulation(); + coordinates_min=(-1.0, -1.0), coordinates_max=(1.0, 1.0), + kwargs...) - e = 1 - offset = 0 - for (e, curve) in enumerate(cutcells) - # map 1D quadrature points to faces - num_cut_faces = length(curve.subcurves) - fids = (1:length(r1D) * num_cut_faces) .+ offset - out = map(x->view(x, fids), (xf.cut, yf.cut, nxJ.cut, nyJ.cut, Jf.cut)) - map_points_to_cut_cell_faces!(out, r1D, curve) - offset += length(fids) - end - - # compute face points + shifting/scaling coefficients for physical frame cut elements. - physical_frame_elements = PhysicalFrame{2}[] # populate this as we iterate through cut cells - - # store cut-cell scaling/shifting coefficients - (; cut_faces_per_cell, cut_face_offsets ) = cutcell_data - num_points_per_face = length(r1D) + # compute intersections of curve with a background Cartesian grid. + vx = LinRange(coordinates_min[1], coordinates_max[1], cells_per_dimension_x + 1) + vy = LinRange(coordinates_min[2], coordinates_max[2], cells_per_dimension_y + 1) - e = 1 - for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y - if is_cut(region_flags[ex, ey]) + return MeshData(rd, objects, vx, vy, quadrature_type; kwargs...) +end - # here, we evaluate a PhysicalFrame basis by shifting and scaling the - # coordinates on an element back to the reference element [-1, 1]^2. - cut_face_node_ids = (1:num_points_per_face * cut_faces_per_cell[e]) .+ - num_points_per_face * cut_face_offsets[e] +num_cartesian_elements(md::MeshData{DIM, <:CutCellMesh}) where {DIM} = size(md.x.cartesian, 2) +num_cut_elements(md::MeshData{DIM, <:CutCellMesh}) where {DIM} = size(md.x.cut, 2) - # store face nodes (extremal) and coordinates of background Cartesian cell - physical_frame_element = - PhysicalFrame(xf.cut[cut_face_node_ids], yf.cut[cut_face_node_ids], - SVector(vx[ex], vx[ex+1]), SVector(vy[ey], vy[ey+1])) +# this is used in src/MeshData.jl in `getproperty` +function num_elements(md::MeshData{DIM, <:CutCellMesh}) where {DIM} + # the number of elements is given by the number of columns of each component of x + return num_cartesian_elements(md) + num_cut_elements(md) +end - push!(physical_frame_elements, physical_frame_element) +""" + construct_cut_volume_quadrature(N, cutcells; target_degree = 2 * N - 1) - e += 1 +Constructs volume quadrature using subtriangulations of cut cells and +Caratheodory pruning. The resulting quadrature is exact for polynomials +of degree `target_degree`. + +We set `target_degree` to `2 * N - 1` by default, which is sufficient to +ensure that ∫du/dx * v is integrated exactly so that integration by parts +holds under the generated cut cell quadrature. +""" +function construct_cut_volume_quadrature(N, cutcells, physical_frame_elements; + target_degree = 2 * N - 1) + + # We make volume quadrature exact for degree N^2 + N(N-1) + 2(N-1) polynomials on the + # reference element, which ensures that ∫du/dx * v * J is integrated exactly over D̂. + # This is because du/dx * v ∈ P^{2N-1}(D^k) and (du/dx * v) ∘ x(r,s) ∈ P^{(2N-1) * N}). + # + # The minimum exactness of volume quadrature is degree N(N-1) + 2N-2 polynomials, + # which ensures Qh * 1 = 0, where Qh = hybridized SBP operator. + + # Integral ∫ v du/dx J over reference element D̂ + # TODO: why can we decrease this and still see exactness? + N_phys_frame_geo = N^2 + N * (N-1) + 2 * (N-1) + + rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, + quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo)) + + Np_target = StartUpDG.Np_cut(target_degree) + xq_pruned, yq_pruned, wJq_pruned = ntuple(_ -> zeros(Np_target, length(cutcells)), 3) + + for (e, cutcell) in enumerate(cutcells) + + xq, yq, wJq = StartUpDG.subtriangulated_cutcell_quadrature(cutcell, rd_tri) + + # if there are not enough quadrature nodes, refine the target polynomial + # degree until there are. + N_refined = N_phys_frame_geo + while length(wJq) < Np_target + N_refined += 1 + rd_tri_refined = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, + quad_rule_vol=NodesAndModes.quad_nodes_tri(N_refined)) + xq, yq, wJq = StartUpDG.subtriangulated_cutcell_quadrature(cutcell, rd_tri_refined) end - end - # interpolation points - x, y = ntuple(_ -> NamedArrayPartition(cartesian=zeros(rd.Np, num_cartesian_cells), - cut=zeros(Np_cut(rd.N), num_cut_cells)), 2) + # perform Caratheodory pruning + Vtarget = vandermonde(physical_frame_elements[e], target_degree, vec(xq), vec(yq)) + w = vec(wJq) + w_pruned, inds = StartUpDG.caratheodory_pruning_qr(Vtarget, w) - # compute interpolation points on cartesian elements - e = 1 - for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y - if is_Cartesian(region_flags[ex, ey]) - vx_element = SVector(vx[ex], vx[ex + 1], vx[ex], vx[ex + 1]) - vy_element = SVector(vy[ey], vy[ey], vy[ey + 1], vy[ey + 1]) - x_element, y_element = map(x -> rd.V1 * x, (vx_element, vy_element)) - view(x.cartesian, :, e) .= x_element - view(y.cartesian, :, e) .= y_element - e = e + 1 + # test exactness of the pruned quadrature rule if applicable + if target_degree >= 2 * N + V = vandermonde(physical_frame_elements[e], N, vec(xq), vec(yq)) + @assert norm(V' * Diagonal(w) * V - V' * Diagonal(w_pruned) * V) < 100 * eps() end + + @. xq_pruned[:, e] = xq[inds] + @. yq_pruned[:, e] = yq[inds] + @. wJq_pruned[:, e] = w_pruned[inds] end - + return xq_pruned, yq_pruned, wJq_pruned +end + +# construct interpolation points using sampling and approximate Fekete +# points (e.g., QR-DEIM). +function construct_cut_interpolation_nodes(N, objects, physical_frame_elements) + + N_sampled = 4 * N # this is arbitrary + num_cut_elements = length(physical_frame_elements) + x, y = ntuple(_ -> zeros(Np_cut(N), num_cut_elements), 2) + # Compute interpolation points on cut elements for e in eachindex(physical_frame_elements) physical_frame_element = physical_frame_elements[e] x_sampled, y_sampled = - generate_sampling_points(objects, physical_frame_element, rd, 2 * Np_cut(rd.N)) - V = vandermonde(physical_frame_element, rd.N, x_sampled, y_sampled) + generate_sampling_points(objects, physical_frame_element, 2 * Np_cut(N), N_sampled) + + V = vandermonde(physical_frame_element, N, x_sampled, y_sampled) # use pivoted QR to find good interpolation points QRfac = qr(V', ColumnNorm()) - ids = QRfac.p[1:Np_cut(rd.N)] + ids = QRfac.p[1:Np_cut(N)] # if the condition number of the VDM is really bad, then increase the # number of sampled points. condV_original = cond(V[ids, :]) condV = condV_original - N_sampled = 20 * rd.N + N_sampled = 8 * N iter, maxiters = 0, 100 while condV > 1e8 && iter < maxiters x_sampled, y_sampled = - generate_sampling_points(objects, physical_frame_element, rd, 2 * Np_cut(rd.N); - N_sampled = N_sampled) - V = vandermonde(physical_frame_element, rd.N, x_sampled, y_sampled) + generate_sampling_points(objects, physical_frame_element, + 2 * Np_cut(N), N_sampled) + V = vandermonde(physical_frame_element, N, x_sampled, y_sampled) # use pivoted QR to find good interpolation points QRfac = qr(V', ColumnNorm()) - ids = QRfac.p[1:Np_cut(rd.N)] + ids = QRfac.p[1:Np_cut(N)] condV = cond(V[ids, :]) - N_sampled += 5 * rd.N + N_sampled += 2 * N if condV < 1e8 - @warn "Conditioning of old VDM for element $e is $condV_original. " * + @info "Conditioning of old VDM for element $e is $condV_original. " * "After recomputing with $(length(x_sampled)) samples: $condV." end iter += 1 end if iter >= maxiters - @warn "Adaptive sampling of cut-cell element $e: maximum number of iterations $maxiters exceeded." + @warn "Adaptive sampling of cut-cell element $e: " * + "maximum number of iterations $maxiters exceeded." end - view(x.cut, :, e) .= x_sampled[ids] - view(y.cut, :, e) .= y_sampled[ids] + view(x, :, e) .= x_sampled[ids] + view(y, :, e) .= y_sampled[ids] end + return x, y +end - # volume geometric terms - rxJ_cartesian = LX / (2 * cells_per_dimension_x) - syJ_cartesian = LY / (2 * cells_per_dimension_y) - J_cartesian = (LX / cells_per_dimension_x) * (LY / cells_per_dimension_y) / 4 +function construct_cartesian_volume_quadrature(vx, vy, region_flags, quad_rule_volume) - # Note: the volume Jacobian for cut elements is 1 since the "reference element" is the - # cut element itself. Similarly, geometric terms should be 1 since `basis` computes - # physical derivatives accounting for element scaling + rq, sq, wq = quad_rule_volume - # Note: we use FillArrays.Fill instead of FillArrays.Ones and FillArrays.Zeros because - # we store `rxJ, sxJ, ryJ, syJ` in a single SMatrix, which assumes one homogeneous - # type for all the entries. Since FillArrays.Ones/Zeros are distinct types, using them - # results in a type instability when accessing entries of md.rstxyzJ - rxJ = NamedArrayPartition(cartesian=Fill(rxJ_cartesian, rd.Np, num_cartesian_cells), - cut=Fill(1.0, Np_cut(rd.N), num_cut_cells)) - syJ = NamedArrayPartition(cartesian=Fill(syJ_cartesian, rd.Np, num_cartesian_cells), - cut=Fill(1.0, Np_cut(rd.N), num_cut_cells)) - sxJ, ryJ = ntuple(_ -> NamedArrayPartition(cartesian=Fill(0.0, rd.Np, num_cartesian_cells), - cut=Fill(0.0, Np_cut(rd.N), num_cut_cells)), 2) - J = NamedArrayPartition(cartesian = Fill(J_cartesian, rd.Np, num_cartesian_cells), - cut = Fill(1.0, Np_cut(rd.N), num_cut_cells)) + num_cartesian_cells = sum(region_flags .== 0) + xq, yq, wJq = ntuple(_ -> zeros(length(wq), num_cartesian_cells), 3) + + # compute quadrature rules for the Cartesian cells + e = 1 + for ex in axes(region_flags, 1), ey in axes(region_flags, 2) + if is_Cartesian(region_flags[ex, ey]) + dx = vx[ex+1] - vx[ex] + dy = vy[ey+1] - vy[ey] + J = dx * dy / sum(wq) + @. xq[:, e] = dx * 0.5 * (1 + rq) + vx[ex] + @. yq[:, e] = dy * 0.5 * (1 + sq) + vy[ey] + @. wJq[:, e] = wq * J + e += 1 + end + end + return xq, yq, wJq +end + +function construct_cartesian_surface_quadrature(vx, vy, region_flags, + quad_rule_face) + + LX, LY = (x -> x[2] - x[1]).(extrema.((vx, vy))) + cells_per_dimension_x, cells_per_dimension_y = size(region_flags) + r1D, w1D = quad_rule_face + + num_cartesian_cells = sum(region_flags .== 0) + Jf = zeros(num_faces(Quad()), num_cartesian_cells) + xf, yf, nxJ, nyJ, wf = ntuple(_ -> zeros(num_faces(Quad()) * length(r1D), + num_cartesian_cells), 5) + + # the face Jacobian involves scaling between mapped and reference domain + # this is precomputed here since it's needed to compute the normals + + # face 1 and 2: ±s faces + face_ids = [eachindex(r1D); eachindex(r1D) .+ length(r1D)] + Jf[SVector(1, 2), :] .= (LY / (cells_per_dimension_y * sum(w1D))) + + # face 3 and 4: ±r faces + face_ids = [eachindex(r1D) .+ 2 * length(r1D); eachindex(r1D) .+ 3 * length(r1D)] + Jf[SVector(3, 4), :] .= (LX / (cells_per_dimension_x * sum(w1D))) + + # compute face data + e = 1 + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + if is_Cartesian(region_flags[ex, ey]) + # face 1: -r + face_index = 1 + face_ids = eachindex(r1D) .+ (face_index - 1) * length(r1D) + xf[face_ids, e] .= vx[ex] + yf[face_ids, e] .= map_to_interval(r1D, vy[ey], vy[ey+1]) + nxJ[face_ids, e] .= -Jf[face_index, e] + nyJ[face_ids, e] .= zero(eltype(nyJ)) + wf[face_ids, e] .= w1D + + # face 2: +r + face_index = 2 + face_ids = eachindex(r1D) .+ (face_index - 1) * length(r1D) + xf[face_ids, e] .= vx[ex+1] + yf[face_ids, e] .= map_to_interval(r1D, vy[ey], vy[ey+1]) + nxJ[face_ids, e] .= Jf[face_index, e] + nyJ[face_ids, e] .= zero(eltype(nyJ)) + wf[face_ids, e] .= w1D + + # face 3: -s + face_index = 3 + face_ids = eachindex(r1D) .+ (face_index - 1) * length(r1D) + xf[face_ids, e] .= map_to_interval(r1D, vx[ex], vx[ex+1]) + yf[face_ids, e] .= vy[ey] + nxJ[face_ids, e] .= zero(eltype(nxJ)) + nyJ[face_ids, e] .= -Jf[face_index, e] + wf[face_ids, e] .= w1D + + # face 3: +s + face_index = 4 + face_ids = eachindex(r1D) .+ (face_index - 1) * length(r1D) + xf[face_ids, e] .= map_to_interval(r1D, vx[ex], vx[ex+1]) + yf[face_ids, e] .= vy[ey+1] + nxJ[face_ids, e] .= zero(eltype(nxJ)) + nyJ[face_ids, e] .= Jf[face_index, e] + wf[face_ids, e] .= w1D + + e = e + 1 + end + end + + return xf, yf, nxJ, nyJ, wf +end + +""" + construct_physical_frame_elements(region_flags, cutcells) + +Computes physical frame shifting and scaling parameters from the +vertices of cut cells and the background cell location. +""" +function construct_physical_frame_elements(region_flags, vx, vy, cutcells) + + physical_frame_elements = typeof(PhysicalFrame())[] # populate this as we iterate through cut cells + e = 1 + for ex in axes(region_flags, 1), ey in axes(region_flags, 2) + if StartUpDG.is_cut(region_flags[ex, ey]) - rstxyzJ = SMatrix{2, 2, typeof(rxJ), 4}(rxJ, sxJ, ryJ, syJ) # pack geometric terms together + # get extremal points (vertices) of the cut cell + cutcell = cutcells[e] + coordinates = cutcell.(cutcell.stop_pts[1:end-1]) - return physical_frame_elements, x, y, rstxyzJ, J, xf, yf, nxJ, nyJ, Jf + # store vertex nodes and coordinates of background Cartesian cell + physical_frame_element = + PhysicalFrame(getindex.(coordinates, 1), getindex.(coordinates, 2), + SVector(vx[ex], vx[ex+1]), SVector(vy[ey], vy[ey+1])) + + push!(physical_frame_elements, physical_frame_element) + + e += 1 + end + end + return physical_frame_elements end """ - connect_mesh(rd, face_centroids, region_flags, cutcells; tol = 1e2 * eps()) + construct_cut_surface_quadrature(N, cutcells, quad_rule_1D = gauss_quad(0, 0, N)) + +Constructs cut surface quadrature using a degree `N` geometric mapping and a reference +quadrature rule `quad_rule_1D`. Returns `xf, yf, nxJ, nyJ, wf` which are vectors, and +`face_node_indices`, which is a `Vector{Vector{Int}}` of global face node indices (which +index into `xf.cut`, `yf.cut`, etc) for each face of each cut element. + +On boundaries of cut cells, the surface quadrature is taken to be exact for degree +`N^2 + (N-1)` polynomials. This ensures satisfaction of a weak GSBP property. +""" +function construct_cut_surface_quadrature(N, cutcells, + quad_rule_1D = gauss_quad(0, 0, N), + quad_rule_boundary = get_boundary_quadrature(N)) + + rd_line = RefElemData(Line(), N, quad_rule_vol = quad_rule_1D) + rq_boundary, wq_boundary = quad_rule_boundary + + interp_to_cut_boundary = vandermonde(Line(), rd_line.N, rq_boundary) / rd_line.VDM + + xf, yf, nxJ, nyJ, wf = ntuple(_ -> Vector{Float64}[], 5) + + face_node_indices = Vector{UnitRange{Int64}}[] + + # recompute normals using isoparametric mapping on cut cells + for cutcell in cutcells + xf_element, yf_element, nxJ_element, nyJ_element, wf_element = + ntuple(_ -> Vector{Float64}[], 5) + + for f in 1:length(cutcell.stop_pts)-1 + + is_boundary_face = + !(cutcell.subcurves[f] isa PathIntersections.PresetGeometries.Line) + + # switch between lower order and higher order face quadrature + # if it's the curved boundary, it may be higher order + interp_to_face_quad_pts = (is_boundary_face) ? interp_to_cut_boundary : rd_line.Vq + + # map 1D interpolation points to curved faces + points = map(s -> cutcell(StartUpDG.map_to_interval(s, cutcell.stop_pts[f:f+1]...)), rd_line.r) + xf_face = getindex.(points, 1) + yf_face = getindex.(points, 2) + + # interp face coordinates to face quad nodes + xfq_face = interp_to_face_quad_pts * xf_face + yfq_face = interp_to_face_quad_pts * yf_face + + # compute tangent vector at quadrature points using + # the polynomial geometric mapping. + dxdr = interp_to_face_quad_pts * rd_line.Dr * xf_face + dydr = interp_to_face_quad_pts * rd_line.Dr * yf_face + scaled_normal = SVector.(-dydr, dxdr) + nxJ_face = @. getindex(scaled_normal, 1) + nyJ_face = @. getindex(scaled_normal, 2) + wf_face = (is_boundary_face) ? wq_boundary : rd_line.wq + + push!(xf_element, xfq_face) + push!(yf_element, yfq_face) + push!(nxJ_element, nxJ_face) + push!(nyJ_element, nyJ_face) + push!(wf_element, wf_face) + end + + # compute indices of nodes on each face before flattening + num_nodes_per_face = length.(xf_element) + face_node_indices_on_this_element = + map((length, offset) -> (1:length) .+ offset, + num_nodes_per_face, [0; cumsum(num_nodes_per_face[1:end-1])]) + push!(face_node_indices, face_node_indices_on_this_element) + + # flatten the collection of quadrature points for each cutcell + # and push them to the arrays xf, yf, ... + map((x,y) -> push!(x, vcat(y...)), + (xf, yf, nxJ, nyJ, wf), + (xf_element, yf_element, nxJ_element, nyJ_element, wf_element)) + end + + # compute global cut face node indices + num_nodes_per_cut_element = map(x -> maximum(maximum.(x)), face_node_indices) + face_node_offsets = [0; cumsum(num_nodes_per_cut_element[1:end-1])] + for e in eachindex(face_node_indices) + for f in eachindex(face_node_indices[e]) + face_node_indices[e][f] = face_node_indices[e][f] .+ face_node_offsets[e] + end + end + + return xf, yf, nxJ, nyJ, wf, face_node_indices +end + +# this assumes Cartesian cells are Quads +function compute_face_centroids(vx, vy, region_flags, + cutcells::AbstractVector{<:PathIntersections.PiecewiseCurve}) + + num_cartesian_cells = sum(region_flags .== 0) + face_centroids_cartesian = + zeros(SVector{2, Float64}, num_faces(Quad()), num_cartesian_cells) -Connects faces of a cut mesh to each other, returns `FToF` such that face -`f` is connected to `FToF[f]`. + num_faces_per_cut_cell = map(cutcell->length(cutcell.subcurves), cutcells) + face_centroids_cut = [zeros(SVector{2, Float64}, num_faces_per_cut_cell[e]) + for e in eachindex(cutcells)] + + # calculate Cartesian centroids + e = 1 + for ex in axes(region_flags, 1), ey in axes(region_flags, 2) + if is_Cartesian(region_flags[ex, ey]) + # quad faces are -r, +r, -s, +s + face_centroids_cartesian[1, e] = SVector(vx[ex], 0.5 * (vy[ey+1] + vy[ey])) + face_centroids_cartesian[2, e] = SVector(vx[ex+1], 0.5 * (vy[ey+1] + vy[ey])) + face_centroids_cartesian[3, e] = SVector(0.5 * (vx[ex] + vx[ex+1]), vy[ey]) + face_centroids_cartesian[4, e] = SVector(0.5 * (vx[ex] + vx[ex+1]), vy[ey+1]) + e += 1 + end + end -Inputs: -- rd::RefElemData -- face_centroids = (face_centroids_x, face_centroids_y), where `face_centroids_x/y` - are vectors of coordinates of face centroids -- `region_flags`, `cutcells` are return arguments from `PathIntersections.define_regions` -The keyword argument `tol` is the tolerance for matches between face centroids. -""" -function connect_mesh(rd, face_centroids, cutcell_data; tol = 1e2 * eps()) + for (e, cutcell) in enumerate(cutcells) + for (f, face_parametrization) in enumerate(cutcell.subcurves) + s_midpoint = 0.5 * sum(cutcell.sub_bounds[f]) + x, y = face_parametrization(s_midpoint) + face_centroids_cut[e][f] = SVector(x, y) + end + end - (; region_flags, cut_faces_per_cell, cut_face_offsets ) = cutcell_data + return face_centroids_cartesian, face_centroids_cut +end + +num_faces(cutcell::PathIntersections.PiecewiseCurve) = length(cutcell.subcurves) + +function connect_mesh(face_centroids, region_flags, + cutcells::Vector{<:PathIntersections.PiecewiseCurve}; + tol = 1e2 * eps()) - cells_per_dimension_x, cells_per_dimension_y = size(region_flags) num_cartesian_cells = sum(region_flags .== 0) + cut_faces_per_cell = num_faces.(cutcells) + cut_face_offsets = [0; cumsum(cut_faces_per_cell)[1:end-1]] num_cut_faces = sum(cut_faces_per_cell) - num_total_faces = num_cut_faces + num_faces(rd.element_type) * num_cartesian_cells + num_total_faces = num_cut_faces + num_faces(Quad()) * num_cartesian_cells # element_indices[ex, ey] returns the global (flattened) element index into # the arrays `xf.cartesian[:, e]` or `xf.cut[:, e]` @@ -424,7 +679,7 @@ function connect_mesh(rd, face_centroids, cutcell_data; tol = 1e2 * eps()) # flat-sided faces. It wouldn't work for meshes with curved interior interfaces. FToF = collect(1:num_total_faces) - for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + for ex in axes(region_flags, 1), ey in axes(region_flags, 2) e = element_indices[ex, ey] @@ -434,7 +689,7 @@ function connect_mesh(rd, face_centroids, cutcell_data; tol = 1e2 * eps()) # `e_nbr` because the ordering of faces is different for cut elements # and Cartesian elements. if is_Cartesian(region_flags[ex, ey]) - face_ids = (1:num_faces(rd.element_type)) .+ (e-1) * num_faces(rd.element_type) + face_ids = (1:num_faces(Quad())) .+ (e-1) * num_faces(Quad()) elseif is_cut(region_flags[ex, ey]) face_ids = (1:cut_faces_per_cell[e]) .+ cut_face_offsets[e] @@ -453,7 +708,7 @@ function connect_mesh(rd, face_centroids, cutcell_data; tol = 1e2 * eps()) # determine face indices of neighboring cells if is_Cartesian(region_flags[ex_nbr, ey_nbr]) - nbr_face_ids = (1:num_faces(rd.element_type)) .+ (e_nbr-1) * num_faces(rd.element_type) + nbr_face_ids = (1:num_faces(Quad())) .+ (e_nbr-1) * num_faces(Quad()) elseif is_cut(region_flags[ex_nbr, ey_nbr]) nbr_face_ids = (1:cut_faces_per_cell[e_nbr]) .+ cut_face_offsets[e_nbr] @@ -471,8 +726,6 @@ function connect_mesh(rd, face_centroids, cutcell_data; tol = 1e2 * eps()) xy_nbr = SVector(face_centroids_x[j], face_centroids_y[j]) if norm(xy - xy_nbr) < tol * max(1, norm(xy), norm(xy_nbr)) FToF[i] = j - # println("match found for f = $f, e=($ex, $ey), - # enbr=($ex_nbr, $ey_nbr)") end end end @@ -485,242 +738,165 @@ function connect_mesh(rd, face_centroids, cutcell_data; tol = 1e2 * eps()) return FToF end -is_Cartesian(flag) = flag==0 ? true : false -is_cut(flag) = flag > 0 - -# returns the 1D quadrature used to build a RefElemData surface quadrature -function get_1d_quadrature(rd::RefElemData{2, Quad}) - rf = reshape(rd.rf, :, num_faces(rd.element_type)) - wf = reshape(rd.wf, :, num_faces(rd.element_type)) - - # face ordering on a quad is -/+ x, -/+ y. face 3 = -y - return rf[:, 3], wf[:, 3] -end - -function calculate_cutcells(vx, vy, objects, ds = 1e-3, arc_tol = 1e-10, corner_tol = 1e-10) - - stop_pts = PathIntersections.find_mesh_intersections((vx, vy), objects, ds, arc_tol, corner_tol, - closed_list=true, closure_tol=1e-12) - - # Calculate cutcells - region_flags, cutcell_indices, cutcells = - PathIntersections.define_regions((vx, vy), objects, stop_pts, binary_regions=false) - - cells_per_dimension_x = length(vx) - 1 - cells_per_dimension_y = length(vy) - 1 - - # sort the vector of cut cells so that they match the ordering when - # iterating through Cartesian mesh indices via (ex, ey). - cutcell_ordering = zeros(Int, length(cutcells)) - e = 1 - for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y - if is_cut(region_flags[ex, ey]) - cutcell_ordering[e] = cutcell_indices[ex, ey] - e += 1 - end - end - permute!(cutcells, cutcell_ordering) - - return region_flags, cutcells +# If we wish to exactly integrate ∫ u * v * (nx * Jf) and ∫ u * v * (ny * Jf), +# the reference quadrature must be exact for degree 2 * N^2 + (N-1) polynomials +# on the reference line. +# +# Note that we don't perform Caratheodory pruning for boundary faces at the moment. +get_boundary_quadrature(rd::RefElemData) = get_boundary_quadrature(rd.N) +function get_boundary_quadrature(N) + # N_boundary = ceil(Int, (N^2 + (N-1) - 1) / 2) # to ensure weak SBP + N_boundary = ceil(Int, (2 * N^2 + (N - 1) - 1) / 2) # to ensure GSBP property + return gauss_quad(0, 0, N_boundary) end """ function MeshData(rd, geometry, vxyz...) -Creates a cut-cell mesh where the boundary is given by `geometry`, which should be a tuple of functions. -These functions can be generated using PathIntersections.PresetGeometries, for example: +Creates a cut-cell mesh where the boundary is given by `geometry`, which should be a +tuple of functions. These functions can be generated using PathIntersections.PresetGeometries, +for example: ```julia julia> geometry = (PresetGeometries.Circle(R=0.33, x0=0, y0=0), ) ``` -Here, `coordinates_min`, `coordinates_max` contain `(smallest value of x, smallest value of y)` and -`(largest value of x, largest value of y)`, and `cells_per_dimension_x/y` is the number of Cartesian grid -cells placed along each dimension. +Here, `coordinates_min`, `coordinates_max` contain `(smallest value of x, smallest value of y)` +and `(largest value of x, largest value of y)`, and `cells_per_dimension_x/y` is the number of +Cartesian grid cells placed along each dimension. """ -MeshData(rd::RefElemData, objects, cells_per_dimension; kwargs...) = - MeshData(rd::RefElemData, objects, cells_per_dimension, cells_per_dimension; kwargs...) function MeshData(rd::RefElemData, objects, - cells_per_dimension_x::Int, cells_per_dimension_y::Int; - quad_rule_face = get_1d_quadrature(rd), - coordinates_min=(-1.0, -1.0), coordinates_max=(1.0, 1.0), - precompute_operators=false) - - # compute intersections of curve with a background Cartesian grid. - vx = LinRange(coordinates_min[1], coordinates_max[1], cells_per_dimension_x + 1) - vy = LinRange(coordinates_min[2], coordinates_max[2], cells_per_dimension_y + 1) - - return MeshData(rd, objects, vx, vy; quad_rule_face, precompute_operators) -end - -function MeshData(rd::RefElemData, objects, - vx::AbstractVector, vy::AbstractVector; + vx::AbstractVector, vy::AbstractVector, + quadrature_type::Subtriangulation; quad_rule_face=get_1d_quadrature(rd), + quad_rule_boundary=get_boundary_quadrature(rd), + cut_quadrature_target_degree = 2 * rd.N-1, precompute_operators=false) cells_per_dimension_x = length(vx) - 1 cells_per_dimension_y = length(vy) - 1 + + LX, LY = (x -> x[2] - x[1]).(extrema.((vx, vy))) # compute mesh intersections and cut cell elements. # `regions` is a matrix of dimensions `(cells_per_dimension_x, cells_per_dimension_y)` with 3 values: # * 1: cut cell # * 0: Cartesian cell # * -1: excluded cells (in the cut-out region) - region_flags, cutcells = calculate_cutcells(vx, vy, objects) + region_flags, cutcells = StartUpDG.calculate_cutcells(vx, vy, objects) # pack useful cut cell information together. num_cartesian_cells = sum(region_flags .== 0) num_cut_cells = sum(region_flags .> 0) - cut_faces_per_cell = count_cut_faces(cutcells) - cut_face_offsets = [0; cumsum(cut_faces_per_cell)[1:end-1]] - cutcell_data = (; objects, region_flags, cutcells, cut_faces_per_cell, cut_face_offsets) + + N = rd.N - # Compute volume, face points, and physical frame element scalings - physical_frame_elements, x, y, rstxyzJ, J, xf, yf, nxJ, nyJ, Jf = - compute_geometric_data(rd, quad_rule_face, vx, vy, cutcell_data) + #################################################### + # Construct cut cells stuff # + #################################################### + + physical_frame_elements = + construct_physical_frame_elements(region_flags, vx, vy, cutcells) + + x_cut, y_cut = + construct_cut_interpolation_nodes(N, objects, physical_frame_elements) + + xq_cut, yq_cut, wJq_cut = + construct_cut_volume_quadrature(N, cutcells, physical_frame_elements; + target_degree = cut_quadrature_target_degree) + + # On the curved cut boundary, we want to exactly integrate degree N^2 + (N-1) + # polynomials since our target integrand is nxJ * u (nxJ is degree N-1 on the + # reference face, and u is degree N on the physical face, but mapped back to + # the reference face it's N^2). quad_rule_boundary defaults to this degree. + xf_cut, yf_cut, nxJ_cut, nyJ_cut, wf_cut, cut_face_node_indices = + construct_cut_surface_quadrature(N, cutcells, quad_rule_face, quad_rule_boundary) + + #################################################### + # Construct Cartesian cells # + #################################################### + + xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wf_cartesian = + construct_cartesian_surface_quadrature(vx, vy, region_flags, quad_rule_face) + + xq_cartesian, yq_cartesian, wJq_cartesian = + construct_cartesian_volume_quadrature(vx, vy, region_flags, + (rd.rq, rd.sq, rd.wq)) + + # reuse quadrature mapping to construct interpolation nodes + x_cartesian, y_cartesian, _ = + construct_cartesian_volume_quadrature(vx, vy, region_flags, + (rd.r, rd.s, ones(size(rd.r)))) + + #################################################### + # Assemble combined cut/Cartesian arrays # + #################################################### + + x, y = map((x, y) -> NamedArrayPartition((; cartesian=x, cut=y)), + (x_cartesian, y_cartesian), (x_cut, y_cut)) + + xq, yq, wJq = map((x, y) -> NamedArrayPartition((; cartesian=x, cut=y)), + (xq_cartesian, yq_cartesian, wJq_cartesian), + (xq_cut, yq_cut, wJq_cut)) + + xf, yf, nxJ, nyJ, wf = + map((x, y) -> NamedArrayPartition((; cartesian=x, cut=y)), + (xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wf_cartesian), + map(x -> vcat(x...), (xf_cut, yf_cut, nxJ_cut, nyJ_cut, wf_cut))) + + Jf = @. sqrt(nxJ^2 + nyJ^2) + wJf = @. wf * Jf - # Compute face-to-face connectivity by matching face centroids - face_centroids = compute_face_centroids(rd, xf, yf, cutcell_data) - FToF = connect_mesh(rd, face_centroids, cutcell_data) - - # Compute node-to-node mappings - # !!! Warning: this only works if the same quadrature rule is used for all faces! - num_total_faces = length(FToF) - num_points_per_face = length(rd.rf) ÷ num_faces(rd.element_type) - mapM = collect(reshape(1:num_points_per_face * num_total_faces, - num_points_per_face, num_total_faces)) + #################################################### + # Mesh connectivity # + #################################################### + + # note: face_centroids_cartesian has dims [4, num_elements] + # face_centroids_cut has dims [num_elements][num_faces_on_this_element] + face_centroids_cartesian, face_centroids_cut = + StartUpDG.compute_face_centroids(vx, vy, region_flags, cutcells) + + face_centroids_x = NamedArrayPartition(cartesian = getindex.(face_centroids_cartesian, 1), + cut = getindex.(vcat(face_centroids_cut...), 1)) + face_centroids_y = NamedArrayPartition(cartesian = getindex.(face_centroids_cartesian, 2), + cut = getindex.(vcat(face_centroids_cut...), 2)) + face_centroids = (face_centroids_x, face_centroids_y) + FToF = StartUpDG.connect_mesh(face_centroids, region_flags, cutcells) + + # create arrays of node indices on each face, e.g., + # xf.cut[cut_face_node_indices_by_face[f]] gives nodes on face f + cartesian_nodes_per_face = rd.Nfq ÷ rd.num_faces + cartesian_face_node_indices_by_face = + [(1:cartesian_nodes_per_face) .+ (f - 1) * cartesian_nodes_per_face + for f in 1:num_cartesian_cells * rd.num_faces] + cut_face_node_indices_by_face = vcat(cut_face_node_indices...) + + # create global index vector (note that we offset cut cell indices + # by num_cartesian nodes). + face_node_indices_by_face = vcat(cartesian_face_node_indices_by_face, + map(x -> x .+ (num_cartesian_cells * rd.Nfq), + cut_face_node_indices_by_face)) + + # create node mappings + mapM_cartesian = num_cartesian_cells > 0 ? + reshape(eachindex(xf.cartesian), :, num_cartesian_cells) : Int64[] + mapM = NamedArrayPartition(cartesian = collect(mapM_cartesian), + cut = collect(eachindex(xf.cut) .+ length(mapM_cartesian))) mapP = copy(mapM) - p = zeros(Int, num_points_per_face) # temp storage for a permutation vector - for f in eachindex(FToF) - idM, idP = view(mapM, :, f), view(mapM, :, FToF[f]) - xyzM = (view(xf, idM), view(yf, idM)) - xyzP = (view(xf, idP), view(yf, idP)) - StartUpDG.match_coordinate_vectors!(p, xyzM, xyzP) - mapP[p, f] .= idP - end - mapB = findall(vec(mapM) .== vec(mapP)) # determine boundary nodes - - # compute cut-cell surface quadrature - _, w1D = quad_rule_face - wJf = similar(Jf) - wJf.cartesian = reshape(Diagonal(w1D) * reshape(Jf.cartesian, length(w1D), :), size(Jf.cartesian)) - wJf.cut = reshape(Diagonal(w1D) * reshape(Jf.cut, length(w1D), :), size(Jf.cut)) - - # The minimum number of cut cell quadrature points is `Np_cut(2 * rd.N)`. However, - # oversampling slightly seems to improve the conditioning of the quadrature weights. - num_cut_quad_points = Np_cut(2 * rd.N) + 1 - xq, yq, wJq = ntuple(_ -> NamedArrayPartition(cartesian=zeros(rd.Nq, num_cartesian_cells), - cut=zeros(num_cut_quad_points, num_cut_cells)), 3) - # compute quadrature rules for the Cartesian cells - e = 1 - for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y - if is_Cartesian(region_flags[ex, ey]) - dx = vx[ex+1] - vx[ex] - dy = vy[ey+1] - vy[ey] - J = dx * dy / sum(rd.wq) - @. xq.cartesian[:, e] = dx * 0.5 * (1 + rd.rq) + vx[ex] - @. yq.cartesian[:, e] = dy * 0.5 * (1 + rd.sq) + vy[ey] - @. wJq.cartesian[:, e] = rd.wq * J - e += 1 + # allocate a vector for points on non-boundary faces + p = zeros(Int, length(first(quad_rule_face))) + for f in eachindex(FToF) + # if it's not a boundary face + if f !== FToF[f] + idM = view(mapM, face_node_indices_by_face[f]) + idP = view(mapM, face_node_indices_by_face[FToF[f]]) + xyzM = (view(xf, idM), view(yf, idM)) + xyzP = (view(xf, idP), view(yf, idP)) + StartUpDG.match_coordinate_vectors!(p, xyzM, xyzP) + @. mapP[idM[p]] = idP end end - - # polynomial antidifferentiation operators for computing volume integrals - Ix, Iy = StartUpDG.antidiff_operators(2 * rd.N) - - # refine the surface rule used to compute the volume quadrature - if length(quad_rule_face[1]) < 3 * rd.N + 1 - r1D, w1D = gauss_quad(0, 0, 3 * rd.N) - else - r1D, w1D = quad_rule_face - end - # compute quadrature rules for the cut cells. integrate exactly degree 2N basis - for e in eachindex(cutcells) - # compute these quantities using a higher accuracy surface quadrature rule - xf_element, yf_element, nxJ_element, nyJ_element, Jf_element = - map_points_to_cut_cell_faces(r1D, cutcells[e]) - nx_element = nxJ_element ./ Jf_element - ny_element = nyJ_element ./ Jf_element - wJf_element = vec(Diagonal(w1D) * reshape(Jf_element, length(w1D), :)) - - # compute volume integrals using numerical Green's theorem and surface integrals - scaling = physical_frame_elements[e].scaling - Vf = vandermonde(physical_frame_elements[e], 2 * rd.N + 1, xf_element, yf_element) - b = 0.5 * ((Vf * Ix)' * (wJf_element .* nx_element) ./ scaling[1] + - (Vf * Iy)' * (wJf_element .* ny_element) ./ scaling[2]) - - # compute degree 2N basis matrix at sampled points - x_sampled, y_sampled = generate_sampling_points(objects, physical_frame_elements[e], - rd, 2 * Np_cut(rd.N); N_sampled = 8 * rd.N) - Vq = vandermonde(physical_frame_elements[e], 2 * rd.N, x_sampled, y_sampled) - - # if the condition number of the VDM matrix is large, select more points until cond(VDM) < 1e8 - QRfac = qr(Vq', ColumnNorm()) - - # if we have enough points, compute the cond number. Otherwise, assume it's Inf - if length(QRfac.p) >= num_cut_quad_points - ids = view(QRfac.p, 1:num_cut_quad_points) - condV_original = cond(Vq[ids, :]) - else - condV_original = Inf - end - condV = condV_original - N_sampled = 20 * rd.N # start with a lot more points - iter, maxiters = 0, 100 - while condV > 1e8 && iter < maxiters - x_sampled, y_sampled = - generate_sampling_points(objects, physical_frame_elements[e], rd, 2 * Np_cut(rd.N); - N_sampled = N_sampled) - - # if the number of sampled points isn't large enough, add more - if length(x_sampled) < num_cut_quad_points - N_sampled += 5 * rd.N - continue - end - - Vq = vandermonde(physical_frame_elements[e], 2 * rd.N, x_sampled, y_sampled) - - # use pivoted QR to find good interpolation points - QRfac = qr(Vq', ColumnNorm()) - ids = view(QRfac.p, 1:num_cut_quad_points) - condV = cond(Vq[ids, :]) - N_sampled += 5 * rd.N - - if condV < 1e8 - @warn "Conditioning of old VDM for element $e is $condV_original. " * - "After recomputing with $(length(x_sampled)) samples: $condV." - end - - iter += 1 - end - - if iter >= maxiters - @warn "Adaptive sampling of cut-cell element $e: maximum number of iterations $maxiters exceeded." - end - - # naive approach to computing quadrature weights; no guarantees of positivity - QRfac = qr(Vq', ColumnNorm()) - ids = view(QRfac.p, 1:num_cut_quad_points) - wq = Vq[ids, :]' \ b - - quadrature_error = norm(Vq[ids, :]' * wq - b) - quadrature_condition_number = sum(abs.(wq)) / sum(wq) - if quadrature_condition_number > 10 || quadrature_error > 1e3 * eps() - @warn "Quadrature error on element $e is $quadrature_error, " * - "quadrature condition number = $quadrature_condition_number. " * - "Condition number of quadrature VDM is $condV." - end - - if e==249 - if isdefined(Main, :Infiltrator) - Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__) - end - end - - view(xq.cut, :, e) .= view(x_sampled, ids) - view(yq.cut, :, e) .= view(y_sampled, ids) - view(wJq.cut, :, e) .= wq - end + mapB = findall(vec(mapM) .== vec(mapP)) # determine boundary nodes # default to non-periodic is_periodic = (false, false) @@ -741,63 +917,153 @@ function MeshData(rd::RefElemData, objects, cut_cell_data = (; cutcells, region_flags, wJf, cartesian_to_linear_element_indices, linear_to_cartesian_element_indices, - vxyz=(vx, vy), cells_per_dimension=(cells_per_dimension_x, cells_per_dimension_y), # background Cartesian grid info - ) + vxyz=(vx, vy), + cells_per_dimension=(cells_per_dimension_x, + cells_per_dimension_y), + cut_face_node_indices_by_elem_by_face = cut_face_node_indices, + ) # background Cartesian grid info - # get indices of cut face nodes - face_ids(e) = (1:(num_points_per_face * cut_faces_per_cell[e])) .+ - cut_face_offsets[e] * num_points_per_face - cut_face_node_ids = [face_ids(e) for e in 1:num_cut_cells] + # get flattened indices of cut face nodes. note that + # `cut_face_node_indices = [[face_1_indices, face_2_indices, ...], ...]` + cut_face_node_indices_by_cell = + map(x -> UnitRange(first(x[1]), last(x[end])), cut_face_node_indices) if precompute_operators == true - # precompute cut-cell operators and store them in the `md.mesh_type.cut_cell_operators` field - cut_face_nodes = cut_face_node_ids - face_interpolation_matrices = Matrix{eltype(x)}[] - lift_matrices = Matrix{eltype(x)}[] + # precompute cut-cell operators and store them in the `md.mesh_type.cut_cell_operators` field. differentiation_matrices = Tuple{Matrix{eltype(x)}, Matrix{eltype(x)}}[] - mass_matrices = Matrix{eltype(x)}[] + volume_interpolation_matrices, face_interpolation_matrices = + ntuple(_ -> Matrix{eltype(x)}[], 3) + + # TODO: remove these? these can be computed from other matrices + mass_matrices, lift_matrices = ntuple(_ -> Matrix{eltype(x)}[], 2) for (e, elem) in enumerate(physical_frame_elements) VDM = vandermonde(elem, rd.N, x.cut[:, e], y.cut[:, e]) - Vq, Vrq, Vsq = map(A -> A / VDM, basis(elem, rd.N, xq.cut[:,e], yq.cut[:, e])) + Vq, Vrq, Vsq = map(A -> A / VDM, + basis(elem, rd.N, view(xq.cut, :, e), view(yq.cut, :, e))) - M = Vq' * diagm(wJq.cut[:, e]) * Vq - Qr = Vq' * diagm(wJq.cut[:, e]) * Vrq - Qs = Vq' * diagm(wJq.cut[:, e]) * Vsq + M = Vq' * Diagonal(wJq.cut[:, e]) * Vq + Qr = Vq' * Diagonal(wJq.cut[:, e]) * Vrq + Qs = Vq' * Diagonal(wJq.cut[:, e]) * Vsq Dx_e, Dy_e = M \ Qr, M \ Qs - Vf = vandermonde(elem, rd.N, xf.cut[cut_face_nodes[e]], yf.cut[cut_face_nodes[e]]) / VDM + xf_e = xf.cut[cut_face_node_indices_by_cell[e]] + yf_e = yf.cut[cut_face_node_indices_by_cell[e]] + Vf = vandermonde(elem, rd.N, xf_e, yf_e) / VDM - # don't include jacobian scaling in LIFT matrix (for consistency with the Cartesian mesh) - _, w1D = quad_rule_face - num_cut_faces = length(cut_face_nodes[e]) ÷ length(w1D) - wf = vec(repeat(w1D, 1, num_cut_faces)) + # don't include jacobian scaling in LIFT matrix (for consistency + # with the Cartesian mesh) + wf = wJf.cut[cut_face_node_indices_by_cell[e]] ./ + Jf.cut[cut_face_node_indices_by_cell[e]] - push!(lift_matrices, M \ (Vf' * diagm(wf))) + push!(volume_interpolation_matrices, Vq) push!(face_interpolation_matrices, Vf) push!(differentiation_matrices, (Dx_e, Dy_e)) push!(mass_matrices, M) + push!(lift_matrices, M \ (Vf' * Diagonal(wf))) end - cut_cell_operators = (; differentiation_matrices, face_interpolation_matrices, + cut_cell_operators = (; volume_interpolation_matrices, + face_interpolation_matrices, + differentiation_matrices, mass_matrices, lift_matrices) else cut_cell_operators = nothing end - - return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_ids, objects, cut_cell_operators, cut_cell_data), + + + ################################################ + # Construct volume geometric terms # + ################################################ + + rxJ_cartesian = LX / (2 * cells_per_dimension_x) + syJ_cartesian = LY / (2 * cells_per_dimension_y) + J_cartesian = (LX / cells_per_dimension_x) * (LY / cells_per_dimension_y) / 4 + + # Note: the volume Jacobian for cut elements is 1 since the "reference element" is the + # cut element itself. Similarly, geometric terms should be 1 since `basis` computes + # physical derivatives accounting for element scaling + + # Note: we use FillArrays.Fill instead of FillArrays.Ones and FillArrays.Zeros because + # we store `rxJ, sxJ, ryJ, syJ` in a single SMatrix, which assumes one homogeneous + # type for all the entries. Since FillArrays.Ones/Zeros are distinct types, using them + # results in a type instability when accessing entries of md.rstxyzJ + rxJ = NamedArrayPartition(cartesian=Fill(rxJ_cartesian, rd.Np, num_cartesian_cells), + cut=Fill(1.0, Np_cut(rd.N), num_cut_cells)) + syJ = NamedArrayPartition(cartesian=Fill(syJ_cartesian, rd.Np, num_cartesian_cells), + cut=Fill(1.0, Np_cut(rd.N), num_cut_cells)) + sxJ, ryJ = ntuple(_ -> NamedArrayPartition(cartesian=Fill(0.0, rd.Np, num_cartesian_cells), + cut=Fill(0.0, Np_cut(rd.N), num_cut_cells)), 2) + J = NamedArrayPartition(cartesian = Fill(J_cartesian, rd.Np, num_cartesian_cells), + cut = Fill(1.0, Np_cut(rd.N), num_cut_cells)) + + # pack geometric terms together + rstxyzJ = SMatrix{2, 2, typeof(rxJ), 4}(rxJ, sxJ, ryJ, syJ) + + return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_indices_by_cell, objects, + cut_cell_operators, cut_cell_data, quadrature_type), FToF, (x, y), (xf, yf), (xq, yq), wJq, mapM, mapP, mapB, rstxyzJ, J, (nxJ, nyJ), Jf, is_periodic) end -num_cartesian_elements(md::MeshData{DIM, <:CutCellMesh}) where {DIM} = size(md.x.cartesian, 2) -num_cut_elements(md::MeshData{DIM, <:CutCellMesh}) where {DIM} = size(md.x.cut, 2) +""" + hybridized_SBP_operators(md::MeshData{2, <:CutCellMesh}) -# this is used in src/MeshData.jl in `getproperty` -function num_elements(md::MeshData{DIM, <:CutCellMesh}) where {DIM} - # the number of elements is given by the number of columns of each component of x - return num_cartesian_elements(md) + num_cut_elements(md) -end +This constructs hybridized SBP operators using the approach taken in Chan (2019), +"Skew-Symmetric Entropy Stable Modal Discontinuous Galerkin Formulations". +[https://doi.org/10.1007/s10915-019-01026-w](https://doi.org/10.1007/s10915-019-01026-w) + +This function returns `hybridized_operators::Vector{Tuple{<:Matrix, <:Matrix}}` and +`project_and_interp_operators, projection_operators, interpolation_operators`, which are +all `Vector{<:Matrix}`, where each entry corresponds to a cut element. +""" +function hybridized_SBP_operators(md::MeshData{2, <:CutCellMesh}) + mt = md.mesh_type + (; volume_interpolation_matrices, + face_interpolation_matrices, + differentiation_matrices, + mass_matrices) = mt.cut_cell_operators + (; wJf) = mt.cut_cell_data + wf = wJf ./ md.Jf + (; x, nxJ, nyJ) = md + + hybridized_operators = Tuple{Matrix{eltype(x)}, Matrix{eltype(x)}}[] + interpolation_operators, projection_operators, project_and_interp_operators = + ntuple(_ -> Matrix{eltype(x)}[], 3) + for e in eachindex(differentiation_matrices) + Dx, Dy = differentiation_matrices[e] + Vq = volume_interpolation_matrices[e] + Vf = face_interpolation_matrices[e] + M = mass_matrices[e] + + Pq = M \ (Vq' * Diagonal(md.wJq.cut[:, e])) + E = Vf * Pq + + fids = mt.cut_face_nodes[e] + Bx = Diagonal(wf.cut[fids] .* nxJ.cut[fids]) + By = Diagonal(wf.cut[fids] .* nyJ.cut[fids]) + + Qx = Pq' * M * Dx * Pq + Qy = Pq' * M * Dy * Pq + Qxh = 0.5 * [Qx - Qx' E' * Bx; + -Bx * E Bx ] + Qyh = 0.5 * [Qy - Qy' E' * By; + -By * E By ] + + # interpolation matrices + Vh = [Vq; Vf] + VhP = Vh * Pq + Ph = M \ Vh' + + push!(hybridized_operators, (Qxh, Qyh)) + push!(interpolation_operators, Vh) + push!(projection_operators, Ph) + push!(project_and_interp_operators, VhP) + end + + return hybridized_operators, project_and_interp_operators, + projection_operators, interpolation_operators +end \ No newline at end of file diff --git a/src/cut_cell_moment_fitting.jl b/src/cut_cell_moment_fitting.jl new file mode 100644 index 00000000..23310fe2 --- /dev/null +++ b/src/cut_cell_moment_fitting.jl @@ -0,0 +1,632 @@ +# !!! Warning: MomentFitting cut cell quadrature is experimental +# !!! and may be removed in a future release. +struct MomentFitting end + +function MeshData(rd::RefElemData, objects, + vx::AbstractVector, vy::AbstractVector, + quadrature_type::MomentFitting; + quad_rule_face=get_1d_quadrature(rd), + precompute_operators=false) + + cells_per_dimension_x = length(vx) - 1 + cells_per_dimension_y = length(vy) - 1 + + # compute mesh intersections and cut cell elements. + # `regions` is a matrix of dimensions `(cells_per_dimension_x, cells_per_dimension_y)` with 3 values: + # * 1: cut cell + # * 0: Cartesian cell + # * -1: excluded cells (in the cut-out region) + region_flags, cutcells = calculate_cutcells(vx, vy, objects) + + # pack useful cut cell information together. + num_cartesian_cells = sum(region_flags .== 0) + num_cut_cells = sum(region_flags .> 0) + cut_faces_per_cell = count_cut_faces(cutcells) + cut_face_offsets = [0; cumsum(cut_faces_per_cell)[1:end-1]] + cutcell_data = (; objects, region_flags, cutcells, cut_faces_per_cell, cut_face_offsets) + + # Compute volume, face points, and physical frame element scalings + physical_frame_elements, x, y, rstxyzJ, J, xf, yf, nxJ, nyJ, Jf = + compute_geometric_data(rd, quad_rule_face, vx, vy, cutcell_data) + + # Compute face-to-face connectivity by matching face centroids + face_centroids = compute_face_centroids(rd, xf, yf, cutcell_data) + FToF = connect_mesh(rd, face_centroids, cutcell_data) + + # Compute node-to-node mappings + # !!! Warning: this only works if the same quadrature rule is used for all faces! + num_total_faces = length(FToF) + num_points_per_face = length(rd.rf) ÷ num_faces(rd.element_type) + mapM = collect(reshape(1:num_points_per_face * num_total_faces, + num_points_per_face, num_total_faces)) + mapP = copy(mapM) + p = zeros(Int, num_points_per_face) # temp storage for a permutation vector + for f in eachindex(FToF) + idM, idP = view(mapM, :, f), view(mapM, :, FToF[f]) + xyzM = (view(xf, idM), view(yf, idM)) + xyzP = (view(xf, idP), view(yf, idP)) + StartUpDG.match_coordinate_vectors!(p, xyzM, xyzP) + mapP[p, f] .= idP + end + mapB = findall(vec(mapM) .== vec(mapP)) # determine boundary nodes + + # compute cut-cell surface quadrature + _, w1D = quad_rule_face + wJf = similar(Jf) + wJf.cartesian = reshape(Diagonal(w1D) * reshape(Jf.cartesian, length(w1D), :), size(Jf.cartesian)) + wJf.cut = reshape(Diagonal(w1D) * reshape(Jf.cut, length(w1D), :), size(Jf.cut)) + + # The minimum number of cut cell quadrature points is `Np_cut(2 * rd.N)`. However, + # oversampling slightly seems to improve the conditioning of the quadrature weights. + num_cut_quad_points = Np_cut(2 * rd.N) + 1 + xq, yq, wJq = ntuple(_ -> NamedArrayPartition(cartesian=zeros(rd.Nq, num_cartesian_cells), + cut=zeros(num_cut_quad_points, num_cut_cells)), 3) + + # compute quadrature rules for the Cartesian cells + e = 1 + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + if is_Cartesian(region_flags[ex, ey]) + dx = vx[ex+1] - vx[ex] + dy = vy[ey+1] - vy[ey] + J = dx * dy / sum(rd.wq) + @. xq.cartesian[:, e] = dx * 0.5 * (1 + rd.rq) + vx[ex] + @. yq.cartesian[:, e] = dy * 0.5 * (1 + rd.sq) + vy[ey] + @. wJq.cartesian[:, e] = rd.wq * J + e += 1 + end + end + + # polynomial antidifferentiation operators for computing volume integrals + Ix, Iy = StartUpDG.antidiff_operators(2 * rd.N) + + # refine the surface rule used to compute the volume quadrature + if length(quad_rule_face[1]) < 3 * rd.N + 1 + r1D, w1D = gauss_quad(0, 0, 3 * rd.N) + else + r1D, w1D = quad_rule_face + end + # compute quadrature rules for the cut cells. integrate exactly degree 2N basis + for e in eachindex(cutcells) + # compute these quantities using a higher accuracy surface quadrature rule + xf_element, yf_element, nxJ_element, nyJ_element, Jf_element = + map_points_to_cut_cell_faces(r1D, cutcells[e]) + nx_element = nxJ_element ./ Jf_element + ny_element = nyJ_element ./ Jf_element + wJf_element = vec(Diagonal(w1D) * reshape(Jf_element, length(w1D), :)) + + # compute volume integrals using numerical Green's theorem and surface integrals + scaling = physical_frame_elements[e].scaling + Vf = vandermonde(physical_frame_elements[e], 2 * rd.N + 1, xf_element, yf_element) + b = 0.5 * ((Vf * Ix)' * (wJf_element .* nx_element) ./ scaling[1] + + (Vf * Iy)' * (wJf_element .* ny_element) ./ scaling[2]) + + # compute degree 2N basis matrix at sampled points + x_sampled, y_sampled = generate_sampling_points(objects, physical_frame_elements[e], + rd, 2 * Np_cut(rd.N); N_sampled = 8 * rd.N) + Vq = vandermonde(physical_frame_elements[e], 2 * rd.N, x_sampled, y_sampled) + + # if the condition number of the VDM matrix is large, select more points until cond(VDM) < 1e8 + QRfac = qr(Vq', ColumnNorm()) + + # if we have enough points, compute the cond number. Otherwise, assume it's Inf + if length(QRfac.p) >= num_cut_quad_points + ids = view(QRfac.p, 1:num_cut_quad_points) + condV_original = cond(Vq[ids, :]) + else + condV_original = Inf + end + condV = condV_original + N_sampled = 20 * rd.N # start with a lot more points + iter, maxiters = 0, 100 + while condV > 1e8 && iter < maxiters + x_sampled, y_sampled = + generate_sampling_points(objects, physical_frame_elements[e], rd, 2 * Np_cut(rd.N); + N_sampled = N_sampled) + + # if the number of sampled points isn't large enough, add more + if length(x_sampled) < num_cut_quad_points + N_sampled += 5 * rd.N + continue + end + + Vq = vandermonde(physical_frame_elements[e], 2 * rd.N, x_sampled, y_sampled) + + # use pivoted QR to find good interpolation points + QRfac = qr(Vq', ColumnNorm()) + ids = view(QRfac.p, 1:num_cut_quad_points) + condV = cond(Vq[ids, :]) + N_sampled += 5 * rd.N + + if condV < 1e8 + @warn "Conditioning of old VDM for element $e is $condV_original. " * + "After recomputing with $(length(x_sampled)) samples: $condV." + end + + iter += 1 + end + + if iter >= maxiters + @warn "Adaptive sampling of cut-cell element $e: maximum number of iterations $maxiters exceeded." + end + + # naive approach to computing quadrature weights; no guarantees of positivity + QRfac = qr(Vq', ColumnNorm()) + ids = view(QRfac.p, 1:num_cut_quad_points) + wq = Vq[ids, :]' \ b + + quadrature_error = norm(Vq[ids, :]' * wq - b) + quadrature_condition_number = sum(abs.(wq)) / sum(wq) + if quadrature_condition_number > 10 || quadrature_error > 1e3 * eps() + @warn "Quadrature error on element $e is $quadrature_error, " * + "quadrature condition number = $quadrature_condition_number. " * + "Condition number of quadrature VDM is $condV." + end + + if e==249 + if isdefined(Main, :Infiltrator) + Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__) + end + end + + view(xq.cut, :, e) .= view(x_sampled, ids) + view(yq.cut, :, e) .= view(y_sampled, ids) + view(wJq.cut, :, e) .= wq + end + + # default to non-periodic + is_periodic = (false, false) + + # compute mapping from linear element indices to Cartesian element indices + cartesian_to_linear_element_indices = compute_element_indices(region_flags) + linear_to_cartesian_element_indices = (; cut=zeros(SVector{2, Int}, num_cut_cells), + cartesian=zeros(SVector{2, Int}, num_cartesian_cells)) + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + e = cartesian_to_linear_element_indices[ex, ey] + if is_cut(region_flags[ex, ey]) + linear_to_cartesian_element_indices.cut[e] = SVector(ex, ey) + elseif is_Cartesian(region_flags[ex, ey]) + linear_to_cartesian_element_indices.cartesian[e] = SVector(ex, ey) + end + end + + cut_cell_data = (; cutcells, region_flags, wJf, + cartesian_to_linear_element_indices, + linear_to_cartesian_element_indices, + vxyz=(vx, vy), cells_per_dimension=(cells_per_dimension_x, cells_per_dimension_y), # background Cartesian grid info + ) + + # get indices of cut face nodes + face_ids(e) = (1:(num_points_per_face * cut_faces_per_cell[e])) .+ + cut_face_offsets[e] * num_points_per_face + cut_face_node_ids = [face_ids(e) for e in 1:num_cut_cells] + + if precompute_operators == true + + # precompute cut-cell operators and store them in the `md.mesh_type.cut_cell_operators` field + cut_face_nodes = cut_face_node_ids + face_interpolation_matrices = Matrix{eltype(x)}[] + lift_matrices = Matrix{eltype(x)}[] + differentiation_matrices = Tuple{Matrix{eltype(x)}, Matrix{eltype(x)}}[] + mass_matrices = Matrix{eltype(x)}[] + for (e, elem) in enumerate(physical_frame_elements) + + VDM = vandermonde(elem, rd.N, x.cut[:, e], y.cut[:, e]) + Vq, Vrq, Vsq = map(A -> A / VDM, basis(elem, rd.N, xq.cut[:,e], yq.cut[:, e])) + + M = Vq' * diagm(wJq.cut[:, e]) * Vq + Qr = Vq' * diagm(wJq.cut[:, e]) * Vrq + Qs = Vq' * diagm(wJq.cut[:, e]) * Vsq + Dx_e, Dy_e = M \ Qr, M \ Qs + + Vf = vandermonde(elem, rd.N, xf.cut[cut_face_nodes[e]], yf.cut[cut_face_nodes[e]]) / VDM + + # don't include jacobian scaling in LIFT matrix (for consistency with the Cartesian mesh) + _, w1D = quad_rule_face + num_cut_faces = length(cut_face_nodes[e]) ÷ length(w1D) + wf = vec(repeat(w1D, 1, num_cut_faces)) + + push!(lift_matrices, M \ (Vf' * diagm(wf))) + push!(face_interpolation_matrices, Vf) + push!(differentiation_matrices, (Dx_e, Dy_e)) + push!(mass_matrices, M) + end + cut_cell_operators = (; differentiation_matrices, face_interpolation_matrices, + mass_matrices, lift_matrices) + + else + + cut_cell_operators = nothing + end + + # added for consistency + J = NamedArrayPartition(cartesian=Fill(J, size(x.cartesian)), + cut=Fill(1.0, size(x.cut))) + return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_ids, objects, + cut_cell_operators, cut_cell_data, quadrature_type), + FToF, (x, y), (xf, yf), (xq, yq), wJq, + mapM, mapP, mapB, rstxyzJ, J, (nxJ, nyJ), Jf, is_periodic) + +end + +# Computes face geometric terms from a RefElemData, `quad_rule_face = (r1D, w1D)`, +# the vectors of the 1D vertex nodes `vx` and `vy`, and named tuple +# `cutcell_data is a NamedTuple containing `objects`, `region_flags`, `stop_pts``, `cutcells`. +function compute_geometric_data(rd::RefElemData{2, Quad}, quad_rule_face, + vx, vy, cutcell_data; tol=100 * eps()) + + # domain size and reference face weights + cells_per_dimension_x, cells_per_dimension_y = length(vx) - 1, length(vy) - 1 + LX, LY = (x -> x[2] - x[1]).(extrema.((vx, vy))) + + r1D, w1D = quad_rule_face + + (; objects, region_flags, cutcells, cut_faces_per_cell ) = cutcell_data + + # count number of cells and cut face nodes + num_cartesian_cells = sum(region_flags .== 0) + num_cut_cells = sum(region_flags .> 0) + nodes_per_face = length(r1D) + num_cut_face_nodes = nodes_per_face * sum(cut_faces_per_cell) + + # compute face data + xf, yf, nxJ, nyJ, Jf = ntuple(_ -> NamedArrayPartition(cartesian=zeros(rd.Nfq, num_cartesian_cells), + cut=zeros(num_cut_face_nodes)), 5) + + # the face Jacobian involves scaling between mapped and reference domain + # this is precomputed here since it's needed to compute the normals + face_ids_left_right = 1:(length(rd.rf) ÷ 2) + face_ids_top_bottom = ((length(rd.rf) ÷ 2) + 1):length(rd.rf) + Jf.cartesian[face_ids_top_bottom, :] .= LX / (cells_per_dimension_x * sum(w1D)) + Jf.cartesian[face_ids_left_right, :] .= LY / (cells_per_dimension_y * sum(w1D)) + + # compute face data + e = 1 + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + if is_Cartesian(region_flags[ex, ey]) + vx_element = SVector(vx[ex], vx[ex + 1], vx[ex], vx[ex + 1]) + vy_element = SVector(vy[ey], vy[ey], vy[ey + 1], vy[ey + 1]) + x_element, y_element = map(x -> rd.V1 * x, (vx_element, vy_element)) + mul!(view(xf.cartesian, :, e), rd.Vf, x_element) + mul!(view(yf.cartesian, :, e), rd.Vf, y_element) + view(nxJ.cartesian, :, e) .= rd.nrJ .* view(Jf.cartesian, :, e) + view(nyJ.cartesian, :, e) .= rd.nsJ .* view(Jf.cartesian, :, e) + e = e + 1 + end + end + + e = 1 + offset = 0 + for (e, curve) in enumerate(cutcells) + # map 1D quadrature points to faces + num_cut_faces = length(curve.subcurves) + fids = (1:length(r1D) * num_cut_faces) .+ offset + out = map(x->view(x, fids), (xf.cut, yf.cut, nxJ.cut, nyJ.cut, Jf.cut)) + map_points_to_cut_cell_faces!(out, r1D, curve) + offset += length(fids) + end + + # compute face points + shifting/scaling coefficients for physical frame cut elements. + physical_frame_elements = PhysicalFrame{2}[] # populate this as we iterate through cut cells + + # store cut-cell scaling/shifting coefficients + (; cut_faces_per_cell, cut_face_offsets ) = cutcell_data + num_points_per_face = length(r1D) + + e = 1 + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + if is_cut(region_flags[ex, ey]) + + # here, we evaluate a PhysicalFrame basis by shifting and scaling the + # coordinates on an element back to the reference element [-1, 1]^2. + cut_face_node_ids = (1:num_points_per_face * cut_faces_per_cell[e]) .+ + num_points_per_face * cut_face_offsets[e] + + # store face nodes (extremal) and coordinates of background Cartesian cell + physical_frame_element = + PhysicalFrame(xf.cut[cut_face_node_ids], yf.cut[cut_face_node_ids], + SVector(vx[ex], vx[ex+1]), SVector(vy[ey], vy[ey+1])) + + push!(physical_frame_elements, physical_frame_element) + + e += 1 + end + end + + # interpolation points + x, y = ntuple(_ -> NamedArrayPartition(cartesian=zeros(rd.Np, num_cartesian_cells), + cut=zeros(Np_cut(rd.N), num_cut_cells)), 2) + + # compute interpolation points on cartesian elements + e = 1 + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + if is_Cartesian(region_flags[ex, ey]) + vx_element = SVector(vx[ex], vx[ex + 1], vx[ex], vx[ex + 1]) + vy_element = SVector(vy[ey], vy[ey], vy[ey + 1], vy[ey + 1]) + x_element, y_element = map(x -> rd.V1 * x, (vx_element, vy_element)) + view(x.cartesian, :, e) .= x_element + view(y.cartesian, :, e) .= y_element + e = e + 1 + end + end + + # Compute interpolation points on cut elements + for e in eachindex(physical_frame_elements) + + physical_frame_element = physical_frame_elements[e] + + x_sampled, y_sampled = + generate_sampling_points(objects, physical_frame_element, rd, 2 * Np_cut(rd.N)) + V = vandermonde(physical_frame_element, rd.N, x_sampled, y_sampled) + + # use pivoted QR to find good interpolation points + QRfac = qr(V', ColumnNorm()) + ids = QRfac.p[1:Np_cut(rd.N)] + + # if the condition number of the VDM is really bad, then increase the + # number of sampled points. + condV_original = cond(V[ids, :]) + condV = condV_original + N_sampled = 20 * rd.N + iter, maxiters = 0, 100 + while condV > 1e8 && iter < maxiters + x_sampled, y_sampled = + generate_sampling_points(objects, physical_frame_element, rd, 2 * Np_cut(rd.N); + N_sampled = N_sampled) + V = vandermonde(physical_frame_element, rd.N, x_sampled, y_sampled) + + # use pivoted QR to find good interpolation points + QRfac = qr(V', ColumnNorm()) + ids = QRfac.p[1:Np_cut(rd.N)] + condV = cond(V[ids, :]) + N_sampled += 5 * rd.N + + if condV < 1e8 + @warn "Conditioning of old VDM for element $e is $condV_original. " * + "After recomputing with $(length(x_sampled)) samples: $condV." + end + + iter += 1 + end + if iter >= maxiters + @warn "Adaptive sampling of cut-cell element $e: maximum number of iterations $maxiters exceeded." + end + + view(x.cut, :, e) .= x_sampled[ids] + view(y.cut, :, e) .= y_sampled[ids] + end + + # volume geometric terms + rxJ_cartesian = LX / (2 * cells_per_dimension_x) + syJ_cartesian = LY / (2 * cells_per_dimension_y) + J_cartesian = (LX / cells_per_dimension_x) * (LY / cells_per_dimension_y) / 4 + + # Note: the volume Jacobian for cut elements is 1 since the "reference element" is the + # cut element itself. Similarly, geometric terms should be 1 since `basis` computes + # physical derivatives accounting for element scaling + + # Note: we use FillArrays.Fill instead of FillArrays.Ones and FillArrays.Zeros because + # we store `rxJ, sxJ, ryJ, syJ` in a single SMatrix, which assumes one homogeneous + # type for all the entries. Since FillArrays.Ones/Zeros are distinct types, using them + # results in a type instability when accessing entries of md.rstxyzJ + rxJ = NamedArrayPartition(cartesian=Fill(rxJ_cartesian, rd.Np, num_cartesian_cells), + cut=Fill(1.0, Np_cut(rd.N), num_cut_cells)) + syJ = NamedArrayPartition(cartesian=Fill(syJ_cartesian, rd.Np, num_cartesian_cells), + cut=Fill(1.0, Np_cut(rd.N), num_cut_cells)) + sxJ, ryJ = ntuple(_ -> NamedArrayPartition(cartesian=Fill(0.0, rd.Np, num_cartesian_cells), + cut=Fill(0.0, Np_cut(rd.N), num_cut_cells)), 2) + J = NamedArrayPartition(cartesian = Fill(J_cartesian, rd.Np, num_cartesian_cells), + cut = Fill(1.0, Np_cut(rd.N), num_cut_cells)) + + rstxyzJ = SMatrix{2, 2, typeof(rxJ), 4}(rxJ, sxJ, ryJ, syJ) # pack geometric terms together + + return physical_frame_elements, x, y, rstxyzJ, J, xf, yf, nxJ, nyJ, Jf +end + +""" + connect_mesh(rd, face_centroids, region_flags, cutcells; tol = 1e2 * eps()) + +Connects faces of a cut mesh to each other, returns `FToF` such that face +`f` is connected to `FToF[f]`. + +Inputs: +- rd::RefElemData +- face_centroids = (face_centroids_x, face_centroids_y), where `face_centroids_x/y` + are vectors of coordinates of face centroids +- `region_flags`, `cutcells` are return arguments from `PathIntersections.define_regions` +The keyword argument `tol` is the tolerance for matches between face centroids. +""" +function connect_mesh(rd::RefElemData, face_centroids, cutcell_data; tol = 1e2 * eps()) + + (; region_flags, cut_faces_per_cell, cut_face_offsets ) = cutcell_data + + cells_per_dimension_x, cells_per_dimension_y = size(region_flags) + num_cartesian_cells = sum(region_flags .== 0) + num_cut_faces = sum(cut_faces_per_cell) + num_total_faces = num_cut_faces + num_faces(rd.element_type) * num_cartesian_cells + + # element_indices[ex, ey] returns the global (flattened) element index into + # the arrays `xf.cartesian[:, e]` or `xf.cut[:, e]` + element_indices = compute_element_indices(region_flags) + + # compute face centroids for making face matches + face_centroids_x, face_centroids_y = face_centroids + + # To determine face-to-face matches, we work with each background Cartesian element + # and search through the 4 neighboring background Cartesian elements for a match in + # the face centroids of the current cell and the face centroids of its neighbors. + # NOTE: this works because the cut cells can only have Cartesian neighbors across + # flat-sided faces. It wouldn't work for meshes with curved interior interfaces. + + FToF = collect(1:num_total_faces) + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + + e = element_indices[ex, ey] + + # Determine face indices of current cell. The face indices are determined + # from the flattened (e.g., not ex, ey) element ordering. + # NOTE: we search for matches between all faces of `e` and all faces of + # `e_nbr` because the ordering of faces is different for cut elements + # and Cartesian elements. + if is_Cartesian(region_flags[ex, ey]) + face_ids = (1:num_faces(rd.element_type)) .+ (e-1) * num_faces(rd.element_type) + elseif is_cut(region_flags[ex, ey]) + face_ids = (1:cut_faces_per_cell[e]) .+ cut_face_offsets[e] + + # we offset by the number of Cartesian faces so we can index globally + # into the arrays `face_centroids_x`, `face_centroids_y`. + face_ids = face_ids .+ length(face_centroids_x.cartesian) + end + + if is_inside_domain(ex, ey, region_flags) + + for f in 1:4 # each Cartesian background element has 4 neighbors + + ex_nbr, ey_nbr = neighbor_across_face(f, ex, ey) + if is_inside_domain(ex_nbr, ey_nbr, region_flags) + e_nbr = element_indices[ex_nbr, ey_nbr] + + # determine face indices of neighboring cells + if is_Cartesian(region_flags[ex_nbr, ey_nbr]) + nbr_face_ids = (1:num_faces(rd.element_type)) .+ (e_nbr-1) * num_faces(rd.element_type) + elseif is_cut(region_flags[ex_nbr, ey_nbr]) + nbr_face_ids = (1:cut_faces_per_cell[e_nbr]) .+ cut_face_offsets[e_nbr] + + # we offset by the number of Cartesian faces so we can index globally + # into the arrays `face_centroids_x`, `face_centroids_y`. + nbr_face_ids = nbr_face_ids .+ length(face_centroids_x.cartesian) + end + + # check for matches in face and neighbor face centroids. + # note: we index into the global `face_centroids_x/y` container + # rather than the `.cut` or `.cartesian subarrays`. + for i in face_ids + xy = SVector(face_centroids_x[i], face_centroids_y[i]) + for j in nbr_face_ids + xy_nbr = SVector(face_centroids_x[j], face_centroids_y[j]) + if norm(xy - xy_nbr) < tol * max(1, norm(xy), norm(xy_nbr)) + FToF[i] = j + end + end + end + + end # if enbr is_inside_domain + end + end # if e is_inside_domain + end + + return FToF +end + +# returns points (xf, yf), scaled normals (nxJ, nyJ), and face Jacobian (Jf) +# for a curve returned from PathIntersections. +# `out` should hold `xf, yf, nxJ, nyJ, Jf = ntuple(_ -> similar(points, (length(points), num_faces)), 5)` +function map_points_to_cut_cell_faces(points, curve) + num_faces = length(curve.subcurves) + out = ntuple(_ -> similar(points, (length(points) * num_faces)), 5) + return map_points_to_cut_cell_faces!(out, points, curve) +end + +function map_points_to_cut_cell_faces!(out, points, curve) + xf, yf, nxJ, nyJ, Jf = out + stop_points = curve.stop_pts + r1D = points + for f in 1:length(stop_points)-1 + for i in eachindex(r1D) + fid = i + (f-1) * length(r1D) + + s = map_to_interval(r1D[i], stop_points[f], stop_points[f+1]) + + # compute tangent vector at a node, face Jacobian, and normals + x_node, y_node = curve(s) + xf[fid], yf[fid] = x_node, y_node + tangent_vector = PathIntersections.ForwardDiff.derivative(curve, s) + + # the face Jacobian involves scaling between mapped and reference face + # reference face = [-1, 1] + scaling = (stop_points[f+1] - stop_points[f]) / 2 + Jf[fid] = norm(tangent_vector) * scaling + + # we have to flip the sign to get the outward normal. + # note: we compute the scaled normal nxJ for consistency with other meshes. + normal_node = SVector{2}(tangent_vector[2], -tangent_vector[1]) + nxJ[fid], nyJ[fid] = (-normal_node / norm(normal_node)) .* Jf[fid] + end + end + return vec.((xf, yf, nxJ, nyJ, Jf)) +end + +# generates at least Np_target sampling points within a cut cell defined by `curve` +# returns both x_sampled, y_sampled (physical points inside the cut cell), as well as +# r_sampled, y_sampled (reference points which correspond to x_sampled, y_sampled). +function generate_sampling_points(objects, elem, rd, Np_target; N_sampled = 4 * rd.N) + + r_sampled, s_sampled = equi_nodes(rd.element_type, N_sampled) # oversampled nodes + + # map sampled points to the background Cartesian cell + x_sampled, y_sampled = map_nodes_to_background_cell(elem, r_sampled, s_sampled) + is_in_domain = fill(true, length(x_sampled)) + for (index, point) in enumerate(zip(x_sampled, y_sampled)) + is_in_domain[index] = !any(map(obj -> PathIntersections.is_contained(obj, point), objects)) + end + + # increase number of background points until we are left with `Np_target` sampling points + while sum(is_in_domain) < Np_target + + N_sampled *= 2 # double degree of sampling + r_sampled, s_sampled = equi_nodes(rd.element_type, N_sampled) # oversampled nodes + x_sampled, y_sampled = map_nodes_to_background_cell(elem, r_sampled, s_sampled) + + # check if all the points are in all the objects + is_in_domain = fill(true, length(x_sampled)) + for (index, point) in enumerate(zip(x_sampled, y_sampled)) + is_in_domain[index] = !any(map(obj -> PathIntersections.is_contained(obj, point), objects)) + end + end + + ids_in_element = findall(is_in_domain) + + return x_sampled[ids_in_element], y_sampled[ids_in_element] +end + +function count_cut_faces(cutcells) + num_cut_faces = zeros(Int, length(cutcells)) + for e in eachindex(cutcells) + curve = cutcells[e] + stop_points = curve.stop_pts + num_cut_faces[e] = length(stop_points) - 1 + end + return num_cut_faces +end + +function compute_face_centroids(rd::RefElemData, xf, yf, cutcell_data) + + (; region_flags, cut_faces_per_cell, cut_face_offsets ) = cutcell_data + num_cut_cells = length(cut_faces_per_cell) + num_cartesian_cells = sum(region_flags .== 0) + num_cut_faces = sum(cut_faces_per_cell) + + num_points_per_face = length(rd.rf) ÷ num_faces(rd.element_type) + + face_centroids_x = NamedArrayPartition(cartesian=zeros(num_faces(rd.element_type), num_cartesian_cells), + cut=zeros(num_cut_faces)) + face_centroids_y = similar(face_centroids_x) + + for e in 1:num_cartesian_cells + xf_element = reshape(view(xf.cartesian, :, e), num_points_per_face, num_faces(rd.element_type)) + yf_element = reshape(view(yf.cartesian, :, e), num_points_per_face, num_faces(rd.element_type)) + view(face_centroids_x.cartesian, :, e) .= vec(sum(xf_element, dims=1) / num_points_per_face) + view(face_centroids_y.cartesian, :, e) .= vec(sum(yf_element, dims=1) / num_points_per_face) + end + + for e in 1:num_cut_cells + face_node_ids = (1:(num_points_per_face * cut_faces_per_cell[e])) .+ cut_face_offsets[e] * num_points_per_face + xf_element = reshape(view(xf.cut, face_node_ids), num_points_per_face, cut_faces_per_cell[e]) + yf_element = reshape(view(yf.cut, face_node_ids), num_points_per_face, cut_faces_per_cell[e]) + + face_ids = (1:cut_faces_per_cell[e]) .+ cut_face_offsets[e] + view(face_centroids_x.cut, face_ids) .= vec(sum(xf_element, dims=1) / num_points_per_face) + view(face_centroids_y.cut, face_ids) .= vec(sum(yf_element, dims=1) / num_points_per_face) + end + + return face_centroids_x, face_centroids_y +end \ No newline at end of file diff --git a/src/explicit_timestep_utils.jl b/src/explicit_timestep_utils.jl index 77382a77..be915750 100644 --- a/src/explicit_timestep_utils.jl +++ b/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() diff --git a/src/geometric_functions.jl b/src/geometric_functions.jl index 2982b339..00e8a251 100644 --- a/src/geometric_functions.jl +++ b/src/geometric_functions.jl @@ -77,7 +77,8 @@ function estimate_h(rd::RefElemData{DIM}, md::MeshData{DIM}) where {DIM} h_e = estimate_h(e, rd, md) hmin = min(hmin, h_e) end - return hmin * compute_domain_size(rd, md)^(1/DIM) + domain_size = sum(rd.M * md.J) + return hmin * domain_size^(1/DIM) end function estimate_h(e, rd::RefElemData, md::MeshData) @@ -106,5 +107,4 @@ face_scaling(rd, f) = 1.0 face_scaling(rd::RefElemData{2, Tri}, f) = f == 3 ? sqrt(2) : 1.0 # Jf incorporates length of long triangle edge face_scaling(rd::RefElemData{3, Tet}, f) = f == 2 ? sqrt(3) : 1.0 # Jf incorporates area of larger triangle face face_scaling(rd::RefElemData{3, Wedge}, f) = f == 2 ? sqrt(2) : 1.0 # Jf incorporates area of larger triangle face -compute_domain_size(rd::RefElemData, md::MeshData) = sum(rd.M * md.J) diff --git a/src/hybrid_meshes.jl b/src/hybrid_meshes.jl index 4bb2e26a..51559bac 100644 --- a/src/hybrid_meshes.jl +++ b/src/hybrid_meshes.jl @@ -142,9 +142,11 @@ end typename(x) = typeof(x).name.name -function RefElemData(element_types::NTuple{N, Union{Tri, Quad}}, args...; kwargs...) where {N} - - rds = NamedTuple((typename(elem) => RefElemData(elem, args...; kwargs...) for elem in element_types)) +function RefElemData(element_types::NTuple{NT, Union{Tri, Quad}}, N::Int; + quad_rule_face=gauss_quad(0, 0, N), kwargs...) where {NT} + + rds = NamedTuple((typename(elem) => + RefElemData(elem, N; quad_rule_face, kwargs...) for elem in element_types)) # check if number of face nodes is the same num_face_nodes = length.(getproperty.(values(rds), :rf)) .÷ num_faces.(getproperty.(values(rds), :element_type)) diff --git a/src/low_order_sbp.jl b/src/low_order_sbp.jl new file mode 100644 index 00000000..7e0d1b2f --- /dev/null +++ b/src/low_order_sbp.jl @@ -0,0 +1,224 @@ +""" + function sparse_low_order_SBP_operators(rd; factor=1.01) + +Constructs sparse low order SBP operators given a `RefElemData`. +Returns operators `Qrst..., E ≈ Vf * Pq` that satisfy a generalized +summation-by-parts (GSBP) property: + + `Q_i + Q_i^T = E' * B_i * E` + +`factor` is a scaling which determines how close a node must be to +another node to be considered a neighbor. +""" +function sparse_low_order_SBP_operators(rd::RefElemData{NDIMS}; factor=1.01) where {NDIMS} + (; Pq, Vf, rstq, wf, nrstJ) = rd + + # if element is a simplex, convert nodes to an equilateral triangle for symmetry + rstq = map_nodes_to_symmetric_element(rd.element_type, rstq...) + + # build distance and adjacency matrices + D = [norm(SVector{NDIMS}(getindex.(rstq, i)) - SVector{NDIMS}(getindex.(rstq, j))) + for i in eachindex(first(rstq)), j in eachindex(first(rstq))] + A = zeros(Int, length(first(rstq)), length(first(rstq))) + for i in axes(D, 1) + # heuristic cutoff - we take NDIMS + 1 neighbors, but the smallest distance = 0, + # so we need to access the first NDIMS + 2 sorted distance entries. + dist = sort(view(D, i, :))[NDIMS + 2] + for j in findall(view(D, i, :) .< factor * dist) + A[i, j] = 1 + end + end + A = (A + A' .> 0) # symmetrize + L = Diagonal(vec(sum(A, dims=2))) - A + sorted_eigvals = sort(eigvals(L)) + + # check that there's only one zero null vector + if sorted_eigvals[2] < 10 * size(A, 1) * eps() + error("Warning: the connectivity between nodes yields a graph with " * + "more than one connected component. Try increasing the value of `factor`.") + end + + # note: this part doesn't do anything for a nodal SBP operator. In that case, E_dense = E = Vf + # and E should just reduce to Vf, which is a face extraction operator. + E_dense = Vf * Pq + E = zeros(size(E_dense)) + for i in axes(E, 1) + # find all j such that E[i,j] ≥ 0.5, e.g., points which positively contribute to at least half of the + # interpolation. These seem to be associated with volume points "j" that are close to face point "i". + ids = findall(E_dense[i, :] .>= 0.5) + E[i, ids] .= E_dense[i, ids] ./ sum(E_dense[i, ids]) # normalize so sum(E, dims=2) = [1, 1, ...] still. + end + Brst = (nJ -> diagm(wf .* nJ)).(nrstJ) + + # compute potential + e = ones(size(L, 2)) + right_hand_sides = map(B -> 0.5 * sum(E' * B, dims=2), Brst) + psi_augmented = map(b -> [L e; e' 0] \ [b; 0], right_hand_sides) + psi = map(x -> x[1:end-1], psi_augmented) + + # construct sparse skew part + function construct_skew_matrix_from_potential(ψ) + S = zeros(length(ψ), length(ψ)) + for i in axes(S, 1), j in axes(S, 2) + if A[i,j] > 0 + S[i,j] = ψ[j] - ψ[i] + end + end + return S + end + + Srst = construct_skew_matrix_from_potential.(psi) + Qrst = map((S, B) -> S + 0.5 * E' * B * E, Srst, Brst) + return sparse.(Qrst), sparse(E) +end + +function sparse_low_order_SBP_1D_operators(rd::RefElemData) + E = zeros(2, rd.N+1) + E[1, 1] = 1 + E[2, end] = 1 + + # create volume operators + Q = diagm(1 => ones(rd.N), -1 => -ones(rd.N)) + Q[1,1] = -1 + Q[end,end] = 1 + Q = 0.5 * Q + + return (sparse(Q),), sparse(E) +end + +sparse_low_order_SBP_operators(rd::RefElemData{1, Line, <:Union{<:SBP, <:Polynomial{Gauss}}}) = + sparse_low_order_SBP_1D_operators(rd) + +function diagonal_1D_mass_matrix(N, ::SBP) + _, w1D = gauss_lobatto_quad(0, 0, N) + return Diagonal(w1D) +end + +function diagonal_1D_mass_matrix(N, ::Polynomial{Gauss}) + _, w1D = gauss_quad(0, 0, N) + return Diagonal(w1D) +end + +function sparse_low_order_SBP_operators(rd::RefElemData{2, Quad, <:Union{<:SBP, <:Polynomial{Gauss}}}) + (Q1D,), E1D = sparse_low_order_SBP_1D_operators(rd) + + # permute face node ordering for the first 2 faces + ids = reshape(1:(rd.N+1) * 2, :, 2) + Er = zeros((rd.N+1) * 2, rd.Np) + Er[vec(ids'), :] .= kron(I(rd.N+1), E1D) + Es = kron(E1D, I(rd.N+1)) + E = vcat(Er, Es) + + M1D = diagonal_1D_mass_matrix(rd.N, rd.approximation_type) + Qr = kron(M1D, Q1D) + Qs = kron(Q1D, M1D) + + return sparse.((Qr, Qs)), sparse(E) +end + +function sparse_low_order_SBP_operators(rd::RefElemData{3, Hex, <:Union{<:SBP, <:Polynomial{Gauss}}}) + (Q1D,), E1D = sparse_low_order_SBP_1D_operators(rd) + + # permute face nodes for the faces in the ±r and ±s directions + ids = reshape(1:(rd.N+1)^2 * 2, rd.N+1, :, 2) + Er, Es, Et = ntuple(_ -> zeros((rd.N+1)^2 * 2, rd.Np), 3) + Er[vec(permutedims(ids, [2, 3, 1])), :] .= kron(I(rd.N+1), E1D, I(rd.N+1)) + Es[vec(permutedims(ids, [3, 2, 1])), :] .= kron(I(rd.N+1), I(rd.N+1), E1D) + Et .= kron(E1D, I(rd.N+1), I(rd.N+1)) + + # create boundary extraction matrix + E = vcat(Er, Es, Et) + + M1D = diagonal_1D_mass_matrix(rd.N, rd.approximation_type) + Qr = kron(M1D, Q1D, M1D) + Qs = kron(M1D, M1D, Q1D) + Qt = kron(Q1D, M1D, M1D) + + return sparse.((Qr, Qs, Qt)), sparse(E) +end + +# builds subcell operators D * R such that for r such that sum(r) = 0, +# D * Diagonal(θ) * R * r is also diagonal for any choice of θ. This is +# useful in subcell limiting (see, for example +# https://doi.org/10.1016/j.compfluid.2022.105627 for a 1D example) +function subcell_limiting_operators(Qr::AbstractMatrix; tol = 100 * eps()) + Qr = sparse(Qr) + Sr = droptol!(Qr - Qr', tol) + Br = droptol!(Qr + Qr', tol) + + num_interior_fluxes = nnz(Sr) ÷ 2 + num_boundary_indices = nnz(Br) + num_fluxes = num_interior_fluxes + num_boundary_indices + + interior_indices = findall(Sr .> tol) + boundary_indices = findall(abs.(Br) .> tol) + + matrix_indices = zeros(Int, size(Sr)) + for (i, ij) in enumerate(boundary_indices) + matrix_indices[ij] = i + matrix_indices[ij.I[2], ij.I[1]] = i + end + + for (i, ij) in enumerate(interior_indices) + matrix_indices[ij] = i + length(boundary_indices) + matrix_indices[ij.I[2], ij.I[1]] = i + length(boundary_indices) + end + + D = zeros(size(Qr, 2), num_fluxes) + for i in axes(matrix_indices, 1), j in axes(matrix_indices, 2) + if matrix_indices[i, j] !== 0 + D[i, matrix_indices[i, j]] = Qr[i, j] + end + end + + d = vec(ones(size(D, 1))' * D) + ids = findall(@. abs(d) > 1e2 * eps()) + + Z = nullspace(D)[ids, :] + A = pinv(D)[ids,:] + + # compute R via linear algebra + m, n = size(A) + z = size(Z, 2) + C = hcat(kron(I(n), Z), kron(ones(n), I(m))) + sol = pinv(C) * -vec(A) + X = reshape(sol[1:n*z], (z, n)) + R = pinv(D) + nullspace(D) * X + + # check if R satisfies our pseudoinverse condition + @assert norm(D * R - I(size(D,1))) < tol * length(D) + + return D, R +end + +""" + Δrst, Rrst = subcell_limiting_operators(rd::RefElemData) + +Returns tuples of subcell limiting operators Drst = (Δr, Δs, ...) and R = (Rr, Rs, ...) +such that for r where sum(r) = 0, sum(D * Diagonal(θ) * R * r) = 0 for any choice of θ. +These operators are useful for conservative subcell limiting techniques (see +https://doi.org/10.1016/j.compfluid.2022.105627 for an example of such an approach on +tensor product elements). + +Sparse SBP operators used in an intermediate step when buidling these subcell +limiting operators; by default, these operators are constructed using +`sparse_low_order_SBP_operators`. To construct subcell limiting operators for a +general SBP operator, one can use the following: + + Δ, R = subcell_limiting_operators(Q::AbstractMatrix; tol = 100 * eps()) +""" +function subcell_limiting_operators(rd::RefElemData) + Qrst, _ = sparse_low_order_SBP_operators(rd) + Δrst, Rrst = subcell_limiting_operators.(Qrst) + return zip(Δrst, Rrst) +end + +# specialize for single dimension to return tuples of operators +function subcell_limiting_operators(rd::RefElemData{1}) + Qrst, _ = sparse_low_order_SBP_operators(rd) + Δrst, Rrst = first(subcell_limiting_operators.(Qrst)) + return (Δrst,), (Rrst,) +end + +# TODO: specialize for quadrilaterals and hex? + diff --git a/src/mesh/gmsh_utilities.jl b/src/mesh/gmsh_utilities.jl index db958d7e..936b2668 100644 --- a/src/mesh/gmsh_utilities.jl +++ b/src/mesh/gmsh_utilities.jl @@ -103,7 +103,7 @@ Notes: the version 4 format has a more detailed block data format this leads to more complicated parser. """ function read_Gmsh_2D_v4(filename::String, options::MeshImportOptions) - @unpack grouping, remap_group_name = options + (; grouping, remap_group_name) = options if !isfile(filename) throw(ArgumentError("file $filename does not exist")) @@ -205,7 +205,6 @@ function read_Gmsh_2D_v4(filename::String, options::MeshImportOptions) block_line_start = block_line_start + num_elem_in_block + 1 end - EToV = EToV[:, vec([1 3 2])] # permute for Gmsh ordering EToV = correct_negative_Jacobians!((VX, VY), EToV) if grouping @@ -320,7 +319,6 @@ function read_Gmsh_2D_v2(filename::String) end end - EToV = EToV[:, vec([1 3 2])] # permute for Gmsh ordering EToV = correct_negative_Jacobians!((VX, VY), EToV) return (VX, VY), EToV @@ -331,7 +329,7 @@ end # Computes the area of a triangle given `tri`, which is a tuple of three points (vectors), # using the [Shoelace_formula](https://en.wikipedia.org/wiki/Shoelace_formula). function compute_triangle_area(tri) - B, A, C = tri + A, B, C = tri return 0.5 * (A[1] * (B[2] - C[2]) + B[1] * (C[2] - A[2]) + C[1] * (A[2] - B[2])) end diff --git a/src/mesh/simple_meshes.jl b/src/mesh/simple_meshes.jl index ffa92075..36aef802 100644 --- a/src/mesh/simple_meshes.jl +++ b/src/mesh/simple_meshes.jl @@ -29,8 +29,8 @@ function uniform_mesh(elem::Tri, Kx, Ky) id2 = id(ex + 1, ey) id3 = id(ex + 1, ey + 1) id4 = id(ex, ey + 1) - EToV[2*sk-1, :] = [id1 id3 id2] - EToV[2*sk, :] = [id3 id1 id4] + EToV[2*sk-1, :] = [id1 id2 id3] + EToV[2*sk, :] = [id3 id4 id1] sk += 1 end end @@ -123,15 +123,15 @@ function uniform_mesh(elem::Hex, Nx, Ny, Nz) end -function uniform_mesh(elem::Wedge, Kx, Ky, Kz) +function uniform_mesh(::Wedge, Kx, Ky, Kz) (VY, VX, VZ) = meshgrid(LinRange(-1, 1, Ky + 1), LinRange(-1, 1, Kx + 1), LinRange(-1, 1, Kz + 1)) sk = 1 - shift = (Kx+1)*(Ky+1) + shift = (Kx + 1) * (Ky + 1) + id(ex, ey, ez) = ex + (ey - 1) * (Kx + 1) + (ez - 1) * (shift) # index function EToV = zeros(Int, 2 * Kx * Ky * Kz, 6) for ey = 1:Ky for ex = 1:Kx for ez = 1:Kz - id(ex, ey, ez) = ex + (ey - 1) * (Kx + 1) + (ez - 1) * (shift)# index function id1 = id(ex, ey, ez) id2 = id(ex + 1, ey, ez) id3 = id(ex, ey + 1, ez) @@ -140,14 +140,13 @@ function uniform_mesh(elem::Wedge, Kx, Ky, Kz) id6 = id(ex + 1, ey, ez + 1) id7 = id(ex, ey + 1, ez + 1) id8 = id(ex + 1, ey + 1, ez + 1) - EToV[2*sk-1, :] = [id1 id5 id4 id8 id2 id6] - EToV[2*sk, :] = [id1 id5 id3 id7 id4 id8] + EToV[2 * sk - 1, :] = [id1 id2 id4 id5 id6 id8] + EToV[2 * sk, :] = [id1 id4 id3 id5 id8 id7] sk += 1 end end end - EToV .= EToV[:, SVector(1, 3, 5, 2, 4, 6)] - return (VX[:], VY[:], VZ[:]), EToV + return vec.((VX, VY, VZ)), EToV end # split each cube into 6 pyramids. Pyramid 1: diff --git a/src/mesh/triangulate_utils.jl b/src/mesh/triangulate_utils.jl index f189af7f..918c16e9 100644 --- a/src/mesh/triangulate_utils.jl +++ b/src/mesh/triangulate_utils.jl @@ -19,7 +19,6 @@ Computes `VX`,`VY`,`EToV` from a `TriangulateIO` object. function triangulateIO_to_VXYEToV(triout::TriangulateIO) VX,VY = (triout.pointlist[i,:] for i = 1:size(triout.pointlist,1)) EToV = permutedims(triout.trianglelist) - Base.swapcols!(EToV,2,3) # to match MeshData ordering return (VX, VY), Matrix{Int}(EToV) end @@ -37,8 +36,8 @@ function get_boundary_face_labels(triout::TriangulateIO, rd::RefElemData{2, Tri} for (f,boundary_face) in enumerate(boundary_faces) element = (boundary_face - 1) ÷ rd.Nfaces + 1 face = (boundary_face - 1) % rd.Nfaces + 1 - vertices_on_face = sort(md.EToV[element,rd.fv[face]]) - tag_id = findfirst(c->view(segmentlist,:,c)==vertices_on_face,axes(segmentlist,2)) + vertices_on_face = sort(md.EToV[element, rd.fv[face]]) + tag_id = findfirst(c -> view(segmentlist,:,c) == vertices_on_face,axes(segmentlist, 2)) boundary_face_tags[f] = triout.segmentmarkerlist[tag_id] end return boundary_face_tags, boundary_faces diff --git a/src/named_array_partition.jl b/src/named_array_partition.jl deleted file mode 100644 index 32e4d947..00000000 --- a/src/named_array_partition.jl +++ /dev/null @@ -1,117 +0,0 @@ - -using RecursiveArrayTools: RecursiveArrayTools, ArrayPartition, npartitions, unpack - -""" - NamedArrayPartition(; kwargs...) - NamedArrayPartition(x::NamedTuple) - -Similar to an `ArrayPartition` but the individual arrays can be accessed via the -constructor-specified names. However, unlike `ArrayPartition`, each individual array -must have the same element type. -""" -struct NamedArrayPartition{T, A<:ArrayPartition{T}, NT<:NamedTuple} <: AbstractVector{T} - array_partition::A - names_to_indices::NT -end -NamedArrayPartition(; kwargs...) = NamedArrayPartition(NamedTuple(kwargs)) -function NamedArrayPartition(x::NamedTuple) - names_to_indices = NamedTuple(Pair(symbol, index) for (index, symbol) in enumerate(keys(x))) - - # enforce homogeneity of eltypes - @assert all(eltype.(values(x)) .== eltype(first(x))) - T = eltype(first(x)) - S = typeof(values(x)) - return NamedArrayPartition(ArrayPartition{T, S}(values(x)), names_to_indices) -end - -# note that overloading `getproperty` means we cannot access `NamedArrayPartition` -# fields except through `getfield` and accessor functions. -ArrayPartition(x::NamedArrayPartition) = getfield(x, :array_partition) - -Base.Array(x::NamedArrayPartition) = Array(ArrayPartition(x)) - -Base.zero(x::NamedArrayPartition{T, S, TN}) where {T, S, TN} = - NamedArrayPartition{T, S, TN}(zero(ArrayPartition(x)), getfield(x, :names_to_indices)) -Base.zero(A::NamedArrayPartition, dims::NTuple{N, Int}) where {N} = zero(A) # ignore dims since named array partitions are vectors - - -Base.propertynames(x::NamedArrayPartition) = propertynames(getfield(x, :names_to_indices)) -Base.getproperty(x::NamedArrayPartition, s::Symbol) = - getindex(ArrayPartition(x).x, getproperty(getfield(x, :names_to_indices), s)) - -# this enables x.s = some_array. -@inline function Base.setproperty!(x::NamedArrayPartition, s::Symbol, v) - index = getproperty(getfield(x, :names_to_indices), s) - ArrayPartition(x).x[index] .= v -end - -# print out NamedArrayPartition as a NamedTuple -Base.summary(x::NamedArrayPartition) = string(typeof(x), " with arrays:") -Base.show(io::IO, m::MIME"text/plain", x::NamedArrayPartition) = - show(io, m, NamedTuple(Pair.(keys(getfield(x, :names_to_indices)), ArrayPartition(x).x))) - -Base.size(x::NamedArrayPartition) = size(ArrayPartition(x)) -Base.length(x::NamedArrayPartition) = length(ArrayPartition(x)) -Base.getindex(x::NamedArrayPartition, args...) = getindex(ArrayPartition(x), args...) - -Base.setindex!(x::NamedArrayPartition, args...) = setindex!(ArrayPartition(x), args...) -Base.map(f, x::NamedArrayPartition) = NamedArrayPartition(map(f, ArrayPartition(x)), getfield(x, :names_to_indices)) -Base.mapreduce(f, op, x::NamedArrayPartition) = mapreduce(f, op, ArrayPartition(x)) -# Base.filter(f, x::NamedArrayPartition) = filter(f, ArrayPartition(x)) - -Base.similar(x::NamedArrayPartition{T, S, NT}) where {T, S, NT} = - NamedArrayPartition{T, S, NT}(similar(ArrayPartition(x)), getfield(x, :names_to_indices)) - -# broadcasting -Base.BroadcastStyle(::Type{<:NamedArrayPartition}) = Broadcast.ArrayStyle{NamedArrayPartition}() -function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{NamedArrayPartition}}, - ::Type{ElType}) where {ElType} - x = find_NamedArrayPartition(bc) - return NamedArrayPartition(similar(ArrayPartition(x)), getfield(x, :names_to_indices)) -end - -# when broadcasting with ArrayPartition + another array type, the output is the other array tupe -Base.BroadcastStyle(::Broadcast.ArrayStyle{NamedArrayPartition}, ::Broadcast.DefaultArrayStyle{1}) = - Broadcast.DefaultArrayStyle{1}() - -# hook into ArrayPartition broadcasting routines -@inline RecursiveArrayTools.npartitions(x::NamedArrayPartition) = npartitions(ArrayPartition(x)) -@inline RecursiveArrayTools.unpack(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{NamedArrayPartition}}, i) = - Broadcast.Broadcasted(bc.f, RecursiveArrayTools.unpack_args(i, bc.args)) -@inline RecursiveArrayTools.unpack(x::NamedArrayPartition, i) = unpack(ArrayPartition(x), i) - -Base.copy(A::NamedArrayPartition{T,S,NT}) where {T,S,NT} = - NamedArrayPartition{T,S,NT}(copy(ArrayPartition(A)), getfield(A, :names_to_indices)) - -@inline NamedArrayPartition(f::F, N, names_to_indices) where F<:Function = - NamedArrayPartition(ArrayPartition(ntuple(f, Val(N))), names_to_indices) - -@inline function Base.copy(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{NamedArrayPartition}}) - N = npartitions(bc) - @inline function f(i) - copy(unpack(bc, i)) - end - x = find_NamedArrayPartition(bc) - NamedArrayPartition(f, N, getfield(x, :names_to_indices)) -end - -@inline function Base.copyto!(dest::NamedArrayPartition, - bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{NamedArrayPartition}}) - N = npartitions(dest, bc) - @inline function f(i) - copyto!(ArrayPartition(dest).x[i], unpack(bc, i)) - end - ntuple(f, Val(N)) - return dest -end - -# `x = find_NamedArrayPartition(x)` returns the first `NamedArrayPartition` among broadcast arguments. -find_NamedArrayPartition(bc::Base.Broadcast.Broadcasted) = find_NamedArrayPartition(bc.args) -find_NamedArrayPartition(args::Tuple) = - find_NamedArrayPartition(find_NamedArrayPartition(args[1]), Base.tail(args)) -find_NamedArrayPartition(x) = x -find_NamedArrayPartition(::Tuple{}) = nothing -find_NamedArrayPartition(x::NamedArrayPartition, rest) = x -find_NamedArrayPartition(::Any, rest) = find_NamedArrayPartition(rest) - - diff --git a/src/physical_frame_basis.jl b/src/physical_frame_basis.jl index e81870aa..dcff6a46 100644 --- a/src/physical_frame_basis.jl +++ b/src/physical_frame_basis.jl @@ -63,11 +63,13 @@ function NodesAndModes.basis(elem::PhysicalFrame{2}, N, x, y) V, Vr, Vs = ntuple(x->zeros(length(r), Np), 3) for j = 0:N P_j = jacobiP(s, 0, 0, j) + dP_j = grad_jacobiP(s, 0, 0, j) for i = 0:N-j P_i = jacobiP(r, 0, 0, i) - V[:, sk] = P_i .* P_j - Vr[:, sk] = grad_jacobiP(r, 0, 0, i) .* P_j * scaling[1] - Vs[:, sk] = P_i .* grad_jacobiP(s, 0, 0, j) * scaling[2] + dP_i = grad_jacobiP(r, 0, 0, i) + @. V[:, sk] = P_i * P_j + @. Vr[:, sk] = dP_i * P_j * scaling[1] + @. Vs[:, sk] = P_i * dP_j * scaling[2] sk += 1 end end @@ -144,4 +146,53 @@ function map_nodes_to_background_cell(elem::PhysicalFrame{2}, r, s) x = @. 0.5 * (1 + r) * dx + vx[1] y = @. 0.5 * (1 + s) * dy + vy[1] return x, y +end + +function triangulate_points(coordinates::AbstractMatrix) + triin=Triangulate.TriangulateIO() + triin.pointlist = coordinates + triout, _ = triangulate("Q", triin) + VX, VY = (triout.pointlist[i,:] for i = 1:size(triout.pointlist,1)) + EToV = permutedims(triout.trianglelist) + return (VX, VY), EToV +end + +""" + caratheodory_pruning_qr(V, w_in) + +This performs Caratheodory pruning using a naive QR-based algorithm. +Returns (w, inds), where `inds` denotes sub-selected indices for a +reduced quadrature rule, and `w` is a vector of reduced positive weights. + +The original Matlab code this was based on was authored by Akil Narayan. +""" +function caratheodory_pruning_qr(V, w_in) + + if length(w_in) <= size(V, 2) + return w_in, eachindex(w_in) + end + w = copy(w_in) + M, N = size(V) + inds = collect(1:M) + m = M-N + Q, _ = qr(V) + for _ in 1:m + kvec = view(Q, :, size(Q, 2)) + + # for subtracting the kernel vector + idp = findall(@. kvec > 0) + alphap, k0p = findmin(view(w, inds[idp]) ./ view(kvec, idp)) + k0p = idp[k0p] + + # for adding the kernel vector + idn = findall(@. kvec < 0) + alphan, k0n = findmax(view(w, inds[idn]) ./ view(kvec, idn)) + k0n = idn[k0n] + + alpha, k0 = abs(alphan) < abs(alphap) ? (alphan, k0n) : (alphap, k0p) + @. w[inds] = w[inds] - alpha * kvec + deleteat!(inds, k0) + Q, _ = qr(V[inds, :]) + end + return w, inds end \ No newline at end of file diff --git a/src/ref_elem_utils.jl b/src/ref_elem_utils.jl index 1d2bb910..1fa75801 100644 --- a/src/ref_elem_utils.jl +++ b/src/ref_elem_utils.jl @@ -36,6 +36,81 @@ function map_face_nodes(::Tet, face_nodes...) return rf, sf, tf end +function init_face_data(elem::Tri; quad_rule_face = gauss_quad(0,0,N)) + r1D, w1D = quad_rule_face + e = ones(size(r1D)) + z = zeros(size(r1D)) + rf, sf = map_face_nodes(elem, r1D) + wf = vec(repeat(w1D, 3, 1)); + nrJ = [z; e; -e] + nsJ = [-e; e; z] + return rf, sf, wf, nrJ, nsJ +end + +function init_face_data(elem::Quad; quad_rule_face=gauss_quad(0, 0, N)) + Nfaces = 4 + r1D, w1D = quad_rule_face + e = ones(size(r1D)) + z = zeros(size(r1D)) + rf, sf = map_face_nodes(elem, r1D) + wf = vec(repeat(w1D, Nfaces, 1)) + nrJ = [-e; e; z; z] + nsJ = [z; z; -e; e] + + return rf, sf, wf, nrJ, nsJ +end + +function init_face_data(elem::Hex; quad_rule_face=quad_nodes(Quad(), N)) + rquad, squad, wquad = quad_rule_face + e = ones(size(rquad)) + zz = zeros(size(rquad)) + rf, sf, tf = map_face_nodes(elem, rquad, squad) + Nfaces = 6 + wf = vec(repeat(wquad, Nfaces, 1)); + nrJ = [-e; e; zz; zz; zz; zz] + nsJ = [zz; zz; -e; e; zz; zz] + ntJ = [zz; zz; zz; zz; -e; e] + return rf, sf, tf, wf, nrJ, nsJ, ntJ +end + +function init_face_data(elem::Tet; quad_rule_face=quad_nodes(Tri(), N)) + rquad, squad, wquad = quad_rule_face + e = ones(size(rquad)) + zz = zeros(size(rquad)) + rf, sf, tf = map_face_nodes(elem, rquad, squad) + Nfaces = 4 + wf = vec(repeat(wquad, Nfaces, 1)); + nrJ = [zz; e; -e; zz] + nsJ = [-e; e; zz; zz] + ntJ = [zz; e; zz; -e] + return rf, sf, tf, wf, nrJ, nsJ, ntJ +end + +# default to doing nothing +map_nodes_to_symmetric_element(element_type, rst...) = rst + +# for triangles, map to an equilateral triangle +function map_nodes_to_symmetric_element(::Tri, r, s) + # biunit right triangular vertices + v1, v2, v3 = SVector{2}.(zip(nodes(Tri(), 1)...)) + + denom = (v2[2] - v3[2]) * (v1[1] - v3[1]) + (v3[1] - v2[1]) * (v1[2] - v3[2]) + L1 = @. ((v2[2] - v3[2]) * (r - v3[1]) + (v3[1] - v2[1]) * (s - v3[2])) / denom + L2 = @. ((v3[2] - v1[2]) * (r - v3[1]) + (v1[1] - v3[1]) * (s - v3[2])) / denom + L3 = @. 1 - L1 - L2 + + # equilateral vertices + v1 = SVector{2}(2 * [-.5, -sqrt(3) / 6]) + v2 = SVector{2}(2 * [.5, -sqrt(3)/6]) + v3 = SVector{2}(2 * [0, sqrt(3)/3]) + + x = @. v1[1] * L1 + v2[1] * L2 + v3[1] * L3 + y = @. v1[2] * L1 + v2[2] * L2 + v3[2] * L3 + + return x, y +end + + """ function inverse_trace_constant(rd::RefElemData) diff --git a/test/MeshData_tests.jl b/test/MeshData_tests.jl index 39f86075..d939f90d 100644 --- a/test/MeshData_tests.jl +++ b/test/MeshData_tests.jl @@ -208,11 +208,11 @@ end end +# Note: we test wedge and pyr elements in a separate file approx_elem_types_to_test = [(Polynomial(), Hex()), (SBP(), Hex()), (Polynomial(), Tet()), - (Polynomial(), Wedge()), - (Polynomial(), Pyr())] + ] @testset "3D MeshData tests" begin @testset "$approximation_type $element_type MeshData initialization" for (approximation_type, element_type) in approx_elem_types_to_test tol = 5e2*eps() @@ -288,90 +288,4 @@ approx_elem_types_to_test = [(Polynomial(), Hex()), @test isempty(md.mapB) end - - @testset "TensorProductWedge MeshData" begin - element_type = Wedge() - tol = 5e2*eps() - @testset "Degree $tri_grad triangle" for tri_grad = [2, 3] - @testset "Degree $line_grad line" for line_grad = [2, 3] - line = RefElemData(Line(), line_grad) - tri = RefElemData(Tri(), tri_grad) - tensor = TensorProductWedge(tri, line) - - rd = RefElemData(element_type, tensor) - K1D = 2 - md = MeshData(uniform_mesh(element_type, K1D)..., rd) - (; wq, Dr, Ds, Dt, Vq, Vf, wf ) = rd - Nfaces = length(rd.fv) - (; x, y, z, xq, yq, zq, wJq, xf, yf, zf, K ) = md - (; rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ, J ) = md - (; nxJ, nyJ, nzJ, sJ ) = md - (; FToF, mapM, mapP, mapB ) = md - - @test StartUpDG._short_typeof(rd.approximation_type) == "TensorProductWedge{Polynomial, Polynomial}" - @test typeof(md.mesh_type) <: StartUpDG.VertexMappedMesh{<:typeof(rd.element_type)} - @test md.x == md.xyz[1] - @test md.y == md.xyz[2] - @test md.z == md.xyz[3] - - # check positivity of Jacobian - # @show J[1,:] - @test all(J .> 0) - h = estimate_h(rd, md) - @test h <= 2 / K1D + tol - - # check differentiation - u = @. x^2 + 2 * x * y - y^2 + x * y * z - dudx_exact = @. (2*x + 2*y + y*z) - dudy_exact = @. 2*x - 2*y + x*z - dudz_exact = @. x*y - dudr,duds,dudt = (D->D*u).((Dr, Ds, Dt)) - dudx = @. (rxJ * dudr + sxJ * duds + txJ * dudt) / J - dudy = @. (ryJ * dudr + syJ * duds + tyJ * dudt) / J - dudz = @. (rzJ * dudr + szJ * duds + tzJ * dudt) / J - @test dudx ≈ dudx_exact - @test dudy ≈ dudy_exact - @test dudz ≈ dudz_exact - - # check volume integration - @test Vq * x ≈ xq - @test Vq * y ≈ yq - @test Vq * z ≈ zq - @test diagm(wq) * (Vq * J) ≈ wJq - @test abs(sum(xq .* wJq)) < tol - @test abs(sum(yq .* wJq)) < tol - @test abs(sum(zq .* wJq)) < tol - - # check surface integration - @test Vf * x ≈ xf - @test Vf * y ≈ yf - @test Vf * z ≈ zf - @test abs(sum(diagm(wf) * nxJ)) < tol - @test abs(sum(diagm(wf) * nyJ)) < tol - @test abs(sum(diagm(wf) * nzJ)) < tol - @test md.nx .* md.Jf ≈ md.nxJ - @test md.ny .* md.Jf ≈ md.nyJ - @test md.nz .* md.Jf ≈ md.nzJ - - # check connectivity and boundary maps - u = @. (1-x) * (1+x) * (1-y) * (1+y) * (1-z) * (1+z) - uf = Vf * u - @test uf ≈ uf[mapP] - @test norm(uf[mapB]) < tol - - # check periodic node connectivity maps - md_periodic = make_periodic(md, (true, true, true)) - @test md_periodic.mapP != md.mapP # check that the node mapping actually changed - - u = @. sin(pi * (.5 + x)) * sin(pi * (.5 + y)) * sin(pi * (.5 + z)) - (; mapP ) = md_periodic - uf = Vf * u - @test uf ≈ uf[mapP] - - md = MeshData(uniform_mesh(rd.element_type, K1D)..., rd; is_periodic=true) - @test isempty(md.mapB) - - end - end - end end diff --git a/test/MeshData_wedge_pyr_tests.jl b/test/MeshData_wedge_pyr_tests.jl new file mode 100644 index 00000000..e3a139c6 --- /dev/null +++ b/test/MeshData_wedge_pyr_tests.jl @@ -0,0 +1,164 @@ +approx_elem_types_to_test = [(Polynomial(), Wedge()), + (Polynomial(), Pyr())] + +@testset "3D MeshData tests for wedges and pyramids" begin + @testset "$approximation_type $element_type MeshData initialization" for (approximation_type, element_type) in approx_elem_types_to_test + tol = 5e2*eps() + + N = 3 + K1D = 2 + rd = RefElemData(element_type, approximation_type, N) + md = MeshData(uniform_mesh(element_type, K1D)..., rd) + (; wq, Dr, Ds, Dt, Vq, Vf, wf ) = rd + Nfaces = length(rd.fv) + (; x, y, z, xq, yq, zq, wJq, xf, yf, zf, K ) = md + (; rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ, J ) = md + (; nxJ, nyJ, nzJ, sJ ) = md + (; FToF, mapM, mapP, mapB ) = md + + @test typeof(md.mesh_type) <: StartUpDG.VertexMappedMesh{<:typeof(rd.element_type)} + @test md.x == md.xyz[1] + + # check positivity of Jacobian + @test all(J .> 0) + h = estimate_h(rd, md) + @test h <= 2 / K1D + tol + + # check differentiation + u = @. x^2 + 2 * x * y - y^2 + x * y * z + dudx_exact = @. 2*x + 2*y + y*z + dudy_exact = @. 2*x - 2*y + x*z + dudz_exact = @. x*y + dudr,duds,dudt = (D->D*u).((Dr, Ds, Dt)) + dudx = @. (rxJ * dudr + sxJ * duds + txJ * dudt) / J + dudy = @. (ryJ * dudr + syJ * duds + tyJ * dudt) / J + dudz = @. (rzJ * dudr + szJ * duds + tzJ * dudt) / J + @test dudx ≈ dudx_exact + @test dudy ≈ dudy_exact + @test dudz ≈ dudz_exact + + # check volume integration + @test Vq * x ≈ xq + @test Vq * y ≈ yq + @test Vq * z ≈ zq + @test diagm(wq) * (Vq * J) ≈ wJq + @test abs(sum(xq .* wJq)) < tol + @test abs(sum(yq .* wJq)) < tol + @test abs(sum(zq .* wJq)) < tol + + # check surface integration + @test Vf * x ≈ xf + @test Vf * y ≈ yf + @test Vf * z ≈ zf + @test abs(sum(diagm(wf) * nxJ)) < tol + @test abs(sum(diagm(wf) * nyJ)) < tol + @test abs(sum(diagm(wf) * nzJ)) < tol + @test md.nx .* md.Jf ≈ md.nxJ + @test md.ny .* md.Jf ≈ md.nyJ + @test md.nz .* md.Jf ≈ md.nzJ + + # check connectivity and boundary maps + u = @. (1-x) * (1+x) * (1-y) * (1+y) * (1-z) * (1+z) + uf = Vf * u + @test uf ≈ uf[mapP] + @test norm(uf[mapB]) < tol + + # check periodic node connectivity maps + md_periodic = make_periodic(md, (true, true, true)) + @test md_periodic.mapP != md.mapP # check that the node mapping actually changed + + u = @. sin(pi * (.5 + x)) * sin(pi * (.5 + y)) * sin(pi * (.5 + z)) + (; mapP ) = md_periodic + uf = Vf * u + @test uf ≈ uf[mapP] + + md = MeshData(uniform_mesh(rd.element_type, K1D)..., rd; is_periodic=true) + @test isempty(md.mapB) + + end + + @testset "TensorProductWedge MeshData" begin + element_type = Wedge() + tol = 5e2*eps() + @testset "Degree $tri_grad triangle" for tri_grad = [2, 3] + @testset "Degree $line_grad line" for line_grad = [2, 3] + line = RefElemData(Line(), line_grad) + tri = RefElemData(Tri(), tri_grad) + tensor = TensorProductWedge(tri, line) + + rd = RefElemData(element_type, tensor) + K1D = 2 + md = MeshData(uniform_mesh(element_type, K1D)..., rd) + (; wq, Dr, Ds, Dt, Vq, Vf, wf ) = rd + Nfaces = length(rd.fv) + (; x, y, z, xq, yq, zq, wJq, xf, yf, zf, K ) = md + (; rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ, J ) = md + (; nxJ, nyJ, nzJ, sJ ) = md + (; FToF, mapM, mapP, mapB ) = md + + @test StartUpDG._short_typeof(rd.approximation_type) == "TensorProductWedge{Polynomial, Polynomial}" + @test typeof(md.mesh_type) <: StartUpDG.VertexMappedMesh{<:typeof(rd.element_type)} + @test md.x == md.xyz[1] + @test md.y == md.xyz[2] + @test md.z == md.xyz[3] + + # check positivity of Jacobian + @test all(J .> 0) + h = estimate_h(rd, md) + @test h <= 2 / K1D + tol + + # check differentiation + u = @. x^2 + 2 * x * y - y^2 + x * y * z + dudx_exact = @. (2*x + 2*y + y*z) + dudy_exact = @. 2*x - 2*y + x*z + dudz_exact = @. x*y + dudr,duds,dudt = (D->D*u).((Dr, Ds, Dt)) + dudx = @. (rxJ * dudr + sxJ * duds + txJ * dudt) / J + dudy = @. (ryJ * dudr + syJ * duds + tyJ * dudt) / J + dudz = @. (rzJ * dudr + szJ * duds + tzJ * dudt) / J + @test dudx ≈ dudx_exact + @test dudy ≈ dudy_exact + @test dudz ≈ dudz_exact + + # check volume integration + @test Vq * x ≈ xq + @test Vq * y ≈ yq + @test Vq * z ≈ zq + @test diagm(wq) * (Vq * J) ≈ wJq + @test abs(sum(xq .* wJq)) < tol + @test abs(sum(yq .* wJq)) < tol + @test abs(sum(zq .* wJq)) < tol + + # check surface integration + @test Vf * x ≈ xf + @test Vf * y ≈ yf + @test Vf * z ≈ zf + @test abs(sum(diagm(wf) * nxJ)) < tol + @test abs(sum(diagm(wf) * nyJ)) < tol + @test abs(sum(diagm(wf) * nzJ)) < tol + @test md.nx .* md.Jf ≈ md.nxJ + @test md.ny .* md.Jf ≈ md.nyJ + @test md.nz .* md.Jf ≈ md.nzJ + + # check connectivity and boundary maps + u = @. (1-x) * (1+x) * (1-y) * (1+y) * (1-z) * (1+z) + uf = Vf * u + @test uf ≈ uf[mapP] + @test norm(uf[mapB]) < tol + + # check periodic node connectivity maps + md_periodic = make_periodic(md, (true, true, true)) + @test md_periodic.mapP != md.mapP # check that the node mapping actually changed + + u = @. sin(pi * (.5 + x)) * sin(pi * (.5 + y)) * sin(pi * (.5 + z)) + (; mapP ) = md_periodic + uf = Vf * u + @test uf ≈ uf[mapP] + + md = MeshData(uniform_mesh(rd.element_type, K1D)..., rd; is_periodic=true) + @test isempty(md.mapB) + + end + end + end +end diff --git a/test/cut_mesh_tests.jl b/test/cut_mesh_tests.jl index 92ccb15a..93c8cf03 100644 --- a/test/cut_mesh_tests.jl +++ b/test/cut_mesh_tests.jl @@ -1,98 +1,149 @@ using StartUpDG: PathIntersections -@testset "Cut meshes" begin - - cells_per_dimension = 4 - cells_per_dimension_x, cells_per_dimension_y = cells_per_dimension, cells_per_dimension - circle = PresetGeometries.Circle(R=0.33, x0=0, y0=0) - - rd = RefElemData(Quad(), N=3) - md = MeshData(rd, (circle, ), cells_per_dimension_x, cells_per_dimension_y; precompute_operators=true) - - @test_throws ErrorException("Face index f = 5 > 4; too large.") StartUpDG.neighbor_across_face(5, nothing, nothing) - - @test StartUpDG.num_cartesian_elements(md) + StartUpDG.num_cut_elements(md) == md.num_elements - - @test (@capture_out Base.show(stdout, MIME"text/plain"(), md)) == "Cut-cell MeshData of dimension 2 with 16 elements (12 Cartesian, 4 cut)" - - # test constructor with only one "cells_per_dimension" argument - @test_nowarn MeshData(rd, (circle, ), cells_per_dimension_x) - - # check the volume of the domain - A = 4 - pi * .33^2 - @test sum(md.wJq) ≈ A - - # check the length of the boundary of the domain - face_weights = reshape(rd.wf, :, num_faces(rd.element_type))[:, 1] - wJf = vec(Diagonal(face_weights) * reshape(md.Jf, length(face_weights), :)) - @test sum(wJf[md.mapB]) ≈ (8 + 2 * pi * .33) - - # check continuity of a function that's in the global polynomial space - (; physical_frame_elements) = md.mesh_type - (; x, y) = md - u = @. x^rd.N - x * y^(rd.N-1) - x^(rd.N-1) * y + y^rd.N - uf = similar(md.xf) - uf.cartesian .= rd.Vf * u.cartesian - for e in 1:size(md.x.cut, 2) - ids = md.mesh_type.cut_face_nodes[e] - Vf = md.mesh_type.cut_cell_operators.face_interpolation_matrices[e] - uf.cut[ids] .= Vf * u.cut[:, e] +@testset "Cut meshes ($(typeof(quadrature_type)))" for quadrature_type = [Subtriangulation(), MomentFitting()] + @testset "Correctness" begin + + cells_per_dimension = 4 + cells_per_dimension_x, cells_per_dimension_y = cells_per_dimension, cells_per_dimension + circle = PresetGeometries.Circle(R=0.33, x0=0, y0=0) + + rd = RefElemData(Quad(), N=3) + md = MeshData(rd, (circle, ), + cells_per_dimension_x, cells_per_dimension_y, + quadrature_type; + precompute_operators=true) + + @test_throws ErrorException("Face index f = 5 > 4; too large.") StartUpDG.neighbor_across_face(5, nothing, nothing) + + @test StartUpDG.num_cartesian_elements(md) + StartUpDG.num_cut_elements(md) == md.num_elements + + @test (@capture_out Base.show(stdout, MIME"text/plain"(), md)) == "Cut-cell MeshData of dimension 2 with 16 elements (12 Cartesian, 4 cut)" + + # test constructor with only one "cells_per_dimension" argument + @test_nowarn MeshData(rd, (circle, ), cells_per_dimension_x, quadrature_type) + + # check the volume of the domain + A = 4 - pi * .33^2 + @test abs(sum(md.wJq) - A) < 7e-5 # should be 6.417e-5 + + # check the length of the boundary of the domain + wJf = md.mesh_type.cut_cell_data.wJf + @test abs(sum(wJf[md.mapB]) - (8 + 2 * pi * .33)) < 8e-5 # should be 7.66198e-5 + + # check continuity of a function that's in the global polynomial space + (; physical_frame_elements) = md.mesh_type + (; x, y) = md + u = @. x^rd.N - x * y^(rd.N-1) - x^(rd.N-1) * y + y^rd.N + uf = similar(md.xf) + uf.cartesian .= rd.Vf * u.cartesian + for e in 1:size(md.x.cut, 2) + ids = md.mesh_type.cut_face_nodes[e] + Vf = md.mesh_type.cut_cell_operators.face_interpolation_matrices[e] + uf.cut[ids] .= Vf * u.cut[:, e] + end + @test all(uf .≈ vec(uf[md.mapP])) + + dudx_exact = @. rd.N * x^(rd.N-1) - y^(rd.N-1) - (rd.N-1) * x^(rd.N-2) * y + dudy_exact = @. -(rd.N-1) * x * y^(rd.N-2) - x^(rd.N-1) + rd.N * y^(rd.N-1) + (; physical_frame_elements, cut_face_nodes) = md.mesh_type + dudx, dudy = similar(md.x), similar(md.x) + dudx.cartesian .= (md.rxJ.cartesian .* (rd.Dr * u.cartesian)) ./ md.J.cartesian + dudy.cartesian .= (md.syJ.cartesian .* (rd.Ds * u.cartesian)) ./ md.J.cartesian + for (e, elem) in enumerate(physical_frame_elements) + VDM = vandermonde(elem, rd.N, x.cut[:, e], y.cut[:, e]) + Vq, Vxq, Vyq = map(A -> A / VDM, basis(elem, rd.N, md.xq.cut[:,e], md.yq.cut[:, e])) + + M = Vq' * diagm(md.wJq.cut[:, e]) * Vq + Qx = Vq' * diagm(md.wJq.cut[:, e]) * Vxq + Qy = Vq' * diagm(md.wJq.cut[:, e]) * Vyq + Dx, Dy = M \ Qx, M \ Qy + # LIFT = M \ (Vf' * diagm(wJf)) + + # TODO: add interface flux terms into test + dudx.cut[:, e] .= Dx * u.cut[:,e] # (md.rxJ.cut[:,e] .* (Dr * u.cut[:,e])) + dudy.cut[:, e] .= Dy * u.cut[:,e] # (md.rxJ.cut[:,e] .* (Dr * u.cut[:,e])) + end + + @test dudx ≈ dudx_exact + @test dudy ≈ dudy_exact + + # test normals are unit and non-zero + @test all(@. md.nx^2 + md.ny^2 ≈ 1) + + # test creation of equispaced nodes on cut cells + x, y = equi_nodes(physical_frame_elements[1], circle, 10) + # shouldn't have more points than equispaced points on a quad + @test 0 < length(x) <= length(first(equi_nodes(Quad(), 10))) + # no points should be contained in the circle + @test !all(PathIntersections.is_contained.(circle, zip(x, y))) + end - @test all(uf .≈ vec(uf[md.mapP])) - - dudx_exact = @. rd.N * x^(rd.N-1) - y^(rd.N-1) - (rd.N-1) * x^(rd.N-2) * y - dudy_exact = @. -(rd.N-1) * x * y^(rd.N-2) - x^(rd.N-1) + rd.N * y^(rd.N-1) - (; physical_frame_elements, cut_face_nodes) = md.mesh_type - dudx, dudy = similar(md.x), similar(md.x) - dudx.cartesian .= (md.rxJ.cartesian .* (rd.Dr * u.cartesian)) ./ md.J - dudy.cartesian .= (md.syJ.cartesian .* (rd.Ds * u.cartesian)) ./ md.J - for (e, elem) in enumerate(physical_frame_elements) - VDM = vandermonde(elem, rd.N, x.cut[:, e], y.cut[:, e]) - Vq, Vxq, Vyq = map(A -> A / VDM, basis(elem, rd.N, md.xq.cut[:,e], md.yq.cut[:, e])) - - M = Vq' * diagm(md.wJq.cut[:, e]) * Vq - Qx = Vq' * diagm(md.wJq.cut[:, e]) * Vxq - Qy = Vq' * diagm(md.wJq.cut[:, e]) * Vyq - Dx, Dy = M \ Qx, M \ Qy - # LIFT = M \ (Vf' * diagm(wJf)) - - # TODO: add interface flux terms into test - dudx.cut[:, e] .= Dx * u.cut[:,e] # (md.rxJ.cut[:,e] .* (Dr * u.cut[:,e])) - dudy.cut[:, e] .= Dy * u.cut[:,e] # (md.rxJ.cut[:,e] .* (Dr * u.cut[:,e])) + + @testset "State redistribution" begin + # test state redistribution + cells_per_dimension = 4 + circle = PresetGeometries.Circle(R=0.6, x0=0, y0=0) + rd = RefElemData(Quad(), N=4) + md = MeshData(rd, (circle, ), cells_per_dimension, quadrature_type) + + srd = StateRedistribution(rd, md) + e = @. 0 * md.x + 1 # constant + u = @. md.x + md.x^3 .* md.y # degree 4 polynomial + ecopy, ucopy = copy.((e, u)) + + # two ways of applying SRD + apply!(u, srd) + srd(e) # in-place application of SRD functor + + # test exactness + @test norm(e .- ecopy) < 1e2 * length(e) * eps() + @test norm(u .- ucopy) < 1e2 * length(u) * eps() end +end - @test dudx ≈ dudx_exact - @test dudy ≈ dudy_exact +####################################### +# test weak SBP property # +####################################### +@testset "Cut mesh weak SBP property" begin - # test normals are unit and non-zero - @test all(@. md.nx^2 + md.ny^2 ≈ 1) + cells_per_dimension = 2 + circle = PresetGeometries.Circle(R=0.66, x0=.1, y0=0) - # test creation of equispaced nodes on cut cells - x, y = equi_nodes(physical_frame_elements[1], circle, 10) - # shouldn't have more points than equispaced points on a quad - @test 0 < length(x) <= length(first(equi_nodes(Quad(), 10))) - # no points should be contained in the circle - @test !all(PathIntersections.is_contained.(circle, zip(x, y))) + rd = RefElemData(Quad(), N=2) -end - -@testset "State redistribution" begin - # test state redistribution - cells_per_dimension = 4 - circle = PresetGeometries.Circle(R=0.6, x0=0, y0=0) - rd = RefElemData(Quad(), N=4) - md = MeshData(rd, (circle, ), cells_per_dimension, cells_per_dimension) - - srd = StateRedistribution(rd, md) - e = @. 0 * md.x + 1 # constant - u = @. md.x + md.x^3 .* md.y # degree 4 polynomial - ecopy, ucopy = copy.((e, u)) - - # two ways of applying SRD - apply!(u, srd) - srd(e) # in-place application of SRD functor - - # test exactness - @test norm(e .- ecopy) < 1e3 * eps() - @test norm(u .- ucopy) < 1e3 * eps() + md = MeshData(rd, (circle, ), cells_per_dimension; + precompute_operators=true) + + mt = md.mesh_type + + (; x, y, xq, yq) = md + u = @. x^2 + x * y + y^2 + uq = @. xq^2 + xq * yq + yq^2 + + (; volume_interpolation_matrices, differentiation_matrices) = md.mesh_type.cut_cell_operators + Qxyh, hybridized_project_interp_matrices, + hybridized_projection_matrices, hybridized_interp_matrices = + hybridized_SBP_operators(md) + + for e in eachindex(Qxyh) + Qxh, Qyh = Qxyh[e] + Dx, Dy = differentiation_matrices[e] + Vq = volume_interpolation_matrices[e] + M = Vq' * diagm(md.wJq.cut[:,e]) * Vq + Pq = M \ (Vq' * diagm(md.wJq.cut[:,e])) + + # test exactness of interpolation, projection + @test norm(Vq * u.cut[:,e] - uq.cut[:,e]) < 100 * eps() + @test norm(Pq * uq.cut[:,e] - u.cut[:,e]) < 100 * eps() + + # test weak SBP property + @test norm(sum(Qxh, dims=2)) + norm(sum(Qyh, dims=2)) < 100 * eps() + + # test accuracy of hybridized SBP operators when rd::RefElemData uses + # full accuracy quadrature. + Ph = hybridized_projection_matrices[e] + Vh = hybridized_interp_matrices[e] + @test norm(Ph * Qxh * Vh - Dx) + + norm(Ph * Qyh * Vh - Dy) < 100 * length(Dx) * eps() + end end \ No newline at end of file diff --git a/test/named_array_partition_tests.jl b/test/named_array_partition_tests.jl deleted file mode 100644 index 79f3fe85..00000000 --- a/test/named_array_partition_tests.jl +++ /dev/null @@ -1,34 +0,0 @@ -@testset "NamedArrayPartition tests" begin - x = NamedArrayPartition(a = ones(10), b = rand(20)) - @test typeof(@. sin(x * x^2 / x - 1)) <: NamedArrayPartition - @test typeof(x.^2) <: NamedArrayPartition - @test x.a ≈ ones(10) - @test typeof(x .+ x[1:end]) <: Vector # test broadcast precedence - @test all(x .== x[1:end]) - y = copy(x) - @test zero(x, (10, 20)) == zero(x) # test that ignoring dims works - @test typeof(zero(x)) <: NamedArrayPartition - @test (y .*= 2).a[1] ≈ 2 # test in-place bcast - - @test length(Array(x))==30 - @test typeof(Array(x)) <: Array - @test propertynames(x) == (:a, :b) - - x = NamedArrayPartition(a = ones(1), b = 2*ones(1)) - @test Base.summary(x) == string(typeof(x), " with arrays:") - @test (@capture_out Base.show(stdout, MIME"text/plain"(), x)) == "(a = [1.0], b = [2.0])" - - using StructArrays - using StaticArrays: SVector - x = NamedArrayPartition(a = StructArray{SVector{2, Float64}}((ones(5), 2*ones(5))), - b = StructArray{SVector{2, Float64}}((3 * ones(2,2), 4*ones(2,2)))) - @test typeof(x.a) <: StructVector{<:SVector{2}} - @test typeof(x.b) <: StructArray{<:SVector{2}, 2} - @test typeof((x->x[1]).(x)) <: NamedArrayPartition - @test typeof(map(x->x[1], x)) <: NamedArrayPartition -end - -# x = NamedArrayPartition(a = ones(10), b = rand(20)) -# x_ap = ArrayPartition(x) -# @btime @. x_ap * x_ap; # 498.836 ns (5 allocations: 2.77 KiB) -# @btime @. x * x; # 2.032 μs (5 allocations: 2.84 KiB) - 5x slower than ArrayPartition diff --git a/test/reference_elem_tests.jl b/test/reference_elem_tests.jl index 7d487087..46f7a419 100644 --- a/test/reference_elem_tests.jl +++ b/test/reference_elem_tests.jl @@ -24,8 +24,8 @@ @test abs(sum(rd.wf)) ≈ 6 @test abs(sum(rd.wf .* rd.nrJ)) + abs(sum(rd.wf .* rd.nsJ)) < tol @test rd.Pq * rd.Vq ≈ I - Vfp = vandermonde(Line(), N, quad_nodes(Line(), N)[1]) / vandermonde(Line(), N, nodes(Line(), N)) - rstf = (x->Vfp * x[reshape(rd.Fmask, rd.Nfq ÷ rd.Nfaces, rd.Nfaces)]).(rd.rst) + Vf = vandermonde(Line(), N, quad_nodes(Line(), N)[1]) / vandermonde(Line(), N, nodes(Line(), N)) + rstf = (x->Vf * x[reshape(rd.Fmask, :, rd.Nfaces)]).(rd.rst) @test all(vec.(rstf) .≈ rd.rstf) @test StartUpDG.eigenvalue_inverse_trace_constant(rd) ≈ inverse_trace_constant(rd) @test propertynames(rd)[1] == :element_type @@ -44,11 +44,11 @@ @test abs(sum(rd.wf)) ≈ 8 @test abs(sum(rd.wf .* rd.nrJ)) + abs(sum(rd.wf .* rd.nsJ)) < tol @test rd.Pq * rd.Vq ≈ I - Vfp = vandermonde(Line(), N, quad_nodes(Line(), N)[1]) / vandermonde(Line(), N, nodes(Line(), N)) - rstf = (x->Vfp * x[reshape(rd.Fmask,rd.Nfq÷rd.Nfaces,rd.Nfaces)]).(rd.rst) + Vfp = vandermonde(Line(), N, quad_nodes(Line(), N+1)[1]) / vandermonde(Line(), N, nodes(Line(), N)) + rstf = (x->Vfp * x[reshape(rd.Fmask,:,rd.Nfaces)]).(rd.rst) @test all(vec.(rstf) .≈ rd.rstf) - @test StartUpDG.eigenvalue_inverse_trace_constant(rd) ≈ inverse_trace_constant(rd) + @test StartUpDG.eigenvalue_inverse_trace_constant(rd) ≈ inverse_trace_constant(rd) @test StartUpDG.num_vertices(Quad()) == 4 @test StartUpDG.num_faces(Quad()) == 4 end @@ -273,7 +273,7 @@ inverse_trace_constant_compare(rd::RefElemData{3, <:Wedge, <:TensorProductWedge} @test abs(sum(rd.wf .* rd.nsJ)) < tol @test abs(sum(rd.wf .* rd.ntJ)) < tol - @unpack node_ids_by_face = rd.element_type + (; node_ids_by_face) = rd.element_type @test sum(rd.wf[node_ids_by_face[1]]) ≈ 4 # Note: this is not the true area of face 2. Because we map # all faces back to the reference face, there is a factor of @@ -304,25 +304,25 @@ inverse_trace_constant_compare(rd::RefElemData{3, <:Wedge, <:TensorProductWedge} end end -@testset "TensorProductQuadrature on Hex" begin - N = 2 - rd = RefElemData(Hex(), TensorProductQuadrature(quad_nodes(Line(), N+1)), N) - rd_ref = RefElemData(Hex(), N; quad_rule_vol=quad_nodes(Hex(), N+1), quad_rule_face=quad_nodes(Quad(), N+1)) +# @testset "TensorProductQuadrature on Hex" begin +# N = 2 +# rd = RefElemData(Hex(), TensorProductQuadrature(quad_nodes(Line(), N+1)), N) +# rd_ref = RefElemData(Hex(), N; quad_rule_vol=quad_nodes(Hex(), N+1), quad_rule_face=quad_nodes(Quad(), N+1)) - @test typeof(rd) == typeof(RefElemData(Hex(), Polynomial(TensorProductQuadrature(quad_nodes(Line(), N+1))), N)) +# @test typeof(rd) == typeof(RefElemData(Hex(), Polynomial(TensorProductQuadrature(quad_nodes(Line(), N+1))), N)) - for prop in [:N, :element_type] - @test getproperty(rd, prop) == getproperty(rd_ref, prop) - end +# for prop in [:N, :element_type] +# @test getproperty(rd, prop) == getproperty(rd_ref, prop) +# end - for prop in [:fv, :rst, :rstp, :rstq, :rstf, :nrstJ, :Drst] - @test all(getproperty(rd, prop) .≈ getproperty(rd_ref, prop)) - end +# for prop in [:fv, :rst, :rstp, :rstq, :rstf, :nrstJ, :Drst] +# @test all(getproperty(rd, prop) .≈ getproperty(rd_ref, prop)) +# end - for prop in [:Fmask, :VDM, :V1, :wq, :Vq, :wf, :Vf, :Vp, :M, :Pq, :LIFT] - @test norm(getproperty(rd, prop) - getproperty(rd_ref, prop)) < 1e4 * eps() - end -end +# for prop in [:Fmask, :VDM, :V1, :wq, :Vq, :wf, :Vf, :Vp, :M, :Pq, :LIFT] +# @test norm(getproperty(rd, prop) - getproperty(rd_ref, prop)) < 1e4 * eps() +# end +# end @testset "Tensor product Gauss collocation" begin N = 3 @@ -336,55 +336,6 @@ end end end -@testset "Sparse SBP operators" begin - tol = 5e2*eps() - N = 3 - @testset "Line" begin - rd = RefElemData(Line(), Polynomial{Gauss}(), N) - (Q,), E = sparse_low_order_SBP_operators(rd) - @test norm(sum(Q, dims=2)) < tol - @test Q + Q' ≈ E' * Diagonal([-1,1]) * E - end - - @testset "Tri" begin - rd = RefElemData(Tri(), N) - (Qr, Qs), E = sparse_low_order_SBP_operators(rd) - @test norm(sum(Qr, dims=2)) < tol - @test norm(sum(Qs, dims=2)) < tol - @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E - @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E - end - - @testset "Quad (SBP)" begin - rd = RefElemData(Quad(), SBP(), N) - (Qr, Qs), E = sparse_low_order_SBP_operators(rd) - @test norm(sum(Qr, dims=2)) < tol - @test norm(sum(Qs, dims=2)) < tol - @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E - @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E - end - - @testset "Quad (Polynomial{Gauss})" begin - rd = RefElemData(Quad(), Gauss(), N) - (Qr, Qs), E = sparse_low_order_SBP_operators(rd) - @test norm(sum(Qr, dims=2)) < tol - @test norm(sum(Qs, dims=2)) < tol - @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E - @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E - end - - @testset "Hex (Polynomial{Gauss})" begin - rd = RefElemData(Hex(), Gauss(), N) - (Qr, Qs, Qt), E = sparse_low_order_SBP_operators(rd) - @test norm(sum(Qr, dims=2)) < tol - @test norm(sum(Qs, dims=2)) < tol - @test norm(sum(Qt, dims=2)) < tol - @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E - @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E - @test Qt + Qt' ≈ E' * Diagonal(rd.wf .* rd.ntJ) * E - end -end - @testset "RefElemData SBP keyword arguments" begin for elem in [Line(), Tri(), Quad(), Hex()] # make sure Nplot is accepted as a kwarg diff --git a/test/runtests.jl b/test/runtests.jl index af4b0d34..cfdeabbd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,12 +9,13 @@ using SummationByPartsOperators using StartUpDG include("write_vtk_tests.jl") -include("named_array_partition_tests.jl") include("triangulate_tests.jl") include("reference_elem_tests.jl") include("multidim_sbp_tests.jl") +include("sparse_SBP_operator_tests.jl") include("SummationByPartsOperatorsExt_tests.jl") include("MeshData_tests.jl") +include("MeshData_wedge_pyr_tests.jl") include("boundary_util_tests.jl") include("hybrid_mesh_tests.jl") include("noncon_mesh_tests.jl") diff --git a/test/sparse_SBP_operator_tests.jl b/test/sparse_SBP_operator_tests.jl new file mode 100644 index 00000000..758cbb72 --- /dev/null +++ b/test/sparse_SBP_operator_tests.jl @@ -0,0 +1,88 @@ +@testset "Sparse SBP operators for approx_type <: Polynomial" begin + tol = 5e2*eps() + N = 3 + @testset "Line" begin + rd = RefElemData(Line(), Polynomial{Gauss}(), N) + (Q,), E = sparse_low_order_SBP_operators(rd) + @test norm(sum(Q, dims=2)) < tol + @test Q + Q' ≈ E' * Diagonal([-1,1]) * E + end + + @testset "Tri" begin + rd = RefElemData(Tri(), N) + (Qr, Qs), E = sparse_low_order_SBP_operators(rd) + @test norm(sum(Qr, dims=2)) < tol + @test norm(sum(Qs, dims=2)) < tol + @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E + @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E + end + + @testset "Quad (SBP)" begin + rd = RefElemData(Quad(), SBP(), N) + (Qr, Qs), E = sparse_low_order_SBP_operators(rd) + @test norm(sum(Qr, dims=2)) < tol + @test norm(sum(Qs, dims=2)) < tol + @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E + @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E + end + + @testset "Quad (Polynomial{Gauss})" begin + rd = RefElemData(Quad(), Gauss(), N) + (Qr, Qs), E = sparse_low_order_SBP_operators(rd) + @test norm(sum(Qr, dims=2)) < tol + @test norm(sum(Qs, dims=2)) < tol + @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E + @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E + end + + @testset "Hex (Polynomial{Gauss})" begin + rd = RefElemData(Hex(), Gauss(), N) + (Qr, Qs, Qt), E = sparse_low_order_SBP_operators(rd) + @test norm(sum(Qr, dims=2)) < tol + @test norm(sum(Qs, dims=2)) < tol + @test norm(sum(Qt, dims=2)) < tol + @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E + @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E + @test Qt + Qt' ≈ E' * Diagonal(rd.wf .* rd.ntJ) * E + end +end + +@testset "Subcell limiting operators for approx_type <: SBP" begin + @testset "$elem_type" for elem_type in [Tri()] + N = 3 + rd = RefElemData(elem_type, SBP(), N) + Qrst, E = sparse_low_order_SBP_operators(rd) + Δrst, Rrst = subcell_limiting_operators(rd) + for dim in eachindex(Rrst) + @test Δrst[dim] * Rrst[dim] ≈ I + + # check conservation for a random set of limiting factors + r = randn(size(Rrst[dim], 2)) + r = r .- sum(r) / length(r) # so that sum(r) = 0 + θ = rand(size(Rrst[dim], 1)) + @test abs(sum(Δrst[dim] * Diagonal(θ) * Rrst[dim] * r)) < 10 * length(r) * eps() + end + end + + @testset "Tensor product elements: $elem_type" for elem_type in [Quad(), Hex()] + N = 2 + rd = RefElemData(elem_type, SBP(), N) + Qrst, E = sparse_low_order_SBP_operators(rd) + Δrst, Rrst = subcell_limiting_operators(rd) + + for dim in eachindex(Rrst) + conservation_vectors = nullspace(Matrix(Qrst[dim])) + Δ = Δrst[dim] + R = Rrst[dim] + + # make random conservative vector + r = randn(length(rd.r)) + r = r - conservation_vectors * (conservation_vectors' * r) + + # create limiting factors + θ = rand(size(R, 1)) + + @test norm(conservation_vectors' * (Δ * (θ .* (R * r)))) < 100 * eps() + end + end +end \ No newline at end of file diff --git a/test/write_vtk_tests.jl b/test/write_vtk_tests.jl index e2bcdde8..d52eb3ed 100644 --- a/test/write_vtk_tests.jl +++ b/test/write_vtk_tests.jl @@ -11,7 +11,9 @@ deg_one_order(::Hex) = [-1.0 -1.0 1.0 1.0 -1.0 -1.0 1.0 1.0; -1.0 1.0 1.0 -1.0 -1.0 1.0 1.0 -1.0; -1.0 -1.0 -1.0 -1.0 1.0 1.0 1.0 1.0] - deg_one_order(::Wedge) = [-1.0 1.0 -1.0 -1.0 1.0 -1.0; -1.0 -1.0 1.0 -1.0 -1.0 1.0; -1.0 -1.0 -1.0 1.0 1.0 1.0] + deg_one_order(::Wedge) = [-1.0 -1.0 1.0 -1.0 -1.0 1.0 + -1.0 1.0 -1.0 -1.0 1.0 -1.0 + -1.0 -1.0 -1.0 1.0 1.0 1.0] deg_zero_order(elem::Union{Quad, Hex, Wedge}) = deg_one_order(elem)