MFEM  v3.3
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 
51  OptionsParser args(argc, argv);
52  args.AddOption(&elem_type, "-e", "--elem",
53  "Type of elements to use: 0 - triangles, 1 - quads.");
54  args.AddOption(&order, "-o", "--order",
55  "Finite element order (polynomial degree).");
56  args.AddOption(&ref_levels, "-r", "--refine",
57  "Number of times to refine the mesh uniformly.");
58  args.AddOption(&amr, "-amr", "--refine-locally",
59  "Additional local (non-conforming) refinement:"
60  " 1 = refine around north pole, 2 = refine randomly.");
61  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
62  "--no-visualization",
63  "Enable or disable GLVis visualization.");
64  args.AddOption(&always_snap, "-snap", "--always-snap", "-no-snap",
65  "--snap-at-the-end",
66  "If true, snap nodes to the sphere initially and after each refinement "
67  "otherwise, snap only after the last refinement");
68  args.Parse();
69  if (!args.Good())
70  {
71  if (myid == 0)
72  {
73  args.PrintUsage(cout);
74  }
75  MPI_Finalize();
76  return 1;
77  }
78  if (myid == 0)
79  {
80  args.PrintOptions(cout);
81  }
82 
83  // 3. Generate an initial high-order (surface) mesh on the unit sphere. The
84  // Mesh object represents a 2D mesh in 3 spatial dimensions. We first add
85  // the elements and the vertices of the mesh, and then make it high-order
86  // by specifying a finite element space for its nodes.
87  int Nvert = 8, Nelem = 6;
88  if (elem_type == 0)
89  {
90  Nvert = 6;
91  Nelem = 8;
92  }
93  Mesh *mesh = new Mesh(2, Nvert, Nelem, 0, 3);
94 
95  if (elem_type == 0) // inscribed octahedron
96  {
97  const double tri_v[6][3] =
98  {
99  { 1, 0, 0}, { 0, 1, 0}, {-1, 0, 0},
100  { 0, -1, 0}, { 0, 0, 1}, { 0, 0, -1}
101  };
102  const int tri_e[8][3] =
103  {
104  {0, 1, 4}, {1, 2, 4}, {2, 3, 4}, {3, 0, 4},
105  {1, 0, 5}, {2, 1, 5}, {3, 2, 5}, {0, 3, 5}
106  };
107 
108  for (int j = 0; j < Nvert; j++)
109  {
110  mesh->AddVertex(tri_v[j]);
111  }
112  for (int j = 0; j < Nelem; j++)
113  {
114  int attribute = j + 1;
115  mesh->AddTriangle(tri_e[j], attribute);
116  }
117  mesh->FinalizeTriMesh(1, 1, true);
118  }
119  else // inscribed cube
120  {
121  const double quad_v[8][3] =
122  {
123  {-1, -1, -1}, {+1, -1, -1}, {+1, +1, -1}, {-1, +1, -1},
124  {-1, -1, +1}, {+1, -1, +1}, {+1, +1, +1}, {-1, +1, +1}
125  };
126  const int quad_e[6][4] =
127  {
128  {3, 2, 1, 0}, {0, 1, 5, 4}, {1, 2, 6, 5},
129  {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}
130  };
131 
132  for (int j = 0; j < Nvert; j++)
133  {
134  mesh->AddVertex(quad_v[j]);
135  }
136  for (int j = 0; j < Nelem; j++)
137  {
138  int attribute = j + 1;
139  mesh->AddQuad(quad_e[j], attribute);
140  }
141  mesh->FinalizeQuadMesh(1, 1, true);
142  }
143 
144  // Set the space for the high-order mesh nodes.
145  H1_FECollection fec(order, mesh->Dimension());
146  FiniteElementSpace nodal_fes(mesh, &fec, mesh->SpaceDimension());
147  mesh->SetNodalFESpace(&nodal_fes);
148 
149  // 4. Refine the mesh while snapping nodes to the sphere. Number of parallel
150  // refinements is fixed to 2.
151  for (int l = 0; l <= ref_levels; l++)
152  {
153  if (l > 0) // for l == 0 just perform snapping
154  {
155  mesh->UniformRefinement();
156  }
157 
158  // Snap the nodes of the refined mesh back to sphere surface.
159  if (always_snap)
160  {
161  SnapNodes(*mesh);
162  }
163  }
164 
165  if (amr == 1)
166  {
167  Vertex target(0.0, 0.0, 1.0);
168  for (int l = 0; l < 3; l++)
169  {
170  mesh->RefineAtVertex(target);
171  }
172  SnapNodes(*mesh);
173  }
174  else if (amr == 2)
175  {
176  for (int l = 0; l < 2; l++)
177  {
178  mesh->RandomRefinement(0.5); // 50% probability
179  }
180  SnapNodes(*mesh);
181  }
182 
183  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
184  delete mesh;
185  {
186  int par_ref_levels = 2;
187  for (int l = 0; l < par_ref_levels; l++)
188  {
189  pmesh->UniformRefinement();
190 
191  // Snap the nodes of the refined mesh back to sphere surface.
192  if (always_snap)
193  {
194  SnapNodes(*pmesh);
195  }
196  }
197  if (!always_snap || par_ref_levels < 1)
198  {
199  SnapNodes(*pmesh);
200  }
201  }
202 
203  if (amr == 1)
204  {
205  Vertex target(0.0, 0.0, 1.0);
206  for (int l = 0; l < 2; l++)
207  {
208  pmesh->RefineAtVertex(target);
209  }
210  SnapNodes(*pmesh);
211  }
212  else if (amr == 2)
213  {
214  for (int l = 0; l < 2; l++)
215  {
216  pmesh->RandomRefinement(0.5); // 50% probability
217  }
218  SnapNodes(*pmesh);
219  }
220 
221  // 5. Define a finite element space on the mesh. Here we use isoparametric
222  // finite elements -- the same as the mesh nodes.
223  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, &fec);
224  HYPRE_Int size = fespace->GlobalTrueVSize();
225  if (myid == 0)
226  {
227  cout << "Number of unknowns: " << size << endl;
228  }
229 
230  // 6. Set up the linear form b(.) which corresponds to the right-hand side of
231  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
232  // the basis functions in the finite element fespace.
233  ParLinearForm *b = new ParLinearForm(fespace);
234  ConstantCoefficient one(1.0);
237  b->AddDomainIntegrator(new DomainLFIntegrator(rhs_coef));
238  b->Assemble();
239 
240  // 7. Define the solution vector x as a finite element grid function
241  // corresponding to fespace. Initialize x with initial guess of zero.
242  ParGridFunction x(fespace);
243  x = 0.0;
244 
245  // 8. Set up the bilinear form a(.,.) on the finite element space
246  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
247  // and Mass domain integrators.
248  ParBilinearForm *a = new ParBilinearForm(fespace);
250  a->AddDomainIntegrator(new MassIntegrator(one));
251 
252  // 9. Assemble the parallel linear system, applying any transformations
253  // such as: parallel assembly, applying conforming constraints, etc.
254  a->Assemble();
255  HypreParMatrix A;
256  Vector B, X;
257  Array<int> empty_tdof_list;
258  a->FormLinearSystem(empty_tdof_list, x, *b, A, X, B);
259 
260  // 10. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
261  // preconditioner from hypre. Extract the parallel grid function x
262  // corresponding to the finite element approximation X. This is the local
263  // solution on each processor.
264  HypreSolver *amg = new HypreBoomerAMG(A);
265  HyprePCG *pcg = new HyprePCG(A);
266  pcg->SetTol(1e-12);
267  pcg->SetMaxIter(200);
268  pcg->SetPrintLevel(2);
269  pcg->SetPreconditioner(*amg);
270  pcg->Mult(B, X);
271  a->RecoverFEMSolution(X, *b, x);
272 
273  delete a;
274  delete b;
275 
276  // 11. Compute and print the L^2 norm of the error.
277  double err = x.ComputeL2Error(sol_coef);
278  if (myid == 0)
279  {
280  cout << "\nL2 norm of error: " << err << endl;
281  }
282 
283  // 12. Save the refined mesh and the solution. This output can be viewed
284  // later using GLVis: "glvis -np <np> -m sphere_refined -g sol".
285  {
286  ostringstream mesh_name, sol_name;
287  mesh_name << "sphere_refined." << setfill('0') << setw(6) << myid;
288  sol_name << "sol." << setfill('0') << setw(6) << myid;
289 
290  ofstream mesh_ofs(mesh_name.str().c_str());
291  mesh_ofs.precision(8);
292  pmesh->Print(mesh_ofs);
293 
294  ofstream sol_ofs(sol_name.str().c_str());
295  sol_ofs.precision(8);
296  x.Save(sol_ofs);
297  }
298 
299  // 13. Send the solution by socket to a GLVis server.
300  if (visualization)
301  {
302  char vishost[] = "localhost";
303  int visport = 19916;
304  socketstream sol_sock(vishost, visport);
305  sol_sock << "parallel " << num_procs << " " << myid << "\n";
306  sol_sock.precision(8);
307  sol_sock << "solution\n" << *pmesh << x << flush;
308  }
309 
310  // 14. Free the used memory.
311  delete pcg;
312  delete amg;
313  delete fespace;
314  delete pmesh;
315 
316  MPI_Finalize();
317 
318  return 0;
319 }
320 
321 double analytic_solution(const Vector &x)
322 {
323  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
324  return x(0)*x(1)/l2;
325 }
326 
327 double analytic_rhs(const Vector &x)
328 {
329  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
330  return 7*x(0)*x(1)/l2;
331 }
332 
333 void SnapNodes(Mesh &mesh)
334 {
335  GridFunction &nodes = *mesh.GetNodes();
336  Vector node(mesh.SpaceDimension());
337  for (int i = 0; i < nodes.FESpace()->GetNDofs(); i++)
338  {
339  for (int d = 0; d < mesh.SpaceDimension(); d++)
340  {
341  node(d) = nodes(nodes.FESpace()->DofToVDof(i, d));
342  }
343 
344  node /= node.Norml2();
345 
346  for (int d = 0; d < mesh.SpaceDimension(); d++)
347  {
348  nodes(nodes.FESpace()->DofToVDof(i, d)) = node(d);
349  }
350  }
351 }
void SetTol(double tol)
Definition: hypre.cpp:1993
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:161
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:104
double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:177
Subclass constant coefficient.
Definition: coefficient.hpp:57
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:667
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:42
void SnapNodes(Mesh &mesh)
Definition: ex7p.cpp:333
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:314
Abstract parallel finite element space.
Definition: pfespace.hpp:28
Data type for vertex.
Definition: vertex.hpp:21
STL namespace.
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2008
void AddVertex(const double *)
Definition: mesh.cpp:910
The BoomerAMG solver in hypre.
Definition: hypre.hpp:770
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:1054
double analytic_solution(const Vector &x)
Definition: ex7p.cpp:321
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6684
void Assemble(int skip_zeros=1)
Assemble the local matrix.
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:6367
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3261
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:293
int Dimension() const
Definition: mesh.hpp:611
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetMaxIter(int max_iter)
Definition: hypre.cpp:1998
int SpaceDimension() const
Definition: mesh.hpp:612
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator.
Definition: linearform.cpp:19
PCG solver in hypre.
Definition: hypre.hpp:630
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)
Definition: optparser.hpp:74
void AddQuad(const int *vi, int attr=1)
Definition: mesh.cpp:931
void RefineAtVertex(const Vertex &vert, double eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex. Uses GeneralRefinement.
Definition: mesh.cpp:6386
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
double analytic_rhs(const Vector &x)
Definition: ex7p.cpp:327
void AddTriangle(const int *vi, int attr=1)
Definition: mesh.cpp:926
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition: mesh.cpp:1086
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:2013
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
Class for parallel bilinear form.
Abstract class for hypre's solvers and preconditioners.
Definition: hypre.hpp:594
class for C-function coefficient
Vector data type.
Definition: vector.hpp:36
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5368
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:146
Class for parallel grid function.
Definition: pgridfunc.hpp:31
Wrapper for hypre's ParCSR matrix class.
Definition: hypre.hpp:174
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre's PCG.
Definition: hypre.cpp:2034
Class for parallel meshes.
Definition: pmesh.hpp:28
int main(int argc, char *argv[])
Definition: ex7p.cpp:35
virtual void Print(std::ostream &out=std::cout) const
Definition: pmesh.cpp:3395
bool Good() const
Definition: optparser.hpp:120