MFEM  v4.5.2
Finite element discretization library
ex0p.cpp
Go to the documentation of this file.
1 // MFEM Example 0 - Parallel Version
2 //
3 // Compile with: make ex0p
4 //
5 // Sample runs: mpirun -np 4 ex0p
6 // mpirun -np 4 ex0p -m ../data/fichera.mesh
7 // mpirun -np 4 ex0p -m ../data/square-disc.mesh -o 2
8 //
9 // Description: This example code demonstrates the most basic parallel usage of
10 // MFEM to define a simple finite element discretization of the
11 // Laplace problem -Delta u = 1 with zero Dirichlet boundary
12 // conditions. General 2D/3D serial mesh files and finite element
13 // polynomial degrees can be specified by command line options.
14 
15 #include "mfem.hpp"
16 #include <fstream>
17 #include <iostream>
18 
19 using namespace std;
20 using namespace mfem;
21 
22 int main(int argc, char *argv[])
23 {
24  // 1. Initialize MPI and HYPRE.
25  Mpi::Init(argc, argv);
26  Hypre::Init();
27 
28  // 2. Parse command line options.
29  const char *mesh_file = "../data/star.mesh";
30  int order = 1;
31 
32  OptionsParser args(argc, argv);
33  args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
34  args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
35  args.ParseCheck();
36 
37  // 3. Read the serial mesh from the given mesh file.
38  Mesh serial_mesh(mesh_file);
39 
40  // 4. Define a parallel mesh by a partitioning of the serial mesh. Refine
41  // this mesh once in parallel to increase the resolution.
42  ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
43  serial_mesh.Clear(); // the serial mesh is no longer needed
44  mesh.UniformRefinement();
45 
46  // 5. Define a finite element space on the mesh. Here we use H1 continuous
47  // high-order Lagrange finite elements of the given order.
48  H1_FECollection fec(order, mesh.Dimension());
49  ParFiniteElementSpace fespace(&mesh, &fec);
50  HYPRE_BigInt total_num_dofs = fespace.GlobalTrueVSize();
51  if (Mpi::Root())
52  {
53  cout << "Number of unknowns: " << total_num_dofs << endl;
54  }
55 
56  // 6. Extract the list of all the boundary DOFs. These will be marked as
57  // Dirichlet in order to enforce zero boundary conditions.
58  Array<int> boundary_dofs;
59  fespace.GetBoundaryTrueDofs(boundary_dofs);
60 
61  // 7. Define the solution x as a finite element grid function in fespace. Set
62  // the initial guess to zero, which also sets the boundary conditions.
63  ParGridFunction x(&fespace);
64  x = 0.0;
65 
66  // 8. Set up the linear form b(.) corresponding to the right-hand side.
67  ConstantCoefficient one(1.0);
68  ParLinearForm b(&fespace);
69  b.AddDomainIntegrator(new DomainLFIntegrator(one));
70  b.Assemble();
71 
72  // 9. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
73  ParBilinearForm a(&fespace);
74  a.AddDomainIntegrator(new DiffusionIntegrator);
75  a.Assemble();
76 
77  // 10. Form the linear system A X = B. This includes eliminating boundary
78  // conditions, applying AMR constraints, parallel assembly, etc.
80  Vector B, X;
81  a.FormLinearSystem(boundary_dofs, x, b, A, X, B);
82 
83  // 11. Solve the system using PCG with hypre's BoomerAMG preconditioner.
84  HypreBoomerAMG M(A);
85  CGSolver cg(MPI_COMM_WORLD);
86  cg.SetRelTol(1e-12);
87  cg.SetMaxIter(2000);
88  cg.SetPrintLevel(1);
89  cg.SetPreconditioner(M);
90  cg.SetOperator(A);
91  cg.Mult(B, X);
92 
93  // 12. Recover the solution x as a grid function and save to file. The output
94  // can be viewed using GLVis as follows: "glvis -np <np> -m mesh -g sol"
95  a.RecoverFEMSolution(X, b, x);
96  x.Save("sol");
97  mesh.Save("mesh");
98 
99  return 0;
100 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
Conjugate gradient method.
Definition: solvers.hpp:493
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
int main(int argc, char *argv[])
Definition: ex0p.cpp:22
int Dimension() const
Definition: mesh.hpp:1047
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:712
Abstract parallel finite element space.
Definition: pfespace.hpp:28
STL namespace.
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
Class for parallel linear form.
Definition: plinearform.hpp:26
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9878
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
void SetRelTol(double rtol)
Definition: solvers.hpp:199
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
HYPRE_Int HYPRE_BigInt
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:873
double a
Definition: lissajous.cpp:41
Class for parallel bilinear form.
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:912
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
Vector data type.
Definition: vector.hpp:60
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Class for parallel meshes.
Definition: pmesh.hpp:32
void Save(const char *fname, int precision=16) const override
Definition: pmesh.cpp:4955
void ParseCheck(std::ostream &out=mfem::out)
Definition: optparser.cpp:252