Flow123d  release_2.2.0-22-g936454a
mesh.h
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 mesh.h
15  * @brief
16  */
17 
18 #ifndef MAKE_MESH_H
19 #define MAKE_MESH_H
20 
21 #include <vector>
22 #include <mpi.h>
23 
24 #include "mesh/mesh_types.hh"
25 
26 #include "mesh/nodes.hh"
27 //#include "mesh/elements.h"
28 //#include "mesh/sides.h"
29 #include "mesh/edges.h"
30 #include "mesh/neighbours.h"
31 #include "mesh/boundaries.h"
32 #include "mesh/intersection.hh"
33 #include "mesh/partitioning.hh"
34 #include "mesh/region_set.hh"
35 
36 
39 #include "system/exceptions.hh"
40 
41 
42 
43 
44 #define ELM 0
45 #define BC 1
46 #define NODE 2
47 
48 /**
49  * This parameter limits volume of elements from below.
50  */
51 #define MESH_CRITICAL_VOLUME 1.0E-12
52 
53 /**
54  * Provides for statement to iterate over the Nodes of the Mesh.
55  * The parameter is FullIter local variable of the cycle, so it need not be declared before.
56  * Macro assume that variable Mesh *mesh; is declared and points to a valid Mesh structure.
57  */
58 #define FOR_NODES(_mesh_, i) \
59  for( NodeFullIter i( (_mesh_)->node_vector.begin() ); \
60  i != (_mesh_)->node_vector.end(); \
61  ++i)
62 
63 /**
64  * Macro for conversion form Iter to FullIter for nodes.
65  */
66 #define NODE_FULL_ITER(_mesh_,i) \
67  (_mesh_)->node_vector.full_iter(i)
68 
69 /**
70  * Macro to get "NULL" ElementFullIter.
71  */
72 #define NODE_FULL_ITER_NULL(_mesh_) \
73  NodeFullIter((_mesh_)->node_vector)
74 
75 /**
76  * Macro for conversion form Iter to FullIter for elements.
77  */
78 #define ELEM_FULL_ITER(_mesh_,i) \
79  (_mesh_)->element.full_iter(i)
80 
81 
82 #define FOR_NODE_ELEMENTS(i,j) for((j)=0;(j)<(i)->n_elements();(j)++)
83 #define FOR_NODE_SIDES(i,j) for((j)=0;(j)<(i)->n_sides;(j)++)
84 
85 
87 public:
89 };
90 
91 
92 
93 //=============================================================================
94 // STRUCTURE OF THE MESH
95 //=============================================================================
96 
97 class Mesh {
98 public:
99  TYPEDEF_ERR_INFO( EI_ElemLast, int);
100  TYPEDEF_ERR_INFO( EI_ElemNew, int);
101  TYPEDEF_ERR_INFO( EI_RegLast, std::string);
102  TYPEDEF_ERR_INFO( EI_RegNew, std::string);
103  DECLARE_EXCEPTION(ExcDuplicateBoundary,
104  << "Duplicate boundary elements! \n"
105  << "Element id: " << EI_ElemLast::val << " on region name: " << EI_RegLast::val << "\n"
106  << "Element id: " << EI_ElemNew::val << " on region name: " << EI_RegNew::val << "\n");
107 
108 
109 
110  static const unsigned int undef_idx=-1;
111  static const Input::Type::Record & get_input_type();
112 
113 
114 
115  /** Labels for coordinate indexes in arma::vec3 representing vectors and points.*/
116  enum {x_coord=0, y_coord=1, z_coord=2};
117 
118  /**
119  * Empty constructor.
120  *
121  * Use only for unit tests!!!
122  */
123  Mesh();
124  /**
125  * Constructor from an input record.
126  * Do not process input record. That is done in init_from_input.
127  */
128  Mesh(Input::Record in_record, MPI_Comm com = MPI_COMM_WORLD);
129  /**
130  * Common part of both previous constructors and way how to reinitialize a mesh from the given input record.
131  */
132  void reinit(Input::Record in_record);
133 
134  /// Destructor.
135  ~Mesh();
136 
137  inline unsigned int n_nodes() const {
138  return node_vector.size();
139  }
140 
141  inline unsigned int n_elements() const {
142  return element.size();
143  }
144 
145  inline unsigned int n_boundaries() const {
146  return boundary_.size();
147  }
148 
149  inline unsigned int n_edges() const {
150  return edges.size();
151  }
152 
153  unsigned int n_corners();
154 
155  inline const RegionDB &region_db() const {
156  return region_db_;
157  }
158 
159  /// Reserve size of node vector
160  inline void reserve_node_size(unsigned int n_nodes) {
161  node_vector.reserve(n_nodes);
162  }
163 
164  /// Reserve size of element vector
165  inline void reserve_element_size(unsigned int n_elements) {
166  element.reserve(n_elements);
167  }
168 
169  /**
170  * Returns pointer to partitioning object. Partitioning is created during setup_topology.
171  */
172  Partitioning *get_part();
173 
175  { return el_ds; }
176 
177  int *get_row_4_el() const
178  { return row_4_el; }
179 
180  int *get_el_4_loc() const
181  { return el_4_loc; }
182 
183  /**
184  * Returns MPI communicator of the mesh.
185  */
186  inline MPI_Comm get_comm() const { return comm_; }
187 
188 
189  void make_intersec_elements();
190 
191  unsigned int n_sides();
192 
193  inline unsigned int n_vb_neighbours() const {
194  return vb_neighbours_.size();
195  }
196 
197  /**
198  * Returns maximal number of sides of one edge, which connects elements of dimension @p dim.
199  * @param dim Dimension of elements sharing the edge.
200  */
201  unsigned int max_edge_sides(unsigned int dim) const { return max_edge_sides_[dim-1]; }
202 
203 
204  /**
205  * Fills vectors of ID numbers of bulk and bc elements.
206  */
207  void elements_id_maps( vector<int> & bulk_elements_id, vector<int> & boundary_elements_id) const;
208 
209 
210  ElementAccessor<3> element_accessor(unsigned int idx, bool boundary=false);
211 
212  /**
213  * Reads elements and their affiliation to regions and region sets defined by user in input file
214  * Format of input record is defined in method RegionSetBase::get_input_type()
215  *
216  * @param region_list Array input AbstractRecords which define regions, region sets and elements
217  */
218  void read_regions_from_input(Input::Array region_list);
219 
220  /**
221  * Returns nodes_elements vector, if doesn't exist creates its.
222  */
223  vector<vector<unsigned int> > const & node_elements();
224 
225  /// Vector of nodes of the mesh.
227  /// Vector of elements of the mesh.
229 
230  /// Vector of boundary sides where is prescribed boundary condition.
231  /// TODO: apply all boundary conditions in the main assembling cycle over elements and remove this Vector.
233  /// vector of boundary elements - should replace 'boundary'
234  /// TODO: put both bulk and bc elements (on zero level) to the same vector or make better map id->element for field inputs that use element IDs
235  /// the avoid usage of ElementVector etc.
237 
238  /// Vector of MH edges, this should not be part of the geometrical mesh
240 
241  //flow::VectorId<int> bcd_group_id; // gives a index of group for an id
242 
243  /**
244  * Vector of individual intersections of two elements.
245  * This is enough for local mortar.
246  */
248 
249  /**
250  * For every element El we have vector of indices into @var intersections array for every intersection in which El is master element.
251  * This is necessary for true mortar.
252  */
254 
255  /**
256  * Vector of compatible neighbourings.
257  */
259 
260  int n_insides; // # of internal sides
261  int n_exsides; // # of external sides
262  int n_sides_; // total number of sides (should be easy to count when we have separated dimensions
263 
264  int n_lines; // Number of line elements
265  int n_triangles; // Number of triangle elements
266  int n_tetrahedras; // Number of tetrahedra elements
267 
268  // Temporary solution for numbering of nodes on sides.
269  // The data are defined in RefElement<dim>::side_nodes,
270  // Mesh::side_nodes can be removed as soon as Element
271  // is templated by dimension.
272  //
273  // for every side dimension D = 0 .. 2
274  // for every element side 0 .. D+1
275  // for every side node 0 .. D
276  // index into element node array
278 
279  /**
280  * Initialize all mesh structures from raw information about nodes and elements (including boundary elements).
281  * Namely: create remaining boundary elements and Boundary objects, find edges and compatible neighborings.
282  */
283  void setup_topology();
284 
285  /**
286  * Check usage of regions, set regions to elements defined by user, close RegionDB
287  */
288  void check_and_finish();
289 
290  const BIHTree &get_bih_tree();
291 
292  /**
293  * Find intersection of element lists given by Mesh::node_elements_ for elements givne by @p nodes_list parameter.
294  * The result is placed into vector @p intersection_element_list. If the @p node_list is empty, and empty intersection is
295  * returned.
296  */
297  void intersect_element_lists(vector<unsigned int> const &nodes_list, vector<unsigned int> &intersection_element_list);
298 
299  /// Add new node of given id and coordinates to mesh
300  void add_node(unsigned int node_id, arma::vec3 coords);
301 
302  /// Add new element of given id to mesh
303  void add_element(unsigned int elm_id, unsigned int dim, unsigned int region_id, unsigned int partition_id,
304  std::vector<unsigned int> node_ids);
305 
306  /// Add new node of given id and coordinates to mesh
307  void add_physical_name(unsigned int dim, unsigned int id, std::string name);
308 
309  /// Return FilePath object representing "mesh_file" input key
310  inline FilePath mesh_file() {
311  return in_record_.val<FilePath>("mesh_file");
312  }
313 
314  /// Number of elements read from input.
315  unsigned int n_all_input_elements_;
316 
317 protected:
318 
319  /**
320  * This replaces read_neighbours() in order to avoid using NGH preprocessor.
321  *
322  * TODO:
323  * - Avoid maps:
324  *
325  * 4) replace EdgeVector by std::vector<Edge> (need not to know the size)
326  *
327  * 5) need not to have temporary array for Edges, only postpone setting pointers in elements and set them
328  * after edges are found; we can temporary save Edge index instead of pointer in Neigbours and elements
329  *
330  * 6) Try replace Edge * by indexes in Neigbours and elements (anyway we have mesh pointer in elements so it is accessible also from Neigbours)
331  *
332  */
333  void make_neighbours_and_edges();
334 
335  /**
336  * On edges sharing sides of many elements it may happen that each side has its nodes ordered in a different way.
337  * This method finds the permutation for each side so as to obtain the ordering of side 0.
338  */
339  void make_edge_permutations();
340  /**
341  * Create element lists for nodes in Mesh::nodes_elements.
342  */
343  void create_node_element_lists();
344  /**
345  * Remove elements with dimension not equal to @p dim from @p element_list. Index of the first element of dimension @p dim-1,
346  * is returned in @p element_idx. If no such element is found the method returns false, if one such element is found the method returns true,
347  * if more elements are found we report an user input error.
348  */
349  bool find_lower_dim_element(ElementVector&elements, vector<unsigned int> &element_list, unsigned int dim, unsigned int &element_idx);
350 
351  /**
352  * Returns true if side @p si has same nodes as in the list @p side_nodes.
353  */
354  bool same_sides(const SideIter &si, vector<unsigned int> &side_nodes);
355 
356 
357  void element_to_neigh_vb();
358 
359  void count_element_types();
360  void count_side_types();
361 
362  /**
363  * Possibly modify region id of elements sets by user in "regions" part of input file.
364  *
365  * TODO: This method needs check in issue 'Review mesh setting'.
366  * Changes have been done during generalized region key and may be causing problems
367  * during the further development.
368  */
369  void modify_element_ids(const RegionDB::MapElementIDToRegionID &map);
370 
371  unsigned int n_bb_neigh, n_vb_neigh;
372 
373  /// Maximal number of sides per one edge in the actual mesh (set in make_neighbours_and_edges()).
374  unsigned int max_edge_sides_[3];
375 
376  /**
377  * Database of regions (both bulk and boundary) of the mesh. Regions are logical parts of the
378  * domain that allows setting of different data and boundary conditions on them.
379  */
381  /**
382  * Mesh partitioning. Created in setup_topology.
383  */
384  std::shared_ptr<Partitioning> part_;
385 
386  /**
387  * BIH Tree for intersection and observe points lookup.
388  */
389  std::shared_ptr<BIHTree> bih_tree_;
390  /**
391  * Accessor to the input record for the mesh.
392  */
394 
395  /**
396  * MPI communicator used for partitioning and ...
397  */
399 
400  // For each node the vector contains a list of elements that use this node
402 
403  friend class RegionSetBase;
404 
405 
406 private:
407 
408  /// Index set assigning to global element index the local index used in parallel vectors.
409  int *row_4_el;
410  /// Index set assigning to local element index its global index.
411  int *el_4_loc;
412  /// Parallel distribution of elements.
414 };
415 
416 
417 #include "mesh/side_impl.hh"
418 #include "mesh/element_impls.hh"
419 #include "mesh/neighbours_impl.hh"
420 
421 /**
422  * Provides for statement to iterate over the Elements of the Mesh.
423  * The parameter is FullIter local variable of the cycle, so it need not be declared before.
424  * Macro assume that variable Mesh *mesh; is declared and points to a valid Mesh structure.
425  */
426 #define FOR_ELEMENTS(_mesh_,__i) \
427  for( ElementFullIter __i( (_mesh_)->element.begin() ); \
428  __i != (_mesh_)->element.end(); \
429  ++__i)
430 
431 /**
432  * Macro for conversion form Iter to FullIter for elements.
433  */
434 #define ELEMENT_FULL_ITER(_mesh_,i) \
435  (_mesh_)->element.full_iter(i)
436 
437 /**
438  * Macro to get "NULL" ElementFullIter.
439  */
440 #define ELEMENT_FULL_ITER_NULL(_mesh_) \
441  ElementFullIter((_mesh_)->element)
442 
443 
444 #define FOR_BOUNDARIES(_mesh_,i) \
445 for( std::vector<Boundary>::iterator i= (_mesh_)->boundary_.begin(); \
446  i != (_mesh_)->boundary_.end(); \
447  ++i)
448 
449 /**
450  * Macro for conversion form Iter to FullIter for boundaries.
451  */
452 #define BOUNDARY_FULL_ITER(_mesh_,i) \
453  (_mesh_)->boundary.full_iter(i)
454 
455 /**
456  * Macro to get "NULL" BoundaryFullIter.
457  */
458 #define BOUNDARY_NULL(_mesh_) \
459  BoundaryFullIter((_mesh_)->boundary)
460 
461 
462 /**
463  * Provides for statement to iterate over the Edges of the Mesh. see FOR_ELEMENTS
464  */
465 #define FOR_EDGES(_mesh_,__i) \
466  for( vector<Edge>::iterator __i = (_mesh_)->edges.begin(); \
467  __i !=(_mesh_)->edges.end(); \
468  ++__i)
469 
470 #define FOR_SIDES(_mesh_, it) \
471  FOR_ELEMENTS((_mesh_), ele) \
472  for(SideIter it = ele->side(0); it->el_idx() < ele->n_sides(); ++it)
473 
474 #define FOR_SIDE_NODES(i,j) for((j)=0;(j)<(i)->n_nodes;(j)++)
475 
476 
477 #define FOR_NEIGHBOURS(_mesh_, it) \
478  for( std::vector<Neighbour>::iterator it = (_mesh_)->vb_neighbours_.begin(); \
479  (it)!= (_mesh_)->vb_neighbours_.end(); ++it)
480 
481 #define FOR_NEIGH_ELEMENTS(i,j) for((j)=0;(j)<(i)->n_elements;(j)++)
482 #define FOR_NEIGH_SIDES(i,j) for((j)=0;(j)<(i)->n_sides;(j)++)
483 
484 
485 #endif
486 //-----------------------------------------------------------------------------
487 // vim: set cindent:
int n_triangles
Definition: mesh.h:265
Distribution * el_ds
Parallel distribution of elements.
Definition: mesh.h:413
Class for the mesh partitioning. This should provide:
Definition: partitioning.hh:39
Accessor to input data conforming to declared Array.
Definition: accessors.hh:567
int MPI_Comm
Definition: mpi.h:141
Nodes of a mesh.
static Input::Type::Record input_type
Definition: mesh.h:88
int n_lines
Definition: mesh.h:264
#define DECLARE_EXCEPTION(ExcName, Format)
Macro for simple definition of exceptions.
Definition: exceptions.hh:158
Definition: mesh.h:97
Input::Record in_record_
Definition: mesh.h:393
int * get_el_4_loc() const
Definition: mesh.h:180
vector< vector< vector< unsigned int > > > side_nodes
Definition: mesh.h:277
int n_tetrahedras
Definition: mesh.h:266
const RegionDB & region_db() const
Definition: mesh.h:155
vector< Boundary > boundary_
Definition: mesh.h:232
int * el_4_loc
Index set assigning to local element index its global index.
Definition: mesh.h:411
int * row_4_el
Index set assigning to global element index the local index used in parallel vectors.
Definition: mesh.h:409
int n_exsides
Definition: mesh.h:261
unsigned int n_vb_neighbours() const
Definition: mesh.h:193
unsigned int n_elements() const
Definition: mesh.h:141
MPI_Comm get_comm() const
Definition: mesh.h:186
unsigned int n_boundaries() const
Definition: mesh.h:145
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
Class for O(log N) lookup for intersections with a set of bounding boxes.
Definition: bih_tree.hh:36
ElementVector bc_elements
Definition: mesh.h:236
void reserve_element_size(unsigned int n_elements)
Reserve size of element vector.
Definition: mesh.h:165
#define TYPEDEF_ERR_INFO(EI_Type, Type)
Macro to simplify declaration of error_info types.
Definition: exceptions.hh:194
std::shared_ptr< BIHTree > bih_tree_
Definition: mesh.h:389
Distribution * get_el_ds() const
Definition: mesh.h:174
vector< vector< unsigned int > > node_elements_
Definition: mesh.h:401
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
vector< vector< unsigned int > > master_elements
Definition: mesh.h:253
vector< Neighbour > vb_neighbours_
Definition: mesh.h:258
unsigned int n_vb_neigh
Definition: mesh.h:371
unsigned int n_all_input_elements_
Number of elements read from input.
Definition: mesh.h:315
int n_sides_
Definition: mesh.h:262
unsigned int max_edge_sides(unsigned int dim) const
Definition: mesh.h:201
FilePath mesh_file()
Return FilePath object representing "mesh_file" input key.
Definition: mesh.h:310
int * get_row_4_el() const
Definition: mesh.h:177
RegionDB region_db_
Definition: mesh.h:380
std::vector< Edge > edges
Vector of MH edges, this should not be part of the geometrical mesh.
Definition: mesh.h:239
vector< Intersection > intersections
Definition: mesh.h:247
MPI_Comm comm_
Definition: mesh.h:398
int n_insides
Definition: mesh.h:260
#define MPI_COMM_WORLD
Definition: mpi.h:123
unsigned int n_nodes() const
Definition: mesh.h:137
Record type proxy class.
Definition: type_record.hh:182
unsigned int n_edges() const
Definition: mesh.h:149
void reserve_node_size(unsigned int n_nodes)
Reserve size of node vector.
Definition: mesh.h:160
std::shared_ptr< Partitioning > part_
Definition: mesh.h:384
NodeVector node_vector
Vector of nodes of the mesh.
Definition: mesh.h:226
ElementVector element
Vector of elements of the mesh.
Definition: mesh.h:228