MFEM  v3.1 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: mpirun -np 4 ex12p -m ../data/beam-tri.mesh
6 // mpirun -np 4 ex12p -m ../data/beam-quad.mesh
7 // mpirun -np 4 ex12p -m ../data/beam-tet.mesh -n 10 -o 2 -elast
8 // mpirun -np 4 ex12p -m ../data/beam-hex.mesh
9 // mpirun -np 4 ex12p -m ../data/beam-tri.mesh -o 2 -sys
10 // mpirun -np 4 ex12p -m ../data/beam-quad.mesh -n 6 -o 3 -elast
11 // mpirun -np 4 ex12p -m ../data/beam-quad-nurbs.mesh
12 // mpirun -np 4 ex12p -m ../data/beam-hex-nurbs.mesh
13 //
14 // Description: This example code solves the linear elasticity eigenvalue
15 // problem for a multi-material cantilever beam.
16 //
17 // Specifically, we compute a number of the lowest eigenmodes by
18 // approximating the weak form of -div(sigma(u)) = lambda u where
20 // tensor corresponding to displacement field u, and lambda and mu
21 // are the material Lame constants. The boundary conditions are
22 // u=0 on the fixed part of the boundary with attribute 1, and
23 // sigma(u).n=f on the remainder. The geometry of the domain is
24 // assumed to be as follows:
25 //
26 // +----------+----------+
27 // boundary --->| material | material |
28 // attribute 1 | 1 | 2 |
29 // (fixed) +----------+----------+
30 //
31 // The example highlights the use of the LOBPCG eigenvalue solver
32 // together with the BoomerAMG preconditioner in HYPRE. Reusing a
33 // single GLVis visualization window for multiple eigenfunctions
34 // is also illustrated.
35 //
36 // We recommend viewing examples 2 and 11 before viewing this
37 // example.
38
39 #include "mfem.hpp"
40 #include <fstream>
41 #include <iostream>
42
43 using namespace std;
44 using namespace mfem;
45
46 int main(int argc, char *argv[])
47 {
48  // 1. Initialize MPI.
49  int num_procs, myid;
50  MPI_Init(&argc, &argv);
51  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
52  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
53
54  // 2. Parse command-line options.
55  const char *mesh_file = "../data/beam-tri.mesh";
56  int order = 1;
57  int nev = 5;
58  bool visualization = 1;
59  bool amg_elast = 0;
60
61  OptionsParser args(argc, argv);
63  "Mesh file to use.");
65  "Finite element order (polynomial degree).");
67  "Number of desired eigenmodes.");
69  "--amg-for-systems",
70  "Use the special AMG elasticity solver (GM/LN approaches), "
71  "or standard AMG for systems (unknown approach).");
73  "--no-visualization",
74  "Enable or disable GLVis visualization.");
75  args.Parse();
76  if (!args.Good())
77  {
78  if (myid == 0)
79  {
80  args.PrintUsage(cout);
81  }
82  MPI_Finalize();
83  return 1;
84  }
85  if (myid == 0)
86  {
87  args.PrintOptions(cout);
88  }
89
90  // 3. Read the (serial) mesh from the given mesh file on all processors. We
91  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
92  // and volume meshes with the same code.
93  Mesh *mesh;
94  ifstream imesh(mesh_file);
95  if (!imesh)
96  {
97  if (myid == 0)
98  {
99  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
100  }
101  MPI_Finalize();
102  return 2;
103  }
104  mesh = new Mesh(imesh, 1, 1);
105  imesh.close();
106  int dim = mesh->Dimension();
107
108  if (mesh->attributes.Max() < 2)
109  {
110  if (myid == 0)
111  cerr << "\nInput mesh should have at least two materials!"
112  << " (See schematic in ex12p.cpp)\n"
113  << endl;
114  MPI_Finalize();
115  return 3;
116  }
117
118  // 4. Select the order of the finite element discretization space. For NURBS
119  // meshes, we increase the order by degree elevation.
120  if (mesh->NURBSext && order > mesh->NURBSext->GetOrder())
121  {
122  mesh->DegreeElevate(order - mesh->NURBSext->GetOrder());
123  }
124
125  // 5. Refine the serial mesh on all processors to increase the resolution. In
126  // this example we do 'ref_levels' of uniform refinement. We choose
127  // 'ref_levels' to be the largest number that gives a final mesh with no
128  // more than 1,000 elements.
129  {
130  int ref_levels =
131  (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
132  for (int l = 0; l < ref_levels; l++)
133  {
134  mesh->UniformRefinement();
135  }
136  }
137
138  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
139  // this mesh further in parallel to increase the resolution. Once the
140  // parallel mesh is defined, the serial mesh can be deleted.
141  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
142  delete mesh;
143  {
144  int par_ref_levels = 1;
145  for (int l = 0; l < par_ref_levels; l++)
146  {
147  pmesh->UniformRefinement();
148  }
149  }
150
151  // 7. Define a parallel finite element space on the parallel mesh. Here we
152  // use vector finite elements, i.e. dim copies of a scalar finite element
153  // space. We use the ordering by vector dimension (the last argument of
154  // the FiniteElementSpace constructor) which is expected in the systems
155  // version of BoomerAMG preconditioner. For NURBS meshes, we use the
156  // (degree elevated) NURBS space associated with the mesh nodes.
158  ParFiniteElementSpace *fespace;
159  const bool use_nodal_fespace = pmesh->NURBSext && !amg_elast;
160  if (use_nodal_fespace)
161  {
162  fec = NULL;
163  fespace = (ParFiniteElementSpace *)pmesh->GetNodes()->FESpace();
164  }
165  else
166  {
167  fec = new H1_FECollection(order, dim);
168  fespace = new ParFiniteElementSpace(pmesh, fec, dim, Ordering::byVDIM);
169  }
170  HYPRE_Int size = fespace->GlobalTrueVSize();
171  if (myid == 0)
172  {
173  cout << "Number of unknowns: " << size << endl
174  << "Assembling: " << flush;
175  }
176
177  // 8. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
178  // element space corresponding to the linear elasticity integrator with
179  // piece-wise constants coefficient lambda and mu, a simple mass matrix
180  // needed on the right hand side of the generalized eigenvalue problem
181  // below. The boundary conditions are implemented by marking only boundary
182  // attribute 1 as essential. We use special values on the diagonal to
183  // shift the Dirichlet eigenvalues out of the computational range. After
184  // serial/parallel assembly we extract the corresponding parallel matrices
185  // A and M.
186  Vector lambda(pmesh->attributes.Max());
187  lambda = 1.0;
188  lambda(0) = lambda(1)*50;
189  PWConstCoefficient lambda_func(lambda);
190  Vector mu(pmesh->attributes.Max());
191  mu = 1.0;
192  mu(0) = mu(1)*50;
193  PWConstCoefficient mu_func(mu);
194
195  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
196  ess_bdr = 0;
197  ess_bdr[0] = 1;
198
199  ParBilinearForm *a = new ParBilinearForm(fespace);
201  if (myid == 0)
202  {
203  cout << "matrix ... " << flush;
204  }
205  a->Assemble();
206  a->EliminateEssentialBCDiag(ess_bdr, 1.0);
207  a->Finalize();
208
209  ParBilinearForm *m = new ParBilinearForm(fespace);
211  m->Assemble();
212  // shift the eigenvalue corresponding to eliminated dofs to a large value
213  m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
214  m->Finalize();
215  if (myid == 0)
216  {
217  cout << "done." << endl;
218  }
219
222
223  delete a;
224  delete m;
225
226  // 9. Define and configure the LOBPCG eigensolver and the BoomerAMG
227  // preconditioner for A to be used within the solver. Set the matrices
228  // which define the generalized eigenproblem A x = lambda M x.
229  HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
230  amg->SetPrintLevel(0);
231  if (amg_elast)
232  {
233  amg->SetElasticityOptions(fespace);
234  }
235  else
236  {
237  amg->SetSystemsOptions(dim);
238  }
239
240  HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
241  lobpcg->SetNumModes(nev);
242  lobpcg->SetPreconditioner(*amg);
243  lobpcg->SetMaxIter(100);
244  lobpcg->SetTol(1e-8);
245  lobpcg->SetPrecondUsageMode(1);
246  lobpcg->SetPrintLevel(1);
247  lobpcg->SetMassMatrix(*M);
248  lobpcg->SetOperator(*A);
249
250  // 10. Compute the eigenmodes and extract the array of eigenvalues. Define a
251  // parallel grid function to represent each of the eigenmodes returned by
252  // the solver.
253  Array<double> eigenvalues;
254  lobpcg->Solve();
255  lobpcg->GetEigenvalues(eigenvalues);
256  ParGridFunction x(fespace);
257
258  // 11. For non-NURBS meshes, make the mesh curved based on the finite element
259  // space. This means that we define the mesh elements through a fespace
260  // based transformation of the reference element. This allows us to save
261  // the displaced mesh as a curved mesh when using high-order finite
262  // element displacement field. We assume that the initial mesh (read from
263  // the file) is not higher order curved mesh compared to the chosen FE
264  // space.
265  if (!use_nodal_fespace)
266  {
267  pmesh->SetNodalFESpace(fespace);
268  }
269
270  // 12. Save the refined mesh and the modes in parallel. This output can be
271  // viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
272  {
273  ostringstream mesh_name, mode_name;
274  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
275
276  ofstream mesh_ofs(mesh_name.str().c_str());
277  mesh_ofs.precision(8);
278  pmesh->Print(mesh_ofs);
279
280  for (int i=0; i<nev; i++)
281  {
282  // convert eigenvector from HypreParVector to ParGridFunction
283  x = lobpcg->GetEigenvector(i);
284
285  mode_name << "mode_" << setfill('0') << setw(2) << i << "."
286  << setfill('0') << setw(6) << myid;
287
288  ofstream mode_ofs(mode_name.str().c_str());
289  mode_ofs.precision(8);
290  x.Save(mode_ofs);
291  mode_name.str("");
292  }
293  }
294
295  // 13. Send the above data by socket to a GLVis server. Use the "n" and "b"
296  // keys in GLVis to visualize the displacements.
297  if (visualization)
298  {
299  char vishost[] = "localhost";
300  int visport = 19916;
301  socketstream mode_sock(vishost, visport);
302
303  for (int i=0; i<nev; i++)
304  {
305  if ( myid == 0 )
306  {
307  cout << "Eigenmode " << i+1 << '/' << nev
308  << ", Lambda = " << eigenvalues[i] << endl;
309  }
310
311  // convert eigenvector from HypreParVector to ParGridFunction
312  x = lobpcg->GetEigenvector(i);
313
314  mode_sock << "parallel " << num_procs << " " << myid << "\n"
315  << "solution\n" << *pmesh << x << flush
316  << "window_title 'Eigenmode " << i+1 << '/' << nev
317  << ", Lambda = " << eigenvalues[i] << "'" << endl;
318
319  char c;
320  if (myid == 0)
321  {
322  cout << "press (q)uit or (c)ontinue --> " << flush;
323  cin >> c;
324  }
325  MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
326
327  if (c != 'c')
328  {
329  break;
330  }
331  }
332  mode_sock.close();
333  }
334
335  // 14. Free the used memory.
336  delete lobpcg;
337  delete amg;
338  delete M;
339  delete A;
340
341  if (fec)
342  {
343  delete fespace;
344  delete fec;
345  }
346  delete pmesh;
347
348  MPI_Finalize();
349
350  return 0;
351 }
HypreParVector & GetEigenvector(unsigned int i)
Extract a single eigenvector.
Definition: hypre.cpp:3102
void DegreeElevate(int t)
Definition: mesh.cpp:3497
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:454
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:306
void SetPreconditioner(Solver &precond)
Definition: hypre.cpp:3025
void SetMassMatrix(Operator &M)
Definition: hypre.cpp:3080
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void SetPrintLevel(int logging)
Definition: hypre.cpp:3010
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void SetTol(double tol)
Definition: hypre.cpp:2998
void SetSystemsOptions(int dim)
Definition: hypre.cpp:2326
The BoomerAMG solver in hypre.
Definition: hypre.hpp:698
int dim
Definition: ex3.cpp:48
void SetMaxIter(int max_iter)
Definition: hypre.cpp:3004
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7316
void SetPrintLevel(int print_level)
Definition: hypre.hpp:739
T Max() const
Definition: array.cpp:90
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void SetElasticityOptions(ParFiniteElementSpace *fespace)
Definition: hypre.cpp:2408
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3719
int Dimension() const
Definition: mesh.hpp:475
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void EliminateEssentialBCDiag(Array< int > &bdr_attr_is_ess, double value)
Perform elimination and set the diagonal entry to the given value.
Array< int > bdr_attributes
Definition: mesh.hpp:140
int main(int argc, char *argv[])
Definition: ex1.cpp:45
void GetEigenvalues(Array< double > &eigenvalues)
Collect the converged eigenvalues.
Definition: hypre.cpp:3090
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
NURBSExtension * NURBSext
Definition: mesh.hpp:142
class for piecewise constant coefficient
Definition: coefficient.hpp:72
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void SetOperator(Operator &A)
Definition: hypre.cpp:3034
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
Class for parallel bilinear form.
Vector data type.
Definition: vector.hpp:33
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5844
void SetPrecondUsageMode(int pcg_mode)
Definition: hypre.cpp:3019
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:54
Class for parallel grid function.
Definition: pgridfunc.hpp:31
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:143
void SetNumModes(int num_eigs)
Definition: hypre.hpp:937
Class for parallel meshes.
Definition: pmesh.hpp:28
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:3108
Array< int > attributes
Definition: mesh.hpp:139
virtual void Print(std::ostream &out=std::cout) const
Definition: pmesh.cpp:2850
bool Good() const
Definition: optparser.hpp:120