MFEM  v3.0
ex2p.cpp
Go to the documentation of this file.
1 // MFEM Example 2 - Parallel Version
2 //
3 // Compile with: make ex2p
4 //
5 // Sample runs: mpirun -np 4 ex2p -m ../data/beam-tri.mesh
6 // mpirun -np 4 ex2p -m ../data/beam-quad.mesh
7 // mpirun -np 4 ex2p -m ../data/beam-tet.mesh
8 // mpirun -np 4 ex2p -m ../data/beam-hex.mesh
9 // mpirun -np 4 ex2p -m ../data/beam-quad-nurbs.mesh
10 // mpirun -np 4 ex2p -m ../data/beam-hex-nurbs.mesh
11 //
12 // Description: This example code solves a simple linear elasticity problem
13 // describing a multi-material cantilever beam.
14 //
15 // Specifically, we approximate the weak form of -div(sigma(u))=0
17 // tensor corresponding to displacement field u, and lambda and mu
18 // are the material Lame constants. The boundary conditions are
19 // u=0 on the fixed part of the boundary with attribute 1, and
20 // sigma(u).n=f on the remainder with f being a constant pull down
21 // vector on boundary elements with attribute 2, and zero
22 // otherwise. The geometry of the domain is assumed to be as
23 // follows:
24 //
25 // +----------+----------+
26 // boundary --->| material | material |<--- boundary
27 // attribute 1 | 1 | 2 | attribute 2
28 // (fixed) +----------+----------+ (pull down)
29 //
30 // The example demonstrates the use of high-order and NURBS vector
31 // finite element spaces with the linear elasticity bilinear form,
32 // meshes with curved elements, and the definition of piece-wise
33 // constant and vector coefficient objects.
34 //
35 // We recommend viewing Example 1 before viewing this example.
36
37 #include "mfem.hpp"
38 #include <fstream>
39 #include <iostream>
40
41 using namespace std;
42 using namespace mfem;
43
44 int main(int argc, char *argv[])
45 {
46  // 1. Initialize MPI.
47  int num_procs, myid;
48  MPI_Init(&argc, &argv);
49  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
50  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
51
52  // 2. Parse command-line options.
53  const char *mesh_file = "../data/beam-tri.mesh";
54  int order = 1;
55  bool visualization = 1;
56
57  OptionsParser args(argc, argv);
59  "Mesh file to use.");
61  "Finite element order (polynomial degree).");
63  "--no-visualization",
64  "Enable or disable GLVis visualization.");
65  args.Parse();
66  if (!args.Good())
67  {
68  if (myid == 0)
69  args.PrintUsage(cout);
70  MPI_Finalize();
71  return 1;
72  }
73  if (myid == 0)
74  args.PrintOptions(cout);
75
76  // 3. Read the (serial) mesh from the given mesh file on all processors. We
77  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
78  // and volume meshes with the same code.
79  Mesh *mesh;
80  ifstream imesh(mesh_file);
81  if (!imesh)
82  {
83  if (myid == 0)
84  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
85  MPI_Finalize();
86  return 2;
87  }
88  mesh = new Mesh(imesh, 1, 1);
89  imesh.close();
90  int dim = mesh->Dimension();
91
92  if (mesh->attributes.Max() < 2 || mesh->bdr_attributes.Max() < 2)
93  {
94  if (myid == 0)
95  cerr << "\nInput mesh should have at least two materials and "
96  << "two boundary attributes! (See schematic in ex2.cpp)\n"
97  << endl;
98  MPI_Finalize();
99  return 3;
100  }
101
102  // 4. Select the order of the finite element discretization space. For NURBS
103  // meshes, we increase the order by degree elevation.
104  if (mesh->NURBSext && order > mesh->NURBSext->GetOrder())
105  mesh->DegreeElevate(order - mesh->NURBSext->GetOrder());
106
107  // 5. Refine the serial mesh on all processors to increase the resolution. In
108  // this example we do 'ref_levels' of uniform refinement. We choose
109  // 'ref_levels' to be the largest number that gives a final mesh with no
110  // more than 1,000 elements.
111  {
112  int ref_levels =
113  (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
114  for (int l = 0; l < ref_levels; l++)
115  mesh->UniformRefinement();
116  }
117
118  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
119  // this mesh further in parallel to increase the resolution. Once the
120  // parallel mesh is defined, the serial mesh can be deleted.
121  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
122  delete mesh;
123  {
124  int par_ref_levels = 1;
125  for (int l = 0; l < par_ref_levels; l++)
126  pmesh->UniformRefinement();
127  }
128
129  // 7. Define a parallel finite element space on the parallel mesh. Here we
130  // use vector finite elements, i.e. dim copies of a scalar finite element
131  // space. We use the ordering by vector dimension (the last argument of
132  // the FiniteElementSpace constructor) which is expected in the systems
133  // version of BoomerAMG preconditioner. For NURBS meshes, we use the
134  // (degree elevated) NURBS space associated with the mesh nodes.
136  ParFiniteElementSpace *fespace;
137  if (pmesh->NURBSext)
138  {
139  fec = NULL;
140  fespace = (ParFiniteElementSpace *)pmesh->GetNodes()->FESpace();
141  }
142  else
143  {
144  fec = new H1_FECollection(order, dim);
145  fespace = new ParFiniteElementSpace(pmesh, fec, dim, Ordering::byVDIM);
146  }
147  int size = fespace->GlobalTrueVSize();
148  if (myid == 0)
149  cout << "Number of unknowns: " << size << endl
150  << "Assembling: " << flush;
151
152  // 8. Set up the parallel linear form b(.) which corresponds to the
153  // right-hand side of the FEM linear system. In this case, b_i equals the
154  // boundary integral of f*phi_i where f represents a "pull down" force on
155  // the Neumann part of the boundary and phi_i are the basis functions in
156  // the finite element fespace. The force is defined by the object f, which
157  // is a vector of Coefficient objects. The fact that f is non-zero on
158  // boundary attribute 2 is indicated by the use of piece-wise constants
159  // coefficient for its last component.
160  VectorArrayCoefficient f(dim);
161  for (int i = 0; i < dim-1; i++)
162  f.Set(i, new ConstantCoefficient(0.0));
163  {
164  Vector pull_force(pmesh->bdr_attributes.Max());
165  pull_force = 0.0;
166  pull_force(1) = -1.0e-2;
167  f.Set(dim-1, new PWConstCoefficient(pull_force));
168  }
169
170  ParLinearForm *b = new ParLinearForm(fespace);
172  if (myid == 0)
173  cout << "r.h.s. ... " << flush;
174  b->Assemble();
175
176  // 9. Define the solution vector x as a parallel finite element grid function
177  // corresponding to fespace. Initialize x with initial guess of zero,
178  // which satisfies the boundary conditions.
179  ParGridFunction x(fespace);
180  x = 0.0;
181
182  // 10. Set up the parallel bilinear form a(.,.) on the finite element space
183  // corresponding to the linear elasticity integrator with piece-wise
184  // constants coefficient lambda and mu. The boundary conditions are
185  // implemented by marking only boundary attribute 1 as essential. After
186  // serial/parallel assembly we extract the corresponding parallel matrix.
187  Vector lambda(pmesh->attributes.Max());
188  lambda = 1.0;
189  lambda(0) = lambda(1)*50;
190  PWConstCoefficient lambda_func(lambda);
191  Vector mu(pmesh->attributes.Max());
192  mu = 1.0;
193  mu(0) = mu(1)*50;
194  PWConstCoefficient mu_func(mu);
195
196  ParBilinearForm *a = new ParBilinearForm(fespace);
198  if (myid == 0)
199  cout << "matrix ... " << flush;
200  a->Assemble();
201  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
202  ess_bdr = 0;
203  ess_bdr[0] = 1;
204  a->EliminateEssentialBC(ess_bdr, x, *b);
205  a->Finalize();
206  if (myid == 0)
207  cout << "done." << endl;
208
209  // 11. Define the parallel (hypre) matrix and vectors representing a(.,.),
210  // b(.) and the finite element approximation.
214
215  delete a;
216  delete b;
217
218  // 12. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
219  // preconditioner from hypre.
220  HypreBoomerAMG *amg = new HypreBoomerAMG(*A);
221  amg->SetSystemsOptions(dim);
222  HyprePCG *pcg = new HyprePCG(*A);
223  pcg->SetTol(1e-8);
224  pcg->SetMaxIter(500);
225  pcg->SetPrintLevel(2);
226  pcg->SetPreconditioner(*amg);
227  pcg->Mult(*B, *X);
228
229  // 13. Extract the parallel grid function corresponding to the finite element
230  // approximation X. This is the local solution on each processor.
231  x = *X;
232
233  // 14. For non-NURBS meshes, make the mesh curved based on the finite element
234  // space. This means that we define the mesh elements through a fespace
235  // based transformation of the reference element. This allows us to save
236  // the displaced mesh as a curved mesh when using high-order finite
237  // element displacement field. We assume that the initial mesh (read from
238  // the file) is not higher order curved mesh compared to the chosen FE
239  // space.
240  if (!pmesh->NURBSext)
241  pmesh->SetNodalFESpace(fespace);
242
243  // 15. Save in parallel the displaced mesh and the inverted solution (which
244  // gives the backward displacements to the original grid). This output
245  // can be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
246  {
247  GridFunction *nodes = pmesh->GetNodes();
248  *nodes += x;
249  x *= -1;
250
251  ostringstream mesh_name, sol_name;
252  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
253  sol_name << "sol." << setfill('0') << setw(6) << myid;
254
255  ofstream mesh_ofs(mesh_name.str().c_str());
256  mesh_ofs.precision(8);
257  pmesh->Print(mesh_ofs);
258
259  ofstream sol_ofs(sol_name.str().c_str());
260  sol_ofs.precision(8);
261  x.Save(sol_ofs);
262  }
263
264  // 16. Send the above data by socket to a GLVis server. Use the "n" and "b"
265  // keys in GLVis to visualize the displacements.
266  if (visualization)
267  {
268  char vishost[] = "localhost";
269  int visport = 19916;
270  socketstream sol_sock(vishost, visport);
271  sol_sock << "parallel " << num_procs << " " << myid << "\n";
272  sol_sock.precision(8);
273  sol_sock << "solution\n" << *pmesh << x << flush;
274  }
275
276  // 17. Free the used memory.
277  delete pcg;
278  delete amg;
279  delete X;
280  delete B;
281  delete A;
282
283  if (fec)
284  {
285  delete fespace;
286  delete fec;
287  }
288  delete pmesh;
289
290  MPI_Finalize();
291
292  return 0;
293 }
void SetTol(double tol)
Definition: hypre.cpp:1339
Vector coefficient defined by an array of scalar coefficients.
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:26
Subclass constant coefficient.
Definition: coefficient.hpp:57
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:34
void DegreeElevate(int t)
Definition: mesh.cpp:2878
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:396
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:250
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:1354
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void SetSystemsOptions(int dim)
Definition: hypre.cpp:1628
The BoomerAMG solver in hypre.
Definition: hypre.hpp:519
Class for parallel linear form.
Definition: plinearform.hpp:26
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6225
T Max() const
Definition: array.cpp:78
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3066
Definition: linearform.cpp:24
int Dimension() const
Definition: mesh.hpp:417
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:385
void SetMaxIter(int max_iter)
Definition: hypre.cpp:1344
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:38
Array< int > bdr_attributes
Definition: mesh.hpp:305
int main(int argc, char *argv[])
Definition: ex1.cpp:39
PCG solver in hypre.
Definition: hypre.hpp:381
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:307
class for piecewise constant coefficient
Definition: coefficient.hpp:72
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:1360
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:266
Class for parallel bilinear form.
Vector data type.
Definition: vector.hpp:29
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:4948
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:52
Class for parallel grid function.
Definition: pgridfunc.hpp:31
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:103
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:1377
Class for parallel meshes.
Definition: pmesh.hpp:27
void ParallelAverage(Vector &tv) const
Returns the vector averaged on the true dofs.
Definition: pgridfunc.cpp:110
void Set(int i, Coefficient *c)
Sets coefficient in the vector.
void EliminateEssentialBC(Array< int > &bdr_attr_is_ess, Vector &sol, Vector &rhs, int d=0)
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
Definition: plinearform.cpp:34
Array< int > attributes
Definition: mesh.hpp:304
virtual void Print(std::ostream &out=std::cout) const
Definition: pmesh.cpp:2444
bool Good() const
Definition: optparser.hpp:120