MFEM  v3.1 Finite element discretization library
ex6p.cpp
Go to the documentation of this file.
1 // MFEM Example 6 - Parallel Version
2 //
3 // Compile with: make ex6p
4 //
5 // Sample runs: mpirun -np 4 ex6p -m ../data/square-disc.mesh -o 1
6 // mpirun -np 4 ex6p -m ../data/square-disc.mesh -o 2
7 // mpirun -np 4 ex6p -m ../data/square-disc-nurbs.mesh -o 2
8 // mpirun -np 4 ex6p -m ../data/star.mesh -o 3
9 // mpirun -np 4 ex6p -m ../data/escher.mesh -o 1
10 // mpirun -np 4 ex6p -m ../data/fichera.mesh -o 2
11 // mpirun -np 4 ex6p -m ../data/disc-nurbs.mesh -o 2
12 // mpirun -np 4 ex6p -m ../data/ball-nurbs.mesh
13 // mpirun -np 4 ex6p -m ../data/pipe-nurbs.mesh
14 // mpirun -np 4 ex6p -m ../data/star-surf.mesh -o 2
15 // mpirun -np 4 ex6p -m ../data/square-disc-surf.mesh -o 2
16 // mpirun -np 4 ex6p -m ../data/amr-quad.mesh
17 //
18 // Description: This is a version of Example 1 with a simple adaptive mesh
19 // refinement loop. The problem being solved is again the Laplace
20 // equation -Delta u = 1 with homogeneous Dirichlet boundary
21 // conditions. The problem is solved on a sequence of meshes which
22 // are locally refined in a conforming (triangles, tetrahedrons)
23 // or non-conforming (quadrilateral, hexahedrons) manner according
24 // to a simple ZZ error estimator.
25 //
26 // The example demonstrates MFEM's capability to work with both
27 // conforming and nonconforming refinements, in 2D and 3D, on
28 // linear, curved and surface meshes. Interpolation of functions
29 // from coarse to fine meshes, as well as persistent GLVis
30 // visualization are also illustrated.
31 //
32 // We recommend viewing Example 1 before viewing this example.
33
34 #include "mfem.hpp"
35 #include <fstream>
36 #include <iostream>
37
38 using namespace std;
39 using namespace mfem;
40
41 int main(int argc, char *argv[])
42 {
43  // 1. Initialize MPI.
44  int num_procs, myid;
45  MPI_Init(&argc, &argv);
46  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
47  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
48
49  // 2. Parse command-line options.
50  const char *mesh_file = "../data/star.mesh";
51  int order = 1;
52  bool visualization = true;
53
54  OptionsParser args(argc, argv);
56  "Mesh file to use.");
58  "Finite element order (polynomial degree).");
60  "--no-visualization",
61  "Enable or disable GLVis visualization.");
62  args.Parse();
63  if (!args.Good())
64  {
65  if (myid == 0)
66  {
67  args.PrintUsage(cout);
68  }
69  MPI_Finalize();
70  return 1;
71  }
72  if (myid == 0)
73  {
74  args.PrintOptions(cout);
75  }
76
77  // 3. Read the (serial) mesh from the given mesh file on all processors. We
78  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
79  // and volume meshes with the same code.
80  Mesh *mesh;
81  ifstream imesh(mesh_file);
82  if (!imesh)
83  {
84  if (myid == 0)
85  {
86  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
87  }
88  MPI_Finalize();
89  return 2;
90  }
91  mesh = new Mesh(imesh, 1, 1);
92  imesh.close();
93  int dim = mesh->Dimension();
94  int sdim = mesh->SpaceDimension();
95
96  // 4. Refine the serial mesh on all processors to increase the resolution.
97  // Also project a NURBS mesh to a piecewise-quadratic curved mesh. Make
98  // sure that the mesh is non-conforming.
99  if (mesh->NURBSext)
100  {
101  mesh->UniformRefinement();
102  mesh->SetCurvature(2);
103  }
104  mesh->EnsureNCMesh();
105
106  // 5. Define a parallel mesh by partitioning the serial mesh.
107  // Once the parallel mesh is defined, the serial mesh can be deleted.
108  ParMesh pmesh(MPI_COMM_WORLD, *mesh);
109  delete mesh;
110
111  MFEM_VERIFY(pmesh.bdr_attributes.Size() > 0,
112  "Boundary attributes required in the mesh.");
113  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
114  ess_bdr = 1;
115
116  // 6. Define a finite element space on the mesh. The polynomial order is
117  // one (linear) by default, but this can be changed on the command line.
118  H1_FECollection fec(order, dim);
119  ParFiniteElementSpace fespace(&pmesh, &fec);
120
121  // 7. As in Example 1p, we set up bilinear and linear forms corresponding to
122  // the Laplace problem -\Delta u = 1. We don't assemble the discrete
123  // problem yet, this will be done in the main loop.
124  ParBilinearForm a(&fespace);
125  ParLinearForm b(&fespace);
126
127  ConstantCoefficient one(1.0);
128
131
132  // 8. The solution vector x and the associated finite element grid function
133  // will be maintained over the AMR iterations. We initialize it to zero.
134  ParGridFunction x(&fespace);
135  x = 0;
136
137  // 9. Connect to GLVis.
138  char vishost[] = "localhost";
139  int visport = 19916;
140
141  socketstream sout;
142  if (visualization)
143  {
144  sout.open(vishost, visport);
145  if (!sout)
146  {
147  if (myid == 0)
148  {
149  cout << "Unable to connect to GLVis server at "
150  << vishost << ':' << visport << endl;
151  cout << "GLVis visualization disabled.\n";
152  }
153  visualization = false;
154  }
155
156  sout.precision(8);
157  }
158
159  // 10. The main AMR loop. In each iteration we solve the problem on the
160  // current mesh, visualize the solution, estimate the error on all
161  // elements, refine the worst elements and update all objects to work
162  // with the new mesh.
163  const int max_dofs = 100000;
164  for (int it = 0; ; it++)
165  {
166  HYPRE_Int global_dofs = fespace.GlobalTrueVSize();
167  if (myid == 0)
168  {
169  cout << "\nIteration " << it << endl;
170  cout << "Number of unknowns: " << global_dofs << endl;
171  }
172
173  // 11. Assemble the stiffness matrix and the right-hand side. Note that
174  // MFEM doesn't care at this point that the mesh is nonconforming
175  // and parallel. The FE space is considered 'cut' along hanging
176  // edges/faces, and also across processor boundaries.
177  a.Assemble();
178  b.Assemble();
179
180  // 12. Set the initial estimate of the solution and the Dirichlet DOFs,
181  // here we just use zero everywhere.
182  x = 0.0;
183
184  // 13. Create the parallel linear system: eliminate boundary conditions,
185  // constrain hanging nodes and nodes across processor boundaries.
186  // The system will be solved for true (unconstrained/unique) DOFs only.
187  Array<int> ess_tdof_list;
188  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
189
190  HypreParMatrix A;
191  Vector B, X;
192  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
193
194  // 14. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
195  // preconditioner from hypre.
196  HypreBoomerAMG amg(A);
197  amg.SetPrintLevel(0);
198  HyprePCG pcg(A);
199  pcg.SetTol(1e-12);
200  pcg.SetMaxIter(200);
201  pcg.SetPrintLevel(0);
202  pcg.SetPreconditioner(amg);
203  pcg.Mult(B, X);
204
205  // 15. Extract the parallel grid function corresponding to the finite element
206  // approximation X. This is the local solution on each processor.
207  a.RecoverFEMSolution(X, b, x);
208
209  // 16. Send the solution by socket to a GLVis server.
210  if (visualization)
211  {
212  sout << "parallel " << num_procs << " " << myid << "\n";
213  sout << "solution\n" << pmesh << x << flush;
214  }
215
216  if (global_dofs > max_dofs)
217  {
218  break;
219  }
220
221  // 17. Estimate element errors using the Zienkiewicz-Zhu error estimator.
222  // The bilinear form integrator must have the 'ComputeElementFlux'
223  // method defined.
224  Vector errors(pmesh.GetNE());
225  {
226  // Space for the discontinuous (original) flux
227  DiffusionIntegrator flux_integrator(one);
228  L2_FECollection flux_fec(order, dim);
229  ParFiniteElementSpace flux_fes(&pmesh, &flux_fec, sdim);
230
231  // Space for the smoothed (conforming) flux
232  double norm_p = 1;
233  RT_FECollection smooth_flux_fec(order-1, dim);
234  ParFiniteElementSpace smooth_flux_fes(&pmesh, &smooth_flux_fec);
235
236  // Another possible set of options for the smoothed flux space:
237  // norm_p = 1;
238  // H1_FECollection smooth_flux_fec(order, dim);
239  // ParFiniteElementSpace smooth_flux_fes(&pmesh, &smooth_flux_fec, dim);
240
241  L2ZZErrorEstimator(flux_integrator, x,
242  smooth_flux_fes, flux_fes, errors, norm_p);
243  }
244  double local_max_err = errors.Max();
245  double global_max_err;
246  MPI_Allreduce(&local_max_err, &global_max_err, 1,
247  MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
248
249  // 18. Make a list of elements whose error is larger than a fraction
250  // of the maximum element error. These elements will be refined.
251  Array<int> ref_list;
252  const double frac = 0.7;
253  double threshold = frac * global_max_err;
254  for (int i = 0; i < errors.Size(); i++)
255  {
256  if (errors[i] >= threshold) { ref_list.Append(i); }
257  }
258
259  // 19. Refine the selected elements. Since we are going to transfer the
260  // grid function x from the coarse mesh to the new fine mesh in the
261  // next step, we need to request the "two-level state" of the mesh.
262  pmesh.GeneralRefinement(ref_list);
263
264  // 20. Inform the space, grid function and also the bilinear and linear
265  // forms that the space has changed.
266  fespace.Update();
267  x.Update();
268  a.Update();
269  b.Update();
270  }
271
272  MPI_Finalize();
273  return 0;
274 }
void SetTol(double tol)
Definition: hypre.cpp:1917
int Size() const
Logical size of the array.
Definition: array.hpp:109
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
Subclass constant coefficient.
Definition: coefficient.hpp:57
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:34
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:454
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:1932
void Update(ParFiniteElementSpace *pf=NULL)
Definition: plinearform.cpp:21
virtual void Update(FiniteElementSpace *nfes=NULL)
The BoomerAMG solver in hypre.
Definition: hypre.hpp:698
Class for parallel linear form.
Definition: plinearform.hpp:26
int dim
Definition: ex3.cpp:48
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:368
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7316
void SetPrintLevel(int print_level)
Definition: hypre.hpp:739
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3736
void L2ZZErrorEstimator(BilinearFormIntegrator &flux_integrator, ParGridFunction &x, ParFiniteElementSpace &smooth_flux_fes, ParFiniteElementSpace &flux_fes, Vector &errors, int norm_p, double solver_tol, int solver_max_it)
Definition: pgridfunc.cpp:535
T Max() const
Definition: array.cpp:90
void Assemble(int skip_zeros=1)
Assemble the local matrix.
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
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:133
int SpaceDimension() const
Definition: mesh.hpp:476
Definition: linearform.cpp:19
Array< int > bdr_attributes
Definition: mesh.hpp:140
int main(int argc, char *argv[])
Definition: ex1.cpp:45
MPI_Comm GetComm() const
Definition: pmesh.hpp:86
PCG solver in hypre.
Definition: hypre.hpp:558
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list)
Definition: pfespace.cpp:545
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 RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
NURBSExtension * NURBSext
Definition: mesh.hpp:142
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
Class for parallel bilinear form.
int open(const char hostname[], int port)
void EnsureNCMesh()
Definition: mesh.cpp:6929
Vector data type.
Definition: vector.hpp:33
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
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:6864
Class for parallel meshes.
Definition: pmesh.hpp:28
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:95
bool Good() const
Definition: optparser.hpp:120