diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 71381be4..ea8b561c 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.7.3","generation_timestamp":"2023-09-16T16:05:20","documenter_version":"1.0.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.9.3","generation_timestamp":"2023-09-19T00:42:19","documenter_version":"1.0.1"}} \ No newline at end of file diff --git a/dev/MeshData/index.html b/dev/MeshData/index.html index 17aadc84..92ba1e36 100644 --- a/dev/MeshData/index.html +++ b/dev/MeshData/index.html @@ -29,4 +29,4 @@ boundary_face_dict = tag_boundary_faces(md, Dict(:bottom => on_bottom_boundary, :top => on_top_boundary)) boundary_node_dict = tag_boundary_nodes(rd, md, Dict(:bottom => on_bottom_boundary, :top => on_top_boundary))

You can also specify a list of boundaries using NamedTuples 

boundary_face_dict = tag_boundary_faces(md, (; :bottom=>on_bottom_boundary,:top=>on_top_boundary))
-boundary_node_dict = tag_boundary_nodes(rd, md, (; :bottom=>on_bottom_boundary,:top=>on_top_boundary))
+boundary_node_dict = tag_boundary_nodes(rd, md, (; :bottom=>on_bottom_boundary,:top=>on_top_boundary)) diff --git a/dev/RefElemData/index.html b/dev/RefElemData/index.html index cbbd047d..a9d6e659 100644 --- a/dev/RefElemData/index.html +++ b/dev/RefElemData/index.html @@ -34,4 +34,4 @@ rd = RefElemData(Tri(), SBP{Kubatko{LobattoFaceNodes}}(), N) rd = RefElemData(Tri(), SBP{Kubatko{LegendreFaceNodes}}(), N)

Quadrature rules of both degree 2*N-1 (up to N = 6) and 2*N (up to N = 4) are supported on triangles. For Line, Quad, and Hex elements, RefElemData(..., SBP(), N) is the same as the RefElemData for a DG-SEM discretization, though some fields are specialized for the SBP type. These SBP-based RefElemData objects can also be used to initialize a mesh (for example, md = MeshData(uniform_mesh(rd.element_type, 4)..., rd)).

On triangles, we have the following SBP types with the following properties:

klobatto4

klegendre4

hicken4

Tensor product RefElemData on wedge elements

Experimental implementation

This is an experimental feature and may change in future releases.

There is experimental support for RefElemDatas created from tensor products of triangular and 1D RefElemData objects. 

line = RefElemData(Line(), N_line)
 tri  = RefElemData(Tri(), N_tri)
-rd   = RefElemData(Wedge(), TensorProductWedge(tri, line))

This new rd::RefElemData can then be used to create a wedge-based MeshData. The individual RefElemData objects can be accessed from rd.approximation_type::TensorProductWedge.

+rd = RefElemData(Wedge(), TensorProductWedge(tri, line))

This new rd::RefElemData can then be used to create a wedge-based MeshData. The individual RefElemData objects can be accessed from rd.approximation_type::TensorProductWedge.

diff --git a/dev/conventions/index.html b/dev/conventions/index.html index aef13d71..f8cfdca2 100644 --- a/dev/conventions/index.html +++ b/dev/conventions/index.html @@ -1,2 +1,2 @@ -Background and conventions · StartUpDG.jl

Background

Most high order finite element methods rely on a decomposition of a domain into a mesh of "elements" (e.g., triangles or quadrilaterals in 2D, hexahedra or tetrahedra in 3D). Each "physical" element in a mesh is assumed to be the image of single "reference" element under some geometric mapping. Using the chain rule and changes of variables, one can evaluate integrals and derivatives using only operations on the reference element and some geometric mapping data. This transformation of operations on all elements to a single reference element make finite element methods efficient.

We use the convention that coordinates on the reference element are $r$ in 1D, $r, s$ in 2D, or $r, s, t$ in 3D. Physical coordinates use the standard conventions $x$, $x, y$, and $x, y, z$ in 1D, 2D, and 3D.

Mapping

Derivatives of reference coordinates with respect to physical coordinates are abbreviated, e.g., $\frac{\partial r}{\partial x} = r_x$. Additionally, $J$ is used to denote the determinant of the Jacobian of the reference-to-physical mapping.

Assumptions

We make a few simplifying assumptions about the mesh:

  • meshes are conforming (e.g., each face of an element is shared with at most one other element).
  • the geometric mapping from reference to physical elements is the same degree polynomial as the approximation space on the reference element (e.g., the mapping is isoparametric).

Initial experimental support for hybrid, cut-cell, and non-conforming meshes in two dimensions is also available. Please see the corresponding test sets test/hybrid_mesh_tests.jl, test/cut_mesh_tests.jl, and noncon_mesh_tests.jl for examples.

Code conventions

StartUpDG.jl exports structs RefElemData{Dim, ElemShape, ...} (which contains data associated with the reference element, such as interpolation points, quadrature rules, face nodes, normals, and differentiation/interpolation/projection matrices) and MeshData{Dim} (which contains geometric data associated with a mesh). These are currently used for evaluating DG formulations in a matrix-free fashion. These structs contain fields similar to those in Globals1D, Globals2D, Globals3D in the NDG book codes.

We use the following code conventions:

  • variables r, s, t and x, y, z correspond to values at nodal interpolation points.
  • variables ending in q (e.g., rq, sq,... and xq, yq,...) correspond to values at volume quadrature points.
  • variables ending in f (e.g., rf, sf,... and xf, yf,...) correspond to values at face quadrature points.
  • variables ending in p (e.g., rp, sp,...) correspond to equispaced plotting nodes.
  • Dr, Ds, Dt matrices are nodal differentiation matrices with respect to the $r, s, t$ coordinates. For example, Dr * f.(r, s) approximates the derivative of $f(r, s)$ at nodal points.
  • V matrices correspond to interpolation matrices from nodal interpolation points. For example, Vq interpolates to volume quadrature points, Vf interpolates to face quadrature points, Vp interpolates to plotting nodes.
  • geometric quantities in MeshData are stored as matrices of dimension $\text{number of points per element } \times \text{number of elements}$.

Differences from the codes of "Nodal Discontinuous Galerkin Methods"

The codes in StartUpDG.jl are based closely on the Matlab codes from the book "Nodal Discontinuous Galerkin Methods" by Hesthaven and Warburton (2008) (which we will refer to as "NDG"). However, there are some differences in order to allow for more general DG discretizations and enforce certain mathematical properties:

  • In NDG, Fmask extracts the interpolation nodes which lie on a face. These nodes are then used to compute interface fluxes. However, in StartUpDG.jl, we interpolate nodal values to values at face quadrature points via rd.Vf * u. These operations are equivalent if the interpolation nodes which lie on a face are co-located with quadrature points. Similarly, in NDG, the LIFT matrix maps face nodal points to volume nodal points. In StartUpDG.jl, the rd.LIFT matrix maps from face quadrature points to volume nodal points.
  • in NDG, there are connectivity arrays vmapM and vmapP, which directly retrieve interface values from arrays of nodal values. In StartUpDG.jl, face interpolation nodes are not guaranteed to be co-located with face quadrature nodes, so we do not provide vmapM and vmapP. Instead, we expect the user to compute face values and use the md.mapM, md.mapP arrays to access interface values.
  • in NDG, the mass matrix is computed exactly using the formula M = inv(VDM * VDM'), where VDM is the generalized Vandermonde matrix evaluated at nodal interpolation points. In StartUpDG.jl, the mass matrix is computed using quadrature. These are equivalent if the quadrature is exact for the integrands in the mass matrix (e.g., degree $2N$ polynomials for triangular or tetrahedral elements).
  • in NDG, the geometric terms rx, sx, ry, sy, ... are computed and stored. In StartUpDG.jl, the scaled geometric terms md.rxJ, md.sxJ, md.ryJ, md.syJ, ... are computed, which enable us to enforce the metric identities on curved meshes. Similarly, NDG provides Fscale = sJ ./ J(Fmask, :), while StartUpDG.jl only provides md.Jf, which is equivalent to sJ. Fscale, as well as the NDG geometric terms and can be recovered by dividing by md.J.

Internally, NDG uses arrays EToE and EToF to compute the interface connectivity array mapP. StartUpDG.jl uses instead a face-to-face connectivity array FToF. However, EToE, EToF, and FToF are not typically required for the matrix-free explicit solvers targeted by this package.

+Background and conventions · StartUpDG.jl

Background

Most high order finite element methods rely on a decomposition of a domain into a mesh of "elements" (e.g., triangles or quadrilaterals in 2D, hexahedra or tetrahedra in 3D). Each "physical" element in a mesh is assumed to be the image of single "reference" element under some geometric mapping. Using the chain rule and changes of variables, one can evaluate integrals and derivatives using only operations on the reference element and some geometric mapping data. This transformation of operations on all elements to a single reference element make finite element methods efficient.

We use the convention that coordinates on the reference element are $r$ in 1D, $r, s$ in 2D, or $r, s, t$ in 3D. Physical coordinates use the standard conventions $x$, $x, y$, and $x, y, z$ in 1D, 2D, and 3D.

Mapping

Derivatives of reference coordinates with respect to physical coordinates are abbreviated, e.g., $\frac{\partial r}{\partial x} = r_x$. Additionally, $J$ is used to denote the determinant of the Jacobian of the reference-to-physical mapping.

Assumptions

We make a few simplifying assumptions about the mesh:

  • meshes are conforming (e.g., each face of an element is shared with at most one other element).
  • the geometric mapping from reference to physical elements is the same degree polynomial as the approximation space on the reference element (e.g., the mapping is isoparametric).

Initial experimental support for hybrid, cut-cell, and non-conforming meshes in two dimensions is also available. Please see the corresponding test sets test/hybrid_mesh_tests.jl, test/cut_mesh_tests.jl, and noncon_mesh_tests.jl for examples.

Code conventions

StartUpDG.jl exports structs RefElemData{Dim, ElemShape, ...} (which contains data associated with the reference element, such as interpolation points, quadrature rules, face nodes, normals, and differentiation/interpolation/projection matrices) and MeshData{Dim} (which contains geometric data associated with a mesh). These are currently used for evaluating DG formulations in a matrix-free fashion. These structs contain fields similar to those in Globals1D, Globals2D, Globals3D in the NDG book codes.

We use the following code conventions:

  • variables r, s, t and x, y, z correspond to values at nodal interpolation points.
  • variables ending in q (e.g., rq, sq,... and xq, yq,...) correspond to values at volume quadrature points.
  • variables ending in f (e.g., rf, sf,... and xf, yf,...) correspond to values at face quadrature points.
  • variables ending in p (e.g., rp, sp,...) correspond to equispaced plotting nodes.
  • Dr, Ds, Dt matrices are nodal differentiation matrices with respect to the $r, s, t$ coordinates. For example, Dr * f.(r, s) approximates the derivative of $f(r, s)$ at nodal points.
  • V matrices correspond to interpolation matrices from nodal interpolation points. For example, Vq interpolates to volume quadrature points, Vf interpolates to face quadrature points, Vp interpolates to plotting nodes.
  • geometric quantities in MeshData are stored as matrices of dimension $\text{number of points per element } \times \text{number of elements}$.

Differences from the codes of "Nodal Discontinuous Galerkin Methods"

The codes in StartUpDG.jl are based closely on the Matlab codes from the book "Nodal Discontinuous Galerkin Methods" by Hesthaven and Warburton (2008) (which we will refer to as "NDG"). However, there are some differences in order to allow for more general DG discretizations and enforce certain mathematical properties:

  • In NDG, Fmask extracts the interpolation nodes which lie on a face. These nodes are then used to compute interface fluxes. However, in StartUpDG.jl, we interpolate nodal values to values at face quadrature points via rd.Vf * u. These operations are equivalent if the interpolation nodes which lie on a face are co-located with quadrature points. Similarly, in NDG, the LIFT matrix maps face nodal points to volume nodal points. In StartUpDG.jl, the rd.LIFT matrix maps from face quadrature points to volume nodal points.
  • in NDG, there are connectivity arrays vmapM and vmapP, which directly retrieve interface values from arrays of nodal values. In StartUpDG.jl, face interpolation nodes are not guaranteed to be co-located with face quadrature nodes, so we do not provide vmapM and vmapP. Instead, we expect the user to compute face values and use the md.mapM, md.mapP arrays to access interface values.
  • in NDG, the mass matrix is computed exactly using the formula M = inv(VDM * VDM'), where VDM is the generalized Vandermonde matrix evaluated at nodal interpolation points. In StartUpDG.jl, the mass matrix is computed using quadrature. These are equivalent if the quadrature is exact for the integrands in the mass matrix (e.g., degree $2N$ polynomials for triangular or tetrahedral elements).
  • in NDG, the geometric terms rx, sx, ry, sy, ... are computed and stored. In StartUpDG.jl, the scaled geometric terms md.rxJ, md.sxJ, md.ryJ, md.syJ, ... are computed, which enable us to enforce the metric identities on curved meshes. Similarly, NDG provides Fscale = sJ ./ J(Fmask, :), while StartUpDG.jl only provides md.Jf, which is equivalent to sJ. Fscale, as well as the NDG geometric terms and can be recovered by dividing by md.J.

Internally, NDG uses arrays EToE and EToF to compute the interface connectivity array mapP. StartUpDG.jl uses instead a face-to-face connectivity array FToF. However, EToE, EToF, and FToF are not typically required for the matrix-free explicit solvers targeted by this package.

diff --git a/dev/ex_dg_deriv/index.html b/dev/ex_dg_deriv/index.html index c075e8b2..dde65470 100644 --- a/dev/ex_dg_deriv/index.html +++ b/dev/ex_dg_deriv/index.html @@ -25,4 +25,4 @@ return dudxJ ./ J end

We can visualize the result as follows:

dudx = dg_deriv_x(u, rd, md)
 uxp = Vp * dudx
-scatter(xp, yp, uxp, zcolor=uxp, msw=0, leg=false, ratio=1, cam=(0,90))

Plots of the polynomial approximation $u(x,y)$ and the DG approximation of $\frac{\partial u}{\partial x}$ are given below

u dudx ⠀

+scatter(xp, yp, uxp, zcolor=uxp, msw=0, leg=false, ratio=1, cam=(0,90))

Plots of the polynomial approximation $u(x,y)$ and the DG approximation of $\frac{\partial u}{\partial x}$ are given below

u dudx ⠀

diff --git a/dev/index.html b/dev/index.html index 5031db1c..8a0f5e1a 100644 --- a/dev/index.html +++ b/dev/index.html @@ -17,4 +17,4 @@ # Compute derivatives using geometric mapping + chain rule (; Dr, Ds ) = rd (; rxJ, sxJ, J ) = md -dudx = (rxJ .* (Dr * u) + sxJ .* (Ds * u)) ./ J +dudx = (rxJ .* (Dr * u) + sxJ .* (Ds * u)) ./ J diff --git a/dev/index_refs/index.html b/dev/index_refs/index.html index ca794369..09dc1780 100644 --- a/dev/index_refs/index.html +++ b/dev/index_refs/index.html @@ -1,5 +1,5 @@ -Reference · StartUpDG.jl

Index

Functions

StartUpDG.BoundaryTagPlotter — Type
BoundaryTagPlotter(triout::TriangulateIO)

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

source
StartUpDG.CurvedMesh — Type
struct CurvedMesh{T}

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

Fields

originalmeshtype :: T

source
StartUpDG.CutCellMesh — Type

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.

source
StartUpDG.MeshData — Type
struct MeshData{Dim, Tv, Ti}

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

Example:

N, K1D = 3, 2
+Reference · StartUpDG.jl
Index
Functions
StartUpDG.BoundaryTagPlotter — Type
BoundaryTagPlotter(triout::TriangulateIO)
Plot recipe to visualize boundary tags by color. Usage: plot(BoundaryTagPlotter(triout))
source
StartUpDG.CurvedMesh — Type
struct CurvedMesh{T}
Mesh type indicating that the mesh has been curved. Stores the original mesh type as a field.
Fields
originalmeshtype :: T
source
StartUpDG.CutCellMesh — Type
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.
source
StartUpDG.MeshData — Type
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)
@@ -52,7 +52,7 @@
             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
source
StartUpDG.VertexMappedMesh — Type
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
source
StartUpDG.VertexMeshPlotter — Type
VertexMeshPlotter((VX, VY), EToV, fv)
-VertexMeshPlotter(triout::TriangulateIO)
Plot recipe to plot a quadrilateral or triangular mesh. Usage: plot(VertexMeshPlotter(...))
source
NodesAndModes.equi_nodes — Method
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.
source
StartUpDG.MeshData_to_vtk — Function
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
source
StartUpDG.MeshData_to_vtk — Method
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
source
StartUpDG.build_node_maps — Method
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))
source
StartUpDG.ck45 — Method
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
+VertexMeshPlotter(triout::TriangulateIO)
Plot recipe to plot a quadrilateral or triangular mesh. Usage: plot(VertexMeshPlotter(...))
source
NodesAndModes.equi_nodes — Method
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.
source
StartUpDG.MeshData_to_vtk — Function
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
source
StartUpDG.MeshData_to_vtk — Method
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
source
StartUpDG.build_node_maps — Method
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))
source
StartUpDG.ck45 — Method
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
source
StartUpDG.connect_mesh — Method
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. 
source
StartUpDG.connect_mesh — Method
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.
source
StartUpDG.estimate_h — Method
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ᵈ⁻¹). 
source
StartUpDG.findline — Method
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.
source
StartUpDG.get_boundary_face_labels — Method
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
source
StartUpDG.get_node_boundary_tags — Method
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.
source
StartUpDG.get_num_elements — Function
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.
source
StartUpDG.hex_vtk_order — Function
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
             /|     /|
