diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 483370cd..007006bc 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-03T13:26:23","documenter_version":"1.5.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-03T13:27:58","documenter_version":"1.5.0"}} \ No newline at end of file diff --git a/dev/MeshData/index.html b/dev/MeshData/index.html index d1aeda0a..6097ee4c 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 8219d2e3..4c9109da 100644 --- a/dev/RefElemData/index.html +++ b/dev/RefElemData/index.html @@ -43,4 +43,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 2b834f23..408ba1ce 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 795984e4..4b70baee 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 8a06cbb5..95328ab3 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 e7556bf1..c244ad6e 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
MeshData(rd::RefElemData, objects, 
+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
MeshData(rd::RefElemData, objects, 
          vx::AbstractVector, vy::AbstractVector,
          quadrature_type=Subtriangulation(); 
          quad_rule_face=get_1d_quadrature(rd), 
@@ -102,4 +102,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 3ed1beee..ded5cdf7 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/objects.inv b/dev/objects.inv index 6dd3977c..f0c8b44a 100644 Binary files a/dev/objects.inv and b/dev/objects.inv differ diff --git a/dev/tstep_usage/index.html b/dev/tstep_usage/index.html index a3b3c56e..6cd24763 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))