Fixing merge error

examples/CMakeLists.txt

@@ -38,6 +38,7 @@ list(APPEND ALL_EXE_SRCS
   ex28.cpp
   ex29.cpp
   ex30_proposed.cpp
+  ex32_proposed.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -73,6 +74,7 @@ if (MFEM_USE_MPI)
     ex29p.cpp
     ex30p_proposed.cpp
     ex31p_proposed.cpp
+    ex32p_proposed.cpp
     )
 endif()
 

examples/ex32_proposed.cpp

@@ -0,0 +1,300 @@
+//                                MFEM Example 32
+//
+// Compile with: make ex32
+//
+// Sample runs:  ex32 -m ../data/hexagon.mesh -o 2
+//               ex32 -m ../data/inline-tri.mesh -o 3
+//               ex32 -m ../data/amr-quad.mesh -o 3
+//
+// Description:  This is a version of Example 30 with a simple adaptive mesh
+//               refinement loop. The problem being solved is again the
+//               electromagnetic diffusion problem with an anisotropic
+//               conductivity coefficient. The problem is solved on a sequence
+//               of 2D meshes which are locally refined in a conforming
+//               (triangles) or non-conforming (quadrilaterals) manner according
+//               to a simple ZZ error estimator.
+//
+//               The example demonstrates MFEM's capability to work with both
+//               conforming and nonconforming refinements on 2D meshes.
+//               Interpolation of functions from coarse to fine meshes, as well
+//               as persistent GLVis visualization are also illustrated.
+//
+//               We recommend viewing Examples 6 and 30 before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static Vector bb_min;
+static Vector bb_max;
+void f_func(const Vector &, Vector &);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int max_amr_its = 100;
+   int max_dofs = 20000;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&max_amr_its, "-mx", "--max-amr-its",
+                  "Maximum number of AMR iterations.");
+   args.AddOption(&max_dofs, "-md", "--max-amr-dofs",
+                  "Maximum number of degrees of freedom.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+
+   // 2. Read the mesh from the given mesh file.  We can handle triangular
+   //    or quadrilateral meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   mesh.GetBoundingBox(bb_min, bb_max);
+
+   MFEM_VERIFY(dim == 2, "This exmaple requires a 2D mesh.");
+
+   // 3. Refine the mesh to increase the resolution. Also project a NURBS mesh
+   //    to a piecewise-quadratic curved mesh. Make sure that the mesh is
+   //    non-conforming.
+   if (mesh.NURBSext)
+   {
+      mesh.UniformRefinement();
+      mesh.SetCurvature(2);
+   }
+
+   // 4. Define a finite element space on the mesh. The polynomial order is
+   //    one (linear) by default, but this can be changed on the command line.
+   ND_R2D_FECollection fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec);
+
+   // 5. Set up the linear form b(.) which corresponds to the right-hand side
+   //    of the FEM linear system, which in this case is (f,phi_i) where f is
+   //    given by the function f_func and phi_i are the basis functions in the
+   //    finite element fespace.
+   VectorFunctionCoefficient f(3, f_func);
+   LinearForm b(&fespace);
+   b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
+
+   // 6. The solution vector x and the associated finite element grid function
+   //    will be maintained over the AMR iterations. We initialize it to zero.
+   Vector zeroVec(3); zeroVec = 0.0;
+   VectorConstantCoefficient zeroCoef(zeroVec);
+   GridFunction sol(&fespace);
+   sol = 0;
+
+   // 7. Set up the bilinear form corresponding to the EM diffusion operator
+   //    curl muinv curl + sigma I, by adding the curl-curl and the mass domain
+   //    integrators.
+   Array<int> ess_bdr(mesh.bdr_attributes.Max());
+   if (ess_bdr.Size() > 0) { ess_bdr = 1; }
+
+   DenseMatrix sigmaMat(3);
+   sigmaMat(0,0) = 2.0; sigmaMat(1,1) = 2.0; sigmaMat(2,2) = 2.0;
+   sigmaMat(0,2) = 0.0; sigmaMat(2,0) = 0.0;
+   sigmaMat(0,1) = M_SQRT1_2; sigmaMat(1,0) = M_SQRT1_2; // 1/sqrt(2) in cmath
+   sigmaMat(1,2) = M_SQRT1_2; sigmaMat(2,1) = M_SQRT1_2;
+
+   ConstantCoefficient muinv(1.0);
+   MatrixConstantCoefficient sigma(sigmaMat);
+   BilinearForm a(&fespace);
+   BilinearFormIntegrator * integ = new CurlCurlIntegrator(muinv);
+   a.AddDomainIntegrator(integ);
+   a.AddDomainIntegrator(new VectorFEMassIntegrator(sigma));
+
+   // 8. Connect to GLVis.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+
+   socketstream xy_sock, z_sock;
+   if (visualization)
+   {
+      xy_sock.open(vishost, visport);
+      z_sock.open(vishost, visport);
+      if (!xy_sock && !z_sock)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         cout << "GLVis visualization disabled.\n";
+
+         visualization = false;
+      }
+
+      xy_sock.precision(8);
+      z_sock.precision(8);
+   }
+
+   // 9. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
+   //    that uses the ComputeElementFlux method of the CurlCurlIntegrator to
+   //    recover a smoothed flux (curl) that is subtracted from the element
+   //    flux to get an error indicator. We need to supply a space for the
+   //    discontinuous flux (RT) and a space for the smoothed flux (H(curl) is
+   //    used here).
+   RT_R2D_FECollection flux_fec(order-1, dim);
+   FiniteElementSpace flux_fes(&mesh, &flux_fec);
+   ND_R2D_FECollection smooth_flux_fec(order, dim);
+   FiniteElementSpace smooth_flux_fes(&mesh, &smooth_flux_fec);
+   L2ZienkiewiczZhuEstimator estimator(*integ, sol, flux_fes, smooth_flux_fes);
+
+   // 10. A refiner selects and refines elements based on a refinement strategy.
+   //     The strategy here is to refine elements with errors larger than a
+   //     fraction of the maximum element error. Other strategies are possible.
+   //     The refiner will call the given error estimator.
+   ThresholdRefiner refiner(estimator);
+   refiner.SetTotalErrorFraction(0.6);
+
+   // 11. The main AMR loop. In each iteration we solve the problem on the
+   //     current mesh, visualize the solution, and refine the mesh.
+   for (int it = 0; it <= max_amr_its; it++)
+   {
+      int cdofs = fespace.GetTrueVSize();
+      cout << "\nAMR iteration " << it << endl;
+      cout << "Number of unknowns: " << cdofs << endl;
+
+      // 12. Assemble the right-hand side and determine the list of true
+      //     essential boundary dofs.
+      Array<int> ess_tdof_list;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      b.Assemble();
+
+      // 13. Assemble the stiffness matrix. Note that MFEM doesn't care at this
+      //     point that the mesh is nonconforming.  The FE space is considered
+      //     'cut' along hanging edges/faces.
+      a.Assemble();
+
+      // 14. Create the linear system: eliminate boundary conditions. The
+      //     system will be solved for true (unconstrained/unique) DOFs only.
+      OperatorPtr A;
+      Vector B, X;
+
+      sol.ProjectBdrCoefficientTangent(zeroCoef, ess_bdr);
+
+      const int copy_interior = 1;
+      a.FormLinearSystem(ess_tdof_list, sol, b, A, X, B, copy_interior);
+
+      // 15. Solve the linear system A X = B.
+#ifndef MFEM_USE_SUITESPARSE
+      // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
+      GSSmoother M((SparseMatrix&)(*A));
+      PCG(*A, M, B, X, 3, 500, 1e-12, 1e-12);
+#else
+      // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
+      UMFPackSolver umf_solver;
+      umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+      umf_solver.SetOperator(*A);
+      umf_solver.Mult(B, X);
+#endif
+
+      // 16. After solving the linear system, reconstruct the solution as a
+      //     finite element GridFunction. Constrained edges are interpolated
+      //     from true DOFs (it may therefore happen that x.Size() >= X.Size()).
+      a.RecoverFEMSolution(X, b, sol);
+
+      // 17. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         DenseMatrix xyMat(2,3); xyMat = 0.0;
+         xyMat(0,0) = 1.0; xyMat(1,1) = 1.0;
+         MatrixConstantCoefficient xyMatCoef(xyMat);
+         Vector zVec(3); zVec = 0.0; zVec(2) = 1;
+         VectorConstantCoefficient zVecCoef(zVec);
+
+         VectorGridFunctionCoefficient solCoef(&sol);
+         MatrixVectorProductCoefficient xyCoef(xyMatCoef, solCoef);
+         InnerProductCoefficient zCoef(zVecCoef, solCoef);
+
+         H1_FECollection fec_h1(order, dim);
+         ND_FECollection fec_nd(order, dim);
+
+         FiniteElementSpace fes_h1(&mesh, &fec_h1);
+         FiniteElementSpace fes_nd(&mesh, &fec_nd);
+
+         GridFunction xyComp(&fes_nd);
+         GridFunction zComp(&fes_h1);
+
+         xyComp.ProjectCoefficient(xyCoef);
+         zComp.ProjectCoefficient(zCoef);
+
+         xy_sock << "solution\n" << mesh << xyComp << flush;
+         if (it == 0)
+         {
+            xy_sock << "keys vvv "
+                    << "window_geometry 0 0 400 350 "
+                    << "window_title 'XY components'\n";
+         }
+
+         z_sock << "solution\n" << mesh << zComp << flush;
+         if (it == 0)
+         {
+            z_sock << "window_geometry 403 0 400 350 "
+                   << "window_title 'Z component'\n";
+         }
+      }
+
+      if (cdofs > max_dofs)
+      {
+         cout << "Reached the maximum number of dofs. Stop." << endl;
+         break;
+      }
+
+      // 20. Call the refiner to modify the mesh. The refiner calls the error
+      //     estimator to obtain element errors, then it selects elements to be
+      //     refined and finally it modifies the mesh. The Stop() method can be
+      //     used to determine if a stopping criterion was met.
+      refiner.Apply(mesh);
+      if (refiner.Stop())
+      {
+         cout << "Stopping criterion satisfied. Stop." << endl;
+         break;
+      }
+
+      // 21. Update the finite element space (recalculate the number of DOFs,
+      //     etc.) and create a grid function update matrix. Apply the matrix
+      //     to any GridFunctions over the space. In this case, the update
+      //     matrix is an interpolation matrix so the updated GridFunction will
+      //     still represent the same function as before refinement.
+      fespace.Update();
+      sol.Update();
+
+      // 22. Inform also the bilinear and linear forms that the space has
+      //     changed.
+      a.Update();
+      b.Update();
+   }
+   if (visualization)
+   {
+      xy_sock.close();
+      z_sock.close();
+   }
+
+   return 0;
+}
+
+void f_func(const Vector &x, Vector &f)
+{
+   double xc = 0.5 * (bb_min[0] + bb_max[0]);
+   double yc = 0.5 * (bb_min[1] + bb_max[1]);
+   double dx = bb_max[0] - bb_min[0];
+   double dy = bb_max[1] - bb_min[1];
+
+   f = 0.0;
+   if (fabs(x[0] - xc) < 0.2 * dx && fabs(x[1] - yc) < 0.2 * dy)
+   {
+      double a = pow(cos(2.5 * M_PI * (x[0] - xc) / dx) *
+                     cos(2.5 * M_PI * (x[1] - yc) / dy), 2);
+      f(0) = a * sin(2.5 * M_PI * (x[1] - yc) / dy);
+      f(1) = a * sin(5.0 * M_PI * (x[0] - xc) / dx);
+      f(2) = a * cos(5.0 * M_PI * (x[0] - xc) / dx);
+   }
+}

