MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
nurbs_ex11p.cpp
Go to the documentation of this file.
1 // MFEM Example 11 - Parallel NURBS Version
2 //
3 // Compile with: make nurbs_ex11p
4 //
5 // Sample runs: mpirun -np 4 nurbs_ex11p -m ../../data/square-disc.mesh
6 // mpirun -np 4 nurbs_ex11p -m ../../data/star.mesh
7 // mpirun -np 4 nurbs_ex11p -m ../../data/escher.mesh
8 // mpirun -np 4 nurbs_ex11p -m ../../data/fichera.mesh
9 // mpirun -np 4 nurbs_ex11p -m ../../data/square-disc-p2.vtk -o 2
10 // mpirun -np 4 nurbs_ex11p -m ../../data/square-disc-p3.mesh -o 3
11 // mpirun -np 4 nurbs_ex11p -m ../../data/square-disc-nurbs.mesh -o -1
12 // mpirun -np 4 nurbs_ex11p -m ../../data/disc-nurbs.mesh -o -1 -n 20
13 // mpirun -np 4 nurbs_ex11p -m ../../data/pipe-nurbs.mesh -o -1
14 // mpirun -np 4 nurbs_ex11p -m ../../data/ball-nurbs.mesh -o 2
15 // mpirun -np 4 nurbs_ex11p -m ../../data/star-surf.mesh
16 // mpirun -np 4 nurbs_ex11p -m ../../data/square-disc-surf.mesh
17 // mpirun -np 4 nurbs_ex11p -m ../../data/inline-segment.mesh
18 // mpirun -np 4 nurbs_ex11p -m ../../data/amr-quad.mesh
19 // mpirun -np 4 nurbs_ex11p -m ../../data/amr-hex.mesh
20 // mpirun -np 4 nurbs_ex11p -m ../../data/mobius-strip.mesh -n 8
21 // mpirun -np 4 nurbs_ex11p -m ../../data/klein-bottle.mesh -n 10
22 //
23 // Description: This example code demonstrates the use of MFEM to solve the
24 // eigenvalue problem -Delta u = lambda u with homogeneous
25 // Dirichlet boundary conditions.
26 //
27 // We compute a number of the lowest eigenmodes by discretizing
28 // the Laplacian and Mass operators using a FE space of the
29 // specified order, or an isoparametric/isogeometric space if
30 // order < 1 (quadratic for quadratic curvilinear mesh, NURBS for
31 // NURBS mesh, etc.)
32 //
33 // The example highlights the use of the LOBPCG eigenvalue solver
34 // together with the BoomerAMG preconditioner in HYPRE, as well as
35 // optionally the SuperLU or STRUMPACK parallel direct solvers.
36 // Reusing a single GLVis visualization window for multiple
37 // eigenfunctions is also illustrated.
38 //
39 // We recommend viewing Example 1 before viewing this 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/star.mesh";
58  int ser_ref_levels = 2;
59  int par_ref_levels = 1;
60  Array<int> order(1);
61  order[0] = 0;
62  int nev = 5;
63  int seed = 75;
64  bool slu_solver = false;
65  bool sp_solver = false;
66  bool visualization = 1;
67 
68  OptionsParser args(argc, argv);
69  args.AddOption(&mesh_file, "-m", "--mesh",
70  "Mesh file to use.");
71  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
72  "Number of times to refine the mesh uniformly in serial.");
73  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
74  "Number of times to refine the mesh uniformly in parallel.");
75  args.AddOption(&order, "-o", "--order",
76  "Finite element order (polynomial degree) or -1 for"
77  " isoparametric space.");
78  args.AddOption(&nev, "-n", "--num-eigs",
79  "Number of desired eigenmodes.");
80  args.AddOption(&seed, "-s", "--seed",
81  "Random seed used to initialize LOBPCG.");
82 #ifdef MFEM_USE_SUPERLU
83  args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
84  "--no-superlu", "Use the SuperLU Solver.");
85 #endif
86 #ifdef MFEM_USE_STRUMPACK
87  args.AddOption(&sp_solver, "-sp", "--strumpack", "-no-sp",
88  "--no-strumpack", "Use the STRUMPACK Solver.");
89 #endif
90  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
91  "--no-visualization",
92  "Enable or disable GLVis visualization.");
93  args.Parse();
94  if (slu_solver && sp_solver)
95  {
96  if (myid == 0)
97  cout << "WARNING: Both SuperLU and STRUMPACK have been selected,"
98  << " please choose either one." << endl
99  << " Defaulting to SuperLU." << endl;
100  sp_solver = false;
101  }
102  // The command line options are also passed to the STRUMPACK
103  // solver. So do not exit if some options are not recognized.
104  if (!sp_solver)
105  {
106  if (!args.Good())
107  {
108  if (myid == 0)
109  {
110  args.PrintUsage(cout);
111  }
112  MPI_Finalize();
113  return 1;
114  }
115  }
116  if (myid == 0)
117  {
118  args.PrintOptions(cout);
119  }
120 
121  // 3. Read the (serial) mesh from the given mesh file on all processors. We
122  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
123  // and volume meshes with the same code.
124  Mesh *mesh = new Mesh(mesh_file, 1, 1);
125  int dim = mesh->Dimension();
126 
127  // 4. Refine the serial mesh on all processors to increase the resolution. In
128  // this example we do 'ref_levels' of uniform refinement (2 by default, or
129  // specified on the command line with -rs).
130  for (int lev = 0; lev < ser_ref_levels; lev++)
131  {
132  mesh->UniformRefinement();
133  }
134 
135  // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
136  // this mesh further in parallel to increase the resolution (1 time by
137  // default, or specified on the command line with -rp). Once the parallel
138  // mesh is defined, the serial mesh can be deleted.
139  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
140  delete mesh;
141  for (int lev = 0; lev < par_ref_levels; lev++)
142  {
143  pmesh->UniformRefinement();
144  }
145 
146  // 6. Define a parallel finite element space on the parallel mesh. Here we
147  // use continuous Lagrange finite elements of the specified order. If
148  // order < 1, we instead use an isoparametric/isogeometric space.
150  NURBSExtension *NURBSext = NULL;
151  int own_fec = 0;
152 
153  if (order[0] == 0) // Isoparametric
154  {
155  if (pmesh->GetNodes())
156  {
157  fec = pmesh->GetNodes()->OwnFEC();
158  own_fec = 0;
159  if (myid == 0)
160  {
161  cout << "Using isoparametric FEs: " << fec->Name() << endl;
162  }
163  }
164  else
165  {
166  if (myid == 0)
167  {
168  cout <<"Mesh does not have FEs --> Assume order 1.\n";
169  }
170  fec = new H1_FECollection(1, dim);
171  own_fec = 1;
172  }
173  }
174  else if (pmesh->NURBSext && (order[0] > 0) ) // Subparametric NURBS
175  {
176  fec = new NURBSFECollection(order[0]);
177  own_fec = 1;
178  int nkv = pmesh->NURBSext->GetNKV();
179 
180  if (order.Size() == 1)
181  {
182  int tmp = order[0];
183  order.SetSize(nkv);
184  order = tmp;
185  }
186  if (order.Size() != nkv ) { mfem_error("Wrong number of orders set."); }
187  NURBSext = new NURBSExtension(pmesh->NURBSext, order);
188  }
189  else
190  {
191  if (order.Size() > 1) { cout <<"Wrong number of orders set, needs one.\n"; }
192  fec = new H1_FECollection(abs(order[0]), dim);
193  own_fec = 1;
194  }
195  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh,NURBSext,fec);
196  HYPRE_Int size = fespace->GlobalTrueVSize();
197  if (myid == 0)
198  {
199  cout << "Number of unknowns: " << size << endl;
200  }
201 
202  // 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
203  // element space. The first corresponds to the Laplacian operator -Delta,
204  // while the second is a simple mass matrix needed on the right hand side
205  // of the generalized eigenvalue problem below. The boundary conditions
206  // are implemented by elimination with special values on the diagonal to
207  // shift the Dirichlet eigenvalues out of the computational range. After
208  // serial and parallel assembly we extract the corresponding parallel
209  // matrices A and M.
210  ConstantCoefficient one(1.0);
211  Array<int> ess_bdr;
212  if (pmesh->bdr_attributes.Size())
213  {
214  ess_bdr.SetSize(pmesh->bdr_attributes.Max());
215  ess_bdr = 1;
216  }
217 
218  ParBilinearForm *a = new ParBilinearForm(fespace);
220  if (pmesh->bdr_attributes.Size() == 0)
221  {
222  // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
223  // closed surface.
224  a->AddDomainIntegrator(new MassIntegrator(one));
225  }
226  a->Assemble();
227  a->EliminateEssentialBCDiag(ess_bdr, 1.0);
228  a->Finalize();
229 
230  ParBilinearForm *m = new ParBilinearForm(fespace);
231  m->AddDomainIntegrator(new MassIntegrator(one));
232  m->Assemble();
233  // shift the eigenvalue corresponding to eliminated dofs to a large value
234  m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
235  m->Finalize();
236 
239 
240 #if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK)
241  Operator * Arow = NULL;
242 #ifdef MFEM_USE_SUPERLU
243  if (slu_solver)
244  {
245  Arow = new SuperLURowLocMatrix(*A);
246  }
247 #endif
248 #ifdef MFEM_USE_STRUMPACK
249  if (sp_solver)
250  {
251  Arow = new STRUMPACKRowLocMatrix(*A);
252  }
253 #endif
254 #endif
255 
256  delete a;
257  delete m;
258 
259  // 8. Define and configure the LOBPCG eigensolver and the BoomerAMG
260  // preconditioner for A to be used within the solver. Set the matrices
261  // which define the generalized eigenproblem A x = lambda M x.
262  Solver * precond = NULL;
263  if (!slu_solver && !sp_solver)
264  {
265  HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
266  amg->SetPrintLevel(0);
267  precond = amg;
268  }
269  else
270  {
271 #ifdef MFEM_USE_SUPERLU
272  if (slu_solver)
273  {
274  SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
275  superlu->SetPrintStatistics(false);
276  superlu->SetSymmetricPattern(true);
278  superlu->SetOperator(*Arow);
279  precond = superlu;
280  }
281 #endif
282 #ifdef MFEM_USE_STRUMPACK
283  if (sp_solver)
284  {
285  STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
286  strumpack->SetPrintFactorStatistics(true);
287  strumpack->SetPrintSolveStatistics(false);
288  strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
289  strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
290  strumpack->DisableMatching();
291  strumpack->SetOperator(*Arow);
292  strumpack->SetFromCommandLine();
293  precond = strumpack;
294  }
295 #endif
296  }
297 
298 
299  HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
300  lobpcg->SetNumModes(nev);
301  lobpcg->SetRandomSeed(seed);
302  lobpcg->SetPreconditioner(*precond);
303  lobpcg->SetMaxIter(200);
304  lobpcg->SetTol(1e-8);
305  lobpcg->SetPrecondUsageMode(1);
306  lobpcg->SetPrintLevel(1);
307  lobpcg->SetMassMatrix(*M);
308  lobpcg->SetOperator(*A);
309 
310  // 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
311  // parallel grid function to represent each of the eigenmodes returned by
312  // the solver.
313  Array<double> eigenvalues;
314  lobpcg->Solve();
315  lobpcg->GetEigenvalues(eigenvalues);
316  ParGridFunction x(fespace);
317 
318  // 10. Save the refined mesh and the modes in parallel. This output can be
319  // viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
320  {
321  ostringstream mesh_name, mode_name;
322  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
323 
324  ofstream mesh_ofs(mesh_name.str().c_str());
325  mesh_ofs.precision(8);
326  pmesh->Print(mesh_ofs);
327 
328  for (int i=0; i<nev; i++)
329  {
330  // convert eigenvector from HypreParVector to ParGridFunction
331  x = lobpcg->GetEigenvector(i);
332 
333  mode_name << "mode_" << setfill('0') << setw(2) << i << "."
334  << setfill('0') << setw(6) << myid;
335 
336  ofstream mode_ofs(mode_name.str().c_str());
337  mode_ofs.precision(8);
338  x.Save(mode_ofs);
339  mode_name.str("");
340  }
341  }
342 
343  // 11. Send the solution by socket to a GLVis server.
344  if (visualization)
345  {
346  char vishost[] = "localhost";
347  int visport = 19916;
348  socketstream mode_sock(vishost, visport);
349  mode_sock.precision(8);
350 
351  for (int i=0; i<nev; i++)
352  {
353  if ( myid == 0 )
354  {
355  cout << "Eigenmode " << i+1 << '/' << nev
356  << ", Lambda = " << eigenvalues[i] << endl;
357  }
358 
359  // convert eigenvector from HypreParVector to ParGridFunction
360  x = lobpcg->GetEigenvector(i);
361 
362  mode_sock << "parallel " << num_procs << " " << myid << "\n"
363  << "solution\n" << *pmesh << x << flush
364  << "window_title 'Eigenmode " << i+1 << '/' << nev
365  << ", Lambda = " << eigenvalues[i] << "'" << endl;
366 
367  char c;
368  if (myid == 0)
369  {
370  cout << "press (q)uit or (c)ontinue --> " << flush;
371  cin >> c;
372  }
373  MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
374 
375  if (c != 'c')
376  {
377  break;
378  }
379  }
380  mode_sock.close();
381  }
382 
383  // 12. Free the used memory.
384  delete lobpcg;
385  delete precond;
386  delete M;
387  delete A;
388 #if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK)
389  delete Arow;
390 #endif
391 
392  delete fespace;
393  if (own_fec)
394  {
395  delete fec;
396  }
397  delete pmesh;
398 
399  MPI_Finalize();
400 
401  return 0;
402 }
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition: fe_coll.hpp:372
int Size() const
Return the logical size of the array.
Definition: array.hpp:124
HypreParVector & GetEigenvector(unsigned int i)
Extract a single eigenvector.
Definition: hypre.cpp:4335
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
int GetNKV() const
Definition: nurbs.hpp:350
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:819
void SetPreconditioner(Solver &precond)
Definition: hypre.cpp:4258
void SetMassMatrix(Operator &M)
Definition: hypre.cpp:4313
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void SetPrintLevel(int logging)
Definition: hypre.cpp:4243
int main(int argc, char *argv[])
Definition: ex1.cpp:66
void SetSymmetricPattern(bool sym)
Definition: superlu.cpp:420
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void SetTol(double tol)
Definition: hypre.cpp:4221
int close()
Close the socketstream.
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: superlu.cpp:576
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1079
void SetPrintFactorStatistics(bool print_stat)
Definition: strumpack.cpp:134
void SetMaxIter(int max_iter)
Definition: hypre.cpp:4237
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
constexpr char vishost[]
void SetColumnPermutation(superlu::ColPerm col_perm)
Definition: superlu.cpp:340
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:153
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:8382
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1119
constexpr int visport
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:4169
void SetPrintSolveStatistics(bool print_stat)
Definition: strumpack.cpp:139
void SetPrintStatistics(bool print_stat)
Definition: superlu.cpp:322
int Dimension() const
Definition: mesh.hpp:788
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:434
void SetReorderingStrategy(strumpack::ReorderingStrategy method)
Definition: strumpack.cpp:149
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:201
void GetEigenvalues(Array< double > &eigenvalues)
Collect the converged eigenvalues.
Definition: hypre.cpp:4323
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
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:654
void SetRandomSeed(int s)
Definition: hypre.hpp:1337
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:203
virtual const char * Name() const
Definition: fe_coll.hpp:61
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
int dim
Definition: ex24.cpp:53
void SetKrylovSolver(strumpack::KrylovSolver method)
Definition: strumpack.cpp:144
void SetOperator(Operator &A)
Definition: hypre.cpp:4267
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:304
Class for parallel bilinear form.
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: strumpack.cpp:228
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6603
void SetPrecondUsageMode(int pcg_mode)
Definition: hypre.cpp:4252
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:159
Base class for solvers.
Definition: operator.hpp:634
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:181
void SetNumModes(int num_eigs)
Definition: hypre.hpp:1335
Class for parallel meshes.
Definition: pmesh.hpp:32
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:4380
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:145