Flow123d  master-f44eb46
mesh_optimizer.hh
Go to the documentation of this file.
1 
2 #ifndef MESH_OPTIMIZER_HH_
3 #define MESH_OPTIMIZER_HH_
4 
5 #include <vector>
6 
7 #include <armadillo>
8 
9 #include "mesh/mesh.h"
10 #include "mesh/accessors.hh"
11 #include "mesh/range_wrapper.hh"
12 #include "bounding_box.hh"
13 
14 struct Permutee {
15  inline Permutee(uint original_index, double curve_value) : original_index_(original_index), curve_value_(curve_value) {}
17  double curve_value_;
18 };
19 
20 struct Normalizer {
21  inline Normalizer() : shift_({0, 0, 0}), scalar_(1) {}
22  inline Normalizer(arma::vec3 shift, double scalar) : shift_(shift), scalar_(scalar) {}
24  return (vec - shift_) / scalar_;
25  }
27  double scalar_;
28 };
29 
30 inline bool operator<(const Permutee& first, const Permutee& second) {
31  return first.curve_value_ < second.curve_value_;
32 }
33 
34 template <uint DIM>
36  static_assert(DIM == 2 || DIM == 3, "DIM must be either 2 or 3.");
37 public:
38  inline MeshOptimizer(Mesh* mesh) : mesh_(mesh) {}
39 
40  inline void calculate_sizes() {
41  BoundingBox bounding_box( *mesh_->node(0) );
42  for (auto nod : mesh_->node_range()) {
43  bounding_box.expand(*nod);
44  }
45  double mesh_longest_size = bounding_box.longest_size();
46  normalizer_ = Normalizer(bounding_box.min(), mesh_longest_size);
47 
48  node_sizes_.resize(mesh_->n_nodes(), INFINITY);
49  element_sizes_.reserve(mesh_->n_elements());
50  for (const ElementAccessor<3>& elm : mesh_->elements_range()) {
51  double elm_norm_size = elm.bounding_box().longest_size() / mesh_longest_size;
52  element_sizes_.push_back(elm_norm_size);
53  for (uint i = 0; i < elm.dim() + 1; ++i) {
54  node_sizes_[elm->node_idx(i)] = std::min({ node_sizes_[elm->node_idx(i)], elm_norm_size });
55  }
56  }
57  }
58 
60  node_refs_.reserve(mesh_->n_nodes());
61  uint i = 0;
62  for (auto nod : mesh_->node_range()) {
63  node_refs_.emplace_back(i, hilbert_value(normalizer_.normalize( *nod ), node_sizes_[i]));
64  i++;
65  }
66  }
67 
69  node_refs_.reserve(mesh_->n_nodes());
70  uint i = 0;
71  for (auto nod : mesh_->node_range()) {
72  node_refs_.emplace_back(i, zcurve_value(normalizer_.normalize( *nod ), node_sizes_[i]));
73  }
74  }
75 
77  node_refs_.reserve(mesh_->n_nodes());
78  for (uint i = 0; i < mesh_->n_elements(); ++i) {
79  auto elm = mesh_->element_accessor(element_refs_[i].original_index_);
80  for (uint j = 0; j < elm->dim() + 1; ++j) {
81  node_refs_.emplace_back(elm->node_idx(j), element_refs_[i].curve_value_);
82  }
83  }
84  }
85 
87  element_refs_.reserve(mesh_->n_elements());
88  uint i = 0;
89  for (auto elm : mesh_->elements_range()) {
90  element_refs_.emplace_back(i, hilbert_value(normalizer_.normalize(elm.centre()), element_sizes_[i]));
91  i++;
92  }
93  }
94 
96  element_refs_.reserve(mesh_->n_elements());
97  uint i = 0;
98  for (auto elm : mesh_->elements_range()) {
99  element_refs_.emplace_back(i, zcurve_value(normalizer_.normalize(elm.centre()), element_sizes_[i]));
100  i++;
101  }
102  }
103 
105  return this->sort(node_refs_, node_permutation);
106  }
107 
109  return this->sort(element_refs_, elem_permutation);
110  }
111 
112 private:
119 
120  inline double hilbert_value(double x, double y, double eps) const {
121  if (eps > 1) {
122  return 0;
123  } else {
124  if (x < 0.5) {
125  if (y < 0.5) {
126  return hilbert_value(2 * y, 2 * x, 4 * eps) / 4;
127  } else {
128  return (1 + hilbert_value(2 * x, 2 * y - 1, 4 * eps)) / 4;
129  }
130  } else {
131  if (y >= 0.5) {
132  return (2 + hilbert_value(2 * x - 1, 2 * y - 1, 4 * eps)) / 4;
133  } else {
134  return (3 + hilbert_value(1 - 2 * y, 2 - 2 * x, 4 * eps)) / 4;
135  }
136  }
137  }
138  }
139 
140  inline double hilbert_value(double x, double y, double z, double eps) const {
141  if (eps > 1) {
142  return 0;
143  } else {
144  if (z < 0.5) {
145  if (x < 0.5) {
146  if (y < 0.5) {
147  return hilbert_value(2 * z, 2 * x, 2 * y, 8 * eps) / 8;
148  } else {
149  return (1 + hilbert_value(2 * y - 1, 2 * z, 2 * x, 8 * eps)) / 8;
150  }
151  } else {
152  if (y >= 0.5) {
153  return (2 + hilbert_value(2 * y - 1, 2 * z, 2 * x - 1, 8 * eps)) / 8;
154  } else {
155  return (3 + hilbert_value(2 - 2 * x, 1 - 2 * y, 2 * z, 8 * eps)) / 8;
156  }
157  }
158  } else {
159  if (x >= 0.5) {
160  if (y < 0.5) {
161  return (4 + hilbert_value(2 - 2 * x, 1 - 2 * y, 2 * z - 1, 8 * eps)) / 8;
162  } else {
163  return (5 + hilbert_value(2 * y - 1, 2 - 2 * z, 2 - 2 * x, 8 * eps)) / 8;
164  }
165  } else {
166  if (y >= 0.5) {
167  return (6 + hilbert_value(2 * y - 1, 2 - 2 * z, 1 - 2 * x, 8 * eps)) / 8;
168  } else {
169  return (7 + hilbert_value(2 - 2 * z, 2 * x, 1 - 2 * y, 8 * eps)) / 8;
170  }
171  }
172  }
173  }
174  }
175 
176  inline double zcurve_value(double x, double y, double eps) const {
177  if (eps > 1) {
178  return 0;
179  } else {
180  if (y < 0.5) {
181  if (x < 0.5) {
182  return zcurve_value(2 * x, 2 * y, 4 * eps) / 4;
183  } else {
184  return (1 + zcurve_value(2 * x - 1, 2 * y, 4 * eps)) / 4;
185  }
186  } else {
187  if (x < 0.5) {
188  return (2 + zcurve_value(2 * x, 2 * y - 1, 4 * eps)) / 4;
189  } else {
190  return (3 + zcurve_value(2 * x - 1, 2 * y - 1, 4 * eps)) / 4;
191  }
192  }
193  }
194  }
195 
196  inline double zcurve_value(double x, double y, double z, double eps) const {
197  if (eps > 1) {
198  return 0;
199  } else {
200  if (z < 0.5) {
201  if (y < 0.5) {
202  if (x < 0.5) {
203  return zcurve_value(2 * x, 2 * y, 2 * z, 8 * eps) / 8;
204  } else {
205  return (1 + zcurve_value(2 * x - 1, 2 * y, 2 * z, 8 * eps)) / 8;
206  }
207  } else {
208  if (x < 0.5) {
209  return (2 + zcurve_value(2 * x, 2 * y - 1, 2 * z, 8 * eps)) / 8;
210  } else {
211  return (3 + zcurve_value(2 * x - 1, 2 * y - 1, 2 * z, 8 * eps)) / 8;
212  }
213  }
214  } else {
215  if (y < 0.5) {
216  if (x < 0.5) {
217  return (4 + zcurve_value(2 * x, 2 * y, 2 * z - 1, 8 * eps)) / 8;
218  } else {
219  return (5 + zcurve_value(2 * x - 1, 2 * y, 2 * z - 1, 8 * eps)) / 8;
220  }
221  } else {
222  if (x < 0.5) {
223  return (6 + zcurve_value(2 * x, 2 * y - 1, 2 * z - 1, 8 * eps)) / 8;
224  } else {
225  return (7 + zcurve_value(2 * x - 1, 2 * y - 1, 2 * z - 1, 8 * eps)) / 8;
226  }
227  }
228  }
229  }
230  }
231 
232  inline double hilbert_value(arma::vec3 vec, double size) const {
233  return hilbert_value(vec[0], vec[1], vec[2], size * size * size);
234  }
235 
236  inline double zcurve_value(arma::vec3 vec, double size) const {
237  return zcurve_value(vec[0], vec[1], vec[2], size * size * size);
238  }
239 
241  ASSERT(refs.size() <= mesh_perm.size());
242  std::sort(refs.begin(), refs.end());
243  std::vector<int> mesh_ids;
244  mesh_ids.reserve(refs.size());
245  for (uint i = 0; i < refs.size(); ++i) {
246  mesh_perm[ refs[i].original_index_ ] = i;
247  mesh_ids.push_back( refs[i].original_index_ );
248  }
249  return mesh_ids;
250  }
251 
252 };
253 
254 template<>
255 inline double MeshOptimizer<2>::hilbert_value(arma::vec3 vec, double size) const {
256  return hilbert_value(vec[0], vec[1], size * size);
257 }
258 
259 template <>
260 inline double MeshOptimizer<2>::zcurve_value(arma::vec3 vec, double size) const {
261  return zcurve_value(vec[0], vec[1], size * size);
262 }
263 
264 #endif /* MESH_OPTIMIZER_HH_ */
#define ASSERT(expr)
Definition: asserts.hh:351
Bounding box in 3d ambient space.
Definition: bounding_box.hh:53
const Point & min() const
double longest_size() const
void expand(const Point &point)
unsigned int n_nodes() const
Definition: mesh.h:171
NodeAccessor< 3 > node(unsigned int idx) const
Create and return NodeAccessor to node of given idx.
Definition: mesh.cc:872
unsigned int n_elements() const
Definition: mesh.h:111
Range< ElementAccessor< 3 > > elements_range() const
Returns range of mesh elements.
Definition: mesh.cc:1174
Range< NodeAccessor< 3 > > node_range() const
Returns range of nodes.
Definition: mesh.cc:1180
ElementAccessor< 3 > element_accessor(unsigned int idx) const
Create and return ElementAccessor to element of given idx.
Definition: mesh.cc:866
void calculate_sizes()
MeshOptimizer(Mesh *mesh)
void calculate_element_curve_values_as_zcurve_of_center()
std::vector< double > node_sizes_
std::vector< int > sort_elements(std::vector< unsigned int > &elem_permutation)
void calculate_node_curve_values_as_zcurve()
std::vector< double > element_sizes_
void calculate_element_curve_values_as_hilbert_of_centers()
std::vector< int > sort_nodes(std::vector< unsigned int > &node_permutation)
double hilbert_value(double x, double y, double z, double eps) const
void calculate_node_curve_values_as_obtained_from_elements()
double hilbert_value(arma::vec3 vec, double size) const
double zcurve_value(double x, double y, double z, double eps) const
void calculate_node_curve_values_as_hilbert()
std::vector< Permutee > element_refs_
double zcurve_value(double x, double y, double eps) const
std::vector< Permutee > node_refs_
std::vector< int > sort(std::vector< Permutee > &refs, std::vector< unsigned int > &mesh_perm)
double zcurve_value(arma::vec3 vec, double size) const
double hilbert_value(double x, double y, double eps) const
Normalizer normalizer_
Definition: mesh.h:362
bool operator<(const Permutee &first, const Permutee &second)
unsigned int uint
ArmaVec< double, N > vec
Definition: armor.hh:885
Implementation of range helper class.
double scalar_
arma::vec3 shift_
arma::vec3 normalize(const arma::vec3 vec)
Normalizer(arma::vec3 shift, double scalar)
uint original_index_
Permutee(uint original_index, double curve_value)
double curve_value_