Merge remote-tracking branch 'origin/master' into div-free-rt-dev

examples/CMakeLists.txt

@@ -28,7 +28,6 @@ list(APPEND ALL_EXE_SRCS
   ex19.cpp
   ex20.cpp
   ex22.cpp
-  ex23.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -54,7 +53,6 @@ if (MFEM_USE_MPI)
     ex19p.cpp
     ex20p.cpp
     ex22p.cpp
-    ex23p.cpp
     )
 endif()
 

examples/ex23.cpp

@@ -1,734 +0,0 @@
-//                               MFEM Example 23
-//
-// Compile with: make ex23
-//
-// Sample runs:
-//    ex23 -m ../data/periodic-segment.mesh -p 0 -s 2 -dt 0.001 -vs 50
-//    ex23 -m ../data/periodic-segment.mesh -p 0 -s 12 -dt 0.01
-//    ex23 -m ../data/periodic-segment.mesh -p 0 -s 22 -dt 0.01
-//    ex23 -m ../data/periodic-segment.mesh -p 0 -s 32 -dt 0.005 -vs 10
-//    ex23 -m ../data/periodic-square.mesh -p 0 -dt 0.01
-//    ex23 -m ../data/periodic-square.mesh -p 0 -s 32 -dt 0.01
-//    ex23 -m ../data/periodic-hexagon.mesh -p 0 -d 0.001 -s 12 -dt 0.02
-//    ex23 -m ../data/periodic-hexagon.mesh -p 0 -d 0.001 -s 32 -dt 0.009 -vs 10
-//    ex23 -m ../data/periodic-square.mesh -p 1 -dt 0.01 -tf 9
-//    ex23 -m ../data/periodic-hexagon.mesh -p 1 -dt 0.01 -tf 9
-//    ex23 -m ../data/amr-quad.mesh -p 1 -dt 0.01 -tf 9 -vs 2
-//    ex23 -m ../data/disc-nurbs.mesh -p 1 -r 3 -dt 0.01 -tf 9
-//    ex23 -m ../data/disc-nurbs.mesh -p 2 -r 3 -dt 0.01 -tf 9
-//    ex23 -m ../data/disc-nurbs.mesh -p 3 -r 3 -dt 0.01 -tf 9 -d 0.02
-//    ex23 -m ../data/periodic-square.mesh -p 3 -r 3 -dt 0.025 -tf 9
-//    ex23 -m ../data/periodic-cube.mesh -p 0 -o 2 -dt 0.025 -tf 8
-//
-// Description:  This example code solves the time-dependent advection-diffusion
-//               equation
-//               du/dt - div(D grad(u)) + v.grad(u) = 0, where
-//               D is a diffusion coefficient,
-//               v is a given fluid velocity, and
-//               u0(x)=u(0,x) is a given initial condition.
-//
-//               The example demonstrates the use of Discontinuous Galerkin (DG)
-//               bilinear forms in MFEM (face integrators), the use of explicit,
-//               implicit, and implicit-explicit ODE time integrators, the
-//               definition of periodic boundary conditions through periodic
-//               meshes, as well as the use of GLVis for persistent
-//               visualization of a time-evolving solution. The saving of
-//               time-dependent data files for external visualization with
-//               VisIt (visit.llnl.gov) is also illustrated.
-//
-//               This example is a merger of examples 9 and 14.
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-// Choice for the problem setup. The fluid velocity, initial condition and
-// boundary condition are chosen based on this parameter.
-int problem;
-
-// Velocity coefficient
-void velocity_function(const Vector &x, Vector &v);
-
-// Initial condition
-double u0_function(const Vector &x);
-
-// Mesh bounding box
-Vector bb_min, bb_max;
-
-
-/** A time-dependent operator for the right-hand side of the ODE for use with
-    explicit ODE solvers. The DG weak form of du/dt = div(D grad(u))-v.grad(u) is
-    M du/dt = - S u + K u + b, where M, S, and K are the mass,
-    stiffness, and advection matrices, and b describes sources and the flow on
-    the boundary.
-    This can be written as a general ODE,
-    du/dt = M^{-1} (-S u + K u + b), and this class is used to compute the RHS
-    and perform the solve for du/dt. */
-class EX_Evolution : public TimeDependentOperator
-{
-private:
-   SparseMatrix &M, &S, &K;
-   const Vector &b;
-
-   DSmoother M_prec;
-   CGSolver M_solver;
-
-   mutable Vector z;
-
-   void initA(double dt);
-
-public:
-   EX_Evolution(SparseMatrix &_M, SparseMatrix &_S, SparseMatrix &_K,
-                const Vector &_b);
-
-   virtual void Mult(const Vector &x, Vector &y) const;
-
-   virtual ~EX_Evolution() {}
-};
-
-/** A time-dependent operator for the right-hand side of the ODE for use with
-    implicit ODE solvers. The DG weak form of du/dt = div(D grad(u))-v.grad(u) is
-    [M + dt (S - K)] du/dt = - S u + K u + b, where M, S, and K are the mass,
-    stiffness, and advection matrices, and b describes sources and the flow on
-    the boundary.
-    This can be written as a general ODE,
-    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
-    used to perform the fully implicit solve for du/dt. */
-class IM_Evolution : public TimeDependentOperator
-{
-private:
-   SparseMatrix &M, &S, &K;
-   SparseMatrix *A;
-   const Vector &b;
-
-   DSmoother M_prec;
-   CGSolver M_solver;
-
-   DSmoother *A_prec;
-   GMRESSolver *A_solver;
-   double dt;
-
-   mutable Vector z;
-
-   void initA(double dt);
-
-public:
-   IM_Evolution(SparseMatrix &_M, SparseMatrix &_S, SparseMatrix &_K,
-                const Vector &_b);
-
-   virtual void Mult(const Vector &x, Vector &y) const;
-
-   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
-
-   virtual ~IM_Evolution() { delete A_solver; delete A_prec; delete A; }
-};
-
-/** A time-dependent operator for the right-hand side of the ODE for use with
-    IMEX (Implicit-Explicit) ODE solvers. The DG weak form of
-    du/dt = div(D grad(u))-v.grad(u) is
-    [M + dt S] du/dt = - S u + K u + b, where M, S, and K are the mass,
-    stiffness, and advection matrices, and b describes sources and the flow on
-    the boundary.
-    This can be written as a general ODE,
-    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
-    used to perform the implicit or explicit solve for du/dt. */
-class IMEX_Evolution : public TimeDependentOperator
-{
-private:
-   SparseMatrix &M, &S, &K;
-   SparseMatrix *A;
-   const Vector &b;
-
-   DSmoother M_prec;
-   CGSolver M_solver;
-
-   DSmoother *A_prec;
-   CGSolver *A_solver;
-   double dt;
-
-   mutable Vector z;
-
-   void initA(double dt);
-
-public:
-   IMEX_Evolution(SparseMatrix &_M, SparseMatrix &_S, SparseMatrix &_K,
-                  const Vector &_b);
-
-   virtual void ExplicitMult(const Vector &x, Vector &y) const;
-   virtual void Mult(const Vector &x, Vector &y) const;
-
-   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
-
-   virtual ~IMEX_Evolution() { delete A_solver; delete A_prec; delete A; }
-};
-
-
-int main(int argc, char *argv[])
-{
-   // 1. Parse command-line options.
-   problem = 0;
-   const char *mesh_file = "../data/periodic-hexagon.mesh";
-   int ref_levels = 2;
-   int order = 3;
-   int ode_solver_type = 12;
-   double t_final = 10.0;
-   double d_coef = 0.01;
-   double dt = 0.01;
-   double sigma = -1.0;
-   double kappa = -1.0;
-   bool visualization = true;
-   bool visit = false;
-   bool binary = false;
-   int vis_steps = 5;
-
-   int precision = 8;
-   cout.precision(precision);
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&problem, "-p", "--problem",
-                  "Problem setup to use. See options in velocity_function().");
-   args.AddOption(&ref_levels, "-r", "--refine",
-                  "Number of times to refine the 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 Euler, 2 - RK2, 3 - RK3 SSP,"
-                  " 4 - RK4, 5 - Generalized Alpha,\n\t"
-                  "11 - Backward Euler, 12 - SDIRK2, 13 - SDIRK3,\n\t"
-                  "22 - Implicit Midpoint, 23 SDIRK23, 24 - SDIRK34,\n\t"
-                  "31 - IMEX BE/FE, 32 - IMEX RK2.");
-   args.AddOption(&t_final, "-tf", "--t-final",
-                  "Final time; start time is 0.");
-   args.AddOption(&dt, "-dt", "--time-step",
-                  "Time step.");
-   args.AddOption(&d_coef, "-d", "--diff-coef",
-                  "Diffusion coefficient.");
-   args.AddOption(&sigma, "-s", "--sigma",
-                  "One of the two DG penalty parameters, typically +1/-1."
-                  " See the documentation of class DGDiffusionIntegrator.");
-   args.AddOption(&kappa, "-k", "--kappa",
-                  "One of the two DG penalty parameters, should be positive."
-                  " Negative values are replaced with (order+1)^2.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
-                  "--no-visit-datafiles",
-                  "Save data files for VisIt (visit.llnl.gov) visualization.");
-   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
-                  "--ascii-datafiles",
-                  "Use binary (Sidre) or ascii format for VisIt data files.");
-   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
-                  "Visualize every n-th timestep.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(cout);
-      return 1;
-   }
-   if (kappa < 0)
-   {
-      kappa = (order+1)*(order+1);
-   }
-   args.PrintOptions(cout);
-
-   // 2. Define the ODE solver used for time integration. Several explicit
-   //    Runge-Kutta methods are available.
-   ODESolver *ode_solver = NULL;
-   switch (ode_solver_type)
-   {
-      // Explicit methods
-      case 1:  ode_solver = new ForwardEulerSolver; break;
-      case 2:  ode_solver = new RK2Solver(0.5); break; // midpoint method
-      case 3:  ode_solver = new RK3SSPSolver; break;
-      case 4:  ode_solver = new RK4Solver; break;
-      case 5:  ode_solver = new GeneralizedAlphaSolver(0.5); break;
-      // Implicit L-stable methods
-      case 11: ode_solver = new BackwardEulerSolver; break;
-      case 12: ode_solver = new SDIRK23Solver(2); break;
-      case 13: ode_solver = new SDIRK33Solver; break;
-      // Implicit A-stable methods (not L-stable)
-      case 22: ode_solver = new ImplicitMidpointSolver; break;
-      case 23: ode_solver = new SDIRK23Solver; break;
-      case 24: ode_solver = new SDIRK34Solver; break;
-      // Implicit-Explicit methods
-      case 31: ode_solver = new IMEX_BE_FE; break;
-      case 32: ode_solver = new IMEXRK2; break;
-      default:
-         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
-         return 3;
-   }
-
-   // 3. Read the serial mesh from the given mesh file on all processors. We can
-   //    handle geometrically periodic meshes in this code.
-   Mesh mesh(mesh_file, 1, 1);
-   int dim = mesh.Dimension();
-
-   // 4. Refine the mesh in serial to increase the resolution. In this example
-   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
-   //    a command-line parameter. If the mesh is of NURBS type, we convert it
-   //    to a (piecewise-polynomial) high-order mesh.
-   for (int lev = 0; lev < ref_levels; lev++)
-   {
-      mesh.UniformRefinement();
-   }
-   if (mesh.NURBSext)
-   {
-      mesh.SetCurvature(max(order, 1));
-   }
-   mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
-
-   // 5. Define the parallel discontinuous DG finite element space on the
-   //    parallel refined mesh of the given polynomial order.
-   DG_FECollection fec(order, dim);
-   FiniteElementSpace fes(&mesh, &fec);
-
-   cout << "Number of unknowns: " << fes.GetVSize() << endl;
-
-   // 6. Set up and assemble the parallel bilinear and linear forms (and the
-   //    parallel hypre matrices) corresponding to the DG discretization. The
-   //    DGTraceIntegrator involves integrals over mesh interior faces.
-   ConstantCoefficient diff_coef(d_coef);
-   VectorFunctionCoefficient velocity(dim, velocity_function);
-   FunctionCoefficient u0(u0_function);
-
-   BilinearForm m(&fes);
-   m.AddDomainIntegrator(new MassIntegrator);
-
-   BilinearForm s(&fes);
-   s.AddDomainIntegrator(new DiffusionIntegrator(diff_coef));
-   s.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma,
-                                                         kappa));
-   s.AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma, kappa));
-
-   BilinearForm k(&fes);
-   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
-   k.AddInteriorFaceIntegrator(
-      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
-   k.AddBdrFaceIntegrator(
-      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
-
-   LinearForm b(&fes);
-   b.AddBdrFaceIntegrator(
-      new DGDirichletLFIntegrator(u0, diff_coef, sigma, kappa));
-
-   int skip_zeros = 0;
-   m.Assemble(skip_zeros);
-   m.Finalize(skip_zeros);
-   s.Assemble(skip_zeros);
-   s.Finalize(skip_zeros);
-   k.Assemble(skip_zeros);
-   k.Finalize(skip_zeros);
-   b.Assemble();
-
-   // 7. Define the initial conditions, save the corresponding grid function to
-   //    a file and (optionally) save data in the VisIt format and initialize
-   //    GLVis visualization.
-   GridFunction u(&fes);
-   u.ProjectCoefficient(u0);
-
-   {
-      ofstream omesh("ex23.mesh");
-      omesh.precision(precision);
-      mesh.Print(omesh);
-      ofstream osol("ex23-init.gf");
-      osol.precision(precision);
-      u.Save(osol);
-   }
-
-   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
-   DataCollection *dc = NULL;
-   if (visit)
-   {
-      if (binary)
-      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example23", &mesh);
-#else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
-#endif
-      }
-      else
-      {
-         dc = new VisItDataCollection("Example23", &mesh);
-         dc->SetPrecision(precision);
-      }
-      dc->RegisterField("solution", &u);
-      dc->SetCycle(0);
-      dc->SetTime(0.0);
-      dc->Save();
-   }
-
-   socketstream sout;
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      sout.open(vishost, visport);
-      if (!sout)
-      {
-         cout << "Unable to connect to GLVis server at "
-              << vishost << ':' << visport << endl;
-         visualization = false;
-         cout << "GLVis visualization disabled.\n";
-      }
-      else
-      {
-         sout.precision(precision);
-         sout << "solution\n" << mesh << u;
-         sout << "pause\n";
-         sout << flush;
-         cout << "GLVis visualization paused."
-              << " Press space (in the GLVis window) to resume it.\n";
-      }
-   }
-
-   // 8. Define the time-dependent evolution operator describing the ODE
-   //    right-hand side, and perform time-integration (looping over the time
-   //    iterations, ti, with a time-step dt).
-
-   TimeDependentOperator *adv = NULL;
-   if (ode_solver_type < 10)
-   {
-      adv = new EX_Evolution(m.SpMat(), s.SpMat(), k.SpMat(), b);
-   }
-   else if (ode_solver_type < 30)
-   {
-      adv = new IM_Evolution(m.SpMat(), s.SpMat(), k.SpMat(), b);
-   }
-   else
-   {
-      adv = new IMEX_Evolution(m.SpMat(), s.SpMat(), k.SpMat(), b);
-   }
-
-   double t = 0.0;
-   adv->SetTime(t);
-   ode_solver->Init(*adv);
-
-   int n_steps = (int)ceil(t_final / dt);
-   double dt_real = t_final / n_steps;
-
-   for (int ti = 0; ti < n_steps; )
-   {
-      ode_solver->Step(u, t, dt_real);
-      ti++;
-
-      if (ti % vis_steps == 0 || ti == n_steps)
-      {
-         cout << "time step: " << ti << ", time: " << t << endl;
-
-         if (visualization)
-         {
-            sout << "solution\n" << mesh << u << flush;
-         }
-
-         if (visit)
-         {
-            dc->SetCycle(ti);
-            dc->SetTime(t);
-            dc->Save();
-         }
-      }
-   }
-
-   // 9. Save the final solution in parallel. This output can be viewed later
-   //    using GLVis: "glvis -np <np> -m ex23-mesh -g ex23-final".
-   {
-      ofstream osol("ex23-final.gf");
-      osol.precision(precision);
-      u.Save(osol);
-   }
-
-   // 10. Free the used memory.
-   delete ode_solver;
-   delete adv;
-   delete dc;
-
-   return 0;
-}
-
-
-// Implementation of class EX_Evolution
-EX_Evolution::EX_Evolution(SparseMatrix &_M, SparseMatrix &_S,
-                           SparseMatrix &_K, const Vector &_b)
-   : TimeDependentOperator(_M.Height()),
-     M(_M), S(_S), K(_K), b(_b), z(_M.Height())
-{
-   M_solver.SetPreconditioner(M_prec);
-   M_solver.SetOperator(M);
-
-   M_solver.iterative_mode = false;
-   M_solver.SetRelTol(1e-9);
-   M_solver.SetAbsTol(0.0);
-   M_solver.SetMaxIter(100);
-   M_solver.SetPrintLevel(0);
-}
-
-void EX_Evolution::Mult(const Vector &x, Vector &y) const
-{
-   // y = M^{-1} (-S x + K x + b)
-   K.Mult(x, z);
-   S.AddMult(x, z, -1.0);
-   z += b;
-   M_solver.Mult(z, y);
-}
-
-// Implementation of class IM_Evolution
-IM_Evolution::IM_Evolution(SparseMatrix &_M, SparseMatrix &_S,
-                           SparseMatrix &_K, const Vector &_b)
-   : TimeDependentOperator(_M.Height()),
-     M(_M), S(_S), K(_K), A(NULL), b(_b),
-     A_prec(NULL), A_solver(NULL), dt(-1.0), z(_M.Height())
-{
-   M_solver.SetPreconditioner(M_prec);
-   M_solver.SetOperator(M);
-
-   M_solver.iterative_mode = false;
-   M_solver.SetRelTol(1e-9);
-   M_solver.SetAbsTol(0.0);
-   M_solver.SetMaxIter(100);
-   M_solver.SetPrintLevel(0);
-}
-
-void IM_Evolution::initA(double _dt)
-{
-   if (fabs(dt - _dt) > 1e-4 * _dt)
-   {
-      delete A_solver;
-      delete A_prec;
-      delete A;
-
-      SparseMatrix * SK = Add(1.0, S, -1.0, K);
-      A = Add(1.0, M, _dt, *SK);
-      delete SK;
-      dt = _dt;
-
-      A_prec = new DSmoother(*A);
-      A_solver = new GMRESSolver;
-      A_solver->SetOperator(*A);
-      A_solver->SetPreconditioner(*A_prec);
-
-      A_solver->iterative_mode = false;
-      A_solver->SetRelTol(1e-9);
-      A_solver->SetAbsTol(0.0);
-      A_solver->SetMaxIter(100);
-      A_solver->SetPrintLevel(0);
-   }
-}
-
-void IM_Evolution::Mult(const Vector &x, Vector &y) const
-{
-   // y = M^{-1} (-S x + K x + b)
-   K.Mult(x, z);
-   S.AddMult(x, z, -1.0);
-   z += b;
-   M_solver.Mult(z, y);
-}
-
-void IM_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
-{
-   this->initA(_dt);
-
-   // y = (M + dt S - dt K)^{-1} (-S x + K x + b)
-   K.Mult(x, z);
-   S.AddMult(x, z, -1.0);
-   z += b;
-   A_solver->Mult(z, y);
-}
-
-// Implementation of class IMEX_Evolution
-IMEX_Evolution::IMEX_Evolution(SparseMatrix &_M, SparseMatrix &_S,
-                               SparseMatrix &_K, const Vector &_b)
-   : TimeDependentOperator(_M.Height()),
-     M(_M), S(_S), K(_K), A(NULL), b(_b),
-     A_prec(NULL), A_solver(NULL), dt(-1.0), z(_M.Height())
-{
-   M_solver.SetPreconditioner(M_prec);
-   M_solver.SetOperator(M);
-
-   M_solver.iterative_mode = false;
-   M_solver.SetRelTol(1e-9);
-   M_solver.SetAbsTol(0.0);
-   M_solver.SetMaxIter(100);
-   M_solver.SetPrintLevel(0);
-}
-
-void IMEX_Evolution::initA(double _dt)
-{
-   if (fabs(dt - _dt) > 1e-4 * _dt)
-   {
-      delete A_solver;
-      delete A_prec;
-      delete A;
-
-      A = Add(_dt, S, 1.0, M); // A = M + dt * S
-      dt = _dt;
-
-      A_prec = new DSmoother(*A);
-      A_solver = new CGSolver;
-      A_solver->SetOperator(*A);
-      A_solver->SetPreconditioner(*A_prec);
-
-      A_solver->iterative_mode = false;
-      A_solver->SetRelTol(1e-9);
-      A_solver->SetAbsTol(0.0);
-      A_solver->SetMaxIter(100);
-      A_solver->SetPrintLevel(0);
-   }
-}
-
-void IMEX_Evolution::Mult(const Vector &x, Vector &y) const
-{
-   // y = M^{-1} (-S x + K x + b)
-   K.Mult(x, z);
-   S.AddMult(x, z, -1.0);
-   z += b;
-   M_solver.Mult(z, y);
-}
-
-void IMEX_Evolution::ExplicitMult(const Vector &x, Vector &y) const
-{
-   // y = M^{-1} (K x + b)
-   K.Mult(x, z);
-   z += b;
-   M_solver.Mult(z, y);
-}
-
-void IMEX_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
-{
-   this->initA(_dt);
-   // y = (M + dt S)^{-1} (-S x + b)
-   S.Mult(x, z);
-   z *= -1.0;
-   z += b;
-   A_solver->Mult(z, y);
-}
-
-// Velocity coefficient
-void velocity_function(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-
-   // map to the reference [-1,1] domain
-   Vector X(dim);
-   for (int i = 0; i < dim; i++)
-   {
-      double center = (bb_min[i] + bb_max[i]) * 0.5;
-      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
-   }
-
-   switch (problem)
-   {
-      case 0:
-      {
-         // Translations in 1D, 2D, and 3D
-         switch (dim)
-         {
-            case 1: v(0) = 1.0; break;
-            case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
-            case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
-               break;
-         }
-         break;
-      }
-      case 1:
-      case 2:
-      {
-         // Clockwise rotation in 2D around the origin
-         const double w = M_PI/2;
-         switch (dim)
-         {
-            case 1: v(0) = 1.0; break;
-            case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
-            case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
-         }
-         break;
-      }
-      case 3:
-      {
-         // Clockwise twisting rotation in 2D around the origin
-         const double w = M_PI/2;
-         double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
-         d = d*d;
-         switch (dim)
-         {
-            case 1: v(0) = 1.0; break;
-            case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
-            case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
-         }
-         break;
-      }
-   }
-}
-
-// Initial condition
-double u0_function(const Vector &x)
-{
-   int dim = x.Size();
-
-   // map to the reference [-1,1] domain
-   Vector X(dim);
-   for (int i = 0; i < dim; i++)
-   {
-      double center = (bb_min[i] + bb_max[i]) * 0.5;
-      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
-   }
-
-   switch (problem)
-   {
-      case 0:
-      case 1:
-      {
-         switch (dim)
-         {
-            case 1:
-               return exp(-40.*pow(X(0)-0.5,2));
-            case 2:
-            case 3:
-            {
-               double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
-               if (dim == 3)
-               {
-                  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
-                  rx *= s;
-                  ry *= s;
-               }
-               return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
-                        erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
-            }
-         }
-      }
-      case 2:
-      {
-         double x_ = X(0), y_ = X(1), rho, phi;
-         rho = hypot(x_, y_);
-         phi = atan2(y_, x_);
-         return pow(sin(M_PI*rho),2)*sin(3*phi);
-      }
-      case 3:
-      {
-         const double f = M_PI;
-         return sin(f*X(0))*sin(f*X(1));
-      }
-   }
-   return 0.0;
-}
-
-// Inflow boundary condition (zero for the problems considered in this example)
-double inflow_function(const Vector &x)
-{
-   switch (problem)
-   {
-      case 0:
-      case 1:
-      case 2:
-      case 3: return 0.0;
-   }
-   return 0.0;
-}