examples/ex32p_proposed.cpp

@@ -0,0 +1,340 @@
+//                       MFEM Example 32 - Parallel Version
+//
+// Compile with: make ex32p
+//
+// Sample runs:  mpirun -np 3 ex32p -m ../data/hexagon.mesh -o 2
+//               mpirun -np 4 ex32p -m ../data/inline-tri.mesh -o 3
+//               mpirun -np 4 ex32p -m ../data/amr-quad.mesh -o 3
+//
+// Description:  This is a version of Example 30 with a simple adaptive mesh
+//               refinement loop. The problem being solved is again the
+//               electromagnetic diffusion problem with an anisotropic
+//               conductivity coefficient. The problem is solved on a sequence
+//               of 2D meshes which are locally refined in a conforming
+//               (triangles) or non-conforming (quadrilaterals) manner according
+//               to a simple ZZ error estimator.
+//
+//               The example demonstrates MFEM's capability to work with both
+//               conforming and nonconforming refinements on 2D meshes.
+//               Interpolation of functions from coarse to fine meshes, as well
+//               as persistent GLVis visualization are also illustrated.
+//
+//               We recommend viewing Examples 6 and 30 before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static Vector bb_min;
+static Vector bb_max;
+void f_func(const Vector &, Vector &);
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   MPI_Session mpi;
+   int num_procs = mpi.WorldSize();
+   int myid = mpi.WorldRank();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/inline-quad.mesh";
+   int order = 1;
+   int max_amr_its = 100;
+   int max_dofs = 20000;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&max_amr_its, "-mx", "--max-amr-its",
+                  "Maximum number of AMR iterations.");
+   args.AddOption(&max_dofs, "-md", "--max-amr-dofs",
+                  "Maximum number of degrees of freedom.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+
+   // 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();
+
+   mesh->GetBoundingBox(bb_min, bb_max);
+
+   MFEM_VERIFY(dim == 2, "This example requires a 2D mesh.");
+
+   // 4. Refine the serial mesh on all processors to increase the resolution.
+   //    Also project a NURBS mesh to a piecewise-quadratic curved mesh. Make
+   //    sure that the mesh is non-conforming.
+   if (mesh->NURBSext)
+   {
+      mesh->UniformRefinement();
+      mesh->SetCurvature(2);
+   }
+   mesh->EnsureNCMesh();
+
+   // 5. Define a parallel mesh by partitioning the serial mesh.
+   //    Once the parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+
+   MFEM_VERIFY(pmesh.bdr_attributes.Size() > 0,
+               "Boundary attributes required in the mesh.");
+   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+   ess_bdr = 1;
+
+   // 6. Define a finite element space on the mesh. The polynomial order is
+   //    one (linear) by default, but this can be changed on the command line.
+   ND_R2D_FECollection fec(order, dim);
+   ParFiniteElementSpace fespace(&pmesh, &fec);
+
+   // 7. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (f,phi_i) where f is given by the function f_func and phi_i are the
+   //    basis functions in the finite element fespace.
+   VectorFunctionCoefficient f(3, f_func);
+   ParLinearForm b(&fespace);
+   b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
+
+   // 8. The solution vector x and the associated finite element grid function
+   //    will be maintained over the AMR iterations. We initialize it to zero.
+   Vector zeroVec(3); zeroVec = 0.0;
+   VectorConstantCoefficient zeroCoef(zeroVec);
+   ParGridFunction sol(&fespace);
+   sol = 0;
+
+   // 9. Set up the parallel bilinear form corresponding to the EM diffusion
+   //    operator curl muinv curl + sigma I, by adding the curl-curl and the
+   //    mass domain integrators.
+   DenseMatrix sigmaMat(3);
+   sigmaMat(0,0) = 2.0; sigmaMat(1,1) = 2.0; sigmaMat(2,2) = 2.0;
+   sigmaMat(0,2) = 0.0; sigmaMat(2,0) = 0.0;
+   sigmaMat(0,1) = M_SQRT1_2; sigmaMat(1,0) = M_SQRT1_2; // 1/sqrt(2) in cmath
+   sigmaMat(1,2) = M_SQRT1_2; sigmaMat(2,1) = M_SQRT1_2;
+
+   ConstantCoefficient muinv(1.0);
+   MatrixConstantCoefficient sigma(sigmaMat);
+   ParBilinearForm a(&fespace);
+   BilinearFormIntegrator * integ = new CurlCurlIntegrator(muinv);
+   a.AddDomainIntegrator(integ);
+   a.AddDomainIntegrator(new VectorFEMassIntegrator(sigma));
+
+   // 10. Connect to GLVis.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+
+   socketstream xy_sock, z_sock;
+   if (visualization)
+   {
+      xy_sock.open(vishost, visport);
+      z_sock.open(vishost, visport);
+      if (!xy_sock && !z_sock)
+      {
+         if (mpi.Root())
+         {
+            cout << "Unable to connect to GLVis server at "
+                 << vishost << ':' << visport << endl;
+            cout << "GLVis visualization disabled.\n";
+         }
+         visualization = false;
+      }
+
+      xy_sock.precision(8);
+      z_sock.precision(8);
+   }
+
+   // 11. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
+   //     with L2 projection in the smoothing step to better handle hanging
+   //     nodes and parallel partitioning. We need to supply a space for the
+   //     discontinuous flux (RT) and a space for the smoothed flux (H(curl) is
+   //     used here).
+   RT_R2D_FECollection flux_fec(order-1, dim);
+   ParFiniteElementSpace flux_fes(&pmesh, &flux_fec);
+   ND_R2D_FECollection smooth_flux_fec(order, dim);
+   ParFiniteElementSpace smooth_flux_fes(&pmesh, &smooth_flux_fec);
+   // Another possible option for the smoothed flux space:
+   // H1_FECollection smooth_flux_fec(order, dim);
+   // ParFiniteElementSpace smooth_flux_fes(&pmesh, &smooth_flux_fec, 3);
+   L2ZienkiewiczZhuEstimator estimator(*integ, sol, flux_fes, smooth_flux_fes);
+
+   // 12. A refiner selects and refines elements based on a refinement strategy.
+   //     The strategy here is to refine elements with errors larger than a
+   //     fraction of the maximum element error. Other strategies are possible.
+   //     The refiner will call the given error estimator.
+   ThresholdRefiner refiner(estimator);
+   refiner.SetTotalErrorFraction(0.6);
+
+   // 13. The main AMR loop. In each iteration we solve the problem on the
+   //     current mesh, visualize the solution, and refine the mesh.
+   for (int it = 0; it <= max_amr_its; it++)
+   {
+      HYPRE_Int global_dofs = fespace.GlobalTrueVSize();
+      if (mpi.Root())
+      {
+         cout << "\nAMR iteration " << it << endl;
+         cout << "Number of unknowns: " << global_dofs << endl;
+      }
+
+      // 14. Assemble the right-hand side and determine the list of true
+      //     (i.e. parallel conforming) essential boundary dofs.
+      Array<int> ess_tdof_list;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      b.Assemble();
+
+      // 15. Assemble the stiffness matrix. Note that MFEM doesn't care at this
+      //     point that the mesh is nonconforming and parallel.  The FE space is
+      //     considered 'cut' along hanging edges/faces, and also across
+      //     processor boundaries.
+      a.Assemble();
+
+      // 16. Create the parallel linear system: eliminate boundary conditions.
+      //     The system will be solved for true (unconstrained/unique) DOFs only.
+      OperatorPtr A;
+      Vector B, X;
+
+      sol.ProjectBdrCoefficientTangent(zeroCoef, ess_bdr);
+
+      const int copy_interior = 0;
+      a.FormLinearSystem(ess_tdof_list, sol, b, A, X, B, copy_interior);
+
+      // 17. Solve the linear system A X = B.
+      HypreAMS ams(*A.As<HypreParMatrix>(), &fespace);
+      ams.SetPrintLevel(0);
+
+      HyprePCG pcg(*A.As<HypreParMatrix>());
+      pcg.SetTol(1e-12);
+      pcg.SetMaxIter(1000);
+      pcg.SetPrintLevel(3);
+      pcg.SetPreconditioner(ams);
+      pcg.Mult(B, X);
+
+      // 18. Switch back to the host and extract the parallel grid function
+      //     corresponding to the finite element approximation X. This is the
+      //     local solution on each processor.
+      a.RecoverFEMSolution(X, b, sol);
+
+      // 19. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         DenseMatrix xyMat(2,3); xyMat = 0.0;
+         xyMat(0,0) = 1.0; xyMat(1,1) = 1.0;
+         MatrixConstantCoefficient xyMatCoef(xyMat);
+         Vector zVec(3); zVec = 0.0; zVec(2) = 1;
+         VectorConstantCoefficient zVecCoef(zVec);
+
+         VectorGridFunctionCoefficient solCoef(&sol);
+         MatrixVectorProductCoefficient xyCoef(xyMatCoef, solCoef);
+         InnerProductCoefficient zCoef(zVecCoef, solCoef);
+
+         H1_FECollection fec_h1(order, dim);
+         ND_FECollection fec_nd(order, dim);
+
+         ParFiniteElementSpace fes_h1(&pmesh, &fec_h1);
+         ParFiniteElementSpace fes_nd(&pmesh, &fec_nd);
+
+         ParGridFunction xyComp(&fes_nd);
+         ParGridFunction zComp(&fes_h1);
+
+         xyComp.ProjectCoefficient(xyCoef);
+         zComp.ProjectCoefficient(zCoef);
+
+         xy_sock << "parallel " << num_procs << " " << myid << "\n";
+         xy_sock << "solution\n" << pmesh << xyComp << flush;
+         if (it == 0)
+         {
+            xy_sock << "keys vvv "
+                    << "window_geometry 0 0 400 350 "
+                    << "window_title 'XY components'\n";
+         }
+
+         z_sock << "parallel " << num_procs << " " << myid << "\n"
+                << "solution\n" << pmesh << zComp << flush;
+         if (it == 0)
+         {
+            z_sock << "window_geometry 403 0 400 350 "
+                   << "window_title 'Z component'\n";
+         }
+      }
+
+      if (global_dofs > max_dofs)
+      {
+         if (mpi.Root())
+         {
+            cout << "Reached the maximum number of dofs. Stop." << endl;
+         }
+         break;
+      }
+
+      // 20. Call the refiner to modify the mesh. The refiner calls the error
+      //     estimator to obtain element errors, then it selects elements to be
+      //     refined and finally it modifies the mesh. The Stop() method can be
+      //     used to determine if a stopping criterion was met.
+      refiner.Apply(pmesh);
+      if (refiner.Stop())
+      {
+         if (mpi.Root())
+         {
+            cout << "Stopping criterion satisfied. Stop." << endl;
+         }
+         break;
+      }
+
+      // 21. Update the finite element space (recalculate the number of DOFs,
+      //     etc.) and create a grid function update matrix. Apply the matrix
+      //     to any GridFunctions over the space. In this case, the update
+      //     matrix is an interpolation matrix so the updated GridFunction will
+      //     still represent the same function as before refinement.
+      fespace.Update();
+      sol.Update();
+
+      // 22. Load balance the mesh, and update the space and solution. Currently
+      //     available only for nonconforming meshes.
+      if (pmesh.Nonconforming())
+      {
+         pmesh.Rebalance();
+
+         // Update the space and the GridFunction. This time the update matrix
+         // redistributes the GridFunction among the processors.
+         fespace.Update();
+         sol.Update();
+      }
+
+      // 23. Inform also the bilinear and linear forms that the space has
+      //     changed.
+      a.Update();
+      b.Update();
+   }
+   if (visualization)
+   {
+      xy_sock.close();
+      z_sock.close();
+   }
+
+   return 0;
+}
+
+void f_func(const Vector &x, Vector &f)
+{
+   double xc = 0.5 * (bb_min[0] + bb_max[0]);
+   double yc = 0.5 * (bb_min[1] + bb_max[1]);
+   double dx = bb_max[0] - bb_min[0];
+   double dy = bb_max[1] - bb_min[1];
+
+   f = 0.0;
+   if (fabs(x[0] - xc) < 0.2 * dx && fabs(x[1] - yc) < 0.2 * dy)
+   {
+      double a = pow(cos(2.5 * M_PI * (x[0] - xc) / dx) *
+                     cos(2.5 * M_PI * (x[1] - yc) / dy), 2);
+      f(0) = a * sin(2.5 * M_PI * (x[1] - yc) / dy);
+      f(1) = a * sin(5.0 * M_PI * (x[0] - xc) / dx);
+      f(2) = a * cos(5.0 * M_PI * (x[0] - xc) / dx);
+   }
+}

