MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
mesh_operators.cpp
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 #include "mesh_operators.hpp"
13 #include "pmesh.hpp"
14 
15 namespace mfem
16 {
17 
19 {
20  // delete in reverse order
21  for (int i = sequence.Size()-1; i >= 0; i--)
22  {
23  delete sequence[i];
24  }
25 }
26 
28 {
29  if (sequence.Size() == 0) { return NONE; }
30 next_step:
31  step = (step + 1) % sequence.Size();
32  bool last = (step == sequence.Size() - 1);
33  int mod = sequence[step]->ApplyImpl(mesh);
34  switch (mod & MASK_ACTION)
35  {
36  case NONE: if (last) { return NONE; } goto next_step;
37  case CONTINUE: return last ? mod : (REPEAT | (mod & MASK_INFO));
38  case STOP: return STOP;
39  case REPEAT: --step; return mod;
40  }
41  return NONE;
42 }
43 
45 {
46  for (int i = 0; i < sequence.Size(); i++)
47  {
48  sequence[i]->Reset();
49  }
50  step = 0;
51 }
52 
53 
55  : estimator(est)
56 {
59  total_err_goal = 0.0;
60  total_fraction = 0.5;
61  local_err_goal = 0.0;
62  max_elements = std::numeric_limits<long long>::max();
63 
64  threshold = 0.0;
65  num_marked_elements = 0LL;
66  current_sequence = -1;
67 
68  non_conforming = -1;
69  nc_limit = 0;
70 }
71 
72 double ThresholdRefiner::GetNorm(const Vector &local_err, Mesh &mesh) const
73 {
74 #ifdef MFEM_USE_MPI
75  ParMesh *pmesh = dynamic_cast<ParMesh*>(&mesh);
76  if (pmesh)
77  {
78  return ParNormlp(local_err, total_norm_p, pmesh->GetComm());
79  }
80 #endif
81  return local_err.Normlp(total_norm_p);
82 }
83 
85 {
86  threshold = 0.0;
87  num_marked_elements = 0LL;
88  marked_elements.SetSize(0);
90 
91  const long long num_elements = mesh.GetGlobalNE();
92  if (num_elements >= max_elements) { return STOP; }
93 
94  const int NE = mesh.GetNE();
95  const Vector &local_err = estimator.GetLocalErrors();
96  MFEM_ASSERT(local_err.Size() == NE, "invalid size of local_err");
97 
98  const double total_err = GetNorm(local_err, mesh);
99  if (total_err <= total_err_goal) { return STOP; }
100 
101  if (total_norm_p < infinity())
102  {
103  threshold = std::max(total_err * total_fraction *
104  std::pow(num_elements, -1.0/total_norm_p),
106  }
107  else
108  {
109  threshold = std::max(total_err * total_fraction, local_err_goal);
110  }
111 
112  for (int el = 0; el < NE; el++)
113  {
114  if (local_err(el) > threshold)
115  {
116  marked_elements.Append(Refinement(el));
117  }
118  }
119 
120  if (aniso_estimator)
121  {
122  const Array<int> &aniso_flags = aniso_estimator->GetAnisotropicFlags();
123  if (aniso_flags.Size() > 0)
124  {
125  for (int i = 0; i < marked_elements.Size(); i++)
126  {
127  Refinement &ref = marked_elements[i];
128  ref.ref_type = aniso_flags[ref.index];
129  }
130  }
131  }
132 
134  if (num_marked_elements == 0LL) { return STOP; }
135 
137  return CONTINUE + REFINED;
138 }
139 
141 {
142  estimator.Reset();
143  current_sequence = -1;
144  num_marked_elements = 0LL;
145  // marked_elements.SetSize(0); // not necessary
146 }
147 
148 
150 {
151  if (mesh.Conforming()) { return NONE; }
152 
153  const Vector &local_err = estimator.GetLocalErrors();
154  bool derefs = mesh.DerefineByError(local_err, threshold, nc_limit, op);
155 
156  return derefs ? CONTINUE + DEREFINED : NONE;
157 }
158 
159 
161 {
162  int max_it = 1;
163  return PreprocessMesh(mesh, max_it);
164 }
165 
167 {
168  int rank = 0;
169  MFEM_VERIFY(max_it > 0, "max_it must be strictly positive")
170 
171  int dim = mesh.Dimension();
172  L2_FECollection l2fec(order, dim);
173  FiniteElementSpace* l2fes = NULL;
174 
175  bool par = false;
176  GridFunction *gf = NULL;
177 
178 #ifdef MFEM_USE_MPI
179  ParMesh* pmesh = dynamic_cast<ParMesh*>(&mesh);
180  if (pmesh && pmesh->Nonconforming())
181  {
182  par = true;
183  l2fes = new ParFiniteElementSpace(pmesh, &l2fec);
184  gf = new ParGridFunction(static_cast<ParFiniteElementSpace*>(l2fes));
185  }
186 #endif
187  if (!par)
188  {
189  l2fes = new FiniteElementSpace(&mesh, &l2fec);
190  gf = new GridFunction(l2fes);
191  }
192 
193  // If custom integration rule has not been set,
194  // then use the default integration rule
195  if (!irs)
196  {
197  int order_quad = 2*order + 3;
198  for (int i=0; i < Geometry::NumGeom; ++i)
199  {
200  ir_default[i] = &(IntRules.Get(i, order_quad));
201  }
202  irs = ir_default;
203  }
204 
205  for (int i = 0; i < max_it; i++)
206  {
207  // Compute number of elements and L2-norm of f.
208  int NE = mesh.GetNE();
209  int globalNE = 0;
210  double norm_of_coeff = 0.0;
211  if (par)
212  {
213 #ifdef MFEM_USE_MPI
214  globalNE = pmesh->GetGlobalNE();
215  norm_of_coeff = ComputeGlobalLpNorm(2.0,*coeff,*pmesh,irs);
216 #endif
217  }
218  else
219  {
220  globalNE = NE;
221  norm_of_coeff = ComputeLpNorm(2.0,*coeff,mesh,irs);
222  }
223 
224  // Compute average L2-norm of f
225  double av_norm_of_coeff = norm_of_coeff / sqrt(globalNE);
226 
227  // Compute element-wise L2-norms of (I - Π) f
228  Vector element_norms_of_fine_scale(NE);
229  gf->SetSpace(l2fes);
231  gf->ComputeElementL2Errors(*coeff,element_norms_of_fine_scale,irs);
232 
233  // Define osc_K(f) := || h ⋅ (I - Π) f ||_K and select elements
234  // for refinement based on threshold. Also record relative osc(f).
235  global_osc = 0.0;
238  element_oscs.SetSize(NE);
239  element_oscs = 0.0;
240  for (int j = 0; j < NE; j++)
241  {
242  double h = mesh.GetElementSize(j);
243  double element_osc = h * element_norms_of_fine_scale(j);
244  if ( element_osc > threshold * av_norm_of_coeff )
245  {
247  }
248  element_oscs(j) = element_osc/(norm_of_coeff + 1e-10);
249  global_osc += element_osc*element_osc;
250  }
251 #ifdef MFEM_USE_MPI
252  if (par)
253  {
254  MPI_Comm comm = pmesh->GetComm();
255  MPI_Allreduce(MPI_IN_PLACE, &global_osc, 1, MPI_DOUBLE, MPI_SUM, comm);
256  MPI_Comm_rank(comm, &rank);
257  }
258 #endif
259  global_osc = sqrt(global_osc)/(norm_of_coeff + 1e-10);
260 
261  // Exit if the global threshold or maximum number of elements is reached.
262  if (global_osc < threshold || globalNE > max_elements)
263  {
264  if (global_osc > threshold && globalNE > max_elements && rank == 0 &&
265  print_level)
266  {
267  MFEM_WARNING("Reached maximum number of elements "
268  "before resolving data to tolerance.");
269  }
270  delete l2fes;
271  delete gf;
272  return STOP;
273  }
274 
275  // Refine elements.
277  l2fes->Update(false);
278  gf->Update();
279 
280  }
281  delete l2fes;
282  delete gf;
283  return CONTINUE + REFINED;
284 
285 }
286 
288 {
290  global_osc = 0.0;
291  coeff = NULL;
292  irs = NULL;
293 }
294 
295 
297 {
298 #ifdef MFEM_USE_MPI
299  ParMesh *pmesh = dynamic_cast<ParMesh*>(&mesh);
300  if (pmesh && pmesh->Nonconforming())
301  {
302  pmesh->Rebalance();
303  return CONTINUE + REBALANCED;
304  }
305 #endif
306  return NONE;
307 }
308 
309 
310 } // namespace mfem
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
const IntegrationRule * ir_default[Geometry::NumGeom]
char ref_type
refinement XYZ bit mask (7 = full isotropic)
Definition: ncmesh.hpp:39
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3334
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
static const int NumGeom
Definition: geom.hpp:42
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh.
bool Conforming() const
Definition: mesh.hpp:1650
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
the mesh was refined
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
virtual int PreprocessMesh(Mesh &mesh, int max_it)
Apply the operator to the mesh max_it times or until tolerance achieved.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
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.
virtual void ComputeElementL2Errors(Coefficient &exsol, Vector &error, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:605
virtual void Reset()
Reset all MeshOperators in the sequence.
int index
Mesh element number.
Definition: ncmesh.hpp:38
virtual int ApplyImpl(Mesh &mesh)
Apply the MeshOperatorSequence.
bit mask for all &quot;action&quot; bits
long GetSequence() const
Definition: mesh.hpp:1664
void Rebalance()
Definition: pmesh.cpp:3989
bool Nonconforming() const
Definition: mesh.hpp:1651
double Normlp(double p) const
Returns the l_p norm of the vector.
Definition: vector.cpp:852
double ComputeLpNorm(double p, Coefficient &coeff, Mesh &mesh, const IntegrationRule *irs[])
Compute the Lp norm of a function f. .
double GetElementSize(ElementTransformation *T, int type=0)
Definition: mesh.cpp:76
virtual void Reset()
Reset the associated estimator.
virtual const Vector & GetLocalErrors()=0
Get a Vector with all element errors.
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:751
const IntegrationRule ** irs
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
int Dimension() const
Definition: mesh.hpp:1006
virtual const Array< int > & GetAnisotropicFlags()=0
Get an Array&lt;int&gt; with anisotropic flags for all mesh elements.
ErrorEstimator & estimator
virtual void Reset()=0
Force recomputation of the estimates on the next call to GetLocalErrors.
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh once.
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
virtual ~MeshOperatorSequence()
Delete all operators from the sequence.
double ComputeGlobalLpNorm(double p, Coefficient &coeff, ParMesh &pmesh, const IntegrationRule *irs[])
Compute the global Lp norm of a function f. .
a stopping criterion was satisfied
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:165
AnisotropicErrorEstimator * aniso_estimator
ErrorEstimator & estimator
the mesh was rebalanced
bool DerefineByError(Array< double > &elem_error, double threshold, int nc_limit=0, int op=1)
Definition: mesh.cpp:9212
long long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:955
int dim
Definition: ex24.cpp:53
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 ParNormlp(const Vector &vec, double p, MPI_Comm comm)
Compute the l_p norm of the Vector which is split without overlap across the given communicator...
Definition: hypre.cpp:424
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2399
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition: gridfunc.cpp:196
Vector data type.
Definition: vector.hpp:60
ThresholdRefiner(ErrorEstimator &est)
Construct a ThresholdRefiner using the given ErrorEstimator.
virtual long long ReduceInt(int value) const
Utility function: sum integers from all processors (Allreduce).
Definition: mesh.hpp:952
virtual void Reset()
Reset.
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Array< MeshOperator * > sequence
MeshOperators sequence, owned by us.
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:9405
Class for parallel meshes.
Definition: pmesh.hpp:32
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:288