Flow123d  master-f44eb46
accessors.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file accessors.hh
15  * @brief
16  */
17 
18 #ifndef ACCESSORS_HH_
19 #define ACCESSORS_HH_
20 
21 #include "mesh/bounding_box.hh"
22 #include "mesh/region.hh"
23 #include "mesh/elements.h"
24 #include "mesh/mesh.h"
25 #include "mesh/bc_mesh.hh"
26 #include "mesh/node_accessor.hh"
27 #include "mesh/ref_element.hh"
28 #include "la/distribution.hh"
29 #include "mesh/point.hh"
30 #include <vector>
31 #include <armadillo>
32 
33 
34 /**
35  * Due to circular dependence of return parameters in mesh accessors,
36  * and the intention to have all methods inlined,
37  * we gathered the accessors into a single file.
38  * This way, it is possible to implement it, see the simple test below.
39  *
40  * Do not include accessors_impl.hh file anywhere but at the end of this file!
41  *
42  * previous include loops:
43  * side -> boundary -> elemen accessor
44  * element accessor <-> side
45  * edge <-> side
46  * boundary -> edge
47  */
48 
49 // Compilable test of class loop dependence of return parameters.
50 // class A;
51 // class B;
52 
53 // class A{
54 // public:
55 // A()
56 // {};
57 
58 // B create_b();
59 // };
60 
61 // class B{
62 // public:
63 // B()
64 // {};
65 
66 // A create_a();
67 // };
68 
69 
70 // B A::create_b()
71 // { return B(); }
72 
73 // A B::create_a()
74 // { return A(); }
75 
76 
77 class Side;
78 class SiderIter;
79 class Edge;
80 class Boundary;
81 
82 /**
83  * Element accessor templated just by dimension of the embedding space, used by Fields.
84  * This should allow algorithms over elements where dimension of particular element is runtime parameter.
85  *
86  * This class suites as interface of Fields to the mesh elements, in particular this accessor knows directly
87  * the region, and also can be used as an accessor that works on the whole region if used by Fields that do not depend on
88  * particular elements as FieldConstant, FiledFormula, and FieldPython.
89  *
90  * TODO:
91  * - make this kind of accessor subclass of FieldCommonBase or at least move it into src/fields
92  * since it has functionality particular for Fields
93  *
94  * Ideas:
95  * need function to calculate intersection (object) of two ElementAccessors, but this definitely should be templated by
96  * dimension of the ref. element (or rather shape of ref. element), here we can have case dispatch
97  *
98  */
99 template <int spacedim>
101 public:
102  /// Default invalid accessor.
103  ElementAccessor();
104 
105  /// Regional accessor.
106  ElementAccessor(const MeshBase *mesh, RegionIdx r_idx);
107 
108  /// Element accessor.
109  ElementAccessor(const MeshBase *mesh, unsigned int idx);
110 
111  /// Incremental function of the Element iterator.
112  void inc();
113 
114  /// Computes the measure of the element.
115  double measure() const;
116 
117  inline bool inverted() const {
118  return element()->inverted;
119  }
120 
121  inline double sign() const {
122  return inverted() ? -1 : 1;
123  }
124 
125  /**
126  * Returns Jacobian of 3D element.
127  * Used by measure and in intersections.
128  */
129  inline double jacobian_S3() const {
130  ASSERT_EQ(dim(), 3).error("Dimension mismatch.");
131  return arma::dot( arma::cross(*( node(1) ) - *( node(0) ),
132  *( node(2) ) - *( node(0) )),
133  *( node(3) ) - *( node(0) )
134  );
135  }
136 
137  /**
138  * Returns Jacobian of 2D element.
139  */
140  inline double jacobian_S2() const {
141  ASSERT_EQ(dim(), 2).error("Dimension mismatch.");
142  return arma::norm(
143  arma::cross(*( node(1) ) - *( node(0) ), *( node(2) ) - *( node(0) )),
144  2
145  );
146  }
147 
148  /**
149  * Returns Jacobian of 1D element.
150  */
151  inline double jacobian_S1() const {
152  ASSERT_EQ(dim(), 1).error("Dimension mismatch.");
153  return arma::norm(*( node(1) ) - *( node(0) ) , 2);
154  }
155 
156  /// Computes the barycenter.
157  arma::vec::fixed<spacedim> centre() const;
158 
159  /**
160  * Quality of the element based on the smooth and scale-invariant quality measures proposed in:
161  * J. R. Schewchuk: What is a Good Linear Element?
162  *
163  * We scale the measure so that is gives value 1 for regular elements. Line 1d elements
164  * have always quality 1.
165  *
166  * We return signed value in order to detect inverted elements.
167  */
168  double quality_measure_smooth() const;
169 
170  SideIter side(const unsigned int loc_index);
171 
172  const SideIter side(const unsigned int loc_index) const;
173 
174 
175 
176  inline bool is_regional() const {
177  return r_idx_.is_valid();
178  }
179 
180  inline bool is_elemental() const {
181  return ( is_valid() && ! is_regional() );
182  }
183 
184  inline bool is_valid() const {
185  return mesh_ != NULL;
186  }
187 
188  inline unsigned int dim() const {
189  ASSERT(! is_regional());
190  return element()->dim();
191  }
192 
193  inline const Element * element() const {
194  return &(mesh_->element(element_idx_));
195  }
196 
197 
198  inline Region region() const
199  { return Region( region_idx(), mesh_->region_db()); }
200 
201  inline RegionIdx region_idx() const {
202  if (r_idx_.is_valid()) {
203  return r_idx_;
204  } else {
205  return element()->region_idx();
206  }
207  }
208 
209  /// We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular mesh
210  //unsigned int region_id() const {
211  // return region().id();
212  //}
213 
214  /// Return local idx of element in boundary / bulk part of element vector
215  unsigned int idx() const {
216  return element_idx_;
217  }
218 
219  /// Return the element ID in the input mesh. Should be only used for special output.
220  unsigned int input_id() const {
221  return (unsigned int)(mesh_->find_elem_id(element_idx_) );
222  }
223 
224  unsigned int proc() const {
225  ASSERT(is_elemental());
227  }
228 
229 
230  NodeAccessor<3> node(unsigned int ni) const {
231  return mesh_->node( element()->node_idx(ni) );
232  }
233 
234  /**
235  * Return bounding box of the element.
236  */
237  BoundingBox bounding_box() const;
238 
239 
240  inline auto &orig_nodes_order() const {
242  }
243 
244  bool operator==(const ElementAccessor<spacedim>& other) const {
245  return (mesh_ == other.mesh_) && (element_idx_ == other.element_idx_);
246  }
247 
248  inline bool operator!=(const ElementAccessor<spacedim>& other) const {
249  return (element_idx_ != other.element_idx_) || (mesh_ != other.mesh_);
250  }
251 
252  /**
253  * -> dereference operator
254  *
255  * Allow simplify calling of element() method. Example:
256  @code
257  ElementAccessor<3> elm_ac(mesh, index);
258  arma::vec centre;
259  centre = elm_ac.element()->node_idx(0); // full format of access to element
260  centre = elm_ac->node_idx(0); // short format with dereference operator
261  @endcode
262  */
263  inline const Element * operator ->() const {
264  return element();
265  }
266 
267 
268 private:
269  /**
270  * When dim_ == undefined_dim_ ; the value of element_idx_ is invalid.
271  * Is used for ElementAccessors for whole region
272  */
273  static const unsigned int undefined_dim_ = 100;
274 
275  /// Pointer to the mesh owning the element.
276  const MeshBase *mesh_;
277 
278  /// Index into Mesh::element_vec_ array.
279  unsigned int element_idx_;
280 
281  // Hack for sorption tables. TODO: remove
282  // undefined for regular elements
284 
285 };
286 
287 
288 
289 
290 
291 //=============================================================================
292 // Edge class
293 //=============================================================================
294 class Edge
295 {
296 public:
297  /// Default invalid edge accessor constructor.
298  Edge();
299 
300  /// Valid edge accessor constructor.
301  Edge(const MeshBase *mesh, unsigned int edge_idx);
302 
303  /// Gets side iterator of the @p i -th side.
304  SideIter side(const unsigned int i) const;
305 
306 
307 
308  bool is_valid() const
309  { return mesh_ != nullptr; }
310 
311  /// Returns edge global index.
312  unsigned int idx() const {
313  ASSERT(is_valid());
314  return edge_idx_;
315  }
316 
317  /// Incremental function of the Edge iterator.
318  void inc() {
319  ASSERT(is_valid()).error("Do not call inc() for invalid accessor!");
320  edge_idx_++;
321  }
322 
323  /// Comparison operator of the iterator.
324  bool operator==(const Edge& other) const{
325  return (mesh_ == other.mesh_ && edge_idx_ == other.edge_idx_);
326  }
327 
328  /// Returns number of sides aligned with the edge.
329  unsigned int n_sides() const
330  { return edge_data()->n_sides;}
331 
332 private:
333  /// Pointer to the mesh owning the node.
334  const MeshBase *mesh_;
335  /// Index into Mesh::edges vector.
336  unsigned int edge_idx_;
337 
338  /// Getter for edge data from mesh.
339  const EdgeData* edge_data() const;
340 };
341 
342 
343 
344 
345 
346 //=============================================================================
347 // Boundary class
348 //=============================================================================
349 class Boundary
350 {
351 public:
352  Boundary();
353  Boundary(BoundaryData* boundary_data);
354 
355  Edge edge();
357  Region region();
358  const Element * element();
359 
360  bool is_valid() const {
361  return boundary_data_ != nullptr;
362  }
363 
364  Mesh* mesh() {
365  ASSERT(is_valid());
366  return boundary_data_->mesh_;
367  }
368 
370  ASSERT(is_valid());
371  return boundary_data_->edge_idx_;
372  }
373 
375  ASSERT(is_valid());
376  return boundary_data_->bc_ele_idx_;
377  }
378 
379 private:
381 };
382 
383 
384 
385 
386 
387 //=============================================================================
388 // Side class
389 //=============================================================================
390 class Side {
391 public:
392  /// Default invalid side accessor constructor.
393  Side();
394 
395  /// Valid edge accessor constructor.
396  Side(const MeshBase * mesh, unsigned int elem_idx, unsigned int set_lnum);
397 
398  double measure() const; ///< Calculate metrics of the side
399  arma::vec3 centre() const; ///< Centre of side
400  arma::vec3 normal() const; ///< Vector of (generalized) normal
401  double diameter() const; ///< Calculate the side diameter.
402 
403  /// Returns dimension of the side, that is dimension of the element minus one.
404  unsigned int dim() const;
405 
406  /// Returns true for all sides either on boundary or connected to vb neigboring.
407  bool is_external() const;
408 
409  /// Returns true for side on the boundary.
410  bool is_boundary() const;
411 
412  /// Returns node for given local index @p i on the side.
413  NodeAccessor<3> node(unsigned int i) const;
414 
415  /// Returns iterator to the element of the side.
416  ElementAccessor<3> element() const;
417 
418  /// Returns global index of the edge connected to the side.
419  unsigned int edge_idx() const;
420 
421  /// Returns pointer to the edge connected to the side.
422  Edge edge() const;
423 
424  /**
425  * Returns boundary condition prescribed on the side.
426  * Fails on assert if side if not on boundary and no BC is prescribed.
427  */
428  Boundary cond() const;
429 
430  /// Returns global index of the prescribed boundary condition.
431  unsigned int cond_idx() const;
432 
433 
434 
435  /// Returns number of nodes of the side.
436  unsigned int n_nodes() const
437  { return dim()+1; }
438 
439  /// Returns pointer to the mesh.
440  const MeshBase * mesh() const
441  { return this->mesh_; }
442 
443  /// Returns local index of the side on the element.
444  unsigned int side_idx() const
445  { return side_idx_; }
446 
447  /// Returns index of element in Mesh::element_vec_.
448  unsigned int elem_idx() const
449  { return elem_idx_; }
450 
451  /// Returns true if the side has assigned element.
452  bool is_valid() const
453  { return mesh_!= nullptr; }
454 
455  /// Iterate over local sides of the element.
456  void inc() {
457  ASSERT(is_valid()).error("Do not call inc() for invalid accessor!");
458  side_idx_++;
459  }
460 
461  bool operator==(const Side &other) const {
462  return (mesh_ == other.mesh_ ) && ( elem_idx_ == other.elem_idx_ )
463  && ( side_idx_ == other.side_idx_ );
464  }
465 
466  bool operator!=(const Side &other) const {
467  return !( *this == other);
468  }
469 
470  /// This is necessary by current DofHandler, should change this
471  //void *make_ptr() const;
472 private:
473 
474  arma::vec3 normal_point() const;
475  arma::vec3 normal_line() const;
476  arma::vec3 normal_triangle() const;
477 
478  // Topology of the mesh
479 
480  const MeshBase * mesh_; ///< Pointer to Mesh to which belonged
481  unsigned int elem_idx_; ///< Index of element in Mesh::element_vec_
482  unsigned int side_idx_; ///< Local # of side in element (to remove it, we heve to remove calc_side_rhs)
483 
484 };
485 
486 
487 /*
488  * Iterator to a side.
489  */
490 class SideIter {
491 public:
493  {}
494 
495  SideIter(const Side &side)
496  : side_(side)
497  {}
498 
499  bool operator==(const SideIter &other) {
500  return (side_.mesh() == other.side_.mesh() ) && ( side_.elem_idx() == other.side_.elem_idx() )
501  && ( side_.side_idx() == other.side_.side_idx() );
502  }
503 
504 
505  bool operator!=(const SideIter &other) {
506  return !( *this == other);
507  }
508 
509  /// * dereference operator
510  const Side & operator *() const
511  { return side_; }
512 
513  /// -> dereference operator
514  const Side * operator ->() const
515  { return &side_; }
516 
517  /// prefix increment iterate only on local element
519  side_.inc();
520  return (*this);
521  }
522 
523 private:
525 };
526 
527 
528 
529 #include "mesh/accessors_impl.hh"
530 
531 #endif /* ACCESSORS_HH_ */
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
unsigned int bc_ele_idx_
Definition: mesh_data.hh:53
Mesh * mesh_
Definition: mesh_data.hh:54
unsigned int edge_idx_
Definition: mesh_data.hh:49
bool is_valid() const
Definition: accessors.hh:360
const Element * element()
uint bc_ele_idx()
Definition: accessors.hh:374
Region region()
BoundaryData * boundary_data_
Definition: accessors.hh:380
ElementAccessor< 3 > element_accessor()
Mesh * mesh()
Definition: accessors.hh:364
uint edge_idx()
Definition: accessors.hh:369
Bounding box in 3d ambient space.
Definition: bounding_box.hh:53
unsigned int get_proc(unsigned int idx) const
get processor of the given index
unsigned int n_sides
Definition: mesh_data.hh:29
const MeshBase * mesh_
Pointer to the mesh owning the node.
Definition: accessors.hh:334
Edge()
Default invalid edge accessor constructor.
unsigned int n_sides() const
Returns number of sides aligned with the edge.
Definition: accessors.hh:329
bool is_valid() const
Definition: accessors.hh:308
void inc()
Incremental function of the Edge iterator.
Definition: accessors.hh:318
unsigned int edge_idx_
Index into Mesh::edges vector.
Definition: accessors.hh:336
SideIter side(const unsigned int i) const
Gets side iterator of the i -th side.
bool operator==(const Edge &other) const
Comparison operator of the iterator.
Definition: accessors.hh:324
unsigned int idx() const
Returns edge global index.
Definition: accessors.hh:312
const EdgeData * edge_data() const
Getter for edge data from mesh.
bool inverted() const
Definition: accessors.hh:117
auto & orig_nodes_order() const
Definition: accessors.hh:240
bool is_elemental() const
Definition: accessors.hh:180
bool operator==(const ElementAccessor< spacedim > &other) const
Definition: accessors.hh:244
unsigned int input_id() const
Return the element ID in the input mesh. Should be only used for special output.
Definition: accessors.hh:220
bool is_valid() const
Definition: accessors.hh:184
static const unsigned int undefined_dim_
Definition: accessors.hh:273
const Element * operator->() const
Definition: accessors.hh:263
NodeAccessor< 3 > node(unsigned int ni) const
Definition: accessors.hh:230
const MeshBase * mesh_
Pointer to the mesh owning the element.
Definition: accessors.hh:276
bool operator!=(const ElementAccessor< spacedim > &other) const
Definition: accessors.hh:248
void inc()
Incremental function of the Element iterator.
unsigned int element_idx_
Index into Mesh::element_vec_ array.
Definition: accessors.hh:279
double sign() const
Definition: accessors.hh:121
double measure() const
Computes the measure of the element.
SideIter side(const unsigned int loc_index)
BoundingBox bounding_box() const
Region region() const
Definition: accessors.hh:198
double quality_measure_smooth() const
double jacobian_S1() const
Definition: accessors.hh:151
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
Definition: accessors.hh:215
unsigned int dim() const
Definition: accessors.hh:188
arma::vec::fixed< spacedim > centre() const
Computes the barycenter.
double jacobian_S3() const
Definition: accessors.hh:129
RegionIdx region_idx() const
Definition: accessors.hh:201
unsigned int proc() const
Definition: accessors.hh:224
double jacobian_S2() const
Definition: accessors.hh:140
RegionIdx r_idx_
Definition: accessors.hh:283
const Element * element() const
Definition: accessors.hh:193
ElementAccessor()
Default invalid accessor.
bool is_regional() const
Definition: accessors.hh:176
RegionIdx region_idx() const
Definition: elements.h:53
bool inverted
Inverted permutation of element nodes, negative Jacobian.
Definition: elements.h:88
uint permutation_
Index of permutation of input nodes.
Definition: elements.h:92
unsigned int dim() const
Definition: elements.h:120
Base class for Mesh and BCMesh.
Definition: mesh.h:96
std::array< std::array< uint, 4 >, 64 > element_nodes_original_
Definition: mesh.h:325
const RegionDB & region_db() const
Definition: mesh.h:175
NodeAccessor< 3 > node(unsigned int idx) const
Create and return NodeAccessor to node of given idx.
Definition: mesh.cc:872
int find_elem_id(unsigned int pos) const
Return element id (in GMSH file) of element of given position in element vector.
Definition: mesh.h:138
LongIdx * get_row_4_el() const
Definition: mesh.h:125
const Element & element(unsigned idx) const
Definition: mesh.h:128
Distribution * get_el_ds() const
Definition: mesh.h:119
Definition: mesh.h:362
bool is_valid() const
Returns false if the region has undefined/invalid value.
Definition: region.hh:77
Side side_
Definition: accessors.hh:524
const Side * operator->() const
-> dereference operator
Definition: accessors.hh:514
SideIter(const Side &side)
Definition: accessors.hh:495
bool operator!=(const SideIter &other)
Definition: accessors.hh:505
bool operator==(const SideIter &other)
Definition: accessors.hh:499
SideIter & operator++()
prefix increment iterate only on local element
Definition: accessors.hh:518
const Side & operator*() const
Definition: accessors.hh:510
double measure() const
Calculate metrics of the side.
Definition: accessors.cc:27
unsigned int elem_idx() const
Returns index of element in Mesh::element_vec_.
Definition: accessors.hh:448
bool operator!=(const Side &other) const
Definition: accessors.hh:466
unsigned int side_idx() const
Returns local index of the side on the element.
Definition: accessors.hh:444
const MeshBase * mesh() const
Returns pointer to the mesh.
Definition: accessors.hh:440
const MeshBase * mesh_
Pointer to Mesh to which belonged.
Definition: accessors.hh:480
bool is_valid() const
Returns true if the side has assigned element.
Definition: accessors.hh:452
Boundary cond() const
unsigned int n_nodes() const
Returns number of nodes of the side.
Definition: accessors.hh:436
arma::vec3 normal() const
Vector of (generalized) normal.
Definition: accessors.cc:48
arma::vec3 normal_triangle() const
Definition: accessors.cc:101
arma::vec3 normal_line() const
Definition: accessors.cc:80
ElementAccessor< 3 > element() const
Returns iterator to the element of the side.
double diameter() const
Calculate the side diameter.
Definition: accessors.cc:132
unsigned int cond_idx() const
Returns global index of the prescribed boundary condition.
bool is_external() const
Returns true for all sides either on boundary or connected to vb neigboring.
Side()
Default invalid side accessor constructor.
unsigned int side_idx_
Local # of side in element (to remove it, we heve to remove calc_side_rhs)
Definition: accessors.hh:482
bool is_boundary() const
Returns true for side on the boundary.
void inc()
Iterate over local sides of the element.
Definition: accessors.hh:456
unsigned int edge_idx() const
Returns global index of the edge connected to the side.
arma::vec3 centre() const
Centre of side.
Definition: accessors.cc:116
bool operator==(const Side &other) const
Definition: accessors.hh:461
unsigned int dim() const
Returns dimension of the side, that is dimension of the element minus one.
unsigned int elem_idx_
Index of element in Mesh::element_vec_.
Definition: accessors.hh:481
Edge edge() const
Returns pointer to the edge connected to the side.
NodeAccessor< 3 > node(unsigned int i) const
Returns node for given local index i on the side.
arma::vec3 normal_point() const
This is necessary by current DofHandler, should change this.
Definition: accessors.cc:64
Support classes for parallel programing.
Implementation of the inline functions of the mesh accessors.
unsigned int uint
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.