examples/makefile

@@ -22,10 +22,10 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
-  ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30_proposed
+  ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30_proposed ex32_proposed
 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_proposed ex31p_proposed
+  ex25p ex26p ex27p ex28p ex29p ex30p_proposed ex31p_proposed ex32p_proposed
 
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)

fem/estimators.cpp

@@ -31,25 +31,39 @@ void ZienkiewiczZhuEstimator::ComputeEstimates()
 }
 
 
-#ifdef MFEM_USE_MPI
-
 void L2ZienkiewiczZhuEstimator::ComputeEstimates()
 {
    flux_space->Update(false);
    smooth_flux_space->Update(false);
 
-   // TODO: move these parameters in the class, and add Set* methods.
-   const double solver_tol = 1e-12;
-   const int solver_max_it = 200;
-   total_error = L2ZZErrorEstimator(*integ, *solution, *smooth_flux_space,
-                                    *flux_space, error_estimates,
-                                    local_norm_p, solver_tol, solver_max_it);
+#ifdef MFEM_USE_MPI
+   // We need to force the compiler to find the parallel implementation
+   if (dist)
+   {
+      ParGridFunction * par_solution =
+         dynamic_cast<ParGridFunction*>(solution);
+      ParFiniteElementSpace * par_smooth_fes =
+         dynamic_cast<ParFiniteElementSpace*>(smooth_flux_space);
+      ParFiniteElementSpace * par_fes =
+         dynamic_cast<ParFiniteElementSpace*>(flux_space);
+      total_error = L2ZZErrorEstimator(*integ,
+                                       *par_solution,
+                                       *par_smooth_fes,
+                                       *par_fes, error_estimates,
+                                       local_norm_p, solver_tol,
+                                       solver_max_it);
+   }
+   else
+#endif
+   {
+      total_error = L2ZZErrorEstimator(*integ, *solution, *smooth_flux_space,
+                                       *flux_space, error_estimates,
+                                       local_norm_p, solver_tol, solver_max_it);
+   }
 
    current_sequence = solution->FESpace()->GetMesh()->GetSequence();
 }
 
