MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
mesh_operators.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#ifndef MFEM_MESH_OPERATORS
13#define MFEM_MESH_OPERATORS
14
15#include "../config/config.hpp"
16#include "../general/array.hpp"
17#include "mesh.hpp"
18#include "../fem/estimators.hpp"
19
20#include <limits>
21
22namespace mfem
23{
24
25/** @brief The MeshOperator class serves as base for mesh manipulation classes.
26
27 The purpose of the class is to provide a common abstraction for various
28 AMR mesh control schemes. The typical use in an AMR loop is illustrated
29 in examples 6/6p and 15/15p.
30
31 A more general loop that also supports sequences of mesh operators with
32 multiple updates looks like this:
33 \code
34 for (...)
35 {
36 // computations on the current mesh ...
37 while (mesh_operator->Apply(mesh))
38 {
39 // update FiniteElementSpaces and interpolate GridFunctions ...
40 if (mesh_operator->Continue()) { break; }
41 }
42 if (mesh_operator->Stop()) { break; }
43 }
44 \endcode
45 */
47{
48private:
49 int mod;
50
51protected:
53
54 /** @brief Implementation of the mesh operation. Invoked by the Apply()
55 public method.
56 @return Combination of ActionInfo constants. */
57 virtual int ApplyImpl(Mesh &mesh) = 0;
58
59 /// Constructor to be used by derived classes.
60 MeshOperator() : mod(NONE) { }
61
62public:
63 /** @brief Action and information constants and masks.
64
65 Combinations of constants are returned by the Apply() virtual method and
66 can be accessed directly with GetActionInfo() or indirectly with methods
67 like Stop(), Continue(), etc. The information bits (MASK_INFO) can be set
68 only when the update bit is set (see MASK_UPDATE). */
69 enum Action
70 {
71 NONE = 0, /**< continue with computations without updating spaces
72 or grid-functions, i.e. the mesh was not modified */
73 CONTINUE = 1, /**< update spaces and grid-functions and continue
74 computations with the new mesh */
75 STOP = 2, ///< a stopping criterion was satisfied
76 REPEAT = 3, /**< update spaces and grid-functions and call the
77 operator Apply() method again */
78 MASK_UPDATE = 1, ///< bit mask for the "update" bit
79 MASK_ACTION = 3 ///< bit mask for all "action" bits
80 };
81
82 enum Info
83 {
84 REFINED = 4*1, ///< the mesh was refined
85 DEREFINED = 4*2, ///< the mesh was de-refined
86 REBALANCED = 4*3, ///< the mesh was rebalanced
87 MASK_INFO = ~3 ///< bit mask for all "info" bits
88 };
89
90 /** @brief Perform the mesh operation.
91 @return true if FiniteElementSpaces and GridFunctions need to be updated.
92 */
93 bool Apply(Mesh &mesh) { return ((mod = ApplyImpl(mesh)) & MASK_UPDATE); }
94
95 /** @brief Check if STOP action is requested, e.g. stopping criterion is
96 satisfied. */
97 bool Stop() const { return ((mod & MASK_ACTION) == STOP); }
98 /** @brief Check if REPEAT action is requested, i.e. FiniteElementSpaces and
99 GridFunctions need to be updated, and Apply() must be called again. */
100 bool Repeat() const { return ((mod & MASK_ACTION) == REPEAT); }
101 /** @brief Check if CONTINUE action is requested, i.e. FiniteElementSpaces
102 and GridFunctions need to be updated and computations should continue. */
103 bool Continue() const { return ((mod & MASK_ACTION) == CONTINUE); }
104
105 /// Check if the mesh was refined.
106 bool Refined() const { return ((mod & MASK_INFO) == REFINED); }
107 /// Check if the mesh was de-refined.
108 bool Derefined() const { return ((mod & MASK_INFO) == DEREFINED); }
109 /// Check if the mesh was rebalanced.
110 bool Rebalanced() const { return ((mod & MASK_INFO) == REBALANCED); }
111
112 /** @brief Get the full ActionInfo value generated by the last call to
113 Apply(). */
114 int GetActionInfo() const { return mod; }
115
116 /// Reset the MeshOperator.
117 virtual void Reset() = 0;
118
119 /// The destructor is virtual.
120 virtual ~MeshOperator() { }
121};
122
123
124/** Composition of MeshOperators into a sequence. Use the Append() method to
125 create the sequence. */
127{
128protected:
129 int step;
130 Array<MeshOperator*> sequence; ///< MeshOperators sequence, owned by us.
131
132 /// Do not allow copy construction, due to assumed ownership.
134
135 /// Do not allow copy assignment, due to assumed ownership.
137
138 /** @brief Apply the MeshOperatorSequence.
139 @return ActionInfo value corresponding to the last applied operator from
140 the sequence. */
141 virtual int ApplyImpl(Mesh &mesh);
142
143public:
144 /// Constructor. Use the Append() method to create the sequence.
146
147 /// Delete all operators from the sequence.
148 virtual ~MeshOperatorSequence();
149
150 /** @brief Add an operator to the end of the sequence.
151 The MeshOperatorSequence assumes ownership of the operator. */
152 void Append(MeshOperator *mc) { sequence.Append(mc); }
153
154 /// Access the underlying sequence.
156
157 /// Reset all MeshOperators in the sequence.
158 virtual void Reset();
159};
160
161
162/** @brief Mesh refinement operator using an error threshold.
163
164 This class uses the given ErrorEstimator to estimate local element errors
165 and then marks for refinement all elements i such that loc_err_i > threshold.
166 The threshold is computed as
167 \code
168 threshold = max(total_err * total_fraction * pow(num_elements,-1.0/p),
169 local_err_goal);
170 \endcode
171 where p (=total_norm_p), total_fraction, and local_err_goal are settable
172 parameters, total_err = (sum_i local_err_i^p)^{1/p}, when p < inf,
173 or total_err = max_i local_err_i, when p = inf.
174*/
176{
177protected:
180
185 long long max_elements;
186
189
192
195
196 real_t GetNorm(const Vector &local_err, Mesh &mesh) const;
197
198 /** @brief Apply the operator to the mesh.
199 @return STOP if a stopping criterion is satisfied or no elements were
200 marked for refinement; REFINED + CONTINUE otherwise. */
201 virtual int ApplyImpl(Mesh &mesh);
202
203public:
204 /// Construct a ThresholdRefiner using the given ErrorEstimator.
206
207 // default destructor (virtual)
208
209 /** @brief Set the exponent, p, of the discrete p-norm used to compute the
210 total error from the local element errors. */
212 { total_norm_p = norm_p; }
213
214 /** @brief Set the total error stopping criterion: stop when
215 total_err <= total_err_goal. The default value is zero. */
216 void SetTotalErrorGoal(real_t err_goal) { total_err_goal = err_goal; }
217
218 /** @brief Set the total fraction used in the computation of the threshold.
219 The default value is 1/2.
220 @note If fraction == 0, total_err is essentially ignored in the threshold
221 computation, i.e. threshold = local error goal. */
222 void SetTotalErrorFraction(real_t fraction) { total_fraction = fraction; }
223
224 /** @brief Set the local stopping criterion: stop when
225 local_err_i <= local_err_goal. The default value is zero.
226 @note If local_err_goal == 0, it is essentially ignored in the threshold
227 computation. */
228 void SetLocalErrorGoal(real_t err_goal) { local_err_goal = err_goal; }
229
230 /** @brief Set the maximum number of elements stopping criterion: stop when
231 the input mesh has num_elements >= max_elem. The default value is
232 LONG_MAX. */
233 void SetMaxElements(long long max_elem) { max_elements = max_elem; }
234
235 /// Use nonconforming refinement, if possible (triangles, quads, hexes).
237
238 /** @brief Use conforming refinement, if possible (triangles, tetrahedra)
239 -- this is the default. */
241
242 /** @brief Set the maximum ratio of refinement levels of adjacent elements
243 (0 = unlimited). */
244 void SetNCLimit(int nc_limit_)
245 {
246 MFEM_ASSERT(nc_limit_ >= 0, "Invalid NC limit");
247 nc_limit = nc_limit_;
248 }
249
250 /// Get the number of marked elements in the last Apply() call.
251 long long GetNumMarkedElements() const { return num_marked_elements; }
252
253 /// Get the threshold used in the last Apply() call.
254 real_t GetThreshold() const { return threshold; }
255
256 /// Reset the associated estimator.
257 virtual void Reset();
258};
259
260// TODO: BulkRefiner to refine a portion of the global error
261
262
263/** @brief De-refinement operator using an error threshold.
264
265 This de-refinement operator marks elements in the hierarchy whose children
266 are leaves and their combined error is below a given threshold. The
267 errors of the children are combined by one of the following operations:
268 - op = 0: minimum of the errors
269 - op = 1: sum of the errors (default)
270 - op = 2: maximum of the errors. */
272{
273protected:
275
278
279 /** @brief Apply the operator to the mesh.
280 @return DEREFINED + CONTINUE if some elements were de-refined; NONE
281 otherwise. */
282 virtual int ApplyImpl(Mesh &mesh);
283
284public:
285 /// Construct a ThresholdDerefiner using the given ErrorEstimator.
287 : estimator(est)
288 {
289 threshold = 0.0;
290 nc_limit = 0;
291 op = 1;
292 }
293
294 // default destructor (virtual)
295
296 /// Set the de-refinement threshold. The default value is zero.
297 void SetThreshold(real_t thresh) { threshold = thresh; }
298
299 void SetOp(int oper) { op = oper; }
300
301 /** @brief Set the maximum ratio of refinement levels of adjacent elements
302 (0 = unlimited). */
303 void SetNCLimit(int nc_limit_)
304 {
305 MFEM_ASSERT(nc_limit_ >= 0, "Invalid NC limit");
306 nc_limit = nc_limit_;
307 }
308
309 /// Reset the associated estimator.
310 virtual void Reset() { estimator.Reset(); }
311};
312
313
314/** @brief Refinement operator to control data oscillation.
315
316 This class computes osc_K(f) := || h ⋅ (I - Π) f ||_K at each element K.
317 Here, Π is the L2-projection and ||⋅||_K is the L2-norm, restricted to the
318 element K. All elements satisfying the inequality
319 \code
320 osc_K(f) > threshold ⋅ ||f|| / sqrt(n_el),
321 \endcode
322 are refined. Here, threshold is a positive parameter, ||⋅|| is the L2-norm
323 over the entire domain Ω, and n_el is the number of elements in the mesh.
324
325 Note that if osc(f) = threshold ⋅ ||f|| / sqrt(n_el) for each K, then
326 \code
327 osc(f) = sqrt(sum_K osc_K^2(f)) = threshold ⋅ ||f||.
328 \endcode
329 This is the reason for the 1/sqrt(n_el) factor.
330*/
332{
333protected:
334 bool print_level = false;
335 int nc_limit = 1;
337 int order;
338 long long max_elements = std::numeric_limits<long long>::max();
345 const IntegrationRule **irs = NULL;
346
347 /** @brief Apply the operator to the mesh once.
348 @return STOP if a stopping criterion is satisfied or no elements were
349 marked for refinement; REFINED + CONTINUE otherwise. */
350 virtual int ApplyImpl(Mesh &mesh);
351
352public:
353 /// Constructor
354 CoefficientRefiner(Coefficient &coeff_, int order_)
355 {
356 // function f
357 coeff = &coeff_;
358
359 // order of the projection Π
360 order = order_;
361 }
362
363 /** @brief Apply the operator to the mesh max_it times or until tolerance
364 * achieved.
365 @return STOP if a stopping criterion is satisfied or no elements were
366 marked for refinement; REFINED + CONTINUE otherwise. */
367 virtual int PreprocessMesh(Mesh &mesh, int max_it);
368
370 {
371 int max_it = 10;
372 return PreprocessMesh(mesh, max_it);
373 }
374
375 /// Set the refinement threshold. The default value is 1.0e-2.
376 void SetThreshold(real_t threshold_) { threshold = threshold_; }
377
378 /** @brief Set the maximum number of elements stopping criterion: stop when
379 the input mesh has num_elements >= max_elem. The default value is
380 LONG_MAX. */
381 void SetMaxElements(long long max_elements_) { max_elements = max_elements_; }
382
383 /// Reset the function f
385 {
387 global_osc = NAN;
388 coeff = &coeff_;
389 }
390
391 /// Reset the oscillation order
392 void SetOrder(int order_) { order = order_; }
393
394 /** @brief Set the maximum ratio of refinement levels of adjacent elements
395 (0 = unlimited). The default value is 1, which helps ensure appropriate
396 refinements in pathological situations where the default quadrature
397 order is too low. */
398 void SetNCLimit(int nc_limit_)
399 {
400 MFEM_ASSERT(nc_limit_ >= 0, "Invalid NC limit");
401 nc_limit = nc_limit_;
402 }
403
404 // Set a custom integration rule
405 void SetIntRule(const IntegrationRule *irs_[]) { irs = irs_; }
406
407 // Set print level
408 void PrintWarnings() { print_level = true; }
409
410 // Return the value of the global relative data oscillation
411 real_t GetOsc() const { return global_osc; }
412
413 // Return the local relative data oscillation errors
414 const Vector & GetLocalOscs() const
415 {
416 MFEM_ASSERT(element_oscs.Size() > 0,
417 "Local oscillations have not been computed yet")
418 return element_oscs;
419 }
420
421 /// Reset
422 virtual void Reset();
423};
424
425
426/** @brief ParMesh rebalancing operator.
427
428 If the mesh is a parallel mesh, perform rebalancing; otherwise, do nothing.
429*/
431{
432protected:
433 /** @brief Rebalance a parallel mesh (only non-conforming parallel meshes are
434 supported).
435 @return CONTINUE + REBALANCE on success, NONE otherwise. */
436 virtual int ApplyImpl(Mesh &mesh);
437
438public:
439 /// Empty.
440 virtual void Reset() { }
441};
442
443} // namespace mfem
444
445#endif // MFEM_MESH_OPERATORS
The AnisotropicErrorEstimator class is the base class for all error estimators that compute one non-n...
Definition: estimators.hpp:65
Refinement operator to control data oscillation.
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh once.
virtual void Reset()
Reset.
const IntegrationRule * ir_default[Geometry::NumGeom]
void SetOrder(int order_)
Reset the oscillation order.
void SetNCLimit(int nc_limit_)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited). The default value is...
void SetMaxElements(long long max_elements_)
Set the maximum number of elements stopping criterion: stop when the input mesh has num_elements >= m...
const Vector & GetLocalOscs() const
void SetIntRule(const IntegrationRule *irs_[])
int PreprocessMesh(Mesh &mesh)
virtual int PreprocessMesh(Mesh &mesh, int max_it)
Apply the operator to the mesh max_it times or until tolerance achieved.
CoefficientRefiner(Coefficient &coeff_, int order_)
Constructor.
const IntegrationRule ** irs
void SetThreshold(real_t threshold_)
Set the refinement threshold. The default value is 1.0e-2.
void ResetCoefficient(Coefficient &coeff_)
Reset the function f.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:42
Base class for all element based error estimators.
Definition: estimators.hpp:42
virtual void Reset()=0
Force recomputation of the estimates on the next call to GetLocalErrors.
static const int NumGeom
Definition: geom.hpp:42
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:100
MeshOperatorSequence(const MeshOperatorSequence &)
Do not allow copy construction, due to assumed ownership.
Array< MeshOperator * > sequence
MeshOperators sequence, owned by us.
Array< MeshOperator * > & GetSequence()
Access the underlying sequence.
virtual void Reset()
Reset all MeshOperators in the sequence.
virtual int ApplyImpl(Mesh &mesh)
Apply the MeshOperatorSequence.
void Append(MeshOperator *mc)
Add an operator to the end of the sequence. The MeshOperatorSequence assumes ownership of the operato...
MeshOperatorSequence & operator=(const MeshOperatorSequence &s)=delete
Do not allow copy assignment, due to assumed ownership.
virtual ~MeshOperatorSequence()
Delete all operators from the sequence.
MeshOperatorSequence()
Constructor. Use the Append() method to create the sequence.
The MeshOperator class serves as base for mesh manipulation classes.
virtual void Reset()=0
Reset the MeshOperator.
bool Derefined() const
Check if the mesh was de-refined.
bool Repeat() const
Check if REPEAT action is requested, i.e. FiniteElementSpaces and GridFunctions need to be updated,...
bool Apply(Mesh &mesh)
Perform the mesh operation.
int GetActionInfo() const
Get the full ActionInfo value generated by the last call to Apply().
virtual int ApplyImpl(Mesh &mesh)=0
Implementation of the mesh operation. Invoked by the Apply() public method.
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
Action
Action and information constants and masks.
@ MASK_ACTION
bit mask for all "action" bits
@ STOP
a stopping criterion was satisfied
@ MASK_UPDATE
bit mask for the "update" bit
bool Rebalanced() const
Check if the mesh was rebalanced.
MeshOperator()
Constructor to be used by derived classes.
@ MASK_INFO
bit mask for all "info" bits
@ REFINED
the mesh was refined
@ REBALANCED
the mesh was rebalanced
@ DEREFINED
the mesh was de-refined
bool Continue() const
Check if CONTINUE action is requested, i.e. FiniteElementSpaces and GridFunctions need to be updated ...
virtual ~MeshOperator()
The destructor is virtual.
bool Refined() const
Check if the mesh was refined.
Mesh data type.
Definition: mesh.hpp:56
ParMesh rebalancing operator.
virtual void Reset()
Empty.
virtual int ApplyImpl(Mesh &mesh)
Rebalance a parallel mesh (only non-conforming parallel meshes are supported).
De-refinement operator using an error threshold.
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh.
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).
ThresholdDerefiner(ErrorEstimator &est)
Construct a ThresholdDerefiner using the given ErrorEstimator.
ErrorEstimator & estimator
virtual void Reset()
Reset the associated estimator.
Mesh refinement operator using an error threshold.
Array< Refinement > marked_elements
void PreferNonconformingRefinement()
Use nonconforming refinement, if possible (triangles, quads, hexes).
void SetTotalErrorNormP(real_t norm_p=infinity())
Set the exponent, p, of the discrete p-norm used to compute the total error from the local element er...
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh.
void SetNCLimit(int nc_limit_)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited).
AnisotropicErrorEstimator * aniso_estimator
void SetMaxElements(long long max_elem)
Set the maximum number of elements stopping criterion: stop when the input mesh has num_elements >= m...
long long GetNumMarkedElements() const
Get the number of marked elements in the last Apply() call.
real_t GetThreshold() const
Get the threshold used in the last Apply() call.
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.
real_t GetNorm(const Vector &local_err, Mesh &mesh) const
void PreferConformingRefinement()
Use conforming refinement, if possible (triangles, tetrahedra) – this is the default.
void SetTotalErrorGoal(real_t err_goal)
Set the total error stopping criterion: stop when total_err <= total_err_goal. The default value is z...
ErrorEstimator & estimator
void SetLocalErrorGoal(real_t err_goal)
Set the local stopping criterion: stop when local_err_i <= local_err_goal. The default value is zero.
Vector data type.
Definition: vector.hpp:80
void Destroy()
Destroy a vector.
Definition: vector.hpp:615
int Size() const
Returns the size of the vector.
Definition: vector.hpp:218
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:45
float real_t
Definition: config.hpp:43
RefCoord s[3]