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