-#endif // MFEM_USE_MPI
-
 KellyErrorEstimator::KellyErrorEstimator(BilinearFormIntegrator& di_,
                                          GridFunction& sol_,
                                          FiniteElementSpace& flux_fespace_,

fem/estimators.hpp

@@ -203,49 +203,55 @@ public:
 };
 
 
-#ifdef MFEM_USE_MPI
-
 /** @brief The L2ZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu
     error estimation procedure where the flux averaging is replaced by a global
     L2 projection (requiring a mass matrix solve).
 
     The required BilinearFormIntegrator must implement the methods
     ComputeElementFlux() and ComputeFluxEnergy().
-
-    Implemented for the parallel case only.
  */
 class L2ZienkiewiczZhuEstimator : public ErrorEstimator
 {
 protected:
    long current_sequence;
    int local_norm_p; ///< Local L_p norm to use, default is 1.
+   int solver_max_it;
+   double solver_tol;
    Vector error_estimates;
    double total_error;
+#ifdef  MFEM_USE_MPI
+   bool dist;
+#endif
 
    BilinearFormIntegrator *integ; ///< Not owned.
-   ParGridFunction *solution; ///< Not owned.
+   GridFunction *solution; ///< Not owned.
 
-   ParFiniteElementSpace *flux_space; /**< @brief Ownership based on the flag
+   FiniteElementSpace *flux_space; /**< @brief Ownership based on the flag
       own_flux_fes. Its Update() method is called automatically by this class
       when needed. */
-   ParFiniteElementSpace *smooth_flux_space; /**< @brief Ownership based on the
+   FiniteElementSpace *smooth_flux_space; /**< @brief Ownership based on the
       flag own_flux_fes. Its Update() method is called automatically by this
       class when needed.*/
    bool own_flux_fes; ///< Ownership flag for flux_space and smooth_flux_space.
 
    /// Initialize with the integrator, solution, and flux finite element spaces.
    void Init(BilinearFormIntegrator &integ,
-             ParGridFunction &sol,
-             ParFiniteElementSpace *flux_fes,
-             ParFiniteElementSpace *smooth_flux_fes)
+             GridFunction &sol,
+             FiniteElementSpace *flux_fes,
+             FiniteElementSpace *smooth_flux_fes)
    {
       current_sequence = -1;
       local_norm_p = 1;
+      solver_max_it = 200;
+      solver_tol = 1e-12;
       total_error = 0.0;
       this->integ = &integ;
       solution = &sol;
       flux_space = flux_fes;
       smooth_flux_space = smooth_flux_fes;
+#ifdef  MFEM_USE_MPI
+      dist = dynamic_cast<ParGridFunction*>(solution) != NULL;
+#endif
    }
 
    /// Check if the mesh of the solution was modified.