@@ -90,4 +90,4 @@
 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).
source
StartUpDG.vtk_order — Method
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.
source
StartUpDG.vtk_order — Method
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
source
StartUpDG.vtk_order — Method
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
source
StartUpDG.vtk_order — Method
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
source
StartUpDG.wedge_vtk_order — Method
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
source
Triangulate.triangulate — Function
function Triangulate.triangulate(triin::TriangulateIO, maxarea, minangle=20)
Convenience routine to avoid writing out @sprintf each time. Returns a TriangulateIO object.
source

+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).

source
StartUpDG.vtk_order — Method
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.

source
StartUpDG.vtk_order — Method
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

source
StartUpDG.vtk_order — Method
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

source
StartUpDG.vtk_order — Method
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

source
StartUpDG.wedge_vtk_order — Method
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

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

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

source
diff --git a/dev/more_meshes/index.html b/dev/more_meshes/index.html index d3427bb5..cf3b9059 100644 --- a/dev/more_meshes/index.html +++ b/dev/more_meshes/index.html @@ -49,4 +49,4 @@ uf.cut[ids] .= Vf * u.cut[:, e] end -uf[md.mapP] # these are "exterior" values for each entry of `uf` +uf[md.mapP] # these are "exterior" values for each entry of `uf` diff --git a/dev/search_index.js b/dev/search_index.js index d07e5d22..3e65bc49 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"MeshData/#MeshData-type","page":"MeshData","title":"MeshData type","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"MeshData includes fields such as","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"xyz::NTuple{Dim, ...}: nodal interpolation points mapped to physical elements. All elements of xyz are N_p times N_rm elements matrices, where N_p are the number of nodal points on each element.\nxyzq::NTuple{Dim, ...}, wJq: volume quadrature points/weights mapped to physical elements. All elements these tuples are N_q times N_rm elements matrices, where N_q is the number of quadrature points on each element.\nxyzf::NTuple{Dim, ...}: face quadrature points mapped to physical elements. All elements of xyz are N_f times N_rm elements matrices, where N_f is the number of face points on each element.\nmapP, mapB: indexing arrays for inter-element node connectivity (mapP) and for extracting boundary nodes from the list of face nodes xyzf (mapB). mapP is a matrix of size N_f times N_rm elements, while the length of mapB is the total number of nodes on the boundary.\nrstxyzJ::SMatrix{Dim, Dim}: volume geometric terms G_ij = fracpartial x_ipartial hatx_j. Each element of rstxyzJ is a matrix of size N_p times N_rm elements.\nJ, Jf: volume and surface Jacobians evaluated at interpolation points and surface quadrature points, respectively. J is a matrix of size N_p times N_rm elements, while Jf is a matrix of size N_f times N_rm elements. \nnxyz::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 N_rm elements. ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"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. ","category":"page"},{"location":"MeshData/#Setting-up-md::MeshData","page":"MeshData","title":"Setting up md::MeshData","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"The MeshData struct contains data for high order DG methods useful for evaluating DG formulations in a matrix-free fashion.","category":"page"},{"location":"MeshData/#Generating-unstructured-meshes","page":"MeshData","title":"Generating unstructured meshes","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"For convenience, simple uniform meshes are included in with StartUpDG.jl via uniform_mesh","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"using StartUpDG\nnum_cells_x, num_cells_y, num_cells_z = 4, 2, 8\n(VX,), EToV = uniform_mesh(Line(), num_cells_x)\n(VX, VY), EToV = uniform_mesh(Tri(), num_cells_x, num_cells_y)\n(VX, VY), EToV = uniform_mesh(Quad(), num_cells_x, num_cells_y)\n(VX, VY, VZ), EToV = uniform_mesh(Tet(), num_cells_x, num_cells_y, num_cells_z)\n(VX, VY, VZ), EToV = uniform_mesh(Pyr(), num_cells_x, num_cells_y, num_cells_z)\n(VX, VY, VZ), EToV = uniform_mesh(Wedge(), num_cells_x, num_cells_y, num_cells_z)\n(VX, VY, VZ), EToV = uniform_mesh(Hex(), num_cells_x, num_cells_y, num_cells_z)","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"The uniform triangular mesh is constructed by creating a uniform quadrilateral mesh then bisecting each quad into two triangles. Wedge meshes are constructed similarly. Tet meshes are constructed by dividing each hexahedron into 5 tetrahedral elements. Pyramid meshes are constructed by dividing each hexahedron into 6 pyramids. ","category":"page"},{"location":"MeshData/#Initializing-high-order-DG-mesh-data","page":"MeshData","title":"Initializing high order DG mesh data","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"Given unstructured mesh information (tuple of vertex coordinates VXYZ and index array EToV) high order DG mesh data can be constructed as follows:","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"md = MeshData(VXYZ, EToV, rd)","category":"page"},{"location":"MeshData/#Enforcing-periodic-boundary-conditions","page":"MeshData","title":"Enforcing periodic boundary conditions","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"Periodic boundary conditions can be enforced via the is_periodic keyword argument ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"md_periodic = MeshData((VX, VY), EToV, rd; is_periodic=true) # periodic in both x and y coordinates\nmd_periodic_x = MeshData((VX, VY), EToV, rd; is_periodic=(true, false)) # periodic in x direction, but not y","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"or by calling make_periodic, which returns a new MeshData instance","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"md = MeshData((VX, VY), EToV, rd)\nmd_periodic = make_periodic(md) # periodic in both x and y coordinates\nmd_periodic_x = make_periodic(md, true, false) # periodic in x direction, but not y","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"In either case, the MeshData indexing arrays fields mapP,mapB, and FToF are modified to account for periodicity.","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"One can check which dimensions are periodic via the is_periodic field of MeshData. For example, the md_periodic_x example above gives","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"julia> md_periodic_x.is_periodic\n(true, false)","category":"page"},{"location":"MeshData/#Creating-curved-meshes","page":"MeshData","title":"Creating curved meshes","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"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 again with new x, y coordinates:","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"md = MeshData((VX, VY), EToV, rd)\n(; x, y ) = md\n# <-- code to modify high order nodes (x,y)\nmd_curved = MeshData(rd, md, x, y)","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"MeshData(rd, md, x, y) and MeshData(rd, md, x, y, z) are implemented for 2D and 3D, though this is not currently implemented in 1D.","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"More generally, one can create a copy of a MeshData with certain fields modified by using @set or setproperties from Setfield.jl.","category":"page"},{"location":"MeshData/#Unstructured-and-pre-defined-triangular-meshes-using-Triangulate.jl","page":"MeshData","title":"Unstructured and pre-defined triangular meshes using Triangulate.jl","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"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:","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"meshIO = triangulate_domain(Scramjet())\n(VX, VY), EToV = triangulateIO_to_VXYEToV(meshIO)\nrd = RefElemData(Tri(), 7)\nmd = MeshData((VX, VY), EToV, rd)","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"A quick plot of the face nodes via ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"using Plots\nscatter(vec.(md.xyzf)..., msw=0, ms=1, aspect_ratio=:equal, ylims=(0,2), leg=false)","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"shows the following figure (Image: u)","category":"page"},{"location":"MeshData/#Unstructured-curved-quadrilateral-and-hexahedral-meshes-using-HOHQMesh.jl","page":"MeshData","title":"Unstructured curved quadrilateral and hexahedral meshes using HOHQMesh.jl","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"StartUpDG.jl also reads in .mesh files generated by HOHQMesh.jl. The following code constructs a MeshData which represents a curved quadrilateral mesh generated by HOHQMesh.jl. ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"using StartUpDG\nrd = RefElemData(Quad(), 4)\nhmd = read_HOHQMesh(\"test/testset_HOHQMesh_meshes/easy_example.mesh\", Quad())\nmd = MeshData(hmd, rd)","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"We can visualize the mesh using ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"using Plots\nplot(rd, md) # can also use `plot(MeshPlotter(rd, md))`","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"which yields the following figure:","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"(Image: u)","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"The boundary faces are also automatically tagged with the labels provided in the HOHQMesh file. Each boundary tag and the faces that lie on it are stored in md.mesh_type.boundary_faces. ","category":"page"},{"location":"MeshData/#Tagging-boundary-faces-and-boundary-nodes","page":"MeshData","title":"Tagging boundary faces and boundary nodes","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"One can \"tag\" boundary faces (or boundary nodes) by specifying boolean functions which evaluate to true if a point is on a given boundary segment. ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"rd = RefElemData(Tri(), N=3)\nmd = MeshData(uniform_mesh(Tri(), 1)..., rd)\non_bottom_boundary(point, tol=1e-13) = abs(point[2] + 1) < tol # point = (x,y)\non_top_boundary(point, tol=1e-13) = abs(point[2] - 1) < tol \n\nboundary_face_dict = tag_boundary_faces(md, Dict(:bottom => on_bottom_boundary, :top => on_top_boundary))\nboundary_node_dict = tag_boundary_nodes(rd, md, Dict(:bottom => on_bottom_boundary, :top => on_top_boundary))","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"You can also specify a list of boundaries using NamedTuples ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"boundary_face_dict = tag_boundary_faces(md, (; :bottom=>on_bottom_boundary,:top=>on_top_boundary))\nboundary_node_dict = tag_boundary_nodes(rd, md, (; :bottom=>on_bottom_boundary,:top=>on_top_boundary))","category":"page"},{"location":"RefElemData/#RefElemData-type","page":"RefElemData","title":"RefElemData type","text":"","category":"section"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"RefElemData contains the following fields","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"element_type::AbstractElemShape: element shape. Line, Tri, Quad, Hex, Wedge, Pyr, Tet currently supported.\napproximation_type: approximation type. Defaults to Polynomial(), but SBP() is also supported (see RefElemData based on SBP finite differences).\nNfaces: number of faces on a given type of reference element.\nfv: list of vertices defining faces, e.g., [1,2], [2,3], [3,1] for a triangle\nFmask: indices of interpolation nodes which lie on the faces\nrst::NTuple{Dim, ...}: tuple of vectors of length N_p, each of which contains coordinates of degree N optimized polynomial interpolation points.\nrstq::NTuple{Dim, ...},wq, Vq: tuple of volume quadrature points, vector of weights, and quadrature interpolation matrix. Each element of rstq and wq are vectors of length N_q, and Vq is a matrix of size N_q times N_p.\nN_{\\rm plot}: the degree which determines the number of plotting points N_prm plot.\nrstp::NTuple{Dim, ...}, Vp: tuple of plotting points and plotting interpolation matrix. Each element of rstp is a vector of length N_prm plot, and Vp is a matrix of size N_prm plot times N_p.\nrstf::NTuple{Dim, ...},wf, Vf: tuple of face quadrature points, weights, and face interpolation matrix. Each element of rstf and wf are vectors of length N_f, and Vf is a matrix of size N_f times N_p.\nnrstJ::NTuple{Dim, ...}: tuple of outward reference normals, scaled by the face Jacobian. Each element is a vector of length N_f.\nM: mass matrix computed using quadrature. Size N_p times N_p\nPq: quadrature-based L^2 projection matrix. Size N_p times N_q.\nDrst::NTuple{Dim, ...}, LIFT: differentiation and lifting matrices. Differentiation matrices are size N_p times N_p while lift matrices are size N_ptimes N_f.","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"This list is incomplete; other fields are stored or accessible but currently only used for internal computations.","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"Mass, differentiation, lifting, and interpolation matrices can be specialized. For example, these matrices are dense Matrix{T} type for lines and triangles, but could also be stored as sparse matrices for quadrilaterals and hexahedra.","category":"page"},{"location":"RefElemData/#Setting-up-rd::RefElemData","page":"RefElemData","title":"Setting up rd::RefElemData","text":"","category":"section"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"The struct rd::RefElemData contains data for a given element type. All common reference elements are supported: Line, Tri, Quad, Tet, Pyr, Wedge, and Hex.","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"To initalize a RefElemData, just specify the element type and polynomial degree.","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"N = 3\n\n# 1D elements \nrd = RefElemData(Line(), N)\n\n# 2D elements\nrd = RefElemData(Tri(), N)\nrd = RefElemData(Quad(), N)\n\n# 3D elements\nrd = RefElemData(Tet(), N)\nrd = RefElemData(Pyr(), N)\nrd = RefElemData(Wedge(), N)\nrd = RefElemData(Hex(), N)","category":"page"},{"location":"RefElemData/#Specifying-different-quadrature-rules.","page":"RefElemData","title":"Specifying different quadrature rules.","text":"","category":"section"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"By default, RefElemData initializes volume and surface quadrature rules to be the minimum rules which exactly integrate the unweighted volume and surface mass matrices. If different quadrature rules are desired, they can be specified as follows: ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"N = 3\n\n# create degree N tensor product Gauss-Lobatto rule\nr1D, w1D = gauss_lobatto_quad(0, 0, N)\nrq, sq = vec.(StartUpDG.meshgrid(r1D))\nwr, ws = vec.(StartUpDG.meshgrid(w1D))\nwq = @. wr * ws\n\nrd = RefElemData(Quad(), N; quad_rule_vol = (rq, sq, wq), \n quad_rule_face = (r1D, w1D))","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"This results in a DG spectral element method (DG-SEM) discretization, with a diagonal lumped mass matrix and differentiation matrices which satisfy a summation-by-parts property. ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"By default, RefElemData is constructed for a nodal basis (in order to facilitate curved meshes, connectivity, etc). The interpolation nodes are computed using an interpolatory version of the warp-and-blend procedure. ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"note: Note\nWhile specifying the quadrature rule changes the discretization, it is not reflected in the RefElemData type and thus cannot be specialized on. The following constructors produce RefElemData where the quadrature structure is reflected in the type parameters:rd = RefElemData(Hex(), Polynomial(TensorProductQuadrature(quad_nodes(Line(), N+1)), N)) # tensor product quadrature rules\nrd = RefElemData(Quad(), Polynomial{Gauss}(), N) # (N+1)^d point tensor product Gauss quadrature","category":"page"},{"location":"RefElemData/#RefElemData-based-on-SBP-finite-differences","page":"RefElemData","title":"RefElemData based on SBP finite differences","text":"","category":"section"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"It is also possible to construct a RefElemData based on multi-dimensional SBP finite difference operators. These utilize nodes constructed by Tianheng Chen and Chi-Wang Shu, Ethan Kubatko, and Jason Hicken.","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"Some examples:","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"N = 3\nrd = RefElemData(Quad(), SBP(), N) # defaults to SBP{TensorProductLobatto}\nrd = RefElemData(Quad(), SBP{TensorProductLobatto}(), N) \n\nrd = RefElemData(Hex(), SBP(), N) # defaults to SBP{TensorProductLobatto}\nrd = RefElemData(Hex(), SBP{TensorProductLobatto}(), N) \n\nrd = RefElemData(Tri(), SBP(), N) # defaults to SBP{Kubatko{LobattoFaceNodes}}\nrd = RefElemData(Tri(), SBP{Hicken}(), N) \nrd = RefElemData(Tri(), SBP{Kubatko{LobattoFaceNodes}}(), N) \nrd = RefElemData(Tri(), SBP{Kubatko{LegendreFaceNodes}}(), N) ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"Quadrature rules of both degree 2*N-1 (up to N = 6) and 2*N (up to N = 4) are supported on triangles. For Line, Quad, and Hex elements, RefElemData(..., SBP(), N) is the same as the RefElemData for a DG-SEM discretization, though some fields are specialized for the SBP type. These SBP-based RefElemData objects can also be used to initialize a mesh (for example, md = MeshData(uniform_mesh(rd.element_type, 4)..., rd)). ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"On triangles, we have the following SBP types with the following properties:","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"SBP{Kubatko{LobattoFaceNodes}}: degree 2N-1 accurate quadrature rules with N+2 Lobatto nodes on each face. Nodes for N=4: ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"(Image: klobatto4)","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"SBP{Kubatko{LegendreFaceNodes}}: degree 2N-1 accurate quadrature rules with N+1 Legendre nodes on each face. For N = 1,...,4, these are the same as the nodes constructed by Chen and Shu. Nodes for N=4:","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"(Image: klegendre4)","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"SBP{Hicken}: degree 2N accurate quadrature rules with N+2 Lobatto nodes on each face. Nodes for N=4:","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"(Image: hicken4)","category":"page"},{"location":"RefElemData/#Tensor-product-RefElemData-on-wedge-elements","page":"RefElemData","title":"Tensor product RefElemData on wedge elements","text":"","category":"section"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"warning: Experimental implementation\nThis is an experimental feature and may change in future releases.","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"There is experimental support for RefElemDatas created from tensor products of triangular and 1D RefElemData objects. ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"line = RefElemData(Line(), N_line)\ntri = RefElemData(Tri(), N_tri)\nrd = RefElemData(Wedge(), TensorProductWedge(tri, line))","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"This new rd::RefElemData can then be used to create a wedge-based MeshData. The individual RefElemData objects can be accessed from rd.approximation_type::TensorProductWedge. ","category":"page"},{"location":"conventions/#Background","page":"Background and conventions","title":"Background","text":"","category":"section"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"Most high order finite element methods rely on a decomposition of a domain into a mesh of \"elements\" (e.g., triangles or quadrilaterals in 2D, hexahedra or tetrahedra in 3D). Each \"physical\" element in a mesh is assumed to be the image of single \"reference\" element under some geometric mapping. Using the chain rule and changes of variables, one can evaluate integrals and derivatives using only operations on the reference element and some geometric mapping data. This transformation of operations on all elements to a single reference element make finite element methods efficient. ","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"We use the convention that coordinates on the reference element are r in 1D, r s in 2D, or r s t in 3D. Physical coordinates use the standard conventions x, x y, and x y z in 1D, 2D, and 3D. ","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"(Image: Mapping)","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"Derivatives of reference coordinates with respect to physical coordinates are abbreviated, e.g., fracpartial rpartial x = r_x. Additionally, J is used to denote the determinant of the Jacobian of the reference-to-physical mapping. ","category":"page"},{"location":"conventions/#Assumptions","page":"Background and conventions","title":"Assumptions","text":"","category":"section"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"We make a few simplifying assumptions about the mesh:","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"meshes are conforming (e.g., each face of an element is shared with at most one other element). \nthe geometric mapping from reference to physical elements is the same degree polynomial as the approximation space on the reference element (e.g., the mapping is isoparametric). ","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"Initial experimental support for hybrid, cut-cell, and non-conforming meshes in two dimensions is also available. Please see the corresponding test sets test/hybrid_mesh_tests.jl, test/cut_mesh_tests.jl, and noncon_mesh_tests.jl for examples. ","category":"page"},{"location":"conventions/#Code-conventions","page":"Background and conventions","title":"Code conventions","text":"","category":"section"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"StartUpDG.jl exports structs RefElemData{Dim, ElemShape, ...} (which contains data associated with the reference element, such as interpolation points, quadrature rules, face nodes, normals, and differentiation/interpolation/projection matrices) and MeshData{Dim} (which contains geometric data associated with a mesh). These are currently used for evaluating DG formulations in a matrix-free fashion. These structs contain fields similar to those in Globals1D, Globals2D, Globals3D in the NDG book codes. ","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"We use the following code conventions:","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"variables r, s, t and x, y, z correspond to values at nodal interpolation points. \nvariables ending in q (e.g., rq, sq,... and xq, yq,...) correspond to values at volume quadrature points. \nvariables ending in f (e.g., rf, sf,... and xf, yf,...) correspond to values at face quadrature points. \nvariables ending in p (e.g., rp, sp,...) correspond to equispaced plotting nodes.\nDr, Ds, Dt matrices are nodal differentiation matrices with respect to the r s t coordinates. For example, Dr * f.(r, s) approximates the derivative of f(r s) at nodal points. \nV matrices correspond to interpolation matrices from nodal interpolation points. For example, Vq interpolates to volume quadrature points, Vf interpolates to face quadrature points, Vp interpolates to plotting nodes. \ngeometric quantities in MeshData are stored as matrices of dimension textnumber of points per element times textnumber of elements.","category":"page"},{"location":"conventions/#Differences-from-the-codes-of-\"Nodal-Discontinuous-Galerkin-Methods\"","page":"Background and conventions","title":"Differences from the codes of \"Nodal Discontinuous Galerkin Methods\"","text":"","category":"section"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"The codes in StartUpDG.jl are based closely on the Matlab codes from the book \"Nodal Discontinuous Galerkin Methods\" by Hesthaven and Warburton (2008) (which we will refer to as \"NDG\"). However, there are some differences in order to allow for more general DG discretizations and enforce certain mathematical properties:","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"In NDG, Fmask extracts the interpolation nodes which lie on a face. These nodes are then used to compute interface fluxes. However, in StartUpDG.jl, we interpolate nodal values to values at face quadrature points via rd.Vf * u. These operations are equivalent if the interpolation nodes which lie on a face are co-located with quadrature points. Similarly, in NDG, the LIFT matrix maps face nodal points to volume nodal points. In StartUpDG.jl, the rd.LIFT matrix maps from face quadrature points to volume nodal points. \nin NDG, there are connectivity arrays vmapM and vmapP, which directly retrieve interface values from arrays of nodal values. In StartUpDG.jl, face interpolation nodes are not guaranteed to be co-located with face quadrature nodes, so we do not provide vmapM and vmapP. Instead, we expect the user to compute face values and use the md.mapM, md.mapP arrays to access interface values. \nin NDG, the mass matrix is computed exactly using the formula M = inv(VDM * VDM'), where VDM is the generalized Vandermonde matrix evaluated at nodal interpolation points. In StartUpDG.jl, the mass matrix is computed using quadrature. These are equivalent if the quadrature is exact for the integrands in the mass matrix (e.g., degree 2N polynomials for triangular or tetrahedral elements).\nin NDG, the geometric terms rx, sx, ry, sy, ... are computed and stored. In StartUpDG.jl, the scaled geometric terms md.rxJ, md.sxJ, md.ryJ, md.syJ, ... are computed, which enable us to enforce the metric identities on curved meshes. Similarly, NDG provides Fscale = sJ ./ J(Fmask, :), while StartUpDG.jl only provides md.Jf, which is equivalent to sJ. Fscale, as well as the NDG geometric terms and can be recovered by dividing by md.J. ","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"Internally, NDG uses arrays EToE and EToF to compute the interface connectivity array mapP. StartUpDG.jl uses instead a face-to-face connectivity array FToF. However, EToE, EToF, and FToF are not typically required for the matrix-free explicit solvers targeted by this package. ","category":"page"},{"location":"index_refs/#Index","page":"Reference","title":"Index","text":"","category":"section"},{"location":"index_refs/","page":"Reference","title":"Reference","text":"","category":"page"},{"location":"index_refs/#Functions","page":"Reference","title":"Functions","text":"","category":"section"},{"location":"index_refs/","page":"Reference","title":"Reference","text":"Modules = [StartUpDG]","category":"page"},{"location":"index_refs/#StartUpDG.BoundaryTagPlotter","page":"Reference","title":"StartUpDG.BoundaryTagPlotter","text":"BoundaryTagPlotter(triout::TriangulateIO)\n\nPlot recipe to visualize boundary tags by color. Usage: plot(BoundaryTagPlotter(triout))\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.CurvedMesh","page":"Reference","title":"StartUpDG.CurvedMesh","text":"struct CurvedMesh{T}\n\nMesh type indicating that the mesh has been curved. Stores the original mesh type as a field.\n\nFields\n\noriginalmeshtype :: T\n\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.CutCellMesh","page":"Reference","title":"StartUpDG.CutCellMesh","text":"CutCellMesh is used in the MeshData field mesh_type for cut cell meshes.\n\nThe 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. \n\nThe 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. \n\nWe assume all cut elements have the same number of volume quadrature points (which is at least the dimension of a degree 2N polynomial space). \n\nThe field objects contains a tuple of the objects used to define the cut region.\n\nThe field cut_cell_operators contains optionally precomputed operators (mass, differntiation, face interpolation, and lifting operators). \n\nThe field cut_cell_data contains additional data from PathIntersections.\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.MeshData","page":"Reference","title":"StartUpDG.MeshData","text":"struct MeshData{Dim, Tv, Ti}\n\nMeshData: contains info for a high order piecewise polynomial discretization on an unstructured mesh. \n\nExample:\n\nN, K1D = 3, 2\nrd = RefElemData(Tri(), N)\nVXY, EToV = uniform_mesh(Tri(), K1D)\nmd = MeshData(VXY, EToV, rd)\n(; x, y ) = md\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.MeshData-Tuple{RefElemData, Any, Any}","page":"Reference","title":"StartUpDG.MeshData","text":"function MeshData(rd, geometry, vxyz...)\n\nCreates 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:\n\njulia> geometry = (PresetGeometries.Circle(R=0.33, x0=0, y0=0), )\n\nHere, 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. \n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.MeshData-Tuple{Tuple{Tuple, Matrix{Int64}}, RefElemData, Vararg{Any}}","page":"Reference","title":"StartUpDG.MeshData","text":"MeshData(VXYZ, EToV, rd::RefElemData)\nMeshData((VXYZ, EToV), rd::RefElemData)\n\nReturns a MeshData struct with high order DG mesh information from the unstructured mesh information (VXYZ..., EToV).\n\nMeshData(rd::RefElemData, md::MeshData, xyz...)\n\nGiven 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.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.MeshImportOptions","page":"Reference","title":"StartUpDG.MeshImportOptions","text":"MeshImportOptions\n\nThis struct allows the user to opt for supported features when importing a Gmsh 4.1 .msh file.\n\nSupport\n\ngrouping::Bool | On import would you like to include physical group assignements of 2D elements?\nremap_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.\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.MeshPlotter","page":"Reference","title":"StartUpDG.MeshPlotter","text":"MeshPlotter(rd::RefElemData, md::RefElemData)\n\nPlot recipe to plot a (possibly curved) quadrilateral or triangular mesh. Usage: plot(MeshPlotter(...))\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.MultipleRefElemData","page":"Reference","title":"StartUpDG.MultipleRefElemData","text":"struct MultipleRefElemData{T <: NamedTuple}\n data::T\nend\n\nHolds multiple RefElemData objects in data where typeof(data) <: NamedTuple.\n\nIndividual RefElemData can be accessed via getproperty, for example rds.Tri. \n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.NamedArrayPartition","page":"Reference","title":"StartUpDG.NamedArrayPartition","text":"NamedArrayPartition(; kwargs...)\nNamedArrayPartition(x::NamedTuple)\n\nSimilar 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. \n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.NonConformingMesh","page":"Reference","title":"StartUpDG.NonConformingMesh","text":"warning: Experimental implementation\nThis is an experimental feature and may change without warning in future releases.\n\nThis is a proof of concept implementation of a non-conforming mesh in StartUpDG.jl. The intended usage is as follows:\n\nrd = RefElemData(Quad(), N=7)\nmd = MeshData(NonConformingQuadMeshExample(), rd)\n\n(; x, y) = md\nu = @. sin(pi * x) * sin(pi * y)\n\n# interpolate to faces\nnum_total_faces = num_faces(rd.element_type) * md.num_elements\nuf = reshape(rd.Vf * u, :, num_total_faces)\n\n# interpolate faces to mortars (`uf` denotes mortar faces for `NonConformingMesh` types)\n(; nonconforming_faces, mortar_interpolation_matrix) = md.mesh_type\n\nu_mortar = reshape(mortar_interpolation_matrix * uf[:, nonconforming_faces], :, 2 * length(nonconforming_faces))\n\n# construct interior (uM = u⁻ \"minus\") values and exterior (uP = u⁺ \"plus\") values\nuM = hcat(uf, u_mortar) # uM = both element faces and mortar faces\nuP = uM[md.mapP]\n\nThe 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.\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.PhysicalFrame","page":"Reference","title":"StartUpDG.PhysicalFrame","text":"`PhysicalFrame{NDIMS} <: AbstractElemShape{NDIMS}`\n\nPhysicalFrame 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. \n\nPhysicalFrame()\nPhysicalFrame(x, y)\nPhysicalFrame(x, y, vx, vy): stores coordinates `vx, vy` of background Cartesian cell\n\nConstructors for a PhysicalFrame object (optionally uses arrays of points x, y on a cut element).\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.Polynomial","page":"Reference","title":"StartUpDG.Polynomial","text":"Polynomial{T}\n\nRepresents 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). \n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.RefElemData","page":"Reference","title":"StartUpDG.RefElemData","text":"struct RefElemData\n\nRefElemData: contains info (interpolation points, volume/face quadrature, operators) for a high order nodal basis on a given reference element. \n\nExample:\n\nN = 3\nrd = RefElemData(Tri(), N)\n(; r, s ) = rd\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"function RefElemData(elem; N, kwargs...)\nfunction RefElemData(elem, approx_type; N, kwargs...)\n\nKeyword argument constructor for RefElemData (to \"label\" N via rd = RefElemData(Line(), N=3))\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Hex, TensorProductQuadrature, Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"RefElemData(elem::Hex, ::TensorProductQuadrature, N; Nplot = 10)\nRefElemData(elem::Hex, approximation_type::Polynomial{<:TensorProductQuadrature}, N; Nplot = 10)\n\nConstructor for hexahedral RefElemData where the quadrature is assumed to have tensor product structure.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Line, Polynomial{StartUpDG.DefaultPolynomialType}, Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"RefElemData(elem::Line, N;\n quad_rule_vol = quad_nodes(elem, N+1))\nRefElemData(elem::Union{Tri, Quad}, N;\n quad_rule_vol = quad_nodes(elem, N),\n quad_rule_face = gauss_quad(0, 0, N))\nRefElemData(elem::Union{Hex, Tet}, N;\n quad_rule_vol = quad_nodes(elem, N),\n quad_rule_face = quad_nodes(Quad(), N))\nRefElemData(elem; N, kwargs...) # version with keyword args\n\nConstructor for RefElemData for different element types.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Line, SBP{TensorProductLobatto}, Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"function RefElemData(elementType::Line, approxType::SBP, N)\nfunction RefElemData(elementType::Quad, approxType::SBP, N)\nfunction RefElemData(elementType::Hex, approxType::SBP, N)\nfunction RefElemData(elementType::Tri, approxType::SBP, N)\n\nSBP reference element data for Quad(), Hex(), and Tri() elements. \n\nFor Line(), Quad(), and Hex(), approxType is SBP{TensorProductLobatto}.\n\nFor Tri(), approxType can be SBP{Kubatko{LobattoFaceNodes}}, SBP{Kubatko{LegendreFaceNodes}}, or SBP{Hicken}. \n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Pyr, Polynomial, Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"RefElemData(elem::Pyr, approximation_type::Polynomial, N;\n quad_rule_vol=quad_nodes(elem, N),\n quad_rule_face_quad=quad_nodes(Quad(), N), \n quad_rule_face_tri=quad_nodes(Tri(), N), \n quad_rule_face=(quad_rule_face_quad, quad_rule_face_tri),\n Nplot=10)\n\nBuilds operators for pyramids.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Union{Hex, Line, Quad}, Polynomial{<:Gauss}, Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"RefElemData(elem::Union{Line, Quad, Hex}, approximation_type::Polynomial{Gauss}, N)\n\nBuilds rd::RefElemData with (N+1)-point Gauss quadrature in each dimension. \n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Wedge, Polynomial, Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"RefElemData(elem::Wedge, approximation_type::Polynomial, N;\n quad_rule_vol=quad_nodes(elem, N),\n quad_rule_face_quad=quad_nodes(Quad(), N), \n quad_rule_face_tri=quad_nodes(Tri(), N), \n quad_rule_face=(quad_rule_face_quad, quad_rule_face_tri),\n Nplot=10)\n\nBuilds operators for prisms/wedges\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.TensorProductQuadrature","page":"Reference","title":"StartUpDG.TensorProductQuadrature","text":"TensorProductQuadrature{T}\n\nA type parameter to Polynomial indicating that \n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.VertexMappedMesh","page":"Reference","title":"StartUpDG.VertexMappedMesh","text":"struct VertexMappedMesh\n\nThe 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.\n\nFields\n\nelement_type :: TE <: AbstractElemShape \nVXYZ :: TV \nEToV :: TEV\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.VertexMeshPlotter","page":"Reference","title":"StartUpDG.VertexMeshPlotter","text":"VertexMeshPlotter((VX, VY), EToV, fv)\nVertexMeshPlotter(triout::TriangulateIO)\n\nPlot recipe to plot a quadrilateral or triangular mesh. Usage: plot(VertexMeshPlotter(...))\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#NodesAndModes.equi_nodes-Tuple{PhysicalFrame{2, Shifting, Scaling} where {Shifting<:(StaticArraysCore.SVector{2}), Scaling<:(StaticArraysCore.SVector{2})}, Any, Any}","page":"Reference","title":"NodesAndModes.equi_nodes","text":"function NodesAndModes.equi_nodes(elem::PhysicalFrame, curve, N)\n\nReturns back Np(N) equally spaced nodes on the background quadrilateral corresponding to elem, with points inside of curve removed.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.MeshData_to_vtk","page":"Reference","title":"StartUpDG.MeshData_to_vtk","text":"MeshDatatovtk(md, rd, dim, data, dataname, datatype, filename, writedata = false, equidist_nodes = true)\n\nTranslate 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\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.MeshData_to_vtk-Union{Tuple{DIM}, Tuple{MeshData, RefElemData{DIM, ElemShape} where ElemShape<:AbstractElemShape{DIM}, Any, Any, Any}, Tuple{MeshData, RefElemData{DIM, ElemShape} where ElemShape<:AbstractElemShape{DIM}, Any, Any, Any, Any}, Tuple{MeshData, RefElemData{DIM, ElemShape} where ElemShape<:AbstractElemShape{DIM}, Any, Any, Any, Any, Any}} where DIM","page":"Reference","title":"StartUpDG.MeshData_to_vtk","text":"MeshData_to_vtk(md, rd, dim, data, dataname, datatype, filename, write_data = false, equi_dist_nodes = true)\n\nTranslate 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\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.SUD_to_vtk_order-Union{Tuple{RefElemData{DIM, ElemShape} where ElemShape<:AbstractElemShape{DIM}}, Tuple{DIM}} where DIM","page":"Reference","title":"StartUpDG.SUD_to_vtk_order","text":"SUD_to_vtk_order(rd::RefElemData, dim)\n\nCompute the permutation of the nodes between StartUpDG and VTK\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.boundary_face_centroids-Tuple{Any}","page":"Reference","title":"StartUpDG.boundary_face_centroids","text":"function boundary_face_centroids(md)\n\nReturns face centroids and boundary_face_ids on the boundaries of the domain given by md::MeshData.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.build_node_maps-Tuple{Any, Any}","page":"Reference","title":"StartUpDG.build_node_maps","text":"build_node_maps(FToF, Xf)\n\nIntialize the connectivity table along all edges and boundary node tables of all elements. mapM - map minus (interior). mapP - map plus (exterior).\n\nXf = (xf, yf, zf) and FToF is size (Nfaces * K) and FToF[face] = face neighbor\n\nmapM, mapP are size Nfp x (Nfaces*K)\n\nExamples\n\njulia> mapM, mapP, mapB = build_node_maps(FToF, (xf, yf))\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.ck45-Tuple{}","page":"Reference","title":"StartUpDG.ck45","text":"ck45()\n\nReturns 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.,\n\nExample\n\nres = rk4a[i]*res + dt*rhs # i = RK stage\n@. u += rk4b[i]*res\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.connect_mesh-Tuple{Any, Any, Any}","page":"Reference","title":"StartUpDG.connect_mesh","text":"connect_mesh(rd, face_centroids, region_flags, cutcells; tol = 1e2 * eps())\n\nConnects faces of a cut mesh to each other, returns FToF such that face f is connected to FToF[f]. \n\nInputs:\n\nrd::RefElemData\nfacecentroids = (facecentroidsx, facecentroidsy), where `facecentroids_x/y` are vectors of coordinates of face centroids\nregion_flags, cutcells are return arguments from PathIntersections.define_regions\n\nThe keyword argument tol is the tolerance for matches between face centroids. \n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.connect_mesh-Tuple{Any, Any}","page":"Reference","title":"StartUpDG.connect_mesh","text":"connect_mesh(EToV,fv)\n\nInputs:\n\nEToV is a num_elements by Nv matrix whose rows identify the Nv vertices\n\nwhich make up one of the num_elements elements.\n\nfv (an array of arrays containing unordered indices of face vertices).\n\nOutput: FToF, an length(fv) by num_elements index array containing face-to-face connectivity.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.estimate_h-Union{Tuple{DIM}, Tuple{RefElemData{DIM, ElemShape} where ElemShape<:AbstractElemShape{DIM}, MeshData{DIM}}} where DIM","page":"Reference","title":"StartUpDG.estimate_h","text":"estimate_h(rd::RefElemData, md::MeshData)\nestimate_h(e, rd::RefElemData, md::MeshData) # e = element index\n\nEstimates the mesh size via min sizeofdomain * |J|/|Jf|, since |J| = O(hᵈ) and |Jf| = O(hᵈ⁻¹). \n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.findline-Tuple{String, Vector{String}}","page":"Reference","title":"StartUpDG.findline","text":"findline(word::String, lines)\n\nOutputs the line number of word in lines. \n\nIt is assumed that the word exists at least once in the file.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.get_boundary_face_labels-Tuple{Triangulate.TriangulateIO, RefElemData{2, Tri}, MeshData{2}}","page":"Reference","title":"StartUpDG.get_boundary_face_labels","text":"function get_boundary_face_labels(triout::TriangulateIO, md::MeshData{2})\n\nFind Triangle segment labels of boundary faces. Returns two arguments:\n\nboundary_face_tags: tags of faces on the boundary\nboundary_faces: list of faces on the boundary of the domain\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.get_node_boundary_tags-Tuple{Triangulate.TriangulateIO, RefElemData{2, Tri}, MeshData{2}}","page":"Reference","title":"StartUpDG.get_node_boundary_tags","text":"function get_node_boundary_tags(triout::TriangulateIO,md::MeshData{2},rd::RefElemData{2,Tri})\n\nComputes node_tags = Nfp x Nfaces * num_elements array where each entry is a Triangulate.jl tag number.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.get_num_elements","page":"Reference","title":"StartUpDG.get_num_elements","text":"returns the number of elements in a .msh file of a specified dimension\n\nNotes: Gmsh includes elements in a .msh file of multiple dimensions. We want a count of how many\n\n2D elements are in our file. This corisponds to the number of elements in our tri mesh.\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.hex_vtk_order","page":"Reference","title":"StartUpDG.hex_vtk_order","text":"hex_vtk_order(corner_verts, order, dim, skip = false)\n\nCompute 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\n\nInspired by: https://github.com/ju-kreber/paraview-scripts/blob/master/node_ordering.py\n\nVTK node numbering of a hexagon:\n\n 8+------+7\n /| /|\n / | / |\n 5+------+6 |\n\nz | 4+–-|–+3 | y | / | / |/ |/ |/ 0–> x 1+–––+2\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.hybridized_SBP_operators-Tuple{Any}","page":"Reference","title":"StartUpDG.hybridized_SBP_operators","text":"function hybridized_SBP_operators(rd::RefElemData{DIMS})\n\nConstructs hybridized SBP operators given a RefElemData. Returns operators Qrsth..., VhP, Ph.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.inverse_trace_constant-Tuple{RefElemData{1, Line, <:Polynomial}}","page":"Reference","title":"StartUpDG.inverse_trace_constant","text":"function inverse_trace_constant(rd::RefElemData)\n\nReturns 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).\n\nCan be used to estimate dependence of maximum stable timestep on degree of approximation.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.make_periodic-Union{Tuple{MeshData{Dim}}, Tuple{Dim}, Tuple{MeshData{Dim}, Bool}} where Dim","page":"Reference","title":"StartUpDG.make_periodic","text":"make_periodic(md::MeshData{Dim}, is_periodic...) where {Dim}\nmake_periodic(md::MeshData{Dim}, is_periodic = ntuple(x->true,Dim)) where {Dim}\nmake_periodic(md::MeshData, is_periodic = true)\n\nReturns 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\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.n_verts_between-Tuple{Any, Any, Any}","page":"Reference","title":"StartUpDG.n_verts_between","text":"n_verts_between(n, from, to, dim)\n\nCompute 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\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.quad_vtk_order","page":"Reference","title":"StartUpDG.quad_vtk_order","text":"quad_vtk_order(corner_verts, order, dim, skip = false)\n\nCompute 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\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.read_Gmsh_2D-Tuple{String, Vararg{Any}}","page":"Reference","title":"StartUpDG.read_Gmsh_2D","text":"read_Gmsh_2D(filename, args...)\n\nReads a 2D triangular Gmsh file. Mesh formats 2.2 and 4.1 supported. Returns (VX, VY), EToV. \n\nExamples\n\nVXY, EToV = read_Gmsh_2D(\"eulerSquareCylinder2D.msh\") # v2.2 file format\nVXY, EToV = read_Gmsh_2D(\"test/testset_Gmsh_meshes/periodicity_mesh_v4.msh\") # v4.1 file format\n\n# if MeshImportOptions.grouping=true, then a third variable `grouping` is returned\nVXY, EToV, grouping = read_Gmsh_2D(\"test/testset_Gmsh_meshes/periodicity_mesh_v4.msh\", MeshImportOptions(true, false))\nVXY, EToV, grouping = read_Gmsh_2D(\"test/testset_Gmsh_meshes/periodicity_mesh_v4.msh\", true) # same as above\n\nSee also\n\nhttps://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format\nhttps://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format-version-2-0028Legacy0029\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.read_Gmsh_2D_v2-Tuple{String}","page":"Reference","title":"StartUpDG.read_Gmsh_2D_v2","text":"read_Gmsh_2D_v2(filename)\n\nReads triangular GMSH 2D file format 2.2 0 8. Returns (VX, VY), EToV.\n\nExamples\n\nVXY, EToV = read_Gmsh_2D_v2(\"eulerSquareCylinder2D.msh\")\n\nhttps://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format-version-2-0028Legacy0029\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.read_Gmsh_2D_v4","page":"Reference","title":"StartUpDG.read_Gmsh_2D_v4","text":"For brevity when grouping is the only supported feature.\n\nexample: VXY, EToV, grouping = read_Gmsh_2D_v4(\"file.msh\",true)\nexample: VXY, EToV = read_Gmsh_2D_v4(\"file.msh\",false)\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.read_Gmsh_2D_v4-Tuple{String, MeshImportOptions}","page":"Reference","title":"StartUpDG.read_Gmsh_2D_v4","text":"function read_Gmsh_2D_v4(filename, options)\n\nreads triangular GMSH 2D .msh files.\n\nOutput\n\nThis depends on if grouping is opted for or not\n\nreturns: (VX, VY), EToV\nreturns: (VX, VY), EToV, grouping\n\nSupported formats and features:\n\nversion 4.1 'physical group support 'remap group ids\n\ngrouping application\n\nWhen 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.\n\nVXY, EToV = read_Gmsh_2D_v4(\"eulerSquareCylinder2D.msh\")\nVXY, EToV = read_Gmsh_2D_v4(\"eulerSquareCylinder2D.msh\",false)\nVXY, EToV, grouping = read_Gmsh_2D_v4(\"eulerSquareCylinder2D.msh\", true)\n\noption = MeshImportOption(true)\nVXY, EToV, grouping = read_Gmsh_2D_v4(\"eulerSquareCylinder2D.msh\", option)\n\nhttps://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format\n\nNotes: the version 4 format has a more detailed block data format this leads to more complicated parser.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.refine","page":"Reference","title":"StartUpDG.refine","text":"function refine(triout, h, href = h/2)\n\nRefinement of a previous mesh given the current mesh size h. Preserves boundary/volume tags.\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.remap_element_grouping-Tuple{Vector{Int64}}","page":"Reference","title":"StartUpDG.remap_element_grouping","text":"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\n\nExample output\n\nremapelementgrouping([16,16,17,17]) -> [1,1,2,2]\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.sort_by_axis-Tuple{Any}","page":"Reference","title":"StartUpDG.sort_by_axis","text":"sort_by_axis(corner_verts)\n\nGiven the points 'corner_verts' sort them in a lexicographical order and return the permutated points. \n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.sparse_low_order_SBP_operators-Union{Tuple{RefElemData{NDIMS, ElemShape} where ElemShape<:AbstractElemShape{NDIMS}}, Tuple{NDIMS}} where NDIMS","page":"Reference","title":"StartUpDG.sparse_low_order_SBP_operators","text":"function sparse_low_order_SBP_operators(rd)\n\nConstructs sparse low order SBP operators given a RefElemData. Returns operators Qrst..., E ≈ Vf * Pq that satisfy a generalized summation-by-parts (GSBP) property:\n\n `Q_i + Q_i^T = E' * B_i * E`\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.tag_boundary_faces-Tuple{Any, Nothing}","page":"Reference","title":"StartUpDG.tag_boundary_faces","text":"function tag_boundary_faces(md, boundary_name::Symbol = :entire_boundary)\nfunction tag_boundary_faces(md, boundary_list::Dict{Symbol, <:Function})\n\nWhen called without arguments, just returns Dict(:entire_boundary => boundary_faces)`.\n\nExample usage: \n\njulia> rd = RefElemData(Tri(), N=1)\njulia> md = MeshData(uniform_mesh(Tri(), 2)..., rd)\njulia> on_bottom_boundary(x, y, tol = 1e-13) = abs(y+1) < tol\njulia> on_top_boundary(x, y, tol = 1e-13) = abs(y-1) < tol\njulia> tag_boundary_faces(Dict(:bottom => on_bottom_boundary,\n :top => on_top_boundary), md)\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.tag_boundary_faces-Tuple{Triangulate.TriangulateIO, RefElemData{2, Tri}, MeshData{2}, NamedTuple}","page":"Reference","title":"StartUpDG.tag_boundary_faces","text":"function tag_boundary_faces(triout::TriangulateIO,\n rd::RefElemData{2,Tri}, md::MeshData{2},\n boundary_list::Union{NamedTuple,Dict{Symbol,Int}})\n\nHere, 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.\n\nExample usage:\n\njulia> using Triangulate, StartUpDG\njulia> triout = scramjet()\njulia> rd = RefElemData(Tri(),N=1)\njulia> md = MeshData(triangulateIO_to_VXYEToV(triout)...,rd)\njulia> tag_boundary_faces(triout,rd,md, Dict(:wall=>1, :inflow=>2, :outflow=>3))\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.triangle_vtk_order","page":"Reference","title":"StartUpDG.triangle_vtk_order","text":"triangle_vtk_order(corner_verts, order, dim, skip = false)\n\nCompute 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\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.triangulateIO_to_VXYEToV-Tuple{Triangulate.TriangulateIO}","page":"Reference","title":"StartUpDG.triangulateIO_to_VXYEToV","text":"function triangulateIO_to_VXYEToV(triout::TriangulateIO)\n\nComputes VX,VY,EToV from a TriangulateIO object.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.type_to_vtk-Tuple{Hex}","page":"Reference","title":"StartUpDG.type_to_vtk","text":"type_to_vtk(elem::Hex)\nreturn the VTK-type\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.type_to_vtk-Tuple{Quad}","page":"Reference","title":"StartUpDG.type_to_vtk","text":"type_to_vtk(elem::Quad)\nreturn the VTK-type\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.type_to_vtk-Tuple{Tri}","page":"Reference","title":"StartUpDG.type_to_vtk","text":"type_to_vtk(elem::Tri)\nreturn the VTK-type\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.type_to_vtk-Tuple{Wedge}","page":"Reference","title":"StartUpDG.type_to_vtk","text":"type_to_vtk(elem::Wedge)\nreturn the VTK-type\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.uniform_mesh-Tuple{Line, Any}","page":"Reference","title":"StartUpDG.uniform_mesh","text":"uniform_mesh(elem::Line,Kx)\nuniform_mesh(elem::Tri,Kx,Ky)\nuniform_mesh(elem::Quad,Kx,Ky)\nuniform_mesh(elem::Hex,Kx,Ky,Kz)\nuniform_mesh(elem, K)\n\nUniform Kx (by Ky by Kz) mesh on -11^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.\n\nK can also be specified using a keyword argument K1D, e.g., uniform_mesh(elem; K1D = 16).\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.vtk_order-Tuple{Hex, Any}","page":"Reference","title":"StartUpDG.vtk_order","text":"vtk_order(elem::Hex, order)\n\nConstruct 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.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.vtk_order-Tuple{Quad, Any}","page":"Reference","title":"StartUpDG.vtk_order","text":"vtk_order(elem::Quad, order)\n\nConstruct 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\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.vtk_order-Tuple{Tri, Any}","page":"Reference","title":"StartUpDG.vtk_order","text":"vtk_order(elem::Tri, order)\n\nConstruct 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\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.vtk_order-Tuple{Wedge, Any}","page":"Reference","title":"StartUpDG.vtk_order","text":"vtk_order(elem::Wedge, order)\n\nConstruct all node-points of a VTKLAGRANGEWEDGE of order order. The corner-nodes are given by the reference-wedge used by StartUpDG\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.wedge_vtk_order-Tuple{Any, Any, Any}","page":"Reference","title":"StartUpDG.wedge_vtk_order","text":"wedge_vtk_order(corner_verts, order, dim)\n\nCompute 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\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#Triangulate.triangulate","page":"Reference","title":"Triangulate.triangulate","text":"function Triangulate.triangulate(triin::TriangulateIO, maxarea, minangle=20)\n\nConvenience routine to avoid writing out @sprintf each time. Returns a TriangulateIO object.\n\n\n\n\n\n","category":"function"},{"location":"tstep_usage/#Time-stepping","page":"Timestepping","title":"Time-stepping","text":"","category":"section"},{"location":"tstep_usage/","page":"Timestepping","title":"Timestepping","text":"For convenience, we include a commonly used explicit Runge-Kutta scheme. For more advanced time-stepping functionality, we recommend OrdinaryDiffEq.jl. ","category":"page"},{"location":"tstep_usage/#Carpenter-and-Kennedy's-(4,5)-method","page":"Timestepping","title":"Carpenter and Kennedy's (4,5) method","text":"","category":"section"},{"location":"tstep_usage/","page":"Timestepping","title":"Timestepping","text":"ck45() returns coefficients for a low-storage Runge-Kutta method.","category":"page"},{"location":"tstep_usage/#Example-usage:","page":"Timestepping","title":"Example usage:","text":"","category":"section"},{"location":"tstep_usage/","page":"Timestepping","title":"Timestepping","text":"using Plots\nusing StartUpDG\n\n# Brusselator\nB = 3\nf(y1, y2) = [1 + y1^2 * y2 - (B+1) * y1, B * y1 - y1^2 * y2]\nf(Q) = f(Q[1], Q[2])\np,q = 1.5, 3.0\nQ = [p; q]\n\ndt = .1\nFinalTime = 20\n\nres = zero.(Q) # init RK residual\nrk4a, rk4b, rk4c = ck45()\nNsteps = ceil(FinalTime / dt)\ndt = FinalTime / Nsteps\n\nplot() # init plot\nfor i = 1:Nsteps\n global res # yes, I know...this is just for simplicty\n for INTRK = 1:5\n rhsQ = f(Q)\n @. res = rk4a[INTRK] * res + dt * rhsQ # i = RK stage\n @. Q = Q + rk4b[INTRK] * res\n end\n scatter!([i*dt;i*dt],Q)\nend\ndisplay(plot!(leg=false))","category":"page"},{"location":"more_meshes/#Additional-mesh-types","page":"Additional mesh types","title":"Additional mesh types","text":"","category":"section"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"In addition to more \"standard\" mesh types, StartUpDG.jl also has experimental support for hybrid and cut-cell meshes. Both are currently restricted to two dimensional domains.","category":"page"},{"location":"more_meshes/#Hybrid-meshes","page":"Additional mesh types","title":"Hybrid meshes","text":"","category":"section"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"warning: Experimental implementation\nThis is an experimental feature and may change in future releases.","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"There is initial support for hybrid meshes in StartUpDG.jl. The following is a short example where we interpolate a polynomial and compute its derivative.","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"rds = RefElemData((Tri(), Quad()), N = 3)\n\n# Simple hybrid mesh for testing\n# 1 7______8______9\n# | | 3 / |\n# | 4 | / 5 |\n# 0 4 ---- 5 ---- 6 \n# | | |\n# | 1 | 2 |\n# -1 1 ---- 2 ---- 3\n# -1 0 1\nVX = [-1; 0; 1; -1; 0; 1; -1; 0; 1]\nVY = [-1; -1; -1; 0; 0; 0; 1; 1; 1]\nEToV = [[1 2 4 5], [2 3 5 6], [5 8 9], [4 5 7 8], [9 6 5]]\n\nmd = MeshData(VX, VY, EToV, rds)\n\n# test that the local derivatives of a polynomial recover the exact derivative\n(; x, y ) = md\nu = @. x^3 - x^2 * y + 2 * y^3\ndudx = @. 3 * x^2 - 2 * x * y\n\n# compute local derivatives\n(; rxJ, sxJ, J ) = md\ndudr, duds = similar(md.x), similar(md.x)\ndudr.Quad .= rds[Quad()].Dr * u.Quad\nduds.Quad .= rds[Quad()].Ds * u.Quad\ndudr.Tri .= rds[Tri()].Dr * u.Tri\nduds.Tri .= rds[Tri()].Ds * u.Tri\n\n@show norm(@. dudx - (rxJ * dudr + sxJ * duds) / J) # should be O(1e-14)","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"The main difference in the representation of hybrid meshes compared with standard MeshData objects is the use of NamedArrayPartition arrays as storage for the geometric coordinates. These arrays have \"fields\" corresponding to the element type, for example","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"md.x.Tri\nmd.x.Quad","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"but can still be indexed as linear arrays. ","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"The mapP field behaves similarly. If we interpolate the values of u for each element type to surface quadrature nodes, we can use mapP to linearly index into the array to find neighbors. ","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"uf = similar(md.xf)\nuf.Quad .= rds.Quad.Vf * u.Quad\nuf.Tri .= rds.Tri.Vf * u.Tri\nuf[md.mapP] # this returns the exterior node values","category":"page"},{"location":"more_meshes/#Cut-Meshes","page":"Additional mesh types","title":"Cut Meshes","text":"","category":"section"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"warning: Experimental implementation\nThis is an experimental feature and may change in future releases.","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"Initial support for cut-cell meshes is available via PathIntersections.jl. By passing in a tuple of curves (defined as parametrized functions of one coordinate, see PathIntersections.jl documentation for more detail), StartUpDG.jl can compute a MeshData for a cut-cell mesh. ","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"circle = PresetGeometries.Circle(R=0.33, x0=0, y0=0)\ncells_per_dimension_x, cells_per_dimension_y = 4, 4\n\nrd = RefElemData(Quad(), N=3)\nmd = MeshData(rd, (circle, ), cells_per_dimension_x, cells_per_dimension_y; precompute_operators=true)","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"The interpolation points on cut cells md.x.cut are determined from sampled points and a pivoted QR decomposition. The quadrature points on cut cells md.xq.cut are determined similarly. However, the cut-cell quadrature weights md.wJq.cut are not currently positive. The optional keyword argument precompute_operators specifies whether to precompute differentiation, face interpolation, mass, and lifting matrices for each cut cell. If precompute_operators=true, these are stored in md.mesh_type.cut_cell_operators. ","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"As with hybrid meshes, the nodal coordinates md.x, md.y are NamedArrayPartitions with cartesian and cut fields. For example, md.x.cartesian and md.x.cut are the x-coordinates of the Cartesian and cut cells, respectively. Likewise, md.mapP indexes linearly into the array of face coordinates and specifies exterior node indices. For example, we can interpolate a function to face nodes and compute exterior values via the following code:","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"(; x, y) = md\nu = @. x^rd.N - x * y^(rd.N-1) - x^(rd.N-1) * y + y^rd.N # some random function \n\n# interpolate the solution to face nodes\nuf = similar(md.xf) \nuf.cartesian .= rd.Vf * u.cartesian\nfor e in 1:size(md.x.cut, 2)\n ids = md.mesh_type.cut_face_nodes[e]\n Vf = md.mesh_type.cut_cell_operators.face_interpolation_matrices[e]\n uf.cut[ids] .= Vf * u.cut[:, e]\nend\n\nuf[md.mapP] # these are \"exterior\" values for each entry of `uf`","category":"page"},{"location":"ex_dg_deriv/#Example:-approximating-derivatives-using-DG","page":"Example: computing DG derivatives","title":"Example: approximating derivatives using DG","text":"","category":"section"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"MeshData can be used to compute DG derivatives. Suppose f is a differentiable function and the domain Omega can be decomposed into non-overlapping elements D^k. The approximation of fracpartial fpartial x can be approximated using the following formulation: find piecewise polynomial u such that for all piecewise polynomials v","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"int_Omega u v = sum_k left( int_D^k fracpartial upartial xv + int_partial D^k frac12 lefturightn_x v right)","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"Here, lefturight = u^+ - u denotes the jump across an element interface, and n_x is the x-component of the outward unit normal on D^k.","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"Discretizing the left-hand side of this formulation yields a mass matrix. Inverting this mass matrix to the right hand side yields the DG derivative. We show how to compute it for a uniform triangular mesh using MeshData and StartUpDG.jl.","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"We first construct the triangular mesh and initialize md::MeshData.","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"using StartUpDG\nusing Plots\n\nN = 3\nK1D = 8\nrd = RefElemData(Tri(), N)\nVXY, EToV = uniform_mesh(Tri(), K1D)\nmd = MeshData(VXY, EToV, rd)","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"We can approximate a function f(x y) using interpolation","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"f(x, y) = exp(-5 * (x^2 + y^2)) * sin(1 + pi*x) * sin(2 + pi*y)\n(; x, y) = md\nu = @. f(x, y)","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"or using quadrature-based projection","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"(; Pq ) = rd\n(; x, y, xq, yq ) = md\nu = Pq * f.(xq, yq)","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"We can use scatter in Plots.jl to quickly visualize the approximation. This is not intended to create a high quality image (see other libraries, e.g., Makie.jl,VTK.jl, or Triplot.jl for publication-quality images).","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"(; Vp ) = rd\nxp, yp, up = Vp * x, Vp * y, Vp * u # interp to plotting points\nscatter(xp, yp, uxp, zcolor=uxp, msw=0, leg=false, ratio=1, cam=(0, 90))","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"Both interpolation and projection create a matrix u of size N_p times K which contains coefficients (nodal values) of the DG polynomial approximation to f(x y). We can approximate the derivative of f(x y) using the DG derivative formulation","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"function dg_deriv_x(u, rd::RefElemData, md::MeshData)\n (; Vf, Dr, Ds, LIFT ) = rd\n (; rxJ, sxJ, J, nxJ, mapP ) = md\n uf = Vf * u\n ujump = uf[mapP] - uf\n\n # derivatives using chain rule + lifted flux terms\n ux = rxJ .* (Dr * u) + sxJ .* (Ds * u) \n dudxJ = ux + LIFT * (.5 * ujump .* nxJ)\n\n return dudxJ ./ J\nend","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"We can visualize the result as follows:","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"dudx = dg_deriv_x(u, rd, md)\nuxp = Vp * dudx\nscatter(xp, yp, uxp, zcolor=uxp, msw=0, leg=false, ratio=1, cam=(0,90))","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"Plots of the polynomial approximation u(xy) and the DG approximation of fracpartial upartial x are given below","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"(Image: u) (Image: dudx) ⠀","category":"page"},{"location":"#Overview","page":"Home","title":"Overview","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package contains routines to initialize reference element operators, physical mesh arrays, and connectivity arrays for nodal DG methods. The codes roughly based on Nodal Discontinuous Galerkin Methods by Hesthaven and Warburton (2007).","category":"page"},{"location":"#A-short-example","page":"Home","title":"A short example","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"using StartUpDG\n\n# polynomial degree and mesh size\nN = 3\nK1D = 8\n\n# init ref element and mesh\nrd = RefElemData(Tri(), N)\nVXY, EToV = uniform_mesh(Tri(), K1D)\nmd = MeshData(VXY, EToV, rd)\n\n# Define a function by interpolation\n(; x, y ) = md\nu = @. exp(-10 * (x^2 + y^2))\n\n# Compute derivatives using geometric mapping + chain rule\n(; Dr, Ds ) = rd\n(; rxJ, sxJ, J ) = md\ndudx = (rxJ .* (Dr * u) + sxJ .* (Ds * u)) ./ J","category":"page"}] +[{"location":"MeshData/#MeshData-type","page":"MeshData","title":"MeshData type","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"MeshData includes fields such as","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"xyz::NTuple{Dim, ...}: nodal interpolation points mapped to physical elements. All elements of xyz are N_p times N_rm elements matrices, where N_p are the number of nodal points on each element.\nxyzq::NTuple{Dim, ...}, wJq: volume quadrature points/weights mapped to physical elements. All elements these tuples are N_q times N_rm elements matrices, where N_q is the number of quadrature points on each element.\nxyzf::NTuple{Dim, ...}: face quadrature points mapped to physical elements. All elements of xyz are N_f times N_rm elements matrices, where N_f is the number of face points on each element.\nmapP, mapB: indexing arrays for inter-element node connectivity (mapP) and for extracting boundary nodes from the list of face nodes xyzf (mapB). mapP is a matrix of size N_f times N_rm elements, while the length of mapB is the total number of nodes on the boundary.\nrstxyzJ::SMatrix{Dim, Dim}: volume geometric terms G_ij = fracpartial x_ipartial hatx_j. Each element of rstxyzJ is a matrix of size N_p times N_rm elements.\nJ, Jf: volume and surface Jacobians evaluated at interpolation points and surface quadrature points, respectively. J is a matrix of size N_p times N_rm elements, while Jf is a matrix of size N_f times N_rm elements. \nnxyz::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 N_rm elements. ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"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. ","category":"page"},{"location":"MeshData/#Setting-up-md::MeshData","page":"MeshData","title":"Setting up md::MeshData","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"The MeshData struct contains data for high order DG methods useful for evaluating DG formulations in a matrix-free fashion.","category":"page"},{"location":"MeshData/#Generating-unstructured-meshes","page":"MeshData","title":"Generating unstructured meshes","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"For convenience, simple uniform meshes are included in with StartUpDG.jl via uniform_mesh","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"using StartUpDG\nnum_cells_x, num_cells_y, num_cells_z = 4, 2, 8\n(VX,), EToV = uniform_mesh(Line(), num_cells_x)\n(VX, VY), EToV = uniform_mesh(Tri(), num_cells_x, num_cells_y)\n(VX, VY), EToV = uniform_mesh(Quad(), num_cells_x, num_cells_y)\n(VX, VY, VZ), EToV = uniform_mesh(Tet(), num_cells_x, num_cells_y, num_cells_z)\n(VX, VY, VZ), EToV = uniform_mesh(Pyr(), num_cells_x, num_cells_y, num_cells_z)\n(VX, VY, VZ), EToV = uniform_mesh(Wedge(), num_cells_x, num_cells_y, num_cells_z)\n(VX, VY, VZ), EToV = uniform_mesh(Hex(), num_cells_x, num_cells_y, num_cells_z)","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"The uniform triangular mesh is constructed by creating a uniform quadrilateral mesh then bisecting each quad into two triangles. Wedge meshes are constructed similarly. Tet meshes are constructed by dividing each hexahedron into 5 tetrahedral elements. Pyramid meshes are constructed by dividing each hexahedron into 6 pyramids. ","category":"page"},{"location":"MeshData/#Initializing-high-order-DG-mesh-data","page":"MeshData","title":"Initializing high order DG mesh data","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"Given unstructured mesh information (tuple of vertex coordinates VXYZ and index array EToV) high order DG mesh data can be constructed as follows:","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"md = MeshData(VXYZ, EToV, rd)","category":"page"},{"location":"MeshData/#Enforcing-periodic-boundary-conditions","page":"MeshData","title":"Enforcing periodic boundary conditions","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"Periodic boundary conditions can be enforced via the is_periodic keyword argument ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"md_periodic = MeshData((VX, VY), EToV, rd; is_periodic=true) # periodic in both x and y coordinates\nmd_periodic_x = MeshData((VX, VY), EToV, rd; is_periodic=(true, false)) # periodic in x direction, but not y","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"or by calling make_periodic, which returns a new MeshData instance","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"md = MeshData((VX, VY), EToV, rd)\nmd_periodic = make_periodic(md) # periodic in both x and y coordinates\nmd_periodic_x = make_periodic(md, true, false) # periodic in x direction, but not y","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"In either case, the MeshData indexing arrays fields mapP,mapB, and FToF are modified to account for periodicity.","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"One can check which dimensions are periodic via the is_periodic field of MeshData. For example, the md_periodic_x example above gives","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"julia> md_periodic_x.is_periodic\n(true, false)","category":"page"},{"location":"MeshData/#Creating-curved-meshes","page":"MeshData","title":"Creating curved meshes","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"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 again with new x, y coordinates:","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"md = MeshData((VX, VY), EToV, rd)\n(; x, y ) = md\n# <-- code to modify high order nodes (x,y)\nmd_curved = MeshData(rd, md, x, y)","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"MeshData(rd, md, x, y) and MeshData(rd, md, x, y, z) are implemented for 2D and 3D, though this is not currently implemented in 1D.","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"More generally, one can create a copy of a MeshData with certain fields modified by using @set or setproperties from Setfield.jl.","category":"page"},{"location":"MeshData/#Unstructured-and-pre-defined-triangular-meshes-using-Triangulate.jl","page":"MeshData","title":"Unstructured and pre-defined triangular meshes using Triangulate.jl","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"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:","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"meshIO = triangulate_domain(Scramjet())\n(VX, VY), EToV = triangulateIO_to_VXYEToV(meshIO)\nrd = RefElemData(Tri(), 7)\nmd = MeshData((VX, VY), EToV, rd)","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"A quick plot of the face nodes via ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"using Plots\nscatter(vec.(md.xyzf)..., msw=0, ms=1, aspect_ratio=:equal, ylims=(0,2), leg=false)","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"shows the following figure (Image: u)","category":"page"},{"location":"MeshData/#Unstructured-curved-quadrilateral-and-hexahedral-meshes-using-HOHQMesh.jl","page":"MeshData","title":"Unstructured curved quadrilateral and hexahedral meshes using HOHQMesh.jl","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"StartUpDG.jl also reads in .mesh files generated by HOHQMesh.jl. The following code constructs a MeshData which represents a curved quadrilateral mesh generated by HOHQMesh.jl. ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"using StartUpDG\nrd = RefElemData(Quad(), 4)\nhmd = read_HOHQMesh(\"test/testset_HOHQMesh_meshes/easy_example.mesh\", Quad())\nmd = MeshData(hmd, rd)","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"We can visualize the mesh using ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"using Plots\nplot(rd, md) # can also use `plot(MeshPlotter(rd, md))`","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"which yields the following figure:","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"(Image: u)","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"The boundary faces are also automatically tagged with the labels provided in the HOHQMesh file. Each boundary tag and the faces that lie on it are stored in md.mesh_type.boundary_faces. ","category":"page"},{"location":"MeshData/#Tagging-boundary-faces-and-boundary-nodes","page":"MeshData","title":"Tagging boundary faces and boundary nodes","text":"","category":"section"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"One can \"tag\" boundary faces (or boundary nodes) by specifying boolean functions which evaluate to true if a point is on a given boundary segment. ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"rd = RefElemData(Tri(), N=3)\nmd = MeshData(uniform_mesh(Tri(), 1)..., rd)\non_bottom_boundary(point, tol=1e-13) = abs(point[2] + 1) < tol # point = (x,y)\non_top_boundary(point, tol=1e-13) = abs(point[2] - 1) < tol \n\nboundary_face_dict = tag_boundary_faces(md, Dict(:bottom => on_bottom_boundary, :top => on_top_boundary))\nboundary_node_dict = tag_boundary_nodes(rd, md, Dict(:bottom => on_bottom_boundary, :top => on_top_boundary))","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"You can also specify a list of boundaries using NamedTuples ","category":"page"},{"location":"MeshData/","page":"MeshData","title":"MeshData","text":"boundary_face_dict = tag_boundary_faces(md, (; :bottom=>on_bottom_boundary,:top=>on_top_boundary))\nboundary_node_dict = tag_boundary_nodes(rd, md, (; :bottom=>on_bottom_boundary,:top=>on_top_boundary))","category":"page"},{"location":"RefElemData/#RefElemData-type","page":"RefElemData","title":"RefElemData type","text":"","category":"section"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"RefElemData contains the following fields","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"element_type::AbstractElemShape: element shape. Line, Tri, Quad, Hex, Wedge, Pyr, Tet currently supported.\napproximation_type: approximation type. Defaults to Polynomial(), but SBP() is also supported (see RefElemData based on SBP finite differences).\nNfaces: number of faces on a given type of reference element.\nfv: list of vertices defining faces, e.g., [1,2], [2,3], [3,1] for a triangle\nFmask: indices of interpolation nodes which lie on the faces\nrst::NTuple{Dim, ...}: tuple of vectors of length N_p, each of which contains coordinates of degree N optimized polynomial interpolation points.\nrstq::NTuple{Dim, ...},wq, Vq: tuple of volume quadrature points, vector of weights, and quadrature interpolation matrix. Each element of rstq and wq are vectors of length N_q, and Vq is a matrix of size N_q times N_p.\nN_{\\rm plot}: the degree which determines the number of plotting points N_prm plot.\nrstp::NTuple{Dim, ...}, Vp: tuple of plotting points and plotting interpolation matrix. Each element of rstp is a vector of length N_prm plot, and Vp is a matrix of size N_prm plot times N_p.\nrstf::NTuple{Dim, ...},wf, Vf: tuple of face quadrature points, weights, and face interpolation matrix. Each element of rstf and wf are vectors of length N_f, and Vf is a matrix of size N_f times N_p.\nnrstJ::NTuple{Dim, ...}: tuple of outward reference normals, scaled by the face Jacobian. Each element is a vector of length N_f.\nM: mass matrix computed using quadrature. Size N_p times N_p\nPq: quadrature-based L^2 projection matrix. Size N_p times N_q.\nDrst::NTuple{Dim, ...}, LIFT: differentiation and lifting matrices. Differentiation matrices are size N_p times N_p while lift matrices are size N_ptimes N_f.","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"This list is incomplete; other fields are stored or accessible but currently only used for internal computations.","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"Mass, differentiation, lifting, and interpolation matrices can be specialized. For example, these matrices are dense Matrix{T} type for lines and triangles, but could also be stored as sparse matrices for quadrilaterals and hexahedra.","category":"page"},{"location":"RefElemData/#Setting-up-rd::RefElemData","page":"RefElemData","title":"Setting up rd::RefElemData","text":"","category":"section"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"The struct rd::RefElemData contains data for a given element type. All common reference elements are supported: Line, Tri, Quad, Tet, Pyr, Wedge, and Hex.","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"To initalize a RefElemData, just specify the element type and polynomial degree.","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"N = 3\n\n# 1D elements \nrd = RefElemData(Line(), N)\n\n# 2D elements\nrd = RefElemData(Tri(), N)\nrd = RefElemData(Quad(), N)\n\n# 3D elements\nrd = RefElemData(Tet(), N)\nrd = RefElemData(Pyr(), N)\nrd = RefElemData(Wedge(), N)\nrd = RefElemData(Hex(), N)","category":"page"},{"location":"RefElemData/#Specifying-different-quadrature-rules.","page":"RefElemData","title":"Specifying different quadrature rules.","text":"","category":"section"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"By default, RefElemData initializes volume and surface quadrature rules to be the minimum rules which exactly integrate the unweighted volume and surface mass matrices. If different quadrature rules are desired, they can be specified as follows: ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"N = 3\n\n# create degree N tensor product Gauss-Lobatto rule\nr1D, w1D = gauss_lobatto_quad(0, 0, N)\nrq, sq = vec.(StartUpDG.meshgrid(r1D))\nwr, ws = vec.(StartUpDG.meshgrid(w1D))\nwq = @. wr * ws\n\nrd = RefElemData(Quad(), N; quad_rule_vol = (rq, sq, wq), \n quad_rule_face = (r1D, w1D))","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"This results in a DG spectral element method (DG-SEM) discretization, with a diagonal lumped mass matrix and differentiation matrices which satisfy a summation-by-parts property. ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"By default, RefElemData is constructed for a nodal basis (in order to facilitate curved meshes, connectivity, etc). The interpolation nodes are computed using an interpolatory version of the warp-and-blend procedure. ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"note: Note\nWhile specifying the quadrature rule changes the discretization, it is not reflected in the RefElemData type and thus cannot be specialized on. The following constructors produce RefElemData where the quadrature structure is reflected in the type parameters:rd = RefElemData(Hex(), Polynomial(TensorProductQuadrature(quad_nodes(Line(), N+1)), N)) # tensor product quadrature rules\nrd = RefElemData(Quad(), Polynomial{Gauss}(), N) # (N+1)^d point tensor product Gauss quadrature","category":"page"},{"location":"RefElemData/#RefElemData-based-on-SBP-finite-differences","page":"RefElemData","title":"RefElemData based on SBP finite differences","text":"","category":"section"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"It is also possible to construct a RefElemData based on multi-dimensional SBP finite difference operators. These utilize nodes constructed by Tianheng Chen and Chi-Wang Shu, Ethan Kubatko, and Jason Hicken.","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"Some examples:","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"N = 3\nrd = RefElemData(Quad(), SBP(), N) # defaults to SBP{TensorProductLobatto}\nrd = RefElemData(Quad(), SBP{TensorProductLobatto}(), N) \n\nrd = RefElemData(Hex(), SBP(), N) # defaults to SBP{TensorProductLobatto}\nrd = RefElemData(Hex(), SBP{TensorProductLobatto}(), N) \n\nrd = RefElemData(Tri(), SBP(), N) # defaults to SBP{Kubatko{LobattoFaceNodes}}\nrd = RefElemData(Tri(), SBP{Hicken}(), N) \nrd = RefElemData(Tri(), SBP{Kubatko{LobattoFaceNodes}}(), N) \nrd = RefElemData(Tri(), SBP{Kubatko{LegendreFaceNodes}}(), N) ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"Quadrature rules of both degree 2*N-1 (up to N = 6) and 2*N (up to N = 4) are supported on triangles. For Line, Quad, and Hex elements, RefElemData(..., SBP(), N) is the same as the RefElemData for a DG-SEM discretization, though some fields are specialized for the SBP type. These SBP-based RefElemData objects can also be used to initialize a mesh (for example, md = MeshData(uniform_mesh(rd.element_type, 4)..., rd)). ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"On triangles, we have the following SBP types with the following properties:","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"SBP{Kubatko{LobattoFaceNodes}}: degree 2N-1 accurate quadrature rules with N+2 Lobatto nodes on each face. Nodes for N=4: ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"(Image: klobatto4)","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"SBP{Kubatko{LegendreFaceNodes}}: degree 2N-1 accurate quadrature rules with N+1 Legendre nodes on each face. For N = 1,...,4, these are the same as the nodes constructed by Chen and Shu. Nodes for N=4:","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"(Image: klegendre4)","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"SBP{Hicken}: degree 2N accurate quadrature rules with N+2 Lobatto nodes on each face. Nodes for N=4:","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"(Image: hicken4)","category":"page"},{"location":"RefElemData/#Tensor-product-RefElemData-on-wedge-elements","page":"RefElemData","title":"Tensor product RefElemData on wedge elements","text":"","category":"section"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"warning: Experimental implementation\nThis is an experimental feature and may change in future releases.","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"There is experimental support for RefElemDatas created from tensor products of triangular and 1D RefElemData objects. ","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"line = RefElemData(Line(), N_line)\ntri = RefElemData(Tri(), N_tri)\nrd = RefElemData(Wedge(), TensorProductWedge(tri, line))","category":"page"},{"location":"RefElemData/","page":"RefElemData","title":"RefElemData","text":"This new rd::RefElemData can then be used to create a wedge-based MeshData. The individual RefElemData objects can be accessed from rd.approximation_type::TensorProductWedge. ","category":"page"},{"location":"conventions/#Background","page":"Background and conventions","title":"Background","text":"","category":"section"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"Most high order finite element methods rely on a decomposition of a domain into a mesh of \"elements\" (e.g., triangles or quadrilaterals in 2D, hexahedra or tetrahedra in 3D). Each \"physical\" element in a mesh is assumed to be the image of single \"reference\" element under some geometric mapping. Using the chain rule and changes of variables, one can evaluate integrals and derivatives using only operations on the reference element and some geometric mapping data. This transformation of operations on all elements to a single reference element make finite element methods efficient. ","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"We use the convention that coordinates on the reference element are r in 1D, r s in 2D, or r s t in 3D. Physical coordinates use the standard conventions x, x y, and x y z in 1D, 2D, and 3D. ","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"(Image: Mapping)","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"Derivatives of reference coordinates with respect to physical coordinates are abbreviated, e.g., fracpartial rpartial x = r_x. Additionally, J is used to denote the determinant of the Jacobian of the reference-to-physical mapping. ","category":"page"},{"location":"conventions/#Assumptions","page":"Background and conventions","title":"Assumptions","text":"","category":"section"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"We make a few simplifying assumptions about the mesh:","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"meshes are conforming (e.g., each face of an element is shared with at most one other element). \nthe geometric mapping from reference to physical elements is the same degree polynomial as the approximation space on the reference element (e.g., the mapping is isoparametric). ","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"Initial experimental support for hybrid, cut-cell, and non-conforming meshes in two dimensions is also available. Please see the corresponding test sets test/hybrid_mesh_tests.jl, test/cut_mesh_tests.jl, and noncon_mesh_tests.jl for examples. ","category":"page"},{"location":"conventions/#Code-conventions","page":"Background and conventions","title":"Code conventions","text":"","category":"section"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"StartUpDG.jl exports structs RefElemData{Dim, ElemShape, ...} (which contains data associated with the reference element, such as interpolation points, quadrature rules, face nodes, normals, and differentiation/interpolation/projection matrices) and MeshData{Dim} (which contains geometric data associated with a mesh). These are currently used for evaluating DG formulations in a matrix-free fashion. These structs contain fields similar to those in Globals1D, Globals2D, Globals3D in the NDG book codes. ","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"We use the following code conventions:","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"variables r, s, t and x, y, z correspond to values at nodal interpolation points. \nvariables ending in q (e.g., rq, sq,... and xq, yq,...) correspond to values at volume quadrature points. \nvariables ending in f (e.g., rf, sf,... and xf, yf,...) correspond to values at face quadrature points. \nvariables ending in p (e.g., rp, sp,...) correspond to equispaced plotting nodes.\nDr, Ds, Dt matrices are nodal differentiation matrices with respect to the r s t coordinates. For example, Dr * f.(r, s) approximates the derivative of f(r s) at nodal points. \nV matrices correspond to interpolation matrices from nodal interpolation points. For example, Vq interpolates to volume quadrature points, Vf interpolates to face quadrature points, Vp interpolates to plotting nodes. \ngeometric quantities in MeshData are stored as matrices of dimension textnumber of points per element times textnumber of elements.","category":"page"},{"location":"conventions/#Differences-from-the-codes-of-\"Nodal-Discontinuous-Galerkin-Methods\"","page":"Background and conventions","title":"Differences from the codes of \"Nodal Discontinuous Galerkin Methods\"","text":"","category":"section"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"The codes in StartUpDG.jl are based closely on the Matlab codes from the book \"Nodal Discontinuous Galerkin Methods\" by Hesthaven and Warburton (2008) (which we will refer to as \"NDG\"). However, there are some differences in order to allow for more general DG discretizations and enforce certain mathematical properties:","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"In NDG, Fmask extracts the interpolation nodes which lie on a face. These nodes are then used to compute interface fluxes. However, in StartUpDG.jl, we interpolate nodal values to values at face quadrature points via rd.Vf * u. These operations are equivalent if the interpolation nodes which lie on a face are co-located with quadrature points. Similarly, in NDG, the LIFT matrix maps face nodal points to volume nodal points. In StartUpDG.jl, the rd.LIFT matrix maps from face quadrature points to volume nodal points. \nin NDG, there are connectivity arrays vmapM and vmapP, which directly retrieve interface values from arrays of nodal values. In StartUpDG.jl, face interpolation nodes are not guaranteed to be co-located with face quadrature nodes, so we do not provide vmapM and vmapP. Instead, we expect the user to compute face values and use the md.mapM, md.mapP arrays to access interface values. \nin NDG, the mass matrix is computed exactly using the formula M = inv(VDM * VDM'), where VDM is the generalized Vandermonde matrix evaluated at nodal interpolation points. In StartUpDG.jl, the mass matrix is computed using quadrature. These are equivalent if the quadrature is exact for the integrands in the mass matrix (e.g., degree 2N polynomials for triangular or tetrahedral elements).\nin NDG, the geometric terms rx, sx, ry, sy, ... are computed and stored. In StartUpDG.jl, the scaled geometric terms md.rxJ, md.sxJ, md.ryJ, md.syJ, ... are computed, which enable us to enforce the metric identities on curved meshes. Similarly, NDG provides Fscale = sJ ./ J(Fmask, :), while StartUpDG.jl only provides md.Jf, which is equivalent to sJ. Fscale, as well as the NDG geometric terms and can be recovered by dividing by md.J. ","category":"page"},{"location":"conventions/","page":"Background and conventions","title":"Background and conventions","text":"Internally, NDG uses arrays EToE and EToF to compute the interface connectivity array mapP. StartUpDG.jl uses instead a face-to-face connectivity array FToF. However, EToE, EToF, and FToF are not typically required for the matrix-free explicit solvers targeted by this package. ","category":"page"},{"location":"index_refs/#Index","page":"Reference","title":"Index","text":"","category":"section"},{"location":"index_refs/","page":"Reference","title":"Reference","text":"","category":"page"},{"location":"index_refs/#Functions","page":"Reference","title":"Functions","text":"","category":"section"},{"location":"index_refs/","page":"Reference","title":"Reference","text":"Modules = [StartUpDG]","category":"page"},{"location":"index_refs/#StartUpDG.BoundaryTagPlotter","page":"Reference","title":"StartUpDG.BoundaryTagPlotter","text":"BoundaryTagPlotter(triout::TriangulateIO)\n\nPlot recipe to visualize boundary tags by color. Usage: plot(BoundaryTagPlotter(triout))\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.CurvedMesh","page":"Reference","title":"StartUpDG.CurvedMesh","text":"struct CurvedMesh{T}\n\nMesh type indicating that the mesh has been curved. Stores the original mesh type as a field.\n\nFields\n\noriginalmeshtype :: T\n\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.CutCellMesh","page":"Reference","title":"StartUpDG.CutCellMesh","text":"CutCellMesh is used in the MeshData field mesh_type for cut cell meshes.\n\nThe 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. \n\nThe 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. \n\nWe assume all cut elements have the same number of volume quadrature points (which is at least the dimension of a degree 2N polynomial space). \n\nThe field objects contains a tuple of the objects used to define the cut region.\n\nThe field cut_cell_operators contains optionally precomputed operators (mass, differntiation, face interpolation, and lifting operators). \n\nThe field cut_cell_data contains additional data from PathIntersections.\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.MeshData","page":"Reference","title":"StartUpDG.MeshData","text":"struct MeshData{Dim, Tv, Ti}\n\nMeshData: contains info for a high order piecewise polynomial discretization on an unstructured mesh. \n\nExample:\n\nN, K1D = 3, 2\nrd = RefElemData(Tri(), N)\nVXY, EToV = uniform_mesh(Tri(), K1D)\nmd = MeshData(VXY, EToV, rd)\n(; x, y ) = md\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.MeshData-Tuple{RefElemData, Any, Any}","page":"Reference","title":"StartUpDG.MeshData","text":"function MeshData(rd, geometry, vxyz...)\n\nCreates 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:\n\njulia> geometry = (PresetGeometries.Circle(R=0.33, x0=0, y0=0), )\n\nHere, 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. \n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.MeshData-Tuple{Tuple{Tuple, Matrix{Int64}}, RefElemData, Vararg{Any}}","page":"Reference","title":"StartUpDG.MeshData","text":"MeshData(VXYZ, EToV, rd::RefElemData)\nMeshData((VXYZ, EToV), rd::RefElemData)\n\nReturns a MeshData struct with high order DG mesh information from the unstructured mesh information (VXYZ..., EToV).\n\nMeshData(rd::RefElemData, md::MeshData, xyz...)\n\nGiven 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.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.MeshImportOptions","page":"Reference","title":"StartUpDG.MeshImportOptions","text":"MeshImportOptions\n\nThis struct allows the user to opt for supported features when importing a Gmsh 4.1 .msh file.\n\nSupport\n\ngrouping::Bool | On import would you like to include physical group assignements of 2D elements?\nremap_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.\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.MeshPlotter","page":"Reference","title":"StartUpDG.MeshPlotter","text":"MeshPlotter(rd::RefElemData, md::RefElemData)\n\nPlot recipe to plot a (possibly curved) quadrilateral or triangular mesh. Usage: plot(MeshPlotter(...))\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.MultipleRefElemData","page":"Reference","title":"StartUpDG.MultipleRefElemData","text":"struct MultipleRefElemData{T <: NamedTuple}\n data::T\nend\n\nHolds multiple RefElemData objects in data where typeof(data) <: NamedTuple.\n\nIndividual RefElemData can be accessed via getproperty, for example rds.Tri. \n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.NamedArrayPartition","page":"Reference","title":"StartUpDG.NamedArrayPartition","text":"NamedArrayPartition(; kwargs...)\nNamedArrayPartition(x::NamedTuple)\n\nSimilar 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. \n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.NonConformingMesh","page":"Reference","title":"StartUpDG.NonConformingMesh","text":"warning: Experimental implementation\nThis is an experimental feature and may change without warning in future releases.\n\nThis is a proof of concept implementation of a non-conforming mesh in StartUpDG.jl. The intended usage is as follows:\n\nrd = RefElemData(Quad(), N=7)\nmd = MeshData(NonConformingQuadMeshExample(), rd)\n\n(; x, y) = md\nu = @. sin(pi * x) * sin(pi * y)\n\n# interpolate to faces\nnum_total_faces = num_faces(rd.element_type) * md.num_elements\nuf = reshape(rd.Vf * u, :, num_total_faces)\n\n# interpolate faces to mortars (`uf` denotes mortar faces for `NonConformingMesh` types)\n(; nonconforming_faces, mortar_interpolation_matrix) = md.mesh_type\n\nu_mortar = reshape(mortar_interpolation_matrix * uf[:, nonconforming_faces], :, 2 * length(nonconforming_faces))\n\n# construct interior (uM = u⁻ \"minus\") values and exterior (uP = u⁺ \"plus\") values\nuM = hcat(uf, u_mortar) # uM = both element faces and mortar faces\nuP = uM[md.mapP]\n\nThe 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.\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.PhysicalFrame","page":"Reference","title":"StartUpDG.PhysicalFrame","text":"`PhysicalFrame{NDIMS} <: AbstractElemShape{NDIMS}`\n\nPhysicalFrame 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. \n\nPhysicalFrame()\nPhysicalFrame(x, y)\nPhysicalFrame(x, y, vx, vy): stores coordinates `vx, vy` of background Cartesian cell\n\nConstructors for a PhysicalFrame object (optionally uses arrays of points x, y on a cut element).\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.Polynomial","page":"Reference","title":"StartUpDG.Polynomial","text":"Polynomial{T}\n\nRepresents 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). \n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.RefElemData","page":"Reference","title":"StartUpDG.RefElemData","text":"struct RefElemData\n\nRefElemData: contains info (interpolation points, volume/face quadrature, operators) for a high order nodal basis on a given reference element. \n\nExample:\n\nN = 3\nrd = RefElemData(Tri(), N)\n(; r, s ) = rd\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"function RefElemData(elem; N, kwargs...)\nfunction RefElemData(elem, approx_type; N, kwargs...)\n\nKeyword argument constructor for RefElemData (to \"label\" N via rd = RefElemData(Line(), N=3))\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Hex, TensorProductQuadrature, Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"RefElemData(elem::Hex, ::TensorProductQuadrature, N; Nplot = 10)\nRefElemData(elem::Hex, approximation_type::Polynomial{<:TensorProductQuadrature}, N; Nplot = 10)\n\nConstructor for hexahedral RefElemData where the quadrature is assumed to have tensor product structure.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Line, Polynomial{StartUpDG.DefaultPolynomialType}, Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"RefElemData(elem::Line, N;\n quad_rule_vol = quad_nodes(elem, N+1))\nRefElemData(elem::Union{Tri, Quad}, N;\n quad_rule_vol = quad_nodes(elem, N),\n quad_rule_face = gauss_quad(0, 0, N))\nRefElemData(elem::Union{Hex, Tet}, N;\n quad_rule_vol = quad_nodes(elem, N),\n quad_rule_face = quad_nodes(Quad(), N))\nRefElemData(elem; N, kwargs...) # version with keyword args\n\nConstructor for RefElemData for different element types.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Line, SBP{TensorProductLobatto}, Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"function RefElemData(elementType::Line, approxType::SBP, N)\nfunction RefElemData(elementType::Quad, approxType::SBP, N)\nfunction RefElemData(elementType::Hex, approxType::SBP, N)\nfunction RefElemData(elementType::Tri, approxType::SBP, N)\n\nSBP reference element data for Quad(), Hex(), and Tri() elements. \n\nFor Line(), Quad(), and Hex(), approxType is SBP{TensorProductLobatto}.\n\nFor Tri(), approxType can be SBP{Kubatko{LobattoFaceNodes}}, SBP{Kubatko{LegendreFaceNodes}}, or SBP{Hicken}. \n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Pyr, Polynomial, Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"RefElemData(elem::Pyr, approximation_type::Polynomial, N;\n quad_rule_vol=quad_nodes(elem, N),\n quad_rule_face_quad=quad_nodes(Quad(), N), \n quad_rule_face_tri=quad_nodes(Tri(), N), \n quad_rule_face=(quad_rule_face_quad, quad_rule_face_tri),\n Nplot=10)\n\nBuilds operators for pyramids.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Union{Hex, Line, Quad}, Polynomial{<:Gauss}, Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"RefElemData(elem::Union{Line, Quad, Hex}, approximation_type::Polynomial{Gauss}, N)\n\nBuilds rd::RefElemData with (N+1)-point Gauss quadrature in each dimension. \n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.RefElemData-Tuple{Wedge, Polynomial, Any}","page":"Reference","title":"StartUpDG.RefElemData","text":"RefElemData(elem::Wedge, approximation_type::Polynomial, N;\n quad_rule_vol=quad_nodes(elem, N),\n quad_rule_face_quad=quad_nodes(Quad(), N), \n quad_rule_face_tri=quad_nodes(Tri(), N), \n quad_rule_face=(quad_rule_face_quad, quad_rule_face_tri),\n Nplot=10)\n\nBuilds operators for prisms/wedges\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.TensorProductQuadrature","page":"Reference","title":"StartUpDG.TensorProductQuadrature","text":"TensorProductQuadrature{T}\n\nA type parameter to Polynomial indicating that \n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.VertexMappedMesh","page":"Reference","title":"StartUpDG.VertexMappedMesh","text":"struct VertexMappedMesh\n\nThe 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.\n\nFields\n\nelement_type :: TE <: AbstractElemShape \nVXYZ :: TV \nEToV :: TEV\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#StartUpDG.VertexMeshPlotter","page":"Reference","title":"StartUpDG.VertexMeshPlotter","text":"VertexMeshPlotter((VX, VY), EToV, fv)\nVertexMeshPlotter(triout::TriangulateIO)\n\nPlot recipe to plot a quadrilateral or triangular mesh. Usage: plot(VertexMeshPlotter(...))\n\n\n\n\n\n","category":"type"},{"location":"index_refs/#NodesAndModes.equi_nodes-Tuple{PhysicalFrame{2, Shifting, Scaling} where {Shifting<:(StaticArraysCore.SVector{2}), Scaling<:(StaticArraysCore.SVector{2})}, Any, Any}","page":"Reference","title":"NodesAndModes.equi_nodes","text":"function NodesAndModes.equi_nodes(elem::PhysicalFrame, curve, N)\n\nReturns back Np(N) equally spaced nodes on the background quadrilateral corresponding to elem, with points inside of curve removed.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.MeshData_to_vtk","page":"Reference","title":"StartUpDG.MeshData_to_vtk","text":"MeshDatatovtk(md, rd, dim, data, dataname, datatype, filename, writedata = false, equidist_nodes = true)\n\nTranslate 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\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.MeshData_to_vtk-Union{Tuple{DIM}, Tuple{MeshData, RefElemData{DIM, ElemShape} where ElemShape<:AbstractElemShape{DIM}, Any, Any, Any}, Tuple{MeshData, RefElemData{DIM, ElemShape} where ElemShape<:AbstractElemShape{DIM}, Vararg{Any, 4}}, Tuple{MeshData, RefElemData{DIM, ElemShape} where ElemShape<:AbstractElemShape{DIM}, Vararg{Any, 5}}} where DIM","page":"Reference","title":"StartUpDG.MeshData_to_vtk","text":"MeshData_to_vtk(md, rd, dim, data, dataname, datatype, filename, write_data = false, equi_dist_nodes = true)\n\nTranslate 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\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.SUD_to_vtk_order-Union{Tuple{RefElemData{DIM, ElemShape} where ElemShape<:AbstractElemShape{DIM}}, Tuple{DIM}} where DIM","page":"Reference","title":"StartUpDG.SUD_to_vtk_order","text":"SUD_to_vtk_order(rd::RefElemData, dim)\n\nCompute the permutation of the nodes between StartUpDG and VTK\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.boundary_face_centroids-Tuple{Any}","page":"Reference","title":"StartUpDG.boundary_face_centroids","text":"function boundary_face_centroids(md)\n\nReturns face centroids and boundary_face_ids on the boundaries of the domain given by md::MeshData.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.build_node_maps-Tuple{Any, Any}","page":"Reference","title":"StartUpDG.build_node_maps","text":"build_node_maps(FToF, Xf)\n\nIntialize the connectivity table along all edges and boundary node tables of all elements. mapM - map minus (interior). mapP - map plus (exterior).\n\nXf = (xf, yf, zf) and FToF is size (Nfaces * K) and FToF[face] = face neighbor\n\nmapM, mapP are size Nfp x (Nfaces*K)\n\nExamples\n\njulia> mapM, mapP, mapB = build_node_maps(FToF, (xf, yf))\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.ck45-Tuple{}","page":"Reference","title":"StartUpDG.ck45","text":"ck45()\n\nReturns 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.,\n\nExample\n\nres = rk4a[i]*res + dt*rhs # i = RK stage\n@. u += rk4b[i]*res\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.connect_mesh-Tuple{Any, Any, Any}","page":"Reference","title":"StartUpDG.connect_mesh","text":"connect_mesh(rd, face_centroids, region_flags, cutcells; tol = 1e2 * eps())\n\nConnects faces of a cut mesh to each other, returns FToF such that face f is connected to FToF[f]. \n\nInputs:\n\nrd::RefElemData\nfacecentroids = (facecentroidsx, facecentroidsy), where `facecentroids_x/y` are vectors of coordinates of face centroids\nregion_flags, cutcells are return arguments from PathIntersections.define_regions\n\nThe keyword argument tol is the tolerance for matches between face centroids. \n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.connect_mesh-Tuple{Any, Any}","page":"Reference","title":"StartUpDG.connect_mesh","text":"connect_mesh(EToV,fv)\n\nInputs:\n\nEToV is a num_elements by Nv matrix whose rows identify the Nv vertices\n\nwhich make up one of the num_elements elements.\n\nfv (an array of arrays containing unordered indices of face vertices).\n\nOutput: FToF, an length(fv) by num_elements index array containing face-to-face connectivity.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.estimate_h-Union{Tuple{DIM}, Tuple{RefElemData{DIM, ElemShape} where ElemShape<:AbstractElemShape{DIM}, MeshData{DIM}}} where DIM","page":"Reference","title":"StartUpDG.estimate_h","text":"estimate_h(rd::RefElemData, md::MeshData)\nestimate_h(e, rd::RefElemData, md::MeshData) # e = element index\n\nEstimates the mesh size via min sizeofdomain * |J|/|Jf|, since |J| = O(hᵈ) and |Jf| = O(hᵈ⁻¹). \n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.findline-Tuple{String, Vector{String}}","page":"Reference","title":"StartUpDG.findline","text":"findline(word::String, lines)\n\nOutputs the line number of word in lines. \n\nIt is assumed that the word exists at least once in the file.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.get_boundary_face_labels-Tuple{Triangulate.TriangulateIO, RefElemData{2, Tri}, MeshData{2}}","page":"Reference","title":"StartUpDG.get_boundary_face_labels","text":"function get_boundary_face_labels(triout::TriangulateIO, md::MeshData{2})\n\nFind Triangle segment labels of boundary faces. Returns two arguments:\n\nboundary_face_tags: tags of faces on the boundary\nboundary_faces: list of faces on the boundary of the domain\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.get_node_boundary_tags-Tuple{Triangulate.TriangulateIO, RefElemData{2, Tri}, MeshData{2}}","page":"Reference","title":"StartUpDG.get_node_boundary_tags","text":"function get_node_boundary_tags(triout::TriangulateIO,md::MeshData{2},rd::RefElemData{2,Tri})\n\nComputes node_tags = Nfp x Nfaces * num_elements array where each entry is a Triangulate.jl tag number.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.get_num_elements","page":"Reference","title":"StartUpDG.get_num_elements","text":"returns the number of elements in a .msh file of a specified dimension\n\nNotes: Gmsh includes elements in a .msh file of multiple dimensions. We want a count of how many\n\n2D elements are in our file. This corisponds to the number of elements in our tri mesh.\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.hex_vtk_order","page":"Reference","title":"StartUpDG.hex_vtk_order","text":"hex_vtk_order(corner_verts, order, dim, skip = false)\n\nCompute 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\n\nInspired by: https://github.com/ju-kreber/paraview-scripts/blob/master/node_ordering.py\n\nVTK node numbering of a hexagon:\n\n 8+------+7\n /| /|\n / | / |\n 5+------+6 |\n\nz | 4+–-|–+3 | y | / | / |/ |/ |/ 0–> x 1+–––+2\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.hybridized_SBP_operators-Tuple{Any}","page":"Reference","title":"StartUpDG.hybridized_SBP_operators","text":"function hybridized_SBP_operators(rd::RefElemData{DIMS})\n\nConstructs hybridized SBP operators given a RefElemData. Returns operators Qrsth..., VhP, Ph.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.inverse_trace_constant-Tuple{RefElemData{1, Line, <:Polynomial}}","page":"Reference","title":"StartUpDG.inverse_trace_constant","text":"function inverse_trace_constant(rd::RefElemData)\n\nReturns 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).\n\nCan be used to estimate dependence of maximum stable timestep on degree of approximation.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.make_periodic-Union{Tuple{MeshData{Dim}}, Tuple{Dim}, Tuple{MeshData{Dim}, Bool}} where Dim","page":"Reference","title":"StartUpDG.make_periodic","text":"make_periodic(md::MeshData{Dim}, is_periodic...) where {Dim}\nmake_periodic(md::MeshData{Dim}, is_periodic = ntuple(x->true,Dim)) where {Dim}\nmake_periodic(md::MeshData, is_periodic = true)\n\nReturns 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\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.n_verts_between-Tuple{Any, Any, Any}","page":"Reference","title":"StartUpDG.n_verts_between","text":"n_verts_between(n, from, to, dim)\n\nCompute 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\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.quad_vtk_order","page":"Reference","title":"StartUpDG.quad_vtk_order","text":"quad_vtk_order(corner_verts, order, dim, skip = false)\n\nCompute 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\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.read_Gmsh_2D-Tuple{String, Vararg{Any}}","page":"Reference","title":"StartUpDG.read_Gmsh_2D","text":"read_Gmsh_2D(filename, args...)\n\nReads a 2D triangular Gmsh file. Mesh formats 2.2 and 4.1 supported. Returns (VX, VY), EToV. \n\nExamples\n\nVXY, EToV = read_Gmsh_2D(\"eulerSquareCylinder2D.msh\") # v2.2 file format\nVXY, EToV = read_Gmsh_2D(\"test/testset_Gmsh_meshes/periodicity_mesh_v4.msh\") # v4.1 file format\n\n# if MeshImportOptions.grouping=true, then a third variable `grouping` is returned\nVXY, EToV, grouping = read_Gmsh_2D(\"test/testset_Gmsh_meshes/periodicity_mesh_v4.msh\", MeshImportOptions(true, false))\nVXY, EToV, grouping = read_Gmsh_2D(\"test/testset_Gmsh_meshes/periodicity_mesh_v4.msh\", true) # same as above\n\nSee also\n\nhttps://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format\nhttps://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format-version-2-0028Legacy0029\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.read_Gmsh_2D_v2-Tuple{String}","page":"Reference","title":"StartUpDG.read_Gmsh_2D_v2","text":"read_Gmsh_2D_v2(filename)\n\nReads triangular GMSH 2D file format 2.2 0 8. Returns (VX, VY), EToV.\n\nExamples\n\nVXY, EToV = read_Gmsh_2D_v2(\"eulerSquareCylinder2D.msh\")\n\nhttps://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format-version-2-0028Legacy0029\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.read_Gmsh_2D_v4","page":"Reference","title":"StartUpDG.read_Gmsh_2D_v4","text":"For brevity when grouping is the only supported feature.\n\nexample: VXY, EToV, grouping = read_Gmsh_2D_v4(\"file.msh\",true)\nexample: VXY, EToV = read_Gmsh_2D_v4(\"file.msh\",false)\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.read_Gmsh_2D_v4-Tuple{String, MeshImportOptions}","page":"Reference","title":"StartUpDG.read_Gmsh_2D_v4","text":"function read_Gmsh_2D_v4(filename, options)\n\nreads triangular GMSH 2D .msh files.\n\nOutput\n\nThis depends on if grouping is opted for or not\n\nreturns: (VX, VY), EToV\nreturns: (VX, VY), EToV, grouping\n\nSupported formats and features:\n\nversion 4.1 'physical group support 'remap group ids\n\ngrouping application\n\nWhen 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.\n\nVXY, EToV = read_Gmsh_2D_v4(\"eulerSquareCylinder2D.msh\")\nVXY, EToV = read_Gmsh_2D_v4(\"eulerSquareCylinder2D.msh\",false)\nVXY, EToV, grouping = read_Gmsh_2D_v4(\"eulerSquareCylinder2D.msh\", true)\n\noption = MeshImportOption(true)\nVXY, EToV, grouping = read_Gmsh_2D_v4(\"eulerSquareCylinder2D.msh\", option)\n\nhttps://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format\n\nNotes: the version 4 format has a more detailed block data format this leads to more complicated parser.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.refine","page":"Reference","title":"StartUpDG.refine","text":"function refine(triout, h, href = h/2)\n\nRefinement of a previous mesh given the current mesh size h. Preserves boundary/volume tags.\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.remap_element_grouping-Tuple{Vector{Int64}}","page":"Reference","title":"StartUpDG.remap_element_grouping","text":"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\n\nExample output\n\nremapelementgrouping([16,16,17,17]) -> [1,1,2,2]\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.sort_by_axis-Tuple{Any}","page":"Reference","title":"StartUpDG.sort_by_axis","text":"sort_by_axis(corner_verts)\n\nGiven the points 'corner_verts' sort them in a lexicographical order and return the permutated points. \n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.sparse_low_order_SBP_operators-Union{Tuple{RefElemData{NDIMS, ElemShape} where ElemShape<:AbstractElemShape{NDIMS}}, Tuple{NDIMS}} where NDIMS","page":"Reference","title":"StartUpDG.sparse_low_order_SBP_operators","text":"function sparse_low_order_SBP_operators(rd)\n\nConstructs sparse low order SBP operators given a RefElemData. Returns operators Qrst..., E ≈ Vf * Pq that satisfy a generalized summation-by-parts (GSBP) property:\n\n `Q_i + Q_i^T = E' * B_i * E`\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.tag_boundary_faces-Tuple{Any, Nothing}","page":"Reference","title":"StartUpDG.tag_boundary_faces","text":"function tag_boundary_faces(md, boundary_name::Symbol = :entire_boundary)\nfunction tag_boundary_faces(md, boundary_list::Dict{Symbol, <:Function})\n\nWhen called without arguments, just returns Dict(:entire_boundary => boundary_faces)`.\n\nExample usage: \n\njulia> rd = RefElemData(Tri(), N=1)\njulia> md = MeshData(uniform_mesh(Tri(), 2)..., rd)\njulia> on_bottom_boundary(x, y, tol = 1e-13) = abs(y+1) < tol\njulia> on_top_boundary(x, y, tol = 1e-13) = abs(y-1) < tol\njulia> tag_boundary_faces(Dict(:bottom => on_bottom_boundary,\n :top => on_top_boundary), md)\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.tag_boundary_faces-Tuple{Triangulate.TriangulateIO, RefElemData{2, Tri}, MeshData{2}, NamedTuple}","page":"Reference","title":"StartUpDG.tag_boundary_faces","text":"function tag_boundary_faces(triout::TriangulateIO,\n rd::RefElemData{2,Tri}, md::MeshData{2},\n boundary_list::Union{NamedTuple,Dict{Symbol,Int}})\n\nHere, 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.\n\nExample usage:\n\njulia> using Triangulate, StartUpDG\njulia> triout = scramjet()\njulia> rd = RefElemData(Tri(),N=1)\njulia> md = MeshData(triangulateIO_to_VXYEToV(triout)...,rd)\njulia> tag_boundary_faces(triout,rd,md, Dict(:wall=>1, :inflow=>2, :outflow=>3))\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.triangle_vtk_order","page":"Reference","title":"StartUpDG.triangle_vtk_order","text":"triangle_vtk_order(corner_verts, order, dim, skip = false)\n\nCompute 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\n\n\n\n\n\n","category":"function"},{"location":"index_refs/#StartUpDG.triangulateIO_to_VXYEToV-Tuple{Triangulate.TriangulateIO}","page":"Reference","title":"StartUpDG.triangulateIO_to_VXYEToV","text":"function triangulateIO_to_VXYEToV(triout::TriangulateIO)\n\nComputes VX,VY,EToV from a TriangulateIO object.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.type_to_vtk-Tuple{Hex}","page":"Reference","title":"StartUpDG.type_to_vtk","text":"type_to_vtk(elem::Hex)\nreturn the VTK-type\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.type_to_vtk-Tuple{Quad}","page":"Reference","title":"StartUpDG.type_to_vtk","text":"type_to_vtk(elem::Quad)\nreturn the VTK-type\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.type_to_vtk-Tuple{Tri}","page":"Reference","title":"StartUpDG.type_to_vtk","text":"type_to_vtk(elem::Tri)\nreturn the VTK-type\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.type_to_vtk-Tuple{Wedge}","page":"Reference","title":"StartUpDG.type_to_vtk","text":"type_to_vtk(elem::Wedge)\nreturn the VTK-type\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.uniform_mesh-Tuple{Line, Any}","page":"Reference","title":"StartUpDG.uniform_mesh","text":"uniform_mesh(elem::Line,Kx)\nuniform_mesh(elem::Tri,Kx,Ky)\nuniform_mesh(elem::Quad,Kx,Ky)\nuniform_mesh(elem::Hex,Kx,Ky,Kz)\nuniform_mesh(elem, K)\n\nUniform Kx (by Ky by Kz) mesh on -11^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.\n\nK can also be specified using a keyword argument K1D, e.g., uniform_mesh(elem; K1D = 16).\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.vtk_order-Tuple{Hex, Any}","page":"Reference","title":"StartUpDG.vtk_order","text":"vtk_order(elem::Hex, order)\n\nConstruct 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.\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.vtk_order-Tuple{Quad, Any}","page":"Reference","title":"StartUpDG.vtk_order","text":"vtk_order(elem::Quad, order)\n\nConstruct 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\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.vtk_order-Tuple{Tri, Any}","page":"Reference","title":"StartUpDG.vtk_order","text":"vtk_order(elem::Tri, order)\n\nConstruct 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\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.vtk_order-Tuple{Wedge, Any}","page":"Reference","title":"StartUpDG.vtk_order","text":"vtk_order(elem::Wedge, order)\n\nConstruct all node-points of a VTKLAGRANGEWEDGE of order order. The corner-nodes are given by the reference-wedge used by StartUpDG\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#StartUpDG.wedge_vtk_order-Tuple{Any, Any, Any}","page":"Reference","title":"StartUpDG.wedge_vtk_order","text":"wedge_vtk_order(corner_verts, order, dim)\n\nCompute 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\n\n\n\n\n\n","category":"method"},{"location":"index_refs/#Triangulate.triangulate","page":"Reference","title":"Triangulate.triangulate","text":"function Triangulate.triangulate(triin::TriangulateIO, maxarea, minangle=20)\n\nConvenience routine to avoid writing out @sprintf each time. Returns a TriangulateIO object.\n\n\n\n\n\n","category":"function"},{"location":"tstep_usage/#Time-stepping","page":"Timestepping","title":"Time-stepping","text":"","category":"section"},{"location":"tstep_usage/","page":"Timestepping","title":"Timestepping","text":"For convenience, we include a commonly used explicit Runge-Kutta scheme. For more advanced time-stepping functionality, we recommend OrdinaryDiffEq.jl. ","category":"page"},{"location":"tstep_usage/#Carpenter-and-Kennedy's-(4,5)-method","page":"Timestepping","title":"Carpenter and Kennedy's (4,5) method","text":"","category":"section"},{"location":"tstep_usage/","page":"Timestepping","title":"Timestepping","text":"ck45() returns coefficients for a low-storage Runge-Kutta method.","category":"page"},{"location":"tstep_usage/#Example-usage:","page":"Timestepping","title":"Example usage:","text":"","category":"section"},{"location":"tstep_usage/","page":"Timestepping","title":"Timestepping","text":"using Plots\nusing StartUpDG\n\n# Brusselator\nB = 3\nf(y1, y2) = [1 + y1^2 * y2 - (B+1) * y1, B * y1 - y1^2 * y2]\nf(Q) = f(Q[1], Q[2])\np,q = 1.5, 3.0\nQ = [p; q]\n\ndt = .1\nFinalTime = 20\n\nres = zero.(Q) # init RK residual\nrk4a, rk4b, rk4c = ck45()\nNsteps = ceil(FinalTime / dt)\ndt = FinalTime / Nsteps\n\nplot() # init plot\nfor i = 1:Nsteps\n global res # yes, I know...this is just for simplicty\n for INTRK = 1:5\n rhsQ = f(Q)\n @. res = rk4a[INTRK] * res + dt * rhsQ # i = RK stage\n @. Q = Q + rk4b[INTRK] * res\n end\n scatter!([i*dt;i*dt],Q)\nend\ndisplay(plot!(leg=false))","category":"page"},{"location":"more_meshes/#Additional-mesh-types","page":"Additional mesh types","title":"Additional mesh types","text":"","category":"section"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"In addition to more \"standard\" mesh types, StartUpDG.jl also has experimental support for hybrid and cut-cell meshes. Both are currently restricted to two dimensional domains.","category":"page"},{"location":"more_meshes/#Hybrid-meshes","page":"Additional mesh types","title":"Hybrid meshes","text":"","category":"section"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"warning: Experimental implementation\nThis is an experimental feature and may change in future releases.","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"There is initial support for hybrid meshes in StartUpDG.jl. The following is a short example where we interpolate a polynomial and compute its derivative.","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"rds = RefElemData((Tri(), Quad()), N = 3)\n\n# Simple hybrid mesh for testing\n# 1 7______8______9\n# | | 3 / |\n# | 4 | / 5 |\n# 0 4 ---- 5 ---- 6 \n# | | |\n# | 1 | 2 |\n# -1 1 ---- 2 ---- 3\n# -1 0 1\nVX = [-1; 0; 1; -1; 0; 1; -1; 0; 1]\nVY = [-1; -1; -1; 0; 0; 0; 1; 1; 1]\nEToV = [[1 2 4 5], [2 3 5 6], [5 8 9], [4 5 7 8], [9 6 5]]\n\nmd = MeshData(VX, VY, EToV, rds)\n\n# test that the local derivatives of a polynomial recover the exact derivative\n(; x, y ) = md\nu = @. x^3 - x^2 * y + 2 * y^3\ndudx = @. 3 * x^2 - 2 * x * y\n\n# compute local derivatives\n(; rxJ, sxJ, J ) = md\ndudr, duds = similar(md.x), similar(md.x)\ndudr.Quad .= rds[Quad()].Dr * u.Quad\nduds.Quad .= rds[Quad()].Ds * u.Quad\ndudr.Tri .= rds[Tri()].Dr * u.Tri\nduds.Tri .= rds[Tri()].Ds * u.Tri\n\n@show norm(@. dudx - (rxJ * dudr + sxJ * duds) / J) # should be O(1e-14)","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"The main difference in the representation of hybrid meshes compared with standard MeshData objects is the use of NamedArrayPartition arrays as storage for the geometric coordinates. These arrays have \"fields\" corresponding to the element type, for example","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"md.x.Tri\nmd.x.Quad","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"but can still be indexed as linear arrays. ","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"The mapP field behaves similarly. If we interpolate the values of u for each element type to surface quadrature nodes, we can use mapP to linearly index into the array to find neighbors. ","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"uf = similar(md.xf)\nuf.Quad .= rds.Quad.Vf * u.Quad\nuf.Tri .= rds.Tri.Vf * u.Tri\nuf[md.mapP] # this returns the exterior node values","category":"page"},{"location":"more_meshes/#Cut-Meshes","page":"Additional mesh types","title":"Cut Meshes","text":"","category":"section"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"warning: Experimental implementation\nThis is an experimental feature and may change in future releases.","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"Initial support for cut-cell meshes is available via PathIntersections.jl. By passing in a tuple of curves (defined as parametrized functions of one coordinate, see PathIntersections.jl documentation for more detail), StartUpDG.jl can compute a MeshData for a cut-cell mesh. ","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"circle = PresetGeometries.Circle(R=0.33, x0=0, y0=0)\ncells_per_dimension_x, cells_per_dimension_y = 4, 4\n\nrd = RefElemData(Quad(), N=3)\nmd = MeshData(rd, (circle, ), cells_per_dimension_x, cells_per_dimension_y; precompute_operators=true)","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"The interpolation points on cut cells md.x.cut are determined from sampled points and a pivoted QR decomposition. The quadrature points on cut cells md.xq.cut are determined similarly. However, the cut-cell quadrature weights md.wJq.cut are not currently positive. The optional keyword argument precompute_operators specifies whether to precompute differentiation, face interpolation, mass, and lifting matrices for each cut cell. If precompute_operators=true, these are stored in md.mesh_type.cut_cell_operators. ","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"As with hybrid meshes, the nodal coordinates md.x, md.y are NamedArrayPartitions with cartesian and cut fields. For example, md.x.cartesian and md.x.cut are the x-coordinates of the Cartesian and cut cells, respectively. Likewise, md.mapP indexes linearly into the array of face coordinates and specifies exterior node indices. For example, we can interpolate a function to face nodes and compute exterior values via the following code:","category":"page"},{"location":"more_meshes/","page":"Additional mesh types","title":"Additional mesh types","text":"(; x, y) = md\nu = @. x^rd.N - x * y^(rd.N-1) - x^(rd.N-1) * y + y^rd.N # some random function \n\n# interpolate the solution to face nodes\nuf = similar(md.xf) \nuf.cartesian .= rd.Vf * u.cartesian\nfor e in 1:size(md.x.cut, 2)\n ids = md.mesh_type.cut_face_nodes[e]\n Vf = md.mesh_type.cut_cell_operators.face_interpolation_matrices[e]\n uf.cut[ids] .= Vf * u.cut[:, e]\nend\n\nuf[md.mapP] # these are \"exterior\" values for each entry of `uf`","category":"page"},{"location":"ex_dg_deriv/#Example:-approximating-derivatives-using-DG","page":"Example: computing DG derivatives","title":"Example: approximating derivatives using DG","text":"","category":"section"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"MeshData can be used to compute DG derivatives. Suppose f is a differentiable function and the domain Omega can be decomposed into non-overlapping elements D^k. The approximation of fracpartial fpartial x can be approximated using the following formulation: find piecewise polynomial u such that for all piecewise polynomials v","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"int_Omega u v = sum_k left( int_D^k fracpartial upartial xv + int_partial D^k frac12 lefturightn_x v right)","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"Here, lefturight = u^+ - u denotes the jump across an element interface, and n_x is the x-component of the outward unit normal on D^k.","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"Discretizing the left-hand side of this formulation yields a mass matrix. Inverting this mass matrix to the right hand side yields the DG derivative. We show how to compute it for a uniform triangular mesh using MeshData and StartUpDG.jl.","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"We first construct the triangular mesh and initialize md::MeshData.","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"using StartUpDG\nusing Plots\n\nN = 3\nK1D = 8\nrd = RefElemData(Tri(), N)\nVXY, EToV = uniform_mesh(Tri(), K1D)\nmd = MeshData(VXY, EToV, rd)","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"We can approximate a function f(x y) using interpolation","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"f(x, y) = exp(-5 * (x^2 + y^2)) * sin(1 + pi*x) * sin(2 + pi*y)\n(; x, y) = md\nu = @. f(x, y)","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"or using quadrature-based projection","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"(; Pq ) = rd\n(; x, y, xq, yq ) = md\nu = Pq * f.(xq, yq)","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"We can use scatter in Plots.jl to quickly visualize the approximation. This is not intended to create a high quality image (see other libraries, e.g., Makie.jl,VTK.jl, or Triplot.jl for publication-quality images).","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"(; Vp ) = rd\nxp, yp, up = Vp * x, Vp * y, Vp * u # interp to plotting points\nscatter(xp, yp, uxp, zcolor=uxp, msw=0, leg=false, ratio=1, cam=(0, 90))","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"Both interpolation and projection create a matrix u of size N_p times K which contains coefficients (nodal values) of the DG polynomial approximation to f(x y). We can approximate the derivative of f(x y) using the DG derivative formulation","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"function dg_deriv_x(u, rd::RefElemData, md::MeshData)\n (; Vf, Dr, Ds, LIFT ) = rd\n (; rxJ, sxJ, J, nxJ, mapP ) = md\n uf = Vf * u\n ujump = uf[mapP] - uf\n\n # derivatives using chain rule + lifted flux terms\n ux = rxJ .* (Dr * u) + sxJ .* (Ds * u) \n dudxJ = ux + LIFT * (.5 * ujump .* nxJ)\n\n return dudxJ ./ J\nend","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"We can visualize the result as follows:","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"dudx = dg_deriv_x(u, rd, md)\nuxp = Vp * dudx\nscatter(xp, yp, uxp, zcolor=uxp, msw=0, leg=false, ratio=1, cam=(0,90))","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"Plots of the polynomial approximation u(xy) and the DG approximation of fracpartial upartial x are given below","category":"page"},{"location":"ex_dg_deriv/","page":"Example: computing DG derivatives","title":"Example: computing DG derivatives","text":"(Image: u) (Image: dudx) ⠀","category":"page"},{"location":"#Overview","page":"Home","title":"Overview","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package contains routines to initialize reference element operators, physical mesh arrays, and connectivity arrays for nodal DG methods. The codes roughly based on Nodal Discontinuous Galerkin Methods by Hesthaven and Warburton (2007).","category":"page"},{"location":"#A-short-example","page":"Home","title":"A short example","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"using StartUpDG\n\n# polynomial degree and mesh size\nN = 3\nK1D = 8\n\n# init ref element and mesh\nrd = RefElemData(Tri(), N)\nVXY, EToV = uniform_mesh(Tri(), K1D)\nmd = MeshData(VXY, EToV, rd)\n\n# Define a function by interpolation\n(; x, y ) = md\nu = @. exp(-10 * (x^2 + y^2))\n\n# Compute derivatives using geometric mapping + chain rule\n(; Dr, Ds ) = rd\n(; rxJ, sxJ, J ) = md\ndudx = (rxJ .* (Dr * u) + sxJ .* (Ds * u)) ./ J","category":"page"}] } diff --git a/dev/tstep_usage/index.html b/dev/tstep_usage/index.html index 4a3fc1d2..86a27f1f 100644 --- a/dev/tstep_usage/index.html +++ b/dev/tstep_usage/index.html @@ -27,4 +27,4 @@ end scatter!([i*dt;i*dt],Q) end -display(plot!(leg=false)) +display(plot!(leg=false))