Weak form modal DG methods on simplicial (and quad/hex meshes) (#647)

* first draft of weak form modal DG methods * fix naming of new folder * adding structarrays and Trixi interface stuff * adding calc_source! * adding sandbox folder for experiments todo: remove later * new VertexMappedMesh type introduces a wrapper around MeshData and boundary tags * adding timers * removing "has_nonconservative_terms" for now * updating scratch elixir * adding draft of BoundaryStateDirichlet/Wall * updating elixir * updating elixir for testing * adding specialization for nodal SBP methods * minor minor formatting and docstring edits, removing an unnecessary type annotation, renaming variables * formatting + fixed BC bug * working elixir with weak form DG on triangles * calc_error_norms returns component-wise errors * removing unused boundary condition code * adding `func` to calc_error_norms * making AbstractMeshData !<: AbstractMesh * minor formatting * adding periodic/non-periodic elixirs * remove unused type alias / annotations * adding integrate * adding analysis callback code adding analysis to Trixi.jl, moving analysis out of dg.jl * removing extra ::AbstractVector annotation * adding tests * adding `2d_tri` tests to CI * force add Manifest.toml * Revert "force add Manifest.toml" This reverts commit 74a9a46. * adding ape triangular example * modified @requires statements make it so that Trixi @requires both StructArrays/StartUpDG for triangular mesh * switch to L2 projection for compute_coefficients * removing cruft "sandbox" testing folder * formatting (2-space tabs, alignment) * renaming folder dg_tri -> dg_triangles * Update src/solvers/dg_triangles/mesh.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update .github/workflows/ci.yml Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update examples/triangular_mesh_2D/elixir_ape_triangular_mesh.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update examples/triangular_mesh_2D/elixir_ape_triangular_mesh.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update src/solvers/dg_triangles/mesh.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update examples/triangular_mesh_2D/elixir_euler_triangular_mesh.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update examples/triangular_mesh_2D/elixir_euler_periodic_triangular_mesh.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update src/solvers/dg_triangles/dg.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update test/test_examples_2d_triangular_mesh.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update test/runtests.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update src/solvers/dg_triangles/dg.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update src/solvers/dg_triangles/dg.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update src/solvers/dg_triangles/dg.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update src/solvers/dg_triangles/dg.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update src/solvers/dg_triangles/mesh.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update src/solvers/dg_triangles/mesh.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * adding new elixirs + tagging as experimental * using `each___` for iteration * minor doc changes * update NEWS.md * compat bounds * fix eltype in analysis * update path dg_tri -> dg_triangles * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * caught mistake in L2 norm computation * fixed elixir name: euler -> ape * adding convergence elixir * adding compat with StartUpDG@0.9.1 new SBP type selector API. minor version bump was to avoid conflict in exporting DGSEM... * remove SBP() from convergence + bump test project.toml versions * renaming folder dg_triangles - > dg_simplices * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * fixing bug in digest_boundary_conditions for triangles * specializing source_terms in calc_sources! * change dg_triangles -> dg_simplices in Trixi.jl * adding Triangulate VertexMappedMesh constructor * adding example using Triangulate * specializing digest_boundary_conditions to avoid ambiguity * adding docs for simplicial meshes * update gitignore for MacOS ignore .DS_Store * updating docs [skip ci] * Apply suggestions from code review [skip ci] Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * updating docs - capping lines at 100 columns - replacing N -> polydeg - replacing md.sJ -> md.Jf - adding more descriptions of quadrature rules and SBP nodes - finishing descriptions of plotting-related fields in RefElemData * replacing N -> polydeg in elixirs * replacing md.sJ -> md.Jf in solver code * bump StartUpDG compat * specializing existing `digest_boundary_conditions` routines - restricting x_neg, x_pos, ... unpacking routines to mesh::TreeMesh * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * removing elixir which duplicates coverage updating elixirs for CI tests set up parameters to reduce runtimes * updating # todos - # todo -> # Todo: simplices, removing old todos. - fixing todo for calc_error_norms updating todos in dg.jl * updating `each___` iteration function names * adding @threaded * adding Base.show for VertexMappedMesh * adding Octavian - adding Octavian to Project.toml - adding `using Octavian` to `Trixi.jl` * updating l2, linf errors for test for triangular elixirs * updating test/Project.toml - temporarily removing Plots.jl to test CI on the PR * trying Octavian: matmul! * Update test/test_examples_2d_triangular_mesh.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * fixing linf::MVector issue - `similar(::SVector)` returns an `MVector`, so I replaced it * selectively using Octavian `matmul!` revert to `mul!` for 5-argument multiply with vector arguments. Octavian hangs (see JuliaLinearAlgebra/Octavian.jl#103) * compat bound for Octavian, reverting OrdinaryDiffEq compat * fixing bugs for SBP-type discretizations - specialize `mul_by!` for UniformScaling - using LinearAlgebra: UniformScaling - dispatch on subtype of SBP * fixing Octavian compat had put Octavian compat in test/Project.toml instead of Project.toml * fixing thread-related bug - adding threading and preallocated thread arrays other minor change: removed unnecessary StructArray conversion in `compute_coefficients` * update `todo` label for Octavian.matmul! * update NEWS.md with simplicial mesh updates * renaming elixirs and test: triangular_mesh -> simplicial_mesh * adding support for 3D tet meshes - made most simplicial routines dimension-agnostic * added `min_max_speed_naive` for CompressibleEulerEquations3D * updating NEWS.md with tet meshes * fixing test naming: 2d_triangles -> simplices * updating docs with news about simplices * adding minor comments * VertexMappedMesh: fix order of type parameters * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com> Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update docs/make.jl Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com> * reorganizing simplicial examples new folders: `simplicial_2d_dg` and `simplicial_3d_dg` * resolving issue with `mul_by!` dispatch for UniformScaling - defaults to `mul!` for `A::UniformScaling` * renaming/reorganizing simplicial tests in ci * updating docs with tentative name for Triangulate.jl elixir * adding Plots.jl v1.16 to test compat * adding convergence test for triangular mesh * updating `eachdim(mesh)`, also `@inline`ing `each___` functions Co-authored-by: Jesse Chan <jesse.chan@rice.edu> Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
trixi-framework · Jul 1, 2021 · b37468f · b37468f
1 parent 5a1c448
commit b37468f
Show file tree

Hide file tree

Showing 23 changed files with 1,021 additions and 6 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -63,6 +63,8 @@ jobs:
           - 2d_part1
           - 2d_part2
           - 2d_part3
+          - 2d_simplices
+          - 3d_simplices          
           - 2d_mpi
           - 2d_threaded
           - 3d_part1

diff --git a/.gitignore b/.gitignore
@@ -21,3 +21,5 @@ coverage_report/
 **/solution_*.txt
 **/restart_*.txt
 .vscode/
+
+.DS_Store
diff --git a/NEWS.md b/NEWS.md
@@ -22,6 +22,7 @@ for human readability.
 - New unstructured, curvilinear, adaptive (non-conforming) mesh type `P4estMesh` in 2D and 3D (experimental)
 - Experimental support for finite difference (FD) summation-by-parts (SBP) methods via
   [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl)
+- New support for modal DG and SBP-DG methods on triangular and tetrahedral meshes via [StartUpDG.jl](https://github.com/jlchan/StartUpDG.jl) (experimental). 
 
 #### Changed
 

diff --git a/Project.toml b/Project.toml
@@ -13,6 +13,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
@@ -36,6 +37,7 @@ LinearMaps = "2.7, 3.0"
 LoopVectorization = "0.12.35"
 MPI = "0.16, 0.17, 0.18"
 OffsetArrays = "1.3"
+Octavian = "0.2.20"
 P4est = "0.2.1"
 Polyester = "0.3"
 RecipesBase = "1.1"

diff --git a/docs/make.jl b/docs/make.jl
@@ -54,6 +54,7 @@ makedocs(
                 "Tree mesh" => joinpath("meshes", "tree_mesh.md"),
                 "Structured mesh" => joinpath("meshes", "structured_mesh.md"),
                 "Unstructured mesh" => joinpath("meshes", "unstructured_quad_mesh.md"),
+                "Simplicial mesh" => joinpath("meshes", "mesh_data_meshes.md"),
             ],
             "Time integration" => "time_integration.md",
             "Callbacks" => "callbacks.md",

diff --git a/docs/src/meshes/mesh_data_meshes.md b/docs/src/meshes/mesh_data_meshes.md
@@ -0,0 +1,123 @@
+# Unstructured simplicial meshes
+
+Support for simplicial meshes uses the [StartUpDG.jl](https://github.com/jlchan/StartUpDG.jl) 
+package. These meshes are constructed by specifying a list of vertex coordinates `VX`, `VY`, `VZ` 
+and a connectivity matrix `EToV` where `EToV[e,:]` gives the vertices which correspond to element `e`. 
+
+We make a few simplifying assumptions about simplicial meshes:
+* meshes consist of a single type of element
+* 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 polynomial (currently, only affine 
+mappings are supported).
+
+## `AbstractMeshData` wrapper type
+
+Simplicial meshes in Trixi are represented using an `AbstractMeshData{Dim, ElemType}` wrapper. 
+For example, `VertexMappedMesh{Dim, Tri} <: AbstractMeshData{Dim, Tri}` describes a mesh whose
+reference-to-physical mapping can be constructed using only the vertex positions (e.g., an affine 
+simplicial mesh or a bilinear quadrilateral mesh). These quantities are generated by most simplicial 
+mesh generators, and `StartUpDG.jl` includes both simple uniform meshes via `uniform_mesh`, as 
+well as support for triangular meshes constructed using 
+[Triangulate.jl](https://github.com/JuliaGeometry/Triangulate.jl), a wrapper around Jonathan Shewchuk's 
+[Triangle](https://www.cs.cmu.edu/~quake/triangle.html) package.
+
+The only fields are `md::MeshData`, which contains geometric terms derived from the mapping between 
+the reference and physical elements, and `boundary_faces`, which contains a `Dict` of boundary 
+segment names (symbols) and list of faces which lie on that boundary segment. 
+
+## Variable naming conventions
+
+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. 
+
+!["Ref-to-physical mapping"](https://jlchan.github.io/StartUpDG.jl/dev/assets/mapping_diagram.png)
+
+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. 
+
+## Variable meanings and conventions in `StartUpDG.jl`
+
+`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 (but not identical) to those in 
+`Globals1D, Globals2D, Globals3D` in the Matlab codes from "Nodal Discontinuous Galerkin Methods" 
+by Hesthaven and Warburton (2007). 
+
+In general, we use the following code conventions:
+* variables such as `r, s,...` and `x, y,...` 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 "plotting" points, which are usually 
+a fine grid of equispaced points.
+* `V` matrices correspond to interpolation matrices from nodal interpolation points, e.g., `Vq` 
+interpolates to volume quadrature points, `Vf` interpolates to face quadrature points. 
+* geometric quantities in `MeshData` are stored as matrices of dimension 
+``\text{number of points per element } \times \text{number of elements}``. 
+
+Quantities in `rd::RefElemData`: 
+* `rd.Np, rd.Nq, rd.Nf`: the number of nodal interpolation points, volume quadrature points, and 
+face quadrature points on the reference element, respectively. 
+* `rd.Vq`: interpolation matrices from values at nodal interpolation points to volume quadrature points
+* `rd.wq`: volume quadrature weights on the reference element
+* `rd.Vf`: interpolation matrices from values at nodal interpolation points to face quadrature points
+* `rd.wf`: a vector containing face quadrature weights on the reference element
+* `rd.M`: the quadrature-based mass matrix, computed via `rd.Vq' * diagm(rd.wq) * rd.Vq`.
+* `rd.Pq`: a quadrature-based ``L^2`` projection matrix `rd.Pq = rd.M \ rd.Vq' * diagm(rd.wq)` 
+which maps between values at quadrature points and values at nodal points. 
+* `Dr, Ds, Dt` matrices are nodal differentiation matrices with respect to the ``r,s,t`` coordinates, 
+e.g., `Dr*f.(r,s)` approximates the derivative of ``f(r,s)`` at nodal points. 
+
+Quantities in `md::MeshData`: 
+* `md.xyz` is a tuple of matrices `md.x`, `md.y`, `md.z`, where column `e` contains coordinates of 
+physical interpolation points. 
+* `md.xyzq` is a tuple of matrices `md.xq`, `md.yq`, `md.zq`, where column `e` contains coordinates 
+of physical quadrature points. 
+* `md.rxJ, md.sxJ, ...` are matrices where column `e` contains values of 
+``J\frac{\partial r}{\partial x}``, ``J\frac{\partial s}{\partial x}``, etc. at nodal interpolation 
+points on the element `e`.
+* `md.J` is a matrix where column `e` contains values of the Jacobian ``J`` at nodal interpolation points.
+* `md.Jf` is a matrix where column `e` contains values of the face Jacobian (e.g., determinant of 
+the geometric mapping between a physical face and a reference face) at face quadrature points.
+* `md.nxJ, md.nyJ, ...` are matrices where column `e` contains values of components of the unit 
+normal scaled by the face Jacobian `md.Jf` at face quadrature points.
+
+For more details, please see the [StartUpDG.jl docs](https://jlchan.github.io/StartUpDG.jl/dev/). 
+## Special options
+
+Trixi solvers on simplicial meshes use `DG` solver types with the basis field set as a `RefElemData` 
+type. By default, `RefElemData(Tri(), polydeg)` will generate data for a degree `polydeg` polynomial 
+approximation on a triangle. There are also several parameters which can be tweaked:
+
+* `RefElemData(Tri(), polydeg, quad_rule_vol = quad_nodes(Tri(), Nq))` will substitute in a volume 
+quadrature rule of degree `Nq` instead of the default (which is a quadrature rule of degree `polydeg`).
+Here, a degree `Nq` rule will be exact for at least degree `2*Nq` integrands (such that the mass 
+matrix is integrated exactly). Quadrature rules of which exactly integrate degree `Nq` integrands 
+may also be specified by setting `quad_rule_vol = StartUpDG.quad_nodes_tri(Nq)`. 
+* `RefElemData(Tri(), polydeg, quad_rule_face = quad_nodes(Line(), Nq))` will use a face quadrature rule 
+of degree `Nq` rather than the default. This rule is also exact for at least degree `2*Nq` integrands. 
+* `RefElemData(Tri(), SBP(), polydeg)` will generate a reference element corresponding to a 
+multi-dimensional SBP discretization. Types of SBP discretizations available include 
+`SBP{Kubatko{LobattoFaceNodes}}()` (the default choice), `SBP{Kubatko{LegendreFaceNodes}}()`, and 
+[`SBP{Hicken}()`](https://doi.org/10.1007/s10915-020-01154-8). For `polydeg = 1, ..., 4`, the 
+`SBP{Kubatko{LegendreFaceNodes}}()` SBP nodes are identical to the SBP nodes of 
+[Chen and Shu](https://doi.org/10.1016/j.jcp.2017.05.025). 
+More detailed descriptions of each SBP node set can be found in the 
+[StartUpDG.jl docs](https://jlchan.github.io/StartUpDG.jl/dev/RefElemData/#RefElemData-based-on-SBP-finite-differences). 
+Trixi will also specialize certain parts of the solver based on the `SBP` approximation type. 
+
+## Trixi elixirs on simplicial meshes
+
+Example elixirs on simplicial meshes can be found in the `examples/simplicial_mesh` folder. 
+Some key elixirs to look at: 
+
+* `elixir_euler_triangular_mesh.jl`: basic weak form DG discretization on a uniform triangular mesh. 
+* `elixir_euler_periodic_triangular_mesh.jl`: same as above, but enforces periodicity in the ``x,y`` directions.
+* `elixir_euler_triangulate_pkg_mesh.jl`: uses a `TriangulateIO` unstructured mesh generated by `Triangulate.jl`. 
+* `elixir_ape_sbp_triangular_mesh.jl`: uses a multi-dimensional SBP discretization in weak form. 
+* `elixir_euler_tet_mesh.jl`: basic weak form DG discretization on a uniform tet mesh.
diff --git a/examples/simplicial_2d_dg/elixir_ape_sbp_triangular_mesh.jl b/examples/simplicial_2d_dg/elixir_ape_sbp_triangular_mesh.jl
@@ -0,0 +1,47 @@
+# !!! warning "Experimental features"
+
+using StartUpDG, StructArrays
+using Trixi, OrdinaryDiffEq
+
+polydeg = 3
+rd = RefElemData(Tri(), SBP(), polydeg)
+dg = DG(rd, nothing #= mortar =#, 
+        SurfaceIntegralWeakForm(FluxLaxFriedrichs()), VolumeIntegralWeakForm())
+
+v_mean_global = (0.25, 0.25)
+c_mean_global = 1.0
+rho_mean_global = 1.0
+equations = AcousticPerturbationEquations2D(v_mean_global, c_mean_global, rho_mean_global)
+
+initial_condition = initial_condition_convergence_test
+source_terms = source_terms_convergence_test
+
+# example where we tag two separate boundary segments of the mesh
+VX, VY, EToV = StartUpDG.uniform_mesh(Tri(), 8)
+mesh = VertexMappedMesh(VX, VY, EToV, rd)
+
+boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = (; :entire_boundary => boundary_condition_convergence_test)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
+                                    source_terms = source_terms, 
+                                    boundary_conditions = boundary_conditions) 
+
+tspan = (0.0, 0.1)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+dt0 = StartUpDG.estimate_h(rd,mesh.md) / StartUpDG.inverse_trace_constant(rd)
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt = .5*dt0, save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
+
+l2,linf = analysis_callback(sol)
diff --git a/examples/simplicial_2d_dg/elixir_euler_periodic_triangular_mesh.jl b/examples/simplicial_2d_dg/elixir_euler_periodic_triangular_mesh.jl
@@ -0,0 +1,37 @@
+# !!! warning "Experimental features"
+
+using StartUpDG, StructArrays
+using Trixi, OrdinaryDiffEq
+
+polydeg = 3
+rd = RefElemData(Tri(), polydeg) # equivalent to a "basis"
+dg = DG(rd, nothing #= mortar =#, 
+        SurfaceIntegralWeakForm(FluxHLL()), VolumeIntegralWeakForm())
+
+equations = CompressibleEulerEquations2D(1.4)
+initial_condition = initial_condition_convergence_test
+source_terms = source_terms_convergence_test
+
+VX, VY, EToV = StartUpDG.uniform_mesh(rd.elementType, 8)
+mesh = VertexMappedMesh(VX, VY, EToV, rd, is_periodic=(true,true))
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
+                                    source_terms = source_terms) 
+
+tspan = (0.0, .1)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+dt0 = StartUpDG.estimate_h(rd,mesh.md) / StartUpDG.inverse_trace_constant(rd)
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt = .5*dt0, save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
+
+l2,linf = analysis_callback(sol)
diff --git a/examples/simplicial_2d_dg/elixir_euler_quadrilateral_mesh.jl b/examples/simplicial_2d_dg/elixir_euler_quadrilateral_mesh.jl
@@ -0,0 +1,47 @@
+# !!! warning "Experimental features"
+
+using StartUpDG, StructArrays
+using Trixi, OrdinaryDiffEq
+
+polydeg = 3
+rd = RefElemData(Quad(), polydeg)
+dg = DG(rd, nothing #= mortar =#, 
+        SurfaceIntegralWeakForm(FluxHLL()), VolumeIntegralWeakForm())
+
+equations = CompressibleEulerEquations2D(1.4)
+initial_condition = initial_condition_convergence_test
+source_terms = source_terms_convergence_test
+
+# example where we tag two separate boundary segments of the mesh
+top_boundary(x,y,tol=50*eps()) = abs(y-1)<tol 
+rest_of_boundary(x,y,tol=50*eps()) = !top_boundary(x,y,tol)
+is_on_boundary = Dict(:top => top_boundary, :rest => rest_of_boundary)
+VX, VY, EToV = StartUpDG.uniform_mesh(rd.elementType, 8)
+mesh = VertexMappedMesh(VX, VY, EToV, rd, is_on_boundary = is_on_boundary)
+
+boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = (; :top => boundary_condition_convergence_test,
+                        :rest => boundary_condition_convergence_test)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
+                                    source_terms = source_terms, 
+                                    boundary_conditions = boundary_conditions) 
+
+tspan = (0.0, .1)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+dt0 = StartUpDG.estimate_h(rd,mesh.md) / StartUpDG.inverse_trace_constant(rd)
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt = .5*dt0, save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
+
+l2,linf = analysis_callback(sol)
diff --git a/examples/simplicial_2d_dg/elixir_euler_triangular_mesh.jl b/examples/simplicial_2d_dg/elixir_euler_triangular_mesh.jl
@@ -0,0 +1,47 @@
+# !!! warning "Experimental features"
+
+using StartUpDG, StructArrays
+using Trixi, OrdinaryDiffEq
+
+polydeg = 3
+rd = RefElemData(Tri(), polydeg)
+dg = DG(rd, nothing #= mortar =#, 
+        SurfaceIntegralWeakForm(FluxHLL()), VolumeIntegralWeakForm())
+
+equations = CompressibleEulerEquations2D(1.4)
+initial_condition = initial_condition_convergence_test
+source_terms = source_terms_convergence_test
+
+# example where we tag two separate boundary segments of the mesh
+top_boundary(x,y,tol=50*eps()) = abs(y-1)<tol 
+rest_of_boundary(x,y,tol=50*eps()) = !top_boundary(x,y,tol)
+is_on_boundary = Dict(:top => top_boundary, :rest => rest_of_boundary)
+VX, VY, EToV = StartUpDG.uniform_mesh(Tri(), 8)
+mesh = VertexMappedMesh(VX, VY, EToV, rd, is_on_boundary = is_on_boundary)
+
+boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = (; :top => boundary_condition_convergence_test,
+                        :rest => boundary_condition_convergence_test)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
+                                    source_terms = source_terms, 
+                                    boundary_conditions = boundary_conditions) 
+
+tspan = (0.0, .1)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+dt0 = StartUpDG.estimate_h(rd,mesh.md) / StartUpDG.inverse_trace_constant(rd)
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt = .5*dt0, save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
+
+l2,linf = analysis_callback(sol)