BoundaryTagPlotter(triout::TriangulateIO)

Plot recipe to visualize boundary tags by color. Usage: plot(BoundaryTagPlotter(triout))

struct CurvedMesh{T}

Mesh type indicating that the mesh has been curved. Stores the original mesh type as a field.

Fields

originalmeshtype :: T

CutCellMesh is used in the MeshData field mesh_type for cut cell meshes.

The field physical_frame_elements is a container with shifting/scaling information for each element. We evaluate the physical basis over each element by applying a shifting and scaling of the physical coordinates. The resulting shifted/scaled coordinates then fall into the reference element and can be used to evaluate a reference element basis.

The field cut_face_nodes is a container whose elements are indices of face nodes for a cut element. In other words, md.xf.cut[cut_face_nodes[1]] returns the face nodes of the first element.

We assume all cut elements have the same number of volume quadrature points (which is at least the dimension of a degree 2N polynomial space).

The field objects contains a tuple of the objects used to define the cut region.

The field cut_cell_operators contains optionally precomputed operators (mass, differntiation, face interpolation, and lifting operators).

The field cut_cell_data contains additional data from PathIntersections.

struct MeshData{Dim, Tv, Ti}

MeshData: contains info for a high order piecewise polynomial discretization on an unstructured mesh.

Example:

N, K1D = 3, 2
rd = RefElemData(Tri(), N)
VXY, EToV = uniform_mesh(Tri(), K1D)
md = MeshData(VXY, EToV, rd)
(; x, y ) = md

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> 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(VXYZ, EToV, rd::RefElemData)
MeshData((VXYZ, EToV), rd::RefElemData)

Returns a MeshData struct with high order DG mesh information from the unstructured mesh information (VXYZ..., EToV).

MeshData(rd::RefElemData, md::MeshData, xyz...)

Given new nodal positions xyz... (e.g., from mesh curving), recomputes geometric terms and outputs a new MeshData struct. Only fields modified are the coordinate-dependent terms xyz, xyzf, xyzq, rstxyzJ, J, nxyzJ, Jf.

MeshImportOptions

This struct allows the user to opt for supported features when importing a Gmsh 4.1 .msh file.

Support

• grouping::Bool | On import would you like to include physical group assignements of 2D elements?
• remap_group_name::Bool | On import would you like to maintain or remap physical group ID? Remap results in groupIds in the range 1:number_group_ids.

MeshPlotter(rd::RefElemData, md::RefElemData)

Plot recipe to plot a (possibly curved) quadrilateral or triangular mesh. Usage: plot(MeshPlotter(...))

struct MultipleRefElemData{T <: NamedTuple}
    data::T
end

Holds multiple RefElemData objects in data where typeof(data) <: NamedTuple.

Individual RefElemData can be accessed via getproperty, for example rds.Tri.

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.

This is a proof of concept implementation of a non-conforming mesh in StartUpDG.jl. The intended usage is as follows:

rd = RefElemData(Quad(), N=7)
md = MeshData(NonConformingQuadMeshExample(), rd)

(; x, y) = md
u = @. sin(pi * x) * sin(pi * y)

# interpolate to faces
num_total_faces = num_faces(rd.element_type) * md.num_elements
uf = reshape(rd.Vf * u, :, num_total_faces)

# interpolate faces to mortars (`uf` denotes mortar faces for `NonConformingMesh` types)
(; nonconforming_faces, mortar_interpolation_matrix) = md.mesh_type

u_mortar = reshape(mortar_interpolation_matrix * uf[:, nonconforming_faces], :, 2 * length(nonconforming_faces))

# construct interior (uM = u⁻ "minus") values and exterior (uP = u⁺ "plus") values
uM = hcat(uf, u_mortar) # uM = both element faces and mortar faces
uP = uM[md.mapP]

The mortar_projection_matrix similarly maps values from 2 mortar faces back to values on the original non-conforming face. These can be used to create DG solvers on non-conforming meshes.