examples/ex23p.cpp

@@ -1,797 +0,0 @@
-//                       MFEM Example 23 - Parallel Version
-//
-// Compile with: make ex23p
-//
-// Sample runs:
-//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -s 2 -dt 0.001 -vs 50
-//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -s 12 -dt 0.01
-//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -s 22 -dt 0.01
-//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -s 32 -dt 0.005 -vs 10
-//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 0 -dt 0.01
-//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 0 -s 32 -dt 0.01
-//    mpirun -np 4 ex23p -m ../data/periodic-hexagon.mesh -p 0 -d 0.001 -s 12 -dt 0.02
-//    mpirun -np 4 ex23p -m ../data/periodic-hexagon.mesh -p 0 -d 0.001 -s 32 -dt 0.009 -vs 10
-//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 1 -dt 0.01 -tf 9
-//    mpirun -np 4 ex23p -m ../data/periodic-hexagon.mesh -p 1 -dt 0.01 -tf 9
-//    mpirun -np 4 ex23p -m ../data/amr-quad.mesh -p 1 -dt 0.01 -tf 9 -vs 2
-//    mpirun -np 4 ex23p -m ../data/disc-nurbs.mesh -p 1 -rp 1 -dt 0.01 -tf 9
-//    mpirun -np 4 ex23p -m ../data/disc-nurbs.mesh -p 2 -rp 1 -dt 0.01 -tf 9
-//    mpirun -np 4 ex23p -m ../data/disc-nurbs.mesh -p 3 -rp 1 -dt 0.01 -tf 9 -d 0.02
-//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 3 -rp 1 -dt 0.025 -tf 9
-//    mpirun -np 4 ex23p -m ../data/periodic-cube.mesh -p 0 -o 2 -dt 0.025 -tf 8
-//
-// Description:  This example code solves the time-dependent advection-diffusion
-//               equation
-//               du/dt - div(D grad(u)) + v.grad(u) = 0, where
-//               D is a diffusion coefficient,
-//               v is a given fluid velocity, and
-//               u0(x)=u(0,x) is a given initial condition.
-//
-//               The example demonstrates the use of Discontinuous Galerkin (DG)
-//               bilinear forms in MFEM (face integrators), the use of explicit,
-//               implicit, and implicit-explicit ODE time integrators, the
-//               definition of periodic boundary conditions through periodic
-//               meshes, as well as the use of GLVis for persistent
-//               visualization of a time-evolving solution. The saving of
-//               time-dependent data files for external visualization with
-//               VisIt (visit.llnl.gov) is also illustrated.
-//
-//               This example is a merger of examples 9 and 14.
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-// Choice for the problem setup. The fluid velocity, initial condition and
-// boundary condition are chosen based on this parameter.
-int problem;
-
-// Velocity coefficient
-void velocity_function(const Vector &x, Vector &v);
-
-// Initial condition
-double u0_function(const Vector &x);
-
-// Mesh bounding box
-Vector bb_min, bb_max;
-
-/** A time-dependent operator for the right-hand side of the ODE for use with
-    explicit ODE solvers. The DG weak form of du/dt = div(D grad(u))-v.grad(u) is
-    M du/dt = - S u + K u + b, where M, S, and K are the mass,
-    stiffness, and advection matrices, and b describes sources and the flow on
-    the boundary.
-    This can be written as a general ODE,
-    du/dt = M^{-1} (-S u + K u + b), and this class is used to compute the RHS
-    and perform the solve for du/dt. */
-class EX_Evolution : public TimeDependentOperator
-{
-private:
-   HypreParMatrix &M, &S, &K;
-   const Vector &b;
-
-   HypreSmoother M_prec;
-   CGSolver M_solver;
-
-   mutable Vector z;
-
-   void initA(double dt);
-
-public:
-   EX_Evolution(HypreParMatrix &_M, HypreParMatrix &_S, HypreParMatrix &_K,
-                const Vector &_b);
-
-   virtual void Mult(const Vector &x, Vector &y) const;
-
-   virtual ~EX_Evolution() {}
-};
-
-/** A time-dependent operator for the right-hand side of the ODE for use with
-    implicit ODE solvers. The DG weak form of du/dt = div(D grad(u))-v.grad(u) is
-    [M + dt (S - K)] du/dt = - S u + K u + b, where M, S, and K are the mass,
-    stiffness, and advection matrices, and b describes sources and the flow on
-    the boundary.
-    This can be written as a general ODE,
-    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
-    used to perform the fully implicit solve for du/dt. */
-class IM_Evolution : public TimeDependentOperator
-{
-private:
-   HypreParMatrix &M, &S, &K;
-   HypreParMatrix *A;
-   const Vector &b;
-
-   HypreSmoother M_prec;
-   CGSolver M_solver;
-
-   HypreBoomerAMG *A_prec;
-   GMRESSolver *A_solver;
-   double dt;
-
-   mutable Vector z;
-
-   void initA(double dt);
-
-public:
-   IM_Evolution(HypreParMatrix &_M, HypreParMatrix &_S, HypreParMatrix &_K,
-                const Vector &_b);
-
-   virtual void Mult(const Vector &x, Vector &y) const;
-
-   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
-
-   virtual ~IM_Evolution() { delete A_solver; delete A_prec; delete A; }
-};
-
-/** A time-dependent operator for the right-hand side of the ODE for use with
-    IMEX (Implicit-Explicit) ODE solvers. The DG weak form of
-    du/dt = div(D grad(u))-v.grad(u) is
-    [M + dt S] du/dt = - S u + K u + b, where M, S, and K are the mass,
-    stiffness, and advection matrices, and b describes sources and the flow on
-    the boundary.
-    This can be written as a general ODE,
-    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
-    used to perform the implicit or explicit solve for du/dt. */
-class IMEX_Evolution : public TimeDependentOperator
-{
-private:
-   HypreParMatrix &M, &S, &K;
-   HypreParMatrix *A;
-   const Vector &b;
-
-   HypreSmoother M_prec;
-   CGSolver M_solver;
-
-   HypreBoomerAMG *A_prec;
-   CGSolver *A_solver;
-   double dt;
-
-   mutable Vector z;
-
-   void initA(double dt);
-
-public:
-   IMEX_Evolution(HypreParMatrix &_M, HypreParMatrix &_S, HypreParMatrix &_K,
-                  const Vector &_b);
-
-   virtual void ExplicitMult(const Vector &x, Vector &y) const;
-   virtual void Mult(const Vector &x, Vector &y) const;
-
-   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
-
-   virtual ~IMEX_Evolution() { delete A_solver; delete A_prec; delete A; }
-};
-
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI.
-   int num_procs, myid;
-   MPI_Init(&argc, &argv);
-   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
-   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
-
-   // 2. Parse command-line options.
-   problem = 0;
-   const char *mesh_file = "../data/periodic-hexagon.mesh";
-   int ser_ref_levels = 2;
-   int par_ref_levels = 0;
-   int order = 3;
-   int ode_solver_type = 12;
-   double t_final = 10.0;
-   double d_coef = 0.01;
-   double dt = 0.01;
-   double sigma = -1.0;
-   double kappa = -1.0;
-   bool visualization = true;
-   bool visit = false;
-   bool binary = false;
-   int vis_steps = 5;
-
-   int precision = 8;
-   cout.precision(precision);
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&problem, "-p", "--problem",
-                  "Problem setup to use. See options in velocity_function().");
-   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   args.AddOption(&order, "-o", "--order",
-                  "Order (degree) of the finite elements.");
-   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
-                  "ODE solver: 1 - Forward Euler, 2 - RK2, 3 - RK3 SSP,"
-                  " 4 - RK4, 5 - Generalized Alpha,\n\t"
-                  "11 - Backward Euler, 12 - SDIRK2, 13 - SDIRK3,\n\t"
-                  "22 - Implicit Midpoint, 23 SDIRK23, 24 - SDIRK34,\n\t"
-                  "31 - IMEX BE/FE, 32 - IMEX RK2.");
-   args.AddOption(&t_final, "-tf", "--t-final",
-                  "Final time; start time is 0.");
-   args.AddOption(&dt, "-dt", "--time-step",
-                  "Time step.");
-   args.AddOption(&d_coef, "-d", "--diff-coef",
-                  "Diffusion coefficient.");
-   args.AddOption(&sigma, "-s", "--sigma",
-                  "One of the two DG penalty parameters, typically +1/-1."
-                  " See the documentation of class DGDiffusionIntegrator.");
-   args.AddOption(&kappa, "-k", "--kappa",
-                  "One of the two DG penalty parameters, should be positive."
-                  " Negative values are replaced with (order+1)^2.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
-                  "--no-visit-datafiles",
-                  "Save data files for VisIt (visit.llnl.gov) visualization.");
-   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
-                  "--ascii-datafiles",
-                  "Use binary (Sidre) or ascii format for VisIt data files.");
-   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
-                  "Visualize every n-th timestep.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      MPI_Finalize();
-      return 1;
-   }
-   if (kappa < 0)
-   {
-      kappa = (order+1)*(order+1);
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-
-   // 3. Define the ODE solver used for time integration. Several explicit,
-   //    implicitit, and implicit-explicit Runge-Kutta methods are available.
-   ODESolver *ode_solver = NULL;
-
-   switch (ode_solver_type)
-   {
-      // Explicit methods
-      case 1:  ode_solver = new ForwardEulerSolver; break;
-      case 2:  ode_solver = new RK2Solver(0.5); break; // midpoint method
-      case 3:  ode_solver = new RK3SSPSolver; break;
-      case 4:  ode_solver = new RK4Solver; break;
-      case 5:  ode_solver = new GeneralizedAlphaSolver(0.5); break;
-      // Implicit L-stable methods
-      case 11: ode_solver = new BackwardEulerSolver; break;
-      case 12: ode_solver = new SDIRK23Solver(2); break;
-      case 13: ode_solver = new SDIRK33Solver; break;
-      // Implicit A-stable methods (not L-stable)
-      case 22: ode_solver = new ImplicitMidpointSolver; break;
-      case 23: ode_solver = new SDIRK23Solver; break;
-      case 24: ode_solver = new SDIRK34Solver; break;
-      // Implicit-Explicit methods
-      case 31: ode_solver = new IMEX_BE_FE; break;
-      case 32: ode_solver = new IMEXRK2; break;
-      default:
-         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
-         return 3;
-   }
-
-   // 4. Read the serial mesh from the given mesh file on all processors. We can
-   //    handle geometrically periodic meshes in this code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   int dim = mesh->Dimension();
-
-   // 5. Refine the mesh in serial to increase the resolution. In this example
-   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
-   //    a command-line parameter. If the mesh is of NURBS type, we convert it
-   //    to a (piecewise-polynomial) high-order mesh.
-   for (int lev = 0; lev < ser_ref_levels; lev++)
-   {
-      mesh->UniformRefinement();
-   }
-   if (mesh->NURBSext)
-   {
-      mesh->SetCurvature(max(order, 1));
-   }
-   mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
-
-   // 6. Define the 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 = new ParMesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
-   for (int lev = 0; lev < par_ref_levels; lev++)
-   {
-      pmesh->UniformRefinement();
-   }
-
-   // 7. Define the parallel discontinuous DG finite element space on the
-   //    parallel refined mesh of the given polynomial order.
-   DG_FECollection fec(order, dim);
-   ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
-
-   HYPRE_Int global_vSize = fes->GlobalTrueVSize();
-   if (myid == 0)
-   {
-      cout << "Number of unknowns: " << global_vSize << endl;
-   }
-
-   // 8. Set up and assemble the parallel bilinear and linear forms (and the
-   //    parallel hypre matrices) corresponding to the DG discretization. The
-   //    DGTraceIntegrator involves integrals over mesh interior faces.
-   ConstantCoefficient diff_coef(d_coef);
-   VectorFunctionCoefficient velocity(dim, velocity_function);
-   FunctionCoefficient u0(u0_function);
-
-   ParBilinearForm *m = new ParBilinearForm(fes);
-   m->AddDomainIntegrator(new MassIntegrator);
-
-   ParBilinearForm *s = new ParBilinearForm(fes);
-   s->AddDomainIntegrator(new DiffusionIntegrator(diff_coef));
-   s->AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma,
-                                                          kappa));
-   s->AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma, kappa));
-
-   ParBilinearForm *k = new ParBilinearForm(fes);
-   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
-   k->AddInteriorFaceIntegrator(
-      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
-   k->AddBdrFaceIntegrator(
-      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
-
-   ParLinearForm *b = new ParLinearForm(fes);
-   b->AddBdrFaceIntegrator(
-      new DGDirichletLFIntegrator(u0, diff_coef, sigma, kappa));
-
-   int skip_zeros = 0;
-   m->Assemble(skip_zeros);
-   m->Finalize(skip_zeros);
-   s->Assemble(skip_zeros);
-   s->Finalize(skip_zeros);
-   k->Assemble(skip_zeros);
-   k->Finalize(skip_zeros);
-   b->Assemble();
-
-   HypreParMatrix *M = m->ParallelAssemble();
-   HypreParMatrix *S = s->ParallelAssemble();
-   HypreParMatrix *K = k->ParallelAssemble();
-   HypreParVector *B = b->ParallelAssemble();
-
-   // 9. Define the initial conditions, save the corresponding grid function to
-   //    a file and (optionally) save data in the VisIt format and initialize
-   //    GLVis visualization.
-   ParGridFunction *u = new ParGridFunction(fes);
-   u->ProjectCoefficient(u0);
-   HypreParVector *U = u->GetTrueDofs();
-
-   {
-      ostringstream mesh_name, sol_name;
-      mesh_name << "ex23-mesh." << setfill('0') << setw(6) << myid;
-      sol_name << "ex23-init." << setfill('0') << setw(6) << myid;
-      ofstream omesh(mesh_name.str().c_str());
-      omesh.precision(precision);
-      pmesh->Print(omesh);
-      ofstream osol(sol_name.str().c_str());
-      osol.precision(precision);
-      u->Save(osol);
-   }
-
-   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
-   DataCollection *dc = NULL;
-   if (visit)
-   {
-      if (binary)
-      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example23-Parallel", pmesh);
-#else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
-#endif
-      }
-      else
-      {
-         dc = new VisItDataCollection("Example23-Parallel", pmesh);
-         dc->SetPrecision(precision);
-         // To save the mesh using MFEM's parallel mesh format:
-         // dc->SetFormat(DataCollection::PARALLEL_FORMAT);
-      }
-      dc->RegisterField("solution", u);
-      dc->SetCycle(0);
-      dc->SetTime(0.0);
-      dc->Save();
-   }
-
-   socketstream sout;
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      sout.open(vishost, visport);
-      if (!sout)
-      {
-         if (myid == 0)
-            cout << "Unable to connect to GLVis server at "
-                 << vishost << ':' << visport << endl;
-         visualization = false;
-         if (myid == 0)
-         {
-            cout << "GLVis visualization disabled.\n";
-         }
-      }
-      else
-      {
-         sout << "parallel " << num_procs << " " << myid << "\n";
-         sout.precision(precision);
-         sout << "solution\n" << *pmesh << *u;
-         sout << "pause\n";
-         sout << flush;
-         if (myid == 0)
-            cout << "GLVis visualization paused."
-                 << " Press space (in the GLVis window) to resume it.\n";
-      }
-   }
-
-   // 10. Define the time-dependent evolution operator describing the ODE
-   //     right-hand side, and perform time-integration (looping over the time
-   //     iterations, ti, with a time-step dt).
-
-   TimeDependentOperator *adv = NULL;
-   if (ode_solver_type < 10)
-   {
-      adv = new EX_Evolution(*M, *S, *K, *B);
-   }
-   else if (ode_solver_type < 30)
-   {
-      adv = new IM_Evolution(*M, *S, *K, *B);
-   }
-   else
-   {
-      adv = new IMEX_Evolution(*M, *S, *K, *B);
-   }
-
-   double t = 0.0;
-   adv->SetTime(t);
-   ode_solver->Init(*adv);
-
-   int n_steps = (int)ceil(t_final / dt);
-   double dt_real = t_final / n_steps;
-
-   for (int ti = 0; ti < n_steps; )
-   {
-      ode_solver->Step(*U, t, dt_real);
-      ti++;
-
-      if (ti % vis_steps == 0 || ti == n_steps)
-      {
-         if (myid == 0)
-         {
-            cout << "time step: " << ti << ", time: " << t << endl;
-         }
-
-         // 11. Extract the parallel grid function corresponding to the finite
-         //     element approximation U (the local solution on each processor).
-         *u = *U;
-
-         if (visualization)
-         {
-            sout << "parallel " << num_procs << " " << myid << "\n";
-            sout << "solution\n" << *pmesh << *u << flush;
-         }
-
-         if (visit)
-         {
-            dc->SetCycle(ti);
-            dc->SetTime(t);
-            dc->Save();
-         }
-      }
-   }
-
-   // 12. Save the final solution in parallel. This output can be viewed later
-   //     using GLVis: "glvis -np <np> -m ex23-mesh -g ex23-final".
-   {
-      *u = *U;
-      ostringstream sol_name;
-      sol_name << "ex23-final." << setfill('0') << setw(6) << myid;
-      ofstream osol(sol_name.str().c_str());
-      osol.precision(precision);
-      u->Save(osol);
-   }
-
-   // 13. Free the used memory.
-   delete U;
-   delete u;
-   delete B;
-   delete b;
-   delete K;
-   delete k;
-   delete S;
-   delete s;
-   delete M;
-   delete m;
-   delete fes;
-   delete pmesh;
-   delete ode_solver;
-   delete adv;
-   delete dc;
-
-   MPI_Finalize();
-   return 0;
-}
-
-
-// Implementation of class EX_Evolution
-EX_Evolution::EX_Evolution(HypreParMatrix &_M, HypreParMatrix &_S,
-                           HypreParMatrix &_K, const Vector &_b)
-   : TimeDependentOperator(_M.Height()),
-     M(_M), S(_S), K(_K), b(_b),
-     M_prec(M), M_solver(M.GetComm()), z(M.Height())
-{
-   M_prec.SetType(HypreSmoother::Jacobi);
-   M_solver.SetPreconditioner(M_prec);
-   M_solver.SetOperator(M);
-
-   M_solver.iterative_mode = false;
-   M_solver.SetRelTol(1e-9);
-   M_solver.SetAbsTol(0.0);
-   M_solver.SetMaxIter(100);
-   M_solver.SetPrintLevel(0);
-}
-
-void EX_Evolution::Mult(const Vector &x, Vector &y) const
-{
-   // y = M^{-1} (-S x + K x + b)
-   S.Mult(-1.0, x, 0.0, z);
-   K.Mult(1.0, x, 1.0, z);
-   z += b;
-   M_solver.Mult(z, y);
-}
-
-// Implementation of class IM_Evolution
-IM_Evolution::IM_Evolution(HypreParMatrix &_M, HypreParMatrix &_S,
-                           HypreParMatrix &_K, const Vector &_b)
-   : TimeDependentOperator(_M.Height()),
-     M(_M), S(_S), K(_K), A(NULL), b(_b),
-     M_prec(M), M_solver(M.GetComm()),
-     A_prec(NULL), A_solver(NULL), dt(-1.0), z(M.Height())
-{
-   M_prec.SetType(HypreSmoother::Jacobi);
-   M_solver.SetPreconditioner(M_prec);
-   M_solver.SetOperator(M);
-
-   M_solver.iterative_mode = false;
-   M_solver.SetRelTol(1e-9);
-   M_solver.SetAbsTol(0.0);
-   M_solver.SetMaxIter(100);
-   M_solver.SetPrintLevel(0);
-}
-
-void IM_Evolution::initA(double _dt)
-{
-   if (fabs(dt - _dt) > 1e-4 * _dt)
-   {
-      delete A_solver;
-      delete A_prec;
-      delete A;
-
-      HypreParMatrix * SK = Add(1.0, S, -1.0, K); // SK = S - K
-      A = Add(_dt, *SK, 1.0, M);                  // A = M + dt * (S - K)
-      delete SK;
-      dt = _dt;
-
-      A_prec = new HypreBoomerAMG(*A);
-      A_solver = new GMRESSolver(A->GetComm());
-      A_solver->SetOperator(*A);
-      A_solver->SetPreconditioner(*A_prec);
-
-      A_solver->iterative_mode = false;
-      A_solver->SetRelTol(1e-9);
-      A_solver->SetAbsTol(0.0);
-      A_solver->SetMaxIter(100);
-      A_solver->SetPrintLevel(0);
-   }
-}
-
-void IM_Evolution::Mult(const Vector &x, Vector &y) const
-{
-   // y = M^{-1} (-S x + K x + b)
-   S.Mult(-1.0, x, 0.0, z);
-   K.Mult(1.0, x, 1.0, z);
-   z += b;
-   M_solver.Mult(z, y);
-}
-
-void IM_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
-{
-   this->initA(_dt);
-
-   // y = (M + dt S - dt K)^{-1} (-S x + K x + b)
-   S.Mult(-1.0, x, 0.0, z);
-   K.Mult(1.0, x, 1.0, z);
-   z += b;
-   A_solver->Mult(z, y);
-}
-
-// Implementation of class IMEX_Evolution
-IMEX_Evolution::IMEX_Evolution(HypreParMatrix &_M, HypreParMatrix &_S,
-                               HypreParMatrix &_K, const Vector &_b)
-   : TimeDependentOperator(_M.Height()),
-     M(_M), S(_S), K(_K), A(NULL), b(_b),
-     M_prec(M), M_solver(M.GetComm()),
-     A_prec(NULL), A_solver(NULL), dt(-1.0), z(M.Height())
-{
-   M_prec.SetType(HypreSmoother::Jacobi);
-   M_solver.SetPreconditioner(M_prec);
-   M_solver.SetOperator(M);
-
-   M_solver.iterative_mode = false;
-   M_solver.SetRelTol(1e-9);
-   M_solver.SetAbsTol(0.0);
-   M_solver.SetMaxIter(100);
-   M_solver.SetPrintLevel(0);
-}
-
-void IMEX_Evolution::initA(double _dt)
-{
-   if (fabs(dt - _dt) > 1e-4 * _dt)
-   {
-      delete A_solver;
-      delete A_prec;
-      delete A;
-
-      A = Add(_dt, S, 1.0, M); // A = M + dt * S
-      dt = _dt;
-
-      A_prec = new HypreBoomerAMG(*A);
-      A_solver = new CGSolver(A->GetComm());
-      A_solver->SetOperator(*A);
-      A_solver->SetPreconditioner(*A_prec);
-
-      A_solver->iterative_mode = false;
-      A_solver->SetRelTol(1e-9);
-      A_solver->SetAbsTol(0.0);
-      A_solver->SetMaxIter(100);
-      A_solver->SetPrintLevel(0);
-   }
-}
-
-void IMEX_Evolution::Mult(const Vector &x, Vector &y) const
-{
-   // y = M^{-1} (-S x + K x + b)
-   S.Mult(-1.0, x, 0.0, z);
-   K.Mult(1.0, x, 1.0, z);
-   z += b;
-   M_solver.Mult(z, y);
-}
-
-void IMEX_Evolution::ExplicitMult(const Vector &x, Vector &y) const
-{
-   // y = M^{-1} (K x + b)
-   K.Mult(1.0, x, 0.0, z);
-   z += b;
-   M_solver.Mult(z, y);
-}
-
-void IMEX_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
-{
-   this->initA(_dt);
-   // y = (M + dt S)^{-1} (-S x + b)
-   S.Mult(-1.0, x, 0.0, z);
-   z += b;
-   A_solver->Mult(z, y);
-}
-
-// Velocity coefficient
-void velocity_function(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-
-   // map to the reference [-1,1] domain
-   Vector X(dim);
-   for (int i = 0; i < dim; i++)
-   {
-      double center = (bb_min[i] + bb_max[i]) * 0.5;
-      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
-   }
-
-   switch (problem)
-   {
-      case 0:
-      {
-         // Translations in 1D, 2D, and 3D
-         switch (dim)
-         {
-            case 1: v(0) = 1.0; break;
-            case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
-            case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
-               break;
-         }
-         break;
-      }
-      case 1:
-      case 2:
-      {
-         // Clockwise rotation in 2D around the origin
-         const double w = M_PI/2;
-         switch (dim)
-         {
-            case 1: v(0) = 1.0; break;
-            case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
-            case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
-         }
-         break;
-      }
-      case 3:
-      {
-         // Clockwise twisting rotation in 2D around the origin
-         const double w = M_PI/2;
-         double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
-         d = d*d;
-         switch (dim)
-         {
-            case 1: v(0) = 1.0; break;
-            case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
-            case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
-         }
-         break;
-      }
-   }
-}
-
-// Initial condition
-double u0_function(const Vector &x)
-{
-   int dim = x.Size();
-
-   // map to the reference [-1,1] domain
-   Vector X(dim);
-   for (int i = 0; i < dim; i++)
-   {
-      double center = (bb_min[i] + bb_max[i]) * 0.5;
-      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
-   }
-
-   switch (problem)
-   {
-      case 0:
-      case 1:
-      {
-         switch (dim)
-         {
-            case 1:
-               return exp(-40.*pow(X(0)-0.5,2));
-            case 2:
-            case 3:
-            {
-               double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
-               if (dim == 3)
-               {
-                  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
-                  rx *= s;
-                  ry *= s;
-               }
-               return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
-                        erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
-            }
-         }
-      }
-      case 2:
-      {
-         double x_ = X(0), y_ = X(1), rho, phi;
-         rho = hypot(x_, y_);
-         phi = atan2(y_, x_);
-         return pow(sin(M_PI*rho),2)*sin(3*phi);
-      }
-      case 3:
-      {
-         const double f = M_PI;
-         return sin(f*X(0))*sin(f*X(1));
-      }
-   }
-   return 0.0;
-}

