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