MFEM  v3.0
 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 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
16 //
17 // Description: This is a version of Example 1 with a simple adaptive mesh
18 // refinement loop. The problem being solved is again the Laplace
19 // equation -Delta u = 1 with homogeneous Dirichlet boundary
20 // conditions. The problem is solved on a sequence of meshes which
21 // are locally refined in a conforming (triangles, tetrahedrons)
22 // or non-conforming (quadrilateral, hexahedrons) manner according
23 // to a simple ZZ error estimator.
24 //
25 // The example demonstrates MFEM's capability to work with both
26 // conforming and nonconforming refinements, in 2D and 3D, on
27 // linear, curved and surface meshes. Interpolation of functions
28 // from coarse to fine meshes, as well as persistent GLVis
29 // visualization are also illustrated.
30 //
31 // We recommend viewing Example 1 before viewing this example.
32 
33 #include "mfem.hpp"
34 #include <fstream>
35 #include <iostream>
36 
37 using namespace std;
38 using namespace mfem;
39 
40 int main(int argc, char *argv[])
41 {
42  // 1. Parse command-line options.
43  const char *mesh_file = "../data/star.mesh";
44  int order = 1;
45  bool visualization = 1;
46 
47  OptionsParser args(argc, argv);
48  args.AddOption(&mesh_file, "-m", "--mesh",
49  "Mesh file to use.");
50  args.AddOption(&order, "-o", "--order",
51  "Finite element order (polynomial degree).");
52  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
53  "--no-visualization",
54  "Enable or disable GLVis visualization.");
55  args.Parse();
56  if (!args.Good())
57  {
58  args.PrintUsage(cout);
59  return 1;
60  }
61  args.PrintOptions(cout);
62 
63  // 2. Read the mesh from the given mesh file. We can handle triangular,
64  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
65  // the same code.
66  ifstream imesh(mesh_file);
67  if (!imesh)
68  {
69  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
70  return 2;
71  }
72  Mesh mesh(imesh, 1, 1);
73  imesh.close();
74  int dim = mesh.Dimension();
75 
76  // 3. Since a NURBS mesh can currently only be refined uniformly, we need to
77  // convert it to a piecewise-polynomial curved mesh. First we refine the
78  // NURBS mesh a bit and then project the curvature to quadratic Nodes.
79  if (mesh.NURBSext)
80  {
81  for (int i = 0; i < 2; i++)
82  mesh.UniformRefinement();
83 
84  FiniteElementCollection* nfec = new H1_FECollection(2, dim);
85  FiniteElementSpace* nfes = new FiniteElementSpace(&mesh, nfec, dim);
86  mesh.SetNodalFESpace(nfes);
87  mesh.GetNodes()->MakeOwner(nfec);
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  Array<int> ess_bdr(mesh.bdr_attributes.Max());
114  ess_bdr = 1;
115 
116  // 8. Connect to GLVis.
117  char vishost[] = "localhost";
118  int visport = 19916;
119  socketstream sol_sock;
120  if (visualization)
121  sol_sock.open(vishost, visport);
122 
123  // 9. The main AMR loop. In each iteration we solve the problem on the
124  // current mesh, visualize the solution, estimate the error on all
125  // elements, refine the worst elements and update all objects to work
126  // with the new mesh.
127  const int max_it = 15;
128  for (int it = 0; it < max_it; it++)
129  {
130  cout << "\nIteration " << it << endl;
131  cout << "Number of unknowns: " << fespace.GetNConformingDofs() << endl;
132 
133  // 10. Assemble the stiffness matrix and the right-hand side. Note that
134  // MFEM doesn't care at this point if the mesh is nonconforming (i.e.,
135  // contains hanging nodes). The FE space is considered 'cut' along
136  // hanging edges/faces.
137  a.Assemble();
138  b.Assemble();
139 
140  x.ProjectBdrCoefficient(zero, ess_bdr);
141 
142  // 11. Take care of nonconforming meshes by applying the interpolation
143  // matrix P to a, b and x, so that slave degrees of freedom get
144  // eliminated from the linear system. The system becomes P'AP x = P'b.
145  // (If the mesh is conforming, P is identity.)
146  a.ConformingAssemble(x, b);
147 
148  // 12. As usual, we also need to eliminate the essential BC from the
149  // system. This needs to be done after ConformingAssemble.
150  a.EliminateEssentialBC(ess_bdr, x, b);
151 
152  const SparseMatrix &A = a.SpMat();
153 #ifndef MFEM_USE_SUITESPARSE
154  // 13. Define a simple symmetric Gauss-Seidel preconditioner and use it to
155  // solve the linear system with PCG.
156  GSSmoother M(A);
157  PCG(A, M, b, x, 1, 200, 1e-12, 0.0);
158 #else
159  // 13. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the
160  // the linear system.
161  UMFPackSolver umf_solver;
162  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
163  umf_solver.SetOperator(A);
164  umf_solver.Mult(b, x);
165 #endif
166 
167  // 14. For nonconforming meshes, bring the solution vector back from
168  // the conforming space to the nonconforming (cut) space, i.e.,
169  // x = Px. Slave DOFs receive the correct values to make the solution
170  // continuous.
172 
173  // 15. Send solution by socket to the GLVis server.
174  if (visualization && sol_sock.good())
175  {
176  sol_sock.precision(8);
177  sol_sock << "solution\n" << mesh << x << flush;
178  }
179 
180  // 16. Estimate element errors using the Zienkiewicz-Zhu error estimator.
181  // The bilinear form integrator must have the 'ComputeElementFlux'
182  // method defined.
183  Vector errors(mesh.GetNE());
184  {
185  FiniteElementSpace flux_fespace(&mesh, &fec, dim);
186  DiffusionIntegrator flux_integrator(one);
187  GridFunction flux(&flux_fespace);
188  ComputeFlux(flux_integrator, x, flux);
189  ZZErrorEstimator(flux_integrator, x, flux, errors, 1);
190  }
191 
192  // 17. Make a list of elements whose error is larger than a fraction (0.7)
193  // of the maximum element error. These elements will be refined.
194  Array<int> ref_list;
195  const double frac = 0.7;
196  // the 'errors' are squared, so we need to square the fraction
197  double threshold = (frac*frac) * errors.Max();
198  for (int i = 0; i < errors.Size(); i++)
199  if (errors[i] >= threshold)
200  ref_list.Append(i);
201 
202  // 18. Refine the selected elements. Since we are going to transfer the
203  // grid function x from the coarse mesh to the new fine mesh in the
204  // next step, we need to request the "two-level state" of the mesh.
205  mesh.UseTwoLevelState(1);
206  mesh.GeneralRefinement(ref_list);
207 
208  // 19. Update the space to reflect the new state of the mesh. Also,
209  // interpolate the solution x so that it lies in the new space but
210  // represents the same function. This saves solver iterations since
211  // we'll have a good initial guess of x in the next step.
212  // The interpolation algorithm needs the mesh to hold some information
213  // about the previous state, which is why the call UseTwoLevelState
214  // above is required.
215  fespace.UpdateAndInterpolate(&x);
216 
217  // Note: If interpolation was not needed, we could just use the following
218  // two calls to update the space and the grid function. (No need to
219  // call UseTwoLevelState in this case.)
220  // fespace.Update();
221  // x.Update();
222 
223  // 20. Inform also the bilinear and linear forms that the space has
224  // changed.
225  a.Update();
226  b.Update();
227  }
228 
229  return 0;
230 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:26
void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &u, GridFunction &flux, int wcoef, int sd)
Definition: gridfunc.cpp:1984
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 ConformingProlongate(const Vector &x)
Definition: gridfunc.cpp:1750
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:34
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:396
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:329
void ZZErrorEstimator(BilinearFormIntegrator &blfi, GridFunction &u, GridFunction &flux, Vector &ErrorEstimates, int wsd)
Definition: gridfunc.cpp:2041
Data type sparse matrix.
Definition: sparsemat.hpp:38
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:330
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6225
T Max() const
Definition: array.cpp:78
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:424
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3066
int Dimension() const
Definition: mesh.hpp:417
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:385
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator.
Definition: linearform.cpp:19
Array< int > bdr_attributes
Definition: mesh.hpp:305
void UseTwoLevelState(int use)
Definition: mesh.hpp:684
int main(int argc, char *argv[])
Definition: ex1.cpp:39
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:340
Abstract finite element space.
Definition: fespace.hpp:61
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:307
virtual void UpdateAndInterpolate(int num_grid_fns,...)
Definition: fespace.cpp:1080
int GetNConformingDofs() const
Definition: fespace.hpp:155
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:266
int open(const char hostname[], int port)
Vector data type.
Definition: vector.hpp:29
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:4948
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:52
Class for linear form - Vector with associated FE space and LFIntegrators.
Definition: linearform.hpp:23
const SparseMatrix & SpMat() const
Returns a reference to the sparse martix.
void EliminateEssentialBC(Array< int > &bdr_attr_is_ess, Vector &sol, Vector &rhs, int d=0)
virtual void Mult(const Vector &b, Vector &x) const
Operator application.
Definition: solvers.cpp:1597
void GeneralRefinement(Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:5894
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
Definition: gridfunc.hpp:146
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:1505
bool Good() const
Definition: optparser.hpp:120