From 576d69312cecde2f07958b5fdd9eb835effa7d13 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Mon, 2 Jan 2023 21:07:40 -0600 Subject: [PATCH] Move `VXYZ` and `EToV` into a `VertexMappedMesh` mesh type (#70) * adding VertexMappedMesh to MeshData * moving VXYZ, EToV into hybrid mesh types * moving VXYZ, EToV into cut cell mesh type * moving VXYZ, EToV into noncon mesh type * remove outdated comment * remove macro for _getproperty, use function with inlining * update docs * update news * add deprecation * cleanup comments * comments * misc formatting --- NEWS.md | 10 +++- docs/src/MeshData.md | 11 ++--- src/MeshData.jl | 91 ++++++++++++++++++++++------------- src/RefElemData_SBP.jl | 18 +++---- src/connectivity_functions.jl | 8 +-- src/cut_cell_meshes.jl | 5 +- src/hybrid_meshes.jl | 10 ++-- src/nonconforming.jl | 10 ++-- test/MeshData_tests.jl | 8 +-- test/hybrid_mesh_tests.jl | 2 +- test/misc_tests.jl | 9 ++-- test/runtests.jl | 2 +- 12 files changed, 110 insertions(+), 74 deletions(-) diff --git a/NEWS.md b/NEWS.md index 94ec101c..f994b266 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,15 +10,21 @@ StartUpDG.jl follows the interpretation of [semantic versioning (semver)](https: #### Changed -* Changes related to element types and `MeshData`: +* Changes related to element types: * upstream change in NodesAndModes.jl: the abstract type `AbstractElemShape` is now parametrized by the dimension, e.g., `Line <: AbstractElemShape{1}`, `Tri <: AbstractElemShape{2}`. * the spatial dimension `Dim` parameter in `MeshData` is now inferred from the element type through `elem <: AbstractElemShape{Dim}` * the `PhysicalFrame` type now has leading type parameter `PhysicalFrame{NDIMS} <: AbstractElemShape{NDIMS}` -* the `Nplot` field has been removed from `RefElemData` * `PhysicalFrame`'s fields are now restricted to be type `SVector` only (instead of `Union{SVector, Tuple}`) +* The `MeshData` fields `VXYZ` and `EToV` have been moved into a `VertexMappedMesh` type. However, they can still be accessed as a property of `MeshData` (e.g., `md.VX` will still work). + +#### Removed + +* the `Nplot` field has been removed from `RefElemData` * all usages of `ComponentArrays` have been replaced by `NamedArrayPartition` +* the deprecated `MeshData(md::MeshData, rd::RefElemData, xyz...)` constructor has been removed #### Deprecated +* the old constructor for `MeshData` with fields `VXYZ` and `EToV` has been deprecated and will be removed in the next breaking release. * the old constructor for `RefElemData` with field `Nplot` has been deprecated and will be removed in the next breaking release. diff --git a/docs/src/MeshData.md b/docs/src/MeshData.md index 4a37fd53..9f392440 100644 --- a/docs/src/MeshData.md +++ b/docs/src/MeshData.md @@ -9,7 +9,8 @@ * `J, Jf`: volume and surface Jacobians evaluated at interpolation points and surface quadrature points, respectively. `J` is a matrix of size ``N_p \times num_elements``, while `Jf` is a matrix of size ``N_f \times num_elements``. * `nxyz::NTuple{Dim, ...}` and `nxyzJ::NTuple{Dim, ...}`: normalized and `Jf` scaled outward normals evaluated at surface quadrature points. Each element of `nxyzJ` is a matrix of size ``N_f\times num_elements``. -These are the main quantities used to construct a DG solver. +These are the main quantities used to construct a DG solver. Information specific to the type of mesh used is +stored in the `md.mesh_type` field. # Setting up `md::MeshData` @@ -51,7 +52,7 @@ julia> md_periodic_x.is_periodic ## Creating curved meshes -It is common to generate curved meshes by first generating a linear mesh, then moving high order nodes on the linear mesh. This can be done by calling [`MeshData`](@ref) again with new `x,y` coordinates: +It is common to generate curved meshes by first generating a linear mesh, then moving high order nodes on the linear mesh. This can be done by calling [`MeshData`](@ref) again with new `x, y` coordinates: ```julia md = MeshData((VX, VY), EToV, rd) @unpack x, y = md @@ -62,11 +63,9 @@ md_curved = MeshData(rd, md, x, y) More generally, one can create a copy of a `MeshData` with certain fields modified by using `@set` or `setproperties` from `Setfield.jl`. -## Unstructured (and pre-defined) triangular meshes using Triangulate +## Unstructured and pre-defined triangular meshes using Triangulate -If `Triangulate` is also loaded, then StartUpDG will include additional utilities for creating and visualizing meshes. - -Several pre-defined geometries are included in StartUpDG.jl. A few examples are `SquareDomain`, `RectangularDomainWithHole`, `Scramjet`, and `CircularDomain`. See `triangulate_example_meshes.jl` for a more complete list and field arguments. These can each be called using `triangulate_domain`, for example the following code will create a mesh of a scramjet: +StartUpDG.jl also includes additional utilities based on Triangulate.jl for creating and visualizing meshes. Several pre-defined geometries are included in StartUpDG.jl. A few examples are `SquareDomain`, `RectangularDomainWithHole`, `Scramjet`, and `CircularDomain`. See `triangulate_example_meshes.jl` for a more complete list and field arguments. These can each be called using `triangulate_domain`, for example the following code will create a mesh of a scramjet: ```julia meshIO = triangulate_domain(Scramjet()) (VX, VY), EToV = triangulateIO_to_VXYEToV(meshIO) diff --git a/src/MeshData.jl b/src/MeshData.jl index baf61910..15e15e7d 100644 --- a/src/MeshData.jl +++ b/src/MeshData.jl @@ -13,19 +13,15 @@ md = MeshElemData(VXY, EToV, rd) @unpack x, y = md ``` """ -Base.@kwdef struct MeshData{Dim, MeshType, VolumeType, FaceType, VolumeQType, - VertexType, EToVType, FToFType, - VolumeWeightType, VolumeGeofacsType, VolumeJType, +Base.@kwdef struct MeshData{Dim, MeshType, + VolumeType, FaceType, VolumeQType, FToFType, + VolumeWeightType, VolumeGeofacsType, VolumeJType, ConnectivityTypeM, ConnectivityTypeP, BoundaryMapType} - # this field defaults to the element shape, but can be - # used to specify cut-cell, hybrid, non-conforming meshes. + # this field defaults to the element shape, but can be used to specify + # different mesh types (e.g., vertex-mapped, cut-cell, hybrid, non-conforming meshes. mesh_type::MeshType - # TODO: move VXYZ, EToV into a `VertexMappedMesh` mesh_type? - VXYZ::NTuple{Dim, VertexType} # vertex coordinates - EToV::EToVType # mesh vertex array - FToF::FToFType # face connectivity xyz::NTuple{Dim, VolumeType} # physical points @@ -33,7 +29,6 @@ Base.@kwdef struct MeshData{Dim, MeshType, VolumeType, FaceType, VolumeQType, xyzq::NTuple{Dim, VolumeQType} # phys quad points, Jacobian-scaled weights wJq::VolumeWeightType - # TODO: move mapP, mapB into a "conforming mesh type"? # arrays of connectivity indices between face nodes mapM::ConnectivityTypeM mapP::ConnectivityTypeP @@ -53,24 +48,24 @@ Base.@kwdef struct MeshData{Dim, MeshType, VolumeType, FaceType, VolumeQType, is_periodic::NTuple{Dim, Bool} end -# MeshData constructor where we do not specify `nxyz` -function MeshData(mesh_type, VXYZ, EToV, FToF, xyz, xyzf, xyzq, wJq, - mapM, mapP, mapB, rstxyzJ, J, nxyzJ, Jf, is_periodic) +# MeshData constructor where we do not specify `nxyz` and instead compute it from `nrstJ` and `Jf` +function MeshData(mesh_type, FToF, xyz, xyzf, xyzq, wJq, mapM, mapP, mapB, rstxyzJ, J, nxyzJ, Jf, is_periodic) nxyz = map(nJ -> nJ ./ Jf, nxyzJ) - return MeshData(mesh_type, VXYZ, EToV, FToF, - xyz, xyzf, xyzq, wJq, mapM, mapP, mapB, - rstxyzJ, J, nxyz, nxyzJ, Jf, is_periodic) + return MeshData(mesh_type, FToF, xyz, xyzf, xyzq, wJq, mapM, mapP, mapB, rstxyzJ, J, nxyz, nxyzJ, Jf, is_periodic) end +# TODO: remove in next breaking release after v0.15. `MeshData` constructor where `VXYZ`, `EToV` are specified +@deprecate MeshData(element_type, VXYZ, EToV, FToF, xyz, xyzf, xyzq, wJq, mapM, mapP, mapB, rstxyzJ, J, nxyzJ, Jf, periodicity) MeshData(VertexMappedMesh(element_type, VXYZ, EToV), FToF, xyz, xyzf, xyzq, wJq, mapM, mapP, mapB, rstxyzJ, J, nxyzJ, Jf, periodicity) + function ConstructionBase.setproperties(md::MeshData, patch::NamedTuple) fields = (haskey(patch, symbol) ? getproperty(patch, symbol) : getproperty(md, symbol) for symbol in fieldnames(typeof(md))) return MeshData(fields...) end ConstructionBase.getproperties(md::MeshData) = - (; mesh_type=md.mesh_type, VXYZ=md.VXYZ, EToV=md.EToV, FToF=md.FToF, xyz=md.xyz, xyzf=md.xyzf, xyzq=md.xyzq, wJq=md.wJq, + (; mesh_type=md.mesh_type, FToF=md.FToF, xyz=md.xyz, xyzf=md.xyzf, xyzq=md.xyzq, wJq=md.wJq, mapM=md.mapM, mapP=md.mapP, mapB=md.mapB, rstxyzJ=md.rstxyzJ, J=md.J, nxyz = md.nxyz, nxyzJ=md.nxyzJ, Jf=md.Jf, is_periodic=md.is_periodic) @@ -81,7 +76,7 @@ end function Base.show(io::IO, ::MIME"text/plain", md::MeshData{DIM, MeshType}) where {DIM, MeshType} @nospecialize md - print(io,"$MeshType MeshData of dimension $DIM with $(md.num_elements) elements") + print(io,"MeshData of dimension $DIM with $(md.num_elements) elements") end function Base.propertynames(x::MeshData{1}, private::Bool = false) @@ -99,17 +94,9 @@ function Base.propertynames(x::MeshData{3}, private::Bool = false) :nx, :ny, :nz, :rxJ, :sxJ, :txJ, :ryJ, :syJ, :tyJ, :rzJ, :szJ, :tzJ) end -# convenience routines for unpacking individual tuple entries -function Base.getproperty(x::MeshData, s::Symbol) - - if s==:VX - return getfield(x, :VXYZ)[1] - elseif s==:VY - return getfield(x, :VXYZ)[2] - elseif s==:VZ - return getfield(x, :VXYZ)[3] - - elseif s==:x +# seems like we need to use @inline to ensure type stability +@inline function meshdata_getproperty(x::MeshData, s::Symbol) + if s==:x return getfield(x, :xyz)[1] elseif s==:y return getfield(x, :xyz)[2] @@ -177,7 +164,44 @@ function Base.getproperty(x::MeshData, s::Symbol) end end -num_elements(md) = size(getfield(md, :EToV), 1) +# generic fallback +Base.getproperty(x::MeshData, s::Symbol) = meshdata_getproperty(x, s) + +""" + struct VertexMappedMesh + +The default `MeshData` mesh type, represents a mesh which is defined purely by +vertex locations and element-to-vertex connectivities. For example, these include +affine triangular meshes or bilinear quadrilateral or trilinear hexahedral meshes. + +# Fields +element_type :: TE <: AbstractElemShape \\ +VXYZ :: TV \\ +EToV :: TEV +""" +struct VertexMappedMesh{TE <: AbstractElemShape, TV, TEV} + element_type::TE + VXYZ::TV + EToV::TEV +end + +# convenience accessor routines for `VertexMappedMesh` types (lets you do `md.VX` instead of `md.mesh_type.VX`) +function Base.getproperty(x::MeshData{Dim, <:VertexMappedMesh}, s::Symbol) where {Dim} + + if s===:VX + return getfield(getfield(x, :mesh_type), :VXYZ)[1] + elseif s===:VY + return getfield(getfield(x, :mesh_type), :VXYZ)[2] + elseif s===:VZ + return getfield(getfield(x, :mesh_type), :VXYZ)[3] + elseif s===:EToV + return getfield(getfield(x, :mesh_type), s) + else + meshdata_getproperty(x, s) + end +end + +num_elements(md::MeshData{Dim, <:VertexMappedMesh}) where {Dim} = size(md.mesh_type.EToV, 1) """ MeshData(VXYZ, EToV, rd::RefElemData) @@ -236,7 +260,7 @@ function MeshData(VX::AbstractVector{Tv}, EToV, rd::RefElemData{1}) where {Tv} is_periodic = (false,) - return MeshData(rd.element_type, tuple(VX), EToV, FToF, + return MeshData(VertexMappedMesh(rd.element_type, tuple(VX), EToV), FToF, tuple(x), tuple(xf), tuple(xq), wJq, collect(mapM), mapP, mapB, SMatrix{1,1}(tuple(rxJ)), J, @@ -276,7 +300,7 @@ function MeshData(VX, VY, EToV, rd::RefElemData{2}) nxJ, nyJ, Jf = compute_normals(rstxyzJ, rd.Vf, rd.nrstJ...) is_periodic = (false, false) - return MeshData(rd.element_type, tuple(VX, VY), EToV, FToF, + return MeshData(VertexMappedMesh(rd.element_type, tuple(VX, VY), EToV), FToF, tuple(x, y), tuple(xf, yf), tuple(xq, yq), wJq, mapM, mapP, mapB, SMatrix{2, 2}(tuple(rxJ, ryJ, sxJ, syJ)), J, @@ -314,13 +338,14 @@ function MeshData(VX, VY, VZ, EToV, rd::RefElemData{3}) nxJ, nyJ, nzJ, Jf = compute_normals(rstxyzJ, rd.Vf, rd.nrstJ...) is_periodic = (false, false, false) - return MeshData(rd.element_type, tuple(VX, VY, VZ), EToV, FToF, + return MeshData(VertexMappedMesh(rd.element_type, tuple(VX, VY, VZ), EToV), FToF, tuple(x, y, z), tuple(xf, yf, zf), tuple(xq, yq, zq), wJq, mapM, mapP, mapB, rstxyzJ, J, tuple(nxJ, nyJ, nzJ), Jf, is_periodic) end +# TODO: remove with breaking change v0.15 @deprecate MeshData(md::MeshData, rd::RefElemData, xyz...) MeshData(rd, md, xyz...) function recompute_geometry(rd::RefElemData{Dim}, xyz) where {Dim} diff --git a/src/RefElemData_SBP.jl b/src/RefElemData_SBP.jl index c20f7469..cba423f0 100644 --- a/src/RefElemData_SBP.jl +++ b/src/RefElemData_SBP.jl @@ -75,24 +75,24 @@ function RefElemData(elementType::Tri, approxType::SBP, N; tol = 100*eps()) (Qrh, Qsh), _ = hybridized_SBP_operators(rd) Nq = length(rd.wq) Vh_sbp = [I(Nq); Ef] - Qr = Vh_sbp'*Qrh*Vh_sbp - Qs = Vh_sbp'*Qsh*Vh_sbp - Dr,Ds = (x->diagm(1 ./ rd.wq) * x).((Qr,Qs)) + Qr = Vh_sbp' * Qrh * Vh_sbp + Qs = Vh_sbp' * Qsh * Vh_sbp + Dr,Ds = (x -> diagm(1 ./ rd.wq) * x).((Qr, Qs)) rd = @set rd.rst = quad_rule_vol[1:2] # set nodes = SBP nodes rd = @set rd.rstq = quad_rule_vol[1:2] # set quad nodes = SBP nodes - rd = @set rd.Drst = (Dr,Ds) + rd = @set rd.Drst = (Dr, Ds) rd = @set rd.Fmask = vec(Fmask) - # TODO: make these more efficient with custom operator? - rd = @set rd.Vf = droptol!(sparse(Ef),tol) - rd = @set rd.LIFT = Diagonal(rd.wq)\(rd.Vf'*Diagonal(rd.wf)) + # TODO: make these more efficient with custom operators? + rd = @set rd.Vf = droptol!(sparse(Ef), tol) + rd = @set rd.LIFT = Diagonal(rd.wq) \ (rd.Vf' * Diagonal(rd.wf)) # make V1 the interpolation matrix from triangle vertices to SBP nodal points - rd = @set rd.V1 = vandermonde(elementType,N,rd.rst...)/rd.VDM * rd.V1 + rd = @set rd.V1 = vandermonde(elementType, N, rd.rst...) / rd.VDM * rd.V1 # Vp operator = projects SBP nodal vector onto degree N polynomial, then interpolates to plotting points - rd = @set rd.Vp = vandermonde(elementType,N,rd.rstp...)/rd.VDM * rd.Pq + rd = @set rd.Vp = vandermonde(elementType, N, rd.rstp...) / rd.VDM * rd.Pq return _convert_RefElemData_fields_to_SBP(rd, approxType) end diff --git a/src/connectivity_functions.jl b/src/connectivity_functions.jl index c9fca77c..012b15a8 100644 --- a/src/connectivity_functions.jl +++ b/src/connectivity_functions.jl @@ -387,14 +387,15 @@ function build_periodic_boundary_maps!(xf, yf, zf, return mapPB[:] end -function compute_boundary_face_centroids(md::MeshData{NDIMS, <:Union{Wedge, Pyr}}) where {NDIMS} - element_type = md.mesh_type +function compute_boundary_face_centroids(md::MeshData{3, <:VertexMappedMesh{<:Union{<:Wedge, <:Pyr}}}) + (; element_type) = md.mesh_type (; node_ids_by_face) = element_type face_nodes_per_element = maximum(maximum.(node_ids_by_face)) boundary_faces = findall(eachindex(md.FToF) .== vec(md.FToF)) # compute face centroids + NDIMS = 3 face_centroids = ntuple(_ -> similar(md.xf, (length(boundary_faces), )), NDIMS) sk = 1 for e in 1:md.num_elements @@ -414,7 +415,8 @@ function compute_boundary_face_centroids(md::MeshData{NDIMS, <:Union{Wedge, Pyr} end # specialize on 3D elements with different types of faces -function make_periodic(md::MeshData{3, <:Union{Wedge, Pyr}}, is_periodic::NTuple{3}; tol=1e-12, kwargs...) +function make_periodic(md::MeshData{3, <:VertexMappedMesh{<:Union{<:Wedge, <:Pyr}}}, + is_periodic::NTuple{3}; tol=1e-12, kwargs...) NDIMS = 3 (; mapM) = md diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index ec0b2af8..00b47958 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -611,9 +611,6 @@ function MeshData(rd::RefElemData, curves, cells_per_dimension_x, cells_per_dime view(wJq.cut, :, e) .= wq end - VXYZ = ntuple(_ -> nothing, 2) - EToV = nothing # dummy field for cut cells - # default to non-periodic is_periodic = (false, false) @@ -642,7 +639,7 @@ function MeshData(rd::RefElemData, curves, cells_per_dimension_x, cells_per_dime cut_face_node_ids = [face_ids(e) for e in 1:num_cut_cells] return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_ids, curves, cut_cell_data), - VXYZ, EToV, FToF, (x, y), (xf, yf), (xq, yq), wJq, + FToF, (x, y), (xf, yf), (xq, yq), wJq, mapM, mapP, mapB, rstxyzJ, J, (nxJ, nyJ), Jf, is_periodic) end diff --git a/src/hybrid_meshes.jl b/src/hybrid_meshes.jl index e4eb7b9a..a69488a8 100644 --- a/src/hybrid_meshes.jl +++ b/src/hybrid_meshes.jl @@ -1,11 +1,13 @@ # mesh_type identifier for hybrid meshes -struct HybridMesh{T} +struct HybridMesh{T, TV, TE} element_types::T + VXYZ::TV + EToV::TE # This is an inner constructor for HybridMesh. It sorts tuples of element types alphabetically # so that `HybridMesh.element_types` always yields the same ordering of element types. - function HybridMesh(element_types::T) where {T <: Tuple} + function HybridMesh(element_types::T, VXYZ, EToV) where {T <: Tuple} p = sortperm(collect(typeof.(element_types)); by=string) - return new{typeof(Tuple(element_types[p]))}(Tuple(element_types[p])) + return new{typeof(Tuple(element_types[p])), typeof(VXYZ), typeof(EToV)}(Tuple(element_types[p]), VXYZ, EToV) end end @@ -170,7 +172,7 @@ function MeshData(VX, VY, EToV_unsorted, rds::LittleDict{<:AbstractElemShape, <: mapM, mapP, mapB = vec.(build_node_maps(FToF, xyzf)) - return MeshData(HybridMesh(Tuple(element_types)), (VX, VY), EToV, FToF, + return MeshData(HybridMesh(Tuple(element_types), (VX, VY), EToV), FToF, xyz, xyzf, xyzq, wJq, mapM, mapP, mapB, rstxyzJ, J, nxyzJ, Jf, diff --git a/src/nonconforming.jl b/src/nonconforming.jl index cb58545d..f96afb28 100644 --- a/src/nonconforming.jl +++ b/src/nonconforming.jl @@ -34,7 +34,9 @@ The `mortar_projection_matrix` similarly maps values from 2 mortar faces back to original non-conforming face. These can be used to create DG solvers on non-conforming meshes. """ -struct NonConformingMesh{CF, NCF, I, P} +struct NonConformingMesh{TV, TE, CF, NCF, I, P} + VXYZ::TV + EToV::TE conforming_faces::CF non_conforming_faces::NCF mortar_interpolation_matrix::I @@ -91,6 +93,8 @@ end # one non-conforming quad face is split into 2 mortar faces num_mortars_per_face(rd::RefElemData{2, Quad}) = 2 +num_elements(md::MeshData{Dim, <:NonConformingMesh}) where {Dim} = size(getproperty(md.mesh_type, :EToV), 1) + function MeshData(mesh::NonConformingQuadMeshExample, rd::RefElemData{2, Quad}) (VX, VY) = mesh.VXY @@ -163,10 +167,10 @@ function MeshData(mesh::NonConformingQuadMeshExample, rd::RefElemData{2, Quad}) is_periodic = (false, false) - mesh_type = NonConformingMesh(conforming_faces, non_conforming_faces, + mesh_type = NonConformingMesh(tuple(VX, VY), EToV, conforming_faces, non_conforming_faces, mortar_interpolation_matrix, mortar_projection_matrix) - return MeshData(mesh_type, tuple(VX, VY), EToV, FToF, + return MeshData(mesh_type, FToF, tuple(x, y), tuple(x_mortar, y_mortar), tuple(xq, yq), wJq, mapM, mapP, mapB, SMatrix{2, 2}(tuple(rxJ, ryJ, sxJ, syJ)), J, diff --git a/test/MeshData_tests.jl b/test/MeshData_tests.jl index 8c9b5929..390bafc1 100644 --- a/test/MeshData_tests.jl +++ b/test/MeshData_tests.jl @@ -13,7 +13,7 @@ @unpack rxJ, J, nxJ, wJq = md @unpack mapM, mapP, mapB = md - @test md.mesh_type==rd.element_type + @test typeof(md.mesh_type) <: StartUpDG.VertexMappedMesh{<:typeof(rd.element_type)} @test md.x == md.xyz[1] @@ -65,7 +65,7 @@ @unpack rxJ, sxJ, ryJ, syJ, J, nxJ, nyJ, sJ, wJq = md @unpack FToF, mapM, mapP, mapB = md - @test md.mesh_type==rd.element_type + @test typeof(md.mesh_type) <: StartUpDG.VertexMappedMesh{<:typeof(rd.element_type)} @test md.x == md.xyz[1] # check positivity of Jacobian @@ -132,7 +132,7 @@ @unpack rxJ, sxJ, ryJ, syJ, J, nxJ, nyJ, sJ, wJq = md @unpack FToF, mapM, mapP, mapB = md - @test md.mesh_type==rd.element_type + @test typeof(md.mesh_type) <: StartUpDG.VertexMappedMesh{<:typeof(rd.element_type)} @test md.x == md.xyz[1] # check positivity of Jacobian @@ -214,7 +214,7 @@ approx_elem_types_to_test = [(Polynomial(), Hex()), @unpack nxJ, nyJ, nzJ, sJ = md @unpack FToF, mapM, mapP, mapB = md - @test md.mesh_type==rd.element_type + @test typeof(md.mesh_type) <: StartUpDG.VertexMappedMesh{<:typeof(rd.element_type)} @test md.x == md.xyz[1] # check positivity of Jacobian diff --git a/test/hybrid_mesh_tests.jl b/test/hybrid_mesh_tests.jl index 6451243b..d63542c9 100644 --- a/test/hybrid_mesh_tests.jl +++ b/test/hybrid_mesh_tests.jl @@ -36,7 +36,7 @@ EToV = [[1 2 4 5], [2 3 5 6], [5 8 9], [4 5 7 8], [9 6 5]] md = MeshData(VX, VY, EToV, rds) - @test md.mesh_type == StartUpDG.HybridMesh((Quad(), Tri())) + @test typeof(md.mesh_type) <: typeof(StartUpDG.HybridMesh((Quad(), Tri()), (VX, VY), EToV)) # test if all nodes on boundary are ±1 @test all(@. abs(max(abs(md.xf[md.mapB]), abs(md.yf[md.mapB])) - 1) < 100 * eps() ) diff --git a/test/misc_tests.jl b/test/misc_tests.jl index 844ff6b5..5ef21004 100644 --- a/test/misc_tests.jl +++ b/test/misc_tests.jl @@ -5,8 +5,8 @@ md = MeshData(uniform_mesh(Tri(),1)...,rd) @test (@capture_out Base.show(stdout,MIME"text/plain"(),rd)) == "RefElemData for a degree 3 Polynomial() approximation on Tri() element." @test (@capture_out Base.show(stdout,rd)) == "RefElemData{N=3,Polynomial(),Tri()}." - @test (@capture_out Base.show(stdout,MIME"text/plain"(),md)) == "Tri MeshData of dimension 2 with 2 elements" - @test (@capture_out Base.show(stdout,md)) == "MeshData{2}" + @test (@capture_out Base.show(stdout, MIME"text/plain"(), md)) == "MeshData of dimension 2 with 2 elements" + @test (@capture_out Base.show(stdout, md)) == "MeshData{2}" # test recipes # see https://discourse.julialang.org/t/how-to-test-plot-recipes/2648/6?u=jlchan @@ -37,10 +37,11 @@ function foo(rd::RefElemData) rd.Nfq, rd.Np, rd.Nq end + @inferred foo(rd) + function foo(md::MeshData) - md.VX, md.VY, md.VZ, md.num_elements, md.K + md.VX, md.VY, md.VZ, md.num_elements end - @inferred foo(rd) @inferred foo(md) # test setproperties diff --git a/test/runtests.jl b/test/runtests.jl index cd771e0f..15ecf96e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,7 +2,7 @@ using Test using Suppressor using LinearAlgebra using RecipesBase -using Triangulate # load before StartUpDG b/c of @require +using Triangulate using StartUpDG include("named_array_partition_tests.jl")