MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
mesh_operators.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 
22 namespace 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 {
48 private:
49  int mod;
50 
51 protected:
52  friend class MeshOperatorSequence;
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 
62 public:
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 {
128 protected:
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 
143 public:
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 {
177 protected:
180 
181  double total_norm_p;
185  long long max_elements;
186 
187  double threshold;
189 
192 
194  int nc_limit;
195 
196  double 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 
203 public:
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. */
211  void SetTotalErrorNormP(double norm_p = infinity())
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(double 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(double 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(double 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  double 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 {
273 protected:
275 
276  double threshold;
277  int nc_limit, op;
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 
284 public:
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(double 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 {
333 protected:
334  bool print_level = false;
335  int nc_limit = 1;
336  int nonconforming = -1;
337  int order;
338  long long max_elements = std::numeric_limits<long long>::max();
339  double threshold = 1.0e-2;
340  double global_osc = NAN;
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 
352 public:
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 
369  int PreprocessMesh(Mesh &mesh)
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(double 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  double 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 */
430 class Rebalancer : public MeshOperator
431 {
432 protected:
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 
438 public:
439  /// Empty.
440  virtual void Reset() { }
441 };
442 
443 } // namespace mfem
444 
445 #endif // MFEM_MESH_OPERATORS
const IntegrationRule * ir_default[Geometry::NumGeom]
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void SetTotalErrorNormP(double norm_p=infinity())
Set the exponent, p, of the discrete p-norm used to compute the total error from the local element er...
bool Continue() const
Check if CONTINUE action is requested, i.e. FiniteElementSpaces and GridFunctions need to be updated ...
MeshOperatorSequence(const MeshOperatorSequence &)
Do not allow copy construction, due to assumed ownership.
virtual void Reset()=0
Reset the MeshOperator.
bool Repeat() const
Check if REPEAT action is requested, i.e. FiniteElementSpaces and GridFunctions need to be updated...
static const int NumGeom
Definition: geom.hpp:42
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).
const Vector & GetLocalOscs() const
the mesh was refined
void ResetCoefficient(Coefficient &coeff_)
Reset the function f.
Action
Action and information constants and masks.
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
void SetIntRule(const IntegrationRule *irs_[])
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
virtual int PreprocessMesh(Mesh &mesh, int max_it)
Apply the operator to the mesh max_it times or until tolerance achieved.
void PreferNonconformingRefinement()
Use nonconforming refinement, if possible (triangles, quads, hexes).
double GetNorm(const Vector &local_err, Mesh &mesh) const
the mesh was de-refined
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh.
int GetActionInfo() const
Get the full ActionInfo value generated by the last call to Apply().
MeshOperatorSequence & operator=(const MeshOperatorSequence &s)=delete
Do not allow copy assignment, due to assumed ownership.
virtual void Reset()
Reset all MeshOperators in the sequence.
void SetMaxElements(long long max_elements_)
Set the maximum number of elements stopping criterion: stop when the input mesh has num_elements &gt;= m...
virtual int ApplyImpl(Mesh &mesh)
Apply the MeshOperatorSequence.
bit mask for all &quot;action&quot; bits
bool Apply(Mesh &mesh)
Perform the mesh operation.
int PreprocessMesh(Mesh &mesh)
long long GetNumMarkedElements() const
Get the number of marked elements in the last Apply() call.
void SetOrder(int order_)
Reset the oscillation order.
virtual void Reset()
Reset the associated estimator.
bit mask for the &quot;update&quot; bit
virtual void Reset()
Reset the associated estimator.
Mesh refinement operator using an error threshold.
const IntegrationRule ** irs
void SetNCLimit(int nc_limit_)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited). The default value is...
MeshOperator()
Constructor to be used by derived classes.
virtual int ApplyImpl(Mesh &mesh)
Rebalance a parallel mesh (only non-conforming parallel meshes are supported).
bit mask for all &quot;info&quot; bits
Base class for all element based error estimators.
Definition: estimators.hpp:41
Array< Refinement > marked_elements
ErrorEstimator & estimator
virtual void Reset()=0
Force recomputation of the estimates on the next call to GetLocalErrors.
bool Derefined() const
Check if the mesh was de-refined.
void SetTotalErrorGoal(double err_goal)
Set the total error stopping criterion: stop when total_err &lt;= total_err_goal. The default value is z...
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh once.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
Refinement operator to control data oscillation.
virtual ~MeshOperator()
The destructor is virtual.
The MeshOperator class serves as base for mesh manipulation classes.
void PreferConformingRefinement()
Use conforming refinement, if possible (triangles, tetrahedra) – this is the default.
virtual ~MeshOperatorSequence()
Delete all operators from the sequence.
a stopping criterion was satisfied
ThresholdDerefiner(ErrorEstimator &est)
Construct a ThresholdDerefiner using the given ErrorEstimator.
AnisotropicErrorEstimator * aniso_estimator
ErrorEstimator & estimator
the mesh was rebalanced
void SetThreshold(double threshold_)
Set the refinement threshold. The default value is 1.0e-2.
void Append(MeshOperator *mc)
Add an operator to the end of the sequence. The MeshOperatorSequence assumes ownership of the operato...
bool Refined() const
Check if the mesh was refined.
Array< MeshOperator * > & GetSequence()
Access the underlying sequence.
MeshOperatorSequence()
Constructor. Use the Append() method to create the sequence.
void Destroy()
Destroy a vector.
Definition: vector.hpp:590
The AnisotropicErrorEstimator class is the base class for all error estimators that compute one non-n...
Definition: estimators.hpp:64
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:46
double GetThreshold() const
Get the threshold used in the last Apply() call.
void SetThreshold(double thresh)
Set the de-refinement threshold. The default value is zero.
bool Rebalanced() const
Check if the mesh was rebalanced.
void SetLocalErrorGoal(double err_goal)
Set the local stopping criterion: stop when local_err_i &lt;= local_err_goal. The default value is zero...
Vector data type.
Definition: vector.hpp:60
ThresholdRefiner(ErrorEstimator &est)
Construct a ThresholdRefiner using the given ErrorEstimator.
virtual void Reset()
Reset.
RefCoord s[3]
void SetMaxElements(long long max_elem)
Set the maximum number of elements stopping criterion: stop when the input mesh has num_elements &gt;= m...
void SetTotalErrorFraction(double fraction)
Set the total fraction used in the computation of the threshold. The default value is 1/2...
virtual void Reset()
Empty.
Array< MeshOperator * > sequence
MeshOperators sequence, owned by us.
CoefficientRefiner(Coefficient &coeff_, int order_)
Constructor.
void SetNCLimit(int nc_limit_)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited).
virtual int ApplyImpl(Mesh &mesh)=0
Implementation of the mesh operation. Invoked by the Apply() public method.
ParMesh rebalancing operator.
De-refinement operator using an error threshold.