Merge remote-tracking branch 'origin/master' into cyl-coords-dev

# Conflicts: # .gitignore # examples/CMakeLists.txt # examples/makefile
2025-10-20 14:08:39 -07:00 · 2021-04-27 18:29:10 -07:00 · 2021-04-27 18:27:25 -07:00 · 2021-04-27 18:26:00 -07:00 · 2021-04-27 10:33:43 -07:00 · 2021-04-18 19:18:51 -07:00
7 changed files with 2281 additions and 1 deletions
@@ -55,8 +55,11 @@ doc/warnings.log

 examples/ex[0-9]
 examples/ex[0-9]p
+examples/ex[0-9]-orth
+examples/ex[0-9]p-orth
 examples/ex1[04-9]
 examples/ex1[0-9]p
+examples/ex1[0-9]p-cyl
 examples/ex2[0-9]
 examples/ex2[0-9]p
 examples/ex3[0-9]
@@ -52,6 +52,7 @@ if (MFEM_USE_MPI)
  list(APPEND ALL_EXE_SRCS
    ex0p.cpp
    ex1p.cpp
+    ex1p-orth.cpp
    ex2p.cpp
    ex3p.cpp
    ex4p.cpp
@@ -62,8 +63,11 @@ if (MFEM_USE_MPI)
    ex9p.cpp
    ex10p.cpp
    ex11p.cpp
+    ex11p-cyl.cpp
    ex12p.cpp
    ex13p.cpp
+    ex13p-cyl.cpp
+    ex13p-cyl-3d.cpp
    ex14p.cpp
    ex15p.cpp
    ex16p.cpp
@@ -0,0 +1,399 @@
+//                       MFEM Example 11-cyl - Parallel Version
+//
+// Compile with: make ex11p-cyl
+//
+// Sample runs:  mpirun -np 4 ex11p-cyl
+//               mpirun -np 4 ex11p-cyl -o 2
+//               mpirun -np 4 ex11p-cyl -o 2 -e 0
+//
+// Description:  This example code demonstrates the use of MFEM to solve PDEs
+//               on an axisymmetric domain.  The eigenvalue problem:
+//                  -Delta u = lambda u
+//               with homogeneous Dirichlet boundary conditions is solved on
+//               a cylindrical domain by meshing only a rectangle in the
+//               rho, z plane.  In cylindrical coordinates the weak form of
+//               the eigenvalue problem is given by:
+//                  (rho Grad(u), Grad(v)) = lambda (rho u, v)
+//
+//               We compute the five lowest eigenmodes by discretizing
+//               the Laplacian and Mass operators using a FE space of the
+//               specified order and compare to the known values.  Because the
+//               eigenvalue spectrum of a domain is unique this provides a
+//               reliable test that the axisymmetric domain is being faithfully
+//               characterized.
+//
+//               The example highlights the use of specialized coefficients
+//               with existing operators to mimic axisymmetric domains.  The
+//               gradient of each eigenmode is also computed and displayed to
+//               ilustrate that no special steps need to be taken to compute
+//               gradients in this coordinate system.
+//
+//               We recommend viewing Example 11 before viewing this example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Zeros of Bessel function J_0
+static double J0z[] = {2.40482555769577,
+                       5.52007811028631,
+                       8.65372791291101,
+                       11.7915344390143
+                      };
+
+// Modes numbers in the rho and z directions for the first five eigenmodes
+static int mode_nums[] = {0, 1,
+                          1, 1,
+                          0, 2,
+                          1, 2,
+                          2, 1
+                         };
+
+double rhoFunc(const Vector &x)
+{
+   return x[0];
+}
+
+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.
+   int nr = 1;
+   int nz = 1;
+   int el_type_flag = 1;
+   Element::Type el_type;
+   int ser_ref_levels = 2;
+   int par_ref_levels = 1;
+   int order = 1;
+   int nev = 5;
+   int seed = 75;
+   bool slu_solver  = false;
+   bool sp_solver = false;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&nz, "-nz", "--num-elements-z",
+                  "Number of elements in z-direction.");
+   args.AddOption(&nr, "-nr", "--num-elements-rho",
+                  "Number of elements in radial direction.");
+   args.AddOption(&el_type_flag, "-e", "--element-type",
+                  "Element type: 0 - Triangle, 1 - Quadrilateral.");
+   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",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&seed, "-s", "--seed",
+                  "Random seed used to initialize LOBPCG.");
+#ifdef MFEM_USE_SUPERLU
+   args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
+                  "--no-superlu", "Use the SuperLU Solver.");
+#endif
+#ifdef MFEM_USE_STRUMPACK
+   args.AddOption(&sp_solver, "-sp", "--strumpack", "-no-sp",
+                  "--no-strumpack", "Use the STRUMPACK Solver.");
+#endif
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (slu_solver && sp_solver)
+   {
+      if (myid == 0)
+         cout << "WARNING: Both SuperLU and STRUMPACK have been selected,"
+              << " please choose either one." << endl
+              << "         Defaulting to SuperLU." << endl;
+      sp_solver = false;
+   }
+   // The command line options are also passed to the STRUMPACK
+   // solver. So do not exit if some options are not recognized.
+   if (!sp_solver)
+   {
+      if (!args.Good())
+      {
+         if (myid == 0)
+         {
+            args.PrintUsage(cout);
+         }
+         MPI_Finalize();
+         return 1;
+      }
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // The output mesh could be quadrilaterals or triangles
+   el_type = (el_type_flag == 0) ? Element::TRIANGLE : Element::QUADRILATERAL;
+   if (el_type != Element::TRIANGLE && el_type != Element::QUADRILATERAL)
+   {
+      cout << "Unsupported element type" << endl;
+      exit(1);
+   }
+
+   // 3. Prepare a rectangular mesh with the desired dimensions and element
+   //    type.  Other 2D meshes could be used but then we couldn't check the
+   //    eigenvalues.
+   Mesh *mesh = new Mesh(nr, nz, el_type);
+   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 (2 by default, or
+   //    specified on the command line with -rs).
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution (1 time by
+   //    default, or specified on the command line with -rp). 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();
+   }
+
+   // 6. Define a parallel finite element space on the parallel mesh. Here we
+   //    use continuous Lagrange finite elements (H1) of the specified order.
+   //    We also create a Nedelec space to represent the gradients of the modes.
+   H1_FECollection fec_h1(order, dim);
+   ND_FECollection fec_nd(order, dim);
+   ParFiniteElementSpace fespace_h1(pmesh, &fec_h1);
+   ParFiniteElementSpace fespace_nd(pmesh, &fec_nd);
+   HYPRE_Int size = fespace_h1.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of unknowns: " << size << endl;
+   }
+
+   // 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
+   //    element space. The first corresponds to the Laplacian operator -Delta,
+   //    while the second is a simple mass matrix needed on the right hand side
+   //    of the generalized eigenvalue problem below. The boundary conditions
+   //    are implemented by elimination with special values on the diagonal to
+   //    shift the Dirichlet eigenvalues out of the computational range. After
+   //    serial and parallel assembly we extract the corresponding parallel
+   //    matrices A and M.
+   FunctionCoefficient rhoCoef(rhoFunc);
+   Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+   ess_bdr = 1;     // Homogeneous Dirichlet BCs everywhere except for
+   ess_bdr[3] = 0;  // attribute 4 which is on the axis of symmetry.
+
+   ParBilinearForm *a = new ParBilinearForm(&fespace_h1);
+   a->AddDomainIntegrator(new DiffusionIntegrator(rhoCoef));
+   a->Assemble();
+   a->EliminateEssentialBCDiag(ess_bdr, 1.0);
+   a->Finalize();
+
+   ParBilinearForm *m = new ParBilinearForm(&fespace_h1);
+   m->AddDomainIntegrator(new MassIntegrator(rhoCoef));
+   m->Assemble();
+   // shift the eigenvalue corresponding to eliminated dofs to a large value
+   m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
+   m->Finalize();
+
+   HypreParMatrix *A = a->ParallelAssemble();
+   HypreParMatrix *M = m->ParallelAssemble();
+
+#if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK)
+   Operator * Arow = NULL;
+#ifdef MFEM_USE_SUPERLU
+   if (slu_solver)
+   {
+      Arow = new SuperLURowLocMatrix(*A);
+   }
+#endif
+#ifdef MFEM_USE_STRUMPACK
+   if (sp_solver)
+   {
+      Arow = new STRUMPACKRowLocMatrix(*A);
+   }
+#endif
+#endif
+
+   delete a;
+   delete m;
+
+   // 8. Define and configure the LOBPCG eigensolver and the BoomerAMG
+   //    preconditioner for A to be used within the solver. Set the matrices
+   //    which define the generalized eigenproblem A x = lambda M x.
+   Solver * precond = NULL;
+   if (!slu_solver && !sp_solver)
+   {
+      HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
+      amg->SetPrintLevel(0);
+      precond = amg;
+   }
+   else
+   {
+#ifdef MFEM_USE_SUPERLU
+      if (slu_solver)
+      {
+         SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
+         superlu->SetPrintStatistics(false);
+         superlu->SetSymmetricPattern(true);
+         superlu->SetColumnPermutation(superlu::PARMETIS);
+         superlu->SetOperator(*Arow);
+         precond = superlu;
+      }
+#endif
+#ifdef MFEM_USE_STRUMPACK
+      if (sp_solver)
+      {
+         STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
+         strumpack->SetPrintFactorStatistics(true);
+         strumpack->SetPrintSolveStatistics(false);
+         strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
+         strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
+         strumpack->DisableMatching();
+         strumpack->SetOperator(*Arow);
+         strumpack->SetFromCommandLine();
+         precond = strumpack;
+      }
+#endif
+   }
+
+
+   HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
+   lobpcg->SetNumModes(nev);
+   lobpcg->SetRandomSeed(seed);
+   lobpcg->SetPreconditioner(*precond);
+   lobpcg->SetMaxIter(200);
+   lobpcg->SetTol(1e-8);
+   lobpcg->SetPrecondUsageMode(1);
+   lobpcg->SetPrintLevel(1);
+   lobpcg->SetMassMatrix(*M);
+   lobpcg->SetOperator(*A);
+
+   // 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
+   //    parallel grid function to represent each of the eigenmodes returned by
+   //    the solver.  Also define a discrete gradient operator.
+   Array<double> eigenvalues;
+   lobpcg->Solve();
+   lobpcg->GetEigenvalues(eigenvalues);
+   ParGridFunction x(&fespace_h1);
+   ParGridFunction dx(&fespace_nd);
+
+   ParDiscreteLinearOperator grad(&fespace_h1, &fespace_nd);
+   grad.AddDomainInterpolator(new GradientInterpolator());
+   grad.Assemble();
+
+   if ( myid == 0 )
+   {
+      // Display the eigenvalues and their relative errors
+      cout << "\nRelative error in eigenvalues:\n";
+      for (int i=0; i<nev; i++)
+      {
+         double lambda =
+            pow(J0z[mode_nums[2*i]], 2) +
+            pow(M_PI * mode_nums[2*i+1], 2);
+         cout << "Lambda " << i+1 << '/' << nev << " = " << eigenvalues[i]
+              << ", rel err = " << fabs(eigenvalues[i] - lambda) / lambda
+              << endl;
+      }
+      cout << endl;
+   }
+
+   // 10. Save the refined mesh and the modes in parallel. This output can be
+   //     viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
+   {
+      ostringstream mesh_name, mode_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      for (int i=0; i<nev; i++)
+      {
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);
+
+         mode_name << "mode_" << setfill('0') << setw(2) << i << "."
+                   << setfill('0') << setw(6) << myid;
+
+         ofstream mode_ofs(mode_name.str().c_str());
+         mode_ofs.precision(8);
+         x.Save(mode_ofs);
+         mode_name.str("");
+      }
+   }
+
+   // 11. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mode_sock(vishost, visport);
+      mode_sock.precision(8);
+      socketstream grad_sock(vishost, visport);
+      grad_sock.precision(8);
+
+      for (int i=0; i<nev; i++)
+      {
+         if ( myid == 0 )
+         {
+            cout << "Eigenmode " << i+1 << '/' << nev
+                 << ", Lambda = " << eigenvalues[i] << endl;
+         }
+
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);
+
+         grad.Mult(x, dx);
+
+         mode_sock << "parallel " << num_procs << " " << myid << "\n"
+                   << "solution\n" << *pmesh << x << flush
+                   << "window_title 'Eigenmode " << i+1 << '/' << nev
+                   << ", Lambda = " << eigenvalues[i] << "'" << endl;
+
+         grad_sock << "parallel " << num_procs << " " << myid << "\n"
+                   << "solution\n" << *pmesh << dx << flush
+                   << "window_title 'Grad of Eigenmode " << i+1 << '/' << nev
+                   << ", Lambda = " << eigenvalues[i] << "'"
+                   << "window_geometry 400 0 400 350" << endl;
+
+         char c;
+         if (myid == 0)
+         {
+            cout << "press (q)uit or (c)ontinue --> " << flush;
+            cin >> c;
+         }
+         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      mode_sock.close();
+   }
+
+   // 12. Free the used memory.
+   delete lobpcg;
+   delete precond;
+   delete M;
+   delete A;
+#if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK)
+   delete Arow;
+#endif
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
@@ -0,0 +1,563 @@
+//                       MFEM Example 13-cyl - Parallel Version
+//
+// Compile with: make ex13p-cyl
+//
+// Sample runs:  mpirun -np 4 ex13p-cyl
+//               mpirun -np 4 ex13p-cyl -o 2
+//               mpirun -np 4 ex13p-cyl -o 2 -e 0
+//
+// Description:  This example code demonstrates the use of MFEM to solve the
+//               Maxwell (electromagnetic) eigenvalue problem on an
+//               axisymmetric domain. The eigenvalue problem:
+//                 curl curl E = lambda E
+//               with homogeneous Dirichlet boundary conditions E x n = 0 is
+//               solved on a cylindrical domain by meshing only a rectangle in
+//               the rho, z plane.  In cylindrical coordinates the weak form of
+//               the eigenvalue problem is given by:
+//                  (rho Curl(u), Curl(v)) = lambda (rho u, v)
+//
+//               We compute the eight lowest nonzero eigenmodes by discretizing
+//               the curl curl operator using a Nedelec FE space of the
+//               specified order and compare to the known values.  Because the
+//               eigenvalue spectrum of a domain is unique this provides a
+//               reliable test that the axisymmetric domain is being faithfully
+//               characterized.
+//
+//               In two dimensions the curl curl operator, with isotropic
+//               material coefficients, splits into two separate PDEs. The rho
+//               and z components form a 2D vector field in the rho-z plane
+//               which is discretized with the 2D Nedelec vector basis
+//               functions. The Maxwell eigenvalue problem for the rho-z field
+//               can be written with cartesian operators as:
+//                   curl (rho curl E_rz) = lambda rho E_rz
+//
+//               The angular component can be discretized with the 2D H1 scalar
+//               basis functions. The angular portion of the eigenvalue problem
+//               can be written with cartesian operators as:
+//                   -div (rho grad E_phi) + (1/rho) E_phi = lambda rho E_phi
+//
+//               The example highlights the use of specialized coefficients
+//               with existing operators to mimic axisymmetric domains.  The
+//               curl of each eigenmode is also computed and displayed to
+//               ilustrate that no special steps need to be taken to compute
+//               the curl in this coordinate system.
+//
+//               We recommend viewing examples 13 and 11-cyl before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Zeros of Bessel function J_0
+static double J0z[] = {2.404825557695773,
+                       5.520078110286312,
+                       8.653727912911013,
+                       11.79153443901428
+                      };
+
+// Zeros of Bessel function J_1
+static double J1z[] = {3.831705970207512,
+                       7.015586669815622,
+                       10.17346813506272,
+                       13.32369193631422
+                      };
+
+// Modes numbers in the rho and z directions for the first five eigenmodes
+static int mode_nums_rz[] = {0, 0,
+                             0, 1,
+                             1, 0,
+                             1, 1,
+                             0, 2,
+			     1, 2,
+			     2, 0,
+			     2, 1
+                            };
+
+// Modes numbers in the phi direction for the first five eigenmodes
+static int mode_nums_phi[] = {0, 1,
+                              0, 2,
+                              1, 1,
+                              1, 2,
+                              0, 3,
+			      2, 1,
+			      1, 3,
+			      2, 2
+                             };
+
+// Mode polarizations for the first several modes; 0 - rz, 1 - phi
+static int mode_type[] = {0,0,1,0,0,0,1,1,0,0,0,1,0,1,1,0,1,0,1};
+
+double rhoFunc(const Vector &x)
+{
+   return x[0];
+}
+
+double rhoInvFunc(const Vector &x)
+{
+   return 1.0 / x[0];
+}
+
+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.
+   int dim = 2;
+   int nr = 1;
+   int nz = 1;
+   int el_type_flag = 1;
+   Element::Type el_type;
+   int ser_ref_levels = 2;
+   int par_ref_levels = 1;
+   int order = 1;
+   int nev = 8;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&nz, "-nz", "--num-elements-z",
+                  "Number of elements in z-direction.");
+   args.AddOption(&nr, "-nr", "--num-elements-rho",
+                  "Number of elements in radial direction.");
+   args.AddOption(&el_type_flag, "-e", "--element-type",
+                  "Element type: 0 - Triangle, 1 - Quadrilateral.");
+   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",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   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);
+   }
+
+   // The output mesh could be quadrilaterals or triangles
+   el_type = (el_type_flag == 0) ? Element::TRIANGLE : Element::QUADRILATERAL;
+   if (el_type != Element::TRIANGLE && el_type != Element::QUADRILATERAL)
+   {
+      cout << "Unsupported element type" << endl;
+      exit(1);
+   }
+
+   // 3. Prepare a rectangular mesh with the desired dimensions and element
+   //    type.  Other 2D meshes could be used but then we couldn't check the
+   //    eigenvalues.
+   ParMesh pmesh;
+   {
+      Mesh mesh = Mesh::MakeCartesian2D(nr, nz, el_type);
+
+      // 4. Refine the serial mesh on all processors to increase the resolution. In
+      //    this example we do 'ref_levels' of uniform refinement (2 by default, or
+      //    specified on the command line with -rs).
+      for (int lev = 0; lev < ser_ref_levels; lev++)
+      {
+         mesh.UniformRefinement();
+      }
+
+      // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+      //    this mesh further in parallel to increase the resolution (1 time by
+      //    default, or specified on the command line with -rp). Once the parallel
+      //    mesh is defined, the serial mesh can be deleted.
+      pmesh = ParMesh(MPI_COMM_WORLD, mesh);
+      for (int lev = 0; lev < par_ref_levels; lev++)
+      {
+         pmesh.UniformRefinement();
+      }
+      pmesh.ReorientTetMesh();
+   }
+
+   // 6. Define a parallel finite element space on the parallel mesh. Here we
+   //    use the Nedelec finite elements (ND) of the specified order. We also
+   //    create an L2 space to represent the z-component of the curl of the
+   //    modes.
+   ND_FECollection fec_nd(order, dim);
+   RT_FECollection fec_rt(order - 1, dim);
+   H1_FECollection fec_ndp(order, dim);
+   L2_FECollection fec_rtp(order - 1, dim,
+                           BasisType::GaussLegendre, FiniteElement::INTEGRAL);
+   L2_FECollection fec_l2(order - 1, dim);
+   ParFiniteElementSpace fespace_nd(&pmesh, &fec_nd);
+   ParFiniteElementSpace fespace_rt(&pmesh, &fec_rt);
+   ParFiniteElementSpace fespace_ndp(&pmesh, &fec_ndp);
+   ParFiniteElementSpace fespace_rtp(&pmesh, &fec_rtp);
+   ParFiniteElementSpace fespace_l2(&pmesh, &fec_l2);
+   HYPRE_Int size = fespace_nd.GlobalTrueVSize() +
+                    fespace_ndp.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of unknowns: " << size << endl;
+   }
+
+   // 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
+   //    element space. The first corresponds to the curl curl, while the second
+   //    is a simple mass matrix needed on the right hand side of the
+   //    generalized eigenvalue problem below. The boundary conditions are
+   //    implemented by marking all the boundary attributes from the mesh as
+   //    essential. The corresponding degrees of freedom are eliminated with
+   //    special values on the diagonal to shift the Dirichlet eigenvalues out
+   //    of the computational range. After serial and parallel assembly we
+   //    extract the corresponding parallel matrices A and M.
+   FunctionCoefficient rhoCoef(rhoFunc);
+   FunctionCoefficient rhoInvCoef(rhoInvFunc);
+   Array<int> ess_bdr_rz(pmesh.bdr_attributes.Size());
+   Array<int> ess_bdr_phi(pmesh.bdr_attributes.Size());
+   ess_bdr_rz = 1;  ess_bdr_rz[3] = 0;
+   ess_bdr_phi = 1;
+
+   ParBilinearForm *arz = new ParBilinearForm(&fespace_nd);
+   arz->AddDomainIntegrator(new CurlCurlIntegrator(rhoCoef));
+   arz->Assemble();
+   arz->EliminateEssentialBCDiag(ess_bdr_rz, 1.0);
+   arz->Finalize();
+
+   ParBilinearForm *mrz = new ParBilinearForm(&fespace_nd);
+   mrz->AddDomainIntegrator(new VectorFEMassIntegrator(rhoCoef));
+   mrz->Assemble();
+   // shift the eigenvalue corresponding to eliminated dofs to a large value
+   mrz->EliminateEssentialBCDiag(ess_bdr_rz, numeric_limits<double>::min());
+   mrz->Finalize();
+
+   ParBilinearForm *aphi = new ParBilinearForm(&fespace_ndp);
+   aphi->AddDomainIntegrator(new MassIntegrator(rhoInvCoef));
+   aphi->AddDomainIntegrator(new DiffusionIntegrator(rhoCoef));
+   aphi->Assemble();
+   aphi->EliminateEssentialBCDiag(ess_bdr_phi, 1.0);
+   aphi->Finalize();
+
+   ParBilinearForm *mphi = new ParBilinearForm(&fespace_ndp);
+   mphi->AddDomainIntegrator(new MassIntegrator(rhoCoef));
+   mphi->Assemble();
+   // shift the eigenvalue corresponding to eliminated dofs to a large value
+   mphi->EliminateEssentialBCDiag(ess_bdr_phi, numeric_limits<double>::min());
+   mphi->Finalize();
+
+   HypreParMatrix *Arz = arz->ParallelAssemble();
+   HypreParMatrix *Mrz = mrz->ParallelAssemble();
+
+   HypreParMatrix *Aphi = aphi->ParallelAssemble();
+   HypreParMatrix *Mphi = mphi->ParallelAssemble();
+
+   delete arz;
+   delete mrz;
+   delete aphi;
+   delete mphi;
+
+   // 8. Define and configure the AME eigensolver and the AMS preconditioner for
+   //    A to be used within the solver. Set the matrices which define the
+   //    generalized eigenproblem A x = lambda M x.
+   HypreAMS *ams = new HypreAMS(*Arz,&fespace_nd);
+   ams->SetPrintLevel(0);
+   ams->SetSingularProblem();
+
+   HypreAME *ame = new HypreAME(MPI_COMM_WORLD);
+   ame->SetNumModes(nev);
+   ame->SetPreconditioner(*ams);
+   ame->SetMaxIter(100);
+   ame->SetTol(1e-8);
+   ame->SetPrintLevel(1);
+   ame->SetMassMatrix(*Mrz);
+   ame->SetOperator(*Arz);
+
+   HypreBoomerAMG *amg = new HypreBoomerAMG(*Aphi);
+   amg->SetPrintLevel(0);
+
+   HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
+   lobpcg->SetNumModes(nev);
+   // lobpcg->SetRandomSeed(seed);
+   lobpcg->SetPreconditioner(*amg);
+   lobpcg->SetMaxIter(200);
+   lobpcg->SetTol(1e-8);
+   lobpcg->SetPrecondUsageMode(1);
+   lobpcg->SetPrintLevel(1);
+   lobpcg->SetMassMatrix(*Mphi);
+   lobpcg->SetOperator(*Aphi);
+
+   // 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
+   //    parallel grid function to represent each of the eigenmodes returned by
+   //    the solver. Also, define a discrete curl operator.
+   Array<double> eigenvalues_rz;
+   Array<double> eigenvalues_phi;
+
+   ame->Solve();
+   ame->GetEigenvalues(eigenvalues_rz);
+
+   lobpcg->Solve();
+   lobpcg->GetEigenvalues(eigenvalues_phi);
+
+   Array<double> eigenvalues;
+
+   ParGridFunction x_rz(&fespace_nd);
+   ParGridFunction x_phi(&fespace_ndp);
+   ParGridFunction dx_rz(&fespace_rt);
+   ParGridFunction dx_phi(&fespace_rtp);
+
+   ParDiscreteLinearOperator curl(&fespace_nd, &fespace_rtp);
+   curl.AddDomainInterpolator(new CurlInterpolator());
+   curl.Assemble();
+   curl.Finalize();
+
+   ParDiscreteLinearOperator curl2(&fespace_ndp, &fespace_rt);
+   curl2.AddDomainInterpolator(new CurlInterpolator());
+   curl2.Assemble();
+   curl2.Finalize();
+
+   // This is one workaround for GLVis limitations
+   ParGridFunction dx_l2(&fespace_l2);
+
+   GridFunctionCoefficient dxCoef(&dx_phi);
+
+   if ( myid == 0 )
+   {
+      cout << "\nRelative error in eigenvalues of RZ modes:\n";
+      for (int i=0; i<nev; i++)
+      {
+         double lambda =
+            pow(J0z[mode_nums_rz[2*i]], 2) +
+            pow(M_PI * mode_nums_rz[2*i+1], 2);
+         cout << "Lambda " << i+1 << '/' << nev << " = " << eigenvalues_rz[i]
+              << ", rel err = " << fabs(eigenvalues_rz[i] - lambda) / lambda
+              << endl;
+      }
+      cout << endl;
+      cout << "\nRelative error in eigenvalues of Phi modes:\n";
+      for (int i=0; i<nev; i++)
+      {
+         double lambda =
+            pow(J1z[mode_nums_phi[2*i]], 2) +
+            pow(M_PI * mode_nums_phi[2*i+1], 2);
+         cout << "Lambda " << i+1 << '/' << nev << " = " << eigenvalues_phi[i]
+              << ", rel err = " << fabs(eigenvalues_phi[i] - lambda) / lambda
+              << endl;
+      }
+      cout << endl;
+   }
+
+   // 10. Save the refined mesh and the modes in parallel. This output can be
+   //     viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
+   {
+      ostringstream mesh_name, mode_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh.Print(mesh_ofs);
+
+      for (int i=0; i<nev; i++)
+      {
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x_rz = ame->GetEigenvector(i);
+
+         mode_name << "mode_rz_" << setfill('0') << setw(2) << i << "."
+                   << setfill('0') << setw(6) << myid;
+
+         ofstream mode_rz_ofs(mode_name.str().c_str());
+         mode_rz_ofs.precision(8);
+         x_rz.Save(mode_rz_ofs);
+         mode_name.str("");
+
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x_phi = lobpcg->GetEigenvector(i);
+
+         mode_name << "mode_phi_" << setfill('0') << setw(2) << i << "."
+                   << setfill('0') << setw(6) << myid;
+
+         ofstream mode_phi_ofs(mode_name.str().c_str());
+         mode_phi_ofs.precision(8);
+         x_phi.Save(mode_phi_ofs);
+         mode_name.str("");
+      }
+   }
+
+   // 11. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mode_rz_sock(vishost, visport);
+      socketstream curl_rz_sock(vishost, visport);
+      socketstream mode_phi_sock(vishost, visport);
+      socketstream curl_phi_sock(vishost, visport);
+      mode_rz_sock.precision(8);
+      curl_rz_sock.precision(8);
+      mode_phi_sock.precision(8);
+      curl_phi_sock.precision(8);
+
+      int irz = 0;
+      int iphi = 0;
+
+      for (int i=0; i<nev; i++)
+      {
+         if (mode_type[i] == 0)
+         {
+            if ( myid == 0 )
+            {
+               cout << "Eigenmode " << i+1 << '/' << nev
+                    << ", Lambda = " << eigenvalues_rz[irz] << endl;
+            }
+
+            // convert eigenvector from HypreParVector to ParGridFunction
+            x_rz = ame->GetEigenvector(irz);
+
+            curl.Mult(x_rz, dx_phi);
+            dx_l2.ProjectCoefficient(dxCoef);
+
+            mode_rz_sock << "parallel " << num_procs << " " << myid << "\n"
+			 << "solution\n" << pmesh << x_rz << flush
+			 << "window_title 'Eigenmode " << i+1 << '/' << nev
+			 << ", Lambda = " << eigenvalues_rz[irz] << "'\n";
+	    if (irz == 0)
+	      {
+		mode_rz_sock << "keys vvv\n";
+	      }
+	    mode_rz_sock << flush;
+	    
+            curl_rz_sock << "parallel " << num_procs << " " << myid << "\n"
+			 << "solution\n" << pmesh << dx_l2 << flush
+			 << "window_title 'Curl of Eigenmode " << i+1
+			 << '/' << nev
+			 << ", Lambda = " << eigenvalues_rz[irz] << "' "
+			 << "window_geometry 400 0 400 350\n" << flush;
+
+            irz++;
+         }
+         else
+         {
+            if ( myid == 0 )
+            {
+               cout << "Eigenmode " << i+1 << '/' << nev
+                    << ", Lambda = " << eigenvalues_phi[iphi] << endl;
+            }
+
+            // convert eigenvector from HypreParVector to ParGridFunction
+            x_phi = lobpcg->GetEigenvector(iphi);
+
+            curl2.Mult(x_phi, dx_rz);
+
+            mode_phi_sock << "parallel " << num_procs << " " << myid << "\n"
+			  << "solution\n" << pmesh << x_phi << flush
+			  << "window_title 'Eigenmode " << i+1 << '/' << nev
+			  << ", Lambda = " << eigenvalues_phi[iphi] << "' "
+			  << "window_geometry 0 375 400 350\n"
+			  << flush;
+
+            curl_phi_sock << "parallel " << num_procs << " " << myid << "\n"
+			  << "solution\n" << pmesh << dx_rz << flush
+			  << "window_title 'Curl of Eigenmode " << i+1
+			  << '/' << nev
+			  << ", Lambda = " << eigenvalues_phi[iphi] << "' "
+			  << "window_geometry 400 375 400 350\n";
+	    if (iphi == 0)
+	    {
+	      curl_phi_sock << "keys vvv\n";
+	    }
+	    curl_phi_sock << flush;
+	    
+            iphi++;
+         }
+         char c;
+         if (myid == 0)
+         {
+            cout << "press (q)uit or (c)ontinue --> " << flush;
+            cin >> c;
+         }
+         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      mode_rz_sock.close();
+      curl_rz_sock.close();
+      mode_phi_sock.close();
+      curl_phi_sock.close();
+   }
+   /*
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mode_sock(vishost, visport);
+      mode_sock.precision(8);
+      socketstream curl_sock(vishost, visport);
+      curl_sock.precision(8);
+
+      for (int i=0; i<nev; i++)
+      {
+         if ( myid == 0 )
+         {
+            cout << "Eigenmode " << i+1 << '/' << nev
+                 << ", Lambda = " << eigenvalues_phi[i] << endl;
+         }
+
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x_phi = lobpcg->GetEigenvector(i);
+
+         curl2.Mult(x_phi, dx_rz);
+    // dx_l2.ProjectCoefficient(dxCoef);
+
+         mode_sock << "parallel " << num_procs << " " << myid << "\n"
+                   << "solution\n" << pmesh << x_phi << flush
+                   << "window_title 'Eigenmode " << i+1 << '/' << nev
+                   << ", Lambda = " << eigenvalues_phi[i] << "' "
+                   << "keys vvv\n" << flush;
+
+         // Limitations in the GridFunction and GLVis prevent this from working
+         curl_sock << "parallel " << num_procs << " " << myid << "\n"
+                   << "solution\n" << pmesh << dx_rz << flush
+                   << "window_title 'Curl of Eigenmode " << i+1 << '/' << nev
+                   << ", Lambda = " << eigenvalues_phi[i] << "' "
+                   << "window_geometry 400 0 400 350\n" << flush;
+
+         char c;
+         if (myid == 0)
+         {
+            cout << "press (q)uit or (c)ontinue --> " << flush;
+            cin >> c;
+         }
+         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      mode_sock.close();
+   }
+   */
+   // 12. Free the used memory.
+   delete ame;
+   delete ams;
+   delete lobpcg;
+   delete amg;
+   delete Mrz;
+   delete Arz;
+   delete Mphi;
+   delete Aphi;
+
+   MPI_Finalize();
+
+   return 0;
+}
@@ -0,0 +1,336 @@
+//                       MFEM Example 13-cyl - Parallel Version
+//
+// Compile with: make ex13p-cyl
+//
+// Sample runs:  mpirun -np 4 ex13p-cyl
+//               mpirun -np 4 ex13p-cyl -o 2
+//               mpirun -np 4 ex13p-cyl -o 2 -e 0
+//
+// Description:  This example code demonstrates the use of MFEM to solve the
+//               Maxwell (electromagnetic) eigenvalue problem on an
+//               axisymmetric domain. The eigenvalue problem:
+//                 curl curl E = lambda E
+//               with homogeneous Dirichlet boundary conditions E x n = 0 is
+//               solved on a cylindrical domain by meshing only a rectangle in
+//               the rho, z plane.  In cylindrical coordinates the weak form of
+//               the eigenvalue problem is given by:
+//                  (rho Curl(u), Curl(v)) = lambda (rho u, v)
+//
+//               We compute the five lowest nonzero eigenmodes by discretizing
+//               the curl curl operator using a Nedelec FE space of the
+//               specified order and compare to the known values.  Because the
+//               eigenvalue spectrum of a domain is unique this provides a
+//               reliable test that the axisymmetric domain is being faithfully
+//               characterized.
+//
+//               The example highlights the use of specialized coefficients
+//               with existing operators to mimic axisymmetric domains.  The
+//               curl of each eigenmode is also computed and displayed to
+//               ilustrate that no special steps need to be taken to compute
+//               the curl in this coordinate system.
+//
+//               We recommend viewing examples 13 and 11-cyl before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Zeros of Bessel function J_0
+static double J0z[] = {2.40482555769577,
+                       5.52007811028631,
+                       8.65372791291101,
+                       11.7915344390143
+                      };
+
+// Modes numbers in the rho and z directions for the first five eigenmodes
+static int mode_nums[] = {0, 0,
+                          0, 1,
+                          1, 0,
+                          1, 1,
+                          0, 2
+                         };
+
+double rhoFunc(const Vector &x)
+{
+   return x[0];
+}
+
+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.
+   int dim = 2;
+   int nr = 1;
+   int nz = 1;
+   int el_type_flag = 1;
+   Element::Type el_type;
+   int ser_ref_levels = 2;
+   int par_ref_levels = 1;
+   int order = 1;
+   int nev = 5;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&nz, "-nz", "--num-elements-z",
+                  "Number of elements in z-direction.");
+   args.AddOption(&nr, "-nr", "--num-elements-rho",
+                  "Number of elements in radial direction.");
+   args.AddOption(&el_type_flag, "-e", "--element-type",
+                  "Element type: 0 - Triangle, 1 - Quadrilateral.");
+   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",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   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);
+   }
+
+   // The output mesh could be quadrilaterals or triangles
+   el_type = (el_type_flag == 0) ? Element::TRIANGLE : Element::QUADRILATERAL;
+   if (el_type != Element::TRIANGLE && el_type != Element::QUADRILATERAL)
+   {
+      cout << "Unsupported element type" << endl;
+      exit(1);
+   }
+
+   // 3. Prepare a rectangular mesh with the desired dimensions and element
+   //    type.  Other 2D meshes could be used but then we couldn't check the
+   //    eigenvalues.
+   ParMesh pmesh;
+   {
+      Mesh mesh = Mesh::MakeCartesian2D(nr, nz, el_type);
+
+      // 4. Refine the serial mesh on all processors to increase the resolution. In
+      //    this example we do 'ref_levels' of uniform refinement (2 by default, or
+      //    specified on the command line with -rs).
+      for (int lev = 0; lev < ser_ref_levels; lev++)
+      {
+         mesh.UniformRefinement();
+      }
+
+      // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+      //    this mesh further in parallel to increase the resolution (1 time by
+      //    default, or specified on the command line with -rp). Once the parallel
+      //    mesh is defined, the serial mesh can be deleted.
+      pmesh = ParMesh(MPI_COMM_WORLD, mesh);
+      for (int lev = 0; lev < par_ref_levels; lev++)
+      {
+         pmesh.UniformRefinement();
+      }
+      pmesh.ReorientTetMesh();
+   }
+
+   // 6. Define a parallel finite element space on the parallel mesh. Here we
+   //    use the Nedelec finite elements (ND) of the specified order. We also
+   //    create an L2 space to represent the z-component of the curl of the
+   //    modes.
+   ND_FECollection fec_nd(order, dim);
+   L2_FECollection fec_rt(order - 1, dim,
+                          BasisType::GaussLegendre, FiniteElement::INTEGRAL);
+   L2_FECollection fec_l2(order - 1, dim);
+   ParFiniteElementSpace fespace_nd(&pmesh, &fec_nd);
+   ParFiniteElementSpace fespace_rt(&pmesh, &fec_rt);
+   ParFiniteElementSpace fespace_l2(&pmesh, &fec_l2);
+   HYPRE_Int size = fespace_nd.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of unknowns: " << size << endl;
+   }
+
+   // 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
+   //    element space. The first corresponds to the curl curl, while the second
+   //    is a simple mass matrix needed on the right hand side of the
+   //    generalized eigenvalue problem below. The boundary conditions are
+   //    implemented by marking all the boundary attributes from the mesh as
+   //    essential. The corresponding degrees of freedom are eliminated with
+   //    special values on the diagonal to shift the Dirichlet eigenvalues out
+   //    of the computational range. After serial and parallel assembly we
+   //    extract the corresponding parallel matrices A and M.
+   FunctionCoefficient rhoCoef(rhoFunc);
+   Array<int> ess_bdr(pmesh.bdr_attributes.Size());
+   ess_bdr = 1;
+   ess_bdr[3] = 0;
+
+   ParBilinearForm *a = new ParBilinearForm(&fespace_nd);
+   a->AddDomainIntegrator(new CurlCurlIntegrator(rhoCoef));
+   a->Assemble();
+   a->EliminateEssentialBCDiag(ess_bdr, 1.0);
+   a->Finalize();
+
+   ParBilinearForm *m = new ParBilinearForm(&fespace_nd);
+   m->AddDomainIntegrator(new VectorFEMassIntegrator(rhoCoef));
+   m->Assemble();
+   // shift the eigenvalue corresponding to eliminated dofs to a large value
+   m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
+   m->Finalize();
+
+   HypreParMatrix *A = a->ParallelAssemble();
+   HypreParMatrix *M = m->ParallelAssemble();
+
+   delete a;
+   delete m;
+
+   // 8. Define and configure the AME eigensolver and the AMS preconditioner for
+   //    A to be used within the solver. Set the matrices which define the
+   //    generalized eigenproblem A x = lambda M x.
+   HypreAMS *ams = new HypreAMS(*A,&fespace_nd);
+   ams->SetPrintLevel(0);
+   ams->SetSingularProblem();
+
+   HypreAME *ame = new HypreAME(MPI_COMM_WORLD);
+   ame->SetNumModes(nev);
+   ame->SetPreconditioner(*ams);
+   ame->SetMaxIter(100);
+   ame->SetTol(1e-8);
+   ame->SetPrintLevel(1);
+   ame->SetMassMatrix(*M);
+   ame->SetOperator(*A);
+
+   // 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
+   //    parallel grid function to represent each of the eigenmodes returned by
+   //    the solver. Also, define a discrete curl operator.
+   Array<double> eigenvalues;
+   ame->Solve();
+   ame->GetEigenvalues(eigenvalues);
+   ParGridFunction x(&fespace_nd);
+   ParGridFunction dx(&fespace_rt);
+
+   ParDiscreteLinearOperator curl(&fespace_nd, &fespace_rt);
+   curl.AddDomainInterpolator(new CurlInterpolator());
+   curl.Assemble();
+   curl.Finalize();
+
+   // This is one workaround for GLVis limitations
+   ParGridFunction dx_l2(&fespace_l2);
+
+   GridFunctionCoefficient dxCoef(&dx);
+
+   if ( myid == 0 )
+   {
+      cout << "\nRelative error in eigenvalues:\n";
+      for (int i=0; i<nev; i++)
+      {
+         double lambda =
+            pow(J0z[mode_nums[2*i]], 2) +
+            pow(M_PI * mode_nums[2*i+1], 2);
+         cout << "Lambda " << i+1 << '/' << nev << " = " << eigenvalues[i]
+              << ", rel err = " << fabs(eigenvalues[i] - lambda) / lambda
+              << endl;
+      }
+      cout << endl;
+   }
+
+   // 10. Save the refined mesh and the modes in parallel. This output can be
+   //     viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
+   {
+      ostringstream mesh_name, mode_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh.Print(mesh_ofs);
+
+      for (int i=0; i<nev; i++)
+      {
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = ame->GetEigenvector(i);
+
+         mode_name << "mode_" << setfill('0') << setw(2) << i << "."
+                   << setfill('0') << setw(6) << myid;
+
+         ofstream mode_ofs(mode_name.str().c_str());
+         mode_ofs.precision(8);
+         x.Save(mode_ofs);
+         mode_name.str("");
+      }
+   }
+
+   // 11. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mode_sock(vishost, visport);
+      mode_sock.precision(8);
+      socketstream curl_sock(vishost, visport);
+      curl_sock.precision(8);
+
+      for (int i=0; i<nev; i++)
+      {
+         if ( myid == 0 )
+         {
+            cout << "Eigenmode " << i+1 << '/' << nev
+                 << ", Lambda = " << eigenvalues[i] << endl;
+         }
+
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = ame->GetEigenvector(i);
+
+         curl.Mult(x, dx);
+         dx_l2.ProjectCoefficient(dxCoef);
+
+         mode_sock << "parallel " << num_procs << " " << myid << "\n"
+                   << "solution\n" << pmesh << x << flush
+                   << "window_title 'Eigenmode " << i+1 << '/' << nev
+                   << ", Lambda = " << eigenvalues[i] << "' "
+                   << "keys vvv\n" << flush;
+
+         // Limitations in the GridFunction and GLVis prevent this from working
+         curl_sock << "parallel " << num_procs << " " << myid << "\n"
+                   << "solution\n" << pmesh << dx_l2 << flush
+                   << "window_title 'Curl of Eigenmode " << i+1 << '/' << nev
+                   << ", Lambda = " << eigenvalues[i] << "' "
+                   << "window_geometry 400 0 400 350\n" << flush;
+
+         char c;
+         if (myid == 0)
+         {
+            cout << "press (q)uit or (c)ontinue --> " << flush;
+            cin >> c;
+         }
+         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      mode_sock.close();
+   }
+
+   // 12. Free the used memory.
+   delete ame;
+   delete ams;
+   delete M;
+   delete A;
+
+   MPI_Finalize();
+
+   return 0;
+}
@@ -0,0 +1,975 @@
+//                       MFEM Example 1 Ortho - Parallel Version
+//
+// Compile with: make ex1p
+//
+// Sample runs:  mpirun -np 4 ex1p-orth
+//               mpirun -np 4 ex1p-orth -c 2
+//               mpirun -np 4 ex1p-orth -c 3
+//               mpirun -np 4 ex1p-orth -c 4 -rs 4
+//               mpirun -np 4 ex1p-orth -c 6
+//               mpirun -np 4 ex1p-orth -c 7
+//               mpirun -np 4 ex1p-orth -c 8
+//               mpirun -np 4 ex1p-orth -c 9
+//               mpirun -np 4 ex1p-orth -c 10 -rs 4
+//               mpirun -np 4 ex1p-orth -c 11 -n2 2 -rs 3
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions in
+//               a variety of orthogonal coordinate systems.  The discretization
+//               is identical to that used in example 1 but here we use
+//               non-trivial coefficients in the Laplace operator and the right-
+//               hand-side vector to mimic a curvilinear coordinate system.  We
+//               also transform the mesh and solve the standard Laplace problem
+//               on the transformed mesh to compare the solutions.
+//
+//               The example highlights the use of standard differential
+//               operators to mimic the behavior of more exotic operators
+//               derived from coordinate transformations.
+//
+//               We recommend viewing Example 1 and Example 11-cyl before
+//               viewing this example.
+//
+//               Note: the notation used in this code comes from the Wikipedia
+//               page https://en.wikipedia.org/wiki/Orthogonal_coordinates.
+//               There are, however, minor differences made to ensure that
+//               we use right-handed coordinate systems in all cases.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Enumeration listing all supported 2D orthogonal coordinate systems
+enum CoordSys {POLAR = 1, PARABOLIC_CYL, ELLIPTIC, BIPOLAR,
+               CYLINDRICAL, SPHERICAL, PARABOLIC, PROLATE_SPHEROIDAL,
+               OBLATE_SPHEROIDAL, TOROIDAL, BISPHERICAL
+              };
+
+static CoordSys coords_ = (CoordSys)1;
+static double q1_min_ = NAN;
+static double q1_max_ = NAN;
+static double q2_min_ = NAN;
+static double q2_max_ = NAN;
+static double a_ = 1.0;
+
+// Set default values for coordinate ranges q1_min_, q1_max_, q2_min_,
+// and q2_max_ based on the selected coordinate system, coords_.
+void SetRanges();
+
+// Shift the mesh so that the origin is at (q1_min_, q2_min_)
+void trans1(const Vector &u, Vector &x)
+{
+   x.SetSize(2);
+   x[0] = u[0] + q1_min_;
+   x[1] = u[1] + q2_min_;
+}
+
+// Apply conformal mapping from cartesian coordinates to the
+// orthogonal coordinate system specified by coords_.
+void trans(const Vector &u, Vector &x);
+
+// Returns one of the three coordinate scale factors h_i describing
+// the orthogonal coordinate system.
+class OrthoCoef : public Coefficient
+{
+private:
+   int ind_;
+
+public:
+   OrthoCoef(int index) : ind_(index) {}
+
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip);
+};
+
+// Integration weight coefficient h_1 * h_2 * h_3
+class OrthoWeightCoef : public Coefficient
+{
+private:
+   Coefficient &h1Coef_;
+   Coefficient &h2Coef_;
+   Coefficient &h3Coef_;
+
+public:
+   OrthoWeightCoef(Coefficient &h1Coef,
+                   Coefficient &h2Coef,
+                   Coefficient &h3Coef)
+      : h1Coef_(h1Coef),
+        h2Coef_(h2Coef),
+        h3Coef_(h3Coef) {}
+
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip);
+};
+
+// Matrix-valued coefficient appearing in the weak form of the
+// Laplacian operator.
+class OrthoMatrixCoef : public MatrixCoefficient
+{
+private:
+   Coefficient &h1Coef_;
+   Coefficient &h2Coef_;
+   Coefficient &h3Coef_;
+
+public:
+   OrthoMatrixCoef(Coefficient &h1Coef,
+                   Coefficient &h2Coef,
+                   Coefficient &h3Coef)
+      : MatrixCoefficient(2),
+        h1Coef_(h1Coef),
+        h2Coef_(h2Coef),
+        h3Coef_(h3Coef)
+   {}
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+
+};
+
+// Radial weight factor to distinguish volumes of revolution from
+// extruded volumes
+class RhoCoef : public Coefficient
+{
+public:
+   RhoCoef() {}
+
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   {
+      double u_data[2];
+      Vector u(u_data, 2);
+      T.Transform(ip, u);
+      return u[0];
+   }
+
+};
+
+// Matrix-valued radial weight factor to distinguish volumes of
+// revolution from extruded volumes within the Laplacian operator
+class RhoMatrixCoef : public MatrixCoefficient
+{
+public:
+   RhoMatrixCoef()
+      : MatrixCoefficient(2)
+   {}
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      double u_data[2];
+      Vector u(u_data, 2);
+      T.Transform(ip, u);
+
+      K.SetSize(2);
+      K(0,0) = u[0];
+      K(0,1) = 0.0;
+      K(1,0) = 0.0;
+      K(1,1) = u[0];
+   }
+
+};
+
+static bool static_cond_ = false;
+static bool pa_ = false;
+
+// Ensure that m >= 3 if a periodic mesh has been selected
+void AdjustDimensions(int &m, int &n, int & rs, int & rp);
+
+// Setup and solve the Poisson problem with boundary conditions
+// appropriate to the selected coordinate system.
+void Poisson(ParMesh &pmesh, ParFiniteElementSpace &fespace,
+             MatrixCoefficient &LCoef, Coefficient &MCoef,
+             ParGridFunction &x);
+
+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.
+   int coords = 1;
+   int n1 = 1;
+   int n2 = 1;
+   int el_type_flag = 1;
+   Element::Type el_type;
+   int ser_ref_levels = 2;
+   int par_ref_levels = 1;
+   int morder = 2;
+   int order = 2;
+   const char *device_config = "cpu";
+   bool comp = true;
+   bool discont = false;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&coords, "-c", "--coord-sys",
+                  "Coordinate system: 1 - POLAR, 2 - PARABOLIC_CYL, "
+                  "3 - ELLIPTIC, 4 - BIPOLAR, 5 - CYLINDRICAL, 6 - SPHERICAL, "
+                  "7 - PARABOLIC, 8 - PROLATE_SPHEROIDAL, "
+                  "9 - OBLATE_SPHEROIDAL, 10 - TOROIDAL, 11 - BISPHERICAL");
+   args.AddOption(&n1, "-n1", "--num-elements-1",
+                  "Number of elements in q1-direction.");
+   args.AddOption(&n2, "-n2", "--num-elements-2",
+                  "Number of elements in q2-direction.");
+   args.AddOption(&q1_min_, "-q1-min", "--q1-minimum-1",
+                  "Minimum value of q1 coordinate.");
+   args.AddOption(&q1_max_, "-q1-max", "--q1-maximum-1",
+                  "Maximum value of q1 coordinate.");
+   args.AddOption(&q2_min_, "-q2-min", "--q2-minimum-1",
+                  "Minimum value of q2 coordinate.");
+   args.AddOption(&q2_max_, "-q2-max", "--q2-maximum-1",
+                  "Maximum value of q2 coordinate.");
+   args.AddOption(&a_, "-a", "--scale-parameter",
+                  "Scale paramter appearing in some of the transformations.");
+   args.AddOption(&el_type_flag, "-e", "--element-type",
+                  "Element type: 0 - Triangle, 1 - Quadrilateral.");
+   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(&morder, "-mo", "--mesh-order",
+                  "Order (polynomial degree) for the mesh geometry.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&comp, "-comp", "--compare", "-no-comp",
+                  "--no-compare", "Compare to standard curved mesh solution.");
+   args.AddOption(&static_cond_, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa_, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   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);
+   }
+
+   // Cast the user input to the enumerated type and set the appropriate
+   // coordinate ranges.
+   coords_ = (CoordSys)coords;
+   SetRanges();
+
+   // The output mesh could be quadrilaterals or triangles
+   el_type = (el_type_flag == 0) ? Element::TRIANGLE : Element::QUADRILATERAL;
+   if (el_type != Element::TRIANGLE && el_type != Element::QUADRILATERAL)
+   {
+      cout << "Unsupported element type" << endl;
+      exit(1);
+   }
+
+   if (coords_ == BIPOLAR || coords_ == TOROIDAL)
+   {
+      AdjustDimensions(n1, n2, ser_ref_levels, par_ref_levels);
+   }
+   else if (coords_ == POLAR || coords_ == ELLIPTIC)
+   {
+      AdjustDimensions(n2, n1, ser_ref_levels, par_ref_levels);
+   }
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (myid == 0) { device.Print(); }
+
+   // 3. Prepare a rectangular mesh with the desired dimensions and element
+   //    type.
+   Mesh *mesh = new Mesh(n1, n2, el_type, false,
+                         q1_max_ - q1_min_, q2_max_ - q2_min_);
+   mesh->Transform(trans1);
+   int dim = mesh->Dimension();
+
+   if (coords_ == POLAR || coords_ == ELLIPTIC || coords_ == BIPOLAR ||
+       coords_ == TOROIDAL)
+   {
+      // 4. Stitch the ends of the mesh together
+      discont = true;
+      mesh->SetCurvature(1, discont, 2, Ordering::byVDIM);
+      Array<int> v2v(mesh->GetNV());
+      for (int i = 0; i < v2v.Size(); i++)
+      {
+         v2v[i] = i;
+      }
+
+      if (coords_ == POLAR || coords_ == ELLIPTIC)
+      {
+         // identify vertices at the extremes of the mesh in the q2 direction
+         for (int i=0; i<n1 + 1; i++)
+         {
+            v2v[v2v.Size() - n1 - 1 + i] = i;
+         }
+      }
+      else if (coords_ == BIPOLAR || coords_ == TOROIDAL)
+      {
+         // identify vertices at the extremes of the mesh in the q1 direction
+         for (int i=0; i<n2 + 1; i++)
+         {
+            v2v[(n1 + 1) * i + n1] = (n1 + 1) * i;
+         }
+      }
+      // renumber elements
+      for (int i = 0; i < mesh->GetNE(); i++)
+      {
+         Element *el = mesh->GetElement(i);
+         int *v = el->GetVertices();
+         int nv = el->GetNVertices();
+         for (int j = 0; j < nv; j++)
+         {
+            v[j] = v2v[v[j]];
+         }
+      }
+      // renumber boundary elements
+      for (int i = 0; i < mesh->GetNBE(); i++)
+      {
+         Element *el = mesh->GetBdrElement(i);
+         int *v = el->GetVertices();
+         int nv = el->GetNVertices();
+         for (int j = 0; j < nv; j++)
+         {
+            v[j] = v2v[v[j]];
+         }
+      }
+      mesh->RemoveUnusedVertices();
+      mesh->RemoveInternalBoundaries();
+      mesh->FinalizeTopology();
+   }
+
+   // 5. Refine the serial mesh on all processors to increase the resolution.
+   {
+      for (int l = 0; l < ser_ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   // 6. 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_ortho = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   {
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh_ortho->UniformRefinement();
+      }
+   }
+
+   // 7. Create a standard curved mesh to describe the same geometry
+   ParMesh *pmesh_curved = new ParMesh(*pmesh_ortho);
+   pmesh_curved->SetCurvature(morder, discont);
+   pmesh_curved->Transform(trans);
+
+   // 8. Define a parallel finite element space on the parallel mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order.
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fespace_ortho(pmesh_ortho, &fec);
+   HYPRE_Int size = fespace_ortho.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 9. Declare the coordinate scaling factors and the coefficients
+   //    needed to form the mass matrix and Laplacian.
+   OrthoCoef h1Coef(0);
+   OrthoCoef h2Coef(1);
+   OrthoCoef h3Coef(2);
+   OrthoWeightCoef WCoef(h1Coef, h2Coef, h3Coef);
+   OrthoMatrixCoef LCoef(h1Coef, h2Coef, h3Coef);
+
+   // 10. Setup and solve the Poisson problem on the cartesian mesh
+   ParGridFunction x_ortho(&fespace_ortho); x_ortho = 0.0;
+   Poisson(*pmesh_ortho, fespace_ortho, LCoef, WCoef, x_ortho);
+
+   // 11. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh_ortho." << setfill('0') << setw(6) << myid;
+      sol_name << "sol_ortho." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh_ortho->Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+      x_ortho.Save(sol_ofs);
+   }
+
+   // 12. Send the solution by socket to a GLVis server.
+   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_ortho << x_ortho << flush
+               << "window_title 'Straight Mesh'"
+               << "keys m\n";
+
+      MPI_Barrier(MPI_COMM_WORLD);
+      socketstream mix_sol_sock(vishost, visport);
+      mix_sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      mix_sol_sock.precision(8);
+      mix_sol_sock << "solution\n" << *pmesh_curved << x_ortho << flush
+                   << "window_title 'Straight Solution - Curved Mesh' "
+                   << "window_geometry 400 0 400 350"
+                   << "keys m\n";
+   }
+
+   // 13. Compare to a solution computed on the corresponding curved mesh.
+   if (comp)
+   {
+      ParFiniteElementSpace fespace_curved(pmesh_curved, &fec);
+      ParGridFunction x_curved(&fespace_curved); x_curved = 0.0;
+
+      double err = -1.0;
+      GridFunctionCoefficient xCoef(&x_ortho);
+
+      // 14. Setup and solve the Poisson problem on the cartesian mesh
+      if (coords_ == POLAR || coords_ == PARABOLIC_CYL ||
+          coords_ == ELLIPTIC || coords_ == BIPOLAR)
+      {
+         // These coordinate systems can be viewed as truly two-dimensional
+         // or simply extruded into the third dimension and so they require
+         // no special coefficients.
+         DenseMatrix OneMat(2);
+         OneMat = 0.0; OneMat(0,0) = 1.0; OneMat(1,1) = 1.0;
+         MatrixConstantCoefficient OneCoef(OneMat);
+         ConstantCoefficient oneCoef(1.0);
+         Poisson(*pmesh_curved, fespace_curved, OneCoef, oneCoef, x_curved);
+
+         // 15a. Measure the difference in the two solutions using an L2 norm.
+         err = x_curved.ComputeL2Error(xCoef);
+      }
+      else
+      {
+         // The remaining coordinate systems are truly three-dimensional and
+         // involve rotation about the second coordinate axis. Consequently,
+         // they require a radial scale factor both in the mass matrix and the
+         // Laplacian operator.
+         RhoCoef rhoCoef;
+         RhoMatrixCoef RhoCoef;
+         Poisson(*pmesh_curved, fespace_curved, RhoCoef, rhoCoef, x_curved);
+
+         // 15b. Measure the difference in the two solutions using an L2 norm.
+         PowerCoefficient sqrtRhoCoef(rhoCoef, 0.5);
+         ProductCoefficient rxoCoef(sqrtRhoCoef, xCoef);
+         GridFunctionCoefficient xcCoef(&x_curved);
+         ProductCoefficient rxcCoef(sqrtRhoCoef, xcCoef);
+         ParGridFunction rx_curved(&fespace_curved);
+         rx_curved.ProjectCoefficient(rxcCoef);
+         err = rx_curved.ComputeL2Error(rxoCoef);
+      }
+
+      if (myid == 0)
+      {
+         cout << "\n|| u_curved - u_ortho ||_{L^2} = " << err << '\n' << endl;
+      }
+
+      // 15. Save the refined mesh and the solution in parallel. This output can
+      //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+      {
+         ostringstream mesh_name, sol_name;
+         mesh_name << "mesh_std." << setfill('0') << setw(6) << myid;
+         sol_name << "sol_std." << setfill('0') << setw(6) << myid;
+
+         ofstream mesh_ofs(mesh_name.str().c_str());
+         mesh_ofs.precision(8);
+         pmesh_curved->Print(mesh_ofs);
+
+         ofstream sol_ofs(sol_name.str().c_str());
+         sol_ofs.precision(8);
+         x_curved.Save(sol_ofs);
+      }
+
+      // 16. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+
+         socketstream cart_sol_sock(vishost, visport);
+         cart_sol_sock << "parallel " << num_procs << " " << myid << "\n";
+         cart_sol_sock.precision(8);
+         cart_sol_sock << "solution\n" << *pmesh_curved << x_curved << flush
+                       << "window_title 'Curved Mesh' "
+                       << "window_geometry 800 0 400 350"
+                       << "keys m\n";
+      }
+
+      delete pmesh_curved;
+   }
+
+   // 17. Free the used memory.
+   delete pmesh_ortho;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+void AdjustDimensions(int &m, int &n, int & rs, int & rp)
+{
+   while (m < 3 && rs + rp > 0)
+   {
+      m *= 2;
+      n *= 2;
+      (rs > 0) ? rs-- : rp--;
+   }
+   if (m < 3) { m = 3; }
+}
+
+void Poisson(ParMesh &pmesh, ParFiniteElementSpace &fespace,
+             MatrixCoefficient &LCoef, Coefficient &MCoef,
+             ParGridFunction &x)
+{
+   // 8. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    by marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet) and converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (pmesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      if (coords_ == CYLINDRICAL)
+      {
+         ess_bdr[3] = 0;
+      }
+      else if (coords_ == SPHERICAL || coords_ == PROLATE_SPHEROIDAL ||
+               coords_ == OBLATE_SPHEROIDAL)
+      {
+         ess_bdr[0] = 0;
+         ess_bdr[2] = 0;
+      }
+      else if (coords_ == PARABOLIC)
+      {
+         ess_bdr[0] = 0;
+         ess_bdr[3] = 0;
+      }
+      else if (coords_ == BISPHERICAL)
+      {
+         ess_bdr[1] = 0;
+         ess_bdr[3] = 0;
+      }
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 9. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (1,phi_i) where phi_i are the basis functions in fespace.
+   ParLinearForm b(&fespace);
+   b.AddDomainIntegrator(new DomainLFIntegrator(MCoef));
+   b.Assemble();
+
+   // 10. Define the solution vector x as a parallel finite element grid function
+   //     corresponding to fespace. Initialize x with initial guess of zero,
+   //     which satisfies the boundary conditions.
+   // ParGridFunction x(fespace);
+   // x = 0.0;
+
+   // 11. Set up the parallel bilinear form a(.,.) on the finite element space
+   //     corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   //     domain integrator.
+   ParBilinearForm a(&fespace);
+   if (pa_) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a.AddDomainIntegrator(new DiffusionIntegrator(LCoef));
+
+   // 12. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation, etc.
+   if (static_cond_) { a.EnableStaticCondensation(); }
+   a.Assemble();
+
+   OperatorPtr A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   // 13. Solve the linear system A X = B.
+   //     * With full assembly, use the BoomerAMG preconditioner from hypre.
+   //     * With partial assembly, use Jacobi smoothing, for now.
+   Solver *prec = NULL;
+   if (pa_)
+   {
+      if (UsesTensorBasis(fespace))
+      {
+         prec = new OperatorJacobiSmoother(a, ess_tdof_list);
+      }
+   }
+   else
+   {
+      prec = new HypreBoomerAMG;
+   }
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetRelTol(1e-12);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(1);
+   if (prec) { cg.SetPreconditioner(*prec); }
+   cg.SetOperator(*A);
+   cg.Mult(B, X);
+   delete prec;
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a.RecoverFEMSolution(X, b, x);
+}
+
+void SetRanges()
+{
+   switch (coords_)
+   {
+      case POLAR:
+         if (isnan(q1_min_)) { q1_min_ = 0.5; }
+         if (isnan(q1_max_)) { q1_max_ = 4.0; }
+         if (isnan(q2_min_)) { q2_min_ = -M_PI; }
+         if (isnan(q2_max_)) { q2_max_ =  M_PI; }
+         break;
+      case PARABOLIC_CYL:
+         if (isnan(q1_min_)) { q1_min_ = -4.0; }
+         if (isnan(q1_max_)) { q1_max_ = 4.0; }
+         if (isnan(q2_min_)) { q2_min_ = 0.5; }
+         if (isnan(q2_max_)) { q2_max_ = 4.0; }
+         break;
+      case ELLIPTIC:
+         if (isnan(q1_min_)) { q1_min_ = 0.5; }
+         if (isnan(q1_max_)) { q1_max_ = 2.0; }
+         if (isnan(q2_min_)) { q2_min_ = -M_PI; }
+         if (isnan(q2_max_)) { q2_max_ =  M_PI; }
+         break;
+      case BIPOLAR:
+         if (isnan(q1_min_)) { q1_min_ = -M_PI; }
+         if (isnan(q1_max_)) { q1_max_ =  M_PI; }
+         if (isnan(q2_min_)) { q2_min_ = 0.5; }
+         if (isnan(q2_max_)) { q2_max_ = 4.0; }
+         break;
+      case CYLINDRICAL:
+         if (isnan(q1_min_)) { q1_min_ = 0.0; }
+         if (isnan(q1_max_)) { q1_max_ = 4.0; }
+         if (isnan(q2_min_)) { q2_min_ = 0.0; }
+         if (isnan(q2_max_)) { q2_max_ = 4.0; }
+         break;
+      case SPHERICAL:
+         if (isnan(q1_min_)) { q1_min_ = 0.5; }
+         if (isnan(q1_max_)) { q1_max_ = 4.0; }
+         if (isnan(q2_min_)) { q2_min_ = 0.0; }
+         if (isnan(q2_max_)) { q2_max_ = M_PI; }
+         break;
+      case PARABOLIC:
+         if (isnan(q1_min_)) { q1_min_ = 0.2; }
+         if (isnan(q1_max_)) { q1_max_ = 4.0; }
+         if (isnan(q2_min_)) { q2_min_ = 0.2; }
+         if (isnan(q2_max_)) { q2_max_ = 4.0; }
+         break;
+      case PROLATE_SPHEROIDAL:
+         if (isnan(q1_min_)) { q1_min_ = 0.5; }
+         if (isnan(q1_max_)) { q1_max_ = 2.0; }
+         if (isnan(q2_min_)) { q2_min_ = 0.0; }
+         if (isnan(q2_max_)) { q2_max_ = M_PI; }
+         break;
+      case OBLATE_SPHEROIDAL:
+         if (isnan(q1_min_)) { q1_min_ = 0.5; }
+         if (isnan(q1_max_)) { q1_max_ = 2.0; }
+         if (isnan(q2_min_)) { q2_min_ = -0.5 * M_PI; }
+         if (isnan(q2_max_)) { q2_max_ =  0.5 * M_PI; }
+         break;
+      case TOROIDAL:
+         if (isnan(q1_min_)) { q1_min_ = -M_PI; }
+         if (isnan(q1_max_)) { q1_max_ =  M_PI; }
+         if (isnan(q2_min_)) { q2_min_ = 0.5; }
+         if (isnan(q2_max_)) { q2_max_ = 4.0; }
+         break;
+      case BISPHERICAL:
+         if (isnan(q1_min_)) { q1_min_ = 0.0; }
+         if (isnan(q1_max_)) { q1_max_ = M_PI; }
+         if (isnan(q2_min_)) { q2_min_ = 0.5; }
+         if (isnan(q2_max_)) { q2_max_ = 4.0; }
+         break;
+   }
+}
+
+void trans(const Vector &u, Vector &x)
+{
+   x.SetSize(2);
+
+   switch (coords_)
+   {
+      case POLAR:
+         x[0] = u[0] * cos(u[1]);
+         x[1] = u[0] * sin(u[1]);
+         break;
+      case PARABOLIC_CYL:
+         x[0] = 0.5 * (u[0] * u[0] - u[1] * u[1]);
+         x[1] = u[0] * u[1];
+         break;
+      case ELLIPTIC:
+         x[0] = a_ * cosh(u[0]) * cos(u[1]);
+         x[1] = a_ * sinh(u[0]) * sin(u[1]);
+         break;
+      case BIPOLAR:
+      {
+         double den = (cosh(u[1]) - cos(u[0]));
+         x[0] = a_ * sinh(u[1]) / den;
+         x[1] = a_ * sin(u[0]) / den;
+      }
+      break;
+      case CYLINDRICAL:
+      {
+         x[0] = u[0];
+         x[1] = u[1];
+      }
+      break;
+      case SPHERICAL:
+      {
+         x[0] = u[0] * sin(u[1]);
+         x[1] = -u[0] * cos(u[1]);
+      }
+      break;
+      case PARABOLIC:
+      {
+         x[0] = u[0] * u[1];
+         x[1] = -0.5 * (u[0] * u[0] - u[1] * u[1]);
+      }
+      break;
+      case PROLATE_SPHEROIDAL:
+      {
+         x[0] = a_ * sinh(u[0]) * sin(u[1]);
+         x[1] = -a_ * cosh(u[0]) * cos(u[1]);
+      }
+      break;
+      case OBLATE_SPHEROIDAL:
+      {
+         x[0] = a_ * cosh(u[0]) * cos(u[1]);
+         x[1] = a_ * sinh(u[0]) * sin(u[1]);
+      }
+      break;
+      case TOROIDAL:
+      {
+         double den = (cosh(u[1]) - cos(u[0]));
+         x[0] = a_ * sinh(u[1]) / den;
+         x[1] = a_ * sin(u[0]) / den;
+      }
+      break;
+      case BISPHERICAL:
+      {
+         double den = (cosh(u[1]) + cos(u[0]));
+         x[0] = a_ * sin(u[0]) / den;
+         x[1] = a_ * sinh(u[1]) / den;
+      }
+      break;
+   }
+}
+
+double OrthoCoef::Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+{
+   double u_data[2];
+   Vector u(u_data, 2);
+   T.Transform(ip, u);
+
+   switch (coords_)
+   {
+      case POLAR:
+         switch (ind_)
+         {
+            case 0:
+               return 1.0;
+            case 1:
+               return u[0];
+            case 2:
+               return 1.0;
+            default:
+               return 0.0;
+         }
+         break;
+      case PARABOLIC_CYL:
+         switch (ind_)
+         {
+            case 0:
+               return sqrt(u[0] * u[0] + u[1] * u[1]);
+            case 1:
+               return sqrt(u[0] * u[0] + u[1] * u[1]);
+            case 2:
+               return 1.0;
+            default:
+               return 0.0;
+         }
+         break;
+      case ELLIPTIC:
+         switch (ind_)
+         {
+            case 0:
+               return a_ * sqrt(pow(sinh(u[0]), 2) + pow(sin(u[1]), 2));
+            case 1:
+               return a_ * sqrt(pow(sinh(u[0]), 2) + pow(sin(u[1]), 2));
+            case 2:
+               return 1.0;
+            default:
+               return 0.0;
+         }
+         break;
+      case BIPOLAR:
+      {
+         double den = (cosh(u[1]) - cos(u[0]));
+
+         switch (ind_)
+         {
+            case 0:
+               return a_ / den;
+            case 1:
+               return a_ / den;
+            case 2:
+               return 1.0;
+            default:
+               return 0.0;
+         }
+         break;
+      }
+      case CYLINDRICAL:
+         switch (ind_)
+         {
+            case 0:
+               return 1.0;
+            case 1:
+               return 1.0;
+            case 2:
+               return u[0];
+            default:
+               return 0.0;
+         }
+         break;
+      case SPHERICAL:
+         switch (ind_)
+         {
+            case 0:
+               return 1.0;
+            case 1:
+               return u[0];
+            case 2:
+               return u[0] * sin(u[1]);
+            default:
+               return 0.0;
+         }
+         break;
+      case PARABOLIC:
+         switch (ind_)
+         {
+            case 0:
+               return sqrt(u[0] * u[0] + u[1] * u[1]);
+            case 1:
+               return sqrt(u[0] * u[0] + u[1] * u[1]);
+            case 2:
+               return u[0] * u[1];
+            default:
+               return 0.0;
+         }
+         break;
+      case PROLATE_SPHEROIDAL:
+         switch (ind_)
+         {
+            case 0:
+               return a_ * sqrt(pow(sinh(u[0]), 2) + pow(sin(u[1]), 2));
+            case 1:
+               return a_ * sqrt(pow(sinh(u[0]), 2) + pow(sin(u[1]), 2));
+            case 2:
+               return a_ * sinh(u[0]) * sin(u[1]);
+            default:
+               return 0.0;
+         }
+         break;
+      case OBLATE_SPHEROIDAL:
+         switch (ind_)
+         {
+            case 0:
+               return a_ * sqrt(pow(sinh(u[0]), 2) + pow(sin(u[1]), 2));
+            case 1:
+               return a_ * sqrt(pow(sinh(u[0]), 2) + pow(sin(u[1]), 2));
+            case 2:
+               return a_ * cosh(u[0]) * cos(u[1]);
+            default:
+               return 0.0;
+         }
+         break;
+      case TOROIDAL:
+      {
+         double den = (cosh(u[1]) - cos(u[0]));
+
+         switch (ind_)
+         {
+            case 0:
+               return a_ / den;
+            case 1:
+               return a_ / den;
+            case 2:
+               return a_ * sinh(u[1]) / den;
+            default:
+               return 0.0;
+         }
+         break;
+      }
+      case BISPHERICAL:
+      {
+         double den = (cosh(u[1]) + cos(u[0]));
+
+         switch (ind_)
+         {
+            case 0:
+               return a_ / den;
+            case 1:
+               return a_ / den;
+            case 2:
+               return a_ * sin(u[0]) / den;
+            default:
+               return 0.0;
+         }
+         break;
+      }
+   }
+   return 0.0;
+}
+
+double OrthoWeightCoef::Eval(ElementTransformation &T,
+                             const IntegrationPoint &ip)
+{
+   double h1 = h1Coef_.Eval(T, ip);
+   double h2 = h2Coef_.Eval(T, ip);
+   double h3 = h3Coef_.Eval(T, ip);
+
+   return h1 * h2 * h3;
+}
+
+void OrthoMatrixCoef::Eval(DenseMatrix &K, ElementTransformation &T,
+                           const IntegrationPoint &ip)
+{
+   double h1 = h1Coef_.Eval(T, ip);
+   double h2 = h2Coef_.Eval(T, ip);
+   double h3 = h3Coef_.Eval(T, ip);
+
+   K.SetSize(2);
+   K(0,0) = h2 * h3 / h1;
+   K(0,1) = 0.0;
+   K(1,0) = 0.0;
+   K(1,1) = h1 * h3 / h2;
+}
@@ -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 ex1p-orth ex11p-cyl ex13p-cyl ex13p-cyl-3d
 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
