MFEM  v4.5.0
Finite element discretization library
ex12p.cpp
Go to the documentation of this file.
1 // MFEM Example 12 - Parallel Version
2 //
3 // Compile with: make ex12p
4 //
5 // Sample runs:
6 // mpirun -np 4 ex12p -m ../data/beam-tri.mesh
7 // mpirun -np 4 ex12p -m ../data/beam-quad.mesh
8 // mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 462 -n 10 -o 2 -elast
9 // mpirun -np 4 ex12p -m ../data/beam-hex.mesh -s 3878
10 // mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 81
11 // mpirun -np 4 ex12p -m ../data/beam-tri.mesh -s 3877 -o 2 -sys
12 // mpirun -np 4 ex12p -m ../data/beam-quad.mesh -s 4544 -n 6 -o 3 -elast
13 // mpirun -np 4 ex12p -m ../data/beam-quad-nurbs.mesh
14 // mpirun -np 4 ex12p -m ../data/beam-hex-nurbs.mesh
15 //
16 // Description: This example code solves the linear elasticity eigenvalue
17 // problem for a multi-material cantilever beam.
18 //
19 // Specifically, we compute a number of the lowest eigenmodes by
20 // approximating the weak form of -div(sigma(u)) = lambda u where
21 // sigma(u)=lambda*div(u)*I+mu*(grad*u+u*grad) is the stress
22 // tensor corresponding to displacement field u, and lambda and mu
23 // are the material Lame constants. The boundary conditions are
24 // u=0 on the fixed part of the boundary with attribute 1, and
25 // sigma(u).n=f on the remainder. The geometry of the domain is
26 // assumed to be as follows:
27 //
28 // +----------+----------+
29 // boundary --->| material | material |
30 // attribute 1 | 1 | 2 |
31 // (fixed) +----------+----------+
32 //
33 // The example highlights the use of the LOBPCG eigenvalue solver
34 // together with the BoomerAMG preconditioner in HYPRE. Reusing a
35 // single GLVis visualization window for multiple eigenfunctions
36 // and optional saving with ADIOS2 (adios2.readthedocs.io) streams
37 // are also illustrated.
38 //
39 // We recommend viewing examples 2 and 11 before viewing this
40 // example.
41 
42 #include "mfem.hpp"
43 #include <fstream>
44 #include <iostream>
45 
46 using namespace std;
47 using namespace mfem;
48 
49 int main(int argc, char *argv[])
50 {
51  // 1. Initialize MPI and HYPRE.
52  Mpi::Init(argc, argv);
53  int num_procs = Mpi::WorldSize();
54  int myid = Mpi::WorldRank();
55  Hypre::Init();
56 
57  // 2. Parse command-line options.
58  const char *mesh_file = "../data/beam-tri.mesh";
59  int order = 1;
60  int nev = 5;
61  int seed = 66;
62  bool visualization = 1;
63  bool amg_elast = 0;
64  bool adios2 = false;
65 
66  OptionsParser args(argc, argv);
67  args.AddOption(&mesh_file, "-m", "--mesh",
68  "Mesh file to use.");
69  args.AddOption(&order, "-o", "--order",
70  "Finite element order (polynomial degree).");
71  args.AddOption(&nev, "-n", "--num-eigs",
72  "Number of desired eigenmodes.");
73  args.AddOption(&seed, "-s", "--seed",
74  "Random seed used to initialize LOBPCG.");
75  args.AddOption(&amg_elast, "-elast", "--amg-for-elasticity", "-sys",
76  "--amg-for-systems",
77  "Use the special AMG elasticity solver (GM/LN approaches), "
78  "or standard AMG for systems (unknown approach).");
79  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
80  "--no-visualization",
81  "Enable or disable GLVis visualization.");
82  args.AddOption(&adios2, "-adios2", "--adios2-streams", "-no-adios2",
83  "--no-adios2-streams",
84  "Save data using adios2 streams.");
85  args.Parse();
86  if (!args.Good())
87  {
88  if (myid == 0)
89  {
90  args.PrintUsage(cout);
91  }
92  return 1;
93  }
94  if (myid == 0)
95  {
96  args.PrintOptions(cout);
97  }
98 
99  // 3. Read the (serial) mesh from the given mesh file on all processors. We
100  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
101  // and volume meshes with the same code.
102  Mesh *mesh = new Mesh(mesh_file, 1, 1);
103  int dim = mesh->Dimension();
104 
105  if (mesh->attributes.Max() < 2)
106  {
107  if (myid == 0)
108  cerr << "\nInput mesh should have at least two materials!"
109  << " (See schematic in ex12p.cpp)\n"
110  << endl;
111  return 3;
112  }
113 
114  // 4. Select the order of the finite element discretization space. For NURBS
115  // meshes, we increase the order by degree elevation.
116  if (mesh->NURBSext)
117  {
118  mesh->DegreeElevate(order, order);
119  }
120 
121  // 5. Refine the serial mesh on all processors to increase the resolution. In
122  // this example we do 'ref_levels' of uniform refinement. We choose
123  // 'ref_levels' to be the largest number that gives a final mesh with no
124  // more than 1,000 elements.
125  {
126  int ref_levels =
127  (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
128  for (int l = 0; l < ref_levels; l++)
129  {
130  mesh->UniformRefinement();
131  }
132  }
133 
134  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
135  // this mesh further in parallel to increase the resolution. Once the
136  // parallel mesh is defined, the serial mesh can be deleted.
137  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
138  delete mesh;
139  {
140  int par_ref_levels = 1;
141  for (int l = 0; l < par_ref_levels; l++)
142  {
143  pmesh->UniformRefinement();
144  }
145  }
146 
147  // 7. Define a parallel finite element space on the parallel mesh. Here we
148  // use vector finite elements, i.e. dim copies of a scalar finite element
149  // space. We use the ordering by vector dimension (the last argument of
150  // the FiniteElementSpace constructor) which is expected in the systems
151  // version of BoomerAMG preconditioner. For NURBS meshes, we use the
152  // (degree elevated) NURBS space associated with the mesh nodes.
154  ParFiniteElementSpace *fespace;
155  const bool use_nodal_fespace = pmesh->NURBSext && !amg_elast;
156  if (use_nodal_fespace)
157  {
158  fec = NULL;
159  fespace = (ParFiniteElementSpace *)pmesh->GetNodes()->FESpace();
160  }
161  else
162  {
163  fec = new H1_FECollection(order, dim);
164  fespace = new ParFiniteElementSpace(pmesh, fec, dim, Ordering::byVDIM);
165  }
166  HYPRE_BigInt size = fespace->GlobalTrueVSize();
167  if (myid == 0)
168  {
169  cout << "Number of unknowns: " << size << endl
170  << "Assembling: " << flush;
171  }
172 
173  // 8. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
174  // element space corresponding to the linear elasticity integrator with
175  // piece-wise constants coefficient lambda and mu, a simple mass matrix
176  // needed on the right hand side of the generalized eigenvalue problem
177  // below. The boundary conditions are implemented by marking only boundary
178  // attribute 1 as essential. We use special values on the diagonal to
179  // shift the Dirichlet eigenvalues out of the computational range. After
180  // serial/parallel assembly we extract the corresponding parallel matrices
181  // A and M.
182  Vector lambda(pmesh->attributes.Max());
183  lambda = 1.0;
184  lambda(0) = lambda(1)*50;
185  PWConstCoefficient lambda_func(lambda);
186  Vector mu(pmesh->attributes.Max());
187  mu = 1.0;
188  mu(0) = mu(1)*50;
189  PWConstCoefficient mu_func(mu);
190 
191  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
192  ess_bdr = 0;
193  ess_bdr[0] = 1;
194 
195  ParBilinearForm *a = new ParBilinearForm(fespace);
196  a->AddDomainIntegrator(new ElasticityIntegrator(lambda_func, mu_func));
197  if (myid == 0)
198  {
199  cout << "matrix ... " << flush;
200  }
201  a->Assemble();
202  a->EliminateEssentialBCDiag(ess_bdr, 1.0);
203  a->Finalize();
204 
205  ParBilinearForm *m = new ParBilinearForm(fespace);
207  m->Assemble();
208  // shift the eigenvalue corresponding to eliminated dofs to a large value
209  m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
210  m->Finalize();
211  if (myid == 0)
212  {
213  cout << "done." << endl;
214  }
215 
216  HypreParMatrix *A = a->ParallelAssemble();
218 
219  delete a;
220  delete m;
221 
222  // 9. Define and configure the LOBPCG eigensolver and the BoomerAMG
223  // preconditioner for A to be used within the solver. Set the matrices
224  // which define the generalized eigenproblem A x = lambda M x.
225  HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
226  amg->SetPrintLevel(0);
227  if (amg_elast)
228  {
229  amg->SetElasticityOptions(fespace);
230  }
231  else
232  {
233  amg->SetSystemsOptions(dim);
234  }
235 
236  HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
237  lobpcg->SetNumModes(nev);
238  lobpcg->SetRandomSeed(seed);
239  lobpcg->SetPreconditioner(*amg);
240  lobpcg->SetMaxIter(100);
241  lobpcg->SetTol(1e-8);
242  lobpcg->SetPrecondUsageMode(1);
243  lobpcg->SetPrintLevel(1);
244  lobpcg->SetMassMatrix(*M);
245  lobpcg->SetOperator(*A);
246 
247  // 10. Compute the eigenmodes and extract the array of eigenvalues. Define a
248  // parallel grid function to represent each of the eigenmodes returned by
249  // the solver.
250  Array<double> eigenvalues;
251  lobpcg->Solve();
252  lobpcg->GetEigenvalues(eigenvalues);
253  ParGridFunction x(fespace);
254 
255  // 11. For non-NURBS meshes, make the mesh curved based on the finite element
256  // space. This means that we define the mesh elements through a fespace
257  // based transformation of the reference element. This allows us to save
258  // the displaced mesh as a curved mesh when using high-order finite
259  // element displacement field. We assume that the initial mesh (read from
260  // the file) is not higher order curved mesh compared to the chosen FE
261  // space.
262  if (!use_nodal_fespace)
263  {
264  pmesh->SetNodalFESpace(fespace);
265  }
266 
267  // 12. Save the refined mesh and the modes in parallel. This output can be
268  // viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
269  {
270  ostringstream mesh_name, mode_name;
271  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
272 
273  ofstream mesh_ofs(mesh_name.str().c_str());
274  mesh_ofs.precision(8);
275  pmesh->Print(mesh_ofs);
276 
277  for (int i=0; i<nev; i++)
278  {
279  // convert eigenvector from HypreParVector to ParGridFunction
280  x = lobpcg->GetEigenvector(i);
281 
282  mode_name << "mode_" << setfill('0') << setw(2) << i << "."
283  << setfill('0') << setw(6) << myid;
284 
285  ofstream mode_ofs(mode_name.str().c_str());
286  mode_ofs.precision(8);
287  x.Save(mode_ofs);
288  mode_name.str("");
289  }
290  }
291 
292  // 13. Optionally output a BP (binary pack) file using ADIOS2. This can be
293  // visualized with the ParaView VTX reader.
294 #ifdef MFEM_USE_ADIOS2
295  if (adios2)
296  {
297  std::string postfix(mesh_file);
298  postfix.erase(0, std::string("../data/").size() );
299  postfix += "_o" + std::to_string(order);
300 
301  adios2stream adios2output("ex12-p-" + postfix + ".bp",
302  adios2stream::openmode::out, MPI_COMM_WORLD);
303  pmesh->Print(adios2output);
304  for (int i=0; i<nev; i++)
305  {
306  x = lobpcg->GetEigenvector(i);
307  // x is a temporary that must be saved immediately
308  x.Save(adios2output, "mode_" + std::to_string(i));
309  }
310  }
311 #endif
312 
313  // 14. Send the above data by socket to a GLVis server. Use the "n" and "b"
314  // keys in GLVis to visualize the displacements.
315  if (visualization)
316  {
317  char vishost[] = "localhost";
318  int visport = 19916;
319  socketstream mode_sock(vishost, visport);
320 
321  for (int i=0; i<nev; i++)
322  {
323  if ( myid == 0 )
324  {
325  cout << "Eigenmode " << i+1 << '/' << nev
326  << ", Lambda = " << eigenvalues[i] << endl;
327  }
328 
329  // convert eigenvector from HypreParVector to ParGridFunction
330  x = lobpcg->GetEigenvector(i);
331 
332  mode_sock << "parallel " << num_procs << " " << myid << "\n"
333  << "solution\n" << *pmesh << x << flush
334  << "window_title 'Eigenmode " << i+1 << '/' << nev
335  << ", Lambda = " << eigenvalues[i] << "'" << endl;
336 
337  char c;
338  if (myid == 0)
339  {
340  cout << "press (q)uit or (c)ontinue --> " << flush;
341  cin >> c;
342  }
343  MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
344 
345  if (c != 'c')
346  {
347  break;
348  }
349  }
350  mode_sock.close();
351  }
352 
353  // 15. Free the used memory.
354  delete lobpcg;
355  delete amg;
356  delete M;
357  delete A;
358 
359  if (fec)
360  {
361  delete fespace;
362  delete fec;
363  }
364  delete pmesh;
365 
366  return 0;
367 }
void GetEigenvalues(Array< double > &eigenvalues) const
Collect the converged eigenvalues.
Definition: hypre.cpp:6229
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
int Dimension() const
Definition: mesh.hpp:1006
const HypreParVector & GetEigenvector(unsigned int i) const
Extract a single eigenvector.
Definition: hypre.cpp:6241
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
void SetPreconditioner(Solver &precond)
Definition: hypre.cpp:6163
void SetMassMatrix(Operator &M)
Definition: hypre.cpp:6219
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void SetPrintLevel(int logging)
Definition: hypre.cpp:6148
STL namespace.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void SetTol(double tol)
Definition: hypre.cpp:6126
int close()
Close the socketstream.
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1579
void SetMaxIter(int max_iter)
Definition: hypre.cpp:6142
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
constexpr char vishost[]
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9802
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1659
constexpr int visport
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 main(int argc, char *argv[])
Definition: ex12p.cpp:49
void SetElasticityOptions(ParFiniteElementSpace *fespace)
Definition: hypre.cpp:5083
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
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)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
HYPRE_Int HYPRE_BigInt
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:873
void SetRandomSeed(int s)
Definition: hypre.hpp:1993
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
A piecewise constant coefficient with the constants keyed off the element attribute numbers...
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
double mu
Definition: ex25.cpp:139
void DegreeElevate(int rel_degree, int degree=16)
Definition: mesh.cpp:5107
int dim
Definition: ex24.cpp:53
void SetOperator(Operator &A)
Definition: hypre.cpp:6172
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Class for parallel bilinear form.
Vector data type.
Definition: vector.hpp:60
void SetPrecondUsageMode(int pcg_mode)
Definition: hypre.cpp:6157
void SetSystemsOptions(int dim, bool order_bynodes=false)
Definition: hypre.cpp:4962
void SetNodalFESpace(FiniteElementSpace *nfes) override
Definition: pmesh.cpp:2076
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7894
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4770
Class for parallel grid function.
Definition: pgridfunc.hpp:32
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
void SetNumModes(int num_eigs)
Definition: hypre.hpp:1991
Class for parallel meshes.
Definition: pmesh.hpp:32
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:6286
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:268