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