MFEM  v3.3
Finite element discretization library
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-tri.mesh -o 2 -sys
10 // mpirun -np 4 ex2p -m ../data/beam-quad.mesh -o 3 -elast
11 // mpirun -np 4 ex2p -m ../data/beam-quad.mesh -o 3 -sc
12 // mpirun -np 4 ex2p -m ../data/beam-quad-nurbs.mesh
13 // mpirun -np 4 ex2p -m ../data/beam-hex-nurbs.mesh
14 //
15 // Description: This example code solves a simple linear elasticity problem
16 // describing a multi-material cantilever beam.
17 //
18 // Specifically, we approximate the weak form of -div(sigma(u))=0
19 // where 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 with f being a constant pull down
24 // vector on boundary elements with attribute 2, and zero
25 // otherwise. The geometry of the domain is assumed to be as
26 // follows:
27 //
28 // +----------+----------+
29 // boundary --->| material | material |<--- boundary
30 // attribute 1 | 1 | 2 | attribute 2
31 // (fixed) +----------+----------+ (pull down)
32 //
33 // The example demonstrates the use of high-order and NURBS vector
34 // finite element spaces with the linear elasticity bilinear form,
35 // meshes with curved elements, and the definition of piece-wise
36 // constant and vector coefficient objects. Static condensation is
37 // 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/beam-tri.mesh";
58  int order = 1;
59  bool static_cond = false;
60  bool visualization = 1;
61  bool amg_elast = 0;
62 
63  OptionsParser args(argc, argv);
64  args.AddOption(&mesh_file, "-m", "--mesh",
65  "Mesh file to use.");
66  args.AddOption(&order, "-o", "--order",
67  "Finite element order (polynomial degree).");
68  args.AddOption(&amg_elast, "-elast", "--amg-for-elasticity", "-sys",
69  "--amg-for-systems",
70  "Use the special AMG elasticity solver (GM/LN approaches), "
71  "or standard AMG for systems (unknown approach).");
72  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
73  "--no-static-condensation", "Enable static condensation.");
74  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
75  "--no-visualization",
76  "Enable or disable GLVis visualization.");
77  args.Parse();
78  if (!args.Good())
79  {
80  if (myid == 0)
81  {
82  args.PrintUsage(cout);
83  }
84  MPI_Finalize();
85  return 1;
86  }
87  if (myid == 0)
88  {
89  args.PrintOptions(cout);
90  }
91 
92  // 3. Read the (serial) mesh from the given mesh file on all processors. We
93  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
94  // and volume meshes with the same code.
95  Mesh *mesh = new Mesh(mesh_file, 1, 1);
96  int dim = mesh->Dimension();
97 
98  if (mesh->attributes.Max() < 2 || mesh->bdr_attributes.Max() < 2)
99  {
100  if (myid == 0)
101  cerr << "\nInput mesh should have at least two materials and "
102  << "two boundary attributes! (See schematic in ex2.cpp)\n"
103  << endl;
104  MPI_Finalize();
105  return 3;
106  }
107 
108  // 4. Select the order of the finite element discretization space. For NURBS
109  // meshes, we increase the order by degree elevation.
110  if (mesh->NURBSext && order > mesh->NURBSext->GetOrder())
111  {
112  mesh->DegreeElevate(order - mesh->NURBSext->GetOrder());
113  }
114 
115  // 5. Refine the serial mesh on all processors to increase the resolution. In
116  // this example we do 'ref_levels' of uniform refinement. We choose
117  // 'ref_levels' to be the largest number that gives a final mesh with no
118  // more than 1,000 elements.
119  {
120  int ref_levels =
121  (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
122  for (int l = 0; l < ref_levels; l++)
123  {
124  mesh->UniformRefinement();
125  }
126  }
127 
128  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
129  // this mesh further in parallel to increase the resolution. Once the
130  // parallel mesh is defined, the serial mesh can be deleted.
131  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
132  delete mesh;
133  {
134  int par_ref_levels = 1;
135  for (int l = 0; l < par_ref_levels; l++)
136  {
137  pmesh->UniformRefinement();
138  }
139  }
140 
141  // 7. Define a parallel finite element space on the parallel mesh. Here we
142  // use vector finite elements, i.e. dim copies of a scalar finite element
143  // space. We use the ordering by vector dimension (the last argument of
144  // the FiniteElementSpace constructor) which is expected in the systems
145  // version of BoomerAMG preconditioner. For NURBS meshes, we use the
146  // (degree elevated) NURBS space associated with the mesh nodes.
148  ParFiniteElementSpace *fespace;
149  const bool use_nodal_fespace = pmesh->NURBSext && !amg_elast;
150  if (use_nodal_fespace)
151  {
152  fec = NULL;
153  fespace = (ParFiniteElementSpace *)pmesh->GetNodes()->FESpace();
154  }
155  else
156  {
157  fec = new H1_FECollection(order, dim);
158  fespace = new ParFiniteElementSpace(pmesh, fec, dim, Ordering::byVDIM);
159  }
160  HYPRE_Int size = fespace->GlobalTrueVSize();
161  if (myid == 0)
162  {
163  cout << "Number of finite element unknowns: " << size << endl
164  << "Assembling: " << flush;
165  }
166 
167  // 8. Determine the list of true (i.e. parallel conforming) essential
168  // boundary dofs. In this example, the boundary conditions are defined by
169  // marking only boundary attribute 1 from the mesh as essential and
170  // converting it to a list of true dofs.
171  Array<int> ess_tdof_list, ess_bdr(pmesh->bdr_attributes.Max());
172  ess_bdr = 0;
173  ess_bdr[0] = 1;
174  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
175 
176  // 9. Set up the parallel linear form b(.) which corresponds to the
177  // right-hand side of the FEM linear system. In this case, b_i equals the
178  // boundary integral of f*phi_i where f represents a "pull down" force on
179  // the Neumann part of the boundary and phi_i are the basis functions in
180  // the finite element fespace. The force is defined by the object f, which
181  // is a vector of Coefficient objects. The fact that f is non-zero on
182  // boundary attribute 2 is indicated by the use of piece-wise constants
183  // coefficient for its last component.
184  VectorArrayCoefficient f(dim);
185  for (int i = 0; i < dim-1; i++)
186  {
187  f.Set(i, new ConstantCoefficient(0.0));
188  }
189  {
190  Vector pull_force(pmesh->bdr_attributes.Max());
191  pull_force = 0.0;
192  pull_force(1) = -1.0e-2;
193  f.Set(dim-1, new PWConstCoefficient(pull_force));
194  }
195 
196  ParLinearForm *b = new ParLinearForm(fespace);
198  if (myid == 0)
199  {
200  cout << "r.h.s. ... " << flush;
201  }
202  b->Assemble();
203 
204  // 10. Define the solution vector x as a parallel finite element grid
205  // function corresponding to fespace. Initialize x with initial guess of
206  // zero, which satisfies the boundary conditions.
207  ParGridFunction x(fespace);
208  x = 0.0;
209 
210  // 11. Set up the parallel bilinear form a(.,.) on the finite element space
211  // corresponding to the linear elasticity integrator with piece-wise
212  // constants coefficient lambda and mu.
213  Vector lambda(pmesh->attributes.Max());
214  lambda = 1.0;
215  lambda(0) = lambda(1)*50;
216  PWConstCoefficient lambda_func(lambda);
217  Vector mu(pmesh->attributes.Max());
218  mu = 1.0;
219  mu(0) = mu(1)*50;
220  PWConstCoefficient mu_func(mu);
221 
222  ParBilinearForm *a = new ParBilinearForm(fespace);
223  a->AddDomainIntegrator(new ElasticityIntegrator(lambda_func, mu_func));
224 
225  // 12. Assemble the parallel bilinear form and the corresponding linear
226  // system, applying any necessary transformations such as: parallel
227  // assembly, eliminating boundary conditions, applying conforming
228  // constraints for non-conforming AMR, static condensation, etc.
229  if (myid == 0) { cout << "matrix ... " << flush; }
230  if (static_cond) { a->EnableStaticCondensation(); }
231  a->Assemble();
232 
233  HypreParMatrix A;
234  Vector B, X;
235  a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
236  if (myid == 0)
237  {
238  cout << "done." << endl;
239  cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
240  }
241 
242  // 13. Define and apply a parallel PCG solver for A X = B with the BoomerAMG
243  // preconditioner from hypre.
244  HypreBoomerAMG *amg = new HypreBoomerAMG(A);
245  if (amg_elast && !a->StaticCondensationIsEnabled())
246  {
247  amg->SetElasticityOptions(fespace);
248  }
249  else
250  {
251  amg->SetSystemsOptions(dim);
252  }
253  HyprePCG *pcg = new HyprePCG(A);
254  pcg->SetTol(1e-8);
255  pcg->SetMaxIter(500);
256  pcg->SetPrintLevel(2);
257  pcg->SetPreconditioner(*amg);
258  pcg->Mult(B, X);
259 
260  // 14. Recover the parallel grid function corresponding to X. This is the
261  // local finite element solution on each processor.
262  a->RecoverFEMSolution(X, *b, x);
263 
264  // 15. For non-NURBS meshes, make the mesh curved based on the finite element
265  // space. This means that we define the mesh elements through a fespace
266  // based transformation of the reference element. This allows us to save
267  // the displaced mesh as a curved mesh when using high-order finite
268  // element displacement field. We assume that the initial mesh (read from
269  // the file) is not higher order curved mesh compared to the chosen FE
270  // space.
271  if (!use_nodal_fespace)
272  {
273  pmesh->SetNodalFESpace(fespace);
274  }
275 
276  // 16. Save in parallel the displaced mesh and the inverted solution (which
277  // gives the backward displacements to the original grid). This output
278  // can be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
279  {
280  GridFunction *nodes = pmesh->GetNodes();
281  *nodes += x;
282  x *= -1;
283 
284  ostringstream mesh_name, sol_name;
285  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
286  sol_name << "sol." << setfill('0') << setw(6) << myid;
287 
288  ofstream mesh_ofs(mesh_name.str().c_str());
289  mesh_ofs.precision(8);
290  pmesh->Print(mesh_ofs);
291 
292  ofstream sol_ofs(sol_name.str().c_str());
293  sol_ofs.precision(8);
294  x.Save(sol_ofs);
295  }
296 
297  // 17. Send the above data by socket to a GLVis server. Use the "n" and "b"
298  // keys in GLVis to visualize the displacements.
299  if (visualization)
300  {
301  char vishost[] = "localhost";
302  int visport = 19916;
303  socketstream sol_sock(vishost, visport);
304  sol_sock << "parallel " << num_procs << " " << myid << "\n";
305  sol_sock.precision(8);
306  sol_sock << "solution\n" << *pmesh << x << flush;
307  }
308 
309  // 18. Free the used memory.
310  delete pcg;
311  delete amg;
312  delete a;
313  delete b;
314  if (fec)
315  {
316  delete fespace;
317  delete fec;
318  }
319  delete pmesh;
320 
321  MPI_Finalize();
322 
323  return 0;
324 }
void SetTol(double tol)
Definition: hypre.cpp:1993
Vector coefficient defined by an array of scalar coefficients.
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
Subclass constant coefficient.
Definition: coefficient.hpp:57
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:42
void DegreeElevate(int t)
Definition: mesh.cpp:3038
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:584
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:314
Abstract parallel finite element space.
Definition: pfespace.hpp:28
STL namespace.
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2008
void SetSystemsOptions(int dim)
Definition: hypre.cpp:2407
int main(int argc, char *argv[])
Definition: ex2p.cpp:48
The BoomerAMG solver in hypre.
Definition: hypre.hpp:770
Class for parallel linear form.
Definition: plinearform.hpp:26
HYPRE_Int GetGlobalNumRows() const
Definition: hypre.hpp:371
int dim
Definition: ex3.cpp:47
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6684
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:108
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void SetElasticityOptions(ParFiniteElementSpace *fespace)
Definition: hypre.cpp:2489
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3261
void AddBoundaryIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Integrator.
Definition: linearform.cpp:24
int Dimension() const
Definition: mesh.hpp:611
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetMaxIter(int max_iter)
Definition: hypre.cpp:1998
bool StaticCondensationIsEnabled() const
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:142
PCG solver in hypre.
Definition: hypre.hpp:630
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list)
Definition: pfespace.cpp:521
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
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:144
class for piecewise constant coefficient
Definition: coefficient.hpp:72
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:2013
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:36
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5368
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:146
Class for parallel grid function.
Definition: pgridfunc.hpp:31
Wrapper for hypre's ParCSR matrix class.
Definition: hypre.hpp:174
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre's PCG.
Definition: hypre.cpp:2034
Class for parallel meshes.
Definition: pmesh.hpp:28
void Set(int i, Coefficient *c)
Sets coefficient in the vector.
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:140
void EnableStaticCondensation()
virtual void Print(std::ostream &out=std::cout) const
Definition: pmesh.cpp:3395
bool Good() const
Definition: optparser.hpp:120