examples/makefile

@@ -22,9 +22,9 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES = ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 ex17\
-  ex18 ex19 ex20 ex22 ex23
+  ex18 ex19 ex20 ex22
 PAR_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p ex12p\
- ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex22p ex23p
+ ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex22p
 
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)
@@ -117,7 +117,7 @@ clean-build:
 
 clean-exec:
 	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh
-	@rm -rf Example5* Example9* Example15* Example16* Example23*
+	@rm -rf Example5* Example9* Example15* Example16*
 	@rm -f sphere_refined.* sol.* sol_u.* sol_p.*
 	@rm -f ex9.mesh ex9-mesh.* ex9-init.* ex9-final.*
 	@rm -f deformed.* velocity.* elastic_energy.* mode_*
@@ -126,4 +126,3 @@ clean-exec:
 	@rm -f deformation.* pressure.*
 	@rm -f ex20.dat ex20p_?????.dat gnuplot_ex20.inp gnuplot_ex20p.inp
 	@rm -f ex22*.mesh ex22*.sol ex22p_*.*
-	@rm -f ex23.mesh ex23-mesh.* ex23-init.* ex23-final.*

linalg/hypre_parcsr.cpp

@@ -1520,14 +1520,10 @@ hypre_ParCSRMatrixSum(hypre_ParCSRMatrix *A,
    hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
    hypre_CSRMatrix *B_diag = hypre_ParCSRMatrixDiag(B);
    hypre_CSRMatrix *B_offd = hypre_ParCSRMatrixOffd(B);
-   HYPRE_Int ncols_B_offd  = hypre_CSRMatrixNumCols(B_offd);
    HYPRE_Int error;
 
    error = hypre_CSRMatrixSum(A_diag, beta, B_diag);
-   if (ncols_B_offd > 0) /* treat B_offd as zero if it has no columns */
-   {
-      error = error ? error : hypre_CSRMatrixSum(A_offd, beta, B_offd);
-   }
+   error = error ? error : hypre_CSRMatrixSum(A_offd, beta, B_offd);
 
    return error;
 }