Author	SHA1	Message	Date
Stowell, Mark L.	33e1ca24b9	Merge remote-tracking branch 'origin/master' into cyl-coords-dev # Conflicts: # .gitignore # examples/CMakeLists.txt # examples/makefile	2025-10-20 14:08:39 -07:00
Stowell, Mark L	1d75ab10d3	Merge remote-tracking branch 'origin/master' into cyl-coords-dev # Conflicts: # examples/CMakeLists.txt # examples/makefile	2021-04-27 18:29:10 -07:00
Stowell, Mark L	a36b54327b	Moving definition of dim	2021-04-27 18:27:25 -07:00
Stowell, Mark L	c004978124	Adding curl curl example	2021-04-27 18:26:00 -07:00
Stowell, Mark L	580d6135b3	Adding viz of curl and fixing deprecated Mesh c'tor	2021-04-27 10:33:43 -07:00
Stowell, Mark L	027cf57246	Merge remote-tracking branch 'origin/master' into cyl-coords-dev # Conflicts: # examples/CMakeLists.txt # examples/makefile	2021-04-18 19:18:51 -07:00
Stowell, Mark L	a61901c706	Adjusting default range for parabolic coordinate system	2021-02-16 19:11:25 -08:00
Stowell, Mark L	85103bf9dd	expanding comments	2021-02-11 16:59:28 -08:00
Stowell, Mark L	c23df0143b	Merge remote-tracking branch 'origin/master' into cyl-coords-dev	2021-02-10 14:23:15 -08:00
Stowell, Mark L	634d91fa0b	Adding two more sample runs	2021-01-15 13:32:40 -08:00
Stowell, Mark L	21c4bdaadd	Merge remote-tracking branch 'origin/master' into cyl-coords-dev # Conflicts: # examples/CMakeLists.txt # examples/makefile	2021-01-14 11:40:12 -08:00
Stowell, Mark L	2abe24725d	make style	2020-04-07 16:58:02 -07:00
Stowell, Mark L	bc922ab30a	Rewriting comments and adding mesh geometry order option	2020-03-29 11:58:25 -07:00
Stowell, Mark L	f4fddd0e16	Merge remote-tracking branch 'origin/master' into cyl-coords-dev # Conflicts: # examples/CMakeLists.txt # examples/makefile	2020-03-27 13:59:16 -07:00
Stowell, Mark L	a086fdca97	Adding new ortho example to the build system	2020-03-27 13:53:25 -07:00
Stowell, Mark L	e894b0782e	Adding Poisson example using various orthogonal coordinate systems	2020-03-27 12:32:11 -07:00
Stowell, Mark L	8e6d39574e	Adding a temporary line to .gitignore to address the Travis tests	2020-03-25 16:13:13 -07:00
Stowell, Mark L	131f12fe24	Adding new examples to buildsystem	2020-03-25 15:48:15 -07:00
Stowell, Mark L	6253b48837	Cleaning up comments	2020-03-25 15:46:21 -07:00
Stowell, Mark L	4a2da44902	Adding example codes to demonstrate operators with axial symmetry	2020-03-25 08:46:40 -07:00