MFEM  v4.0 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 2
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 // Device sample runs:
19 // mpirun -np 4 ex6p -pa -d cuda
20 // mpirun -np 4 ex6p -pa -d occa-cuda
21 // mpirun -np 4 ex6p -pa -d raja-omp
22 //
23 // Description: This is a version of Example 1 with a simple adaptive mesh
24 // refinement loop. The problem being solved is again the Laplace
25 // equation -Delta u = 1 with homogeneous Dirichlet boundary
26 // conditions. The problem is solved on a sequence of meshes which
27 // are locally refined in a conforming (triangles, tetrahedrons)
28 // or non-conforming (quadrilaterals, hexahedra) manner according
29 // to a simple ZZ error estimator.
30 //
31 // The example demonstrates MFEM's capability to work with both
32 // conforming and nonconforming refinements, in 2D and 3D, on
33 // linear, curved and surface meshes. Interpolation of functions
34 // from coarse to fine meshes, as well as persistent GLVis
35 // visualization are also illustrated.
36 //
37 // We recommend viewing Example 1 before viewing this example.
38
39 #include "mfem.hpp"
40 #include <fstream>
41 #include <iostream>
42
43 using namespace std;
44 using namespace mfem;
45
46 int main(int argc, char *argv[])
47 {
48  // 1. Initialize MPI.
49  int num_procs, myid;
50  MPI_Init(&argc, &argv);
51  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
52  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
53
54  // 2. Parse command-line options.
55  const char *mesh_file = "../data/star.mesh";
56  int order = 1;
57  bool pa = false;
58  const char *device_config = "cpu";
59  bool visualization = true;
60
61  OptionsParser args(argc, argv);
63  "Mesh file to use.");
65  "Finite element order (polynomial degree).");
67  "--no-partial-assembly", "Enable Partial Assembly.");
69  "Device configuration string, see Device::Configure().");
71  "--no-visualization",
72  "Enable or disable GLVis visualization.");
73  args.Parse();
74  if (!args.Good())
75  {
76  if (myid == 0)
77  {
78  args.PrintUsage(cout);
79  }
80  MPI_Finalize();
81  return 1;
82  }
83  if (myid == 0)
84  {
85  args.PrintOptions(cout);
86  }
87
88  // 3. Enable hardware devices such as GPUs, and programming models such as
89  // CUDA, OCCA, RAJA and OpenMP based on command line options.
90  Device device(device_config);
91  if (myid == 0) { device.Print(); }
92
93  // 4. Read the (serial) mesh from the given mesh file on all processors. We
94  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
95  // and volume meshes with the same code.
96  Mesh *mesh = new Mesh(mesh_file, 1, 1);
97  int dim = mesh->Dimension();
98  int sdim = mesh->SpaceDimension();
99
100  // 5. Refine the serial mesh on all processors to increase the resolution.
101  // Also project a NURBS mesh to a piecewise-quadratic curved mesh. Make
102  // sure that the mesh is non-conforming.
103  if (mesh->NURBSext)
104  {
105  mesh->UniformRefinement();
106  mesh->SetCurvature(2);
107  }
108  mesh->EnsureNCMesh();
109
110  // 6. Define a parallel mesh by partitioning the serial mesh.
111  // Once the parallel mesh is defined, the serial mesh can be deleted.
112  ParMesh pmesh(MPI_COMM_WORLD, *mesh);
113  delete mesh;
114
115  MFEM_VERIFY(pmesh.bdr_attributes.Size() > 0,
116  "Boundary attributes required in the mesh.");
117  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
118  ess_bdr = 1;
119
120  // 7. Define a finite element space on the mesh. The polynomial order is
121  // one (linear) by default, but this can be changed on the command line.
122  H1_FECollection fec(order, dim);
123  ParFiniteElementSpace fespace(&pmesh, &fec);
124
125  // 8. As in Example 1p, we set up bilinear and linear forms corresponding to
126  // the Laplace problem -\Delta u = 1. We don't assemble the discrete
127  // problem yet, this will be done in the main loop.
128  ParBilinearForm a(&fespace);
129  if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
130  ParLinearForm b(&fespace);
131
132  ConstantCoefficient one(1.0);
133
137
138  // 9. The solution vector x and the associated finite element grid function
139  // will be maintained over the AMR iterations. We initialize it to zero.
140  ParGridFunction x(&fespace);
141  x = 0;
142
143  // 10. Connect to GLVis.
144  char vishost[] = "localhost";
145  int visport = 19916;
146
147  socketstream sout;
148  if (visualization)
149  {
150  sout.open(vishost, visport);
151  if (!sout)
152  {
153  if (myid == 0)
154  {
155  cout << "Unable to connect to GLVis server at "
156  << vishost << ':' << visport << endl;
157  cout << "GLVis visualization disabled.\n";
158  }
159  visualization = false;
160  }
161
162  sout.precision(8);
163  }
164
165  // 11. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
166  // with L2 projection in the smoothing step to better handle hanging
167  // nodes and parallel partitioning. We need to supply a space for the
168  // discontinuous flux (L2) and a space for the smoothed flux (H(div) is
169  // used here).
170  L2_FECollection flux_fec(order, dim);
171  ParFiniteElementSpace flux_fes(&pmesh, &flux_fec, sdim);
172  RT_FECollection smooth_flux_fec(order-1, dim);
173  ParFiniteElementSpace smooth_flux_fes(&pmesh, &smooth_flux_fec);
174  // Another possible option for the smoothed flux space:
175  // H1_FECollection smooth_flux_fec(order, dim);
176  // ParFiniteElementSpace smooth_flux_fes(&pmesh, &smooth_flux_fec, dim);
177  L2ZienkiewiczZhuEstimator estimator(*integ, x, flux_fes, smooth_flux_fes);
178
179  // 12. A refiner selects and refines elements based on a refinement strategy.
180  // The strategy here is to refine elements with errors larger than a
181  // fraction of the maximum element error. Other strategies are possible.
182  // The refiner will call the given error estimator.
183  ThresholdRefiner refiner(estimator);
184  refiner.SetTotalErrorFraction(0.7);
185
186  // 13. The main AMR loop. In each iteration we solve the problem on the
187  // current mesh, visualize the solution, and refine the mesh.
188  const int max_dofs = 100000;
189  for (int it = 0; ; it++)
190  {
191  HYPRE_Int global_dofs = fespace.GlobalTrueVSize();
192  if (myid == 0)
193  {
194  cout << "\nAMR iteration " << it << endl;
195  cout << "Number of unknowns: " << global_dofs << endl;
196  }
197
198  // 14. Assemble the right-hand side and determine the list of true
199  // (i.e. parallel conforming) essential boundary dofs.
200  Array<int> ess_tdof_list;
201  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
202  b.Assemble();
203
204  // 15. Assemble the stiffness matrix. Note that MFEM doesn't care at this
205  // point that the mesh is nonconforming and parallel. The FE space is
206  // considered 'cut' along hanging edges/faces, and also across
207  // processor boundaries.
208  a.Assemble();
209
210  // 16. Create the parallel linear system: eliminate boundary conditions.
211  // The system will be solved for true (unconstrained/unique) DOFs only.
212  OperatorPtr A;
213  Vector B, X;
214
215  const int copy_interior = 1;
216  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
217
218  // 17. Solve the linear system A X = B.
219  // * With full assembly, use the BoomerAMG preconditioner from hypre.
220  // * With partial assembly, use no preconditioner, for now.
221  HypreBoomerAMG *amg = NULL;
222  if (!pa) { amg = new HypreBoomerAMG; amg->SetPrintLevel(0); }
223  CGSolver cg(MPI_COMM_WORLD);
224  cg.SetRelTol(1e-6);
225  cg.SetMaxIter(2000);
226  cg.SetPrintLevel(3); // print the first and the last iterations only
227  if (amg) { cg.SetPreconditioner(*amg); }
228  cg.SetOperator(*A);
229  cg.Mult(B, X);
230  delete amg;
231
232  // 18. Switch back to the host and extract the parallel grid function
233  // corresponding to the finite element approximation X. This is the
234  // local solution on each processor.
235  a.RecoverFEMSolution(X, b, x);
236
237  // 19. Send the solution by socket to a GLVis server.
238  if (visualization)
239  {
240  sout << "parallel " << num_procs << " " << myid << "\n";
241  sout << "solution\n" << pmesh << x << flush;
242  }
243
244  if (global_dofs > max_dofs)
245  {
246  if (myid == 0)
247  {
248  cout << "Reached the maximum number of dofs. Stop." << endl;
249  }
250  break;
251  }
252
253  // 20. Call the refiner to modify the mesh. The refiner calls the error
254  // estimator to obtain element errors, then it selects elements to be
255  // refined and finally it modifies the mesh. The Stop() method can be
256  // used to determine if a stopping criterion was met.
257  refiner.Apply(pmesh);
258  if (refiner.Stop())
259  {
260  if (myid == 0)
261  {
262  cout << "Stopping criterion satisfied. Stop." << endl;
263  }
264  break;
265  }
266
267  // 21. Update the finite element space (recalculate the number of DOFs,
268  // etc.) and create a grid function update matrix. Apply the matrix
269  // to any GridFunctions over the space. In this case, the update
270  // matrix is an interpolation matrix so the updated GridFunction will
271  // still represent the same function as before refinement.
272  fespace.Update();
273  x.Update();
274
275  // 22. Load balance the mesh, and update the space and solution. Currently
276  // available only for nonconforming meshes.
277  if (pmesh.Nonconforming())
278  {
279  pmesh.Rebalance();
280
281  // Update the space and the GridFunction. This time the update matrix
282  // redistributes the GridFunction among the processors.
283  fespace.Update();
284  x.Update();
285  }
286
287  // 23. Inform also the bilinear and linear forms that the space has
288  // changed.
289  a.Update();
290  b.Update();
291  }
292
293  MPI_Finalize();
294  return 0;
295 }
int Size() const
Logical size of the array.
Definition: array.hpp:118
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:93
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:737
Definition: solvers.hpp:111
Subclass constant coefficient.
Definition: coefficient.hpp:67
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(...
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:290
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:2754
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:79
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:98
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::FULL.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
int main(int argc, char *argv[])
Definition: ex1.cpp:58
The L2ZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu error estimation procedure where t...
Definition: estimators.hpp:194
void Update(ParFiniteElementSpace *pf=NULL)
Update the object according to the given new FE space *pf.
Definition: plinearform.cpp:21
virtual void Update(FiniteElementSpace *nfes=NULL)
The BoomerAMG solver in hypre.
Definition: hypre.hpp:873
bool Apply(Mesh &mesh)
Perform the mesh operation.
void Rebalance()
Load balance the mesh. NC meshes only.
Definition: pmesh.cpp:3085
bool Nonconforming() const
Definition: mesh.hpp:1136
Class for parallel linear form.
Definition: plinearform.hpp:26
int dim
Definition: ex3.cpp:48
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:67
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:80
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7591
void SetPrintLevel(int print_level)
Definition: hypre.hpp:914
Mesh refinement operator using an error threshold.
void SetMaxIter(int max_it)
Definition: solvers.hpp:63
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3644
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
HYPRE_Int GlobalTrueVSize() const
Definition: pfespace.hpp:251
void Assemble(int skip_zeros=1)
Assemble the local matrix.
int Dimension() const
Definition: mesh.hpp:713
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:180
int SpaceDimension() const
Definition: mesh.hpp:714
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:23
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:179
void SetRelTol(double rtol)
Definition: solvers.hpp:61
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:76
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:181
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
Class for parallel bilinear form.
void EnsureNCMesh(bool triangles_nonconforming=false)
Definition: mesh.cpp:7253
int open(const char hostname[], int port)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:125
Vector data type.
Definition: vector.hpp:48
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:88
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void SetTotalErrorFraction(double fraction)
Set the total fraction used in the computation of the threshold. The default value is 1/2...
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:96
Class for parallel meshes.
Definition: pmesh.hpp:32
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:134
bool Good() const
Definition: optparser.hpp:122