MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
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
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
46using namespace std;
47using namespace mfem;
48
49int 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();
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;
65
66 OptionsParser args(argc, argv);
68 "Mesh file to use.");
70 "Finite element order (polynomial degree).");
72 "Number of desired eigenmodes.");
74 "Random seed used to initialize LOBPCG.");
76 "--amg-for-systems",
77 "Use the special AMG elasticity solver (GM/LN approaches), "
78 "or standard AMG for systems (unknown approach).");
80 "--no-visualization",
81 "Enable or disable GLVis visualization.");
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);
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<real_t>::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 {
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<real_t> 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.
296 {
297 std::string postfix(mesh_file);
298 postfix.erase(0, std::string("../data/").size() );
299 postfix += "_o" + std::to_string(order);
300
304 for (int i=0; i<nev; i++)
305 {
306 x = lobpcg->GetEigenvector(i);
307 // x is a temporary that must be saved immediately
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}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
Adds new Domain Integrator. Assumes ownership of bfi.
void EliminateEssentialBCDiag(const Array< int > &bdr_attr_is_ess, real_t value)
Perform elimination and set the diagonal entry to the given value.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetSystemsOptions(int dim, bool order_bynodes=false)
Definition hypre.cpp:5111
void SetElasticityOptions(ParFiniteElementSpace *fespace, bool interp_refine=true)
Definition hypre.cpp:5238
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
void SetMassMatrix(Operator &M)
Definition hypre.cpp:6343
void SetPrintLevel(int logging)
Definition hypre.cpp:6272
void SetPreconditioner(Solver &precond)
Definition hypre.cpp:6287
void GetEigenvalues(Array< real_t > &eigenvalues) const
Collect the converged eigenvalues.
Definition hypre.cpp:6353
void SetTol(real_t tol)
Definition hypre.cpp:6250
void SetOperator(Operator &A)
Definition hypre.cpp:6296
void SetNumModes(int num_eigs)
Definition hypre.hpp:2103
void Solve()
Solve the eigenproblem.
Definition hypre.cpp:6410
void SetPrecondUsageMode(int pcg_mode)
Definition hypre.cpp:6281
void SetRandomSeed(int s)
Definition hypre.hpp:2105
void SetMaxIter(int max_iter)
Definition hypre.cpp:6266
const HypreParVector & GetEigenvector(unsigned int i) const
Extract a single eigenvector.
Definition hypre.cpp:6365
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:290
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
void DegreeElevate(int rel_degree, int degree=16)
Definition mesh.cpp:5779
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
A piecewise constant coefficient with the constants keyed off the element attribute numbers.
Class for parallel bilinear form.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Save(std::ostream &out) const override
Class for parallel meshes.
Definition pmesh.hpp:34
void SetNodalFESpace(FiniteElementSpace *nfes) override
Definition pmesh.cpp:2028
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Vector data type.
Definition vector.hpp:80
int close()
Close the socketstream.
int dim
Definition ex24.cpp:53
real_t mu
Definition ex25.cpp:140
int main()
HYPRE_Int HYPRE_BigInt
real_t a
Definition lissajous.cpp:41
const int visport
const char vishost[]