From a8d2becd5c742ffbcbcc7f74daf6b745c3230e8f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 29 Apr 2024 17:20:54 -0500 Subject: [PATCH] add MomentFitting dispatch for old cut cell MeshData d --- src/StartUpDG.jl | 1 + src/cut_cell_meshes.jl | 31 ++++++++++++++++--------------- test/cut_mesh_tests.jl | 12 +++++++----- 3 files changed, 24 insertions(+), 20 deletions(-) diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index 11acfb57..f85c028a 100644 --- a/src/StartUpDG.jl +++ b/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! diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index fc975a4c..b4c900fe 100644 --- a/src/cut_cell_meshes.jl +++ b/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) diff --git a/test/cut_mesh_tests.jl b/test/cut_mesh_tests.jl index 92ccb15a..07922814 100644 --- a/test/cut_mesh_tests.jl +++ b/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