MFEM  v3.3
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 {
58  total_norm_p = std::numeric_limits<double>::infinity();
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  double total_err = GetNorm(local_err, mesh);
99  if (total_err <= total_err_goal) { return STOP; }
100 
101  threshold = std::max(total_err * total_fraction *
102  std::pow(num_elements, -1.0/total_norm_p),
104 
105  for (int el = 0; el < NE; el++)
106  {
107  if (local_err(el) > threshold)
108  {
109  marked_elements.Append(Refinement(el));
110  }
111  }
112 
113  if (aniso_estimator)
114  {
115  const Array<int> &aniso_flags = aniso_estimator->GetAnisotropicFlags();
116  if (aniso_flags.Size() > 0)
117  {
118  for (int i = 0; i < marked_elements.Size(); i++)
119  {
120  Refinement &ref = marked_elements[i];
121  ref.ref_type = aniso_flags[ref.index];
122  }
123  }
124  }
125 
127  if (num_marked_elements == 0) { return STOP; }
128 
130  return CONTINUE + REFINED;
131 }
132 
134 {
135  estimator.Reset();
136  current_sequence = -1;
138  // marked_elements.SetSize(0); // not necessary
139 }
140 
141 
143 {
144  if (mesh.Conforming()) { return NONE; }
145 
146  const Vector &local_err = estimator.GetLocalErrors();
147  bool derefs = mesh.DerefineByError(local_err, threshold, nc_limit, op);
148 
149  return derefs ? CONTINUE + DEREFINED : NONE;
150 }
151 
152 
154 {
155 #ifdef MFEM_USE_MPI
156  ParMesh *pmesh = dynamic_cast<ParMesh*>(&mesh);
157  if (pmesh && pmesh->Nonconforming())
158  {
159  pmesh->Rebalance();
160  return CONTINUE + REBALANCED;
161  }
162 #endif
163  return NONE;
164 }
165 
166 
167 } // namespace mfem
int Size() const
Logical size of the array.
Definition: array.hpp:109
char ref_type
refinement XYZ bit mask (7 = full isotropic)
Definition: ncmesh.hpp:36
long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:602
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh.
bool Conforming() const
Definition: mesh.hpp:967
the mesh was refined
int Size() const
Returns the size of the vector.
Definition: vector.hpp:106
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:584
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:599
virtual void Reset()
Reset all MeshOperators in the sequence.
int index
Mesh element number.
Definition: ncmesh.hpp:35
virtual int ApplyImpl(Mesh &mesh)
Apply the MeshOperatorSequence.
bit mask for all &quot;action&quot; bits
long GetSequence() const
Definition: mesh.hpp:981
void Rebalance()
Load balance the mesh. NC meshes only.
Definition: pmesh.cpp:2713
bool Nonconforming() const
Definition: mesh.hpp:968
double Normlp(double p) const
Returns the l_p norm of the vector.
Definition: vector.cpp:692
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:116
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:6089
The AnisotropicErrorEstimator class is the base class for all error estimators that compute one non-n...
Definition: estimators.hpp:56
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:210
Vector data type.
Definition: vector.hpp:36
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:6285
Class for parallel meshes.
Definition: pmesh.hpp:28