MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex15.cpp
Go to the documentation of this file.
1 // MFEM Example 15
2 //
3 // Compile with: make ex15
4 //
5 // Sample runs: ex15
6 // ex15 -o 1 -y 0.4
7 // ex15 -o 4 -y 0.1
8 // ex15 -n 5
9 // ex15 -p 1 -n 3
10 //
11 // Other meshes:
12 //
13 // ex15 -m ../data/square-disc-nurbs.mesh
14 // ex15 -m ../data/disc-nurbs.mesh
15 // ex15 -m ../data/fichera.mesh -tf 0.3
16 // ex15 -m ../data/ball-nurbs.mesh -tf 0.3
17 // ex15 -m ../data/mobius-strip.mesh
18 // ex15 -m ../data/amr-quad.mesh
19 // ex15 -m ../data/square-disc.mesh
20 // ex15 -m ../data/escher.mesh -r 2 -tf 0.3
21 //
22 // Description: Building on Example 6, this example demonstrates dynamic AMR.
23 // The mesh is adapted to a time-dependent solution by refinement
24 // as well as by derefinement. For simplicity, the solution is
25 // prescribed and no time integration is done. However, the error
26 // estimation and refinement/derefinement decisions are realistic.
27 //
28 // At each outer iteration the right hand side function is changed
29 // to mimic a time dependent problem. Within each inner iteration
30 // the problem is solved on a sequence of meshes which are locally
31 // refined according to a simple ZZ error estimator. At the end
32 // of the inner iteration the error estimates are also used to
33 // identify any elements which may be over-refined and a single
34 // derefinement step is performed.
35 //
36 // The example demonstrates MFEM's capability to refine and
37 // derefine nonconforming meshes, in 2D and 3D, and on linear,
38 // curved and surface meshes. Interpolation of functions between
39 // coarse and fine meshes, persistent GLVis visualization, and
40 // saving of time-dependent fields for external visualization with
41 // VisIt (visit.llnl.gov) are also illustrated.
42 //
43 // We recommend viewing Examples 1, 6 and 9 before viewing this
44 // example.
45 
46 #include "mfem.hpp"
47 #include <fstream>
48 #include <iostream>
49 
50 using namespace std;
51 using namespace mfem;
52 
53 // Choices for the problem setup. Affect bdr_func and rhs_func.
54 int problem;
56 
57 // Prescribed time-dependent boundary and right-hand side functions.
58 double bdr_func(const Vector &pt, double t);
59 double rhs_func(const Vector &pt, double t);
60 
61 // Update the finite element space, interpolate the solution and perform
62 // parallel load balancing.
63 void UpdateProblem(Mesh &mesh, FiniteElementSpace &fespace,
65 
66 
67 int main(int argc, char *argv[])
68 {
69  // 1. Parse command-line options.
70  problem = 0;
71  nfeatures = 1;
72  const char *mesh_file = "../data/star-hilbert.mesh";
73  int order = 2;
74  double t_final = 1.0;
75  double max_elem_error = 5.0e-3;
76  double hysteresis = 0.15; // derefinement safety coefficient
77  int ref_levels = 0;
78  int nc_limit = 3; // maximum level of hanging nodes
79  bool visualization = true;
80  bool visit = false;
81 
82  OptionsParser args(argc, argv);
83  args.AddOption(&mesh_file, "-m", "--mesh",
84  "Mesh file to use.");
85  args.AddOption(&problem, "-p", "--problem",
86  "Problem setup to use: 0 = spherical front, 1 = ball.");
87  args.AddOption(&nfeatures, "-n", "--nfeatures",
88  "Number of solution features (fronts/balls).");
89  args.AddOption(&order, "-o", "--order",
90  "Finite element order (polynomial degree).");
91  args.AddOption(&max_elem_error, "-e", "--max-err",
92  "Maximum element error");
93  args.AddOption(&hysteresis, "-y", "--hysteresis",
94  "Derefinement safety coefficient.");
95  args.AddOption(&ref_levels, "-r", "--ref-levels",
96  "Number of initial uniform refinement levels.");
97  args.AddOption(&nc_limit, "-l", "--nc-limit",
98  "Maximum level of hanging nodes.");
99  args.AddOption(&t_final, "-tf", "--t-final",
100  "Final time; start time is 0.");
101  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
102  "--no-visualization",
103  "Enable or disable GLVis visualization.");
104  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
105  "--no-visit-datafiles",
106  "Save data files for VisIt (visit.llnl.gov) visualization.");
107  args.Parse();
108  if (!args.Good())
109  {
110  args.PrintUsage(cout);
111  return 1;
112  }
113  args.PrintOptions(cout);
114 
115  // 2. Read the mesh from the given mesh file on all processors. We can handle
116  // triangular, quadrilateral, tetrahedral, hexahedral, surface and volume
117  // meshes with the same code.
118  Mesh mesh(mesh_file, 1, 1);
119  int dim = mesh.Dimension();
120  int sdim = mesh.SpaceDimension();
121 
122  // 3. Project a NURBS mesh to a piecewise-quadratic curved mesh. Make sure
123  // that the mesh is non-conforming if it has quads or hexes and refine it.
124  if (mesh.NURBSext)
125  {
126  mesh.UniformRefinement();
127  if (ref_levels > 0) { ref_levels--; }
128  mesh.SetCurvature(2);
129  }
130  mesh.EnsureNCMesh(true);
131  for (int l = 0; l < ref_levels; l++)
132  {
133  mesh.UniformRefinement();
134  }
135  // Make sure tet-only meshes are marked for local refinement.
136  mesh.Finalize(true);
137 
138  // 4. 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  // 5. Define a finite element space on the mesh. The polynomial order is one
145  // (linear) by default, but this can be changed on the command line.
146  H1_FECollection fec(order, dim);
147  FiniteElementSpace fespace(&mesh, &fec);
148 
149  // 6. As in Example 1p, we set up bilinear and linear forms corresponding to
150  // the Laplace problem -\Delta u = 1. We don't assemble the discrete
151  // problem yet, this will be done in the inner loop.
152  BilinearForm a(&fespace);
153  LinearForm b(&fespace);
154 
155  ConstantCoefficient one(1.0);
158 
160  a.AddDomainIntegrator(integ);
162 
163  // 7. The solution vector x and the associated finite element grid function
164  // will be maintained over the AMR iterations.
165  GridFunction x(&fespace);
166 
167  // 8. Connect to GLVis. Prepare for VisIt output.
168  char vishost[] = "localhost";
169  int visport = 19916;
170 
171  socketstream sout;
172  if (visualization)
173  {
174  sout.open(vishost, visport);
175  if (!sout)
176  {
177  cout << "Unable to connect to GLVis server at "
178  << vishost << ':' << visport << endl;
179  cout << "GLVis visualization disabled.\n";
180  visualization = false;
181  }
182  sout.precision(8);
183  }
184 
185  VisItDataCollection visit_dc("Example15", &mesh);
186  visit_dc.RegisterField("solution", &x);
187  int vis_cycle = 0;
188 
189  // 9. As in Example 6, we set up a Zienkiewicz-Zhu estimator that will be
190  // used to obtain element error indicators. The integrator needs to
191  // provide the method ComputeElementFlux. The smoothed flux space is a
192  // vector valued H1 space here.
193  FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
194  ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
195 
196  // 10. As in Example 6, we also need a refiner. This time the refinement
197  // strategy is based on a fixed threshold that is applied locally to each
198  // element. The global threshold is turned off by setting the total error
199  // fraction to zero. We also enforce a maximum refinement ratio between
200  // adjacent elements.
201  ThresholdRefiner refiner(estimator);
202  refiner.SetTotalErrorFraction(0.0); // use purely local threshold
203  refiner.SetLocalErrorGoal(max_elem_error);
204  refiner.PreferConformingRefinement();
205  refiner.SetNCLimit(nc_limit);
206 
207  // 11. A derefiner selects groups of elements that can be coarsened to form
208  // a larger element. A conservative enough threshold needs to be set to
209  // prevent derefining elements that would immediately be refined again.
210  ThresholdDerefiner derefiner(estimator);
211  derefiner.SetThreshold(hysteresis * max_elem_error);
212  derefiner.SetNCLimit(nc_limit);
213 
214  // 12. The outer time loop. In each iteration we update the right hand side,
215  // solve the problem on the current mesh, visualize the solution and
216  // refine the mesh as many times as necessary. Then we derefine any
217  // elements which have very small errors.
218  x = 0.0;
219  for (double time = 0.0; time < t_final + 1e-10; time += 0.01)
220  {
221  cout << "\nTime " << time << "\n\nRefinement:" << endl;
222 
223  // Set the current time in the coefficients.
224  bdr.SetTime(time);
225  rhs.SetTime(time);
226 
227  // Make sure errors will be recomputed in the following.
228  refiner.Reset();
229  derefiner.Reset();
230 
231  // 13. The inner refinement loop. At the end we want to have the current
232  // time step resolved to the prescribed tolerance in each element.
233  for (int ref_it = 1; ; ref_it++)
234  {
235  cout << "Iteration: " << ref_it << ", number of unknowns: "
236  << fespace.GetVSize() << endl;
237 
238  // 14. Recompute the field on the current mesh: assemble the stiffness
239  // matrix and the right-hand side.
240  a.Assemble();
241  b.Assemble();
242 
243  // 15. Project the exact solution to the essential boundary DOFs.
244  x.ProjectBdrCoefficient(bdr, ess_bdr);
245 
246  // 16. Create and solve the linear system.
247  Array<int> ess_tdof_list;
248  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
249 
250  SparseMatrix A;
251  Vector B, X;
252  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
253 
254 #ifndef MFEM_USE_SUITESPARSE
255  GSSmoother M(A);
256  PCG(A, M, B, X, 0, 500, 1e-12, 0.0);
257 #else
258  UMFPackSolver umf_solver;
259  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
260  umf_solver.SetOperator(A);
261  umf_solver.Mult(B, X);
262 #endif
263 
264  // 17. Extract the local solution on each processor.
265  a.RecoverFEMSolution(X, b, x);
266 
267  // 18. Send the solution by socket to a GLVis server and optionally
268  // save it in VisIt format.
269  if (visualization)
270  {
271  sout.precision(8);
272  sout << "solution\n" << mesh << x << flush;
273  }
274  if (visit)
275  {
276  visit_dc.SetCycle(vis_cycle++);
277  visit_dc.SetTime(time);
278  visit_dc.Save();
279  }
280 
281  // 19. Apply the refiner on the mesh. The refiner calls the error
282  // estimator to obtain element errors, then it selects elements to
283  // be refined and finally it modifies the mesh. The Stop() method
284  // determines if all elements satisfy the local threshold.
285  refiner.Apply(mesh);
286  if (refiner.Stop())
287  {
288  break;
289  }
290 
291  // 20. Update the space and interpolate the solution.
292  UpdateProblem(mesh, fespace, x, a, b);
293  }
294 
295  // 21. Use error estimates from the last inner iteration to check for
296  // possible derefinements. The derefiner works similarly as the
297  // refiner. The errors are not recomputed because the mesh did not
298  // change (and also the estimator was not Reset() at this time).
299  if (derefiner.Apply(mesh))
300  {
301  cout << "\nDerefined elements." << endl;
302 
303  // 22. Update the space and interpolate the solution.
304  UpdateProblem(mesh, fespace, x, a, b);
305  }
306 
307  a.Update();
308  b.Update();
309  }
310 
311  return 0;
312 }
313 
314 
315 void UpdateProblem(Mesh &mesh, FiniteElementSpace &fespace,
317 {
318  // Update the space: recalculate the number of DOFs and construct a matrix
319  // that will adjust any GridFunctions to the new mesh state.
320  fespace.Update();
321 
322  // Interpolate the solution on the new mesh by applying the transformation
323  // matrix computed in the finite element space. Multiple GridFunctions could
324  // be updated here.
325  x.Update();
326 
327  // Free any transformation matrices to save memory.
328  fespace.UpdatesFinished();
329 
330  // Inform the linear and bilinear forms that the space has changed.
331  a.Update();
332  b.Update();
333 }
334 
335 
336 const double alpha = 0.02;
337 
338 // Spherical front with a Gaussian cross section and radius t
339 double front(double x, double y, double z, double t, int)
340 {
341  double r = sqrt(x*x + y*y + z*z);
342  return exp(-0.5*pow((r - t)/alpha, 2));
343 }
344 
345 double front_laplace(double x, double y, double z, double t, int dim)
346 {
347  double x2 = x*x, y2 = y*y, z2 = z*z, t2 = t*t;
348  double r = sqrt(x2 + y2 + z2);
349  double a2 = alpha*alpha, a4 = a2*a2;
350  return -exp(-0.5*pow((r - t)/alpha, 2)) / a4 *
351  (-2*t*(x2 + y2 + z2 - (dim-1)*a2/2)/r + x2 + y2 + z2 + t2 - dim*a2);
352 }
353 
354 // Smooth spherical step function with radius t
355 double ball(double x, double y, double z, double t, int)
356 {
357  double r = sqrt(x*x + y*y + z*z);
358  return -atan(2*(r - t)/alpha);
359 }
360 
361 double ball_laplace(double x, double y, double z, double t, int dim)
362 {
363  double x2 = x*x, y2 = y*y, z2 = z*z, t2 = 4*t*t;
364  double r = sqrt(x2 + y2 + z2);
365  double a2 = alpha*alpha;
366  double den = pow(-a2 - 4*(x2 + y2 + z2 - 2*r*t) - t2, 2.0);
367  return (dim == 2) ? 2*alpha*(a2 + t2 - 4*x2 - 4*y2)/r/den
368  /* */ : 4*alpha*(a2 + t2 - 4*r*t)/r/den;
369 }
370 
371 // Composes several features into one function
372 template<typename F0, typename F1>
373 double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
374 {
375  int dim = pt.Size();
376  double x = pt(0), y = pt(1), z = 0.0;
377  if (dim == 3) { z = pt(2); }
378 
379  if (problem == 0)
380  {
381  if (nfeatures <= 1)
382  {
383  return f0(x, y, z, t, dim);
384  }
385  else
386  {
387  double sum = 0.0;
388  for (int i = 0; i < nfeatures; i++)
389  {
390  double x0 = 0.5*cos(2*M_PI * i / nfeatures);
391  double y0 = 0.5*sin(2*M_PI * i / nfeatures);
392  sum += f0(x - x0, y - y0, z, t, dim);
393  }
394  return sum;
395  }
396  }
397  else
398  {
399  double sum = 0.0;
400  for (int i = 0; i < nfeatures; i++)
401  {
402  double x0 = 0.5*cos(2*M_PI * i / nfeatures + M_PI*t);
403  double y0 = 0.5*sin(2*M_PI * i / nfeatures + M_PI*t);
404  sum += f1(x - x0, y - y0, z, 0.25, dim);
405  }
406  return sum;
407  }
408 }
409 
410 // Exact solution, used for the Dirichlet BC.
411 double bdr_func(const Vector &pt, double t)
412 {
413  return composite_func(pt, t, front, ball);
414 }
415 
416 // Laplace of the exact solution, used for the right hand side.
417 double rhs_func(const Vector &pt, double t)
418 {
419  return composite_func(pt, t, front_laplace, ball_laplace);
420 }
421 
double front_laplace(double x, double y, double z, double t, int dim)
Definition: ex15.cpp:345
int Size() const
Logical size of the array.
Definition: array.hpp:124
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:93
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:381
double bdr_func(const Vector &pt, double t)
Definition: ex15.cpp:411
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
virtual void Update(bool want_transform=true)
Definition: fespace.cpp:2057
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
Subclass constant coefficient.
Definition: coefficient.hpp:67
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:79
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().
void SetNCLimit(int nc_limit)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited).
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: fespace.cpp:366
virtual void Save()
Save the collection and a VisIt root file.
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
int main(int argc, char *argv[])
Definition: ex1.cpp:62
double rhs_func(const Vector &pt, double t)
Definition: ex15.cpp:417
void UpdateProblem(Mesh &mesh, FiniteElementSpace &fespace, GridFunction &x, BilinearForm &a, LinearForm &b)
Definition: ex15.cpp:315
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:580
bool Apply(Mesh &mesh)
Perform the mesh operation.
virtual void Reset()
Reset the associated estimator.
Data type sparse matrix.
Definition: sparsemat.hpp:40
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:7645
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7982
virtual void Reset()
Reset the associated estimator.
Mesh refinement operator using an error threshold.
Data collection with VisIt I/O routines.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3924
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
int nfeatures
Definition: ex15.cpp:55
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:529
int Dimension() const
Definition: mesh.hpp:744
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetTime(double t)
Set physical time (for time-dependent simulations)
void SetTime(double t)
Definition: coefficient.hpp:39
void SetNCLimit(int nc_limit)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited).
int SpaceDimension() const
Definition: mesh.hpp:745
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:24
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:189
double ball_laplace(double x, double y, double z, double t, int dim)
Definition: ex15.cpp:361
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:591
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
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:76
int problem
Definition: ex15.cpp:54
void PreferConformingRefinement()
Use conforming refinement, if possible (triangles, tetrahedra) – this is the default.
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:152
void Update()
Update the object according to the associated FE space fes.
Definition: linearform.hpp:166
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:191
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
Definition: mesh.cpp:2316
int dim
Definition: ex24.cpp:43
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
double ball(double x, double y, double z, double t, int)
Definition: ex15.cpp:355
int open(const char hostname[], int port)
void SetThreshold(double thresh)
Set the de-refinement threshold. The default value is zero.
The ZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu error estimation procedure...
Definition: estimators.hpp:72
void SetLocalErrorGoal(double err_goal)
Set the local stopping criterion: stop when local_err_i &lt;= local_err_goal. The default value is zero...
const double alpha
Definition: ex15.cpp:336
class for C-function coefficient
Vector data type.
Definition: vector.hpp:48
virtual void UpdatesFinished()
Free the GridFunction update operator (if any), to save memory.
Definition: fespace.hpp:654
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
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...
double front(double x, double y, double z, double t, int)
Definition: ex15.cpp:339
double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
Definition: ex15.cpp:373
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:2382
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:287
De-refinement operator using an error threshold.
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:2285
bool Good() const
Definition: optparser.hpp:125