MFEM  v4.0 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);
61  "Mesh file to use.");
63  "Number of times to refine the mesh uniformly in serial.");
65  "Number of times to refine the mesh uniformly in parallel.");
67  "Finite element order (polynomial degree) or -1 for"
68  " isoparametric space.");
70  "Number of desired eigenmodes.");
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  pmesh->ReorientTetMesh();
114
115  // 6. Define a parallel finite element space on the parallel mesh. Here we
116  // use the Nedelec finite elements of the specified order.
117  FiniteElementCollection *fec = new ND_FECollection(order, dim);
118  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
119  HYPRE_Int size = fespace->GlobalTrueVSize();
120  if (myid == 0)
121  {
122  cout << "Number of unknowns: " << size << endl;
123  }
124
125  // 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
126  // element space. The first corresponds to the curl curl, while the second
127  // is a simple mass matrix needed on the right hand side of the
128  // generalized eigenvalue problem below. The boundary conditions are
129  // implemented by marking all the boundary attributes from the mesh as
130  // essential. The corresponding degrees of freedom are eliminated with
131  // special values on the diagonal to shift the Dirichlet eigenvalues out
132  // of the computational range. After serial and parallel assembly we
133  // extract the corresponding parallel matrices A and M.
134  ConstantCoefficient one(1.0);
135  Array<int> ess_bdr;
136  if (pmesh->bdr_attributes.Size())
137  {
138  ess_bdr.SetSize(pmesh->bdr_attributes.Max());
139  ess_bdr = 1;
140  }
141
142  ParBilinearForm *a = new ParBilinearForm(fespace);
144  if (pmesh->bdr_attributes.Size() == 0)
145  {
146  // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
147  // closed surface.
149  }
150  a->Assemble();
151  a->EliminateEssentialBCDiag(ess_bdr, 1.0);
152  a->Finalize();
153
154  ParBilinearForm *m = new ParBilinearForm(fespace);
156  m->Assemble();
157  // shift the eigenvalue corresponding to eliminated dofs to a large value
158  m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
159  m->Finalize();
160
163
164  delete a;
165  delete m;
166
167  // 8. Define and configure the AME eigensolver and the AMS preconditioner for
168  // A to be used within the solver. Set the matrices which define the
169  // generalized eigenproblem A x = lambda M x.
170  HypreAMS *ams = new HypreAMS(*A,fespace);
171  ams->SetPrintLevel(0);
172  ams->SetSingularProblem();
173
174  HypreAME *ame = new HypreAME(MPI_COMM_WORLD);
175  ame->SetNumModes(nev);
176  ame->SetPreconditioner(*ams);
177  ame->SetMaxIter(100);
178  ame->SetTol(1e-8);
179  ame->SetPrintLevel(1);
180  ame->SetMassMatrix(*M);
181  ame->SetOperator(*A);
182
183  // 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
184  // parallel grid function to represent each of the eigenmodes returned by
185  // the solver.
186  Array<double> eigenvalues;
187  ame->Solve();
188  ame->GetEigenvalues(eigenvalues);
189  ParGridFunction x(fespace);
190
191  // 10. Save the refined mesh and the modes in parallel. This output can be
192  // viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
193  {
194  ostringstream mesh_name, mode_name;
195  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
196
197  ofstream mesh_ofs(mesh_name.str().c_str());
198  mesh_ofs.precision(8);
199  pmesh->Print(mesh_ofs);
200
201  for (int i=0; i<nev; i++)
202  {
203  // convert eigenvector from HypreParVector to ParGridFunction
204  x = ame->GetEigenvector(i);
205
206  mode_name << "mode_" << setfill('0') << setw(2) << i << "."
207  << setfill('0') << setw(6) << myid;
208
209  ofstream mode_ofs(mode_name.str().c_str());
210  mode_ofs.precision(8);
211  x.Save(mode_ofs);
212  mode_name.str("");
213  }
214  }
215
216  // 11. Send the solution by socket to a GLVis server.
217  if (visualization)
218  {
219  char vishost[] = "localhost";
220  int visport = 19916;
221  socketstream mode_sock(vishost, visport);
222  mode_sock.precision(8);
223
224  for (int i=0; i<nev; i++)
225  {
226  if ( myid == 0 )
227  {
228  cout << "Eigenmode " << i+1 << '/' << nev
229  << ", Lambda = " << eigenvalues[i] << endl;
230  }
231
232  // convert eigenvector from HypreParVector to ParGridFunction
233  x = ame->GetEigenvector(i);
234
235  mode_sock << "parallel " << num_procs << " " << myid << "\n"
236  << "solution\n" << *pmesh << x << flush
237  << "window_title 'Eigenmode " << i+1 << '/' << nev
238  << ", Lambda = " << eigenvalues[i] << "'" << endl;
239
240  char c;
241  if (myid == 0)
242  {
243  cout << "press (q)uit or (c)ontinue --> " << flush;
244  cin >> c;
245  }
246  MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
247
248  if (c != 'c')
249  {
250  break;
251  }
252  }
253  mode_sock.close();
254  }
255
256  // 12. Free the used memory.
257  delete ame;
258  delete ams;
259  delete M;
260  delete A;
261
262  delete fespace;
263  delete fec;
264  delete pmesh;
265
266  MPI_Finalize();
267
268  return 0;
269 }
int Size() const
Logical size of the array.
Definition: array.hpp:118
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:936
Subclass constant coefficient.
Definition: coefficient.hpp:67
virtual void ReorientTetMesh()
See the remarks for the serial version in mesh.hpp.
Definition: pmesh.cpp:2402
Integrator for (curl u, curl v) for Nedelec elements.
void GetEigenvalues(Array< double > &eigenvalues)
Collect the converged eigenvalues.
Definition: hypre.cpp:3679
void SetPreconditioner(HypreSolver &precond)
Definition: hypre.cpp:3645
void SetTol(double tol)
Definition: hypre.cpp:3614
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:488
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void SetOperator(HypreParMatrix &A)
Definition: hypre.cpp:3651
int main(int argc, char *argv[])
Definition: ex1.cpp:58
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
int dim
Definition: ex3.cpp:48
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7591
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
HYPRE_Int GlobalTrueVSize() const
Definition: pfespace.hpp:251
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.
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:3863
int Dimension() const
Definition: mesh.hpp:713
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:179
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2948
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:76
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:618
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:3673
void SetNumModes(int num_eigs)
Definition: hypre.cpp:3606
void SetMaxIter(int max_iter)
Definition: hypre.cpp:3630
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void SetMassMatrix(HypreParMatrix &M)
Definition: hypre.cpp:3666
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
void SetPrintLevel(int logging)
Definition: hypre.cpp:3636
Class for parallel bilinear form.
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:240
HypreParVector & GetEigenvector(unsigned int i)
Extract a single eigenvector.
Definition: hypre.cpp:3715
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:187
Class for parallel meshes.
Definition: pmesh.hpp:32
void SetSingularProblem()
Set this option when solving a curl-curl problem with zero mass term.
Definition: hypre.hpp:954
Integrator for (Q u, v) for VectorFiniteElements.
bool Good() const
Definition: optparser.hpp:122