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