linalg/ode.cpp

@@ -561,102 +561,6 @@ void GeneralizedAlphaSolver::Step(Vector &x, double &t, double &dt)
 }
 
 
-void IMEX_BE_FE::Init(TimeDependentOperator &_f)
-{
-   ODESolver::Init(_f);
-   k_imp.SetSize(f->Width());
-   y.SetSize(f->Width());
-   k_exp.SetSize(f->Width());
-}
-
-void IMEX_BE_FE::Step(Vector &x, double &t, double &dt)
-{
-   f->ExplicitMult(x, k_exp);
-   add(x, dt, k_exp, y);
-
-   f->SetTime(t + dt);
-   f->ImplicitSolve(dt, y, k_imp);
-
-   x.Add(dt, k_exp);
-   x.Add(dt, k_imp);
-   t += dt;
-}
-
-
-void IMEXRK2::Init(TimeDependentOperator &_f)
-{
-   ODESolver::Init(_f);
-   f = ODESolver::f;
-   k_imp.SetSize(f->Width());
-   k_exp.SetSize(f->Width());
-   y.SetSize(f->Width());
-   z.SetSize(f->Width());
-}
-
-void IMEXRK2::Step(Vector &x, double &t, double &dt)
-{
-   double gamma = 1 - sqrt(2)/2;
-   double delta = -2*sqrt(2)/3;
-
-   // The method is given by
-   // k1_exp = f(u)
-   // k1_imp = g(u + gamma*dt*k1_exp + gamma*dt*k1_imp)
-   // k2_exp = f(u + gamma*dt*k1_exp + gamma*dt*k1_imp)
-   // k2_imp = g(u + delta*dt*k1_exp + (1-gamma)*dt*k1_imp
-   //              + (1-delta)*dt*k2_exp + gamma*dt*k2_imp)
-   // k3_exp = f(u + delta*dt*k1_exp + (1-gamma)*dt*k1_imp
-   //              + (1-delta)*dt*k2_exp + gamma*dt*k2_imp)
-   // u_new = u + dt*((1-gamma)*k1_imp + (1-gamma)*k2_exp
-   //                 + gamma*k2_imp + gamma*k3_exp)
-
-   // Take first explicit step
-   // k1_exp = f(u)
-   f->ExplicitMult(x, k_exp);
-   // b corresponding to this stage is zero, so don't add to solution
-
-   // Solve first implicit step
-   // y = u + gamma*dt*k1_exp
-   add(x, gamma*dt, k_exp, y);
-   // Solve x1_imp = g(u + gamma*dt*k1_exp + gamma*dt*k1_imp)
-   f->SetTime(t + gamma*dt);
-   f->ImplicitSolve(gamma*dt, y, k_imp);
-   // x = u + (1-gamma)*dt*k1_imp
-   x.Add((1-gamma)*dt, k_imp);
-
-   // Begin setting up rhs for second solve
-   // z = u + (1-gamma)*dt*k_imp + delta*dt*k_exp
-   add(x, delta*dt, k_exp, z);
-
-   // Take second explicit step
-   // y = x + gamma*dt*k1_exp + gamma*dt*k1_imp
-   y.Add(gamma*dt, k_imp);
-   // k2_exp = f(x + gamma*dt*k1_exp + gamma*dt*k1_imp)
-   f->ExplicitMult(y, k_exp);
-   // x = u + (1-gamma)*dt*k1_imp + (1-gamma)*dt*k2_exp
-   x.Add((1-gamma)*dt, k_exp);
-
-   // Finish formoing rhs
-   // z = x + (1-gamma)*dt*k1_imp + delta*dt*k1_exp + (1-delta)*dt*k2_exp
-   z.Add((1-delta)*dt, k_exp);
-
-   // Solve second implicit step for k2_imp
-   f->SetTime(t + dt);
-   f->ImplicitSolve(gamma*dt, z, k_imp);
-   // x = u + (1-gamma)*dt*k1_imp + (1-gamma)*dt*k2_exp + gamma*dt*k2_imp
-   x.Add(gamma*dt, k_imp);
-
-   // Take final explicit step for k3_exp
-   z.Add(gamma*dt, k_imp);
-   f->ExplicitMult(z, k_exp);
-
-   // x = u + (1-gamma)*dt*k1_imp + (1-gamma)*dt*k2_exp + gamma*dt*k2_imp
-   //       + gamma*dt*k3_exp
-   x.Add(gamma*dt, k_exp);
-
-   t += dt;
-}
-
-
 void
 SIASolver::Init(Operator &P, TimeDependentOperator & F)
 {

linalg/ode.hpp

@@ -305,33 +305,6 @@ public:
 };
 
 
