MFEM  v4.5.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 and HYPRE.
38  Mpi::Init(argc, argv);
39  int num_procs = Mpi::WorldSize();
40  int myid = Mpi::WorldRank();
41  Hypre::Init();
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);
53  args.AddOption(&elem_type, "-e", "--elem",
54  "Type of elements to use: 0 - triangles, 1 - quads.");
55  args.AddOption(&order, "-o", "--order",
56  "Finite element order (polynomial degree).");
57  args.AddOption(&ref_levels, "-r", "--refine",
58  "Number of times to refine the mesh uniformly.");
59  args.AddOption(&amr, "-amr", "--refine-locally",
60  "Additional local (non-conforming) refinement:"
61  " 1 = refine around north pole, 2 = refine randomly.");
62  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
63  "--no-visualization",
64  "Enable or disable GLVis visualization.");
65  args.AddOption(&always_snap, "-snap", "--always-snap", "-no-snap",
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");
69  args.AddOption(&device_config, "-d", "--device",
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  return 1;
79  }
80  if (myid == 0)
81  {
82  args.PrintOptions(cout);
83  }
84 
85  // 3. Enable hardware devices such as GPUs, and programming models such as
86  // CUDA, OCCA, RAJA and OpenMP based on command line options.
87  Device device(device_config);
88  if (myid == 0) { device.Print(); }
89 
90  // 4. Generate an initial high-order (surface) mesh on the unit sphere. The
91  // Mesh object represents a 2D mesh in 3 spatial dimensions. We first add
92  // the elements and the vertices of the mesh, and then make it high-order
93  // by specifying a finite element space for its nodes.
94  int Nvert = 8, Nelem = 6;
95  if (elem_type == 0)
96  {
97  Nvert = 6;
98  Nelem = 8;
99  }
100  Mesh *mesh = new Mesh(2, Nvert, Nelem, 0, 3);
101 
102  if (elem_type == 0) // inscribed octahedron
103  {
104  const double tri_v[6][3] =
105  {
106  { 1, 0, 0}, { 0, 1, 0}, {-1, 0, 0},
107  { 0, -1, 0}, { 0, 0, 1}, { 0, 0, -1}
108  };
109  const int tri_e[8][3] =
110  {
111  {0, 1, 4}, {1, 2, 4}, {2, 3, 4}, {3, 0, 4},
112  {1, 0, 5}, {2, 1, 5}, {3, 2, 5}, {0, 3, 5}
113  };
114 
115  for (int j = 0; j < Nvert; j++)
116  {
117  mesh->AddVertex(tri_v[j]);
118  }
119  for (int j = 0; j < Nelem; j++)
120  {
121  int attribute = j + 1;
122  mesh->AddTriangle(tri_e[j], attribute);
123  }
124  mesh->FinalizeTriMesh(1, 1, true);
125  }
126  else // inscribed cube
127  {
128  const double quad_v[8][3] =
129  {
130  {-1, -1, -1}, {+1, -1, -1}, {+1, +1, -1}, {-1, +1, -1},
131  {-1, -1, +1}, {+1, -1, +1}, {+1, +1, +1}, {-1, +1, +1}
132  };
133  const int quad_e[6][4] =
134  {
135  {3, 2, 1, 0}, {0, 1, 5, 4}, {1, 2, 6, 5},
136  {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}
137  };
138 
139  for (int j = 0; j < Nvert; j++)
140  {
141  mesh->AddVertex(quad_v[j]);
142  }
143  for (int j = 0; j < Nelem; j++)
144  {
145  int attribute = j + 1;
146  mesh->AddQuad(quad_e[j], attribute);
147  }
148  mesh->FinalizeQuadMesh(1, 1, true);
149  }
150 
151  // Set the space for the high-order mesh nodes.
152  H1_FECollection fec(order, mesh->Dimension());
153  FiniteElementSpace nodal_fes(mesh, &fec, mesh->SpaceDimension());
154  mesh->SetNodalFESpace(&nodal_fes);
155 
156  // 5. Refine the mesh while snapping nodes to the sphere. Number of parallel
157  // refinements is fixed to 2.
158  for (int l = 0; l <= ref_levels; l++)
159  {
160  if (l > 0) // for l == 0 just perform snapping
161  {
162  mesh->UniformRefinement();
163  }
164 
165  // Snap the nodes of the refined mesh back to sphere surface.
166  if (always_snap)
167  {
168  SnapNodes(*mesh);
169  }
170  }
171 
172  if (amr == 1)
173  {
174  Vertex target(0.0, 0.0, 1.0);
175  for (int l = 0; l < 3; l++)
176  {
177  mesh->RefineAtVertex(target);
178  }
179  SnapNodes(*mesh);
180  }
181  else if (amr == 2)
182  {
183  for (int l = 0; l < 2; l++)
184  {
185  mesh->RandomRefinement(0.5); // 50% probability
186  }
187  SnapNodes(*mesh);
188  }
189 
190  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
191  delete mesh;
192  {
193  int par_ref_levels = 2;
194  for (int l = 0; l < par_ref_levels; l++)
195  {
196  pmesh->UniformRefinement();
197 
198  // Snap the nodes of the refined mesh back to sphere surface.
199  if (always_snap)
200  {
201  SnapNodes(*pmesh);
202  }
203  }
204  if (!always_snap || par_ref_levels < 1)
205  {
206  SnapNodes(*pmesh);
207  }
208  }
209 
210  if (amr == 1)
211  {
212  Vertex target(0.0, 0.0, 1.0);
213  for (int l = 0; l < 2; l++)
214  {
215  pmesh->RefineAtVertex(target);
216  }
217  SnapNodes(*pmesh);
218  }
219  else if (amr == 2)
220  {
221  for (int l = 0; l < 2; l++)
222  {
223  pmesh->RandomRefinement(0.5); // 50% probability
224  }
225  SnapNodes(*pmesh);
226  }
227 
228  // 6. Define a finite element space on the mesh. Here we use isoparametric
229  // finite elements -- the same as the mesh nodes.
230  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, &fec);
231  HYPRE_BigInt size = fespace->GlobalTrueVSize();
232  if (myid == 0)
233  {
234  cout << "Number of unknowns: " << size << endl;
235  }
236 
237  // 7. Set up the linear form b(.) which corresponds to the right-hand side of
238  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
239  // the basis functions in the finite element fespace.
240  ParLinearForm *b = new ParLinearForm(fespace);
241  ConstantCoefficient one(1.0);
244  b->AddDomainIntegrator(new DomainLFIntegrator(rhs_coef));
245  b->Assemble();
246 
247  // 8. Define the solution vector x as a finite element grid function
248  // corresponding to fespace. Initialize x with initial guess of zero.
249  ParGridFunction x(fespace);
250  x = 0.0;
251 
252  // 9. Set up the bilinear form a(.,.) on the finite element space
253  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
254  // and Mass domain integrators.
255  ParBilinearForm *a = new ParBilinearForm(fespace);
256  a->AddDomainIntegrator(new DiffusionIntegrator(one));
257  a->AddDomainIntegrator(new MassIntegrator(one));
258 
259  // 10. Assemble the parallel linear system, applying any transformations
260  // such as: parallel assembly, applying conforming constraints, etc.
261  a->Assemble();
262  HypreParMatrix A;
263  Vector B, X;
264  Array<int> empty_tdof_list;
265  a->FormLinearSystem(empty_tdof_list, x, *b, A, X, B);
266 
267  // 11. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
268  // preconditioner from hypre. Extract the parallel grid function x
269  // corresponding to the finite element approximation X. This is the local
270  // solution on each processor.
271  HypreSolver *amg = new HypreBoomerAMG(A);
272  HyprePCG *pcg = new HyprePCG(A);
273  pcg->SetTol(1e-12);
274  pcg->SetMaxIter(200);
275  pcg->SetPrintLevel(2);
276  pcg->SetPreconditioner(*amg);
277  pcg->Mult(B, X);
278  a->RecoverFEMSolution(X, *b, x);
279 
280  delete a;
281  delete b;
282 
283  // 12. Compute and print the L^2 norm of the error.
284  double error = x.ComputeL2Error(sol_coef);
285  if (myid == 0)
286  {
287  cout << "\nL2 norm of error: " << error << endl;
288  }
289 
290  // 13. Save the refined mesh and the solution. This output can be viewed
291  // later using GLVis: "glvis -np <np> -m sphere_refined -g sol".
292  {
293  ostringstream mesh_name, sol_name;
294  mesh_name << "sphere_refined." << setfill('0') << setw(6) << myid;
295  sol_name << "sol." << setfill('0') << setw(6) << myid;
296 
297  ofstream mesh_ofs(mesh_name.str().c_str());
298  mesh_ofs.precision(8);
299  pmesh->Print(mesh_ofs);
300 
301  ofstream sol_ofs(sol_name.str().c_str());
302  sol_ofs.precision(8);
303  x.Save(sol_ofs);
304  }
305 
306  // 14. Send the solution by socket to a GLVis server.
307  if (visualization)
308  {
309  char vishost[] = "localhost";
310  int visport = 19916;
311  socketstream sol_sock(vishost, visport);
312  sol_sock << "parallel " << num_procs << " " << myid << "\n";
313  sol_sock.precision(8);
314  sol_sock << "solution\n" << *pmesh << x << flush;
315  }
316 
317  // 15. Free the used memory.
318  delete pcg;
319  delete amg;
320  delete fespace;
321  delete pmesh;
322 
323  return 0;
324 }
325 
326 double analytic_solution(const Vector &x)
327 {
328  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
329  return x(0)*x(1)/l2;
330 }
331 
332 double analytic_rhs(const Vector &x)
333 {
334  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
335  return 7*x(0)*x(1)/l2;
336 }
337 
338 void SnapNodes(Mesh &mesh)
339 {
340  GridFunction &nodes = *mesh.GetNodes();
341  Vector node(mesh.SpaceDimension());
342  for (int i = 0; i < nodes.FESpace()->GetNDofs(); i++)
343  {
344  for (int d = 0; d < mesh.SpaceDimension(); d++)
345  {
346  node(d) = nodes(nodes.FESpace()->DofToVDof(i, d));
347  }
348 
349  node /= node.Norml2();
350 
351  for (int d = 0; d < mesh.SpaceDimension(); d++)
352  {
353  nodes(nodes.FESpace()->DofToVDof(i, d)) = node(d);
354  }
355  }
356  if (mesh.Nonconforming())
357  {
358  // Snap hanging nodes to the master side.
359  Vector tnodes;
360  nodes.GetTrueDofs(tnodes);
361  nodes.SetFromTrueDofs(tnodes);
362  }
363 }
void SetTol(double tol)
Definition: hypre.cpp:3995
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Definition: mesh.cpp:1674
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
int Dimension() const
Definition: mesh.hpp:1006
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:4043
int AddTriangle(int v1, int v2, int v3, int attr=1)
Definition: mesh.cpp:1660
bool Nonconforming() const
Definition: mesh.hpp:1651
void SnapNodes(Mesh &mesh)
Definition: ex7p.cpp:338
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:584
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
Abstract parallel finite element space.
Definition: pfespace.hpp:28
Data type for vertex.
Definition: vertex.hpp:22
STL namespace.
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:4015
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1579
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction.
Definition: gridfunc.cpp:367
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1613
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:1939
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
constexpr char vishost[]
double analytic_solution(const Vector &x)
Definition: ex7p.cpp:326
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9802
constexpr int visport
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:9483
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:5285
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:668
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: pgridfunc.hpp:281
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
void SetMaxIter(int max_iter)
Definition: hypre.cpp:4005
PCG solver in hypre.
Definition: hypre.hpp:1204
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
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:9502
int SpaceDimension() const
Definition: mesh.hpp:1007
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:873
double analytic_rhs(const Vector &x)
Definition: ex7p.cpp:332
double a
Definition: lissajous.cpp:41
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition: mesh.cpp:1968
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:4020
Class for parallel bilinear form.
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:812
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:1092
A general function coefficient.
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:251
Vector data type.
Definition: vector.hpp:60
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7894
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4770
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:343
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: gridfunc.cpp:382
Class for parallel meshes.
Definition: pmesh.hpp:32
int main(int argc, char *argv[])
Definition: ex7p.cpp:35