`PhysicalFrame{NDIMS} <: AbstractElemShape{NDIMS}`

PhysicalFrame element type. Uses a total degree N approximation space, but is computed with a tensor product Legendre basis as opposed to a triangular PKDO basis. Stores fields shifting and scaling to shift/scale physical coordinates so that they are on the reference element.

PhysicalFrame()
PhysicalFrame(x, y)
PhysicalFrame(x, y, vx, vy): stores coordinates `vx, vy` of background Cartesian cell

Constructors for a PhysicalFrame object (optionally uses arrays of points x, y on a cut element).

Polynomial{T}

Represents polynomial approximation types (as opposed to finite differences). By default, Polynomial() constructs a Polynomial{StartUpDG.DefaultPolynomialType}. Specifying a type parameters allows for dispatch on additional structure within a polynomial approximation (e.g., collocation, tensor product quadrature, etc).

struct RefElemData

RefElemData: contains info (interpolation points, volume/face quadrature, operators) for a high order nodal basis on a given reference element.

Example:

N = 3
rd = RefElemData(Tri(), N)
(; r, s ) = rd

function RefElemData(elem; N, kwargs...)
function RefElemData(elem, approx_type; N, kwargs...)

Keyword argument constructor for RefElemData (to "label" N via rd = RefElemData(Line(), N=3))

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::Line, 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

Constructor for RefElemData for different element types.

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)

SBP reference element data for Quad(), Hex(), and Tri() elements.

For Line(), Quad(), and Hex(), approxType is SBP{TensorProductLobatto}.

For Tri(), approxType can be SBP{Kubatko{LobattoFaceNodes}}, SBP{Kubatko{LegendreFaceNodes}}, or SBP{Hicken}.

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), 
            quad_rule_face=(quad_rule_face_quad, quad_rule_face_tri),
            Nplot=10)

Builds operators for pyramids.

RefElemData(elem::Union{Line, Quad, Hex}, approximation_type::Polynomial{Gauss}, N)

Builds rd::RefElemData with (N+1)-point Gauss quadrature in each dimension.

RefElemData(elem::Wedge, 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), 
            quad_rule_face=(quad_rule_face_quad, quad_rule_face_tri),
            Nplot=10)

Builds operators for prisms/wedges

TensorProductQuadrature{T}

A type parameter to Polynomial indicating that

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

VertexMeshPlotter((VX, VY), EToV, fv)
VertexMeshPlotter(triout::TriangulateIO)

Plot recipe to plot a quadrilateral or triangular mesh. Usage: plot(VertexMeshPlotter(...))

function NodesAndModes.equi_nodes(elem::PhysicalFrame, curve, N)

Returns back Np(N) equally spaced nodes on the background quadrilateral corresponding to elem, with points inside of curve removed.

MeshDatatovtk(md, rd, dim, data, dataname, datatype, filename, writedata = false, equidist_nodes = true)

Translate the given mesh into a vtk-file. md holds a MeshData object rd holds a reference element data/RefElemData of a TensorProductWedge data holds an array of arrays (of size num_nodes by num_elements) with plotting data dataname is an array of strings with name of the associated data write_data, flag if data should be written or not (e.g., if data is not written, only the mesh will be saved as output) equi_dist_nodes flag if points should be interpolated to equidstant nodes

MeshData_to_vtk(md, rd, dim, data, dataname, datatype, filename, write_data = false, equi_dist_nodes = true)

Translate the given mesh into a vtk-file. md holds a MeshData object rd holds a reference element data/RefElemData object. data holds an array of arrays (of size num_nodes by num_elements) with plotting data dataname is an array of strings with name of the associated data write_data, flag if data should be written or not (e.g., if data is not written, only the mesh will be saved as output) equi_dist_nodes flag if points should be interpolated to equidstant nodes

SUD_to_vtk_order(rd::RefElemData, dim)

Compute the permutation of the nodes between StartUpDG and VTK

function boundary_face_centroids(md)

Returns face centroids and boundary_face_ids on the boundaries of the domain given by md::MeshData.

