Refactor RefElemData

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,38 @@ 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. 
+"""
+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. 
 """
 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)
+
+"""
+    TensorProductGaussCollocation
 
-# Polynomial{Gauss} type indicates (N+1)-point Gauss quadrature on tensor product elements
-struct Gauss end 
-Polynomial{Gauss}() = Polynomial(Gauss())
+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 +243,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 +288,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}"

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?

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