MFEM  v3.3
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
16 // ex6 -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. Parse command-line options.
44  const char *mesh_file = "../data/star.mesh";
45  int order = 1;
46  bool visualization = 1;
47 
48  OptionsParser args(argc, argv);
49  args.AddOption(&mesh_file, "-m", "--mesh",
50  "Mesh file to use.");
51  args.AddOption(&order, "-o", "--order",
52  "Finite element order (polynomial degree).");
53  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
54  "--no-visualization",
55  "Enable or disable GLVis visualization.");
56  args.Parse();
57  if (!args.Good())
58  {
59  args.PrintUsage(cout);
60  return 1;
61  }
62  args.PrintOptions(cout);
63 
64  // 2. Read the mesh from the given mesh file. We can handle triangular,
65  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
66  // the same code.
67  Mesh mesh(mesh_file, 1, 1);
68  int dim = mesh.Dimension();
69  int sdim = mesh.SpaceDimension();
70 
71  // 3. Since a NURBS mesh can currently only be refined uniformly, we need to
72  // convert it to a piecewise-polynomial curved mesh. First we refine the
73  // NURBS mesh a bit more and then project the curvature to quadratic Nodes.
74  if (mesh.NURBSext)
75  {
76  for (int i = 0; i < 2; i++)
77  {
78  mesh.UniformRefinement();
79  }
80  mesh.SetCurvature(2);
81  }
82 
83  // 4. Define a finite element space on the mesh. The polynomial order is
84  // one (linear) by default, but this can be changed on the command line.
85  H1_FECollection fec(order, dim);
86  FiniteElementSpace fespace(&mesh, &fec);
87 
88  // 5. As in Example 1, we set up bilinear and linear forms corresponding to
89  // the Laplace problem -\Delta u = 1. We don't assemble the discrete
90  // problem yet, this will be done in the main loop.
91  BilinearForm a(&fespace);
92  LinearForm b(&fespace);
93 
94  ConstantCoefficient one(1.0);
95  ConstantCoefficient zero(0.0);
96 
98  a.AddDomainIntegrator(integ);
100 
101  // 6. The solution vector x and the associated finite element grid function
102  // will be maintained over the AMR iterations. We initialize it to zero.
103  GridFunction x(&fespace);
104  x = 0.0;
105 
106  // 7. All boundary attributes will be used for essential (Dirichlet) BC.
107  MFEM_VERIFY(mesh.bdr_attributes.Size() > 0,
108  "Boundary attributes required in the mesh.");
109  Array<int> ess_bdr(mesh.bdr_attributes.Max());
110  ess_bdr = 1;
111 
112  // 8. Connect to GLVis.
113  char vishost[] = "localhost";
114  int visport = 19916;
115  socketstream sol_sock;
116  if (visualization)
117  {
118  sol_sock.open(vishost, visport);
119  }
120 
121  // 9. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
122  // that uses the ComputeElementFlux method of the DiffusionIntegrator to
123  // recover a smoothed flux (gradient) that is subtracted from the element
124  // flux to get an error indicator. We need to supply the space for the
125  // smoothed flux: an (H1)^sdim (i.e., vector-valued) space is used here.
126  FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
127  ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
128  estimator.SetAnisotropic();
129 
130  // 10. A refiner selects and refines elements based on a refinement strategy.
131  // The strategy here is to refine elements with errors larger than a
132  // fraction of the maximum element error. Other strategies are possible.
133  // The refiner will call the given error estimator.
134  ThresholdRefiner refiner(estimator);
135  refiner.SetTotalErrorFraction(0.7);
136 
137  // 11. The main AMR loop. In each iteration we solve the problem on the
138  // current mesh, visualize the solution, and refine the mesh.
139  const int max_dofs = 50000;
140  for (int it = 0; ; it++)
141  {
142  int cdofs = fespace.GetTrueVSize();
143  cout << "\nAMR iteration " << it << endl;
144  cout << "Number of unknowns: " << cdofs << endl;
145 
146  // 12. Assemble the stiffness matrix and the right-hand side.
147  a.Assemble();
148  b.Assemble();
149 
150  // 13. Set Dirichlet boundary values in the GridFunction x.
151  // Determine the list of Dirichlet true DOFs in the linear system.
152  Array<int> ess_tdof_list;
153  x.ProjectBdrCoefficient(zero, ess_bdr);
154  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
155 
156  // 14. Create the linear system: eliminate boundary conditions, constrain
157  // hanging nodes and possibly apply other transformations. The system
158  // will be solved for true (unconstrained) DOFs only.
159  SparseMatrix A;
160  Vector B, X;
161  const int copy_interior = 1;
162  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
163 
164 #ifndef MFEM_USE_SUITESPARSE
165  // 15. Define a simple symmetric Gauss-Seidel preconditioner and use it to
166  // solve the linear system with PCG.
167  GSSmoother M(A);
168  PCG(A, M, B, X, 3, 200, 1e-12, 0.0);
169 #else
170  // 15. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the
171  // the linear system.
172  UMFPackSolver umf_solver;
173  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
174  umf_solver.SetOperator(A);
175  umf_solver.Mult(B, X);
176 #endif
177 
178  // 16. After solving the linear system, reconstruct the solution as a
179  // finite element GridFunction. Constrained nodes are interpolated
180  // from true DOFs (it may therefore happen that x.Size() >= X.Size()).
181  a.RecoverFEMSolution(X, b, x);
182 
183  // 17. Send solution by socket to the GLVis server.
184  if (visualization && sol_sock.good())
185  {
186  sol_sock.precision(8);
187  sol_sock << "solution\n" << mesh << x << flush;
188  }
189 
190  if (cdofs > max_dofs)
191  {
192  cout << "Reached the maximum number of dofs. Stop." << endl;
193  break;
194  }
195 
196  // 18. Call the refiner to modify the mesh. The refiner calls the error
197  // estimator to obtain element errors, then it selects elements to be
198  // refined and finally it modifies the mesh. The Stop() method can be
199  // used to determine if a stopping criterion was met.
200  refiner.Apply(mesh);
201  if (refiner.Stop())
202  {
203  cout << "Stopping criterion satisfied. Stop." << endl;
204  break;
205  }
206 
207  // 19. Update the space to reflect the new state of the mesh. Also,
208  // interpolate the solution x so that it lies in the new space but
209  // represents the same function. This saves solver iterations later
210  // since we'll have a good initial guess of x in the next step.
211  // Internally, FiniteElementSpace::Update() calculates an
212  // interpolation matrix which is then used by GridFunction::Update().
213  fespace.Update();
214  x.Update();
215 
216  // 20. Inform also the bilinear and linear forms that the space has
217  // changed.
218  a.Update();
219  b.Update();
220  }
221 
222  return 0;
223 }
int Size() const
Logical size of the array.
Definition: array.hpp:109
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
void SetAnisotropic(bool aniso=true)
Enable/disable anisotropic estimates. To enable this option, the BilinearFormIntegrator must support ...
Definition: estimators.hpp:139
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
virtual void Update(bool want_transform=true)
Definition: fespace.cpp:1457
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
Subclass constant coefficient.
Definition: coefficient.hpp:57
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:42
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list)
Definition: fespace.cpp:295
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
STL namespace.
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:337
bool Apply(Mesh &mesh)
Perform the mesh operation.
int dim
Definition: ex3.cpp:47
Data type sparse matrix.
Definition: sparsemat.hpp:38
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6684
Mesh refinement operator using an error threshold.
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3278
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:108
virtual void Update(FiniteElementSpace *nfes=NULL)
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:456
int Dimension() const
Definition: mesh.hpp:611
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
int SpaceDimension() const
Definition: mesh.hpp:612
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator.
Definition: linearform.cpp:19
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:22
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:142
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:348
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 Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:144
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:144
virtual int GetTrueVSize()
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:166
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
int main(int argc, char *argv[])
Definition: ex6.cpp:41
int open(const char hostname[], int port)
The ZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu error estimation procedure...
Definition: estimators.hpp:72
Vector data type.
Definition: vector.hpp:36
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:146
Class for linear form - Vector with associated FE space and LFIntegrators.
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...
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, SparseMatrix &A, Vector &X, Vector &B, int copy_interior=0)
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1774
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
Definition: gridfunc.hpp:197
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:1674
bool Good() const
Definition: optparser.hpp:120