@@ -272,9 +278,9 @@ public:
                        FiniteElementSpace and will call its Update() method when
                        needed. */
    L2ZienkiewiczZhuEstimator(BilinearFormIntegrator &integ,
-                             ParGridFunction &sol,
-                             ParFiniteElementSpace *flux_fes,
-                             ParFiniteElementSpace *smooth_flux_fes)
+                             GridFunction &sol,
+                             FiniteElementSpace *flux_fes,
+                             FiniteElementSpace *smooth_flux_fes)
    { Init(integ, sol, flux_fes, smooth_flux_fes); own_flux_fes = true; }
 
    /** @brief Construct a new L2ZienkiewiczZhuEstimator object.
@@ -289,15 +295,21 @@ public:
                        of this FiniteElementSpace; will call its Update() method
                        when needed. */
    L2ZienkiewiczZhuEstimator(BilinearFormIntegrator &integ,
-                             ParGridFunction &sol,
-                             ParFiniteElementSpace &flux_fes,
-                             ParFiniteElementSpace &smooth_flux_fes)
+                             GridFunction &sol,
+                             FiniteElementSpace &flux_fes,
+                             FiniteElementSpace &smooth_flux_fes)
    { Init(integ, sol, &flux_fes, &smooth_flux_fes); own_flux_fes = false; }
 
    /** @brief Set the exponent, p, of the Lp norm used for computing the local
        element errors. Default value is 1. */
    void SetLocalErrorNormP(int p) { local_norm_p = p; }
 
