MFEM  v3.1 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);
53  "Type of elements to use: 0 - triangles, 1 - quads.");
55  "Finite element order (polynomial degree).");
57  "Number of times to refine the mesh uniformly.");
60  " 1 = refine around north pole, 2 = refine randomly.");
62  "--no-visualization",
63  "Enable or disable GLVis visualization.");
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  {
111  }
112  for (int j = 0; j < Nelem; j++)
113  {
114  int attribute = j + 1;
116  }
117  mesh->FinalizeTriMesh(1, 1, true);
118  }
119  else // inscribed cube
120  {
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  };
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  {
135  }
136  for (int j = 0; j < Nelem; j++)
137  {
138  int attribute = j + 1;
140  }
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  mesh->RefineAtVertex(Vertex(0, 0, 1), 3); // 3 levels
168  SnapNodes(*mesh);
169  }
170  else if (amr == 2)
171  {
172  mesh->RandomRefinement(2, 2); // 2 levels, 1/2 of the elements
173  SnapNodes(*mesh);
174  }
175
176  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
177  delete mesh;
178  {
179  int par_ref_levels = 2;
180  for (int l = 0; l < par_ref_levels; l++)
181  {
182  pmesh->UniformRefinement();
183
184  // Snap the nodes of the refined mesh back to sphere surface.
185  if (always_snap)
186  {
187  SnapNodes(*pmesh);
188  }
189  }
190  if (!always_snap || par_ref_levels < 1)
191  {
192  SnapNodes(*pmesh);
193  }
194  }
195
196  if (amr == 1)
197  {
198  pmesh->RefineAtVertex(Vertex(0, 0, 1), 2); // 2 levels
199  SnapNodes(*pmesh);
200  }
201  else if (amr == 2)
202  {
203  pmesh->RandomRefinement(2, 2); // 2 levels, 1/2 of the elements
204  SnapNodes(*pmesh);
205  }
206
207  // 5. Define a finite element space on the mesh. Here we use isoparametric
208  // finite elements -- the same as the mesh nodes.
209  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, &fec);
210  HYPRE_Int size = fespace->GlobalTrueVSize();
211  if (myid == 0)
212  {
213  cout << "Number of unknowns: " << size << endl;
214  }
215
216  // 6. Set up the linear form b(.) which corresponds to the right-hand side of
217  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
218  // the basis functions in the finite element fespace.
219  ParLinearForm *b = new ParLinearForm(fespace);
220  ConstantCoefficient one(1.0);
224  b->Assemble();
225
226  // 7. Define the solution vector x as a finite element grid function
227  // corresponding to fespace. Initialize x with initial guess of zero.
228  ParGridFunction x(fespace);
229  x = 0.0;
230
231  // 8. Set up the bilinear form a(.,.) on the finite element space
232  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
233  // and Mass domain integrators.
234  ParBilinearForm *a = new ParBilinearForm(fespace);
237
238  // 9. Assemble the parallel linear system, applying any transformations
239  // such as: parallel assembly, applying conforming constraints, etc.
240  a->Assemble();
241  HypreParMatrix A;
242  Vector B, X;
243  Array<int> empty_tdof_list;
244  a->FormLinearSystem(empty_tdof_list, x, *b, A, X, B);
245
246  // 10. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
247  // preconditioner from hypre. Extract the parallel grid function x
248  // corresponding to the finite element approximation X. This is the local
249  // solution on each processor.
250  HypreSolver *amg = new HypreBoomerAMG(A);
251  HyprePCG *pcg = new HyprePCG(A);
252  pcg->SetTol(1e-12);
253  pcg->SetMaxIter(200);
254  pcg->SetPrintLevel(2);
255  pcg->SetPreconditioner(*amg);
256  pcg->Mult(B, X);
257  a->RecoverFEMSolution(X, *b, x);
258
259  delete a;
260  delete b;
261
262  // 11. Compute and print the L^2 norm of the error.
263  double err = x.ComputeL2Error(sol_coef);
264  if (myid == 0)
265  {
266  cout << "\nL2 norm of error: " << err << endl;
267  }
268
269  // 12. Save the refined mesh and the solution. This output can be viewed
270  // later using GLVis: "glvis -np <np> -m sphere_refined -g sol".
271  {
272  ostringstream mesh_name, sol_name;
273  mesh_name << "sphere_refined." << setfill('0') << setw(6) << myid;
274  sol_name << "sol." << setfill('0') << setw(6) << myid;
275
276  ofstream mesh_ofs(mesh_name.str().c_str());
277  mesh_ofs.precision(8);
278  pmesh->Print(mesh_ofs);
279
280  ofstream sol_ofs(sol_name.str().c_str());
281  sol_ofs.precision(8);
282  x.Save(sol_ofs);
283  }
284
285  // 13. Send the solution by socket to a GLVis server.
286  if (visualization)
287  {
288  char vishost[] = "localhost";
289  int visport = 19916;
290  socketstream sol_sock(vishost, visport);
291  sol_sock << "parallel " << num_procs << " " << myid << "\n";
292  sol_sock.precision(8);
293  sol_sock << "solution\n" << *pmesh << x << flush;
294  }
295
296  // 14. Free the used memory.
297  delete pcg;
298  delete amg;
299  delete fespace;
300  delete pmesh;
301
302  MPI_Finalize();
303
304  return 0;
305 }
306
307 double analytic_solution(const Vector &x)
308 {
309  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
310  return x(0)*x(1)/l2;
311 }
312
313 double analytic_rhs(const Vector &x)
314 {
315  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
316  return 7*x(0)*x(1)/l2;
317 }
318
319 void SnapNodes(Mesh &mesh)
320 {
321  GridFunction &nodes = *mesh.GetNodes();
322  Vector node(mesh.SpaceDimension());
323  for (int i = 0; i < nodes.FESpace()->GetNDofs(); i++)
324  {
325  for (int d = 0; d < mesh.SpaceDimension(); d++)
326  {
327  node(d) = nodes(nodes.FESpace()->DofToVDof(i, d));
328  }
329
330  node /= node.Norml2();
331
332  for (int d = 0; d < mesh.SpaceDimension(); d++)
333  {
334  nodes(nodes.FESpace()->DofToVDof(i, d)) = node(d);
335  }
336  }
337 }
void SetTol(double tol)
Definition: hypre.cpp:1917
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:162
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:116
double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:148
Subclass constant coefficient.
Definition: coefficient.hpp:57
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:644
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:34
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:306
Abstract parallel finite element space.
Definition: pfespace.hpp:28
Data type for vertex.
Definition: vertex.hpp:21
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:1932
void RefineAtVertex(const Vertex &vert, int levels, double eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex, &#39;levels&#39; times.
Definition: mesh.cpp:6961
Definition: mesh.cpp:845
The BoomerAMG solver in hypre.
Definition: hypre.hpp:698
Class for parallel linear form.
Definition: plinearform.hpp:26
void FinalizeTriMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Definition: mesh.cpp:992
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7316
void Assemble(int skip_zeros=1)
Assemble the local matrix.
double analytic_solution(const Vector &x)
Definition: ex7.cpp:245
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3719
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:270
int Dimension() const
Definition: mesh.hpp:475
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetMaxIter(int max_iter)
Definition: hypre.cpp:1922
double analytic_rhs(const Vector &x)
Definition: ex7.cpp:251
int SpaceDimension() const
Definition: mesh.hpp:476
Definition: linearform.cpp:19
int main(int argc, char *argv[])
Definition: ex1.cpp:45
PCG solver in hypre.
Definition: hypre.hpp:558
Abstract finite element space.
Definition: fespace.hpp:62
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
Definition: mesh.cpp:866
void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
void RandomRefinement(int levels, int frac=2, bool aniso=false, int nonconforming=-1, int nc_limit=-1, int seed=0)
Refine each element with 1/frac probability, repeat &#39;levels&#39; times.
Definition: mesh.cpp:6938
void AddTriangle(const int *vi, int attr=1)
Definition: mesh.cpp:861
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Definition: mesh.cpp:1020
void FormLinearSystem(Array< int > &ess_tdof_list, Vector &x, Vector &b, HypreParMatrix &A, Vector &X, Vector &B, int copy_interior=0)
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:1937
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
void SnapNodes(Mesh &mesh)
Definition: ex7.cpp:257
Class for parallel bilinear form.
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:522
class for C-function coefficient
Vector data type.
Definition: vector.hpp:33
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5844
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:54
Class for parallel grid function.
Definition: pgridfunc.hpp:31
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:143
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:1958
Class for parallel meshes.
Definition: pmesh.hpp:28
virtual void Print(std::ostream &out=std::cout) const
Definition: pmesh.cpp:2850
bool Good() const
Definition: optparser.hpp:120