MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex1p.cpp
Go to the documentation of this file.
1 // MFEM Example 1 - Parallel Version
2 //
3 // Compile with: make ex1p
4 //
5 // Sample runs: mpirun -np 4 ex1p -m ../data/square-disc.mesh
6 // mpirun -np 4 ex1p -m ../data/star.mesh
7 // mpirun -np 4 ex1p -m ../data/star-mixed.mesh
8 // mpirun -np 4 ex1p -m ../data/escher.mesh
9 // mpirun -np 4 ex1p -m ../data/fichera.mesh
10 // mpirun -np 4 ex1p -m ../data/fichera-mixed.mesh
11 // mpirun -np 4 ex1p -m ../data/toroid-wedge.mesh
12 // mpirun -np 4 ex1p -m ../data/periodic-annulus-sector.msh
13 // mpirun -np 4 ex1p -m ../data/periodic-torus-sector.msh
14 // mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
15 // mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
16 // mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
17 // mpirun -np 4 ex1p -m ../data/star-mixed-p2.mesh -o 2
18 // mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
19 // mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
20 // mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
21 // mpirun -np 4 ex1p -m ../data/fichera-mixed-p2.mesh -o 2
22 // mpirun -np 4 ex1p -m ../data/star-surf.mesh
23 // mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
24 // mpirun -np 4 ex1p -m ../data/inline-segment.mesh
25 // mpirun -np 4 ex1p -m ../data/amr-quad.mesh
26 // mpirun -np 4 ex1p -m ../data/amr-hex.mesh
27 // mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
28 // mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
29 //
30 // Device sample runs:
31 // mpirun -np 4 ex1p -pa -d cuda
32 // mpirun -np 4 ex1p -pa -d occa-cuda
33 // mpirun -np 4 ex1p -pa -d raja-omp
34 // mpirun -np 4 ex1p -pa -d ceed-cpu
35 // * mpirun -np 4 ex1p -pa -d ceed-cuda
36 // * mpirun -np 4 ex1p -pa -d ceed-hip
37 // mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
38 // mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
39 //
40 // Description: This example code demonstrates the use of MFEM to define a
41 // simple finite element discretization of the Laplace problem
42 // -Delta u = 1 with homogeneous Dirichlet boundary conditions.
43 // Specifically, we discretize using a FE space of the specified
44 // order, or if order < 1 using an isoparametric/isogeometric
45 // space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
46 // NURBS mesh, etc.)
47 //
48 // The example highlights the use of mesh refinement, finite
49 // element grid functions, as well as linear and bilinear forms
50 // corresponding to the left-hand side and right-hand side of the
51 // discrete linear system. We also cover the explicit elimination
52 // of essential boundary conditions, static condensation, and the
53 // optional connection to the GLVis tool for visualization.
54 
55 #include "mfem.hpp"
56 #include <fstream>
57 #include <iostream>
58 
59 using namespace std;
60 using namespace mfem;
61 
62 int main(int argc, char *argv[])
63 {
64  // 1. Initialize MPI.
65  int num_procs, myid;
66  MPI_Init(&argc, &argv);
67  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
68  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
69 
70  // 2. Parse command-line options.
71  const char *mesh_file = "../data/star.mesh";
72  int order = 1;
73  bool static_cond = false;
74  bool pa = false;
75  const char *device_config = "cpu";
76  bool visualization = true;
77 
78  OptionsParser args(argc, argv);
79  args.AddOption(&mesh_file, "-m", "--mesh",
80  "Mesh file to use.");
81  args.AddOption(&order, "-o", "--order",
82  "Finite element order (polynomial degree) or -1 for"
83  " isoparametric space.");
84  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
85  "--no-static-condensation", "Enable static condensation.");
86  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
87  "--no-partial-assembly", "Enable Partial Assembly.");
88  args.AddOption(&device_config, "-d", "--device",
89  "Device configuration string, see Device::Configure().");
90  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
91  "--no-visualization",
92  "Enable or disable GLVis visualization.");
93  args.Parse();
94  if (!args.Good())
95  {
96  if (myid == 0)
97  {
98  args.PrintUsage(cout);
99  }
100  MPI_Finalize();
101  return 1;
102  }
103  if (myid == 0)
104  {
105  args.PrintOptions(cout);
106  }
107 
108  // 3. Enable hardware devices such as GPUs, and programming models such as
109  // CUDA, OCCA, RAJA and OpenMP based on command line options.
110  Device device(device_config);
111  if (myid == 0) { device.Print(); }
112 
113  // 4. Read the (serial) mesh from the given mesh file on all processors. We
114  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
115  // and volume meshes with the same code.
116  Mesh mesh(mesh_file, 1, 1);
117  int dim = mesh.Dimension();
118 
119  // 5. Refine the serial mesh on all processors to increase the resolution. In
120  // this example we do 'ref_levels' of uniform refinement. We choose
121  // 'ref_levels' to be the largest number that gives a final mesh with no
122  // more than 10,000 elements.
123  {
124  int ref_levels =
125  (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
126  for (int l = 0; l < ref_levels; l++)
127  {
128  mesh.UniformRefinement();
129  }
130  }
131 
132  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
133  // this mesh further in parallel to increase the resolution. Once the
134  // parallel mesh is defined, the serial mesh can be deleted.
135  ParMesh pmesh(MPI_COMM_WORLD, mesh);
136  mesh.Clear();
137  {
138  int par_ref_levels = 2;
139  for (int l = 0; l < par_ref_levels; l++)
140  {
141  pmesh.UniformRefinement();
142  }
143  }
144 
145  // 7. Define a parallel finite element space on the parallel mesh. Here we
146  // use continuous Lagrange finite elements of the specified order. If
147  // order < 1, we instead use an isoparametric/isogeometric space.
149  bool delete_fec;
150  if (order > 0)
151  {
152  fec = new H1_FECollection(order, dim);
153  delete_fec = true;
154  }
155  else if (pmesh.GetNodes())
156  {
157  fec = pmesh.GetNodes()->OwnFEC();
158  delete_fec = false;
159  if (myid == 0)
160  {
161  cout << "Using isoparametric FEs: " << fec->Name() << endl;
162  }
163  }
164  else
165  {
166  fec = new H1_FECollection(order = 1, dim);
167  delete_fec = true;
168  }
169  ParFiniteElementSpace fespace(&pmesh, fec);
170  HYPRE_Int size = fespace.GlobalTrueVSize();
171  if (myid == 0)
172  {
173  cout << "Number of finite element unknowns: " << size << endl;
174  }
175 
176  // 8. Determine the list of true (i.e. parallel conforming) essential
177  // boundary dofs. In this example, the boundary conditions are defined
178  // by marking all the boundary attributes from the mesh as essential
179  // (Dirichlet) and converting them to a list of true dofs.
180  Array<int> ess_tdof_list;
181  if (pmesh.bdr_attributes.Size())
182  {
183  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
184  ess_bdr = 1;
185  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
186  }
187 
188  // 9. Set up the parallel linear form b(.) which corresponds to the
189  // right-hand side of the FEM linear system, which in this case is
190  // (1,phi_i) where phi_i are the basis functions in fespace.
191  ParLinearForm b(&fespace);
192  ConstantCoefficient one(1.0);
194  b.Assemble();
195 
196  // 10. Define the solution vector x as a parallel finite element grid function
197  // corresponding to fespace. Initialize x with initial guess of zero,
198  // which satisfies the boundary conditions.
199  ParGridFunction x(&fespace);
200  x = 0.0;
201 
202  // 11. Set up the parallel bilinear form a(.,.) on the finite element space
203  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
204  // domain integrator.
205  ParBilinearForm a(&fespace);
206  if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
208 
209  // 12. Assemble the parallel bilinear form and the corresponding linear
210  // system, applying any necessary transformations such as: parallel
211  // assembly, eliminating boundary conditions, applying conforming
212  // constraints for non-conforming AMR, static condensation, etc.
213  if (static_cond) { a.EnableStaticCondensation(); }
214  a.Assemble();
215 
216  OperatorPtr A;
217  Vector B, X;
218  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
219 
220  // 13. Solve the linear system A X = B.
221  // * With full assembly, use the BoomerAMG preconditioner from hypre.
222  // * With partial assembly, use Jacobi smoothing, for now.
223  Solver *prec = NULL;
224  if (pa)
225  {
226  if (UsesTensorBasis(fespace))
227  {
228  prec = new OperatorJacobiSmoother(a, ess_tdof_list);
229  }
230  }
231  else
232  {
233  prec = new HypreBoomerAMG;
234  }
235  CGSolver cg(MPI_COMM_WORLD);
236  cg.SetRelTol(1e-12);
237  cg.SetMaxIter(2000);
238  cg.SetPrintLevel(1);
239  if (prec) { cg.SetPreconditioner(*prec); }
240  cg.SetOperator(*A);
241  cg.Mult(B, X);
242  delete prec;
243 
244  // 14. Recover the parallel grid function corresponding to X. This is the
245  // local finite element solution on each processor.
246  a.RecoverFEMSolution(X, b, x);
247 
248  // 15. Save the refined mesh and the solution in parallel. This output can
249  // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
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 solution by socket to a GLVis server.
265  if (visualization)
266  {
267  char vishost[] = "localhost";
268  int visport = 19916;
269  socketstream sol_sock(vishost, visport);
270  sol_sock << "parallel " << num_procs << " " << myid << "\n";
271  sol_sock.precision(8);
272  sol_sock << "solution\n" << pmesh << x << flush;
273  }
274 
275  // 17. Free the used memory.
276  if (delete_fec)
277  {
278  delete fec;
279  }
280  MPI_Finalize();
281 
282  return 0;
283 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:124
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:93
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:775
Conjugate gradient method.
Definition: solvers.hpp:258
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:535
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:79
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:737
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:819
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:261
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
bool UsesTensorBasis(const FiniteElementSpace &fes)
Definition: fespace.hpp:983
int main(int argc, char *argv[])
Definition: ex1.cpp:66
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1079
Class for parallel linear form.
Definition: plinearform.hpp:26
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:70
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[]
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:128
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:8382
constexpr int visport
void SetMaxIter(int max_it)
Definition: solvers.hpp:98
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.
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4169
int Dimension() const
Definition: mesh.hpp:788
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:434
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:201
void SetRelTol(double rtol)
Definition: solvers.hpp:96
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
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
double a
Definition: lissajous.cpp:41
virtual const char * Name() const
Definition: fe_coll.hpp:61
int dim
Definition: ex24.cpp:53
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 Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:719
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:272
Vector data type.
Definition: vector.hpp:51
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6603
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:91
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
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:118
Class for parallel meshes.
Definition: pmesh.hpp:32
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:145