+   /// Set maximum iteration count for global solve
+   void SetSolverMaxIt(int max_it) { solver_max_it = max_it; }
+
+   /// Set tolerance for global solve
+   void SetSolverTol(double tol) { solver_tol = tol; }
+
    /// Return the total error from the last error estimate.
    virtual double GetTotalError() const override { return total_error; }
 
@@ -319,8 +331,6 @@ public:
    }
 };
 
-#endif // MFEM_USE_MPI
-
 
 /** @brief The LpErrorEstimator class compares the solution to a known
     coefficient.

fem/gridfunc.cpp

@@ -19,6 +19,7 @@
 #include "pfespace.hpp"
 #endif
 
+#include "fem.hpp"
 #include <limits>
 #include <cstring>
 #include <string>
@@ -3934,6 +3935,94 @@ double ZZErrorEstimator(BilinearFormIntegrator &blfi,
 }
 
 
+double L2ZZErrorEstimator(BilinearFormIntegrator &flux_integrator,
+                          const GridFunction &x,
+                          FiniteElementSpace &smooth_flux_fes,
+                          FiniteElementSpace &flux_fes,
+                          Vector &errors,
+                          int norm_p, double solver_tol, int solver_max_it)
+{
+   // Compute fluxes in discontinuous space
+   GridFunction flux(&flux_fes);
+   flux = 0.0;
+
+   const FiniteElementSpace *xfes = x.FESpace();
+   Array<int> xdofs, fdofs;
+   Vector el_x, el_f;
+
+   for (int i = 0; i < xfes->GetNE(); i++)
+   {
+      xfes->GetElementVDofs(i, xdofs);
+      x.GetSubVector(xdofs, el_x);
+
+      ElementTransformation *Transf = xfes->GetElementTransformation(i);
+      flux_integrator.ComputeElementFlux(*xfes->GetFE(i), *Transf, el_x,
+                                         *flux_fes.GetFE(i), el_f, false);
+
+      flux_fes.GetElementVDofs(i, fdofs);
+      flux.AddElementVector(fdofs, el_f);
+   }
+
+   // Assemble the linear system for L2 projection into the "smooth" space
+   BilinearForm a(&smooth_flux_fes);
+   LinearForm b(&smooth_flux_fes);
+   VectorGridFunctionCoefficient f(&flux);
+
+   if (xfes->GetNE())
+   {
+      MFEM_VERIFY(smooth_flux_fes.GetFE(0) != NULL,
+                  "Could not obtain FE of smooth flux space.");
+
+      if (smooth_flux_fes.GetFE(0)->GetRangeType() == FiniteElement::SCALAR)
+      {
+         VectorMassIntegrator *vmass = new VectorMassIntegrator;
+         vmass->SetVDim(smooth_flux_fes.GetVDim());
+         a.AddDomainIntegrator(vmass);
+         b.AddDomainIntegrator(new VectorDomainLFIntegrator(f));
+      }
+      else
+      {
+         a.AddDomainIntegrator(new VectorFEMassIntegrator);
+         b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
+      }
+   }
+
+   b.Assemble();
+   a.Assemble();
+   a.Finalize();
+
+   // The destination of the projected discontinuous flux
+   GridFunction smooth_flux(&smooth_flux_fes);
+   smooth_flux = 0.0;
+
+   OperatorPtr A;
+   Vector B, X;
+   Array<int> ess_tdof_list(0);
+   a.FormLinearSystem(ess_tdof_list, smooth_flux, b, A, X, B);
+
+   // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
+   GSSmoother M((SparseMatrix&)(*A));
+   PCG(*A, M, B, X, 0, 200, 1e-12, 0.0);
+
+   // Extract the parallel grid function corresponding to the finite element
+   // approximation X. This is the local solution on each processor.
+   a.RecoverFEMSolution(X, b, smooth_flux);
+
+   // Proceed through the elements one by one, and find the Lp norm differences
+   // between the flux as computed per element and the flux projected onto the
+   // smooth_flux_fes space.
+   double total_error = 0.0;
+   errors.SetSize(xfes->GetNE());
+   for (int i = 0; i < xfes->GetNE(); i++)
+   {
+      errors(i) = ComputeElementLpDistance(norm_p, i, smooth_flux, flux);
+      total_error += pow(errors(i), norm_p);
+   }
+
+   return pow(total_error, 1.0/norm_p);
+}
+
+
 double ComputeElementLpDistance(double p, int i,
                                 GridFunction& gf1, GridFunction& gf2)
 {

fem/gridfunc.hpp

@@ -897,6 +897,19 @@ double ZZErrorEstimator(BilinearFormIntegrator &blfi,
                         int with_subdomains = 1,
                         bool with_coeff = false);
 
+/** Performs a global L2 projection (through a mass matrix solve) of flux
+    from supplied discontinuous space into supplied smooth (continuous, or at
+    least conforming) space, and computes the Lp norms of the differences
+    between them on each element. This is one approach to handling conforming
+    and non-conforming elements. Returns the total error estimate. */
+double L2ZZErrorEstimator(BilinearFormIntegrator &flux_integrator,
+                          const GridFunction &x,
+                          FiniteElementSpace &smooth_flux_fes,
+                          FiniteElementSpace &flux_fes,
+                          Vector &errors, int norm_p = 2,
+                          double solver_tol = 1e-12,
+                          int solver_max_it = 200);
+
 /// Compute the Lp distance between two grid functions on the given element.
 double ComputeElementLpDistance(double p, int i,
                                 GridFunction& gf1, GridFunction& gf2);

fem/pgridfunc.hpp

@@ -443,7 +443,8 @@ double L2ZZErrorEstimator(BilinearFormIntegrator &flux_integrator,
                           const ParGridFunction &x,
                           ParFiniteElementSpace &smooth_flux_fes,
                           ParFiniteElementSpace &flux_fes,
-                          Vector &errors, int norm_p = 2, double solver_tol = 1e-12,
+                          Vector &errors, int norm_p = 2,
+                          double solver_tol = 1e-12,
                           int solver_max_it = 200);
 
 }