build_node_maps(FToF, Xf)

Intialize the connectivity table along all edges and boundary node tables of all elements. mapM - map minus (interior). mapP - map plus (exterior).

Xf = (xf, yf, zf) and FToF is size (Nfaces * K) and FToF[face] = face neighbor

mapM, mapP are size Nfp x (Nfaces*K)

Examples

julia> mapM, mapP, mapB = build_node_maps(FToF, (xf, yf))

ck45()

Returns coefficients rka,rkb,rkc for the 4th order 5-stage low storage Carpenter/Kennedy Runge Kutta method. Coefficients evolve the residual, solution, and local time, e.g.,

Example

res = rk4a[i]*res + dt*rhs # i = RK stage
@. u += rk4b[i]*res

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
• facecentroids = (facecentroidsx, facecentroidsy), where `facecentroids_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.

connect_mesh(EToV,fv)

Inputs:

• EToV is a num_elements by Nv matrix whose rows identify the Nv vertices

which make up one of the num_elements elements.

• fv (an array of arrays containing unordered indices of face vertices).

Output: FToF, an length(fv) by num_elements index array containing face-to-face connectivity.

estimate_h(rd::RefElemData, md::MeshData)
estimate_h(e, rd::RefElemData, md::MeshData) # e = element index

Estimates the mesh size via min sizeofdomain * |J|/|Jf|, since |J| = O(hᵈ) and |Jf| = O(hᵈ⁻¹).

findline(word::String, lines)

Outputs the line number of word in lines.

It is assumed that the word exists at least once in the file.

function get_boundary_face_labels(triout::TriangulateIO, md::MeshData{2})

Find Triangle segment labels of boundary faces. Returns two arguments:

• boundary_face_tags: tags of faces on the boundary
• boundary_faces: list of faces on the boundary of the domain

function get_node_boundary_tags(triout::TriangulateIO,md::MeshData{2},rd::RefElemData{2,Tri})

Computes node_tags = Nfp x Nfaces * num_elements array where each entry is a Triangulate.jl tag number.

returns the number of elements in a .msh file of a specified dimension

Notes: Gmsh includes elements in a .msh file of multiple dimensions. We want a count of how many

2D elements are in our file. This corisponds to the number of elements in our tri mesh.

hex_vtk_order(corner_verts, order, dim, skip = false)

Compute the coordinates of a VTKLAGRANGEHEXAHEDRON of a hex of order order defined by the coordinates of the vertices given in corner_verts. dim is the dimension of the coordinates given. If skip is set to true, the coordinates of the vertex- and edge-points aren't computed, which can be used to compute points of a VTK_LAGRANGE_HEXHEDRON

Inspired by: https://github.com/ju-kreber/paraview-scripts/blob/master/node_ordering.py

VTK node numbering of a hexagon:

            8+------+7
            /|     /|
           / |    / |
         5+------+6 |

z | 4+–-|–+3 | y | / | / |/ |/ |/ 0–> x 1+–––+2

function hybridized_SBP_operators(rd::RefElemData{DIMS})

Constructs hybridized SBP operators given a RefElemData. Returns operators Qrsth..., VhP, Ph.

function inverse_trace_constant(rd::RefElemData)

Returns the degree-dependent constant in the inverse trace equality over the reference element (as reported in "GPU-accelerated dG methods on hybrid meshes" by Chan, Wang, Modave, Remacle, Warburton 2016).

Can be used to estimate dependence of maximum stable timestep on degree of approximation.

make_periodic(md::MeshData{Dim}, is_periodic...) where {Dim}
make_periodic(md::MeshData{Dim}, is_periodic = ntuple(x->true,Dim)) where {Dim}
make_periodic(md::MeshData, is_periodic = true)

Returns new MeshData such that the node maps mapP and face maps FToF are now periodic. Here, is_periodic is a tuple of Bool indicating whether or not to impose periodic BCs in the x,y, or z coordinate.

n_verts_between(n, from, to, dim)

Compute the coordinates of n equally distributed points between the points given by from and to. dim is the dimension of from and to. Inspired by: https://github.com/ju-kreber/paraview-scripts/blob/master/node_ordering.py

quad_vtk_order(corner_verts, order, dim, skip = false)

Compute the coordinates of a VTKLAGRANGEQUAD of a quad of order order defined by the coordinates of the vertices given in corner_verts. dim is the dimension of the coordinates given. If skip is set to true, the coordinates of the vertex- and edge-points aren't computed, which can be used to compute points of a VTK_LAGRANGE_QUAD Inspired by: https://github.com/ju-kreber/paraview-scripts/blob/master/node_ordering.py

read_Gmsh_2D(filename, args...)

Reads a 2D triangular Gmsh file. Mesh formats 2.2 and 4.1 supported. Returns (VX, VY), EToV.

Examples

VXY, EToV = read_Gmsh_2D("eulerSquareCylinder2D.msh") # v2.2 file format
VXY, EToV = read_Gmsh_2D("test/testset_Gmsh_meshes/periodicity_mesh_v4.msh") # v4.1 file format

# if MeshImportOptions.grouping=true, then a third variable `grouping` is returned
VXY, EToV, grouping = read_Gmsh_2D("test/testset_Gmsh_meshes/periodicity_mesh_v4.msh", MeshImportOptions(true, false))
VXY, EToV, grouping = read_Gmsh_2D("test/testset_Gmsh_meshes/periodicity_mesh_v4.msh", true) # same as above

See also

https://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
https://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format-version-2-0028Legacy0029

read_Gmsh_2D_v2(filename)

Reads triangular GMSH 2D file format 2.2 0 8. Returns (VX, VY), EToV.

Examples

VXY, EToV = read_Gmsh_2D_v2("eulerSquareCylinder2D.msh")

https://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format-version-2-0028Legacy0029

For brevity when grouping is the only supported feature.

example: VXY, EToV, grouping = read_Gmsh_2D_v4("file.msh",true)
example: VXY, EToV = read_Gmsh_2D_v4("file.msh",false)

function read_Gmsh_2D_v4(filename, options)

reads triangular GMSH 2D .msh files.

Output

This depends on if grouping is opted for or not

• returns: (VX, VY), EToV
• returns: (VX, VY), EToV, grouping

Supported formats and features:

• version 4.1 'physical group support 'remap group ids

grouping application

When modeling the wave equation you might want wave speeds to vary across your domain. By assigning physical groups in Gmsh we can maintain such groupings upon importing the .msh file. Each imported element will be a member of a phyical group.

VXY, EToV = read_Gmsh_2D_v4("eulerSquareCylinder2D.msh")
VXY, EToV = read_Gmsh_2D_v4("eulerSquareCylinder2D.msh",false)
VXY, EToV, grouping = read_Gmsh_2D_v4("eulerSquareCylinder2D.msh", true)

option = MeshImportOption(true)
VXY, EToV, grouping = read_Gmsh_2D_v4("eulerSquareCylinder2D.msh", option)

https://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format

Notes: the version 4 format has a more detailed block data format this leads to more complicated parser.

function refine(triout, h, href = h/2)

Refinement of a previous mesh given the current mesh size h. Preserves boundary/volume tags.

remapelementgrouping!(eg::Vector{Int}) GMSH uses integers for naming conventions. This function remaps the Gmsh ids to a list of ids 1:numGroups. This just cleans up a little after Gmsh

Example output

remapelementgrouping([16,16,17,17]) -> [1,1,2,2]

sort_by_axis(corner_verts)

Given the points 'corner_verts' sort them in a lexicographical order and return the permutated points.

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 tag_boundary_faces(md, boundary_name::Symbol = :entire_boundary)
function tag_boundary_faces(md, boundary_list::Dict{Symbol, <:Function})

When called without arguments, just returns Dict(:entire_boundary => boundary_faces)`.