-/// IMEX Backward-Forward Euler ODE solver
-class IMEX_BE_FE : public ODESolver
-{
-protected:
-   Vector k_exp, k_imp, y;
-
-public:
-   virtual void Init(TimeDependentOperator &_f);
-
-   virtual void Step(Vector &x, double &t, double &dt);
-};
-
-/** Second-order IMEX (2,3,2) method, from "Implicit-explicit Runge-Kutta
-    methods for time-dependent partial differential equations" by Ascher, Ruuth
-    and Spiteri, Applied Numerical Mathematics (1997). */
-class IMEXRK2 : public ODESolver
-{
-protected:
-   Vector k_exp, k_imp, y, z;
-
-public:
-   virtual void Init(TimeDependentOperator &_f);
-
-   virtual void Step(Vector &x, double &t, double &dt);
-};
-
-
 /// The SIASolver class is based on the Symplectic Integration Algorithm
 /// described in "A Symplectic Integration Algorithm for Separable Hamiltonian
 /// Functions" by J. Candy and W. Rozmus, Journal of Computational Physics,

miniapps/common/pfem_extras.cpp

@@ -94,13 +94,13 @@ ParDiscreteDivOperator::ParDiscreteDivOperator(ParFiniteElementSpace *dfes,
    this->AddDomainInterpolator(new DivergenceInterpolator);
 }
 
