Merge branch 'master' into ex13p-amr-dev

Conflicts:
	examples/makefile

examples/ex13p_amr.cpp

@@ -0,0 +1,375 @@
+//                       MFEM Example 13 - Parallel, AMR Version
+//
+// Compile with: make ex13p_amr
+//
+// Sample runs:  mpirun -np 4 ex13p_amr -m ../data/beam-tet.mesh
+//               mpirun -np 4 ex13p_amr -m ../data/beam-hex.mesh
+//               mpirun -np 4 ex13p_amr -m ../data/escher.mesh
+//               mpirun -np 4 ex13p_amr -m ../data/fichera.mesh -rs 2
+//               mpirun -np 4 ex13p_amr -m ../data/fichera-q2.vtk
+//               mpirun -np 4 ex13p_amr -m ../data/fichera-q3.mesh -rs 2
+//               mpirun -np 4 ex13p_amr -m ../data/beam-hex-nurbs.mesh
+//               mpirun -np 4 ex13p_amr -m ../data/amr-hex.mesh
+//
+// Description:  This example code solves the Maxwell (electromagnetic)
+//               eigenvalue problem curl curl E = lambda E with homogeneous
+//               Dirichlet boundary conditions E x n = 0.
+//
+//               We compute a number of the lowest nonzero eigenmodes
+//               by discretizing the curl curl operator using a
+//               Nedelec FE space of the specified order in 3D.
+//
+//               The example highlights the use of the AME subspace
+//               eigenvalue solver from HYPRE, which uses LOBPCG and
+//               AMS internally, in the context of adaptive mesh
+//               refinement (AMR).  Reusing a single GLVis
+//               visualization window for multiple eigenfunctions is
+//               also illustrated.
+//
+//               We recommend viewing example 13 before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int ser_ref_levels = 1;
+   int par_ref_levels = 0;
+   int max_amr_levels = 5;
+   int order = 1;
+   int nev = 5;
+   bool visualization = 1;
+   double tol = 1.0e-3;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   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(&max_amr_levels, "-ra", "--refine-adaptively",
+                  "Number of times to adaptively refine the mesh.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&nev, "-n", "--num-eigs",
+                  "Number of desired eigenmodes.");
+   args.AddOption(&tol, "-t", "--tolerance",
+                  "Desired eigenvalue tolerance.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Read the (serial) mesh from the given mesh file on all processors. We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
+   //    specified on the command line with -rs).
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+   if (mesh->NURBSext)
+   {
+      mesh->SetCurvature(2);
+   }
+   mesh->EnsureNCMesh();
+
+   // 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 the Nedelec finite elements of the specified order.
+   FiniteElementCollection *fec = new ND_FECollection(order, dim);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->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.
+   ConstantCoefficient one(1.0);
+   Array<int> ess_bdr;
+   if (pmesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+   }
+
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   a->AddDomainIntegrator(new CurlCurlIntegrator(one));
+   if (pmesh->bdr_attributes.Size() == 0)
+   {
+      // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
+      // closed surface.
+      a->AddDomainIntegrator(new VectorFEMassIntegrator(one));
+   }
+   a->Assemble();
+
+   ParBilinearForm *m = new ParBilinearForm(fespace);
+   m->AddDomainIntegrator(new VectorFEMassIntegrator(one));
+
+   ParGridFunction x(fespace);
+
+   // 8. Setup the AMR loop which will stop after a set number of
+   //    iterations or when the computed eigenvalues change by less than
+   //    some toleratonce.
+   HypreAMS *ams = NULL;
+   HypreAME *ame = NULL;
+
+   Array<double> eigenvalues;
+   Array<double> prev_eigenvalues(nev);
+   Vector eigenvalue_diff(nev);
+   prev_eigenvalues = 0.0;
+
+   double err = 2.0*tol;
+   int it = 0;
+
+   while ( err > tol && it < max_amr_levels )
+   {
+      // 8a. If this is not the first time through the loop, estimate
+      //     the error in the eigenvectors and use this error to mark
+      //     elements for refinement.
+      if ( it > 0 )
+      {
+         Vector errors(pmesh->GetNE());
+         Vector errors_i(pmesh->GetNE());
+
+         // Space for the discontinuous (original) flux
+         CurlCurlIntegrator flux_integrator(one);
+         RT_FECollection flux_fec(order-1, pmesh->SpaceDimension());
+         ParFiniteElementSpace flux_fes(pmesh, &flux_fec);
+
+         // Space for the smoothed (conforming) flux
+         ND_FECollection smooth_flux_fec(order, pmesh->Dimension());
+         ParFiniteElementSpace smooth_flux_fes(pmesh, &smooth_flux_fec);
+
+         double norm_p = 1;
+         errors = 0.0;
+         for (int i=0; i<nev; i++)
+         {
+            // convert eigenvector from HypreParVector to ParGridFunction
+            x = ame->GetEigenvector(i);
+
+            L2ZZErrorEstimator(flux_integrator, x,
+                               smooth_flux_fes, flux_fes, errors_i, norm_p);
+
+            for (int j=0; j<errors.Size(); j++)
+            {
+               errors[j] += pow(errors_i[j], norm_p);
+            }
+         }
+         for (int j=0; j<errors.Size(); j++)
+         {
+            errors[j] = pow(errors[j], 1.0/norm_p);
+         }
+
+         double local_max_err = errors.Max();
+         double global_max_err;
+         MPI_Allreduce(&local_max_err, &global_max_err, 1,
+                       MPI_DOUBLE, MPI_MAX, pmesh->GetComm());
+
+         // Refine the elements whose error is larger than a fraction of the
+         // maximum element error.
+         const double frac = 0.5;
+         double threshold = frac * global_max_err;
+         if ( myid == 0 )
+         {
+            cout << "AMR iteration " << it+1 << ", refining from "
+                 << pmesh->GetNE() << " elements" << flush;
+         }
+         pmesh->RefineByError(errors, threshold);
+         if ( myid == 0 )
+         {
+            cout << " to " << pmesh->GetNE() << "." << endl;
+         }
+
+         fespace->Update();
+
+         a->Update();
+         m->Update();
+         x.Update();
+      }
+      a->Assemble();
+      a->EliminateEssentialBCDiag(ess_bdr, 1.0);
+      a->Finalize();
+
+      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();
+
+      // 8b. 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.
+      delete ams;
+      ams = new HypreAMS(*A,fespace);
+      ams->SetPrintLevel(0);
+      ams->SetSingularProblem();
+
+      delete ame;
+      ame = new HypreAME(MPI_COMM_WORLD);
+      ame->SetNumModes(nev);
+      ame->SetPreconditioner(*ams);
+      ame->SetMaxIter(100);
+      ame->SetTol(1e-2 * tol);
+      ame->SetPrintLevel(1);
+      ame->SetMassMatrix(*M);
+      ame->SetOperator(*A);
+
+      // 8c. Compute the eigenmodes and extract the array of
+      //     eigenvalues.  Define a parallel grid function to
+      //     represent each of the eigenmodes returned by the solver.
+      ame->Solve();
+      ame->GetEigenvalues(eigenvalues);
+
+      delete M;
+      delete A;
+
+      double avg = 0.0;
+      for (int i=0; i<nev; i++)
+      {
+         eigenvalue_diff[i] = eigenvalues[i] - prev_eigenvalues[i];
+         avg += 0.25 * pow(eigenvalues[i] + prev_eigenvalues[i], 2.0);
+      }
+      avg = sqrt(avg);
+      err = eigenvalue_diff.Norml2();
+      if ( avg > tol ) { err /= avg * nev; }
+
+      it++;
+   }
+
+   delete a;
+   delete m;
+
+
+   // 9. 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("");
+      }
+   }
+
+   // 10. 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);
+
+      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);
+
+         mode_sock << "parallel " << num_procs << " " << myid << "\n"
+                   << "solution\n" << *pmesh << x << flush
+                   << "window_title 'Eigenmode " << i+1 << '/' << nev
+                   << ", Lambda = " << eigenvalues[i] << "'" << 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();
+   }
+
+   // 11. Free the used memory.
+   delete ame;
+   delete ams;
+
+   delete fespace;
+   delete fec;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}

examples/makefile

@@ -25,7 +25,7 @@ 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 \
   ex31 ex33
 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 \
+  ex12p ex13p ex13p_amr ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
   ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p
 SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex22 ex24 ex25 ex26
 PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex22p \