add MomentFitting dispatch for old cut cell MeshData

d

src/StartUpDG.jl

@@ -63,6 +63,7 @@ export num_faces, num_vertices, HybridMeshExample
 include("physical_frame_basis.jl")
 include("cut_cell_meshes.jl")
 export PhysicalFrame, equi_nodes
+export MomentFitting
 
 include("state_redistribution.jl")
 export StateRedistribution, apply!

src/cut_cell_meshes.jl

@@ -627,23 +627,24 @@ function subtriangulated_cutcell_quadrature(cutcell, rd_tri,
     return xq, yq, wJq
 end
 
+struct MomentFitting 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:
-```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. 
+    MeshData(quadrature_type::MomentFitting, rd::RefElemData, 
+             objects, vx::AbstractVector, vy::AbstractVector; 
+             quad_rule_face=get_1d_quadrature(rd), 
+             precompute_operators=false)
+
+Constructor for MeshData utilizing moment fitting. Does not guarantee positive
+quadrature weights, and is slower due to the use of adaptive sampling
+to construct 
+
+!!! Warning: this may be deprecated or removed in future versions. 
 """
-MeshData(rd::RefElemData, objects, cells_per_dimension;  kwargs...) = 
-    MeshData(rd::RefElemData, objects, cells_per_dimension, cells_per_dimension;  kwargs...)
+MeshData(quadrature_type::MomentFitting, rd::RefElemData, objects, cells_per_dimension;  kwargs...) = 
+    MeshData(quadrature_type, rd::RefElemData, objects, cells_per_dimension, cells_per_dimension;  kwargs...)
 
-function MeshData(rd::RefElemData, objects, 
+function MeshData(quadrature_type::MomentFitting, 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),
@@ -653,10 +654,10 @@ function MeshData(rd::RefElemData, objects,
     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)
+    return MeshData(quadrature_type, rd, objects, vx, vy; quad_rule_face, precompute_operators)
 end
 
-function MeshData(rd::RefElemData, objects, 
+function MeshData(::MomentFitting, rd::RefElemData, objects, 
                   vx::AbstractVector, vy::AbstractVector; 
                   quad_rule_face=get_1d_quadrature(rd), 
                   precompute_operators=false)

test/cut_mesh_tests.jl

@@ -7,7 +7,9 @@ using StartUpDG: PathIntersections
     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)
+    md = MeshData(MomentFitting(), 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)
 
@@ -16,15 +18,14 @@ using StartUpDG: PathIntersections
     @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)
+    @test_nowarn MeshData(MomentFitting(), 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), :))
+    wJf = md.mesh_type.cut_cell_data.wJf
     @test sum(wJf[md.mapB]) ≈ (8 + 2 * pi * .33)
 
     # check continuity of a function that's in the global polynomial space
@@ -81,7 +82,8 @@ end
     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)
+    md = MeshData(MomentFitting(), rd, (circle, ), 
+                  cells_per_dimension, cells_per_dimension)
 
     srd = StateRedistribution(rd, md)
     e = @. 0 * md.x + 1 # constant