MFEM  v4.3.0 Finite element discretization library
ex7p.cpp
Go to the documentation of this file.
1 // MFEM Example 7 - Parallel Version
2 //
3 // Compile with: make ex7p
4 //
5 // Sample runs: mpirun -np 4 ex7p -e 0 -o 2 -r 4
6 // mpirun -np 4 ex7p -e 1 -o 2 -r 4 -snap
7 // mpirun -np 4 ex7p -e 0 -amr 1
8 // mpirun -np 4 ex7p -e 1 -amr 2 -o 2
9 //
10 // Description: This example code demonstrates the use of MFEM to define a
11 // triangulation of a unit sphere and a simple isoparametric
12 // finite element discretization of the Laplace problem with mass
13 // term, -Delta u + u = f.
14 //
15 // The example highlights mesh generation, the use of mesh
16 // refinement, high-order meshes and finite elements, as well as
17 // surface-based linear and bilinear forms corresponding to the
18 // left-hand side and right-hand side of the discrete linear
19 // system. Simple local mesh refinement is also demonstrated.
20 //
21 // We recommend viewing Example 1 before viewing this example.
22
23 #include "mfem.hpp"
24 #include <fstream>
25 #include <iostream>
26
27 using namespace std;
28 using namespace mfem;
29
30 // Exact solution and r.h.s., see below for implementation.
31 double analytic_solution(const Vector &x);
32 double analytic_rhs(const Vector &x);
33 void SnapNodes(Mesh &mesh);
34
35 int main(int argc, char *argv[])
36 {
37  // 1. Initialize MPI.
38  int num_procs, myid;
39  MPI_Init(&argc, &argv);
40  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
41  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
42
43  // 2. Parse command-line options.
44  int elem_type = 1;
45  int ref_levels = 2;
46  int amr = 0;
47  int order = 2;
48  bool always_snap = false;
49  bool visualization = 1;
50  const char *device_config = "cpu";
51
52  OptionsParser args(argc, argv);
54  "Type of elements to use: 0 - triangles, 1 - quads.");
56  "Finite element order (polynomial degree).");
58  "Number of times to refine the mesh uniformly.");
61  " 1 = refine around north pole, 2 = refine randomly.");
63  "--no-visualization",
64  "Enable or disable GLVis visualization.");
66  "--snap-at-the-end",
67  "If true, snap nodes to the sphere initially and after each refinement "
68  "otherwise, snap only after the last refinement");
70  "Device configuration string, see Device::Configure().");
71  args.Parse();
72  if (!args.Good())
73  {
74  if (myid == 0)
75  {
76  args.PrintUsage(cout);
77  }
78  MPI_Finalize();
79  return 1;
80  }
81  if (myid == 0)
82  {
83  args.PrintOptions(cout);
84  }
85
86  // 3. Enable hardware devices such as GPUs, and programming models such as
87  // CUDA, OCCA, RAJA and OpenMP based on command line options.
88  Device device(device_config);
89  if (myid == 0) { device.Print(); }
90
91  // 4. Generate an initial high-order (surface) mesh on the unit sphere. The
92  // Mesh object represents a 2D mesh in 3 spatial dimensions. We first add
93  // the elements and the vertices of the mesh, and then make it high-order
94  // by specifying a finite element space for its nodes.
95  int Nvert = 8, Nelem = 6;
96  if (elem_type == 0)
97  {
98  Nvert = 6;
99  Nelem = 8;
100  }
101  Mesh *mesh = new Mesh(2, Nvert, Nelem, 0, 3);
102
103  if (elem_type == 0) // inscribed octahedron
104  {
105  const double tri_v[6][3] =
106  {
107  { 1, 0, 0}, { 0, 1, 0}, {-1, 0, 0},
108  { 0, -1, 0}, { 0, 0, 1}, { 0, 0, -1}
109  };
110  const int tri_e[8][3] =
111  {
112  {0, 1, 4}, {1, 2, 4}, {2, 3, 4}, {3, 0, 4},
113  {1, 0, 5}, {2, 1, 5}, {3, 2, 5}, {0, 3, 5}
114  };
115
116  for (int j = 0; j < Nvert; j++)
117  {
119  }
120  for (int j = 0; j < Nelem; j++)
121  {
122  int attribute = j + 1;
124  }
125  mesh->FinalizeTriMesh(1, 1, true);
126  }
127  else // inscribed cube
128  {
130  {
131  {-1, -1, -1}, {+1, -1, -1}, {+1, +1, -1}, {-1, +1, -1},
132  {-1, -1, +1}, {+1, -1, +1}, {+1, +1, +1}, {-1, +1, +1}
133  };
135  {
136  {3, 2, 1, 0}, {0, 1, 5, 4}, {1, 2, 6, 5},
137  {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}
138  };
139
140  for (int j = 0; j < Nvert; j++)
141  {
143  }
144  for (int j = 0; j < Nelem; j++)
145  {
146  int attribute = j + 1;
148  }
150  }
151
152  // Set the space for the high-order mesh nodes.
153  H1_FECollection fec(order, mesh->Dimension());
154  FiniteElementSpace nodal_fes(mesh, &fec, mesh->SpaceDimension());
155  mesh->SetNodalFESpace(&nodal_fes);
156
157  // 5. Refine the mesh while snapping nodes to the sphere. Number of parallel
158  // refinements is fixed to 2.
159  for (int l = 0; l <= ref_levels; l++)
160  {
161  if (l > 0) // for l == 0 just perform snapping
162  {
163  mesh->UniformRefinement();
164  }
165
166  // Snap the nodes of the refined mesh back to sphere surface.
167  if (always_snap)
168  {
169  SnapNodes(*mesh);
170  }
171  }
172
173  if (amr == 1)
174  {
175  Vertex target(0.0, 0.0, 1.0);
176  for (int l = 0; l < 3; l++)
177  {
178  mesh->RefineAtVertex(target);
179  }
180  SnapNodes(*mesh);
181  }
182  else if (amr == 2)
183  {
184  for (int l = 0; l < 2; l++)
185  {
186  mesh->RandomRefinement(0.5); // 50% probability
187  }
188  SnapNodes(*mesh);
189  }
190
191  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
192  delete mesh;
193  {
194  int par_ref_levels = 2;
195  for (int l = 0; l < par_ref_levels; l++)
196  {
197  pmesh->UniformRefinement();
198
199  // Snap the nodes of the refined mesh back to sphere surface.
200  if (always_snap)
201  {
202  SnapNodes(*pmesh);
203  }
204  }
205  if (!always_snap || par_ref_levels < 1)
206  {
207  SnapNodes(*pmesh);
208  }
209  }
210
211  if (amr == 1)
212  {
213  Vertex target(0.0, 0.0, 1.0);
214  for (int l = 0; l < 2; l++)
215  {
216  pmesh->RefineAtVertex(target);
217  }
218  SnapNodes(*pmesh);
219  }
220  else if (amr == 2)
221  {
222  for (int l = 0; l < 2; l++)
223  {
224  pmesh->RandomRefinement(0.5); // 50% probability
225  }
226  SnapNodes(*pmesh);
227  }
228
229  // 6. Define a finite element space on the mesh. Here we use isoparametric
230  // finite elements -- the same as the mesh nodes.
231  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, &fec);
232  HYPRE_BigInt size = fespace->GlobalTrueVSize();
233  if (myid == 0)
234  {
235  cout << "Number of unknowns: " << size << endl;
236  }
237
238  // 7. Set up the linear form b(.) which corresponds to the right-hand side of
239  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
240  // the basis functions in the finite element fespace.
241  ParLinearForm *b = new ParLinearForm(fespace);
242  ConstantCoefficient one(1.0);
246  b->Assemble();
247
248  // 8. Define the solution vector x as a finite element grid function
249  // corresponding to fespace. Initialize x with initial guess of zero.
250  ParGridFunction x(fespace);
251  x = 0.0;
252
253  // 9. Set up the bilinear form a(.,.) on the finite element space
254  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
255  // and Mass domain integrators.
256  ParBilinearForm *a = new ParBilinearForm(fespace);
259
260  // 10. Assemble the parallel linear system, applying any transformations
261  // such as: parallel assembly, applying conforming constraints, etc.
262  a->Assemble();
263  HypreParMatrix A;
264  Vector B, X;
265  Array<int> empty_tdof_list;
266  a->FormLinearSystem(empty_tdof_list, x, *b, A, X, B);
267
268  // 11. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
269  // preconditioner from hypre. Extract the parallel grid function x
270  // corresponding to the finite element approximation X. This is the local
271  // solution on each processor.
272  HypreSolver *amg = new HypreBoomerAMG(A);
273  HyprePCG *pcg = new HyprePCG(A);
274  pcg->SetTol(1e-12);
275  pcg->SetMaxIter(200);
276  pcg->SetPrintLevel(2);
277  pcg->SetPreconditioner(*amg);
278  pcg->Mult(B, X);
279  a->RecoverFEMSolution(X, *b, x);
280
281  delete a;
282  delete b;
283
284  // 12. Compute and print the L^2 norm of the error.
285  double err = x.ComputeL2Error(sol_coef);
286  if (myid == 0)
287  {
288  cout << "\nL2 norm of error: " << err << endl;
289  }
290
291  // 13. Save the refined mesh and the solution. This output can be viewed
292  // later using GLVis: "glvis -np <np> -m sphere_refined -g sol".
293  {
294  ostringstream mesh_name, sol_name;
295  mesh_name << "sphere_refined." << setfill('0') << setw(6) << myid;
296  sol_name << "sol." << setfill('0') << setw(6) << myid;
297
298  ofstream mesh_ofs(mesh_name.str().c_str());
299  mesh_ofs.precision(8);
300  pmesh->Print(mesh_ofs);
301
302  ofstream sol_ofs(sol_name.str().c_str());
303  sol_ofs.precision(8);
304  x.Save(sol_ofs);
305  }
306
307  // 14. Send the solution by socket to a GLVis server.
308  if (visualization)
309  {
310  char vishost[] = "localhost";
311  int visport = 19916;
312  socketstream sol_sock(vishost, visport);
313  sol_sock << "parallel " << num_procs << " " << myid << "\n";
314  sol_sock.precision(8);
315  sol_sock << "solution\n" << *pmesh << x << flush;
316  }
317
318  // 15. Free the used memory.
319  delete pcg;
320  delete amg;
321  delete fespace;
322  delete pmesh;
323
324  MPI_Finalize();
325
326  return 0;
327 }
328
329 double analytic_solution(const Vector &x)
330 {
331  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
332  return x(0)*x(1)/l2;
333 }
334
335 double analytic_rhs(const Vector &x)
336 {
337  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
338  return 7*x(0)*x(1)/l2;
339 }
340
341 void SnapNodes(Mesh &mesh)
342 {
343  GridFunction &nodes = *mesh.GetNodes();
344  Vector node(mesh.SpaceDimension());
345  for (int i = 0; i < nodes.FESpace()->GetNDofs(); i++)
346  {
347  for (int d = 0; d < mesh.SpaceDimension(); d++)
348  {
349  node(d) = nodes(nodes.FESpace()->DofToVDof(i, d));
350  }
351
352  node /= node.Norml2();
353
354  for (int d = 0; d < mesh.SpaceDimension(); d++)
355  {
356  nodes(nodes.FESpace()->DofToVDof(i, d)) = node(d);
357  }
358  }
359  if (mesh.Nonconforming())
360  {
361  // Snap hanging nodes to the master side.
362  Vector tnodes;
363  nodes.GetTrueDofs(tnodes);
364  nodes.SetFromTrueDofs(tnodes);
365  }
366 }
void SetTol(double tol)
Definition: hypre.cpp:3569
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:537
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:233
Definition: mesh.cpp:1324
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(...
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:783
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:275
int AddTriangle(int v1, int v2, int v3, int attr=1)
Definition: mesh.cpp:1310
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
Abstract parallel finite element space.
Definition: pfespace.hpp:28
Data type for vertex.
Definition: vertex.hpp:22
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:3589
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1387
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1263
bool Nonconforming() const
Definition: mesh.hpp:1367
Class for parallel linear form.
Definition: plinearform.hpp:26
void FinalizeTriMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a triangular Mesh.
Definition: mesh.cpp:1556
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[]
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
constexpr int visport
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4382
double analytic_solution(const Vector &x)
Definition: ex7.cpp:252
void RandomRefinement(double prob, bool aniso=false, int nonconforming=-1, int nc_limit=0)
Refine each element with given probability. Uses GeneralRefinement.
Definition: mesh.cpp:8824
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction.
Definition: gridfunc.cpp:332
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:4843
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:629
int Dimension() const
Definition: mesh.hpp:911
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
void SetMaxIter(int max_iter)
Definition: hypre.cpp:3579
double analytic_rhs(const Vector &x)
Definition: ex7.cpp:258
int SpaceDimension() const
Definition: mesh.hpp:912
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
PCG solver in hypre.
Definition: hypre.hpp:1060
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
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
void RefineAtVertex(const Vertex &vert, double eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex. Uses GeneralRefinement.
Definition: mesh.cpp:8843
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
double a
Definition: lissajous.cpp:41
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:273
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition: mesh.cpp:1585
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:3594
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
void SnapNodes(Mesh &mesh)
Definition: ex7.cpp:264
Class for parallel bilinear form.
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:970
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7343
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
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
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:277
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: gridfunc.cpp:347
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:3617
Class for parallel meshes.
Definition: pmesh.hpp:32
int main()
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150