MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
tmop_amr.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_TMOP_AMR_HPP
13 #define MFEM_TMOP_AMR_HPP
14 
15 #include "tmop_tools.hpp"
16 #include "nonlinearform.hpp"
17 #include "pnonlinearform.hpp"
18 #include "estimators.hpp"
19 #include "../mesh/mesh_operators.hpp"
20 
21 namespace mfem
22 {
23 
25 {
26 protected:
27  Mesh *mesh; // not owned
28  NonlinearForm *nlf; // not owned
29  int order;
30  int amrmetric;
35  // An element is refined only if
36  // [mean TMOPEnergy(children)]*energy_scaling_factor < TMOPEnergy(parent)
38  GridFunction *spat_gf; // If specified, can be used to specify the
39  double spat_gf_critical; // region where hr-adaptivity is done.
40 
41  /// Check if the mesh of the solution was modified.
43  {
44  long mesh_sequence = mesh->GetSequence();
45  MFEM_ASSERT(mesh_sequence >= current_sequence, "");
46  return (mesh_sequence > current_sequence);
47  }
48 
49  /// Compute the element error estimates. For an element E in the mesh,
50  /// error(E) = TMOPEnergy(E)*energy_scaling_factor-Mean(TMOPEnergy(ChildofE)),
51  /// where TMOPEnergy of Children of E is obtained by assuming the element E
52  /// is refined using the refinement type being considered based on the TMOP
53  /// mesh quality metric.
54  void ComputeEstimates();
55 
56  /// Construct the integration rules to model how each element type is split
57  /// using different refinement types. ref_type = 0 is the original element
58  /// and reftype \ in [1, 7] represent different refinement type based on
59  /// NCMesh class.
60  void SetQuadIntRules(); // supports ref_type = 1 to 3.
61  void SetTriIntRules(); // currently supports only isotropic refinement.
62  void SetHexIntRules(); // currently supports only isotropic refinement.
63  void SetTetIntRules(); // currently supports only isotropic refinement.
64 
65  /// Get TMOP energy for each element corresponding to the refinement type
66  /// specified.
67  void GetTMOPRefinementEnergy(int reftype, Vector &el_energy_vec);
68 
69  /// Use a mesh to setup an integration rule that will mimic the different
70  /// refinement types.
72 public:
73  TMOPRefinerEstimator(Mesh &mesh_, NonlinearForm &nlf_, int order_,
74  int amrmetric_) :
75  mesh(&mesh_), nlf(&nlf_), order(order_), amrmetric(amrmetric_),
79  {
80  if (mesh->Dimension() == 2)
81  {
84  }
85  else
86  {
89  }
90  }
91 
93  {
94  for (int i = 0; i < QuadIntRule.Size(); i++) { delete QuadIntRule[i]; }
95  for (int i = 0; i < TriIntRule.Size(); i++) { delete TriIntRule[i]; }
96  for (int i = 0; i < HexIntRule.Size(); i++) { delete HexIntRule[i]; }
97  for (int i = 0; i < TetIntRule.Size(); i++) { delete TetIntRule[i]; }
98  }
99 
100  /// Get TMOP-based errors for each element in the mesh computed based on the
101  /// refinement types being considered.
102  virtual const Vector &GetLocalErrors()
103  {
104  if (MeshIsModified()) { ComputeEstimates(); }
105  return error_estimates;
106  }
107  /// For anisotropic refinements, get the refinement type (e.g., x or y)
109  {
110  if (MeshIsModified()) { ComputeEstimates(); }
111  return aniso_flags;
112  }
113 
114  /// Scaling factor for the TMOP refinement energy. An element is refined if
115  /// [mean TMOPEnergy(children)]*energy_scaling_factor < TMOPEnergy(parent)
116  void SetEnergyScalingFactor(double scale) { energy_scaling_factor = scale; }
117 
118  /// Spatial indicator function (eta) that can be used to prevent elements
119  /// from being refined even if the energy criterion is met. Using this,
120  /// an element E is not refined if mean(@a spat_gf(E)) < @a spat_gf_critical.
122  double spat_gf_critical_ = 0.5)
123  { spat_gf = &spat_gf_; spat_gf_critical = spat_gf_critical_; }
124  void SetSpatialIndicatorCritical(double val_) { spat_gf_critical = val_; }
125 
126  /// Reset the error estimator.
127  virtual void Reset() { current_sequence = -1; }
128 };
129 
131 {
132 protected:
135 #ifdef MFEM_USE_MPI
138 #endif
139  int order;
143  bool serial;
144 
145  /// Check if the mesh of the solution was modified.
147  {
148  long mesh_sequence = mesh->GetSequence();
149  MFEM_ASSERT(mesh_sequence >= current_sequence, "");
150  return (mesh_sequence > current_sequence);
151  }
152 
153  /// Compute the element error estimates. For a given element E in the mesh,
154  /// error(E) = TMOPEnergy(parent_of_E)-TMOPEnergy(E). Children element of an
155  /// element are derefined if the mean TMOP energy of children is greater than
156  /// the TMOP energy associated with their parent.
157  void ComputeEstimates();
158 
159  void GetTMOPDerefinementEnergy(Mesh &cmesh,
160  TMOP_Integrator &tmopi,
161  Vector &el_energy_vec);
162 
164  Vector &fine_energy);
165 public:
167  mesh(&mesh_), nlf(&nlf_),
168  current_sequence(-1), error_estimates(), serial(true) { }
169 #ifdef MFEM_USE_MPI
171  mesh(&pmesh_), nlf(&pnlf_), pmesh(&pmesh_), pnlf(&pnlf_),
172  current_sequence(-1), error_estimates(), serial(false) { }
173 #endif
174 
176 
177  virtual const Vector &GetLocalErrors()
178  {
179  if (MeshIsModified()) { ComputeEstimates(); }
180  return error_estimates;
181  }
182 
183  /// Reset the error estimator.
184  virtual void Reset() { current_sequence = -1; }
185 };
186 
187 // hr-adaptivity using TMOP.
188 // If hr-adaptivity is disabled, r-adaptivity is done once using the
189 // TMOPNewtonSolver.
190 // Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
191 // "h_per_r_iter" iterations of h-adaptivity after each r-adaptivity iteration.
192 // The solver terminates early if an h-adaptivity iteration does not
193 // refine/derefine any element in the mesh.
195 {
196 protected:
205 #ifdef MFEM_USE_MPI
210 #endif
211  bool serial;
212 
213  // All are owned.
218 
220 
221  void Update();
222 #ifdef MFEM_USE_MPI
223  void ParUpdate();
224 #endif
226 
227 #ifdef MFEM_USE_MPI
228  // Rebalance ParMesh such that all the children elements are moved to the same
229  // MPI rank where the parent will be if the mesh were to be derefined.
230  void RebalanceParNCMesh();
231 #endif
232 
233 public:
234  TMOPHRSolver(Mesh &mesh_, NonlinearForm &nlf_,
235  TMOPNewtonSolver &tmopns_, GridFunction &x_,
236  bool move_bnd_, bool hradaptivity_,
237  int mesh_poly_deg_, int amr_metric_id_,
238  int hr_iter_ = 5, int h_per_r_iter_ = 1);
239 #ifdef MFEM_USE_MPI
240  TMOPHRSolver(ParMesh &pmesh_, ParNonlinearForm &pnlf_,
241  TMOPNewtonSolver &tmopns_, ParGridFunction &x_,
242  bool move_bnd_, bool hradaptivity_,
243  int mesh_poly_deg_, int amr_metric_id_,
244  int hr_iter_ = 5, int h_per_r_iter_ = 1);
245 #endif
246 
247  void Mult();
248 
249  /// These are used to update spaces and functions that are not owned by the
250  /// TMOPIntegrator or DiscreteAdaptTC. The owned ones are updated in the
251  /// functions UpdateAfterMeshTopologyChange() of both classes.
254 
255 #ifdef MFEM_USE_MPI
257  {
258  pgridfuncarr.Append(pgf_);
259  }
261  {
262  pfespacearr.Append(pfes_);
263  }
264 #endif
265 
267  {
268  if (!hradaptivity) { return; }
269  delete tmop_dr;
270  delete tmop_dr_est;
271  delete tmop_r;
272  delete tmop_r_est;
273  }
274 
275  /// Total number of hr-adaptivity iterations. At each iteration, we do an
276  /// r-adaptivity iteration followed by a number of h-adaptivity iterations.
277  void SetHRAdaptivityIterations(int iter) { hr_iter = iter; }
278 
279  /// Total number of h-adaptivity iterations per r-adaptivity iteration.
280  void SetHAdaptivityIterations(int iter) { h_per_r_iter = iter; }
281 };
282 
283 }
284 #endif
virtual const Vector & GetLocalErrors()
Get a Vector with all element errors.
Definition: tmop_amr.hpp:177
ThresholdRefiner * tmop_r
Definition: tmop_amr.hpp:215
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void GetTMOPRefinementEnergy(int reftype, Vector &el_energy_vec)
Definition: tmop_amr.cpp:87
virtual const Vector & GetLocalErrors()
Definition: tmop_amr.hpp:102
bool MeshIsModified()
Check if the mesh of the solution was modified.
Definition: tmop_amr.hpp:146
bool MeshIsModified()
Check if the mesh of the solution was modified.
Definition: tmop_amr.hpp:42
Array< int > aniso_flags
Definition: tmop_amr.hpp:34
Parallel non-linear operator on the true dofs.
void AddFESpaceForUpdate(ParFiniteElementSpace *pfes_)
Definition: tmop_amr.hpp:260
Abstract parallel finite element space.
Definition: pfespace.hpp:28
TMOPHRSolver(Mesh &mesh_, NonlinearForm &nlf_, TMOPNewtonSolver &tmopns_, GridFunction &x_, bool move_bnd_, bool hradaptivity_, int mesh_poly_deg_, int amr_metric_id_, int hr_iter_=5, int h_per_r_iter_=1)
Definition: tmop_amr.cpp:535
void SetSpatialIndicator(GridFunction &spat_gf_, double spat_gf_critical_=0.5)
Definition: tmop_amr.hpp:121
void SetHAdaptivityIterations(int iter)
Total number of h-adaptivity iterations per r-adaptivity iteration.
Definition: tmop_amr.hpp:280
void AddGridFunctionForUpdate(ParGridFunction *pgf_)
Definition: tmop_amr.hpp:256
Array< FiniteElementSpace * > fespacearr
Definition: tmop_amr.hpp:202
NonlinearForm * nlf
Definition: tmop_amr.hpp:28
long GetSequence() const
Definition: mesh.hpp:1664
void SetEnergyScalingFactor(double scale)
Definition: tmop_amr.hpp:116
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:751
ThresholdDerefiner * tmop_dr
Definition: tmop_amr.hpp:217
Mesh refinement operator using an error threshold.
Array< ParGridFunction * > pgridfuncarr
Definition: tmop_amr.hpp:208
Array< IntegrationRule * > QuadIntRule
Definition: tmop_amr.hpp:31
TMOPRefinerEstimator(Mesh &mesh_, NonlinearForm &nlf_, int order_, int amrmetric_)
Definition: tmop_amr.hpp:73
virtual void Reset()
Reset the error estimator.
Definition: tmop_amr.hpp:184
GridFunction * spat_gf
Definition: tmop_amr.hpp:38
virtual const Array< int > & GetAnisotropicFlags()
For anisotropic refinements, get the refinement type (e.g., x or y)
Definition: tmop_amr.hpp:108
Base class for all element based error estimators.
Definition: estimators.hpp:41
void AddFESpaceForUpdate(FiniteElementSpace *fes)
Definition: tmop_amr.hpp:253
Array< IntegrationRule * > TetIntRule
Definition: tmop_amr.hpp:31
void GetTMOPDerefinementEnergy(Mesh &cmesh, TMOP_Integrator &tmopi, Vector &el_energy_vec)
Definition: tmop_amr.cpp:509
int Dimension() const
Definition: mesh.hpp:1006
void AddGridFunctionForUpdate(GridFunction *gf)
Definition: tmop_amr.hpp:252
Array< GridFunction * > gridfuncarr
Definition: tmop_amr.hpp:201
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
void SetSpatialIndicatorCritical(double val_)
Definition: tmop_amr.hpp:124
void RebalanceParNCMesh()
Definition: tmop_amr.cpp:729
const int mesh_poly_deg
Definition: tmop_amr.hpp:204
TMOPRefinerEstimator * tmop_r_est
Definition: tmop_amr.hpp:214
Array< IntegrationRule * > HexIntRule
Definition: tmop_amr.hpp:31
ParNonlinearForm * pnlf
Definition: tmop_amr.hpp:137
ParNonlinearForm * pnlf
Definition: tmop_amr.hpp:207
const int amr_metric_id
Definition: tmop_amr.hpp:204
IntegrationRule * SetIntRulesFromMesh(Mesh &meshsplit)
Definition: tmop_amr.cpp:317
TMOPDeRefinerEstimator * tmop_dr_est
Definition: tmop_amr.hpp:216
The AnisotropicErrorEstimator class is the base class for all error estimators that compute one non-n...
Definition: estimators.hpp:64
Array< IntegrationRule * > TriIntRule
Definition: tmop_amr.hpp:31
Array< ParFiniteElementSpace * > pfespacearr
Definition: tmop_amr.hpp:209
GridFunction * x
Definition: tmop_amr.hpp:200
void UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf)
Definition: tmop_amr.cpp:840
Vector data type.
Definition: vector.hpp:60
TMOPDeRefinerEstimator(Mesh &mesh_, NonlinearForm &nlf_)
Definition: tmop_amr.hpp:166
NonlinearForm * nlf
Definition: tmop_amr.hpp:198
Class for parallel grid function.
Definition: pgridfunc.hpp:32
TMOPDeRefinerEstimator(ParMesh &pmesh_, ParNonlinearForm &pnlf_)
Definition: tmop_amr.hpp:170
Class for parallel meshes.
Definition: pmesh.hpp:32
virtual void Reset()
Reset the error estimator.
Definition: tmop_amr.hpp:127
TMOPNewtonSolver * tmopns
Definition: tmop_amr.hpp:199
void SetHRAdaptivityIterations(int iter)
Definition: tmop_amr.hpp:277
bool GetDerefineEnergyForIntegrator(TMOP_Integrator &tmopi, Vector &fine_energy)
Definition: tmop_amr.cpp:357
De-refinement operator using an error threshold.
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:1565