MFEM  v4.0
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, 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 #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>::max();
63 
64  threshold = 0.0;
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;
88  marked_elements.SetSize(0);
90 
91  const 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 == 0) { return STOP; }
135 
137  return CONTINUE + REFINED;
138 }
139 
141 {
142  estimator.Reset();
143  current_sequence = -1;
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 #ifdef MFEM_USE_MPI
163  ParMesh *pmesh = dynamic_cast<ParMesh*>(&mesh);
164  if (pmesh && pmesh->Nonconforming())
165  {
166  pmesh->Rebalance();
167  return CONTINUE + REBALANCED;
168  }
169 #endif
170  return NONE;
171 }
172 
173 
174 } // namespace mfem
int Size() const
Logical size of the array.
Definition: array.hpp:118
char ref_type
refinement XYZ bit mask (7 = full isotropic)
Definition: ncmesh.hpp:37
long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:694
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh.
bool Conforming() const
Definition: mesh.hpp:1135
the mesh was refined
int Size() const
Returns the size of the vector.
Definition: vector.hpp:150
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:676
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 long ReduceInt(int value) const
Utility function: sum integers from all processors (Allreduce).
Definition: mesh.hpp:691
virtual void Reset()
Reset all MeshOperators in the sequence.
int index
Mesh element number.
Definition: ncmesh.hpp:36
virtual int ApplyImpl(Mesh &mesh)
Apply the MeshOperatorSequence.
bit mask for all &quot;action&quot; bits
long GetSequence() const
Definition: mesh.hpp:1149
void Rebalance()
Load balance the mesh. NC meshes only.
Definition: pmesh.cpp:3085
bool Nonconforming() const
Definition: mesh.hpp:1136
double Normlp(double p) const
Returns the l_p norm of the vector.
Definition: vector.cpp:766
virtual void Reset()
Reset the associated estimator.
virtual const Vector & GetLocalErrors()=0
Get a Vector with all element errors.
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
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.
MPI_Comm GetComm() const
Definition: pmesh.hpp:223
virtual ~MeshOperatorSequence()
Delete all operators from the sequence.
a stopping criterion was satisfied
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:7000
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 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:264
Vector data type.
Definition: vector.hpp:48
ThresholdRefiner(ErrorEstimator &est)
Construct a ThresholdRefiner using the given ErrorEstimator.
Array< MeshOperator * > sequence
MeshOperators sequence, owned by us.
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:7185
Class for parallel meshes.
Definition: pmesh.hpp:32