hmm

examples/CMakeLists.txt

@@ -89,6 +89,7 @@ if (MFEM_USE_MPI)
     ex37p.cpp
     ex39p.cpp
     ex40p.cpp
+    swe.cpp
   )
 endif()
 

examples/makefile

@@ -26,7 +26,7 @@ SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
 PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \
    ex12p ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
    ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p ex34p ex35p ex36p \
-   ex37p ex39p ex40p
+   ex37p ex39p ex40p swe
 SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
 PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p \
    ex22p ex24p ex25p ex26p ex34p ex35p
@@ -102,6 +102,7 @@ ifeq ($(MFEM_USE_MPI),YES)
 ex18p: $(SRC)ex18.hpp
 ex33p: $(SRC)ex33.hpp
 ex37p: $(SRC)ex37.hpp
+swe: $(SRC)swe.hpp
 endif
 
 MFEM_TESTS = EXAMPLES

examples/swe.cpp

@@ -0,0 +1,365 @@
+//                                MFEM Example 18 - Parallel Version
+//
+// Compile with: make ex18p
+//
+// Sample runs:
+//
+//       mpirun -np 4 ex18p -p 1 -rs 2 -rp 1 -o 1 -s 3
+//       mpirun -np 4 ex18p -p 1 -rs 1 -rp 1 -o 3 -s 4
+//       mpirun -np 4 ex18p -p 1 -rs 1 -rp 1 -o 5 -s 6
+//       mpirun -np 4 ex18p -p 2 -rs 1 -rp 1 -o 1 -s 3 -mf
+//       mpirun -np 4 ex18p -p 2 -rs 1 -rp 1 -o 3 -s 3 -mf
+//
+// Description:  This example code solves the compressible swe system of
+//               equations, a model nonlinear hyperbolic PDE, with a
+//               discontinuous Galerkin (DG) formulation in parallel.
+//
+//                (u_t, v)_T - (F(u), ∇ v)_T + <F̂(u,n), [[v]]>_F = 0
+//
+//               where (⋅,⋅)_T is volume integration, and <⋅,⋅>_F is face
+//               integration, F is the swe flux function, and F̂ is the
+//               numerical flux.
+//
+//               Specifically, it solves for an exact solution of the equations
+//               whereby a vortex is transported by a uniform flow. Since all
+//               boundaries are periodic here, the method's accuracy can be
+//               assessed by measuring the difference between the solution and
+//               the initial condition at a later time when the vortex returns
+//               to its initial location.
+//
+//               Note that as the order of the spatial discretization increases,
+//               the timestep must become smaller. This example currently uses a
+//               simple estimate derived by Cockburn and Shu for the 1D RKDG
+//               method. An additional factor can be tuned by passing the --cfl
+//               (or -c shorter) flag.
+//
+//               The example demonstrates usage of DGHyperbolicConservationLaws
+//               that wraps NonlinearFormIntegrators containing element and face
+//               integration schemes. In this case the system also involves an
+//               external approximate Riemann solver for the DG interface flux.
+//               By default, weak-divergence is pre-assembled in element-wise
+//               manner, which corresponds to (I_h(F(u_h)), ∇ v). This yields
+//               better performance and similar accuracy for the included test
+//               problems. This can be turned off and use nonlinear assembly
+//               similar to matrix-free assembly when -mf flag is provided.
+//               It also demonstrates how to use GLVis for in-situ visualization
+//               of vector grid function and how to set top-view.
+//
+//               We recommend viewing examples 9, 14 and 17 before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include "swe.hpp"
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 0. Parallel setup
+   Mpi::Init(argc, argv);
+   const int numProcs = Mpi::WorldSize();
+   const int myRank = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 1. Parse command-line options.
+   int problem = 1;
+   const real_t specific_heat_ratio = 1.4;
+   const real_t gas_constant = 1.0;
+
+   string mesh_file = "";
+   int IntOrderOffset = 1;
+   int ser_ref_levels = 0;
+   int par_ref_levels = 1;
+   int order = 3;
+   int ode_solver_type = 4;
+   real_t t_final = 1.5;
+   real_t dt = -0.01;
+   real_t cfl = 0.3;
+   bool visualization = true;
+   bool preassembleWeakDiv = false;
+   int vis_steps = 50;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use. If not provided, then a periodic square"
+                  " mesh will be used.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See swe().");
+   args.AddOption(&ser_ref_levels, "-rs", "--serial-refine",
+                  "Number of times to refine the serial mesh uniformly.");
+   args.AddOption(&par_ref_levels, "-rp", "--parallel-refine",
+                  "Number of times to refine the parallel mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Forward swe,\n\t"
+                  "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
+   args.AddOption(&t_final, "-tf", "--t-final", "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step. Positive number skips CFL timestep calculation.");
+   args.AddOption(&cfl, "-c", "--cfl-number",
+                  "CFL number for timestep calculation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&preassembleWeakDiv, "-ea", "--element-assembly-divergence",
+                  "-mf", "--matrix-free-divergence",
+                  "Weak divergence assembly level\n"
+                  "    ea - Element assembly with interpolated F\n"
+                  "    mf - Nonlinear assembly in matrix-free manner");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.ParseCheck();
+
+   // 2. Read the mesh from the given mesh file. When the user does not provide
+   //    mesh file, use the default mesh file for the problem.
+   Mesh mesh = mesh_file.empty() ? SWEMesh(problem) : Mesh(mesh_file);
+   const int dim = mesh.Dimension();
+   const int num_equations = dim + 1;
+
+   // Refine the mesh to increase the resolution. In this example we do
+   // 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is a
+   // command-line parameter.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // Define a parallel mesh by a partitioning of the serial mesh. Refine this
+   // mesh further in parallel to increase the resolution. Once the parallel
+   // mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh = ParMesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   // Refine the mesh to increase the resolution. In this example we do
+   // 'par_ref_levels' of uniform refinement, where 'par_ref_levels' is a
+   // command-line parameter.
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh.UniformRefinement();
+   }
+
+   // 3. Define the ODE solver used for time integration. Several explicit
+   //    Runge-Kutta methods are available.
+   ODESolver *ode_solver = NULL;
+   switch (ode_solver_type)
+   {
+      case 1: ode_solver = new ForwardEulerSolver; break;
+      case 2: ode_solver = new RK2Solver(1.0); break;
+      case 3: ode_solver = new RK3SSPSolver; break;
+      case 4: ode_solver = new RK4Solver; break;
+      case 6: ode_solver = new RK6Solver; break;
+      default:
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
+   }
+
+   // 4. Define the discontinuous DG finite element space of the given
+   //    polynomial order on the refined mesh.
+   DG_FECollection fec(order, dim);
+   // Finite element space for a scalar (thermodynamic quantity)
+   ParFiniteElementSpace fes(&pmesh, &fec);
+   // Finite element space for a mesh-dim vector quantity (momentum)
+   ParFiniteElementSpace dfes(&pmesh, &fec, dim, Ordering::byNODES);
+   // Finite element space for all variables together (total thermodynamic state)
+   ParFiniteElementSpace vfes(&pmesh, &fec, num_equations, Ordering::byNODES);
+
+   // This example depends on this ordering of the space.
+   MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES, "");
+
+   HYPRE_BigInt glob_size = vfes.GlobalTrueVSize();
+   if (Mpi::Root())
+   {
+      cout << "Number of unknowns: " << glob_size << endl;
+   }
+
+   // 5. Define the initial conditions, save the corresponding mesh and grid
+   //    functions to files. These can be opened with GLVis using:
+   //    "glvis -np 4  -m swe-mesh -g swe-1-init" (for x-momentum).
+
+   // Initialize the state.
+   VectorFunctionCoefficient u0 = SWEInitialCondition(problem,
+                                                      specific_heat_ratio,
+                                                      gas_constant);
+   ParGridFunction sol(&vfes);
+   sol.ProjectCoefficient(u0);
+   ParGridFunction mom(&dfes, sol.GetData() + fes.GetTrueVSize());
+   ParGridFunction height(&fes, sol.GetData());
+   // Output the initial solution.
+   {
+      ostringstream mesh_name;
+      mesh_name << "swe-mesh." << setfill('0') << setw(6) << Mpi::WorldRank();
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(precision);
+      mesh_ofs << pmesh;
+
+      for (int k = 0; k < num_equations; k++)
+      {
+         ParGridFunction uk(&fes, sol.GetData() + k * fes.GetNDofs());
+         ostringstream sol_name;
+         sol_name << "swe-" << k << "-init." << setfill('0') << setw(6)
+                  << Mpi::WorldRank();
+         ofstream sol_ofs(sol_name.str().c_str());
+         sol_ofs.precision(precision);
+         sol_ofs << uk;
+      }
+   }
+
+   // 6. Set up the nonlinear form with swe flux and numerical flux
+   ShallowWaterFlux flux(dim, specific_heat_ratio);
+   RusanovFlux numericalFlux(flux);
+   DGHyperbolicConservationLaws swe(
+      vfes, std::unique_ptr<HyperbolicFormIntegrator>(
+         new HyperbolicFormIntegrator(numericalFlux, IntOrderOffset)),
+      preassembleWeakDiv);
+
+   // 7. Visualize momentum with its magnitude
+   socketstream height_sock, mom_sock;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int visport = 19916;
+
+      height_sock.open(vishost, visport);
+      height_sock.precision(precision);
+      // Plot magnitude of vector-valued momentum
+      height_sock << "parallel " << numProcs << " " << myRank << "\n";
+      height_sock << "solution\n" << pmesh << height;
+      height_sock << "window_title 'height, t = 0'\n";
+      height_sock << "autoscale off\n valuerange 0.5 2.5\n";
+      height_sock << "view 0 0\n";  // view from top
+      height_sock << "keys jlm\n";  // turn off perspective and light, show mesh
+      height_sock << flush;
+      MPI_Barrier(pmesh.GetComm());
+
+      mom_sock.open(vishost, visport);
+      mom_sock.precision(precision);
+      // Plot magnitude of vector-valued momentum
+      mom_sock << "parallel " << numProcs << " " << myRank << "\n";
+      mom_sock << "solution\n" << pmesh << mom;
+      mom_sock << "window_title 'mom, t = 0'\n";
+      mom_sock << "autoscale off\n valuerange 0.0 0.5\n";
+      mom_sock << "view 0 0\n";  // view from top
+      mom_sock << "keys jlm\n";  // turn off perspective and light, show mesh
+      mom_sock << flush;
+      MPI_Barrier(pmesh.GetComm());
+   }
+
+   // 8. Time integration
+
+   // When dt is not specified, use CFL condition.
+   // Compute h_min and initial maximum characteristic speed
+   real_t hmin = infinity();
+   if (cfl > 0)
+   {
+      for (int i = 0; i < pmesh.GetNE(); i++)
+      {
+         hmin = min(pmesh.GetElementSize(i, 1), hmin);
+      }
+      MPI_Allreduce(MPI_IN_PLACE, &hmin, 1,  MPITypeMap<real_t>::mpi_type, MPI_MIN,
+                    pmesh.GetComm());
+      // Find a safe dt, using a temporary vector. Calling Mult() computes the
+      // maximum char speed at all quadrature points on all faces (and all
+      // elements with -mf).
+      Vector z(sol.Size());
+      swe.Mult(sol, z);
+
+      real_t max_char_speed = swe.GetMaxCharSpeed();
+      MPI_Allreduce(MPI_IN_PLACE, &max_char_speed, 1,  MPITypeMap<real_t>::mpi_type,
+                    MPI_MAX,
+                    pmesh.GetComm());
+      dt = cfl * hmin / max_char_speed / (2 * order + 1);
+   }
+
+   // Start the timer.
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   // Init time integration
+   real_t t = 0.0;
+   swe.SetTime(t);
+   ode_solver->Init(swe);
+
+   // Integrate in time.
+   bool done = false;
+   for (int ti = 0; !done;)
+   {
+      real_t dt_real = min(dt, t_final - t);
+
+      ode_solver->Step(sol, t, dt_real);
+      if (cfl > 0) // update time step size with CFL
+      {
+         real_t max_char_speed = swe.GetMaxCharSpeed();
+         MPI_Allreduce(MPI_IN_PLACE, &max_char_speed, 1,  MPITypeMap<real_t>::mpi_type,
+                       MPI_MAX,
+                       pmesh.GetComm());
+         dt = cfl * hmin / max_char_speed / (2 * order + 1);
+      }
+      ti++;
+
+      done = (t >= t_final - 1e-8 * dt);
+      if (done || ti % vis_steps == 0)
+      {
+         if (Mpi::Root())
+         {
+            cout << "time step: " << ti << ", time: " << t << endl;
+         }
+         if (visualization)
+         {
+            height_sock << "window_title 'height, t = " << t << "'\n";
+            height_sock << "parallel " << numProcs << " " << myRank << "\n";
+            height_sock << "solution\n" << pmesh << height << flush;
+            mom_sock << "window_title 'mom, t = " << t << "'\n";
+            mom_sock << "parallel " << numProcs << " " << myRank << "\n";
+            mom_sock << "solution\n" << pmesh << mom << flush;
+         }
+      }
+   }
+
+   tic_toc.Stop();
+   if (Mpi::Root())
+   {
+      cout << " done, " << tic_toc.RealTime() << "s." << endl;
+   }
+
+   // 9. Save the final solution. This output can be viewed later using GLVis:
+   //    "glvis -np 4 -m swe-mesh-final -g swe-1-final" (for x-momentum).
+   {
+      ostringstream mesh_name;
+      mesh_name << "swe-mesh-final." << setfill('0') << setw(6)
+                << Mpi::WorldRank();
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(precision);
+      mesh_ofs << pmesh;
+
+      for (int k = 0; k < num_equations; k++)
+      {
+         ParGridFunction uk(&fes, sol.GetData() + k * fes.GetNDofs());
+         ostringstream sol_name;
+         sol_name << "swe-" << k << "-final." << setfill('0') << setw(6)
+                  << Mpi::WorldRank();
+         ofstream sol_ofs(sol_name.str().c_str());
+         sol_ofs.precision(precision);
+         sol_ofs << uk;
+      }
+   }
+
+   // 10. Compute the L2 solution error summed for all components.
+   const real_t error = sol.ComputeLpError(2, u0);
+   if (Mpi::Root())
+   {
+      cout << "Solution error: " << error << endl;
+   }
+
+   // Free the used memory.
+   delete ode_solver;
+
+   return 0;
+}

examples/swe.hpp

