MFEM  v3.1 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 1
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 //
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);
50  "Mesh file to use.");
52  "Finite element order (polynomial degree).");
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  ifstream imesh(mesh_file);
68  if (!imesh)
69  {
70  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
71  return 2;
72  }
73  Mesh mesh(imesh, 1, 1);
74  imesh.close();
75  int dim = mesh.Dimension();
76  int sdim = mesh.SpaceDimension();
77
78  // 3. Since a NURBS mesh can currently only be refined uniformly, we need to
79  // convert it to a piecewise-polynomial curved mesh. First we refine the
80  // NURBS mesh a bit more and then project the curvature to quadratic Nodes.
81  if (mesh.NURBSext)
82  {
83  for (int i = 0; i < 2; i++)
84  {
85  mesh.UniformRefinement();
86  }
87  mesh.SetCurvature(2);
88  }
89
90  // 4. Define a finite element space on the mesh. The polynomial order is
91  // one (linear) by default, but this can be changed on the command line.
92  H1_FECollection fec(order, dim);
93  FiniteElementSpace fespace(&mesh, &fec);
94
95  // 5. As in Example 1, we set up bilinear and linear forms corresponding to
96  // the Laplace problem -\Delta u = 1. We don't assemble the discrete
97  // problem yet, this will be done in the main loop.
98  BilinearForm a(&fespace);
99  LinearForm b(&fespace);
100
101  ConstantCoefficient one(1.0);
102  ConstantCoefficient zero(0.0);
103
106
107  // 6. The solution vector x and the associated finite element grid function
108  // will be maintained over the AMR iterations. We initialize it to zero.
109  GridFunction x(&fespace);
110  x = 0;
111
112  // 7. All boundary attributes will be used for essential (Dirichlet) BC.
113  MFEM_VERIFY(mesh.bdr_attributes.Size() > 0,
114  "Boundary attributes required in the mesh.");
115  Array<int> ess_bdr(mesh.bdr_attributes.Max());
116  ess_bdr = 1;
117
118  // 8. Connect to GLVis.
119  char vishost[] = "localhost";
120  int visport = 19916;
121  socketstream sol_sock;
122  if (visualization)
123  {
124  sol_sock.open(vishost, visport);
125  }
126
127  // 9. The main AMR loop. In each iteration we solve the problem on the
128  // current mesh, visualize the solution, estimate the error on all
129  // elements, refine the worst elements and update all objects to work
130  // with the new mesh.
131  const int max_dofs = 50000;
132  for (int it = 0; ; it++)
133  {
134  int cdofs = fespace.GetTrueVSize();
135  cout << "\nIteration " << it << endl;
136  cout << "Number of unknowns: " << cdofs << endl;
137
138  // 10. Assemble the stiffness matrix and the right-hand side. Note that
139  // MFEM doesn't care at this point if the mesh is nonconforming (i.e.,
140  // contains hanging nodes). The FE space is considered 'cut' along
141  // hanging edges/faces.
142  a.Assemble();
143  b.Assemble();
144
145  // 11. Set Dirichlet boundary values in the GridFunction x.
146  // Determine the list of Dirichlet true DOFs in the linear system.
147  Array<int> ess_tdof_list;
148  x.ProjectBdrCoefficient(zero, ess_bdr);
149  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
150
151  // 12. Create the linear system: eliminate boundary conditions, constrain
152  // hanging nodes and possibly apply other transformations. The system
153  // will be solved for true (unconstrained) DOFs only.
154  SparseMatrix A;
155  Vector B, X;
156  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, 1);
157
158 #ifndef MFEM_USE_SUITESPARSE
159  // 13. Define a simple symmetric Gauss-Seidel preconditioner and use it to
160  // solve the linear system with PCG.
161  GSSmoother M(A);
162  PCG(A, M, B, X, 2, 200, 1e-12, 0.0);
163 #else
164  // 13. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the
165  // the linear system.
166  UMFPackSolver umf_solver;
167  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
168  umf_solver.SetOperator(A);
169  umf_solver.Mult(B, X);
170 #endif
171
172  // 14. After solving the linear system, reconstruct the solution as a finite
173  // element grid function. Constrained nodes are interpolated from true
174  // DOFs (it may therefore happen that dim(x) >= dim(X)).
175  a.RecoverFEMSolution(X, b, x);
176
177  // 15. Send solution by socket to the GLVis server.
178  if (visualization && sol_sock.good())
179  {
180  sol_sock.precision(8);
181  sol_sock << "solution\n" << mesh << x << flush;
182  }
183
184  if (cdofs > max_dofs)
185  {
186  break;
187  }
188
189  // 16. Estimate element errors using the Zienkiewicz-Zhu error estimator.
190  // The bilinear form integrator must have the 'ComputeElementFlux'
191  // method defined.
192  Vector errors(mesh.GetNE());
193  Array<int> aniso_flags;
194  {
195  DiffusionIntegrator flux_integrator(one);
196  FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
197  GridFunction flux(&flux_fespace);
198  ZZErrorEstimator(flux_integrator, x, flux, errors, &aniso_flags);
199  }
200
201  // 17. Make a list of elements whose error is larger than a fraction (0.7)
202  // of the maximum element error. These elements will be refined.
203  Array<Refinement> ref_list;
204  const double frac = 0.7;
205  // the 'errors' are squared, so we need to square the fraction
206  double threshold = (frac*frac) * errors.Max();
207  for (int i = 0; i < errors.Size(); i++)
208  {
209  if (errors[i] >= threshold)
210  {
211  ref_list.Append(Refinement(i, aniso_flags[i]));
212  }
213  }
214
215  // 18. Refine the selected elements. Since we are going to transfer the
216  // grid function x from the coarse mesh to the new fine mesh in the
217  // next step, we need to request the "two-level state" of the mesh.
218  mesh.UseTwoLevelState(1);
219  mesh.GeneralRefinement(ref_list);
220
221  // 19. Update the space to reflect the new state of the mesh. Also,
222  // interpolate the solution x so that it lies in the new space but
223  // represents the same function. This saves solver iterations since
224  // we'll have a good initial guess of x in the next step.
225  // The interpolation algorithm needs the mesh to hold some information
226  // about the previous state, which is why the call UseTwoLevelState
227  // above is required.
228  fespace.UpdateAndInterpolate(&x);
229
230  // Note: If interpolation was not needed, we could just use the following
231  // two calls to update the space and the grid function. (No need to
232  // call UseTwoLevelState in this case.)
233  // fespace.Update();
234  // x.Update();
235
236  // 20. Inform also the bilinear and linear forms that the space has
237  // changed.
238  a.Update();
239  b.Update();
240  }
241
242  return 0;
243 }
int Size() const
Logical size of the array.
Definition: array.hpp:109
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
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:34
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list)
Definition: fespace.cpp:469
void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:454
void FormLinearSystem(Array< int > &ess_tdof_list, Vector &x, Vector &b, SparseMatrix &A, Vector &X, Vector &B, int copy_interior=0)
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:332
int dim
Definition: ex3.cpp:48
Data type sparse matrix.
Definition: sparsemat.hpp:38
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:368
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7316
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3736
T Max() const
Definition: array.cpp:90
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:440
int Dimension() const
Definition: mesh.hpp:475
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void ZZErrorEstimator(BilinearFormIntegrator &blfi, GridFunction &u, GridFunction &flux, Vector &error_estimates, Array< int > *aniso_flags, int with_subdomains)
Definition: gridfunc.cpp:2350
int SpaceDimension() const
Definition: mesh.hpp:476
Definition: linearform.cpp:19
Array< int > bdr_attributes
Definition: mesh.hpp:140
void UseTwoLevelState(int use)
Definition: mesh.hpp:766
int main(int argc, char *argv[])
Definition: ex1.cpp:45
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:343
Abstract finite element space.
Definition: fespace.hpp:62
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
NURBSExtension * NURBSext
Definition: mesh.hpp:142
virtual int GetTrueVSize()
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:167
virtual void UpdateAndInterpolate(int num_grid_fns,...)
Definition: fespace.cpp:1531
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
int open(const char hostname[], int port)
Vector data type.
Definition: vector.hpp:33
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:54
Class for linear form - Vector with associated FE space and LFIntegrators.
Definition: linearform.hpp:23
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:6864
virtual void Mult(const Vector &b, Vector &x) const
Operator application.
Definition: solvers.cpp:1725
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
Definition: gridfunc.hpp:170
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:1625
bool Good() const
Definition: optparser.hpp:120