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