MFEM  v4.0
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
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  /** @brief Apply the MeshOperatorSequence.
136  @return ActionInfo value corresponding to the last applied operator from
137  the sequence. */
138  virtual int ApplyImpl(Mesh &mesh);
139 
140 public:
141  /// Constructor. Use the Append() method to create the sequence.
143 
144  /// Delete all operators from the sequence.
145  virtual ~MeshOperatorSequence();
146 
147  /** @brief Add an operator to the end of the sequence.
148  The MeshOperatorSequence assumes ownership of the operator. */
149  void Append(MeshOperator *mc) { sequence.Append(mc); }
150 
151  /// Access the underlying sequence.
153 
154  /// Reset all MeshOperators in the sequence.
155  virtual void Reset();
156 };
157 
158 
159 /** @brief Mesh refinement operator using an error threshold.
160 
161  This class uses the given ErrorEstimator to estimate local element errors
162  and then marks for refinement all elements i such that loc_err_i > threshold.
163  The threshold is computed as
164  \code
165  threshold = max(total_err * total_fraction * pow(num_elements,-1.0/p),
166  local_err_goal);
167  \endcode
168  where p (=total_norm_p), total_fraction, and local_err_goal are settable
169  parameters, total_err = (sum_i local_err_i^p)^{1/p}, when p < inf,
170  or total_err = max_i local_err_i, when p = inf.
171 */
173 {
174 protected:
177 
178  double total_norm_p;
183 
184  double threshold;
186 
189 
191  int nc_limit;
192 
193  double GetNorm(const Vector &local_err, Mesh &mesh) const;
194 
195  /** @brief Apply the operator to the mesh.
196  @return STOP if a stopping criterion is satisfied or no elements were
197  marked for refinement; REFINED + CONTINUE otherwise. */
198  virtual int ApplyImpl(Mesh &mesh);
199 
200 public:
201  /// Construct a ThresholdRefiner using the given ErrorEstimator.
203 
204  // default destructor (virtual)
205 
206  /** @brief Set the exponent, p, of the discrete p-norm used to compute the
207  total error from the local element errors. */
208  void SetTotalErrorNormP(double norm_p = infinity())
209  { total_norm_p = norm_p; }
210 
211  /** @brief Set the total error stopping criterion: stop when
212  total_err <= total_err_goal. The default value is zero. */
213  void SetTotalErrorGoal(double err_goal) { total_err_goal = err_goal; }
214 
215  /** @brief Set the total fraction used in the computation of the threshold.
216  The default value is 1/2.
217  @note If fraction == 0, total_err is essentially ignored in the threshold
218  computation, i.e. threshold = local error goal. */
219  void SetTotalErrorFraction(double fraction) { total_fraction = fraction; }
220 
221  /** @brief Set the local stopping criterion: stop when
222  local_err_i <= local_err_goal. The default value is zero.
223  @note If local_err_goal == 0, it is essentially ignored in the threshold
224  computation. */
225  void SetLocalErrorGoal(double err_goal) { local_err_goal = err_goal; }
226 
227  /** @brief Set the maximum number of elements stopping criterion: stop when
228  the input mesh has num_elements >= max_elem. The default value is
229  LONG_MAX. */
230  void SetMaxElements(long max_elem) { max_elements = max_elem; }
231 
232  /// Use nonconforming refinement, if possible (triangles, quads, hexes).
234 
235  /** @brief Use conforming refinement, if possible (triangles, tetrahedra)
236  -- this is the default. */
238 
239  /** @brief Set the maximum ratio of refinement levels of adjacent elements
240  (0 = unlimited). */
242  {
243  MFEM_ASSERT(nc_limit >= 0, "Invalid NC limit");
244  this->nc_limit = nc_limit;
245  }
246 
247  /// Get the number of marked elements in the last Apply() call.
248  long GetNumMarkedElements() const { return num_marked_elements; }
249 
250  /// Get the threshold used in the last Apply() call.
251  double GetThreshold() const { return threshold; }
252 
253  /// Reset the associated estimator.
254  virtual void Reset();
255 };
256 
257 // TODO: BulkRefiner to refine a portion of the global error
258 
259 
260 /** @brief De-refinement operator using an error threshold.
261 
262  This de-refinement operator marks elements in the hierarchy whose children
263  are leaves and their combined error is below a given threshold. The
264  errors of the children are combined by one of the following operations:
265  - op = 0: minimum of the errors
266  - op = 1: sum of the errors (default)
267  - op = 2: maximum of the errors. */
269 {
270 protected:
272 
273  double threshold;
274  int nc_limit, op;
275 
276  /** @brief Apply the operator to the mesh.
277  @return DEREFINED + CONTINUE if some elements were de-refined; NONE
278  otherwise. */
279  virtual int ApplyImpl(Mesh &mesh);
280 
281 public:
282  /// Construct a ThresholdDerefiner using the given ErrorEstimator.
284  : estimator(est)
285  {
286  threshold = 0.0;
287  nc_limit = 0;
288  op = 1;
289  }
290 
291  // default destructor (virtual)
292 
293  /// Set the de-refinement threshold. The default value is zero.
294  void SetThreshold(double thresh) { threshold = thresh; }
295 
296  void SetOp(int op) { this->op = op; }
297 
298  /** @brief Set the maximum ratio of refinement levels of adjacent elements
299  (0 = unlimited). */
301  {
302  MFEM_ASSERT(nc_limit >= 0, "Invalid NC limit");
303  this->nc_limit = nc_limit;
304  }
305 
306  /// Reset the associated estimator.
307  virtual void Reset() { estimator.Reset(); }
308 };
309 
310 
311 /** @brief ParMesh rebalancing operator.
312 
313  If the mesh is a parallel mesh, perform rebalancing; otherwise, do nothing.
314 */
315 class Rebalancer : public MeshOperator
316 {
317 protected:
318  /** @brief Rebalance a parallel mesh (only non-conforming parallel meshes are
319  supported).
320  @return CONTINUE + REBALANCE on success, NONE otherwise. */
321  virtual int ApplyImpl(Mesh &mesh);
322 
323 public:
324  /// Empty.
325  virtual void Reset() { }
326 };
327 
328 } // namespace mfem
329 
330 #endif // MFEM_MESH_OPERATORS
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...
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh.
the mesh was refined
Action
Action and information constants and masks.
void SetNCLimit(int nc_limit)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited).
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
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().
virtual void Reset()
Reset all MeshOperators in the sequence.
virtual int ApplyImpl(Mesh &mesh)
Apply the MeshOperatorSequence.
bit mask for all &quot;action&quot; bits
long GetNumMarkedElements() const
Get the number of marked elements in the last Apply() call.
bool Apply(Mesh &mesh)
Perform the mesh operation.
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.
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:39
Array< Refinement > marked_elements
ErrorEstimator & estimator
virtual void Reset()=0
Force recomputation of the estimates on the next call to GetLocalErrors.
void SetNCLimit(int nc_limit)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited).
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 ~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 Append(MeshOperator *mc)
Add an operator to the end of the sequence. The MeshOperatorSequence assumes ownership of the operato...
void SetMaxElements(long max_elem)
Set the maximum number of elements stopping criterion: stop when the input mesh has num_elements &gt;= m...
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.
The AnisotropicErrorEstimator class is the base class for all error estimators that compute one non-n...
Definition: estimators.hpp:56
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:42
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:48
ThresholdRefiner(ErrorEstimator &est)
Construct a ThresholdRefiner using the given ErrorEstimator.
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.
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.