Example usage:

julia> rd = RefElemData(Tri(), N=1)
julia> md = MeshData(uniform_mesh(Tri(), 2)..., rd)
julia> on_bottom_boundary(x, y, tol = 1e-13) = abs(y+1) < tol
julia> on_top_boundary(x, y, tol = 1e-13) = abs(y-1) < tol
julia> tag_boundary_faces(Dict(:bottom => on_bottom_boundary,
                               :top    => on_top_boundary), md)

function tag_boundary_faces(triout::TriangulateIO,
                            rd::RefElemData{2,Tri}, md::MeshData{2},
                            boundary_list::Union{NamedTuple,Dict{Symbol,Int}})

Here, boundary_list is a Dict (or NamedTuple) whose values are the boundary tags for a TriangulateIO mesh format. The output is a Dict or NamedTuple with keys given by boundary_list and values equal to vectors of faces on that given boundary.

Example usage:

julia> using Triangulate, StartUpDG
julia> triout = scramjet()
julia> rd = RefElemData(Tri(),N=1)
julia> md = MeshData(triangulateIO_to_VXYEToV(triout)...,rd)
julia> tag_boundary_faces(triout,rd,md, Dict(:wall=>1, :inflow=>2, :outflow=>3))

triangle_vtk_order(corner_verts, order, dim, skip = false)