@@ -0,0 +1,346 @@
+//                  MFEM Example 18 - Serial/Parallel Shared Code
+//                      (Implementation of Time-dependent DG Operator)
+//
+// This code provide example problems for the SWE equations and implements
+// the time-dependent DG operator given by the equation:
+//
+//            (u_t, v)_T - (F(u), ∇ v)_T + <F̂(u, n), [[v]]>_F = 0.
+//
+// This operator is designed for explicit time stepping methods. Specifically,
+// the function DGHyperbolicConservationLaws::Mult implements the following
+// transformation:
+//
+//                             u ↦ M⁻¹(-DF(u) + NF(u))
+//
+// where M is the mass matrix, DF is the weak divergence of flux, and NF is the
+// interface flux. The inverse of the mass matrix is computed element-wise by
+// leveraging the block-diagonal structure of the DG mass matrix. Additionally,
+// the flux-related terms are computed using the HyperbolicFormIntegrator.
+//
+// The maximum characteristic speed is determined for each time step. For more
+// details, refer to the documentation of DGHyperbolicConservationLaws::Mult.
+//
+
+#include <functional>
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+/// @brief Time dependent DG operator for hyperbolic conservation laws
+class DGHyperbolicConservationLaws : public TimeDependentOperator
+{
+private:
+   const int num_equations; // the number of equations
+   const int dim;
+   FiniteElementSpace &vfes; // vector finite element space
+   // Element integration form. Should contain ComputeFlux
+   std::unique_ptr<HyperbolicFormIntegrator> formIntegrator;
+   // Base Nonlinear Form
+   std::unique_ptr<NonlinearForm> nonlinearForm;
+   // element-wise inverse mass matrix
+   std::vector<DenseMatrix> invmass; // local scalar inverse mass
+   std::vector<DenseMatrix> weakdiv; // local weak divergence (trial space ByDim)
+   // global maximum characteristic speed. Updated by form integrators
+   mutable real_t max_char_speed;
+   // auxiliary variable used in Mult
+   mutable Vector z;
+
+   // Compute element-wise inverse mass matrix
+   void ComputeInvMass();
+   // Compute element-wise weak-divergence matrix
+   void ComputeWeakDivergence();
+
+public:
+   /**
+    * @brief Construct a new DGHyperbolicConservationLaws object
+    *
+    * @param vfes_ vector finite element space. Only tested for DG [Pₚ]ⁿ
+    * @param formIntegrator_ integrator (F(u,x), grad v)
+    * @param preassembleWeakDivergence preassemble weak divergence for faster
+    *                                  assembly
+    */
+   DGHyperbolicConservationLaws(
+      FiniteElementSpace &vfes_,
+      std::unique_ptr<HyperbolicFormIntegrator> formIntegrator_,
+      bool preassembleWeakDivergence=true);
+   /**
+    * @brief Apply nonlinear form to obtain M⁻¹(DIVF + JUMP HAT(F))
+    *
+    * @param x current solution vector
+    * @param y resulting dual vector to be used in an EXPLICIT solver
+    */
+   void Mult(const Vector &x, Vector &y) const override;
+   // get global maximum characteristic speed to be used in CFL condition
+   // where max_char_speed is updated during Mult.
+   real_t GetMaxCharSpeed() { return max_char_speed; }
+   void Update();
+
+};
+
+//////////////////////////////////////////////////////////////////
+///        HYPERBOLIC CONSERVATION LAWS IMPLEMENTATION         ///
+//////////////////////////////////////////////////////////////////
+
+// Implementation of class DGHyperbolicConservationLaws
+DGHyperbolicConservationLaws::DGHyperbolicConservationLaws(
+   FiniteElementSpace &vfes_,
+   std::unique_ptr<HyperbolicFormIntegrator> formIntegrator_,
+   bool preassembleWeakDivergence)
+   : TimeDependentOperator(vfes_.GetTrueVSize()),
+     num_equations(formIntegrator_->num_equations),
+     dim(vfes_.GetMesh()->SpaceDimension()),
+     vfes(vfes_),
+     formIntegrator(std::move(formIntegrator_)),
+     z(vfes_.GetTrueVSize())
+{
+   // Standard local assembly and inversion for energy mass matrices.
+   ComputeInvMass();
+#ifndef MFEM_USE_MPI
+   nonlinearForm.reset(new NonlinearForm(&vfes));
+#else
+   ParFiniteElementSpace *pvfes = dynamic_cast<ParFiniteElementSpace *>(&vfes);
+   if (pvfes)
+   {
+      nonlinearForm.reset(new ParNonlinearForm(pvfes));
+   }
+   else
+   {
+      nonlinearForm.reset(new NonlinearForm(&vfes));
+   }
+#endif
+   if (preassembleWeakDivergence)
+   {
+      ComputeWeakDivergence();
+   }
+   else
+   {
+      nonlinearForm->AddDomainIntegrator(formIntegrator.get());
+   }
+   nonlinearForm->AddInteriorFaceIntegrator(formIntegrator.get());
+   nonlinearForm->UseExternalIntegrators();
+
+}
+
+void DGHyperbolicConservationLaws::ComputeInvMass()
+{
+   InverseIntegrator inv_mass(new MassIntegrator());
+
+   invmass.resize(vfes.GetNE());
+   for (int i=0; i<vfes.GetNE(); i++)
+   {
+      int dof = vfes.GetFE(i)->GetDof();
+      invmass[i].SetSize(dof);
+      inv_mass.AssembleElementMatrix(*vfes.GetFE(i),
+                                     *vfes.GetElementTransformation(i),
+                                     invmass[i]);
+   }
+}
+
+void DGHyperbolicConservationLaws::ComputeWeakDivergence()
+{
+   TransposeIntegrator weak_div(new GradientIntegrator());
+   DenseMatrix weakdiv_bynodes;
+
+   weakdiv.resize(vfes.GetNE());
+   for (int i=0; i<vfes.GetNE(); i++)
+   {
+      int dof = vfes.GetFE(i)->GetDof();
+      weakdiv_bynodes.SetSize(dof, dof*dim);
+      weak_div.AssembleElementMatrix2(*vfes.GetFE(i), *vfes.GetFE(i),
+                                      *vfes.GetElementTransformation(i),
+                                      weakdiv_bynodes);
+      weakdiv[i].SetSize(dof, dof*dim);
+      // Reorder so that trial space is ByDim.
+      // This makes applying weak divergence to flux value simpler.
+      for (int j=0; j<dof; j++)
+      {
+         for (int d=0; d<dim; d++)
+         {
+            weakdiv[i].SetCol(j*dim + d, weakdiv_bynodes.GetColumn(d*dof + j));
+         }
+      }
+
+   }
+}
+
+
+void DGHyperbolicConservationLaws::Mult(const Vector &x, Vector &y) const
+{
+   // 0. Reset wavespeed computation before operator application.
+   formIntegrator->ResetMaxCharSpeed();
+   // 1. Apply Nonlinear form to obtain an auxiliary result
+   //         z = - <F̂(u_h,n), [[v]]>_e
+   //    If weak-divergence is not preassembled, we also have weak-divergence
+   //         z = - <F̂(u_h,n), [[v]]>_e + (F(u_h), ∇v)
+   nonlinearForm->Mult(x, z);
+   if (!weakdiv.empty()) // if weak divergence is pre-assembled
+   {
+      // Apply weak divergence to F(u_h), and inverse mass to z_loc + weakdiv_loc
+      Vector current_state; // view of current state at a node
+      DenseMatrix current_flux; // flux of current state
+      DenseMatrix flux; // element flux value. Whose column is ordered by dim.
+      DenseMatrix current_xmat; // view of current states in an element, dof x num_eq
+      DenseMatrix current_zmat; // view of element auxiliary result, dof x num_eq
+      DenseMatrix current_ymat; // view of element result, dof x num_eq
+      const FluxFunction &fluxFunction = formIntegrator->GetFluxFunction();
+      Array<int> vdofs;
+      Vector xval, zval;
+      for (int i=0; i<vfes.GetNE(); i++)
+      {
+         ElementTransformation* Tr = vfes.GetElementTransformation(i);
+         int dof = vfes.GetFE(i)->GetDof();
+         vfes.GetElementVDofs(i, vdofs);
+         x.GetSubVector(vdofs, xval);
+         current_xmat.UseExternalData(xval.GetData(), dof, num_equations);
+         flux.SetSize(num_equations, dim*dof);
+         for (int j=0; j<dof; j++) // compute flux for all nodes in the element
+         {
+            current_xmat.GetRow(j, current_state);
+            current_flux.UseExternalData(flux.GetData() + num_equations*dim*j,
+                                         num_equations, dof);
+            fluxFunction.ComputeFlux(current_state, *Tr, current_flux);
+         }
+         // Compute weak-divergence and add it to auxiliary result, z
+         // Recalling that weakdiv is reordered by dim, we can apply
+         // weak-divergence to the transpose of flux.
+         z.GetSubVector(vdofs, zval);
+         current_zmat.UseExternalData(zval.GetData(), dof, num_equations);
+         mfem::AddMult_a_ABt(1.0, weakdiv[i], flux, current_zmat);
+         // Apply inverse mass to auxiliary result to obtain the final result
+         current_ymat.SetSize(dof, num_equations);
+         mfem::Mult(invmass[i], current_zmat, current_ymat);
+         y.SetSubVector(vdofs, current_ymat.GetData());
+      }
+   }
+   else
+   {
+      // Apply block inverse mass
+      Vector zval; // z_loc, dof*num_eq
+
+      DenseMatrix current_zmat; // view of element auxiliary result, dof x num_eq
+      DenseMatrix current_ymat; // view of element result, dof x num_eq
+      Array<int> vdofs;
+      for (int i=0; i<vfes.GetNE(); i++)
+      {
+         int dof = vfes.GetFE(i)->GetDof();
+         vfes.GetElementVDofs(i, vdofs);
+         z.GetSubVector(vdofs, zval);
+         current_zmat.UseExternalData(zval.GetData(), dof, num_equations);
+         current_ymat.SetSize(dof, num_equations);
+         mfem::Mult(invmass[i], current_zmat, current_ymat);
+         y.SetSubVector(vdofs, current_ymat.GetData());
+      }
+   }
+   max_char_speed = formIntegrator->GetMaxCharSpeed();
+}
+
+void DGHyperbolicConservationLaws::Update()
+{
+   nonlinearForm->Update();
+   height = nonlinearForm->Height();
+   width = height;
+   z.SetSize(height);
+
+   ComputeInvMass();
+   if (!weakdiv.empty()) {ComputeWeakDivergence();}
+}
+
+std::function<void(const Vector&, Vector&)> GetMovingVortexInit(
+   const real_t radius, const real_t Minf, const real_t beta,
+   const real_t gas_constant, const real_t specific_heat_ratio)
+{
+   return [specific_heat_ratio,
+           gas_constant, Minf, radius, beta](const Vector &x, Vector &y)
+   {
+      MFEM_ASSERT(x.Size() == 2, "");
+
+      const real_t xc = 0.0, yc = 0.0;
+
+      // Nice units
+      const real_t vel_inf = 1.;
+      const real_t den_inf = 1.;
+
+      // Derive remainder of background state from this and Minf
+      const real_t pres_inf = (den_inf / specific_heat_ratio) *
+                              (vel_inf / Minf) * (vel_inf / Minf);
+      const real_t temp_inf = pres_inf / (den_inf * gas_constant);
+
+      real_t r2rad = 0.0;
+      r2rad += (x(0) - xc) * (x(0) - xc);
+      r2rad += (x(1) - yc) * (x(1) - yc);
+      r2rad /= (radius * radius);
+
+      const real_t shrinv1 = 1.0 / (specific_heat_ratio - 1.);
+
+      const real_t velX =
+         vel_inf * (1 - beta * (x(1) - yc) / radius * std::exp(-0.5 * r2rad));
+      const real_t velY =
+         vel_inf * beta * (x(0) - xc) / radius * std::exp(-0.5 * r2rad);
+      const real_t vel2 = velX * velX + velY * velY;
+
+      const real_t specific_heat =
+         gas_constant * specific_heat_ratio * shrinv1;
+      const real_t temp = temp_inf - 0.5 * (vel_inf * beta) *
+                          (vel_inf * beta) / specific_heat *
+                          std::exp(-r2rad);
+
+      const real_t den = den_inf * std::pow(temp / temp_inf, shrinv1);
+      const real_t pres = den * gas_constant * temp;
+      const real_t energy = shrinv1 * pres / den + 0.5 * vel2;
+
+      y(0) = den;
+      y(1) = den * velX;
+      y(2) = den * velY;
+      y(3) = den * energy;
+   };
+}
+
+Mesh SWEMesh(const int problem)
+{
+   switch (problem)
+   {
+      case 1:
+      case 2:
+      case 3:
+      {
+         Mesh mesh("../data/periodic-square.mesh");
+         return mesh;
+         break;
+      }
+      case 4:
+         return Mesh("../data/periodic-segment.mesh");
+         break;
+      default:
+         MFEM_ABORT("Problem Undefined");
+   }
+}
+
+// Initial condition
+VectorFunctionCoefficient SWEInitialCondition(const int problem,
+                                              const real_t specific_heat_ratio,
+                                              const real_t gas_constant)
+{
+   switch (problem)
+   {
+      case 1: // steady flow
+         return VectorFunctionCoefficient(3, [](const Vector &x, Vector &y)
+         {
+            y = 0.0;
+            y[0] = 1.0;
+         });
+      case 2:
+         return VectorFunctionCoefficient(3, [](const Vector &x, Vector &y)
+         {
+            const real_t hmin = 1.0;
+            const real_t hmax = 2.0;
+            const real_t sigma = 0.2;
+            y = 0.0;
+            y[0]= hmin + (hmax - hmin) * std::exp(-(x*x)/(2*sigma*sigma));
+         });
+      default:
+         MFEM_ABORT("Problem Undefined");
+   }
+}
+
+} // namespace mfem

fem/bilininteg.cpp

@@ -774,16 +774,17 @@ void GradientIntegrator::AssembleElementMatrix2(
    ElementTransformation &Trans,  DenseMatrix &elmat)
 {
    dim = test_fe.GetDim();
+   const int sdim = Trans.GetSpaceDim();
    int trial_dof = trial_fe.GetDof();
    int test_dof = test_fe.GetDof();
    real_t c;
    Vector d_col;
 
    dshape.SetSize(trial_dof, dim);
-   gshape.SetSize(trial_dof, dim);
-   Jadj.SetSize(dim);
+   gshape.SetSize(trial_dof, sdim);
+   Jadj.SetSize(dim, sdim);
    shape.SetSize(test_dof);
-   elmat.SetSize(dim * test_dof, trial_dof);
+   elmat.SetSize(sdim * test_dof, trial_dof);
 
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(trial_fe, test_fe,
                                                             Trans);
@@ -799,9 +800,7 @@ void GradientIntegrator::AssembleElementMatrix2(
       CalcAdjugate(Trans.Jacobian(), Jadj);
 
       test_fe.CalcPhysShape(Trans, shape);
-      trial_fe.CalcDShape(ip, dshape);
-
-      Mult(dshape, Jadj, gshape);
+      trial_fe.CalcPhysDShape(Trans, gshape);
 
       c = ip.weight;
       if (Q)
@@ -810,7 +809,7 @@ void GradientIntegrator::AssembleElementMatrix2(
       }
       shape *= c;
 
-      for (int d = 0; d < dim; ++d)
+      for (int d = 0; d < sdim; ++d)
       {
          gshape.GetColumnReference(d, d_col);
          MultVWt(shape, d_col, elmat_comp);

fem/fe/fe_l2.cpp

@@ -28,6 +28,7 @@ L2_SegmentElement::L2_SegmentElement(const int p, const int btype)
 #ifndef MFEM_THREAD_SAFE
    shape_x.SetSize(p + 1);
    dshape_x.SetDataAndSize(NULL, p + 1);
+   d2shape_x.SetSize(p+1);
 #endif
 
    for (int i = 0; i <= p; i++)
@@ -54,6 +55,24 @@ void L2_SegmentElement::CalcDShape(const IntegrationPoint &ip,
    basis1d.ScaleIntegrated(map_type == VALUE);
    basis1d.Eval(ip.x, shape_x, dshape_x);
 }
+void L2_SegmentElement::CalcHessian(const IntegrationPoint &ip,
+                                    DenseMatrix &Hessian) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(p+1), dshape_x(p+1), d2shape_x(p+1);
+#endif
+
+   basis1d.Eval(ip.x, shape_x, dshape_x, d2shape_x);
+
+   Hessian(0,0) = d2shape_x(0);
+   Hessian(1,0) = d2shape_x(p);
+   for (int i = 1; i < p; i++)
+   {
+      Hessian(i+1,0) = d2shape_x(i);
+   }
+}
 
 void L2_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
 {
@@ -89,6 +108,8 @@ L2_QuadrilateralElement::L2_QuadrilateralElement(const int p, const int btype)
    shape_y.SetSize(p + 1);
    dshape_x.SetSize(p + 1);
    dshape_y.SetSize(p + 1);
+   d2shape_x.SetSize(p + 1);
+   d2shape_y.SetSize(p + 1);
 #endif
 
    for (int o = 0, j = 0; j <= p; j++)
@@ -139,6 +160,30 @@ void L2_QuadrilateralElement::CalcDShape(const IntegrationPoint &ip,
       }
 }
 
+void L2_QuadrilateralElement::CalcHessian(const IntegrationPoint &ip,
+                                          DenseMatrix &Hessian) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1),
+          d2shape_x(p+1), d2shape_y(p+1);
+#endif
+
+   basis1d.Eval(ip.x, shape_x, dshape_x, d2shape_x);
+   basis1d.Eval(ip.y, shape_y, dshape_y, d2shape_y);
+
+   for (int o = 0, j = 0; j <= p; j++)
+   {
+      for (int i = 0; i <= p; i++)
+      {
+         Hessian(o,0) = d2shape_x(i)*  shape_y(j);
+         Hessian(o,1) =  dshape_x(i)* dshape_y(j);
+         Hessian(o,2) =   shape_x(i)*d2shape_y(j);  o++;
+      }
+   }
+}
+
 void L2_QuadrilateralElement::ProjectDelta(int vertex, Vector &dofs) const
 {
    const int p = order;
@@ -307,6 +352,9 @@ L2_HexahedronElement::L2_HexahedronElement(const int p, const int btype)
    dshape_x.SetSize(p + 1);
    dshape_y.SetSize(p + 1);
    dshape_z.SetSize(p + 1);
+   d2shape_x.SetSize(p + 1);
+   d2shape_y.SetSize(p + 1);
+   d2shape_z.SetSize(p + 1);
 #endif
 
    for (int o = 0, k = 0; k <= p; k++)
@@ -364,6 +412,35 @@ void L2_HexahedronElement::CalcDShape(const IntegrationPoint &ip,
          }
 }
 
