MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
mesh_operators.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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 int ApplyImpl(Mesh &mesh) override;
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 void Reset() override;
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 int ApplyImpl(Mesh &mesh) override;
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 void Reset() override;
258
259 /// Set the array @a refinements of elements to refine, without refining.
260 int MarkWithoutRefining(Mesh & mesh, Array<Refinement> & refinements);
261};
262
263// TODO: BulkRefiner to refine a portion of the global error
264
265
266/** @brief De-refinement operator using an error threshold.
267
268 This de-refinement operator marks elements in the hierarchy whose children
269 are leaves and their combined error is below a given threshold. The
270 errors of the children are combined by one of the following operations:
271 - op = 0: minimum of the errors
272 - op = 1: sum of the errors (default)
273 - op = 2: maximum of the errors. */
275{
276protected:
278
281
282 /** @brief Apply the operator to the mesh.
283 @return DEREFINED + CONTINUE if some elements were de-refined; NONE
284 otherwise. */
285 int ApplyImpl(Mesh &mesh) override;
286
287public:
288 /// Construct a ThresholdDerefiner using the given ErrorEstimator.
290 : estimator(est)
291 {
292 threshold = 0.0;
293 nc_limit = 0;
294 op = 1;
295 }
296
297 // default destructor (virtual)
298
299 /// Set the de-refinement threshold. The default value is zero.
300 void SetThreshold(real_t thresh) { threshold = thresh; }
301
302 void SetOp(int oper) { op = oper; }
303
304 /** @brief Set the maximum ratio of refinement levels of adjacent elements
305 (0 = unlimited). */
306 void SetNCLimit(int nc_limit_)
307 {
308 MFEM_ASSERT(nc_limit_ >= 0, "Invalid NC limit");
309 nc_limit = nc_limit_;
310 }
311
312 /// Reset the associated estimator.
313 void Reset() override { estimator.Reset(); }
314};
315
316
317/** @brief Refinement operator to control data oscillation.
318
319 This class computes osc_K(f) := || h ⋅ (I - Π) f ||_K at each element K.
320 Here, Π is the L2-projection and ||⋅||_K is the L2-norm, restricted to the
321 element K. All elements satisfying the inequality
322 \code
323 osc_K(f) > threshold ⋅ ||f|| / sqrt(n_el),
324 \endcode
325 are refined. Here, threshold is a positive parameter, ||⋅|| is the L2-norm
326 over the entire domain Ω, and n_el is the number of elements in the mesh.
327
328 Note that if osc(f) = threshold ⋅ ||f|| / sqrt(n_el) for each K, then
329 \code
330 osc(f) = sqrt(sum_K osc_K^2(f)) = threshold ⋅ ||f||.
331 \endcode
332 This is the reason for the 1/sqrt(n_el) factor.
333*/
335{
336protected:
337 bool print_level = false;
338 int nc_limit = 1;
340 int order;
341 long long max_elements = std::numeric_limits<long long>::max();
348 const IntegrationRule **irs = NULL;
349
350 /** @brief Apply the operator to the mesh once.
351 @return STOP if a stopping criterion is satisfied or no elements were
352 marked for refinement; REFINED + CONTINUE otherwise. */
353 int ApplyImpl(Mesh &mesh) override;
354
355public:
356 /// Constructor
357 CoefficientRefiner(Coefficient &coeff_, int order_)
358 {
359 // function f
360 coeff = &coeff_;
361
362 // order of the projection Π
363 order = order_;
364 }
365
366 /** @brief Apply the operator to the mesh max_it times or until tolerance
367 * achieved.
368 @return STOP if a stopping criterion is satisfied or no elements were
369 marked for refinement; REFINED + CONTINUE otherwise. */
370 virtual int PreprocessMesh(Mesh &mesh, int max_it);
371
373 {
374 int max_it = 10;
375 return PreprocessMesh(mesh, max_it);
376 }
377
378 /// Set the refinement threshold. The default value is 1.0e-2.
379 void SetThreshold(real_t threshold_) { threshold = threshold_; }
380
381 /** @brief Set the maximum number of elements stopping criterion: stop when
382 the input mesh has num_elements >= max_elem. The default value is
383 LONG_MAX. */
384 void SetMaxElements(long long max_elements_) { max_elements = max_elements_; }
385
386 /// Reset the function f
388 {
390 global_osc = NAN;
391 coeff = &coeff_;
392 }
393
394 /// Reset the oscillation order
395 void SetOrder(int order_) { order = order_; }
396
397 /** @brief Set the maximum ratio of refinement levels of adjacent elements
398 (0 = unlimited). The default value is 1, which helps ensure appropriate
399 refinements in pathological situations where the default quadrature
400 order is too low. */
401 void SetNCLimit(int nc_limit_)
402 {
403 MFEM_ASSERT(nc_limit_ >= 0, "Invalid NC limit");
404 nc_limit = nc_limit_;
405 }
406
407 // Set a custom integration rule
408 void SetIntRule(const IntegrationRule *irs_[]) { irs = irs_; }
409
410 // Set print level
411 void PrintWarnings() { print_level = true; }
412
413 // Return the value of the global relative data oscillation
414 real_t GetOsc() const { return global_osc; }
415
416 // Return the local relative data oscillation errors
417 const Vector & GetLocalOscs() const
418 {
419 MFEM_ASSERT(element_oscs.Size() > 0,
420 "Local oscillations have not been computed yet")
421 return element_oscs;
422 }
423
424 /// Reset
425 void Reset() override;
426};
427
428
429/** @brief ParMesh rebalancing operator.
430
431 If the mesh is a parallel mesh, perform rebalancing; otherwise, do nothing.
432*/
434{
435protected:
436 /** @brief Rebalance a parallel mesh (only non-conforming parallel meshes are
437 supported).
438 @return CONTINUE + REBALANCE on success, NONE otherwise. */
439 int ApplyImpl(Mesh &mesh) override;
440
441public:
442 /// Empty.
443 void Reset() override { }
444};
445
446} // namespace mfem
447
448#endif // MFEM_MESH_OPERATORS
The AnisotropicErrorEstimator class is the base class for all error estimators that compute one non-n...
Refinement operator to control data oscillation.
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_[])
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 Reset() override
Reset.
void SetThreshold(real_t threshold_)
Set the refinement threshold. The default value is 1.0e-2.
int ApplyImpl(Mesh &mesh) override
Apply the operator to the mesh once.
void ResetCoefficient(Coefficient &coeff_)
Reset the function f.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Base class for all element based error estimators.
virtual void Reset()=0
Force recomputation of the estimates on the next call to GetLocalErrors.
static const int NumGeom
Definition geom.hpp:46
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.
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.
int ApplyImpl(Mesh &mesh) override
Apply the MeshOperatorSequence.
virtual ~MeshOperatorSequence()
Delete all operators from the sequence.
void Reset() override
Reset all MeshOperators in 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:64
ParMesh rebalancing operator.
int ApplyImpl(Mesh &mesh) override
Rebalance a parallel mesh (only non-conforming parallel meshes are supported).
void Reset() override
Empty.
De-refinement operator using an error threshold.
void Reset() override
Reset the associated estimator.
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.
int ApplyImpl(Mesh &mesh) override
Apply the operator to the mesh.
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...
ThresholdRefiner(ErrorEstimator &est)
Construct a ThresholdRefiner using the given ErrorEstimator.
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.
void Reset() override
Reset the associated estimator.
int MarkWithoutRefining(Mesh &mesh, Array< Refinement > &refinements)
Set the array refinements of elements to refine, without refining.
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
int ApplyImpl(Mesh &mesh) override
Apply the operator to the mesh.
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:82
void Destroy()
Destroy a vector.
Definition vector.hpp:635
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition vector.hpp:47
float real_t
Definition config.hpp:43