Compute the coordinates of a VTK_LAGRANGE_TRIANGLE of a triangle or order order defined by the coordinates of the vertices given in corner_verts. dim is the dimension of the coordinates given. If skip is set to true, the coordinates of the vertex- and edge-points aren't computed, which can be used to compute points of a VTK_LAGRANGE_TRIANGLE Inspired by: https://github.com/ju-kreber/paraview-scripts/blob/master/node_ordering.py

function triangulateIO_to_VXYEToV(triout::TriangulateIO)

Computes VX,VY,EToV from a TriangulateIO object.

type_to_vtk(elem::Hex)
return the VTK-type

type_to_vtk(elem::Quad)
return the VTK-type

type_to_vtk(elem::Tri)
return the VTK-type

type_to_vtk(elem::Wedge)
return the VTK-type

uniform_mesh(elem::Line,Kx)
uniform_mesh(elem::Tri,Kx,Ky)
uniform_mesh(elem::Quad,Kx,Ky)
uniform_mesh(elem::Hex,Kx,Ky,Kz)
uniform_mesh(elem, K)

Uniform Kx (by Ky by Kz) mesh on $[-1,1]^d$, where d is the spatial dimension. Returns (VX,VY,VZ), EToV. When only one K is specified, it assumes a uniform mesh with K elements in each coordinate direction.

K can also be specified using a keyword argument K1D, e.g., uniform_mesh(elem; K1D = 16).

vtk_order(elem::Hex, order)

Construct all node-points of a VTKLAGRANGEHEXAHEDRON of order order. The corner-nodes are given by the reference hexahedron used by StartUpDG in the order defined by vtk.

vtk_order(elem::Quad, order)

Construct all node-points of a VTKLAGRANGEQUAD of order order. The corner-nodes are given by the reference quadrilateral used by StartUpDG in the order defined by vtk

vtk_order(elem::Tri, order)

Construct all node-points of a VTK_LAGRANGE_TRIANGLE of order order. The corner-nodes are given by the reference-triangle used by StartUpDG in the order defined by vtk

vtk_order(elem::Wedge, order)

Construct all node-points of a VTKLAGRANGEWEDGE of order order. The corner-nodes are given by the reference-wedge used by StartUpDG

wedge_vtk_order(corner_verts, order, dim)

Compute the coordinates of a VTKLAGRANGEWEDGE of order order defined by the coordinates of the vertices given in corner_verts. dim is the dimension of the coordinates given. Inspired by: https://github.com/ju-kreber/paraview-scripts/blob/master/node_ordering.py

function Triangulate.triangulate(triin::TriangulateIO, maxarea, minangle=20)

Convenience routine to avoid writing out @sprintf each time. Returns a TriangulateIO object.