+void L2_HexahedronElement::CalcHessian(const IntegrationPoint &ip,
+                                       DenseMatrix &Hessian) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(p+1),  shape_y(p+1),  shape_z(p+1);
+   Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
+   Vector d2shape_x(p+1), d2shape_y(p+1), d2shape_z(p+1);
+#endif
+
+   basis1d.Eval(ip.x, shape_x, dshape_x, d2shape_x);
+   basis1d.Eval(ip.y, shape_y, dshape_y, d2shape_y);
+   basis1d.Eval(ip.z, shape_z, dshape_z, d2shape_z);
+
+   for (int o = 0, k = 0; k <= p; k++)
+      for (int j = 0; j <= p; j++)
+         for (int i = 0; i <= p; i++)
+         {
+            Hessian(o,0) = d2shape_x(i)*  shape_y(j)*  shape_z(k);
+            Hessian(o,1) =  dshape_x(i)* dshape_y(j)*  shape_z(k);
+            Hessian(o,2) =  dshape_x(i)*  shape_y(j)* dshape_z(k);
+            Hessian(o,3) =   shape_x(i)*d2shape_y(j)*  shape_z(k);
+            Hessian(o,4) =   shape_x(i)* dshape_y(j)* dshape_z(k);
+            Hessian(o,5) =   shape_x(i)*  shape_y(j)*d2shape_z(k);
+            o++;
+         }
+}
+
 void L2_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
 {
    const int p = order;
@@ -580,8 +657,12 @@ L2_TriangleElement::L2_TriangleElement(const int p, const int btype)
    dshape_x.SetSize(p + 1);
    dshape_y.SetSize(p + 1);
    dshape_l.SetSize(p + 1);
+   ddshape_x.SetSize(p + 1);
+   ddshape_y.SetSize(p + 1);
+   ddshape_l.SetSize(p + 1);
    u.SetSize(dof);
    du.SetSize(dof, dim);
+   ddu.SetSize(dof, (dim * (dim + 1)) / 2 );
 #else
    Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
 #endif
@@ -663,6 +744,38 @@ void L2_TriangleElement::CalcDShape(const IntegrationPoint &ip,
    Ti.Mult(du, dshape);
 }
 
+void L2_TriangleElement::CalcHessian(const IntegrationPoint &ip,
+                                     DenseMatrix &ddshape) const
+{
+   const int p = order;
+#ifdef MFEM_THREAD_SAFE
+   Vector   shape_x(p + 1),   shape_y(p + 1),   shape_l(p + 1);
+   Vector  dshape_x(p + 1),  dshape_y(p + 1),  dshape_l(p + 1);
+   Vector ddshape_x(p + 1), ddshape_y(p + 1), ddshape_l(p + 1);
+   DenseMatrix ddu(dof, dim);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x, ddshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y, ddshape_y);
+   poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l, dshape_l, ddshape_l);
+
+   for (int o = 0, j = 0; j <= p; j++)
+      for (int i = 0; i + j <= p; i++)
+      {
+         int k = p - i - j;
+         // u_xx, u_xy, u_yy
+         ddu(o,0) = ((ddshape_x(i) * shape_l(k)) - 2. * (dshape_x(i) * dshape_l(k)) +
+                     (shape_x(i) * ddshape_l(k))) * shape_y(j);
+         ddu(o,1) = (((shape_x(i) * ddshape_l(k)) - dshape_x(i) * dshape_l(k)) * shape_y(
+                        j)) + (((dshape_x(i) * shape_l(k)) - (shape_x(i) * dshape_l(k))) * dshape_y(j));
+         ddu(o,2) = ((ddshape_y(j) * shape_l(k)) - 2. * (dshape_y(j) * dshape_l(k)) +
+                     (shape_y(j) * ddshape_l(k))) * shape_x(i);
+         o++;
+      }
+
+   Ti.Mult(ddu, ddshape);
+}
+
 void L2_TriangleElement::ProjectDelta(int vertex, Vector &dofs) const
 {
    switch (vertex)
@@ -707,8 +820,13 @@ L2_TetrahedronElement::L2_TetrahedronElement(const int p, const int btype)
    dshape_y.SetSize(p + 1);
    dshape_z.SetSize(p + 1);
    dshape_l.SetSize(p + 1);
+   ddshape_x.SetSize(p + 1);
+   ddshape_y.SetSize(p + 1);
+   ddshape_z.SetSize(p + 1);
+   ddshape_l.SetSize(p + 1);
    u.SetSize(dof);
    du.SetSize(dof, dim);
+   ddu.SetSize(dof, (dim * (dim + 1)) / 2);
 #else
    Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
 #endif
@@ -800,6 +918,52 @@ void L2_TetrahedronElement::CalcDShape(const IntegrationPoint &ip,
    Ti.Mult(du, dshape);
 }
 
+void L2_TetrahedronElement::CalcHessian(const IntegrationPoint &ip,
+                                        DenseMatrix &ddshape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector   shape_x(p + 1),   shape_y(p + 1),   shape_z(p + 1),   shape_l(p + 1);
+   Vector  dshape_x(p + 1),  dshape_y(p + 1),  dshape_z(p + 1),  dshape_l(p + 1);
+   Vector ddshape_x(p + 1), ddshape_y(p + 1), ddshape_z(p + 1), ddshape_l(p + 1);
+   DenseMatrix ddu(dof, ((dim + 1) * dim) / 2);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x, ddshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y, ddshape_y);
+   poly1d.CalcBasis(p, ip.z, shape_z, dshape_z, ddshape_z);
+   poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l, ddshape_l);
+
+   for (int o = 0, k = 0; k <= p; k++)
+      for (int j = 0; j + k <= p; j++)
+         for (int i = 0; i + j + k <= p; i++)
+         {
+            // u_xx, u_xy, u_xz, u_yy, u_yz, u_zz
+            int l = p - i - j - k;
+            ddu(o,0) = ((ddshape_x(i) * shape_l(l)) - 2. * (dshape_x(i) * dshape_l(l)) +
+                        (shape_x(i) * ddshape_l(l))) * shape_y(j) * shape_z(k);
+            ddu(o,1) = ((dshape_y(j) * ((dshape_x(i) * shape_l(l)) -
+                                        (shape_x(i) * dshape_l(l)))) +
+                        (shape_y(j) * ((ddshape_l(l) * shape_x(i)) -
+                                       (dshape_x(i) * dshape_l(l)))))* shape_z(k);
+            ddu(o,2) = ((dshape_z(k) * ((dshape_x(i) * shape_l(l)) -
+                                        (shape_x(i) * dshape_l(l)))) +
+                        (shape_z(k) * ((ddshape_l(l) * shape_x(i)) -
+                                       (dshape_x(i) * dshape_l(l)))))* shape_y(j);
+            ddu(o,3) = ((ddshape_y(j) * shape_l(l)) - 2. * (dshape_y(j) * dshape_l(l)) +
+                        (shape_y(j) * ddshape_l(l))) * shape_x(i) * shape_z(k);
+            ddu(o,4) = ((dshape_z(k) * ((dshape_y(j) * shape_l(l)) -
+                                        (shape_y(j)*dshape_l(l))) ) +
+                        (shape_z(k)* ((ddshape_l(l)*shape_y(j)) -
+                                      (dshape_y(j) * dshape_l(l)) ) ) )* shape_x(i);
+            ddu(o,5) = ((ddshape_z(k) * shape_l(l)) - 2. * (dshape_z(k) * dshape_l(l)) +
+                        (shape_z(k) * ddshape_l(l))) * shape_y(j) * shape_x(i);
+            o++;
+         }
+   Ti.Mult(ddu, ddshape);
+}
+
 void L2_TetrahedronElement::ProjectDelta(int vertex, Vector &dofs) const
 {
    switch (vertex)

fem/fe/fe_l2.hpp

@@ -22,7 +22,7 @@ class L2_SegmentElement : public NodalTensorFiniteElement
 {
 private:
 #ifndef MFEM_THREAD_SAFE
-   mutable Vector shape_x, dshape_x;
+   mutable Vector shape_x, dshape_x, d2shape_x;
 #endif
 
 public:
@@ -31,6 +31,8 @@ public:
    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &Hessian) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
 
    virtual void GetLocalRestriction(ElementTransformation &Trans,
@@ -45,7 +47,7 @@ class L2_QuadrilateralElement : public NodalTensorFiniteElement
 {
 private:
 #ifndef MFEM_THREAD_SAFE
-   mutable Vector shape_x, shape_y, dshape_x, dshape_y;
+   mutable Vector shape_x, shape_y, dshape_x, dshape_y, d2shape_x, d2shape_y;
 #endif
 
 public:
@@ -55,6 +57,8 @@ public:
    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &Hessian) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
    virtual void ProjectCurl(const FiniteElement &fe,
                             ElementTransformation &Trans,
@@ -79,7 +83,8 @@ class L2_HexahedronElement : public NodalTensorFiniteElement
 {
 private:
 #ifndef MFEM_THREAD_SAFE
-   mutable Vector shape_x, shape_y, shape_z, dshape_x, dshape_y, dshape_z;
+   mutable Vector shape_x, shape_y, shape_z, dshape_x, dshape_y, dshape_z,
+           d2shape_x, d2shape_y, d2shape_z;
 #endif
 
 public:
@@ -89,6 +94,8 @@ public:
    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &Hessian) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
 
    virtual void GetLocalRestriction(ElementTransformation &Trans,
@@ -110,7 +117,8 @@ class L2_TriangleElement : public NodalFiniteElement
 private:
 #ifndef MFEM_THREAD_SAFE
    mutable Vector shape_x, shape_y, shape_l, dshape_x, dshape_y, dshape_l, u;
-   mutable DenseMatrix du;
+   mutable Vector ddshape_x, ddshape_y, ddshape_l;
+   mutable DenseMatrix du, ddu;
 #endif
    DenseMatrixInverse Ti;
 
@@ -121,6 +129,8 @@ public:
    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &ddshape) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
    virtual void ProjectCurl(const FiniteElement &fe,
                             ElementTransformation &Trans,
@@ -141,7 +151,8 @@ private:
 #ifndef MFEM_THREAD_SAFE
    mutable Vector shape_x, shape_y, shape_z, shape_l;
    mutable Vector dshape_x, dshape_y, dshape_z, dshape_l, u;
-   mutable DenseMatrix du;
+   mutable Vector ddshape_x, ddshape_y, ddshape_z, ddshape_l;
+   mutable DenseMatrix du, ddu;
 #endif
    DenseMatrixInverse Ti;
 
@@ -152,6 +163,8 @@ public:
    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &ddshape) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
 
    virtual void GetLocalRestriction(ElementTransformation &Trans,

fem/hyperbolic.cpp

@@ -189,7 +189,7 @@ HyperbolicFormIntegrator::HyperbolicFormIntegrator(
 
 real_t FluxFunction::ComputeFluxDotN(const Vector &U,
                                      const Vector &normal,
-                                     FaceElementTransformations &Tr,
+                                     ElementTransformation &Tr,
                                      Vector &FUdotN) const
 {
 #ifdef MFEM_THREAD_SAFE
@@ -208,8 +208,8 @@ real_t RusanovFlux::Eval(const Vector &state1, const Vector &state2,
 #ifdef MFEM_THREAD_SAFE
    Vector fluxN1(fluxFunction.num_equations), fluxN2(fluxFunction.num_equations);
 #endif
-   const real_t speed1 = fluxFunction.ComputeFluxDotN(state1, nor, Tr, fluxN1);
-   const real_t speed2 = fluxFunction.ComputeFluxDotN(state2, nor, Tr, fluxN2);
+   const real_t speed1 = fluxFunction.ComputeFluxDotN(state1, nor, *Tr.Elem1, fluxN1);
+   const real_t speed2 = fluxFunction.ComputeFluxDotN(state2, nor, *Tr.Elem2, fluxN2);
    // NOTE: nor in general is not a unit normal
    const real_t maxE = std::max(speed1, speed2);
    // here, std::sqrt(nor*nor) is multiplied to match the scale with fluxN
@@ -253,7 +253,7 @@ real_t ShallowWaterFlux::ComputeFlux(const Vector &U,
 
    const real_t energy = 0.5 * g * (height * height);
 
-   MFEM_ASSERT(height >= 0, "Negative Height");
+   MFEM_ASSERT(height >= 0, "Negative Height: " << height);
 
    for (int d = 0; d < dim; d++)
    {
@@ -274,7 +274,7 @@ real_t ShallowWaterFlux::ComputeFlux(const Vector &U,
 
 real_t ShallowWaterFlux::ComputeFluxDotN(const Vector &U,
                                          const Vector &normal,
-                                         FaceElementTransformations &Tr,
+                                         ElementTransformation &Tr,
                                          Vector &FUdotN) const
 {
    const real_t height = U(0);
@@ -282,7 +282,7 @@ real_t ShallowWaterFlux::ComputeFluxDotN(const Vector &U,
 
    const real_t energy = 0.5 * g * (height * height);
 
-   MFEM_ASSERT(height >= 0, "Negative Height");
+   MFEM_ASSERT(height >= 0, "Negative Height: " << height);
    FUdotN(0) = h_vel * normal;
    const real_t normal_vel = FUdotN(0) / height;
    for (int i = 0; i < dim; i++)
@@ -348,7 +348,7 @@ real_t EulerFlux::ComputeFlux(const Vector &U,
 
 real_t EulerFlux::ComputeFluxDotN(const Vector &x,
                                   const Vector &normal,
-                                  FaceElementTransformations &Tr,
+                                  ElementTransformation &Tr,
                                   Vector &FUdotN) const
 {
    // 1. Get states

fem/hyperbolic.hpp

@@ -92,7 +92,7 @@ public:
     * @return real_t maximum (normal) characteristic velocity
     */
    virtual real_t ComputeFluxDotN(const Vector &state, const Vector &normal,
-                                  FaceElementTransformations &Tr,
+                                  ElementTransformation &Tr,
                                   Vector &fluxDotN) const;
 
    /**
@@ -375,7 +375,7 @@ public:
     * @return real_t maximum characteristic speed, |u| + √(γp/ρ)
     */
    real_t ComputeFluxDotN(const Vector &state, const Vector &normal,
-                          FaceElementTransformations &Tr,
+                          ElementTransformation &Tr,
                           Vector &fluxN) const override;
 };
 
@@ -418,7 +418,7 @@ public:
     * @return real_t maximum characteristic speed, |u| + √(γp/ρ)
     */
    real_t ComputeFluxDotN(const Vector &x, const Vector &normal,
-                          FaceElementTransformations &Tr,
+                          ElementTransformation &Tr,
                           Vector &fluxN) const override;
 };
 

fem/lininteg.cpp

@@ -274,6 +274,7 @@ void VectorDomainLFIntegrator::AssembleRHSElementVect(
    real_t val,cf;
 
    shape.SetSize(dof);       // vector of size dof
+   Qvec.SetSize(vdim);
 
    elvect.SetSize(dof * vdim);
    elvect = 0.0;

linalg/densemat.cpp

@@ -92,6 +92,15 @@ DenseMatrix::DenseMatrix(const DenseMatrix &mat, char ch)
    }
 }
 
+void DenseMatrix::Resize(int h, int w)
+{
+   MFEM_ASSERT(h*w == height*width,
+               "New size (" << h << " x " << w << ") is incompatible with the original size ("
+               << height << " x " << width << ").");
+   height = h;
+   width = w;
+}
+
 void DenseMatrix::SetSize(int h, int w)
 {
    MFEM_ASSERT(h >= 0 && w >= 0,
@@ -1283,8 +1292,8 @@ int DenseMatrix::Rank(real_t tol) const
 
 real_t DenseMatrix::CalcSingularvalue(const int i) const
 {
-   MFEM_ASSERT(Height() == Width() && Height() > 0 && Height() < 4,
-               "The matrix must be square and sized 1, 2, or 3 to compute the"
+   MFEM_ASSERT(Height() >= Width() && Height() > 0 && Height() < 4,
+               "The matrix must be tall and sized 1, 2, or 3 to compute the"
                " singular values."
                << "  Height() = " << Height()
                << ", Width() = " << Width());
@@ -1292,18 +1301,28 @@ real_t DenseMatrix::CalcSingularvalue(const int i) const
    const int n = Height();
    const real_t *d = data;
 
-   if (n == 1)
+   if (Height() == Width())
    {
-      return d[0];
-   }
-   else if (n == 2)
-   {
-      return kernels::CalcSingularvalue<2>(d,i);
+      if (n == 1)
+      {
+         return d[0];
+      }
+      else if (n == 2)
+      {
+         return kernels::CalcSingularvalue<2>(d,i);
+      }
+      else
+      {
+         return kernels::CalcSingularvalue<3>(d,i);
+      }
    }
    else
    {
-      return kernels::CalcSingularvalue<3>(d,i);
+      DenseMatrix AtA(Width());
+      MultAtB(*this, *this, AtA);
+      return std::sqrt(AtA.CalcSingularvalue(i));
    }
+
 }
 
 void DenseMatrix::CalcEigenvalues(real_t *lambda, real_t *vec) const
@@ -1443,9 +1462,7 @@ void DenseMatrix::Transpose()
       for (i = 0; i < Height(); i++)
          for (j = i+1; j < Width(); j++)
          {
-            t = (*this)(i,j);
-            (*this)(i,j) = (*this)(j,i);
-            (*this)(j,i) = t;
+            std::swap((*this)(j,i), (*this)(i,j));
          }
    }
    else
@@ -3292,6 +3309,8 @@ void AddMult_a_VWt(const real_t a, const Vector &v, const Vector &w,
 #ifdef MFEM_DEBUG
    if (VWt.Height() != m || VWt.Width() != n)
    {
+      out << "Expected VWt size: " << m << " x " << n << std::endl;  
+      out << "Provided VWt size: " << VWt.Height() << " x " << VWt.Width() << std::endl;  
       mfem_error("AddMult_a_VWt(...): dimension mismatch");
    }
 #endif

linalg/densemat.hpp

@@ -104,6 +104,9 @@ public:
    /// Change the size of the DenseMatrix to s x s.
    void SetSize(int s) { SetSize(s, s); }
 
+   /// Change the size while keeping the data
+   void Resize(int h, int w);
+
    /// Change the size of the DenseMatrix to h x w.
    void SetSize(int h, int w);
 

makefile

@@ -125,7 +125,7 @@ EXAMPLE_TEST_DIRS := examples
 
 MINIAPP_SUBDIRS = common electromagnetics meshing navier performance tools \
  toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff hooke \
- multidomain dpg hdiv-linear-solver spde
+ multidomain dpg hdiv-linear-solver spde mani-swe
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \

miniapps/CMakeLists.txt

@@ -37,3 +37,4 @@ add_subdirectory(tribol)
 add_subdirectory(hooke)
 add_subdirectory(dpg)
 add_subdirectory(hdiv-linear-solver)
+add_subdirectory(mani-swe)

miniapps/mani-swe/CMakeLists.txt

@@ -0,0 +1,34 @@
+# Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+if (MFEM_USE_MPI)
+  add_mfem_miniapp(test_sphere
+    MAIN test_sphere.cpp
+    EXTRA_SOURCES manihyp.cpp
+    EXTRA_HEADERS manihyp.hpp
+    LIBRARIES mfem)
+  add_mfem_miniapp(test_square
+    MAIN test_square.cpp
+    EXTRA_SOURCES manihyp.cpp
+    EXTRA_HEADERS manihyp.hpp
+    LIBRARIES mfem)
+  add_mfem_miniapp(test_maniflux
+    MAIN test_maniflux.cpp
+    EXTRA_SOURCES manihyp.cpp
+    EXTRA_HEADERS manihyp.hpp
+    LIBRARIES mfem)
+  add_mfem_miniapp(test_integrator
+    MAIN test_integrator.cpp
+    EXTRA_SOURCES manihyp.cpp
+    EXTRA_HEADERS manihyp.hpp
+    LIBRARIES mfem)
+endif()
+file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/data DESTINATION ${CMAKE_CURRENT_BINARY_DIR})

miniapps/mani-swe/data/icosahedron.mesh

@@ -0,0 +1,119 @@
+MFEM mesh v1.0
+
+dimension
+2
+
+elements
+20
+   1    2    0    4    1 
+   1    2    1    7    0 
+   1    2    4    0    8 
+   1    2    7    9    0 
+   1    2    0    9    8 
+   1    2    4   11    1 
+   1    2    7    1   10 
+   1    2    1   11   10 
+   1    2    3    2    5 
+   1    2    6    2    3 
+   1    2    5    2   11 
+   1    2    2    6   10 
+   1    2   11    2   10 
+   1    2    3    5    8 
+   1    2    6    3    9 
+   1    2    9    3    8 
+   1    2    8    5    4 
+   1    2   11    4    5 
+   1    2    6    9    7 
+   1    2    6    7   10 
+
+boundary
+0
+
+vertices
+13
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: L2_T1_2D_P1
+VDim: 3
+Ordering: 1
+
+0.850651 0.525731 0.000000 
+0.525731 0.000000 0.850651 
+0.850651 -0.525731 0.000000 
+
+0.850651 -0.525731 0.000000 
+0.525731 0.000000 -0.850651 
+0.850651 0.525731 0.000000 
+
+0.525731 0.000000 0.850651 
+0.850651 0.525731 0.000000 
+0.000000 0.850651 0.525731 
+
+0.525731 0.000000 -0.850651 
+0.000000 0.850651 -0.525731 
+0.850651 0.525731 0.000000 
+
+0.850651 0.525731 0.000000 
+0.000000 0.850651 -0.525731 
+0.000000 0.850651 0.525731 
+
+0.525731 0.000000 0.850651 
+0.000000 -0.850651 0.525731 
+0.850651 -0.525731 0.000000 
+
+0.525731 0.000000 -0.850651 
+0.850651 -0.525731 0.000000 
+0.000000 -0.850651 -0.525731 
+
+0.850651 -0.525731 0.000000 
+0.000000 -0.850651 0.525731 
+0.000000 -0.850651 -0.525731 
+
+-0.850651 0.525731 0.000000 
+-0.850651 -0.525731 0.000000 
+-0.525731 0.000000 0.850651 
+
+-0.525731 0.000000 -0.850651 
+-0.850651 -0.525731 0.000000 
+-0.850651 0.525731 0.000000 
+
+-0.525731 0.000000 0.850651 
+-0.850651 -0.525731 0.000000 
+0.000000 -0.850651 0.525731 
+
+-0.850651 -0.525731 0.000000 
+-0.525731 0.000000 -0.850651 
+0.000000 -0.850651 -0.525731 
+
+0.000000 -0.850651 0.525731 
+-0.850651 -0.525731 0.000000 
+0.000000 -0.850651 -0.525731 
+
+-0.850651 0.525731 0.000000 
+-0.525731 0.000000 0.850651 
+0.000000 0.850651 0.525731 
+
+-0.525731 0.000000 -0.850651 
+-0.850651 0.525731 0.000000 
+0.000000 0.850651 -0.525731 
+
+0.000000 0.850651 -0.525731 
+-0.850651 0.525731 0.000000 
+0.000000 0.850651 0.525731 
+
+0.000000 0.850651 0.525731 
+-0.525731 0.000000 0.850651 
+0.525731 0.000000 0.850651 
+
+0.000000 -0.850651 0.525731 
+0.525731 0.000000 0.850651 
+-0.525731 0.000000 0.850651 
+
+-0.525731 0.000000 -0.850651 
+0.000000 0.850651 -0.525731 
+0.525731 0.000000 -0.850651 
+
+-0.525731 0.000000 -0.850651 
+0.525731 0.000000 -0.850651 
+0.000000 -0.850651 -0.525731 

miniapps/mani-swe/data/periodic-square-3d.mesh

@@ -0,0 +1,127 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+2
+
+# format: <attribute> <geometry type> <vertex 0> <vertex 1> ...
+elements
+16
+1 3  0  1  5  4
+1 3  1  2  6  5
+1 3  2  3  7  6
+1 3  3  0  4  7
+1 3  4  5  9  8
+1 3  5  6 10  9
+1 3  6  7 11 10
+1 3  7  4  8 11
+1 3  8  9 13 12
+1 3  9 10 14 13
+1 3 10 11 15 14
+1 3 11  8 12 15
+1 3 12 13  1  0
+1 3 13 14  2  1
+1 3 14 15  3  2
+1 3 15 12  0  3
+
+boundary
+0
+
+vertices
+16
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: L2_T1_2D_P1
+VDim: 3
+Ordering: 1
+
+-1.0 -1.0 1
+-0.5 -1.0 1
+-1.0 -0.5 1
+-0.5 -0.5 1
+
+-0.5 -1.0 1
+   0 -1.0 1
+-0.5 -0.5 1
+   0 -0.5 1
+
+   0 -1.0 1
+ 0.5 -1.0 1
+   0 -0.5 1
+ 0.5 -0.5 1
+
+ 0.5 -1.0 1
+ 1.0 -1.0 1
+ 0.5 -0.5 1
+ 1.0 -0.5 1
+
+-1.0 -0.5 1
+-0.5 -0.5 1
+-1.0    0 1
+-0.5    0 1
+
+-0.5 -0.5 1
+   0 -0.5 1
+-0.5    0 1
+   0    0 1
+
+   0 -0.5 1
+ 0.5 -0.5 1
+   0    0 1
+ 0.5    0 1
+
+ 0.5 -0.5 1
+ 1.0 -0.5 1
+ 0.5    0 1
+ 1.0    0 1
+
+-1.0    0 1
+-0.5    0 1
+-1.0  0.5 1
+-0.5  0.5 1
+
+-0.5    0 1
+   0    0 1
+-0.5  0.5 1
+   0  0.5 1
+
+   0    0 1
+ 0.5    0 1
+   0  0.5 1
+ 0.5  0.5 1
+
+ 0.5    0 1
+ 1.0    0 1
+ 0.5  0.5 1
+ 1.0  0.5 1
+
+-1.0  0.5 1
+-0.5  0.5 1
+-1.0  1.0 1
+-0.5  1.0 1
+
+-0.5  0.5 1
+   0  0.5 1
+-0.5  1.0 1
+   0  1.0 1
+
+   0  0.5 1
+ 0.5  0.5 1
+   0  1.0 1
+ 0.5  1.0 1
+
+ 0.5  0.5 1
+ 1.0  0.5 1
+ 0.5  1.0 1
+ 1.0  1.0 1

miniapps/mani-swe/manihyp.cpp

@@ -0,0 +1,517 @@
+#include "manihyp.hpp"
+
+namespace mfem
+{
+
+void sphere(const Vector &x, Vector &y, const real_t r)
+{
+   y = x; y *= r/std::sqrt(y*y);
+}
+
+void CalcOrtho(const DenseMatrix &faceJ, const DenseMatrix &elemJ, Vector &n)
+{
+   const int sdim = faceJ.Height();
+   const int dim = elemJ.Width();
+   MFEM_ASSERT(sdim == 3 && dim == 2, "Only supports 2D manifold in 3D");
+   MFEM_ASSERT(faceJ.Width() == 1, "FaceJ is not a vector");
+
+   Vector tangent(faceJ.GetData(), sdim);
+   Vector normal1(elemJ.GetData(), sdim);
+   Vector normal2(elemJ.GetData() + sdim, sdim);
+
+   Vector surfaceNormal(sdim);
+   normal1.cross3D(normal2, surfaceNormal);
+
+   tangent.cross3D(surfaceNormal, n);
+   n *= std::sqrt((tangent*tangent)/(n*n));
+}
+
+void ManifoldCoord::convertElemState(ElementTransformation &Tr,
+                                     const int nrScalar, const int nrVector,
+                                     const Vector &state, Vector &phys_state) const
+{
+   for (int i=0; i<nrScalar; i++)
+   {
+      phys_state[i] = state[i];
+   }
+   const DenseMatrix &J = Tr.Jacobian();
+   mani_vec_state.UseExternalData(state.GetData()+nrScalar, dim, nrVector);
+   phys_vec_state.UseExternalData(phys_state.GetData()+nrScalar, sdim, nrVector);
+   Mult(J, mani_vec_state, phys_vec_state);
+}
+
+void ManifoldCoord::convertFaceState(FaceElementTransformations &Tr,
+                                     const int nrScalar, const int nrVector,
+                                     const Vector &stateL, const Vector &stateR,
+                                     Vector &normalL, Vector &normalR,
+                                     Vector &stateL_L, Vector &stateR_L,
+                                     Vector &stateL_R, Vector &stateR_R) const
+{
+   // face Jacobian
+   const DenseMatrix fJ = Tr.Jacobian();
+
+   // element Jacobians
+   const DenseMatrix &J1 = Tr.Elem1->Jacobian();
+   const DenseMatrix &J2 = Tr.Elem2 ? Tr.Elem2->Jacobian() : Tr.Elem1->Jacobian();
+
+   normal_comp.SetSize(nrVector);
+
+   // Compute interface normal vectors at each element
+   CalcOrtho(fJ, J1, normalL);
+   CalcOrtho(fJ, J2, normalR);
+   const real_t tangent_norm = Tr.Weight();
+
+   // copy scalar states
+   for (int i=0; i<nrScalar; i++)
+   {
+      stateL_R[i] = stateL[i];
+      stateR_L[i] = stateR[i];
+   }
+
+   // Convert Left element vector states to physical states
+   mani_vec_state.UseExternalData(stateL.GetData() + nrScalar, dim, nrVector);
+   phys_vec_state.UseExternalData(stateL_R.GetData() + nrScalar, sdim, nrVector);
+   Mult(J1, mani_vec_state, phys_vec_state);
+   stateL_L = stateL_R; // Left to Left done
+   for (int i=0; i<nrVector; i++)
+   {
+      phys_vec_state.GetColumnReference(i, phys_vec);
+      const real_t normal_comp = phys_vec*normalL/(tangent_norm*tangent_norm);
+      phys_vec.Add(-normal_comp, normalL).Add(normal_comp, normalR);
+   }
+
+   // Convert Right element vector states to physical states
+   mani_vec_state.UseExternalData(stateR.GetData() + nrScalar, dim, nrVector);
+   phys_vec_state.UseExternalData(stateR_L.GetData() + nrScalar, sdim, nrVector);
+   Mult(J2, mani_vec_state, phys_vec_state);
+   stateR_R = stateR_L;
+   for (int i=0; i<nrVector; i++)
+   {
+      phys_vec_state.GetColumnReference(i, phys_vec);
+      const real_t normal_comp = phys_vec*normalR/(tangent_norm*tangent_norm);
+      phys_vec.Add(-normal_comp, normalR).Add(normal_comp, normalL);
+   }
+}
+
+const IntegrationRule &ManifoldVectorMassIntegrator::GetRule(
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &Trans)
+{
+   return IntRules.Get(Trans.GetGeometryType(),
+                       trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderJ()*2 + Trans.OrderW());
+}
+
+void ManifoldVectorMassIntegrator::AssembleElementMatrix(
+   const FiniteElement &el, ElementTransformation &Trans,
+   DenseMatrix &elmat )
+{
+   int dof = el.GetDof();
+   dim = el.GetDim();
+   sdim = Trans.GetSpaceDim();
+   real_t w;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape;
+#endif
+   elmat.SetSize(dof*dim);
+   elmat_comp.SetSize(dof);
+   elmat_comp_weighted.SetSize(dof);
+   JtJ.SetSize(dim);
+   shape.SetSize(dof);
+
+   const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el, Trans);
+
+   elmat = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Trans.SetIntPoint (&ip);
+
+      el.CalcPhysShape(Trans, shape);
+      const DenseMatrix &J = Trans.Jacobian();
+      MultAtB(J, J, JtJ);
+
+      w = Trans.Weight() * ip.weight;
+      MultVVt(shape, elmat_comp);
+      for (int col=0; col<dim; col++)
+      {
+         for (int row=0; row<dim; row++)
+         {
+            elmat_comp_weighted = elmat_comp;
+            elmat_comp_weighted *= w*JtJ(row, col);
+            elmat.AddSubMatrix(dof*row, dof*col, elmat_comp_weighted);
+         }
+      }
+   }
+}
+
+void ManifoldVectorGradientIntegrator::AssembleElementMatrix(
+   const FiniteElement &el,
+   ElementTransformation &Trans,
+   DenseMatrix &elmat)
+{
+}
+
+const IntegrationRule &ManifoldVectorGradientIntegrator::GetRule(
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &Trans)
+{
+   return IntRules.Get(Trans.GetGeometryType(),
+                       trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderJ()*2 + Trans.OrderW());
+}
+
+
+
+real_t ManifoldFlux::ComputeFlux(const Vector &state, ElementTransformation &Tr,
+                                 DenseMatrix &flux) const
+{
+   phys_state.SetSize(nrScalar + coord.sdim*nrVector);
+   coord.convertElemState(Tr, nrScalar, nrVector, state, phys_state);
+   return org_flux.ComputeFlux(phys_state, Tr, flux);
+}
+
+real_t ManifoldFlux::ComputeNormalFluxes(const Vector &stateL,
+                                         const Vector &stateR,
+                                         FaceElementTransformations &Tr,
+                                         Vector &normalL, Vector &normalR,
+                                         Vector &stateL_L, Vector &stateR_L,
+                                         Vector &fluxL_L, Vector &fluxR_L,
+                                         Vector &stateL_R, Vector &stateR_R,
+                                         Vector &fluxL_R, Vector &fluxR_R) const
+{
+   coord.convertFaceState(Tr, nrScalar, nrVector,
+                          stateL, stateR,
+                          normalL, normalR,
+                          stateL_L, stateR_L,
+                          stateL_R, stateR_R);
+   ElementTransformation *Tr1 = Tr.Elem1;
+   ElementTransformation *Tr2 = Tr.Elem2 ? Tr.Elem2 : Tr.Elem1;
+
+   real_t mcs = org_flux.ComputeFluxDotN(
+                   stateL_L, normalL, *Tr1, fluxL_L);
+   mcs = std::max(mcs, org_flux.ComputeFluxDotN(
+                     stateR_L, normalL, *Tr1, fluxR_L));
+   mcs = std::max(mcs, org_flux.ComputeFluxDotN(
+                     stateL_R, normalR, *Tr2, fluxL_R));
+   mcs = std::max(mcs, org_flux.ComputeFluxDotN(
+                     stateR_R, normalR, *Tr2, fluxR_R));
+   return mcs;
+}
+
+ManifoldHyperbolicFormIntegrator::ManifoldHyperbolicFormIntegrator(
+   const ManifoldNumericalFlux &flux, const IntegrationRule *ir)
+   :numFlux(flux), maniFlux(flux.GetManifoldFluxFunction()),
+    coord(maniFlux.GetCoordinate()), intrule(ir), dg_fec(0, coord.dim)
+{
+   switch (maniFlux.GetCoordinate().dim)
+   {
+      case 2:
+      {
+         hess_map.SetSize(4);
+         hess_map[0] = 0;
+         hess_map[1] = 1;
+         hess_map[2] = 1;
+         hess_map[3] = 2;
+         break;
+      }
+      default:
+      {
+         MFEM_ABORT("Only support 2D manifold");
+      }
+   }
+}
+
+
+void ManifoldHyperbolicFormIntegrator::AssembleElementVector(
+   const FiniteElement &el, ElementTransformation &Tr,
+   const Vector &elfun, Vector &elvect)
+{
+   // element info
+   const int dof = el.GetDof();
+   const int dim = el.GetDim();
+   const int sdim = Tr.GetSpaceDim();
+   const int nrScalar = maniFlux.GetNumScalars();
+   const int nrVector = (maniFlux.num_equations - nrScalar)/sdim;
+   const int vdim = elfun.Size()/dof;
+   MFEM_ASSERT((vdim - nrScalar) / dim == nrVector,
+               "The number of equations and vector dimension disagree");
+
+   elvect.SetSize(vdim*dof);
+   elvect = 0.0;
+
+   shape.SetSize(dof);
+   dshape.SetSize(dof, dim);
+   gshape.SetSize(dof, sdim);
+   vector_gshape.SetSize(dof*dim, sdim*sdim);
+   vector_gshape_comp.SetSize(dof, sdim*sdim);
+   state.SetSize(vdim);
+   phys_flux.SetSize(maniFlux.num_equations, sdim);
+   phys_flux_scalars.SetSize(nrScalar, sdim);
+   phys_flux_vectors.SetSize(nrVector*sdim, sdim);
+   gradJ.SetSize(sdim, sdim);
+   Vector gradJ_vectorview(gradJ.GetData(), sdim*sdim);
+   HessMat.SetSize(sdim, dim, dim);
+   Vector gshape_vectorview(gshape.GetData(), dof*sdim);
+   Vector J_comp;
+
+   const DenseMatrix u(elfun.GetData(), dof, vdim);
+   const DenseMatrix u_scalars(elfun.GetData(), dof, nrScalar);
+   const DenseMatrix u_vectors(elfun.GetData() + nrScalar*dof, dim*dof, nrVector);
+
+   DenseMatrix divflux_scalars(elvect.GetData(), dof, nrScalar);
+   DenseMatrix divflux_vectors(elvect.GetData() + nrScalar*dof, dim*dof, nrVector);
+
+   const IntegrationRule &ir = intrule ? *intrule : GetRule(el, el, Tr);
+   for (int i=0; i<ir.GetNPoints(); i++)
+   {
+      const IntegrationPoint & ip = ir.IntPoint(i);
+      Tr.SetIntPoint(&ip);
+
+      const DenseMatrix &J = Tr.Jacobian();
+      const DenseMatrix &invJ = Tr.InverseJacobian();
+      const DenseMatrix &Hess = Tr.Hessian();
+      for (int d1 = 0; d1 < dim; d1++)
+         for (int d2 = 0; d2 < dim; d2++)
+            for (int sd = 0; sd < sdim; sd++)
+            {
+               HessMat(sd, d1, d2) = Hess(sd, hess_map[d1*dim + d2]);
+            }
+
+      el.CalcShape(ip, shape);
+      el.CalcPhysDShape(Tr, gshape);
+
+      u.MultTranspose(shape, state);
+
+      max_char_speed = std::max(max_char_speed, maniFlux.ComputeFlux(state, Tr,
+                                                                     phys_flux));
+      phys_flux.GetSubMatrix(0, nrScalar, 0, sdim, phys_flux_scalars);
+      AddMult_a_ABt(ip.weight*Tr.Weight(), gshape, phys_flux_scalars,
+                    divflux_scalars);
+
+      // prepare physical vector flux for integration
+      phys_flux.GetSubMatrix(nrScalar, maniFlux.num_equations, 0, sdim,
+                             phys_flux_vectors);
+      phys_flux_vectors.Transpose();
+      phys_flux_vectors.Resize(sdim*sdim, nrVector);
+
+      for (int d = 0; d < dim; d++)
+      {
+         Mult(HessMat(d), invJ, gradJ);
+         MultVWt(shape, gradJ_vectorview, vector_gshape_comp);
+         J_comp.SetDataAndSize(J.GetData() + sdim*d, sdim);
+         vector_gshape_comp.Resize(sdim*dof, sdim);
+         AddMult_a_VWt(1.0, gshape_vectorview, J_comp, vector_gshape_comp);
+         vector_gshape_comp.Resize(dof, sdim*sdim);
+         vector_gshape.SetSubMatrix(dof*d, 0, vector_gshape_comp);
+      }
+      AddMult_a(ip.weight*Tr.Weight(), vector_gshape, phys_flux_vectors,
+                divflux_vectors);
+      phys_flux_vectors.Resize(sdim*nrVector, sdim);
+   }
+}
+
+void ManifoldHyperbolicFormIntegrator::AssembleFaceVector(
+   const FiniteElement &el1, const FiniteElement &el2,
+   FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect)
+{
+   // element info
+   const int dof1 = el1.GetDof();
+   const int dof2 = el2.GetDof();
+   const int dim = el1.GetDim();
+   const int sdim = Tr.GetSpaceDim();
+   const int nrScalar = maniFlux.GetNumScalars();
+   const int nrVector = (maniFlux.num_equations - nrScalar)/sdim;
+   const int vdim = elfun.Size()/(dof1 + dof2);
+
+   elvect.SetSize((dof1+dof2)*vdim);
+   elvect = 0.0;
+   shape1.SetSize(dof1);
+   shape2.SetSize(dof2);
+   stateL.SetSize(vdim);
+   stateR.SetSize(vdim);
+   phys_hatFL.SetSize(maniFlux.num_equations);
+   phys_hatFR.SetSize(maniFlux.num_equations);
+   hatFL.SetSize(vdim);
+   hatFR.SetSize(vdim);
+   const DenseMatrix phys_hatFL_vectors(phys_hatFL.GetData() + nrScalar, sdim,
+                                        nrVector);
+   const DenseMatrix phys_hatFR_vectors(phys_hatFR.GetData() + nrScalar, sdim,
+                                        nrVector);
+   DenseMatrix hatFL_vectors(hatFL.GetData() + nrScalar, dim, nrVector);
+   DenseMatrix hatFR_vectors(hatFR.GetData() + nrScalar, dim, nrVector);
+   const DenseMatrix u1(elfun.GetData(), dof1, vdim);
+   const DenseMatrix u2(elfun.GetData() + dof1*vdim, dof2, vdim);
+   DenseMatrix jumpflux1(elvect.GetData(), dof1, vdim);
+   DenseMatrix jumpflux2(elvect.GetData()+dof1*vdim, dof2, vdim);
+   const IntegrationRule &ir = intrule ? *intrule : GetRule(el1, el2, Tr);
+   for (int i=0; i < ir.GetNPoints(); i++)
+   {
+      const IntegrationPoint & ip = ir.IntPoint(i);
+      Tr.SetAllIntPoints(&ip);
+      el1.CalcShape(Tr.GetElement1IntPoint(), shape1);
+      el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
+
+      u1.MultTranspose(shape1, stateL);
+      u2.MultTranspose(shape2, stateR);
+
+      const DenseMatrix &J1 = Tr.Elem1->Jacobian();
+      const DenseMatrix &J2 = Tr.Elem2->Jacobian();
+      max_char_speed = std::max(max_char_speed, numFlux.Eval(stateL, stateR, Tr,
+                                                             phys_hatFL, phys_hatFR));
+      for (int j=0; j<nrScalar; j++) { hatFL[j] = phys_hatFL[j]; hatFR[j] = phys_hatFR[j]; }
+      MultAtB(J1, phys_hatFL_vectors, hatFL_vectors);
+      MultAtB(J2, phys_hatFR_vectors, hatFR_vectors);
+      AddMult_a_VWt(-ip.weight, shape1, hatFL, jumpflux1);
+      AddMult_a_VWt(+ip.weight, shape2, hatFR, jumpflux2);
+   }
+}
+
+ManifoldDGHyperbolicConservationLaws::ManifoldDGHyperbolicConservationLaws(
+   FiniteElementSpace &vfes,
+   ManifoldHyperbolicFormIntegrator &formIntegrator,
+   const int nrScalar, const int int_offset)
+   : TimeDependentOperator(vfes.GetTrueVSize()),
+     vfes(vfes),
+     dim(vfes.GetMesh()->Dimension()),
+     sdim(vfes.GetMesh()->SpaceDimension()),
+     nrScalar(nrScalar),
+     nrVector((vfes.GetVDim()-nrScalar)/dim),
+     formIntegrator(formIntegrator),
+     z(vfes.GetTrueVSize()),
+     int_offset(int_offset)
+{
+   ComputeInvMass();
+   // ComputeWeakDivergence();
+#ifdef MFEM_USE_MPI
+   ParFiniteElementSpace *pvfes = dynamic_cast<ParFiniteElementSpace *>(&vfes);
+   if (pvfes)
+   {
+      parallel = true;
+      comm = pvfes->GetComm();
+   }
+#endif
+
+   if (parallel)
+   {
+#ifdef MFEM_USE_MPI
+      nonlinearForm.reset(new ParNonlinearForm(pvfes));
+#endif
+   }
+   else
+   {
+      nonlinearForm.reset(new NonlinearForm(&vfes));
+   }
+   nonlinearForm->AddDomainIntegrator(&formIntegrator);
+   nonlinearForm->AddInteriorFaceIntegrator(&formIntegrator);
+   nonlinearForm->UseExternalIntegrators();
+}
+
+void ManifoldDGHyperbolicConservationLaws::ComputeInvMass()
+{
+   InverseIntegrator inv_mass(new MassIntegrator());
+   InverseIntegrator inv_vec_mass(new ManifoldVectorMassIntegrator());
+
+   invmass.resize(vfes.GetNE());
+   invmass_vec.resize(vfes.GetNE());
+   for (int i=0; i<vfes.GetNE(); i++)
+   {
+      int dof = vfes.GetFE(i)->GetDof();
+      invmass[i].SetSize(dof);
+      inv_mass.AssembleElementMatrix(*vfes.GetFE(i),
+                                     *vfes.GetElementTransformation(i),
+                                     invmass[i]);
+      invmass_vec[i].SetSize(dim*dof);
+      inv_vec_mass.AssembleElementMatrix(*vfes.GetFE(i),
+                                         *vfes.GetElementTransformation(i),
+                                         invmass_vec[i]);
+   }
+}
+
+void ManifoldDGHyperbolicConservationLaws::ComputeWeakDivergence()
+{
+   TransposeIntegrator div(new GradientIntegrator());
+   TransposeIntegrator div_vec(new ManifoldVectorGradientIntegrator());
+
+   weakdiv.resize(vfes.GetNE());
+   weakdiv_vec.resize(vfes.GetNE());
+   DG_FECollection dg_fec(0, dim, BasisType::GaussLegendre);
+   for (int i=0; i<vfes.GetNE(); i++)
+   {
+      int dof = vfes.GetFE(i)->GetDof();
+      const FiniteElement * fe = dg_fec.GetFE(vfes.GetElementTransformation(
+                                                 i)->GetGeometryType(),
+                                              vfes.GetOrder(i) + int_offset);
+
+      weakdiv[i].SetSize(dof, sdim*fe->GetNodes().GetNPoints());
+      div.SetIntegrationRule(fe->GetNodes());
+      div.AssembleElementMatrix2(*vfes.GetFE(i), *fe,
+                                 *vfes.GetElementTransformation(i),
+                                 weakdiv[i]);
+
+      weakdiv_vec[i].SetSize(dim*dof, sdim*sdim*fe->GetNodes().GetNPoints());
+      div_vec.SetIntegrationRule(fe->GetNodes());
+      div_vec.AssembleElementMatrix2(*vfes.GetFE(i), *fe,
+                                     *vfes.GetElementTransformation(i),
+                                     weakdiv_vec[i]);
+   }
+}
+
+void ManifoldDGHyperbolicConservationLaws::Mult(const Vector &x,
+                                                Vector &y) const
+{
+   // 0. Reset wavespeed computation before operator application.
+   formIntegrator.ResetMaxCharSpeed();
+   // 1. Apply Nonlinear form to obtain an auxiliary result
+   //         z = - <F̂(u_h,n), [[v]]>_e
+   //    If weak-divergence is not preassembled, we also have weak-divergence
+   //         z = - <F̂(u_h,n), [[v]]>_e + (F(u_h), ∇v)
+   if (force)
+   {
+      force->Assemble();
+      nonlinearForm->AddMult(x,z);
+   }
+   else
+   {
+      nonlinearForm->Mult(x, z);
+   }
+   // Apply block inverse mass
+   Vector zval; // z_loc, dof*num_eq
+
+   DenseMatrix current_zmat; // view of element auxiliary result, dof x num_eq
+   DenseMatrix current_ymat; // view of element result, dof x num_eq
+   Array<int> vdofs;
+   Array<int> vdofs_scalars;
+   Array<int> vdofs_vectors;
+   for (int i=0; i<vfes.GetNE(); i++)
+   {
+      int dof = vfes.GetFE(i)->GetDof();
+      vfes.GetElementVDofs(i, vdofs);
+
+      // Scalar mass inversion
+      vdofs_scalars.MakeRef(vdofs.GetData(), nrScalar*dof, false);
+      z.GetSubVector(vdofs_scalars, zval);
+      current_zmat.UseExternalData(zval.GetData(), dof, nrScalar);
+      current_ymat.SetSize(dof, nrScalar);
+      mfem::Mult(invmass[i], current_zmat, current_ymat);
+      y.SetSubVector(vdofs_scalars, current_ymat.GetData());
+
+      // Vector mass inversion
+      vdofs_vectors.MakeRef(vdofs.GetData() + nrScalar*dof, nrVector*dof*dim, false);
+      z.GetSubVector(vdofs_vectors, zval);
+      current_zmat.UseExternalData(zval.GetData(), dof*dim, nrVector);
+      current_ymat.SetSize(dof*dim, nrVector);
+      mfem::Mult(invmass_vec[i], current_zmat, current_ymat);
+      y.SetSubVector(vdofs_vectors, current_ymat.GetData());
+   }
+   max_char_speed = formIntegrator.GetMaxCharSpeed();
+   if (parallel)
+   {
+#ifdef MFEM_USE_MPI
+      MPI_Allreduce(MPI_IN_PLACE, &max_char_speed, 1, MFEM_MPI_REAL_T, MPI_MAX, comm);
+#endif
+   }
+}
+
+} // end of namespace mfem

miniapps/mani-swe/manihyp.hpp

@@ -0,0 +1,539 @@
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+void sphere(const Vector &x, Vector &y, const real_t r=1.0);
+
+/**
+ * @brief Extract orthogornal vector from B not belonging to A.
+ *
+ * @param A sub-subspace
+ * @param B subspace
+ * @param n a unit vector in B orthogornal to column space of A.
+ */
+void CalcOrtho(const DenseMatrix &A, const DenseMatrix &B, Vector &n);
+
+class ManifoldCoord
+{
+   // attributes
+private:
+   ElementTransformation *curr_el;
+   mutable DenseMatrix mani_vec_state;
+   mutable DenseMatrix phys_vec_state;
+   mutable Vector normal_comp;
+   mutable Vector phys_vec;
+protected:
+public:
+   const int dim;
+   const int sdim;
+
+   // methods
+private:
+protected:
+public:
+   ManifoldCoord(const int dim, const int sdim):dim(dim), sdim(sdim) {}
+
+   /**
+    * @brief Convert manifold state to physical state
+    *
+    * @param el Target element
+    * @param state  Current state value
+    * @param phys_state Current Physical state value
+    */
+   void convertElemState(ElementTransformation &Tr,
+                         const int nrScalar, const int nrVector,
+                         const Vector &state, Vector &phys_state) const;
+
+   /**
+    * @brief Convert left and right states to physical states
+    *
+    * In physical space, scalar spaces are the same as manifold state
+    * However, vector states are translated in local coordinates
+    * Basic conversion is v -> Jv
+    * To incoporate discontinuity of local coordinates along the interface,
+    * we convert state from one element (left) to another element (right) by
+    * v -> J1v -> J1v + (n1 dot Jv) (n2 - n1)
+    * where J1, n1 are Jacobian and normal vector from one element,
+    * and n2 is the normal vector from another element.
+    *
+    * @param Tr Interface transformation
+    * @param nrScalar The number of scalar states
+    * @param nrVector The number of vector states
+    * @param stateL Input left state
+    * @param stateR Input right state
+    * @param normalL **outward** normal from the left element
+    * @param normalR **inward** normal from the right element
+    * @param stateL_L **left** state in the __left__ coordinate system
+    * @param stateR_L **right** state in the __left__ coordinate system
+    * @param stateL_R **left** state in the __right__ coordinate system
+    * @param stateR_R **right** state in the __right__ coordinate system
+    */
+   void convertFaceState(FaceElementTransformations &Tr,
+                         const int nrScalar, const int nrVector,
+                         const Vector &stateL, const Vector &stateR,
+                         Vector &normalL, Vector &normalR,
+                         Vector &stateL_L, Vector &stateR_L,
+                         Vector &stateL_R, Vector &stateR_R) const;
+
+};
+
+class ManifoldFlux : public FluxFunction
+{
+   // attributes
+private:
+   FluxFunction &org_flux;
+   const ManifoldCoord &coord;
+   int nrScalar;
+   int nrVector;
+   mutable Vector phys_state;
+protected:
+public:
+
+   // methods
+private:
+protected:
+public:
+   ManifoldFlux(FluxFunction &flux, ManifoldCoord &coord, int nrScalar)
+      : FluxFunction(flux.num_equations, flux.dim), org_flux(flux),
+        coord(coord), nrScalar(nrScalar)
+   {
+      nrVector = (org_flux.num_equations - nrScalar)/coord.sdim;
+      phys_state.SetSize(nrScalar + nrVector*coord.sdim);
+   }
+
+   const ManifoldCoord &GetCoordinate() const {return coord;}
+
+   /**
+    * @brief Compute physical flux from manifold state
+    *
+    * @param state manifold state
+    * @param Tr local element transformation
+    * @param flux physical state
+    * @return maximum characteristic speed
+    */
+
+   real_t ComputeFlux(const Vector &state, ElementTransformation &Tr,
+                      DenseMatrix &flux) const override final;
+
+   real_t ComputeFluxDotN(const Vector &state, const Vector &normal,
+                          ElementTransformation &Tr,
+                          Vector &fluxDotN) const override final
+   {
+      MFEM_ABORT("Use ComputeNormalFluxes.");
+   }
+
+   real_t ComputeNormalFluxes(const Vector &stateL,
+                              const Vector &stateR,
+                              FaceElementTransformations &Tr,
+                              Vector &normalL, Vector &normalR,
+                              Vector &stateL_L, Vector &stateR_L,
+                              Vector &fluxL_L, Vector &fluxR_L,
+                              Vector &stateL_R, Vector &stateR_R,
+                              Vector &fluxL_R, Vector &fluxR_R) const;
+
+   int GetNumScalars() const { return nrScalar; }
+};
+
+
+class ManifoldNumericalFlux : public RiemannSolver
+{
+   // attributes
+private:
+protected:
+   const ManifoldFlux &maniflux;
+   mutable Vector fluxL_L, fluxR_L, fluxL_R, fluxR_R;
+   mutable Vector stateL_L, stateR_L, stateL_R, stateR_R;
+   mutable Vector normalL, normalR;
+public:
+
+   // methods
+private:
+protected:
+public:
+   ManifoldNumericalFlux(const ManifoldFlux &flux):RiemannSolver(flux),
+      maniflux(flux)
+   {
+      fluxL_L.SetSize(maniflux.num_equations);
+      fluxR_L.SetSize(maniflux.num_equations);
+      fluxL_R.SetSize(maniflux.num_equations);
+      fluxR_R.SetSize(maniflux.num_equations);
+      stateL_L.SetSize(maniflux.num_equations);
+      stateR_L.SetSize(maniflux.num_equations);
+      stateL_R.SetSize(maniflux.num_equations);
+      stateR_R.SetSize(maniflux.num_equations);
+      normalL.SetSize(maniflux.GetCoordinate().sdim);
+      normalR.SetSize(maniflux.GetCoordinate().sdim);
+   }
+   real_t Eval(const Vector &state1, const Vector &state2,
+               const Vector &nor, FaceElementTransformations &Tr,
+               Vector &flux) const final { MFEM_ABORT("Use the other Eval function") };
+   virtual real_t Eval(const Vector &stateL, const Vector &stateR,
+                       FaceElementTransformations &Tr,
+                       Vector &hatFL, Vector &hatFR) const = 0;
+
+   const ManifoldFlux &GetManifoldFluxFunction() const {return maniflux;}
+   const ManifoldCoord &GetCoordinate() const {return maniflux.GetCoordinate();}
+};
+
+class ManifoldRusanovFlux : public ManifoldNumericalFlux
+{
+   // attributes
+private:
+protected:
+public:
+
+   // methods
+private:
+protected:
+public:
+   ManifoldRusanovFlux(const ManifoldFlux &flux):ManifoldNumericalFlux(flux) {}
+   virtual real_t Eval(const Vector &stateL, const Vector &stateR,
+                       FaceElementTransformations &Tr,
+                       Vector &hatFL, Vector &hatFR) const override
+   {
+#ifdef MFEM_THREAD_SAFE
+      Vector fluxN1(fluxFunction.num_equations), fluxN2(fluxFunction.num_equations);
+#endif
+      const real_t maxE = maniflux.ComputeNormalFluxes(stateL, stateR, Tr,
+                                                       normalL, normalR,
+                                                       stateL_L, stateR_L, fluxL_L, fluxR_L,
+                                                       stateL_R, stateR_R, fluxL_R, fluxR_R);
+      // here, std::sqrt(nor*nor) is multiplied to match the scale with fluxN
+      const real_t scaledMaxE = maxE*Tr.Weight();
+      for (int i=0; i<maniflux.num_equations; i++)
+      {
+         hatFL[i] = 0.5*(scaledMaxE*(stateL_L[i] - stateR_L[i]) +
+                         (fluxL_L[i] + fluxR_L[i]));
+         hatFR[i] = 0.5*(scaledMaxE*(stateL_R[i] - stateR_R[i]) +
+                         (fluxL_R[i] + fluxR_R[i]));
+      }
+      return maxE;
+   }
+};
+
+class ManifoldVectorMassIntegrator : public BilinearFormIntegrator
+{
+protected:
+#ifndef MFEM_THREAD_SAFE
+   Vector shape, te_shape;
+   DenseMatrix elmat_comp, elmat_comp_weighted;
+   DenseMatrix JtJ;
+#endif
+   // PA extension
+   const FiniteElementSpace *fespace;
+   int dim, sdim, ne, nq;
+
+public:
+   ManifoldVectorMassIntegrator(const IntegrationRule *ir = NULL)
+      : BilinearFormIntegrator(ir) { }
+
+   /** Given a particular Finite Element computes the element mass matrix
+       elmat. */
+   virtual void AssembleElementMatrix(const FiniteElement &el,
+                                      ElementTransformation &Trans,
+                                      DenseMatrix &elmat);
+
+   static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
+                                         const FiniteElement &test_fe,
+                                         ElementTransformation &Trans);
+
+};
+class ManifoldVectorGradientIntegrator : public BilinearFormIntegrator
+{
+protected:
+#ifndef MFEM_THREAD_SAFE
+   Vector shape, te_shape;
+   DenseMatrix elmat_comp, elmat_comp_weighted;
+   DenseMatrix JtJ;
+#endif
+   // PA extension
+   const FiniteElementSpace *fespace;
+   int dim, sdim, ne, nq;
+
+public:
+   ManifoldVectorGradientIntegrator(const IntegrationRule *ir = NULL)
+      : BilinearFormIntegrator(ir) { }
+
+   /** Given a particular Finite Element computes the element mass matrix
+       elmat. */
+   virtual void AssembleElementMatrix(const FiniteElement &el,
+                                      ElementTransformation &Trans,
+                                      DenseMatrix &elmat);
+
+   static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
+                                         const FiniteElement &test_fe,
+                                         ElementTransformation &Trans);
+
+};
+
+class ManifoldHyperbolicFormIntegrator : public NonlinearFormIntegrator
+{
+   // attributes
+private:
+protected:
+   const ManifoldNumericalFlux &numFlux;
+   const ManifoldFlux &maniFlux;
+   const ManifoldCoord &coord;
+   real_t max_char_speed=0.0;
+   Vector state, phys_state;
+   Vector stateL, stateR;
+   Vector phys_stateL, phys_stateR;
+   Vector phys_hatFL, phys_hatFR;
+   Vector hatFL, hatFR;
+   Vector shape;
+   Vector shape1, shape2;
+   // DenseMatrix adjJ;
+   DenseMatrix dshape;
+   DenseMatrix gshape, vector_gshape, vector_gshape_comp;
+   DenseMatrix hess_shape;
+   DenseTensor HessMat;
+   DenseMatrix gradJ;
+   Vector x_nodes;
+   DenseMatrix phys_flux;
+   DenseMatrix phys_flux_scalars, phys_flux_vectors;
+   const IntegrationRule *intrule;
+   Array<int> hess_map;
+   DG_FECollection dg_fec;
+public:
+
+   // methods
+private:
+   const int GetElementIntegratioOrder(ElementTransformation &Trans,
+                                       const int order)
+   {
+      return Trans.OrderJ()+Trans.OrderW()+order;
+   }
+
+   const int GetFaceIntegratioOrder(FaceElementTransformations &Trans,
+                                    const int orderL, const int orderR)
+   {
+      return std::max(Trans.Elem1->OrderJ(),
+                      Trans.Elem2->Order())+Trans.OrderW() + std::max(orderL, orderR);
+   }
+   const IntegrationRule &GetRule(const FiniteElement &el1,
+                                  const FiniteElement &el2, FaceElementTransformations &Tr)
+   {
+      return IntRules.Get(Tr.GetGeometryType(),
+                          el1.GetOrder() + el2.GetOrder() + Tr.Elem1->OrderJ() + Tr.Elem2->OrderJ() +
+                          Tr.OrderW());
+   }
+   const IntegrationRule &GetRule(const FiniteElement &trial_fe,
+                                  const FiniteElement &test_fe, ElementTransformation &Trans)
+   {
+      const int order = trial_fe.GetOrder() + trial_fe.GetOrder() + Trans.OrderW() +
+                        Trans.OrderJ()*2;
+      return IntRules.Get(trial_fe.GetGeomType(), order);
+   }
+protected:
+public:
+   /**
+    * @brief Integrator of (F(u), grad v) - <\hat{F}(u), [v]> with given numerical flux.
+    * numerical flux both implements F(u) and numerical flux \hat{F}(u).
+    *
+    * @param flux Numerical flux
+    * @param ir Optionally chosen integration rule
+    */
+   ManifoldHyperbolicFormIntegrator(const ManifoldNumericalFlux &flux,
+                                    const IntegrationRule *ir=nullptr);
+
+   // Compute (F(u), grad v)
+   void AssembleElementVector(const FiniteElement &el, ElementTransformation &Tr,
+                              const Vector &elfun, Vector &elvect) override;
+
+   // Compute -<\hat{F}(u), [v]> with given numerical flux
+   void AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2,
+                           FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect) override;
+   // Get maximum characteristic speed for each processor.
+   // For parallel assembly, you need to use MPI_Allreduce to synchronize.
+   real_t GetMaxCharSpeed() { return max_char_speed; }
+
+   // Set max_char_speed to 0
+   void ResetMaxCharSpeed() { max_char_speed=0.0;}
+
+};
+
+class ManifoldStateCoefficient : public VectorCoefficient
+{
+   // attributes
+private:
+   VectorCoefficient &phys_cf;
+   Vector phys_state;
+   DenseMatrix mani_vecs, phys_vecs;
+   const int nrScalar, nrVector, dim, sdim;
+protected:
+public:
+
+   // methods
+private:
+protected:
+public:
+   ManifoldStateCoefficient(VectorCoefficient &phys_cf,
+                            const int nrScalar, const int nrVector, const int dim)
+      :VectorCoefficient(nrScalar + nrVector*dim), phys_cf(phys_cf),
+       phys_state(phys_cf.GetVDim()), nrScalar(nrScalar), nrVector(nrVector), dim(dim),
+       sdim((phys_cf.GetVDim()-nrScalar)/nrVector)
+   {}
+   virtual void Eval(Vector &mani_state, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      phys_cf.Eval(phys_state, T, ip);
+      for (int i=0; i<nrScalar; i++)
+      {
+         mani_state[i] = phys_state[i];
+      }
+      const DenseMatrix& invJ = T.InverseJacobian();
+      mani_vecs.UseExternalData(mani_state.GetData() + nrScalar, dim, nrVector);
+      phys_vecs.UseExternalData(phys_state.GetData() + nrScalar, sdim, nrVector);
+      Mult(invJ, phys_vecs, mani_vecs);
+   }
+
+};
+
+class ManifoldPhysVectorCoefficient : public VectorCoefficient
+{
+   // attributes
+private:
+   GridFunction &gf;
+   const int vid, dim, sdim;
+   Vector val;
+   Vector val_view;
+protected:
+public:
+
+   // methods
+private:
+protected:
+public:
+   ManifoldPhysVectorCoefficient(GridFunction &gf,
+                                 const int vid, const int dim, const int sdim)
+      :VectorCoefficient(sdim), gf(gf), vid(vid), dim(dim), sdim(sdim)
+   {
+      val.SetSize(gf.VectorDim());
+   }
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      gf.GetVectorValue(T, ip, val);
+      val_view.SetDataAndSize(val.GetData() + vid, dim);
+      T.Jacobian().Mult(val_view, V);
+   }
+
+};
+
+
+class CoriolisForce : public VectorCoefficient
+{
+private:
+   const real_t omega; // Coriolis parameter
+   ManifoldPhysVectorCoefficient &mom_cf;
+   Vector mom; // local momentum
+   Vector normal; // surface normal
+   Vector x;
+   Vector V_phys;
+   Vector V_mani;
+public:
+   CoriolisForce(ManifoldPhysVectorCoefficient &mom_cf,
+                 const real_t omega):VectorCoefficient(3), omega(omega), mom_cf(mom_cf),
+      mom(mom_cf.GetVDim()), normal(mom_cf.GetVDim()), x(mom_cf.GetVDim()),
+      V_phys(mom_cf.GetVDim()), V_mani(mom_cf.GetVDim()-1)
+   {
+      MFEM_ASSERT(mom_cf.GetVDim() == 3, "Momentum should be 3D");
+   }
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      V = 0.0;
+      T.Transform(T.GetIntPoint(), x);
+      const real_t theta = std::acos(x[2] / std::sqrt(x*x));
+      const real_t f = 2*omega*std::sin(theta);
+      mom_cf.Eval(mom, T, ip);
+      CalcOrtho(T.Jacobian(), normal);
+      normal /= std::sqrt(normal*normal);
+      normal.cross3D(mom, V_phys);
+      V_mani.SetData(V.GetData() + 1);
+      T.Jacobian().MultTranspose(V_phys, V_mani);
+      V_mani *= f;
+   }
+};
+
+class ManifoldDGHyperbolicConservationLaws : public TimeDependentOperator
+{
+private:
+   const int dim, sdim, nrScalar, nrVector;
+   FiniteElementSpace &vfes; // vector finite element space
+   // Element integration form. Should contain ComputeFlux
+   ManifoldHyperbolicFormIntegrator &formIntegrator;
+   // Base Nonlinear Form
+   std::unique_ptr<NonlinearForm> nonlinearForm;
+   std::unique_ptr<LinearForm> force;
+   // element-wise inverse mass matrix
+   int int_offset;
+   std::vector<DenseMatrix> invmass; // local scalar inverse mass
+   std::vector<DenseMatrix> invmass_vec; // local scalar inverse mass
+   std::vector<DenseMatrix> weakdiv; // local weak divergence
+   std::vector<DenseMatrix> weakdiv_vec; // local weak divergence
+   // global maximum characteristic speed. Updated by form integrators
+   mutable real_t max_char_speed;
+   // auxiliary variable used in Mult
+   mutable Vector z;
+
+   // Compute element-wise inverse mass matrix
+   void ComputeInvMass();
+   // Compute element-wise weak-divergence matrix
+   void ComputeWeakDivergence();
+   bool parallel = false;
+#ifdef MFEM_USE_MPI
+   MPI_Comm comm;
+#endif
+public:
+   /**
+    * @brief Construct a new DGHyperbolicConservationLaws object
+    *
+    * @param vfes_ vector finite element space. Only tested for DG [Pₚ]ⁿ
+    * @param formIntegrator_ integrator (F(u,x), grad v)
+    * @param preassembleWeakDivergence preassemble weak divergence for faster
+    *                                  assembly
+    */
+   ManifoldDGHyperbolicConservationLaws(
+      FiniteElementSpace &vfes,
+      ManifoldHyperbolicFormIntegrator &formIntegrator,
+      const int nrScalar,
+      const int order_offset=3);
+   /**
+    * @brief Apply nonlinear form to obtain M⁻¹(DIVF + JUMP HAT(F))
+    *
+    * @param x current solution vector
+    * @param y resulting dual vector to be used in an EXPLICIT solver
+    */
+   void Mult(const Vector &x, Vector &y) const override;
+   // get global maximum characteristic speed to be used in CFL condition
+   // where max_char_speed is updated during Mult.
+   real_t GetMaxCharSpeed() { return max_char_speed; }
+   void Update();
+   void AddForce(LinearFormIntegrator *lfdi)
+   {
+      if (force)
+      {
+         force->AddDomainIntegrator(lfdi);
+      }
+      else
+      {
+         if (parallel)
+         {
+#ifdef MFEM_USE_MPI
+            ParFiniteElementSpace *pvfes = static_cast<ParFiniteElementSpace*>(&vfes);
+            force.reset(new ParLinearForm(pvfes, z.GetData()));
+#endif
+         }
+         else
+         {
+            force.reset(new LinearForm(&vfes, z.GetData()));
+         }
+         force->AddDomainIntegrator(lfdi);
+      }
+   }
+};
+
+
+} // end of namespace mfem

miniapps/mani-swe/test_integrator.cpp

@@ -0,0 +1,117 @@
+#include "mfem.hpp"
+#include "manihyp.hpp"
+
+using namespace mfem;
+void gaussian_initial(const Vector &x, Vector &u)
+{
+   // const real_t theta = std::acos(x[2]/std::sqrt(x*x));
+   // const real_t hmin = 1;
+   // const real_t hmax = 2;
+   // const real_t sigma = 0.2;
+   // u = 0.0;
+   // u[0] = hmin + (hmax - hmin)*std::exp(-theta*theta/(2*sigma*sigma));
+
+   const real_t hmin = 1;
+   const real_t hmax = 1;
+   const real_t sigma = 0.2;
+   const real_t r2 = x[0]*x[0] + x[1]*x[1];
+   u = 1.0;
+   u[0] = hmin + (hmax - hmin)*std::exp(-r2/(2*sigma*sigma));
+   if (x.Size() == 3) { u[3] = 0.0; }
+}
+
+int main(int argc, char *argv[])
+{
+   const int order = 4;
+   const int dim = 2;
+   const int sdim = 3;
+   DG_FECollection dg_fec(order, dim);
+   const real_t h = 0.01;
+
+
+   const int nrElem = 1;
+   const int nrVert = 4;
+   const int nrBdr = 4;
+   VectorFunctionCoefficient u0_quad(dim + 1, gaussian_initial);
+   VectorFunctionCoefficient u0_mani_phys(sdim + 1, gaussian_initial);
+   ManifoldStateCoefficient u0_mani(u0_mani_phys, 1, 1, dim);
+
+   // 1. 2D Quad mesh
+   Mesh quad("./data/periodic-square-2d.mesh");
+   // Mesh quad(dim, nrVert, nrElem, nrBdr, dim);
+   // quad.AddVertex(0, 0);
+   // quad.AddVertex(1, 0);
+   // quad.AddVertex(1, 1);
+   // quad.AddVertex(0, 1);
+   //
+   // quad.AddQuad(0,1,2,3);
+   // 
+   // quad.AddBdrSegment(0, 1);
+   // quad.AddBdrSegment(1, 2);
+   // quad.AddBdrSegment(2, 3);
+   // quad.AddBdrSegment(3, 0);
+   //
+   // quad.FinalizeQuadMesh();
+   // quad.Save("quad.mesh");
+
+   FiniteElementSpace fes_quad(&quad, &dg_fec, dim+1);
+   GridFunction x_quad(&fes_quad), y_quad(&fes_quad);
+
+   NonlinearForm form_quad(&fes_quad);
+   ShallowWaterFlux swe_flux_quad(dim);
+
+   RusanovFlux swe_quad_numer(swe_flux_quad);
+   HyperbolicFormIntegrator swe_integ_quad(swe_quad_numer);
+
+   form_quad.AddDomainIntegrator(&swe_integ_quad);
+   form_quad.UseExternalIntegrators();
+   x_quad.ProjectCoefficient(u0_quad);
+   form_quad.Mult(x_quad, y_quad);
+
+   // 2. 3D Quad mesh
+   Mesh mani("./data/periodic-square-3d.mesh");
+   // Mesh mani(dim, nrVert, nrElem, nrBdr, sdim);
+   // mani.AddVertex(0, 0, 10.0);
+   // mani.AddVertex(1, 0, 10.0);
+   // mani.AddVertex(1, 1, 10.0);
+   // mani.AddVertex(0, 1, 10.0);
+   //
+   // mani.AddQuad(0,1,2,3);
+   //
+   // mani.AddBdrSegment(0, 1);
+   // mani.AddBdrSegment(1, 2);
+   // mani.AddBdrSegment(2, 3);
+   // mani.AddBdrSegment(3, 0);
+   //
+   // mani.FinalizeQuadMesh();
+   // mani.SetCurvature(order, true);
+   // mani.Save("mani.mesh");
+
+   FiniteElementSpace fes_mani(&mani, &dg_fec, dim+1);
+   GridFunction x_mani(&fes_mani), y_mani(&fes_mani);
+
+   NonlinearForm form_mani(&fes_mani);
+   ShallowWaterFlux swe_flux_phys(sdim);
+   ManifoldCoord coord(dim, sdim);
+   ManifoldFlux swe_flux_mani(swe_flux_phys, coord, 1);
+
+   ManifoldRusanovFlux swe_mani_numer(swe_flux_mani);
+   ManifoldHyperbolicFormIntegrator swe_integ_mani(swe_mani_numer);
+
+   form_mani.AddDomainIntegrator(&swe_integ_mani);
+   form_mani.UseExternalIntegrators();
+   x_mani.ProjectCoefficient(u0_mani);
+   form_mani.Mult(x_mani, y_mani);
+
+
+   // Logging
+   out << "X in quad: " << std::endl;
+   x_quad.Print(out, x_quad.Size());
+   out << "X in mani: " << std::endl;
+   x_mani.Print(out, x_mani.Size());
+   out << "Y in quad: " << std::endl;
+   y_quad.Print(out, y_quad.Size());
+   out << "Y in mani: " << std::endl;
+   y_mani.Print(out, y_mani.Size());
+   out << y_mani.DistanceTo(y_quad) << std::endl;
+}

miniapps/mani-swe/test_maniflux.cpp

@@ -0,0 +1,262 @@
+#include "mfem.hpp"
+#include "manihyp.hpp"
+#include <cmath>
+
+using namespace mfem;
+
+void gaussian_initial(const Vector &x, Vector &u)
+{
+   // const real_t theta = std::acos(x[2]/std::sqrt(x*x));
+   // const real_t hmin = 1;
+   // const real_t hmax = 2;
+   // const real_t sigma = 0.2;
+   // u = 0.0;
+   // u[0] = hmin + (hmax - hmin)*std::exp(-theta*theta/(2*sigma*sigma));
+
+   const real_t hmin = 1;
+   const real_t hmax = 2;
+   const real_t sigma = 0.2;
+   const real_t r2 = x[0]*x[0] + x[1]*x[1];
+   u = 0.0;
+   u[0] = hmin + (hmax - hmin)*std::exp(-r2/(2*sigma*sigma));
+}
+
+void UniformSpherRefinement(ParMesh &pmesh, int ref_level)
+{
+   ParGridFunction &x = static_cast<ParGridFunction&>(*pmesh.GetNodes());
+   VectorFunctionCoefficient sphere_cf(3, [](const Vector& x, Vector &y) {sphere(x,y,1.0);});
+   x.ProjectCoefficient(sphere_cf);
+   for (int i=0; i<ref_level; i++)
+   {
+      pmesh.UniformRefinement();
+      x.ProjectCoefficient(sphere_cf);
+   }
+}
+
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   const int numProcs = Mpi::WorldSize();
+   const int myRank = Mpi::WorldRank();
+   Hypre::Init();
+
+   int order = 3;
+   int refinement_level = 4;
+   bool visualization = true;
+   bool paraview = true;
+   real_t cfl = 0.3;
+   real_t tF = 1.5;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&refinement_level, "-r", "--refine",
+                  "Mesh refinement level");
+   args.ParseCheck();
+
+   std::unique_ptr<ParMesh> pmesh;
+   {
+      Mesh mesh("./data/periodic-square-3d.mesh");
+      mesh.SetCurvature(order, true);
+      pmesh.reset(new ParMesh(MPI_COMM_WORLD, mesh));
+      mesh.Clear();
+      for (int i=0; i<refinement_level; i++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+   std::unique_ptr<ParMesh> pmesh_visual;
+   {
+      Mesh mesh("./data/periodic-square-2d.mesh");
+      mesh.SetCurvature(order, true);
+      pmesh_visual.reset(new ParMesh(MPI_COMM_WORLD, mesh));
+      mesh.Clear();
+      for (int i=0; i<refinement_level; i++)
+      {
+         pmesh_visual->UniformRefinement();
+      }
+   }
+   {
+      // Mesh mesh("./data/icosahedron.mesh");
+      // mesh.SetCurvature(order);
+      // pmesh.reset(new ParMesh(MPI_COMM_WORLD, mesh));
+      // mesh.Clear();
+      // UniformSpherRefinement(*pmesh, refinement_level);
+   }
+
+   const int dim = pmesh->Dimension();
+   const int sdim = pmesh->SpaceDimension();
+   const int num_equations = dim + 1;
+   const int phys_num_equations = sdim + 1;
+
+   ManifoldCoord coord(dim, sdim);
+   ShallowWaterFlux swe_phys(sdim);
+   ManifoldFlux swe_mani(swe_phys, coord, 1);
+   ManifoldRusanovFlux rusanovFlux(swe_mani);
+   ManifoldHyperbolicFormIntegrator swe_integ(rusanovFlux);
+
+   std::unique_ptr<ODESolver> ode_solver;
+   ode_solver.reset(new RK4Solver());
+
+   DG_FECollection dg_fec(order, dim);
+   // FE Space for state
+   ParFiniteElementSpace vfes(pmesh.get(), &dg_fec, num_equations,
+                              Ordering::byNODES);
+   // FE space for manifold vector
+   ParFiniteElementSpace dfes(pmesh.get(), &dg_fec, dim, Ordering::byNODES);
+   // FE space for physical vector
+   ParFiniteElementSpace sfes(pmesh.get(), &dg_fec, sdim, Ordering::byNODES);
+   // FE space for scalar
+   ParFiniteElementSpace fes(pmesh.get(), &dg_fec);
+
+   // State
+   ParGridFunction u(&vfes);
+   // Height for visualization
+   ParGridFunction height(&fes, u.GetData());
+   ParGridFunction mom(&sfes);
+   ParGridFunction mom_x(&fes, mom.GetData() + 0*fes.GetTrueVSize());
+   ParGridFunction mom_y(&fes, mom.GetData() + 1*fes.GetTrueVSize());
+   ParGridFunction mom_z(&fes, mom.GetData() + 2*fes.GetTrueVSize());
+   ManifoldPhysVectorCoefficient mom_cf(u, 1, dim, sdim);
+
+   VectorFunctionCoefficient u0_phys(phys_num_equations, gaussian_initial);
+   ManifoldStateCoefficient u0_mani(u0_phys, 1, 1, dim);
+   u.ProjectCoefficient(u0_mani);
+
+   ManifoldDGHyperbolicConservationLaws swe(vfes, swe_integ, 1);
+   swe.SetTime(0.0);
+   real_t hmin=infinity();
+   {
+      for (int i=0; i<pmesh->GetNE(); i++)
+      {
+         hmin = std::min(pmesh->GetElementSize(i, 1), hmin);
+      }
+      MPI_Allreduce(MPI_IN_PLACE, &hmin, 1, MFEM_MPI_REAL_T, MPI_MIN,
+                    pmesh->GetComm());
+      Vector z(vfes.GetTrueVSize());
+      swe.Mult(u,z);
+   }
+
+   socketstream height_sock, mom_x_sock, mom_y_sock, mom_z_sock;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      MPI_Barrier(MPI_COMM_WORLD);
+      height_sock.open(vishost, visport);
+      height_sock.precision(8);
+      // Plot height
+      height_sock << "parallel " << numProcs << " " << myRank << "\n";
+      height_sock << "solution\n" << *pmesh_visual << height;
+      height_sock << "window_title 'momentum, t = 0'\n";
+      height_sock << "view 0 0\n";  // view from top
+      height_sock << "autoscale off\n valuerange 0.5 2.5\n";
+      height_sock << "keys jm\n";  // turn off perspective and light, show mesh
+      height_sock << std::flush;
+      MPI_Barrier(MPI_COMM_WORLD);
+
+      // MPI_Barrier(MPI_COMM_WORLD);
+      // mom.ProjectCoefficient(mom_cf);
+      // mom_x_sock.open(vishost, visport);
+      // mom_x_sock.precision(8);
+      // // Plot magnitude of vector-valued momentum
+      // mom_x_sock << "parallel " << numProcs << " " << myRank << "\n";
+      // mom_x_sock << "solution\n" << *pmesh << mom_x;
+      // mom_x_sock << "window_title 'momentum_x, t = 0'\n";
+      // mom_x_sock << "view 0 0\n";  // view from top
+      // mom_x_sock << "keys jm\n";  // turn off perspective and light, show mesh
+      // mom_x_sock << std::flush;
+      // MPI_Barrier(MPI_COMM_WORLD);
+      //
+      // MPI_Barrier(MPI_COMM_WORLD);
+      // mom_y_sock.open(vishost, visport);
+      // mom_y_sock.precision(8);
+      // // Plot magnitude of vector-valued momentum
+      // mom_y_sock << "parallel " << numProcs << " " << myRank << "\n";
+      // mom_y_sock << "solution\n" << *pmesh << mom_y;
+      // mom_y_sock << "window_title 'momentum_y, t = 0'\n";
+      // mom_y_sock << "view 0 0\n";  // view from top
+      // mom_y_sock << "keys jm\n";  // turn off perspective and light, show mesh
+      // mom_y_sock << std::flush;
+      // MPI_Barrier(MPI_COMM_WORLD);
+      //
+      // MPI_Barrier(MPI_COMM_WORLD);
+      // mom_z_sock.open(vishost, visport);
+      // mom_z_sock.precision(8);
+      // // Plot magnitude of vector-valued momentum
+      // mom_z_sock << "parallel " << numProcs << " " << myRank << "\n";
+      // mom_z_sock << "solution\n" << *pmesh << mom_z;
+      // mom_z_sock << "window_title 'momentum_z, t = 0'\n";
+      // mom_z_sock << "view 0 0\n";  // view from top
+      // mom_z_sock << "keys jm\n";  // turn off perspective and light, show mesh
+      // mom_z_sock << std::flush;
+      // MPI_Barrier(MPI_COMM_WORLD);
+   }
+   std::unique_ptr<ParaViewDataCollection> dacol;
+   if (paraview)
+   {
+      dacol.reset(new ParaViewDataCollection("ParaViewSWE", pmesh_visual.get()));
+      dacol->SetLevelsOfDetail(order);
+      dacol->RegisterField("Height", &height);
+      dacol->RegisterField("Momentum", &mom);
+      dacol->SetTime(0.0);
+      dacol->SetCycle(0);
+      dacol->Save();
+   }
+
+   real_t t = 0.0;
+   real_t dt; // dt will be computed using CFL condition
+   swe.SetTime(t);
+   ode_solver->Init(swe);
+   bool done = false;
+
+   for (int ti = 0; !done; ti++)
+   {
+      // CFL condition
+      dt = cfl * hmin / swe.GetMaxCharSpeed() / real_t(2*order+1);
+      // Adjust dt so that t + dt <= tF
+      real_t dt_real = std::min(dt, tF - t);
+      if (Mpi::Root())
+      {
+         out << "time step: " << ti << ", time: " << t << std::endl;
+         out << "\tMaxChar: " << swe.GetMaxCharSpeed() << std::endl;
+      }
+
+      // ODE step
+      ode_solver->Step(u, t, dt_real);
+      done = (t >= tF - 1e-8 * dt);
+
+
+      // Visualize
+      // mom.ProjectCoefficient(mom_cf);
+      if (height_sock.is_open() && height_sock.good())
+      {
+         height_sock << "parallel " << numProcs << " " << myRank << "\n";
+         height_sock << "solution\n" << *pmesh_visual << height;
+         height_sock << "window_title 'height, t = " << t << "'\n";
+      }
+      // if (mom_x_sock.is_open() && mom_x_sock.good())
+      // {
+      //    mom_x_sock << "parallel " << numProcs << " " << myRank << "\n";
+      //    mom_x_sock << "solution\n" << *pmesh << mom_x;
+      //    mom_x_sock << "window_title 'momentum_x, t = " << t << "'\n";
+      // }
+      // if (mom_y_sock.is_open() && mom_y_sock.good())
+      // {
+      //    mom_y_sock << "parallel " << numProcs << " " << myRank << "\n";
+      //    mom_y_sock << "solution\n" << *pmesh << mom_y;
+      //    mom_y_sock << "window_title 'momentum_y, t = " << t << "'\n";
+      // }
+      // if (mom_z_sock.is_open() && mom_z_sock.good())
+      // {
+      //    mom_z_sock << "parallel " << numProcs << " " << myRank << "\n";
+      //    mom_z_sock << "solution\n" << *pmesh << mom_z;
+      //    mom_z_sock << "window_title 'momentum_z, t = " << t << "'\n";
+      // }
+      if (dacol)
+      {
+         dacol->SetTime(t);
+         dacol->SetCycle(ti);
+         dacol->Save();
+      }
+   }
+}

miniapps/mani-swe/test_sphere.cpp

@@ -0,0 +1,332 @@
+#include "mfem.hpp"
+#include "manihyp.hpp"
+#include <cmath>
+
+using namespace mfem;
+
+void gaussian_initial(const Vector &x, Vector &u)
+{
+   const real_t theta = std::acos(x[2]/std::sqrt(x*x));
+   const real_t hmin = 1;
+   const real_t hmax = 2;
+   const real_t sigma = 0.2;
+   u = 0.0;
+   u[0] = hmin + (hmax - hmin)*std::exp(-theta*theta/(2*sigma*sigma));
+}
+
+inline int sgn(const real_t x){ return x >= 0 ? 1 : -1; }
+
+void williamsTest1(const Vector &x, Vector &u)
+{
+   const real_t meter = 1.0 / 6.37122e6;
+   const real_t hour = 1;
+   const real_t second = hour / 3600.0;
+   const real_t R = 6.37122e6 * meter;
+   const real_t omega = 7.292e-05 / second;
+   const real_t g = 9.80616 * meter / second / second;
+   const real_t H = 1e4 * meter;
+
+   const real_t latitude = std::acos(x[2] / std::sqrt(x*x));
+   const real_t longitude = std::acos(x[0] / std::sqrt(x[0]*x[0] + x[1]*x[1]))*sgn(x[1]);
+
+   // const real_t u0 = 2*meter / second;
+   // const real_t h0 = 1e4 * meter;
+   // const real_t dh = 100 * meter;
+   //
+   // u[0] = h0 + dh * std::pow(std::cos(latitude), 2.0) * std::cos(longitude);
+   // u[1] =      u0 * std::pow(std::cos(latitude), 2.0) * std::cos(longitude);
+   // u[2] =      u0 * std::pow(std::cos(latitude), 2.0) * std::sin(longitude);
+   // u[3] =      u0 * std::pow(std::cos(latitude), 1.0) * std::sin(latitude);
+
+   const real_t lat_m = M_PI / 9.0;
+   const real_t h0 = 2000 * meter;
+   const real_t dh = 200 * meter;
+   const real_t rad2degree = 180 * M_1_PI;
+   u = 0.0;
+   u[0] = h0 + (std::fabs(latitude) < lat_m ? dh*std::exp(-latitude*latitude/(lat_m*lat_m*(lat_m*lat_m - latitude*latitude)*rad2degree*rad2degree)) : 0.0);
+   // const real_t lat0 = M_PI / 7.0;
+   // const real_t lat1 = M_PI / 2.0 - lat0;
+   // const real_t lat2 = M_PI / 4.0;
+   //
+   // const real_t lat = M_PI / 2.0 - theta;
+   //
+   // const real_t hpert = 120.0 * meter;
+   // const real_t alpha = 1.0 / 3.0;
+   // const real_t beta = 1.0 / 15.0;
+   //
+   // const real_t umax = 80.0 * meter / second;
+   // const real_t en = std::exp(-4.0 / std::pow(lat1-lat0, 2.0));
+   // bool jet = (lat0 <= lat) && (lat <= lat1);
+   // const real_t u_jet = umax / en * std::exp(1.0 / (lat - lat0) / (lat - lat1));
+   //
+   // const real_t dtheta = M_PI/100;
+   // real_t theta2 = 0.0;
+}
+
+class SphericalHeight : public VectorCoefficient
+{
+private:
+   GridFunctionCoefficient h;
+   const real_t scale;
+   Vector normal;
+public:
+   SphericalHeight(GridFunction &h, const real_t scale=1.0):VectorCoefficient(3),
+      h(&h), scale(scale), normal(3) { }
+   virtual void Eval(Vector &node, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      T.Transform(T.GetIntPoint(), node);
+      node /= std::sqrt(node*node);
+      normal = node;
+      normal *= scale * h.Eval(T, T.GetIntPoint()) / std::sqrt(normal*normal);
+      node += normal;
+   }
+};
+
+void UniformSpherRefinement(ParMesh &pmesh, int ref_level)
+{
+   ParGridFunction &x = static_cast<ParGridFunction&>(*pmesh.GetNodes());
+   VectorFunctionCoefficient sphere_cf(3, [](const Vector& x, Vector &y) {sphere(x,y,1.0);});
+   x.ProjectCoefficient(sphere_cf);
+   for (int i=0; i<ref_level; i++)
+   {
+      pmesh.UniformRefinement();
+      x.ProjectCoefficient(sphere_cf);
+   }
+}
+
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   const int numProcs = Mpi::WorldSize();
+   const int myRank = Mpi::WorldRank();
+   Hypre::Init();
+
+   const real_t meter = 1.0 / 6.37122e6;
+   const real_t hour = 1;
+   const real_t second = hour / 3600.0;
+   const real_t R = 6.37122e6 * meter;
+   const real_t omega = 7.292e-05 / second;
+   const real_t g = 9.80616 * meter / second;
+   const real_t H = 1e4 * meter;
+
+   int order = 3;
+   int refinement_level = 4;
+   int vis_step = 50;
+   bool visualization = true;
+   bool paraview = true;
+   real_t cfl = 0.2;
+   real_t tF = 360*hour;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&refinement_level, "-r", "--refine",
+                  "Mesh refinement level");
+   args.ParseCheck();
+
+   std::unique_ptr<ParMesh> pmesh;
+   {
+      Mesh mesh("./data/icosahedron.mesh");
+      mesh.SetCurvature(order, true);
+      pmesh.reset(new ParMesh(MPI_COMM_WORLD, mesh));
+      mesh.Clear();
+      UniformSpherRefinement(*pmesh, refinement_level);
+   }
+   std::unique_ptr<ParMesh> pmesh_visualize;
+   {
+      Mesh mesh("./data/icosahedron.mesh");
+      mesh.SetCurvature(order+4, true);
+      pmesh_visualize.reset(new ParMesh(MPI_COMM_WORLD, mesh));
+      mesh.Clear();
+      UniformSpherRefinement(*pmesh_visualize, refinement_level);
+   }
+   GridFunction *nodes = pmesh_visualize->GetNodes();
+
+   const int dim = pmesh->Dimension();
+   const int sdim = pmesh->SpaceDimension();
+   const int num_equations = dim + 1;
+   const int phys_num_equations = sdim + 1;
+
+   ManifoldCoord coord(dim, sdim);
+   ShallowWaterFlux swe_phys(sdim, g);
+   ManifoldFlux swe_mani(swe_phys, coord, 1);
+   ManifoldRusanovFlux rusanovFlux(swe_mani);
+   ManifoldHyperbolicFormIntegrator swe_integ(rusanovFlux);
+
+   std::unique_ptr<ODESolver> ode_solver;
+   ode_solver.reset(new RK4Solver());
+
+   DG_FECollection dg_fec(order, dim);
+   // FE Space for state
+   ParFiniteElementSpace vfes(pmesh.get(), &dg_fec, num_equations,
+                              Ordering::byNODES);
+   // FE space for manifold vector
+   ParFiniteElementSpace dfes(pmesh.get(), &dg_fec, dim, Ordering::byNODES);
+   // FE space for physical vector
+   ParFiniteElementSpace sfes(pmesh.get(), &dg_fec, sdim, Ordering::byNODES);
+   // FE space for scalar
+   ParFiniteElementSpace fes(pmesh.get(), &dg_fec);
+
+   // State
+   ParGridFunction u(&vfes);
+   // Height for visualization
+   ParGridFunction height(&fes, u.GetData());
+   ParGridFunction mom(&sfes);
+   ParGridFunction mom_x(&fes, mom.GetData() + 0*fes.GetTrueVSize());
+   ParGridFunction mom_y(&fes, mom.GetData() + 1*fes.GetTrueVSize());
+   ParGridFunction mom_z(&fes, mom.GetData() + 2*fes.GetTrueVSize());
+   ManifoldPhysVectorCoefficient mom_cf(u, 1, dim, sdim);
+   SphericalHeight deform_cf(height);
+
+   VectorFunctionCoefficient u0_phys(phys_num_equations, gaussian_initial);
+   ManifoldStateCoefficient u0_mani(u0_phys, 1, 1, dim);
+   u.ProjectCoefficient(u0_mani);
+
+   ManifoldDGHyperbolicConservationLaws swe(vfes, swe_integ, 1);
+   CoriolisForce force(mom_cf, omega);
+   // swe.AddForce(new VectorDomainLFIntegrator(force));
+   swe.SetTime(0.0);
+   real_t hmin=infinity();
+   {
+      for (int i=0; i<pmesh->GetNE(); i++)
+      {
+         hmin = std::min(pmesh->GetElementSize(i, 1), hmin);
+      }
+      MPI_Allreduce(MPI_IN_PLACE, &hmin, 1, MFEM_MPI_REAL_T, MPI_MIN,
+                    pmesh->GetComm());
+      Vector z(vfes.GetTrueVSize());
+      swe.Mult(u,z);
+   }
+
+   socketstream height_sock, mom_x_sock, mom_y_sock, mom_z_sock;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      MPI_Barrier(MPI_COMM_WORLD);
+      height_sock.open(vishost, visport);
+      height_sock.precision(8);
+      // Plot height
+      nodes->ProjectCoefficient(deform_cf);
+      height_sock << "parallel " << numProcs << " " << myRank << "\n";
+      height_sock << "solution\n" << *pmesh_visualize << height;
+      height_sock << "window_title 'momentum, t = 0'\n";
+      height_sock << "view 0 0\n";  // view from top
+      // height_sock << "autoscale off\n valuerange 0.5 2.5\n";
+      height_sock << "keys jm\n";  // turn off perspective and light, show mesh
+      height_sock << std::flush;
+      MPI_Barrier(MPI_COMM_WORLD);
+
+      // MPI_Barrier(MPI_COMM_WORLD);
+      // mom.ProjectCoefficient(mom_cf);
+      // mom_x_sock.open(vishost, visport);
+      // mom_x_sock.precision(8);
+      // // Plot magnitude of vector-valued momentum
+      // mom_x_sock << "parallel " << numProcs << " " << myRank << "\n";
+      // mom_x_sock << "solution\n" << *pmesh << mom_x;
+      // mom_x_sock << "window_title 'momentum_x, t = 0'\n";
+      // mom_x_sock << "view 0 0\n";  // view from top
+      // mom_x_sock << "keys jm\n";  // turn off perspective and light, show mesh
+      // mom_x_sock << std::flush;
+      // MPI_Barrier(MPI_COMM_WORLD);
+      //
+      // MPI_Barrier(MPI_COMM_WORLD);
+      // mom_y_sock.open(vishost, visport);
+      // mom_y_sock.precision(8);
+      // // Plot magnitude of vector-valued momentum
+      // mom_y_sock << "parallel " << numProcs << " " << myRank << "\n";
+      // mom_y_sock << "solution\n" << *pmesh << mom_y;
+      // mom_y_sock << "window_title 'momentum_y, t = 0'\n";
+      // mom_y_sock << "view 0 0\n";  // view from top
+      // mom_y_sock << "keys jm\n";  // turn off perspective and light, show mesh
+      // mom_y_sock << std::flush;
+      // MPI_Barrier(MPI_COMM_WORLD);
+      //
+      // MPI_Barrier(MPI_COMM_WORLD);
+      // mom_z_sock.open(vishost, visport);
+      // mom_z_sock.precision(8);
+      // // Plot magnitude of vector-valued momentum
+      // mom_z_sock << "parallel " << numProcs << " " << myRank << "\n";
+      // mom_z_sock << "solution\n" << *pmesh << mom_z;
+      // mom_z_sock << "window_title 'momentum_z, t = 0'\n";
+      // mom_z_sock << "view 0 0\n";  // view from top
+      // mom_z_sock << "keys jm\n";  // turn off perspective and light, show mesh
+      // mom_z_sock << std::flush;
+      // MPI_Barrier(MPI_COMM_WORLD);
+   }
+   std::unique_ptr<ParaViewDataCollection> dacol;
+   if (paraview)
+   {
+      std::stringstream paraviewname;
+      paraviewname << "ParaViewSWE_WO_Coriolis_" << refinement_level;
+      dacol.reset(new ParaViewDataCollection(paraviewname.str().c_str(), pmesh.get()));
+      dacol->SetLevelsOfDetail(order);
+      dacol->RegisterField("Height", &height);
+      dacol->RegisterField("Momentum", &mom);
+      dacol->SetTime(0.0);
+      dacol->SetCycle(0);
+      dacol->Save();
+   }
+
+   real_t t = 0.0;
+   real_t dt; // dt will be computed using CFL condition
+   swe.SetTime(t);
+   ode_solver->Init(swe);
+   bool done = false;
+
+   for (int ti = 1; !done; ti++)
+   {
+      // CFL condition
+      dt = cfl * hmin / swe.GetMaxCharSpeed() / real_t(2*order+1);
+      // Adjust dt so that t + dt <= tF
+      real_t dt_real = std::min(dt, tF - t);
+      if (Mpi::Root())
+      {
+         out << "time step: " << ti << ", time: " << t << std::endl;
+         out << "\tMaxChar: " << swe.GetMaxCharSpeed() << std::endl;
+      }
+
+      // ODE step
+      ode_solver->Step(u, t, dt_real);
+      done = (t >= tF - 1e-8 * dt);
+
+
+      // Visualize
+      if (ti % vis_step == 0 || done)
+      {
+         mom.ProjectCoefficient(mom_cf);
+         if (height_sock.is_open() && height_sock.good())
+         {
+            nodes->ProjectCoefficient(deform_cf);
+            height_sock << "parallel " << numProcs << " " << myRank << "\n";
+            height_sock << "solution\n" << *pmesh_visualize << height;
+            height_sock << "window_title 'height, t = " << t << "'\n";
+         }
+         // if (mom_x_sock.is_open() && mom_x_sock.good())
+         // {
+         //    mom_x_sock << "parallel " << numProcs << " " << myRank << "\n";
+         //    mom_x_sock << "solution\n" << *pmesh << mom_x;
+         //    mom_x_sock << "window_title 'momentum_x, t = " << t << "'\n";
+         // }
+         // if (mom_y_sock.is_open() && mom_y_sock.good())
+         // {
+         //    mom_y_sock << "parallel " << numProcs << " " << myRank << "\n";
+         //    mom_y_sock << "solution\n" << *pmesh << mom_y;
+         //    mom_y_sock << "window_title 'momentum_y, t = " << t << "'\n";
+         // }
+         // if (mom_z_sock.is_open() && mom_z_sock.good())
+         // {
+         //    mom_z_sock << "parallel " << numProcs << " " << myRank << "\n";
+         //    mom_z_sock << "solution\n" << *pmesh << mom_z;
+         //    mom_z_sock << "window_title 'momentum_z, t = " << t << "'\n";
+         // }
+         if (dacol)
+         {
+            dacol->SetTime(t);
+            dacol->SetCycle(ti);
+            dacol->Save();
+         }
+
+      }
+   }
+}

miniapps/mani-swe/test_square.cpp

@@ -0,0 +1,51 @@
+#include "mfem.hpp"
+#include "manihyp.hpp"
+
+using namespace mfem;
+
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   int order = 8;
+   int refinement_level = 3;
+
+   Hypre::Init();
+   std::unique_ptr<ParMesh> pmesh;
+   {
+      Mesh mesh("./data/periodic-square-3d.mesh");
+      pmesh.reset(new ParMesh(MPI_COMM_WORLD, mesh));
+      mesh.Clear();
+   }
+
+   const int dim = pmesh->Dimension();
+   const int sdim = pmesh->SpaceDimension();
+
+
+   pmesh->SetCurvature(order, true);
+   for (int i=0; i<refinement_level; i++)
+   {
+      pmesh->UniformRefinement();
+   }
+   bool visualization = true;
+   DG_FECollection dg_fec(order, dim);
+   ParFiniteElementSpace pfes(pmesh.get(), &dg_fec);
+   ParGridFunction x(&pfes);
+   x = 1.0;
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *pmesh << x
+               << "keys 'mj'"
+               << std::flush;
+   }
+}