MFEM  v4.3.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-cpu -o 4 -a
36 // * mpirun -np 4 ex1p -pa -d ceed-cuda
37 // * mpirun -np 4 ex1p -pa -d ceed-hip
38 // mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
39 // mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
40 //
41 // Description: This example code demonstrates the use of MFEM to define a
42 // simple finite element discretization of the Laplace problem
43 // -Delta u = 1 with homogeneous Dirichlet boundary conditions.
44 // Specifically, we discretize using a FE space of the specified
45 // order, or if order < 1 using an isoparametric/isogeometric
46 // space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
47 // NURBS mesh, etc.)
48 //
49 // The example highlights the use of mesh refinement, finite
50 // element grid functions, as well as linear and bilinear forms
51 // corresponding to the left-hand side and right-hand side of the
52 // discrete linear system. We also cover the explicit elimination
53 // of essential boundary conditions, static condensation, and the
54 // optional connection to the GLVis tool for visualization.
55 
56 #include "mfem.hpp"
57 #include <fstream>
58 #include <iostream>
59 
60 using namespace std;
61 using namespace mfem;
62 
63 int main(int argc, char *argv[])
64 {
65  // 1. Initialize MPI.
66  MPI_Session mpi;
67  int num_procs = mpi.WorldSize();
68  int myid = mpi.WorldRank();
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  bool algebraic_ceed = false;
78 
79  OptionsParser args(argc, argv);
80  args.AddOption(&mesh_file, "-m", "--mesh",
81  "Mesh file to use.");
82  args.AddOption(&order, "-o", "--order",
83  "Finite element order (polynomial degree) or -1 for"
84  " isoparametric space.");
85  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
86  "--no-static-condensation", "Enable static condensation.");
87  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
88  "--no-partial-assembly", "Enable Partial Assembly.");
89  args.AddOption(&device_config, "-d", "--device",
90  "Device configuration string, see Device::Configure().");
91 #ifdef MFEM_USE_CEED
92  args.AddOption(&algebraic_ceed, "-a", "--algebraic",
93  "-no-a", "--no-algebraic",
94  "Use algebraic Ceed solver");
95 #endif
96  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
97  "--no-visualization",
98  "Enable or disable GLVis visualization.");
99  args.Parse();
100  if (!args.Good())
101  {
102  if (myid == 0)
103  {
104  args.PrintUsage(cout);
105  }
106  return 1;
107  }
108  if (myid == 0)
109  {
110  args.PrintOptions(cout);
111  }
112 
113  // 3. Enable hardware devices such as GPUs, and programming models such as
114  // CUDA, OCCA, RAJA and OpenMP based on command line options.
115  Device device(device_config);
116  if (myid == 0) { device.Print(); }
117 
118  // 4. Read the (serial) mesh from the given mesh file on all processors. We
119  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
120  // and volume meshes with the same code.
121  Mesh mesh(mesh_file, 1, 1);
122  int dim = mesh.Dimension();
123 
124  // 5. Refine the serial mesh on all processors to increase the resolution. In
125  // this example we do 'ref_levels' of uniform refinement. We choose
126  // 'ref_levels' to be the largest number that gives a final mesh with no
127  // more than 10,000 elements.
128  {
129  int ref_levels =
130  (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
131  for (int l = 0; l < ref_levels; l++)
132  {
133  mesh.UniformRefinement();
134  }
135  }
136 
137  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
138  // this mesh further in parallel to increase the resolution. Once the
139  // parallel mesh is defined, the serial mesh can be deleted.
140  ParMesh pmesh(MPI_COMM_WORLD, mesh);
141  mesh.Clear();
142  {
143  int par_ref_levels = 2;
144  for (int l = 0; l < par_ref_levels; l++)
145  {
146  pmesh.UniformRefinement();
147  }
148  }
149 
150  // 7. Define a parallel finite element space on the parallel mesh. Here we
151  // use continuous Lagrange finite elements of the specified order. If
152  // order < 1, we instead use an isoparametric/isogeometric space.
154  bool delete_fec;
155  if (order > 0)
156  {
157  fec = new H1_FECollection(order, dim);
158  delete_fec = true;
159  }
160  else if (pmesh.GetNodes())
161  {
162  fec = pmesh.GetNodes()->OwnFEC();
163  delete_fec = false;
164  if (myid == 0)
165  {
166  cout << "Using isoparametric FEs: " << fec->Name() << endl;
167  }
168  }
169  else
170  {
171  fec = new H1_FECollection(order = 1, dim);
172  delete_fec = true;
173  }
174  ParFiniteElementSpace fespace(&pmesh, fec);
175  HYPRE_BigInt size = fespace.GlobalTrueVSize();
176  if (myid == 0)
177  {
178  cout << "Number of finite element unknowns: " << size << endl;
179  }
180 
181  // 8. Determine the list of true (i.e. parallel conforming) essential
182  // boundary dofs. In this example, the boundary conditions are defined
183  // by marking all the boundary attributes from the mesh as essential
184  // (Dirichlet) and converting them to a list of true dofs.
185  Array<int> ess_tdof_list;
186  if (pmesh.bdr_attributes.Size())
187  {
188  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
189  ess_bdr = 1;
190  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
191  }
192 
193  // 9. Set up the parallel linear form b(.) which corresponds to the
194  // right-hand side of the FEM linear system, which in this case is
195  // (1,phi_i) where phi_i are the basis functions in fespace.
196  ParLinearForm b(&fespace);
197  ConstantCoefficient one(1.0);
199  b.Assemble();
200 
201  // 10. Define the solution vector x as a parallel finite element grid
202  // function corresponding to fespace. Initialize x with initial guess of
203  // zero, which satisfies the boundary conditions.
204  ParGridFunction x(&fespace);
205  x = 0.0;
206 
207  // 11. Set up the parallel bilinear form a(.,.) on the finite element space
208  // corresponding to the Laplacian operator -Delta, by adding the
209  // Diffusion domain integrator.
210  ParBilinearForm a(&fespace);
211  if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
213 
214  // 12. Assemble the parallel bilinear form and the corresponding linear
215  // system, applying any necessary transformations such as: parallel
216  // assembly, eliminating boundary conditions, applying conforming
217  // constraints for non-conforming AMR, static condensation, etc.
218  if (static_cond) { a.EnableStaticCondensation(); }
219  a.Assemble();
220 
221  OperatorPtr A;
222  Vector B, X;
223  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
224 
225  // 13. Solve the linear system A X = B.
226  // * With full assembly, use the BoomerAMG preconditioner from hypre.
227  // * With partial assembly, use Jacobi smoothing, for now.
228  Solver *prec = NULL;
229  if (pa)
230  {
231  if (UsesTensorBasis(fespace))
232  {
233  if (algebraic_ceed)
234  {
235  prec = new ceed::AlgebraicSolver(a, ess_tdof_list);
236  }
237  else
238  {
239  prec = new OperatorJacobiSmoother(a, ess_tdof_list);
240  }
241  }
242  }
243  else
244  {
245  prec = new HypreBoomerAMG;
246  }
247  CGSolver cg(MPI_COMM_WORLD);
248  cg.SetRelTol(1e-12);
249  cg.SetMaxIter(2000);
250  cg.SetPrintLevel(1);
251  if (prec) { cg.SetPreconditioner(*prec); }
252  cg.SetOperator(*A);
253  cg.Mult(B, X);
254  delete prec;
255 
256  // 14. Recover the parallel grid function corresponding to X. This is the
257  // local finite element solution on each processor.
258  a.RecoverFEMSolution(X, b, x);
259 
260  // 15. Save the refined mesh and the solution in parallel. This output can
261  // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
262  {
263  ostringstream mesh_name, sol_name;
264  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
265  sol_name << "sol." << setfill('0') << setw(6) << myid;
266 
267  ofstream mesh_ofs(mesh_name.str().c_str());
268  mesh_ofs.precision(8);
269  pmesh.Print(mesh_ofs);
270 
271  ofstream sol_ofs(sol_name.str().c_str());
272  sol_ofs.precision(8);
273  x.Save(sol_ofs);
274  }
275 
276  // 16. Send the solution by socket to a GLVis server.
277  if (visualization)
278  {
279  char vishost[] = "localhost";
280  int visport = 19916;
281  socketstream sol_sock(vishost, visport);
282  sol_sock << "parallel " << num_procs << " " << myid << "\n";
283  sol_sock.precision(8);
284  sol_sock << "solution\n" << pmesh << x << flush;
285  }
286 
287  // 17. Free the used memory.
288  if (delete_fec)
289  {
290  delete fec;
291  }
292 
293  return 0;
294 }
int WorldSize() const
Return MPI_COMM_WORLD&#39;s size.
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:790
Conjugate gradient method.
Definition: solvers.hpp:316
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:616
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:275
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:841
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
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:957
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1387
Class for parallel linear form.
Definition: plinearform.hpp:26
A simple convenience class that calls MPI_Init() at construction and MPI_Finalize() at destruction...
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:71
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:130
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
constexpr int visport
void SetMaxIter(int max_it)
Definition: solvers.hpp:100
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4382
int Dimension() const
Definition: mesh.hpp:911
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
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:204
void SetRelTol(double rtol)
Definition: solvers.hpp:98
int WorldRank() const
Return MPI_COMM_WORLD&#39;s rank.
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
HYPRE_Int HYPRE_BigInt
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:327
Class for parallel bilinear form.
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:828
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:330
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7343
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:92
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
Base class for solvers.
Definition: operator.hpp:648
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:121
Class for parallel meshes.
Definition: pmesh.hpp:32
int main()
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:150