MFEM  v3.3
Finite element discretization library
ex13p.cpp
Go to the documentation of this file.
1 // MFEM Example 13 - Parallel Version
2 //
3 // Compile with: make ex13p
4 //
5 // Sample runs: mpirun -np 4 ex13p -m ../data/star.mesh
6 // mpirun -np 4 ex13p -m ../data/square-disc.mesh -o 2
7 // mpirun -np 4 ex13p -m ../data/beam-tet.mesh
8 // mpirun -np 4 ex13p -m ../data/beam-hex.mesh
9 // mpirun -np 4 ex13p -m ../data/escher.mesh
10 // mpirun -np 4 ex13p -m ../data/fichera.mesh
11 // mpirun -np 4 ex13p -m ../data/fichera-q2.vtk
12 // mpirun -np 4 ex13p -m ../data/fichera-q3.mesh
13 // mpirun -np 4 ex13p -m ../data/square-disc-nurbs.mesh
14 // mpirun -np 4 ex13p -m ../data/beam-hex-nurbs.mesh
15 // mpirun -np 4 ex13p -m ../data/amr-quad.mesh -o 2
16 // mpirun -np 4 ex13p -m ../data/amr-hex.mesh
17 // mpirun -np 4 ex13p -m ../data/mobius-strip.mesh -n 8 -o 2
18 // mpirun -np 4 ex13p -m ../data/klein-bottle.mesh -n 10 -o 2
19 //
20 // Description: This example code solves the Maxwell (electromagnetic)
21 // eigenvalue problem curl curl E = lambda E with homogeneous
22 // Dirichlet boundary conditions E x n = 0.
23 //
24 // We compute a number of the lowest nonzero eigenmodes by
25 // discretizing the curl curl operator using a Nedelec FE space of
26 // the specified order in 2D or 3D.
27 //
28 // The example highlights the use of the AME subspace eigenvalue
29 // solver from HYPRE, which uses LOBPCG and AMS internally.
30 // Reusing a single GLVis visualization window for multiple
31 // eigenfunctions is also illustrated.
32 //
33 // We recommend viewing examples 3 and 11 before viewing this
34 // example.
35 
36 #include "mfem.hpp"
37 #include <fstream>
38 #include <iostream>
39 
40 using namespace std;
41 using namespace mfem;
42 
43 int main(int argc, char *argv[])
44 {
45  // 1. Initialize MPI.
46  int num_procs, myid;
47  MPI_Init(&argc, &argv);
48  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
49  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
50 
51  // 2. Parse command-line options.
52  const char *mesh_file = "../data/beam-tet.mesh";
53  int ser_ref_levels = 2;
54  int par_ref_levels = 1;
55  int order = 1;
56  int nev = 5;
57  bool visualization = 1;
58 
59  OptionsParser args(argc, argv);
60  args.AddOption(&mesh_file, "-m", "--mesh",
61  "Mesh file to use.");
62  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
63  "Number of times to refine the mesh uniformly in serial.");
64  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
65  "Number of times to refine the mesh uniformly in parallel.");
66  args.AddOption(&order, "-o", "--order",
67  "Finite element order (polynomial degree) or -1 for"
68  " isoparametric space.");
69  args.AddOption(&nev, "-n", "--num-eigs",
70  "Number of desired eigenmodes.");
71  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
72  "--no-visualization",
73  "Enable or disable GLVis visualization.");
74  args.Parse();
75  if (!args.Good())
76  {
77  if (myid == 0)
78  {
79  args.PrintUsage(cout);
80  }
81  MPI_Finalize();
82  return 1;
83  }
84  if (myid == 0)
85  {
86  args.PrintOptions(cout);
87  }
88 
89  // 3. Read the (serial) mesh from the given mesh file on all processors. We
90  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
91  // and volume meshes with the same code.
92  Mesh *mesh = new Mesh(mesh_file, 1, 1);
93  int dim = mesh->Dimension();
94 
95  // 4. Refine the serial mesh on all processors to increase the resolution. In
96  // this example we do 'ref_levels' of uniform refinement (2 by default, or
97  // specified on the command line with -rs).
98  for (int lev = 0; lev < ser_ref_levels; lev++)
99  {
100  mesh->UniformRefinement();
101  }
102 
103  // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
104  // this mesh further in parallel to increase the resolution (1 time by
105  // default, or specified on the command line with -rp). Once the parallel
106  // mesh is defined, the serial mesh can be deleted.
107  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
108  delete mesh;
109  for (int lev = 0; lev < par_ref_levels; lev++)
110  {
111  pmesh->UniformRefinement();
112  }
113 
114  // 6. Define a parallel finite element space on the parallel mesh. Here we
115  // use the Nedelec finite elements of the specified order.
116  FiniteElementCollection *fec = new ND_FECollection(order, dim);
117  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
118  HYPRE_Int size = fespace->GlobalTrueVSize();
119  if (myid == 0)
120  {
121  cout << "Number of unknowns: " << size << endl;
122  }
123 
124  // 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
125  // element space. The first corresponds to the curl curl, while the second
126  // is a simple mass matrix needed on the right hand side of the
127  // generalized eigenvalue problem below. The boundary conditions are
128  // implemented by marking all the boundary attributes from the mesh as
129  // essential. The corresponding degrees of freedom are eliminated with
130  // special values on the diagonal to shift the Dirichlet eigenvalues out
131  // of the computational range. After serial and parallel assembly we
132  // extract the corresponding parallel matrices A and M.
133  ConstantCoefficient one(1.0);
134  Array<int> ess_bdr;
135  if (pmesh->bdr_attributes.Size())
136  {
137  ess_bdr.SetSize(pmesh->bdr_attributes.Max());
138  ess_bdr = 1;
139  }
140 
141  ParBilinearForm *a = new ParBilinearForm(fespace);
143  if (pmesh->bdr_attributes.Size() == 0)
144  {
145  // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
146  // closed surface.
148  }
149  a->Assemble();
150  a->EliminateEssentialBCDiag(ess_bdr, 1.0);
151  a->Finalize();
152 
153  ParBilinearForm *m = new ParBilinearForm(fespace);
155  m->Assemble();
156  // shift the eigenvalue corresponding to eliminated dofs to a large value
157  m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
158  m->Finalize();
159 
162 
163  delete a;
164  delete m;
165 
166  // 8. Define and configure the AME eigensolver and the AMS preconditioner for
167  // A to be used within the solver. Set the matrices which define the
168  // generalized eigenproblem A x = lambda M x.
169  HypreAMS *ams = new HypreAMS(*A,fespace);
170  ams->SetPrintLevel(0);
171  ams->SetSingularProblem();
172 
173  HypreAME *ame = new HypreAME(MPI_COMM_WORLD);
174  ame->SetNumModes(nev);
175  ame->SetPreconditioner(*ams);
176  ame->SetMaxIter(100);
177  ame->SetTol(1e-8);
178  ame->SetPrintLevel(1);
179  ame->SetMassMatrix(*M);
180  ame->SetOperator(*A);
181 
182  // 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
183  // parallel grid function to represent each of the eigenmodes returned by
184  // the solver.
185  Array<double> eigenvalues;
186  ame->Solve();
187  ame->GetEigenvalues(eigenvalues);
188  ParGridFunction x(fespace);
189 
190  // 10. Save the refined mesh and the modes in parallel. This output can be
191  // viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
192  {
193  ostringstream mesh_name, mode_name;
194  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
195 
196  ofstream mesh_ofs(mesh_name.str().c_str());
197  mesh_ofs.precision(8);
198  pmesh->Print(mesh_ofs);
199 
200  for (int i=0; i<nev; i++)
201  {
202  // convert eigenvector from HypreParVector to ParGridFunction
203  x = ame->GetEigenvector(i);
204 
205  mode_name << "mode_" << setfill('0') << setw(2) << i << "."
206  << setfill('0') << setw(6) << myid;
207 
208  ofstream mode_ofs(mode_name.str().c_str());
209  mode_ofs.precision(8);
210  x.Save(mode_ofs);
211  mode_name.str("");
212  }
213  }
214 
215  // 11. Send the solution by socket to a GLVis server.
216  if (visualization)
217  {
218  char vishost[] = "localhost";
219  int visport = 19916;
220  socketstream mode_sock(vishost, visport);
221  mode_sock.precision(8);
222 
223  for (int i=0; i<nev; i++)
224  {
225  if ( myid == 0 )
226  {
227  cout << "Eigenmode " << i+1 << '/' << nev
228  << ", Lambda = " << eigenvalues[i] << endl;
229  }
230 
231  // convert eigenvector from HypreParVector to ParGridFunction
232  x = ame->GetEigenvector(i);
233 
234  mode_sock << "parallel " << num_procs << " " << myid << "\n"
235  << "solution\n" << *pmesh << x << flush
236  << "window_title 'Eigenmode " << i+1 << '/' << nev
237  << ", Lambda = " << eigenvalues[i] << "'" << endl;
238 
239  char c;
240  if (myid == 0)
241  {
242  cout << "press (q)uit or (c)ontinue --> " << flush;
243  cin >> c;
244  }
245  MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
246 
247  if (c != 'c')
248  {
249  break;
250  }
251  }
252  mode_sock.close();
253  }
254 
255  // 12. Free the used memory.
256  delete ame;
257  delete ams;
258  delete M;
259  delete A;
260 
261  delete fespace;
262  delete fec;
263  delete pmesh;
264 
265  MPI_Finalize();
266 
267  return 0;
268 }
int Size() const
Logical size of the array.
Definition: array.hpp:109
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:833
Subclass constant coefficient.
Definition: coefficient.hpp:57
Integrator for (curl u, curl v) for Nedelec elements.
void GetEigenvalues(Array< double > &eigenvalues)
Collect the converged eigenvalues.
Definition: hypre.cpp:3388
void SetPreconditioner(HypreSolver &precond)
Definition: hypre.cpp:3354
void SetTol(double tol)
Definition: hypre.cpp:3333
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:314
Abstract parallel finite element space.
Definition: pfespace.hpp:28
STL namespace.
void SetOperator(HypreParMatrix &A)
Definition: hypre.cpp:3360
int main(int argc, char *argv[])
Definition: ex13p.cpp:43
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
int dim
Definition: ex3.cpp:47
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6684
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:108
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void EliminateEssentialBCDiag(const Array< int > &bdr_attr_is_ess, double value)
Perform elimination and set the diagonal entry to the given value.
int Dimension() const
Definition: mesh.hpp:611
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:142
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2730
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Definition: optparser.hpp:74
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:349
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:3382
void SetNumModes(int num_eigs)
Definition: hypre.cpp:3325
void SetMaxIter(int max_iter)
Definition: hypre.cpp:3339
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void SetMassMatrix(HypreParMatrix &M)
Definition: hypre.cpp:3375
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
void SetPrintLevel(int logging)
Definition: hypre.cpp:3345
Class for parallel bilinear form.
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:297
HypreParVector & GetEigenvector(unsigned int i)
Extract a single eigenvector.
Definition: hypre.cpp:3424
Class for parallel grid function.
Definition: pgridfunc.hpp:31
Wrapper for hypre's ParCSR matrix class.
Definition: hypre.hpp:174
Class for parallel meshes.
Definition: pmesh.hpp:28
void SetSingularProblem()
Set this option when solving a curl-curl problem with zero mass term.
Definition: hypre.hpp:851
Integrator for (Q u, v) for VectorFiniteElements.
virtual void Print(std::ostream &out=std::cout) const
Definition: pmesh.cpp:3395
bool Good() const
Definition: optparser.hpp:120