MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
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
16 // ex6 -m ../data/amr-quad.mesh
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 LSZZ = false;
59  bool visualization = true;
60 
61  OptionsParser args(argc, argv);
62  args.AddOption(&mesh_file, "-m", "--mesh",
63  "Mesh file to use.");
64  args.AddOption(&order, "-o", "--order",
65  "Finite element order (polynomial degree).");
66  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
67  "--no-partial-assembly", "Enable Partial Assembly.");
68  args.AddOption(&device_config, "-d", "--device",
69  "Device configuration string, see Device::Configure().");
70  args.AddOption(&max_dofs, "-md", "--max-dofs",
71  "Stop after reaching this many degrees of freedom.");
72  args.AddOption(&LSZZ, "-ls", "--ls-zz", "-no-ls",
73  "--no-ls-zz",
74  "Switch to least-squares ZZ estimator.");
75  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
76  "--no-visualization",
77  "Enable or disable GLVis visualization.");
78  args.Parse();
79  if (!args.Good())
80  {
81  args.PrintUsage(cout);
82  return 1;
83  }
84  args.PrintOptions(cout);
85 
86  // 2. Enable hardware devices such as GPUs, and programming models such as
87  // CUDA, OCCA, RAJA and OpenMP based on command line options.
88  Device device(device_config);
89  device.Print();
90 
91  // 3. Read the mesh from the given mesh file. We can handle triangular,
92  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
93  // the same code.
94  Mesh mesh(mesh_file, 1, 1);
95  int dim = mesh.Dimension();
96  int sdim = mesh.SpaceDimension();
97 
98  // 4. Since a NURBS mesh can currently only be refined uniformly, we need to
99  // convert it to a piecewise-polynomial curved mesh. First we refine the
100  // NURBS mesh a bit more and then project the curvature to quadratic Nodes.
101  if (mesh.NURBSext)
102  {
103  for (int i = 0; i < 2; i++)
104  {
105  mesh.UniformRefinement();
106  }
107  mesh.SetCurvature(2);
108  }
109 
110  // 5. Define a finite element space on the mesh. The polynomial order is
111  // one (linear) by default, but this can be changed on the command line.
112  H1_FECollection fec(order, dim);
113  FiniteElementSpace fespace(&mesh, &fec);
114 
115  // 6. As in Example 1, we set up bilinear and linear forms corresponding to
116  // the Laplace problem -\Delta u = 1. We don't assemble the discrete
117  // problem yet, this will be done in the main loop.
118  BilinearForm a(&fespace);
119  if (pa)
120  {
121  a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
122  a.SetDiagonalPolicy(Operator::DIAG_ONE);
123  }
124  LinearForm b(&fespace);
125 
126  ConstantCoefficient one(1.0);
127  ConstantCoefficient zero(0.0);
128 
130  a.AddDomainIntegrator(integ);
132 
133  // 7. The solution vector x and the associated finite element grid function
134  // will be maintained over the AMR iterations. We initialize it to zero.
135  GridFunction x(&fespace);
136  x = 0.0;
137 
138  // 8. All boundary attributes will be used for essential (Dirichlet) BC.
139  MFEM_VERIFY(mesh.bdr_attributes.Size() > 0,
140  "Boundary attributes required in the mesh.");
141  Array<int> ess_bdr(mesh.bdr_attributes.Max());
142  ess_bdr = 1;
143 
144  // 9. Connect to GLVis.
145  char vishost[] = "localhost";
146  int visport = 19916;
147  socketstream sol_sock;
148  if (visualization)
149  {
150  sol_sock.open(vishost, visport);
151  }
152 
153  // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
154  // that uses the ComputeElementFlux method of the DiffusionIntegrator to
155  // recover a smoothed flux (gradient) that is subtracted from the element
156  // flux to get an error indicator. We need to supply the space for the
157  // smoothed flux: an (H1)^sdim (i.e., vector-valued) space is used here.
158  ErrorEstimator *estimator{nullptr};
159 
160  if (LSZZ)
161  {
162  estimator = new LSZienkiewiczZhuEstimator(*integ, x);
163  if (dim == 3 && mesh.GetElementType(0) != Element::HEXAHEDRON)
164  {
165  dynamic_cast<LSZienkiewiczZhuEstimator *>
166  (estimator)->SetTichonovRegularization();
167  }
168  }
169  else
170  {
171  auto flux_fes = new FiniteElementSpace(&mesh, &fec, sdim);
172  estimator = new ZienkiewiczZhuEstimator(*integ, x, flux_fes);
173  dynamic_cast<ZienkiewiczZhuEstimator *>(estimator)->SetAnisotropic();
174  }
175 
176  // 11. A refiner selects and refines elements based on a refinement strategy.
177  // The strategy here is to refine elements with errors larger than a
178  // fraction of the maximum element error. Other strategies are possible.
179  // The refiner will call the given error estimator.
180  ThresholdRefiner refiner(*estimator);
181  refiner.SetTotalErrorFraction(0.7);
182 
183  // 12. The main AMR loop. In each iteration we solve the problem on the
184  // current mesh, visualize the solution, and refine the mesh.
185  for (int it = 0; ; it++)
186  {
187  int cdofs = fespace.GetTrueVSize();
188  cout << "\nAMR iteration " << it << endl;
189  cout << "Number of unknowns: " << cdofs << endl;
190 
191  // 13. Assemble the right-hand side.
192  b.Assemble();
193 
194  // 14. Set Dirichlet boundary values in the GridFunction x.
195  // Determine the list of Dirichlet true DOFs in the linear system.
196  Array<int> ess_tdof_list;
197  x.ProjectBdrCoefficient(zero, ess_bdr);
198  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
199 
200  // 15. Assemble the stiffness matrix.
201  a.Assemble();
202 
203  // 16. Create the linear system: eliminate boundary conditions, constrain
204  // hanging nodes and possibly apply other transformations. The system
205  // will be solved for true (unconstrained) DOFs only.
206  OperatorPtr A;
207  Vector B, X;
208 
209  const int copy_interior = 1;
210  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
211 
212  // 17. Solve the linear system A X = B.
213  if (!pa)
214  {
215 #ifndef MFEM_USE_SUITESPARSE
216  // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
217  GSSmoother M((SparseMatrix&)(*A));
218  PCG(*A, M, B, X, 3, 200, 1e-12, 0.0);
219 #else
220  // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
221  UMFPackSolver umf_solver;
222  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
223  umf_solver.SetOperator(*A);
224  umf_solver.Mult(B, X);
225 #endif
226  }
227  else // Diagonal preconditioning in partial assembly mode.
228  {
229  OperatorJacobiSmoother M(a, ess_tdof_list);
230  PCG(*A, M, B, X, 3, 2000, 1e-12, 0.0);
231  }
232 
233  // 18. After solving the linear system, reconstruct the solution as a
234  // finite element GridFunction. Constrained nodes are interpolated
235  // from true DOFs (it may therefore happen that x.Size() >= X.Size()).
236  a.RecoverFEMSolution(X, b, x);
237 
238  // 19. Send solution by socket to the GLVis server.
239  if (visualization && sol_sock.good())
240  {
241  sol_sock.precision(8);
242  sol_sock << "solution\n" << mesh << x << flush;
243  }
244 
245  if (cdofs > max_dofs)
246  {
247  cout << "Reached the maximum number of dofs. Stop." << endl;
248  break;
249  }
250 
251  // 20. Call the refiner to modify the mesh. The refiner calls the error
252  // estimator to obtain element errors, then it selects elements to be
253  // refined and finally it modifies the mesh. The Stop() method can be
254  // used to determine if a stopping criterion was met.
255  refiner.Apply(mesh);
256  if (refiner.Stop())
257  {
258  cout << "Stopping criterion satisfied. Stop." << endl;
259  break;
260  }
261 
262  // 21. Update the space to reflect the new state of the mesh. Also,
263  // interpolate the solution x so that it lies in the new space but
264  // represents the same function. This saves solver iterations later
265  // since we'll have a good initial guess of x in the next step.
266  // Internally, FiniteElementSpace::Update() calculates an
267  // interpolation matrix which is then used by GridFunction::Update().
268  fespace.Update();
269  x.Update();
270 
271  // 22. Inform also the bilinear and linear forms that the space has
272  // changed.
273  a.Update();
274  b.Update();
275  }
276 
277  delete estimator;
278  return 0;
279 }
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:108
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:3334
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:84
Element::Type GetElementType(int i) const
Returns the type of element i.
Definition: mesh.cpp:6250
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:168
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.
The LSZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu error estimation procedure [1...
Definition: estimators.hpp:241
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:50
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:9816
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:5343
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.
Base class for all element based error estimators.
Definition: estimators.hpp:41
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:590
int Dimension() const
Definition: mesh.hpp:1006
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
int SpaceDimension() const
Definition: mesh.hpp:1007
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:35
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:96
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:165
void Update()
Update the object according to the associated FE space fes.
Definition: linearform.cpp:344
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
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
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:220
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
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:3189
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:449
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:3094
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150