 MFEM  v4.4.0 Finite element discretization library
ex6.cpp
Go to the documentation of this file.
1 // MFEM Example 6
2 //
3 // Compile with: make ex6
4 //
5 // Sample runs: ex6 -m ../data/square-disc.mesh -o 1
6 // ex6 -m ../data/square-disc.mesh -o 2
7 // ex6 -m ../data/square-disc-nurbs.mesh -o 2
8 // ex6 -m ../data/star.mesh -o 3
9 // ex6 -m ../data/escher.mesh -o 2
10 // ex6 -m ../data/fichera.mesh -o 2
11 // ex6 -m ../data/disc-nurbs.mesh -o 2
12 // ex6 -m ../data/ball-nurbs.mesh
13 // ex6 -m ../data/pipe-nurbs.mesh
14 // ex6 -m ../data/star-surf.mesh -o 2
15 // ex6 -m ../data/square-disc-surf.mesh -o 2
17 // ex6 -m ../data/inline-segment.mesh -o 1 -md 100
18 //
19 // Device sample runs:
20 // ex6 -pa -d cuda
21 // ex6 -pa -d occa-cuda
22 // ex6 -pa -d raja-omp
23 // ex6 -pa -d ceed-cpu
24 // * ex6 -pa -d ceed-cuda
25 // ex6 -pa -d ceed-cuda:/gpu/cuda/shared
26 //
27 // Description: This is a version of Example 1 with a simple adaptive mesh
28 // refinement loop. The problem being solved is again the Laplace
29 // equation -Delta u = 1 with homogeneous Dirichlet boundary
30 // conditions. The problem is solved on a sequence of meshes which
31 // are locally refined in a conforming (triangles, tetrahedrons)
32 // or non-conforming (quadrilaterals, hexahedra) manner according
33 // to a simple ZZ error estimator.
34 //
35 // The example demonstrates MFEM's capability to work with both
36 // conforming and nonconforming refinements, in 2D and 3D, on
37 // linear, curved and surface meshes. Interpolation of functions
38 // from coarse to fine meshes, as well as persistent GLVis
39 // visualization are also illustrated.
40 //
41 // We recommend viewing Example 1 before viewing this example.
42
43 #include "mfem.hpp"
44 #include <fstream>
45 #include <iostream>
46
47 using namespace std;
48 using namespace mfem;
49
50 int main(int argc, char *argv[])
51 {
52  // 1. Parse command-line options.
53  const char *mesh_file = "../data/star.mesh";
54  int order = 1;
55  bool pa = false;
56  const char *device_config = "cpu";
57  int max_dofs = 50000;
58  bool visualization = true;
59
60  OptionsParser args(argc, argv);
62  "Mesh file to use.");
64  "Finite element order (polynomial degree).");
66  "--no-partial-assembly", "Enable Partial Assembly.");
68  "Device configuration string, see Device::Configure().");
70  "Stop after reaching this many degrees of freedom.");
72  "--no-visualization",
73  "Enable or disable GLVis visualization.");
74  args.Parse();
75  if (!args.Good())
76  {
77  args.PrintUsage(cout);
78  return 1;
79  }
80  args.PrintOptions(cout);
81
82  // 2. Enable hardware devices such as GPUs, and programming models such as
83  // CUDA, OCCA, RAJA and OpenMP based on command line options.
84  Device device(device_config);
85  device.Print();
86
87  // 3. Read the mesh from the given mesh file. We can handle triangular,
88  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
89  // the same code.
90  Mesh mesh(mesh_file, 1, 1);
91  int dim = mesh.Dimension();
92  int sdim = mesh.SpaceDimension();
93
94  // 4. Since a NURBS mesh can currently only be refined uniformly, we need to
95  // convert it to a piecewise-polynomial curved mesh. First we refine the
96  // NURBS mesh a bit more and then project the curvature to quadratic Nodes.
97  if (mesh.NURBSext)
98  {
99  for (int i = 0; i < 2; i++)
100  {
101  mesh.UniformRefinement();
102  }
103  mesh.SetCurvature(2);
104  }
105
106  // 5. Define a finite element space on the mesh. The polynomial order is
107  // one (linear) by default, but this can be changed on the command line.
108  H1_FECollection fec(order, dim);
109  FiniteElementSpace fespace(&mesh, &fec);
110
111  // 6. As in Example 1, we set up bilinear and linear forms corresponding to
112  // the Laplace problem -\Delta u = 1. We don't assemble the discrete
113  // problem yet, this will be done in the main loop.
114  BilinearForm a(&fespace);
115  if (pa)
116  {
117  a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
118  a.SetDiagonalPolicy(Operator::DIAG_ONE);
119  }
120  LinearForm b(&fespace);
121
122  ConstantCoefficient one(1.0);
123  ConstantCoefficient zero(0.0);
124
128
129  // 7. The solution vector x and the associated finite element grid function
130  // will be maintained over the AMR iterations. We initialize it to zero.
131  GridFunction x(&fespace);
132  x = 0.0;
133
134  // 8. All boundary attributes will be used for essential (Dirichlet) BC.
135  MFEM_VERIFY(mesh.bdr_attributes.Size() > 0,
136  "Boundary attributes required in the mesh.");
137  Array<int> ess_bdr(mesh.bdr_attributes.Max());
138  ess_bdr = 1;
139
140  // 9. Connect to GLVis.
141  char vishost[] = "localhost";
142  int visport = 19916;
143  socketstream sol_sock;
144  if (visualization)
145  {
146  sol_sock.open(vishost, visport);
147  }
148
149  // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
150  // that uses the ComputeElementFlux method of the DiffusionIntegrator to
151  // recover a smoothed flux (gradient) that is subtracted from the element
152  // flux to get an error indicator. We need to supply the space for the
153  // smoothed flux: an (H1)^sdim (i.e., vector-valued) space is used here.
154  FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
155  ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
156  estimator.SetAnisotropic();
157
158  // 11. A refiner selects and refines elements based on a refinement strategy.
159  // The strategy here is to refine elements with errors larger than a
160  // fraction of the maximum element error. Other strategies are possible.
161  // The refiner will call the given error estimator.
162  ThresholdRefiner refiner(estimator);
163  refiner.SetTotalErrorFraction(0.7);
164
165  // 12. The main AMR loop. In each iteration we solve the problem on the
166  // current mesh, visualize the solution, and refine the mesh.
167  for (int it = 0; ; it++)
168  {
169  int cdofs = fespace.GetTrueVSize();
170  cout << "\nAMR iteration " << it << endl;
171  cout << "Number of unknowns: " << cdofs << endl;
172
173  // 13. Assemble the right-hand side.
174  b.Assemble();
175
176  // 14. Set Dirichlet boundary values in the GridFunction x.
177  // Determine the list of Dirichlet true DOFs in the linear system.
178  Array<int> ess_tdof_list;
179  x.ProjectBdrCoefficient(zero, ess_bdr);
180  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
181
182  // 15. Assemble the stiffness matrix.
183  a.Assemble();
184
185  // 16. Create the linear system: eliminate boundary conditions, constrain
186  // hanging nodes and possibly apply other transformations. The system
187  // will be solved for true (unconstrained) DOFs only.
188  OperatorPtr A;
189  Vector B, X;
190
191  const int copy_interior = 1;
192  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
193
194  // 17. Solve the linear system A X = B.
195  if (!pa)
196  {
197 #ifndef MFEM_USE_SUITESPARSE
198  // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
199  GSSmoother M((SparseMatrix&)(*A));
200  PCG(*A, M, B, X, 3, 200, 1e-12, 0.0);
201 #else
202  // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
203  UMFPackSolver umf_solver;
204  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
205  umf_solver.SetOperator(*A);
206  umf_solver.Mult(B, X);
207 #endif
208  }
209  else // Diagonal preconditioning in partial assembly mode.
210  {
211  OperatorJacobiSmoother M(a, ess_tdof_list);
212  PCG(*A, M, B, X, 3, 2000, 1e-12, 0.0);
213  }
214
215  // 18. After solving the linear system, reconstruct the solution as a
216  // finite element GridFunction. Constrained nodes are interpolated
217  // from true DOFs (it may therefore happen that x.Size() >= X.Size()).
218  a.RecoverFEMSolution(X, b, x);
219
220  // 19. Send solution by socket to the GLVis server.
221  if (visualization && sol_sock.good())
222  {
223  sol_sock.precision(8);
224  sol_sock << "solution\n" << mesh << x << flush;
225  }
226
227  if (cdofs > max_dofs)
228  {
229  cout << "Reached the maximum number of dofs. Stop." << endl;
230  break;
231  }
232
233  // 20. Call the refiner to modify the mesh. The refiner calls the error
234  // estimator to obtain element errors, then it selects elements to be
235  // refined and finally it modifies the mesh. The Stop() method can be
236  // used to determine if a stopping criterion was met.
237  refiner.Apply(mesh);
238  if (refiner.Stop())
239  {
240  cout << "Stopping criterion satisfied. Stop." << endl;
241  break;
242  }
243
244  // 21. Update the space to reflect the new state of the mesh. Also,
245  // interpolate the solution x so that it lies in the new space but
246  // represents the same function. This saves solver iterations later
247  // since we'll have a good initial guess of x in the next step.
248  // Internally, FiniteElementSpace::Update() calculates an
249  // interpolation matrix which is then used by GridFunction::Update().
250  fespace.Update();
251  x.Update();
252
253  // 22. Inform also the bilinear and linear forms that the space has
254  // changed.
255  a.Update();
256  b.Update();
257  }
258
259  return 0;
260 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
void SetAnisotropic(bool aniso=true)
Enable/disable anisotropic estimates. To enable this option, the BilinearFormIntegrator must support ...
Definition: estimators.hpp:165
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3336
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:102
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
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 RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition: fespace.cpp:564
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:279
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
void SetDiagonalPolicy(DiagonalPolicy policy)
Sets diagonal policy used upon construction of the linear system.
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:1025
bool Apply(Mesh &mesh)
Perform the mesh operation.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
Data type sparse matrix.
Definition: sparsemat.hpp:46
constexpr char vishost[]
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:274
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
constexpr int visport
Mesh refinement operator using an error threshold.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:5312
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
virtual void Update(FiniteElementSpace *nfes=NULL)
Update the FiniteElementSpace and delete all data associated with the old one.
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition: solvers.cpp:904
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:566
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
int SpaceDimension() const
Definition: mesh.hpp:1000
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:34
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:1036
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
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)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:164
void Update()
Update the object according to the associated FE space fes.
Definition: linearform.hpp:188
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
int dim
Definition: ex24.cpp:53
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
The ZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu error estimation procedure...
Definition: estimators.hpp:88
Vector data type.
Definition: vector.hpp:60
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:23
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:121
int main()
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3179
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
Definition: gridfunc.hpp:447
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:3084
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150