-IrrotationalProjector
-::IrrotationalProjector(ParFiniteElementSpace   & H1FESpace,
-                        ParFiniteElementSpace   & HCurlFESpace,
-                        const int               & irOrder,
-                        ParBilinearForm         * s0,
-                        ParMixedBilinearForm    * weakDiv,
-                        ParDiscreteGradOperator * grad)
+IrrotationalNDProjector
+::IrrotationalNDProjector(ParFiniteElementSpace   & H1FESpace,
+                          ParFiniteElementSpace   & HCurlFESpace,
+                          const int               & irOrder,
+                          ParBilinearForm         * s0,
+                          ParMixedBilinearForm    * weakDiv,
+                          ParDiscreteGradOperator * grad)
    : H1FESpace_(&H1FESpace),
      HCurlFESpace_(&HCurlFESpace),
      s0_(s0),
@@ -115,10 +115,14 @@ IrrotationalProjector
      ownsWeakDiv_(weakDiv == NULL),
      ownsGrad_(grad == NULL)
 {
+   /*
    ess_bdr_.SetSize(H1FESpace_->GetParMesh()->bdr_attributes.Max());
    ess_bdr_ = 1;
    H1FESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
-
+   */
+   ess_bdr_tdofs_.SetSize(1);
+   ess_bdr_tdofs_ = 0;
+   
    int geom = H1FESpace_->GetFE(0)->GetGeomType();
    const IntegrationRule * ir = &IntRules.Get(geom, irOrder);
 
@@ -152,7 +156,7 @@ IrrotationalProjector
    xDiv_ = new ParGridFunction(H1FESpace_);
 }
 
-IrrotationalProjector::~IrrotationalProjector()
+IrrotationalNDProjector::~IrrotationalNDProjector()
 {
    delete psi_;
    delete xDiv_;
@@ -167,7 +171,7 @@ IrrotationalProjector::~IrrotationalProjector()
 }
 
 void
-IrrotationalProjector::InitSolver() const
+IrrotationalNDProjector::InitSolver() const
 {
    delete pcg_;
    delete amg_;
@@ -182,7 +186,7 @@ IrrotationalProjector::InitSolver() const
 }
 
 void
-IrrotationalProjector::Mult(const Vector &x, Vector &y) const
+IrrotationalNDProjector::Mult(const Vector &x, Vector &y) const
 {
    // Compute the divergence of x
    weakDiv_->Mult(x,*xDiv_); *xDiv_ *= -1.0;
@@ -203,7 +207,7 @@ IrrotationalProjector::Mult(const Vector &x, Vector &y) const
 }
 
 void
-IrrotationalProjector::Update()
+IrrotationalNDProjector::Update()
 {
    delete pcg_; pcg_ = NULL;
    delete amg_; amg_ = NULL;
@@ -234,31 +238,213 @@ IrrotationalProjector::Update()
    H1FESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
 }
 
-DivergenceFreeProjector
-::DivergenceFreeProjector(ParFiniteElementSpace   & H1FESpace,
-                          ParFiniteElementSpace   & HCurlFESpace,
-                          const int               & irOrder,
-                          ParBilinearForm         * s0,
-                          ParMixedBilinearForm    * weakDiv,
-                          ParDiscreteGradOperator * grad)
-   : IrrotationalProjector(H1FESpace,HCurlFESpace, irOrder, s0, weakDiv, grad)
+DivergenceFreeNDProjector
+::DivergenceFreeNDProjector(ParFiniteElementSpace   & H1FESpace,
+                            ParFiniteElementSpace   & HCurlFESpace,
+                            const int               & irOrder,
+                            ParBilinearForm         * s0,
+                            ParMixedBilinearForm    * weakDiv,
+                            ParDiscreteGradOperator * grad)
+   : IrrotationalNDProjector(H1FESpace,HCurlFESpace, irOrder, s0, weakDiv, grad)
 {}
 
-DivergenceFreeProjector::~DivergenceFreeProjector()
+DivergenceFreeNDProjector::~DivergenceFreeNDProjector()
 {}
 
 void
-DivergenceFreeProjector::Mult(const Vector &x, Vector &y) const
+DivergenceFreeNDProjector::Mult(const Vector &x, Vector &y) const
 {
-   this->IrrotationalProjector::Mult(x, y);
+   this->IrrotationalNDProjector::Mult(x, y);
    y  -= x;
    y *= -1.0;
 }
 
 void
-DivergenceFreeProjector::Update()
+DivergenceFreeNDProjector::Update()
 {
-   this->IrrotationalProjector::Update();
+   this->IrrotationalNDProjector::Update();
+}
+
+DivergenceFreeRTProjector
+::DivergenceFreeRTProjector(ParFiniteElementSpace   & HCurlFESpace,
+                            ParFiniteElementSpace   & HDivFESpace,
+                            const int               & irOrder,
+                            ParBilinearForm         * s1,
+                            ParMixedBilinearForm    * weakCurl,
+                            ParDiscreteCurlOperator * curl)
+   : HCurlFESpace_(&HCurlFESpace),
+     HDivFESpace_(&HDivFESpace),
+     s1_(s1),
+     weakCurl_(weakCurl),
+     curl_(curl),
+     psi_(NULL),
+     xCurl_(NULL),
+     S1_(NULL),
+     pc_(NULL),
+     pcg_(NULL),
+     dim_(HCurlFESpace_->GetFE(0)->GetDim()),
+     ownsS1_(s1 == NULL),
+     ownsWeakCurl_(weakCurl == NULL),
+     ownsCurl_(curl == NULL)
+{
+   ess_bdr_.SetSize(HCurlFESpace_->GetParMesh()->bdr_attributes.Max());
+   ess_bdr_ = 1;
+   HCurlFESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
+
+   int geom = HCurlFESpace_->GetFE(0)->GetGeomType();
+   const IntegrationRule * ir = &IntRules.Get(geom, irOrder);
+
+   if ( s1 == NULL )
+   {
+      s1_ = new ParBilinearForm(HCurlFESpace_);
+      BilinearFormIntegrator * ccInteg =
+         (dim_==2) ?
+         dynamic_cast<BilinearFormIntegrator*>(new DiffusionIntegrator) :
+         dynamic_cast<BilinearFormIntegrator*>(new CurlCurlIntegrator);
+      ccInteg->SetIntRule(ir);
+      s1_->AddDomainIntegrator(ccInteg);
+      s1_->Assemble();
+      s1_->Finalize();
+      S1_ = new HypreParMatrix;
+   }
+   if ( weakCurl_ == NULL )
+   {
+      weakCurl_ = new ParMixedBilinearForm(HDivFESpace_, HCurlFESpace_);
+      BilinearFormIntegrator * wcurlInteg = new MixedVectorWeakCurlIntegrator;
+      wcurlInteg->SetIntRule(ir);
+      weakCurl_->AddDomainIntegrator(wcurlInteg);
+      weakCurl_->Assemble();
+      weakCurl_->Finalize();
+   }
+   if ( curl_ == NULL )
+   {
+      curl_ = new ParDiscreteCurlOperator(HCurlFESpace_, HDivFESpace_);
+      curl_->Assemble();
+      curl_->Finalize();
+   }
+
+   psi_   = new ParGridFunction(HCurlFESpace_);
+   xCurl_ = new ParGridFunction(HCurlFESpace_);
+}
+
+DivergenceFreeRTProjector::~DivergenceFreeRTProjector()
+{
+   delete psi_;
+   delete xCurl_;
+
+   delete pc_;
+   delete pcg_;
+
+   delete S1_;
+
+   delete s1_;
+   delete weakCurl_;
+}
+
+void
+DivergenceFreeRTProjector::InitSolver() const
+{
+   delete pcg_;
+   delete pc_;
+
+   if (dim_ == 2)
+   {
+      HypreBoomerAMG * amg = new HypreBoomerAMG(*S1_);
+      amg->SetPrintLevel(0);
+      pc_ = amg;
+   }
+   else
+   {
+      HypreAMS * ams = new HypreAMS(*S1_, HCurlFESpace_);
+      ams->SetPrintLevel(0);
+      pc_ = ams;
+   }
+   pcg_ = new HyprePCG(*S1_);
+   pcg_->SetTol(1e-14);
+   pcg_->SetMaxIter(200);
+   pcg_->SetPrintLevel(0);
+   pcg_->SetPreconditioner(*pc_);
+}
+
+void
+DivergenceFreeRTProjector::Mult(const Vector &x, Vector &y) const
+{
+   // Compute the curl of x
+   weakCurl_->Mult(x,*xCurl_);
+
+   // Apply essential BC and form linear system
+   *psi_ = 0.0;
+   s1_->FormLinearSystem(ess_bdr_tdofs_, *psi_, *xCurl_, *S1_, Psi_, RHS_);
+
+   // Solve the linear system for Psi
+   if ( pcg_ == NULL ) { this->InitSolver(); }
+   pcg_->Mult(RHS_, Psi_);
+
+   // Compute the parallel grid function correspoinding to Psi
+   s1_->RecoverFEMSolution(Psi_, *xCurl_, *psi_);
+
+   // Compute the divergence free portion of x
+   curl_->Mult(*psi_, y);
+}
+
+void
+DivergenceFreeRTProjector::Update()
+{
+   delete pcg_; pcg_ = NULL;
+   delete pc_;  pc_  = NULL;
+   delete S1_;  S1_  = new HypreParMatrix;
+
+   psi_->Update();
+   xCurl_->Update();
+
+   if ( ownsS1_ )
+   {
+      s1_->Update();
+      s1_->Assemble();
+      s1_->Finalize();
+   }
+   if ( ownsWeakCurl_ )
+   {
+      weakCurl_->Update();
+      weakCurl_->Assemble();
+      weakCurl_->Finalize();
+   }
+   if ( ownsCurl_ )
+   {
+      curl_->Update();
+      curl_->Assemble();
+      curl_->Finalize();
+   }
+
+   HCurlFESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
+}
+
+IrrotationalRTProjector
+::IrrotationalRTProjector(ParFiniteElementSpace   & HCurlFESpace,
+                          ParFiniteElementSpace   & HDivFESpace,
+                          const int               & irOrder,
+                          ParBilinearForm         * s1,
+                          ParMixedBilinearForm    * weakCurl,
+                          ParDiscreteCurlOperator * curl)
+   : DivergenceFreeRTProjector(HCurlFESpace, HDivFESpace, irOrder,
+                               s1, weakCurl, curl)
+{}
+
+IrrotationalRTProjector::~IrrotationalRTProjector()
+{}
+
+void
+IrrotationalRTProjector::Mult(const Vector &x, Vector &y) const
+{
+   this->DivergenceFreeRTProjector::Mult(x, y);
+   y  -= x;
+   y *= -1.0;
+}
+
+void
+IrrotationalRTProjector::Update()
+{
+   this->DivergenceFreeRTProjector::Update();
 }
 
 void VisualizeMesh(socketstream &sock, const char *vishost, int visport,

miniapps/common/pfem_extras.hpp

@@ -115,16 +115,16 @@ public:
 /// This class computes the irrotational portion of a vector field.
 /// This vector field must be discretized using Nedelec basis
 /// functions.
-class IrrotationalProjector : public Operator
+class IrrotationalNDProjector : public Operator
 {
 public:
-   IrrotationalProjector(ParFiniteElementSpace   & H1FESpace,
-                         ParFiniteElementSpace   & HCurlFESpace,
-                         const int               & irOrder,
-                         ParBilinearForm         * s0 = NULL,
-                         ParMixedBilinearForm    * weakDiv = NULL,
-                         ParDiscreteGradOperator * grad = NULL);
-   virtual ~IrrotationalProjector();
+   IrrotationalNDProjector(ParFiniteElementSpace   & H1FESpace,
+                           ParFiniteElementSpace   & HCurlFESpace,
+                           const int               & irOrder,
+                           ParBilinearForm         * s0 = NULL,
+                           ParMixedBilinearForm    * weakDiv = NULL,
+                           ParDiscreteGradOperator * grad = NULL);
+   virtual ~IrrotationalNDProjector();
 
    // Given a GridFunction 'x' of Nedelec DoFs for an arbitrary vector field,
    // compute the Nedelec DoFs of the irrotational portion, 'y', of
@@ -164,16 +164,16 @@ private:
 /// This class computes the divergence free portion of a vector field.
 /// This vector field must be discretized using Nedelec basis
 /// functions.
-class DivergenceFreeProjector : public IrrotationalProjector
+class DivergenceFreeNDProjector : public IrrotationalNDProjector
 {
 public:
-   DivergenceFreeProjector(ParFiniteElementSpace   & H1FESpace,
-                           ParFiniteElementSpace   & HCurlFESpace,
-                           const int               & irOrder,
-                           ParBilinearForm         * s0 = NULL,
-                           ParMixedBilinearForm    * weakDiv = NULL,
-                           ParDiscreteGradOperator * grad = NULL);
-   virtual ~DivergenceFreeProjector();
+   DivergenceFreeNDProjector(ParFiniteElementSpace   & H1FESpace,
+                             ParFiniteElementSpace   & HCurlFESpace,
+                             const int               & irOrder,
+                             ParBilinearForm         * s0 = NULL,
+                             ParMixedBilinearForm    * weakDiv = NULL,
+                             ParDiscreteGradOperator * grad = NULL);
+   virtual ~DivergenceFreeNDProjector();
 
    // Given a vector 'x' of Nedelec DoFs for an arbitrary vector field,
    // compute the Nedelec DoFs of the divergence free portion, 'y', of
@@ -184,6 +184,79 @@ public:
    void Update();
 };
 
+/// This class computes the divergence free portion of a vector field.
+/// This vector field must be discretized using Raviart-Thomas basis
+/// functions.
+class DivergenceFreeRTProjector : public Operator
+{
+public:
+   DivergenceFreeRTProjector(ParFiniteElementSpace   & HCurlFESpace,
+                             ParFiniteElementSpace   & HDivFESpace,
+                             const int               & irOrder,
+                             ParBilinearForm         * s1 = NULL,
+                             ParMixedBilinearForm    * weakCurl = NULL,
+                             ParDiscreteCurlOperator * curl = NULL);
+   virtual ~DivergenceFreeRTProjector();
+
+   // Given a GridFunction 'x' of Raviart-Thomas DoFs for an arbitrary vector
+   // field, compute the Raviart-Thomas DoFs of the divergence free portion,
+   // 'y', of this vector field.  The resulting GridFunction will satisfy
+   // Div y = 0 to machine precision.
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   void Update();
+
+private:
+   void InitSolver() const;
+
+   ParFiniteElementSpace * HCurlFESpace_;
+   ParFiniteElementSpace * HDivFESpace_;
+
+   ParBilinearForm         * s1_;
+   ParMixedBilinearForm    * weakCurl_;
+   ParDiscreteCurlOperator * curl_;
+
+   ParGridFunction * psi_;
+   ParGridFunction * xCurl_;
+
+   HypreParMatrix * S1_;
+   mutable Vector Psi_;
+   mutable Vector RHS_;
+
+   mutable HypreSolver * pc_;
+   mutable HyprePCG    * pcg_;
+
+   Array<int> ess_bdr_, ess_bdr_tdofs_;
+
+   int dim_;
+   bool ownsS1_;
+   bool ownsWeakCurl_;
+   bool ownsCurl_;
+};
+
+/// This class computes the irrotational portion of a vector field.
+/// This vector field must be discretized using Nedelec basis
+/// functions.
+class IrrotationalRTProjector : public DivergenceFreeRTProjector
+{
+public:
+   IrrotationalRTProjector(ParFiniteElementSpace   & HCurlFESpace,
+                           ParFiniteElementSpace   & HDivFESpace,
+                           const int               & irOrder,
+                           ParBilinearForm         * s1 = NULL,
+                           ParMixedBilinearForm    * weakCurl = NULL,
+                           ParDiscreteCurlOperator * curl = NULL);
+   virtual ~IrrotationalRTProjector();
+
+   // Given a GridFunction 'x' of Raviart-Thomas DoFs for an arbitrary vector
+   // field, compute the Raviart-Thomas DoFs of the irrotational portion,
+   // 'y', of this vector field.  The resulting GridFunction will satisfy
+   // Curl y = 0 to machine precision.
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   void Update();
+};
+
 
 /// Visualize the given parallel mesh object, using a GLVis server on the
 /// specified host and port. Set the visualization window title, and optionally,

miniapps/electromagnetics/tesla_solver.cpp

@@ -152,8 +152,8 @@ TeslaSolver::TeslaSolver(ParMesh & pmesh, int order,
    {
       jr_          = new ParGridFunction(HCurlFESpace_);
       j_           = new ParGridFunction(HCurlFESpace_);
-      DivFreeProj_ = new DivergenceFreeProjector(*H1FESpace_, *HCurlFESpace_,
-                                                 irOrder, NULL, NULL, grad_);
+      DivFreeProj_ = new DivergenceFreeNDProjector(*H1FESpace_, *HCurlFESpace_,
+                                                   irOrder, NULL, NULL, grad_);
    }
 
    if ( kbcs.Size() > 0 )

miniapps/electromagnetics/tesla_solver.hpp

@@ -28,7 +28,7 @@ using miniapps::ND_ParFESpace;
 using miniapps::RT_ParFESpace;
 using miniapps::ParDiscreteGradOperator;
 using miniapps::ParDiscreteCurlOperator;
-using miniapps::DivergenceFreeProjector;
+using miniapps::DivergenceFreeNDProjector;
 
 namespace electromagnetics
 {
@@ -99,8 +99,8 @@ private:
    ParGridFunction * bd_; // Dual of B (HCurl)
    ParGridFunction * jd_; // Dual of J, the rhs vector (HCurl)
 
-   DivergenceFreeProjector * DivFreeProj_;
-   SurfaceCurrent          * SurfCur_;
+   DivergenceFreeNDProjector * DivFreeProj_;
+   SurfaceCurrent            * SurfCur_;
 
    Coefficient       * muInvCoef_; // Dia/Paramagnetic Material Coefficient
    VectorCoefficient * aBCCoef_;   // Vector Potential BC Function

miniapps/tools/hodge-decomp.cpp

@@ -0,0 +1,395 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+//
+//    -------------------------------------------------------------------
+//    Hodge Decomposition Miniapp: Split vector fields into
+//    -------------------------------------------------------------------
+
+#include "../common/pfem_extras.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace mfem::miniapps;
+
+using miniapps::H1_ParFESpace;
+using miniapps::ND_ParFESpace;
+using miniapps::RT_ParFESpace;
+//using miniapps::DivergenceFreeNDProjector;
+using miniapps::DivergenceFreeRTProjector;
+using miniapps::IrrotationalNDProjector;
+//using miniapps::IrrotationalFreeRTProjector;
+
+static int nr_   = 1;
+static int nphi_ = 0;
+
+static double r_ = 0.4;
+static double R_ = 1.1;
+
+void w_exact(const Vector &, Vector &);
+
+double a_exact(const Vector &);
+void da_exact(const Vector &, Vector &);
+
+void b_exact(const Vector &, Vector &);
+void db_exact(const Vector &, Vector &);
+
+void c_exact(const Vector &, Vector &);
+
+int main(int argc, char ** argv)
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/toroid-hex.mesh";
+   int order = 1;
+   int serial_ref_levels = 2;
+   int parallel_ref_levels = 0;
+   bool visualization = 1;
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+
+   int Wx = 0, Wy = 0; // window position
+   int Ww = 350, Wh = 350; // window size
+   int offx = Ww+3, offy = Wh+25; // window offsets
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&serial_ref_levels, "-rs", "--serial-ref-levels",
+                  "Number of serial refinement levels.");
+   args.AddOption(&parallel_ref_levels, "-rp", "--parallel-ref-levels",
+                  "Number of parallel refinement levels.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   for (int l = 0; l < serial_ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. 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(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   int par_ref_levels = parallel_ref_levels;
+   for (int l = 0; l < par_ref_levels; l++)
+   {
+      pmesh.UniformRefinement();
+   }
+
+   H1_ParFESpace fespace_h1(&pmesh, order, pmesh.Dimension());
+   ND_ParFESpace fespace_nd(&pmesh, order, pmesh.Dimension());
+   RT_ParFESpace fespace_rt(&pmesh, order, pmesh.Dimension());
+
+   ParDiscreteGradOperator Grad(&fespace_h1, &fespace_nd);
+   Grad.Assemble();
+   Grad.Finalize();
+
+   ParDiscreteCurlOperator Curl(&fespace_nd, &fespace_rt);
+   Curl.Assemble();
+   Curl.Finalize();
+
+   ParGridFunction a_h1(&fespace_h1);
+   ParGridFunction da_nd(&fespace_nd);
+
+   ParGridFunction b_nd(&fespace_nd);
+   ParGridFunction db_rt(&fespace_rt);
+
+   ParGridFunction w_nd(&fespace_nd);
+   ParGridFunction irr_w_nd(&fespace_nd);
+
+   ParGridFunction w_rt(&fespace_rt);
+   ParGridFunction df_w_rt(&fespace_rt);
+   ParGridFunction w_c_rt(&fespace_rt);
+
+   FunctionCoefficient aCoef(a_exact);
+   VectorFunctionCoefficient daCoef(pmesh.SpaceDimension(), da_exact);
+
+   VectorFunctionCoefficient bCoef(pmesh.SpaceDimension(), b_exact);
+   VectorFunctionCoefficient dbCoef(pmesh.SpaceDimension(), db_exact);
+
+   VectorFunctionCoefficient cCoef(pmesh.SpaceDimension(), c_exact);
+
+   VectorFunctionCoefficient wCoef(pmesh.SpaceDimension(), w_exact);
+
+   VectorGridFunctionCoefficient irr_w_Coef(&irr_w_nd);
+
+   a_h1.ProjectCoefficient(aCoef);
+   b_nd.ProjectCoefficient(bCoef);
+
+   Grad.Mult(a_h1, da_nd);
+   Curl.Mult(b_nd, db_rt);
+
+   double err_a_h1 = a_h1.ComputeL2Error(aCoef);
+   double err_da_nd = da_nd.ComputeL2Error(daCoef);
+
+   double err_b_nd = b_nd.ComputeL2Error(bCoef);
+   double err_db_rt = db_rt.ComputeL2Error(dbCoef);
+
+   if (myid == 0)
+   {
+      cout << "Error in a (H1):  " << err_a_h1 << endl;
+      cout << "Error in da (ND): " << err_da_nd << endl;
+      cout << "Error in b (ND):  " << err_b_nd << endl;
+      cout << "Error in db (RT): " << err_db_rt << endl;
+   }
+
+   w_nd.ProjectCoefficient(wCoef);
+   w_rt.ProjectCoefficient(wCoef);
+
+   double err_w_nd = w_nd.ComputeL2Error(wCoef);
+   double err_w_rt = w_rt.ComputeL2Error(wCoef);
+
+   if (myid == 0)
+   {
+      cout << "Error in w (ND): " << err_w_nd << endl;
+      cout << "Error in w (RT): " << err_w_rt << endl;
+   }
+
+   map<string, socketstream*> socks;
+   {
+      socks["w_nd"] = new socketstream;
+      socks["w_nd"]->precision(8);
+      VisualizeField(*socks["w_nd"], vishost, visport,
+                     w_nd, "w ND", Wx, Wy, Ww, Wh);
+      Wy += offy;
+
+      socks["w_rt"] = new socketstream;
+      socks["w_rt"]->precision(8);
+      VisualizeField(*socks["w_rt"], vishost, visport,
+                     w_rt, "w RT", Wx, Wy, Ww, Wh);
+   }
+
+   IrrotationalNDProjector irr_nd(fespace_h1, fespace_nd, 2 * order + 1);
+   irr_nd.Mult(w_nd, irr_w_nd);
+
+   double err_irr_w_nd = irr_w_nd.ComputeL2Error(daCoef);
+
+   if (myid == 0)
+   {
+      cout << "Error in da (ND): " << err_da_nd << endl;
+      cout << "Error in irr w (ND): " << err_irr_w_nd << endl;
+   }
+   {
+      Wy -= offy;
+      Wx += offx;
+
+      socks["da_nd"] = new socketstream;
+      socks["da_nd"]->precision(8);
+      VisualizeField(*socks["da_nd"], vishost, visport,
+                     irr_w_nd, "irr w ND", Wx, Wy, Ww, Wh);
+   }
+
+   DivergenceFreeRTProjector df_rt(fespace_nd, fespace_rt, 2 * order + 1);
+
+   df_rt.Mult(w_rt, df_w_rt);
+
+   double err_df_w_rt = df_w_rt.ComputeL2Error(dbCoef);
+
+   if (myid == 0)
+   {
+      cout << "Error in df w (RT): " << err_df_w_rt << endl;
+   }
+
+   {
+      Wy += offy;
+
+      socks["db_rt"] = new socketstream;
+      socks["db_rt"]->precision(8);
+      VisualizeField(*socks["db_rt"], vishost, visport,
+                     df_w_rt, "df w RT", Wx, Wy, Ww, Wh);
+   }
+
+   w_c_rt.ProjectCoefficient(irr_w_Coef);
+   w_c_rt += df_w_rt;
+   w_c_rt *= -1.0;
+   w_c_rt += w_rt;
+
+   double err_w_c_rt = w_c_rt.ComputeL2Error(cCoef);
+   if (myid == 0)
+   {
+      cout << "Error in c (RT): " << err_w_c_rt << endl;
+   }
+   {
+      Wx += offx;
+
+      socks["w_c_rt"] = new socketstream;
+      socks["w_c_rt"]->precision(8);
+      VisualizeField(*socks["w_c_rt"], vishost, visport,
+                     w_c_rt, "w c RT", Wx, Wy, Ww, Wh);
+   }
+}
+
+void w_exact(const Vector &x, Vector &w)
+{
+   w.SetSize(3);
+
+   double da_data[3];
+   double db_data[3];
+   Vector da(da_data, 3);
+   Vector db(db_data, 3);
+
+   da_exact(x, da);
+   db_exact(x, db);
+   c_exact(x, w);
+   w += da;
+   w += db;
+}
+
+double a_exact(const Vector &x)
+{
+   double r = sqrt(x[0] * x[0] + x[1] * x[1]);
+   double phi = atan2(x[1], x[0]);
+
+   double ar = 0.5 * M_PI * nr_ * (r - R_) / r_;
+   double ap = phi * nphi_;
+   double az = 0.5 * M_PI * nr_ * x[2] / r_;
+
+   return (2.0 * r_ / (M_PI * nr_)) * cos(ar) * cos(ap) * cos(az);
+}
+
+void da_exact(const Vector &x, Vector &da)
+{
+   da.SetSize(3);
+
+   double r = sqrt(x[0] * x[0] + x[1] * x[1]);
+   double phi = atan2(x[1], x[0]);
+
+   double ar = 0.5 * M_PI * nr_ * (r - R_) / r_;
+   double ap = phi * nphi_;
+   double az = 0.5 * M_PI * nr_ * x[2] / r_;
+
+   double drdx = x[0] / r;
+   double drdy = x[1] / r;
+   double dpdx = -x[1] / (r * r);
+   double dpdy =  x[0] / (r * r);
+
+   double dardr = 0.5 * M_PI * nr_ / r_;
+   double dapdp = (double)nphi_;
+   double dazdz = 0.5 * M_PI * nr_ / r_;
+
+   da(0) = -(dardr * drdx * sin(ar) * cos(ap) +
+             dapdp * dpdx * cos(ar) * sin(ap)
+            ) * cos(az);
+   da(1) = -(dardr * drdy * sin(ar) * cos(ap) +
+             dapdp * dpdy * cos(ar) * sin(ap)
+            ) * cos(az);
+   da(2) = -dazdz * cos(ar) * cos(ap) * sin(az);
+   da *= (2.0 * r_ / (M_PI * nr_));
+}
+
+void b_exact(const Vector &x, Vector &b)
+{
+   b.SetSize(3);
+
+   double r = sqrt(x[0] * x[0] + x[1] * x[1]);
+   double phi = atan2(x[1], x[0]);
+
+   double ar = 0.5 * M_PI * nr_ * (r - R_) / r_;
+   double ap = phi * nphi_;
+   double az = 0.5 * M_PI * nr_ * x[2] / r_;
+
+   double cp = x[0] / r;
+   double sp = x[1] / r;
+
+   b(0) = cp * cos(ap) * cos(az);
+   b(1) = sp * cos(ap) * cos(az);
+   b(2) = cos(ar) * cos(ap);
+
+   b *= r_ / (M_PI * nr_);
+}
+
+void db_exact(const Vector &x, Vector &db)
+{
+   db.SetSize(3);
+
+   double r = sqrt(x[0] * x[0] + x[1] * x[1]);
+   double phi = atan2(x[1], x[0]);
+
+   double ar = 0.5 * M_PI * nr_ * (r - R_) / r_;
+   double ap = phi * nphi_;
+   double az = 0.5 * M_PI * nr_ * x[2] / r_;
+
+   double cp = x[0] / r;
+   double sp = x[1] / r;
+
+   double drdx = x[0] / r;
+   double drdy = x[1] / r;
+   double dpdx = -x[1] / (r * r);
+   double dpdy =  x[0] / (r * r);
+
+   double dardr = 0.5 * M_PI * nr_ / r_;
+   double dapdp = (double)nphi_;
+   double dazdz = 0.5 * M_PI * nr_ / r_;
+
+   double dcpdy = -drdy * cp / r;
+   double dspdx = -drdx * sp / r;
+
+   db(0) = -(dardr * drdy * sin(ar) * cos(ap) +
+             dapdp * dpdy * cos(ar) * sin(ap) -
+             dazdz * sp * cos(ap) * sin(az));
+   db(1) = (dardr * drdx * sin(ar) * cos(ap) +
+            dapdp * dpdx * cos(ar) * sin(ap) -
+            dazdz * cp * cos(ap) * sin(az));
+   db(2) = (dspdx * cos(ap) * cos(az) - dapdp * dpdx * sp * sin(ap)
+            -dcpdy * cos(ap) * cos(az) + dapdp * dpdy * cp * sin(ap)
+           ) * cos(az);
+
+   db *= r_ / (M_PI * nr_);
+}
+
+void c_exact(const Vector &x, Vector &c)
+{
+   c.SetSize(3);
+   c = 0.0;
+   c(0) = -x[1];
+   c(1) = x[0];
+
+   double r2 = x[0] * x[0] + x[1] * x[1];
+   c *= 0.7 / r2;
+   // c(0) = sin(kappa * x[0]) * (cos(kappa * x[1]) - cos(kappa * x[2]));
+   // c(1) = sin(kappa * x[1]) * (cos(kappa * x[2]) - cos(kappa * x[0]));
+   // c(2) = sin(kappa * x[2]) * (cos(kappa * x[0]) - cos(kappa * x[1]));
+}

miniapps/tools/makefile

@@ -23,7 +23,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 
 SEQ_MINIAPPS = display-basis load-dc convert-dc lor-transfer
-PAR_MINIAPPS =
+PAR_MINIAPPS = hodge-decomp
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -35,7 +35,7 @@ endif
 .PHONY: all clean clean-build clean-exec
 .PRECIOUS: %.o
 
-COMMON_O=../common/fem_extras.o ../common/mesh_extras.o
+COMMON_O=../common/fem_extras.o ../common/pfem_extras.o ../common/mesh_extras.o
 
 all: $(MINIAPPS)
 
@@ -52,6 +52,9 @@ display-basis: %: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
 	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<)
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_O) $(MFEM_LIBS)
 
+hodge-decomp: %: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(COMMON_O) $(MFEM_LIBS)
+
 # Rules for compiling dependencies
 $(COMMON_O): %.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
 	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@