Flow123d  JS_before_hm-1804-gf2ad740aa
mesh.cc
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.cc
15  * @ingroup mesh
16  * @brief Mesh construction
17  */
18 
19 
20 #include <unistd.h>
21 #include <set>
22 #include <unordered_map>
23 
24 #include "system/system.hh"
25 #include "system/exceptions.hh"
26 #include "system/index_types.hh"
28 #include "input/input_type.hh"
29 #include "input/accessors.hh"
30 #include "system/sys_profiler.hh"
31 #include "la/distribution.hh"
32 
33 #include "mesh/mesh.h"
34 #include "mesh/bc_mesh.hh"
35 #include "mesh/ref_element.hh"
36 #include "mesh/region_set.hh"
37 #include "mesh/range_wrapper.hh"
38 
39 // think about following dependencies
40 #include "mesh/accessors.hh"
41 #include "mesh/node_accessor.hh"
42 #include "mesh/partitioning.hh"
43 #include "mesh/neighbours.h"
44 
45 
46 #include "mesh/bih_tree.hh"
47 #include "mesh/duplicate_nodes.h"
48 #include "mesh/mesh_optimizer.hh"
49 
51 
52 
53 
54 //TODO: sources, concentrations, initial condition and similarly boundary conditions should be
55 // instances of a Element valued field
56 // concentrations is in fact reimplemented in transport REMOVE it HERE
57 
58 // After removing non-geometrical things from mesh, this should be part of mash initializing.
59 #include "mesh/region.hh"
60 
61 #define NDEF -1
62 
63 namespace IT = Input::Type;
64 
66  return Input::Type::Selection("Types of search algorithm for finding intersection candidates.")
67  .add_value(Mesh::BIHsearch, "BIHsearch",
68  "Use BIH for finding initial candidates, then continue by prolongation.")
69  .add_value(Mesh::BIHonly, "BIHonly",
70  "Use BIH for finding all candidates.")
71  .add_value(Mesh::BBsearch, "BBsearch",
72  "Use bounding boxes for finding initial candidates, then continue by prolongation.")
73  .close();
74 }
75 
77  return IT::Record("Mesh","Record with mesh related data." )
78  .allow_auto_conversion("mesh_file")
80  "Input file with mesh description.")
82  "List of additional region and region set definitions not contained in the mesh. "
83  "There are three region sets implicitly defined:\n\n"
84  "- ALL (all regions of the mesh)\n"
85  "- .BOUNDARY (all boundary regions)\n"
86  "- BULK (all bulk regions)")
87  .declare_key("partitioning", Partitioning::get_input_type(), IT::Default("\"any_neighboring\""), "Parameters of mesh partitioning algorithms.\n" )
88  .declare_key("print_regions", IT::Bool(), IT::Default("true"), "If true, print table of all used regions.")
89  .declare_key("intersection_search", Mesh::get_input_intersection_variant(),
90  IT::Default("\"BIHsearch\""), "Search algorithm for element intersections.")
91  .declare_key("global_snap_radius", IT::Double(0.0), IT::Default("1E-3"),
92  "Maximal snapping distance from the mesh in various search operations. In particular, it is used "
93  "to find the closest mesh element of an observe point; and in FieldFormula to find closest surface "
94  "element in plan view (Z projection).")
96  "Output file with neighboring data from mesh.")
97  .declare_key("optimize_mesh", IT::Bool(), IT::Default("true"), "If true, permute nodes and elements in order to increase cache locality. "
98  "This will speed up the calculations. GMSH output preserves original ordering but is slower. All variants of VTK output use the permuted.")
99  .close();
100 }
101 
103 : tree(nullptr),
104  comm_(MPI_COMM_WORLD),
105  bulk_size_(0),
106  nodes_(3, 1, 0),
107  row_4_el(nullptr),
108  el_4_loc(nullptr),
109  el_ds(nullptr),
110  node_4_loc_(nullptr),
111  node_ds_(nullptr),
112  bc_mesh_(nullptr)
113 
114 {init();}
115 
116 
117 
119 : tree(nullptr),
120  optimize_memory_locality(true),
121  in_record_(in_record),
122  comm_(com),
123  bulk_size_(0),
124  nodes_(3, 1, 0),
125  row_4_el(nullptr),
126  el_4_loc(nullptr),
127  el_ds(nullptr),
128  node_4_loc_(nullptr),
129  node_ds_(nullptr),
130  bc_mesh_(nullptr)
131 {
132 
133  init();
134 }
135 
136 
138  : tree(nullptr),
139  optimize_memory_locality(other.optimize_memory_locality),
140  in_record_(other.in_record_),
141  comm_(other.comm_),
142  bulk_size_(0),
143  nodes_(3, 1, 0),
144  row_4_el(nullptr),
145  el_4_loc(nullptr),
146  el_ds(nullptr),
147  node_4_loc_(nullptr),
148  node_ds_(nullptr),
149  bc_mesh_(nullptr)
150 {
151  init();
152 }
153 
154 
155 
157 {
158  return in_record_.val<Mesh::IntersectionSearch>("intersection_search");
159 }
160 
161 
163 {
164  // set in_record_, if input accessor is empty
165  if (in_record_.is_empty()) {
166  istringstream is("{mesh_file=\"\", optimize_mesh=false}");
167  Input::ReaderToStorage reader;
168  IT::Record &in_rec = const_cast<IT::Record &>(Mesh::get_input_type());
169  in_rec.finish();
170  reader.read_stream(is, in_rec, Input::FileFormat::format_JSON);
172  }
173 
174  optimize_memory_locality = in_record_.val<bool>("optimize_mesh");
175 
176  n_insides = NDEF;
177  n_exsides = NDEF;
178  n_sides_ = NDEF;
179 
180  for (int d=0; d<3; d++) max_edge_sides_[d] = 0;
181 
182  // Initialize numbering of nodes on sides.
183  // This is temporary solution, until class Element is templated
184  // by dimension. Then we can replace Mesh::side_nodes by
185  // RefElement<dim>::side_nodes.
186 
187  // indices of side nodes in element node array
188  // Currently this is made ad libitum
189  // with some ordering here we can get sides with correct orientation.
190  // This speedup normal calculation.
191 
192  side_nodes.resize(3); // three side dimensions
193  for(int dim=0; dim < 3; dim++) {
194  side_nodes[dim].resize(dim+2); // number of sides
195  for(int i_side=0; i_side < dim+2; i_side++)
196  side_nodes[dim][i_side].resize(dim+1);
197  }
198 
199  for (unsigned int sid=0; sid<RefElement<1>::n_sides; sid++)
200  for (unsigned int nid=0; nid<RefElement<1>::n_nodes_per_side; nid++)
201  side_nodes[0][sid][nid] = RefElement<1>::interact(Interaction<0,0>(sid))[nid];
202 
203  for (unsigned int sid=0; sid<RefElement<2>::n_sides; sid++)
204  for (unsigned int nid=0; nid<RefElement<2>::n_nodes_per_side; nid++)
205  side_nodes[1][sid][nid] = RefElement<2>::interact(Interaction<0,1>(sid))[nid];
206 
207  for (unsigned int sid=0; sid<RefElement<3>::n_sides; sid++)
208  for (unsigned int nid=0; nid<RefElement<3>::n_nodes_per_side; nid++)
209  side_nodes[2][sid][nid] = RefElement<3>::interact(Interaction<0,2>(sid))[nid];
210 }
211 
212 
214  for(EdgeData &edg : this->edges)
215  if (edg.side_) delete[] edg.side_;
216 
217  for (unsigned int idx=0; idx < bulk_size_; idx++) {
218  Element *ele=&(element_vec_[idx]);
219  if (ele->boundary_idx_) delete[] ele->boundary_idx_;
220  if (ele->neigh_vb) delete[] ele->neigh_vb;
221  }
222 
223  for(unsigned int idx=bulk_size_; idx < element_vec_.size(); idx++) {
224  Element *ele=&(element_vec_[idx]);
225  if (ele->boundary_idx_) delete[] ele->boundary_idx_;
226  }
227 
228  if (row_4_el != nullptr) delete[] row_4_el;
229  if (el_4_loc != nullptr) delete[] el_4_loc;
230  if (el_ds != nullptr) delete el_ds;
231  if (node_4_loc_ != nullptr) delete[] node_4_loc_;
232  if (node_ds_ != nullptr) delete node_ds_;
233  if (bc_mesh_ != nullptr) delete bc_mesh_;
234  if (tree != nullptr) delete tree;
235 }
236 
237 
238 unsigned int Mesh::n_sides() const
239 {
240  if (n_sides_ == NDEF) {
241  n_sides_=0;
242  for (auto ele : this->elements_range()) n_sides_ += ele->n_sides();
243  }
244  return n_sides_;
245 }
246 
247 unsigned int Mesh::n_vb_neighbours() const {
248  return vb_neighbours_.size();
249  }
250 
251 
252 unsigned int Mesh::n_corners() {
253  unsigned int li, count = 0;
254  for (auto ele : this->elements_range()) {
255  for (li=0; li<ele->n_nodes(); li++) {
256  count++;
257  }
258  }
259  return count;
260 }
261 
262 Edge Mesh::edge(uint edge_idx) const
263 {
264  ASSERT_LT_DBG(edge_idx, edges.size());
265  return Edge(this, edge_idx);
266 }
267 
269 {
270  ASSERT_LT_DBG(bc_idx, boundary_.size());
271  return Boundary(&boundary_[bc_idx]);
272 }
273 
275  return part_.get();
276 }
277 
279  return (LongIdx*)this->get_part()->get_loc_part();
280 }
281 
282 
283 
284 
286 
287  // get dim of the first element in the map, if it exists
288  uint dim_to_check = RegionDB::undefined_dim;
289  std::string reg_name = "UndefinedRegion";
290  if(map.size() > 0){
291  Element &ele = element_vec_[ elem_index(map.begin()->first) ];
292  dim_to_check = ele.dim();
293  reg_name = region_db_.find_id(map.begin()->second).label();
294  }
295 
296  for (auto elem_to_region : map) {
297  Element &ele = element_vec_[ elem_index(elem_to_region.first) ];
298 
299  if( ele.dim() != dim_to_check){
300  THROW(ExcRegionElmDiffDim() << EI_Region(reg_name) << EI_RegIdx(elem_to_region.second) << EI_Dim(dim_to_check)
301  << EI_DimOther(ele.dim()) << EI_ElemId(elem_to_region.first) );
302  }
303 
304  ele.region_idx_ = region_db_.get_region( elem_to_region.second, ele.dim() );
306  }
307 }
308 
309 
311  std::vector<uint> nodes_new_idx( this->n_nodes(), undef_idx );
312 
313  // check element quality and flag used nodes
314  for (auto ele : this->elements_range()) {
315  // element quality
316  double quality = ele.quality_measure_smooth();
317  if (quality < 0) {
318  ASSERT_LT_DBG(ele.jacobian_S3(), 0);
319  element_vec_[ele.mesh_idx()].inverted = true;
320  quality = -quality;
321  }
322  if (quality < 4*std::numeric_limits<double>::epsilon())
323  THROW( ExcBadElement() << EI_Quality(quality) << EI_ElemId(ele.idx()) );
324  if ( quality< 0.001)
325  WarningOut().fmt("Bad quality element ID={}, ({}<0.001).\n", ele.idx(), quality);
326 
327  // flag used nodes
328  for (uint ele_node=0; ele_node<ele->n_nodes(); ele_node++) {
329  uint inode = ele->node_idx(ele_node);
330  nodes_new_idx[inode] = inode;
331  }
332  }
333 
334  // possibly build new node ids map
335  BidirectionalMap<int> new_node_ids_;
336  new_node_ids_.reserve(node_ids_.size());
337 
338  // remove unused nodes from the mesh
339  uint inode_new = 0;
340  for(uint inode = 0; inode < nodes_new_idx.size(); inode++) {
341  if(nodes_new_idx[inode] == undef_idx){
342  WarningOut().fmt("A node {} does not belong to any element "
343  " and will be removed.",
344  find_node_id(inode));
345  }
346  else{
347  // map new node numbering
348  nodes_new_idx[inode] = inode_new;
349 
350  // possibly move the nodes
351  nodes_.vec<3>(inode_new) = nodes_.vec<3>(inode);
352  new_node_ids_.add_item(node_ids_[inode]);
353 
354  inode_new++;
355  }
356  }
357 
358  uint n_nodes_new = inode_new;
359 
360  // if some node erased, update node ids in elements
361  if(n_nodes_new < nodes_new_idx.size()){
362 
363  DebugOut() << "Updating node-element numbering due to unused nodes: "
364  << print_var(n_nodes_new) << print_var(nodes_new_idx.size()) << "\n";
365 
366  // throw away unused nodes
367  nodes_.resize(n_nodes_new);
368  node_ids_ = new_node_ids_;
369 
370  // update node-element numbering
371  for (auto ele : this->elements_range()) {
372  for (uint ele_node=0; ele_node<ele->n_nodes(); ele_node++) {
373  uint inode_orig = ele->node_idx(ele_node);
374  uint inode = nodes_new_idx[inode_orig];
375  ASSERT_DBG(inode != undef_idx);
376  const_cast<Element*>(ele.element())->nodes_[ele_node] = inode;
377  }
378  }
379  }
380 }
381 
382 
383 //void Mesh::array_sort(std::array<uint, 4> &nodes) {
384 // // TODO: use templated insert sort with recursion over length of array so that compiler can
385 // // optimize for the small array size.
386 //
387 // std::sort(nodes.begin(), nodes.end());
388 //}
389 
391  // element_vec_ still contains both bulk and boundary elements
392  for (uint i_el=0; i_el < element_vec_.size(); i_el++) {
393  Element &ele = element_vec_[i_el];
394  std::sort(ele.nodes_.begin(), ele.nodes_.end());
395  }
396 
397 }
398 
401  START_TIMER("MESH - optimizer");
402  this->optimize();
403  END_TIMER("MESH - optimizer");
404  }
405 
406  START_TIMER("MESH - setup topology");
407 
408  canonical_faces();
410 
411 
415 
416  tree = new DuplicateNodes(this);
417 
418  part_ = std::make_shared<Partitioning>(this, in_record_.val<Input::Record>("partitioning") );
419 
420  // create parallel distribution and numbering of elements
421  LongIdx *id_4_old = new LongIdx[n_elements()];
422  int i = 0;
423  for (auto ele : this->elements_range())
424  id_4_old[i++] = ele.idx();
425  part_->id_maps(n_elements(), id_4_old, el_ds, el_4_loc, row_4_el);
426 
427  delete[] id_4_old;
428 
429  this->distribute_nodes();
430 
432 }
433 
434 
436  MeshOptimizer<3> mo(this);
437  mo.calculate_sizes();
440 
441  this->sort_permuted_nodes_elements( mo.sort_nodes(this->node_permutation_), mo.sort_elements(this->elem_permutation_) );
442 }
443 
444 
446  BidirectionalMap<int> node_ids_backup = this->node_ids_;
447  this->node_ids_.clear();
448  this->node_ids_.reserve(this->n_nodes());
449  Armor::Array<double> nodes_backup = this->nodes_;
450  for (uint i = 0; i < this->element_vec_.size(); ++i) {
451  for (uint j = 0; j < this->element_vec_[i].dim() + 1; ++j) {
452  this->element_vec_[i].nodes_[j] = this->node_permutation_[this->element_vec_[i].nodes_[j]];
453  }
454  }
455  for (uint i = 0; i < this->n_nodes(); ++i) {
456  this->nodes_.set(node_permutation_[i]) = nodes_backup.vec<3>(i);
457  this->node_ids_.add_item( node_ids_backup[new_node_ids[i]] );
458  }
459 
460  BidirectionalMap<int> elem_ids_backup = this->element_ids_;
461  this->element_ids_.clear();
463  std::vector<Element> elements_backup = this->element_vec_;
464  for (uint i = 0; i < bulk_size_; ++i) {
465  this->element_vec_[elem_permutation_[i]] = elements_backup[i];
466  this->element_ids_.add_item( elem_ids_backup[new_elem_ids[i]] );
467  }
468 }
469 
470 
471 //
473 {
474 
475  n_insides = 0;
476  n_exsides = 0;
477  for (auto ele : this->elements_range())
478  for(SideIter sde = ele.side(0); sde->side_idx() < ele->n_sides(); ++sde) {
479  if (sde->is_external()) n_exsides++;
480  else n_insides++;
481  }
482 }
483 
484 
485 
487  // for each node we make a list of elements that use this node
488  node_elements_.resize( this->n_nodes() );
489 
490  for (auto ele : this->elements_range())
491  for (unsigned int n=0; n<ele->n_nodes(); n++)
492  node_elements_[ele.node(n).idx()].push_back(ele.idx());
493 
494  for (vector<vector<unsigned int> >::iterator n=node_elements_.begin(); n!=node_elements_.end(); n++)
495  stable_sort(n->begin(), n->end());
496 }
497 
498 
499 void Mesh::intersect_element_lists(vector<unsigned int> const &nodes_list, vector<unsigned int> &intersection_element_list)
500 {
501  if (node_elements_.size() == 0) {
503  }
504 
505  if (nodes_list.size() == 0) {
506  intersection_element_list.clear();
507  } else if (nodes_list.size() == 1) {
508  intersection_element_list = node_elements_[ nodes_list[0] ];
509  } else {
510  vector<unsigned int>::const_iterator it1=nodes_list.begin();
512  intersection_element_list.resize( node_elements_[*it1].size() ); // make enough space
513 
514  it1=set_intersection(
515  node_elements_[*it1].begin(), node_elements_[*it1].end(),
516  node_elements_[*it2].begin(), node_elements_[*it2].end(),
517  intersection_element_list.begin());
518  intersection_element_list.resize(it1-intersection_element_list.begin()); // resize to true size
519 
520  for(;it2<nodes_list.end();++it2) {
521  it1=set_intersection(
522  intersection_element_list.begin(), intersection_element_list.end(),
523  node_elements_[*it2].begin(), node_elements_[*it2].end(),
524  intersection_element_list.begin());
525  intersection_element_list.resize(it1-intersection_element_list.begin()); // resize to true size
526  }
527  }
528 }
529 
530 
531 bool Mesh::find_lower_dim_element( vector<unsigned int> &element_list, unsigned int dim, unsigned int &element_idx) {
532  bool is_neighbour = false;
533 
534  vector<unsigned int>::iterator e_dest=element_list.begin();
535  for( vector<unsigned int>::iterator ele = element_list.begin(); ele!=element_list.end(); ++ele) {
536  //DebugOut() << "Eid: " << this->elem_index(*ele)
537  // << format(element_vec_[*ele].nodes_);
538 
539  if (element_vec_[*ele].dim() == dim) { // keep only indexes of elements of same dimension
540  *e_dest=*ele;
541  ++e_dest;
542  } else if (element_vec_[*ele].dim() == dim-1) { // get only first element of lower dimension
543  if (is_neighbour) THROW(ExcTooMatchingIds() << EI_ElemId(this->elem_index(*ele)) << EI_ElemIdOther(this->elem_index(element_idx)) );
544 
545  is_neighbour = true;
546  element_idx = *ele;
547  }
548  }
549  element_list.resize( e_dest - element_list.begin());
550  return is_neighbour;
551 }
552 
553 bool Mesh::same_sides(const SideIter &si, vector<unsigned int> &side_nodes) {
554  // check if nodes lists match (this is slow and will be faster only when we convert whole mesh into hierarchical design like in deal.ii)
555  unsigned int ni=0;
556  while ( ni < si->n_nodes()
557  && find(side_nodes.begin(), side_nodes.end(), si->node(ni).idx() ) != side_nodes.end() ) ni++;
558  return ( ni == si->n_nodes() );
559 }
560 
561 /**
562  * TODO:
563  * - use std::is_any for setting is_neigbour
564  * - possibly make appropriate constructors for Edge and Neighbour
565  * - check side!=-1 when searching neigbouring element
566  * - process boundary elements first, there should be no Neigh, but check it
567  * set Edge and boundary there
568  */
569 
571 {
572  ASSERT(bc_element_tmp_.size()==0)
573  .error("Temporary structure of boundary element data is not empty. Did you call create_boundary_elements?");
574 
575  Neighbour neighbour;
576  EdgeData *edg = nullptr;
577  unsigned int ngh_element_idx;
578  unsigned int last_edge_idx = undef_idx;
579 
580  neighbour.mesh_ = this;
581 
583 
584  // pointers to created edges
585  //vector<Edge *> tmp_edges;
586  edges.resize(0); // be sure that edges are empty
587 
589  vector<unsigned int> intersection_list; // list of elements in intersection of node element lists
590 
591  for( unsigned int i=bulk_size_; i<element_vec_.size(); ++i) {
592 
593  ElementAccessor<3> bc_ele = this->element_accessor(i);
594  ASSERT(bc_ele.region().is_boundary());
595  // Find all elements that share this side.
596  side_nodes.resize(bc_ele->n_nodes());
597  for (unsigned n=0; n<bc_ele->n_nodes(); n++) side_nodes[n] = bc_ele->node_idx(n);
598  intersect_element_lists(side_nodes, intersection_list);
599  bool is_neighbour = find_lower_dim_element(intersection_list, bc_ele->dim() +1, ngh_element_idx);
600  if (is_neighbour) {
601  THROW( ExcBdrElemMatchRegular() << EI_ElemId(bc_ele.idx()) << EI_ElemIdOther(this->elem_index(ngh_element_idx)) );
602  } else {
603  if (intersection_list.size() == 0) {
604  // no matching dim+1 element found
605  WarningOut().fmt("Lonely boundary element, id: {}, region: {}, dimension {}.\n",
606  bc_ele.idx(), bc_ele.region().id(), bc_ele->dim());
607  continue; // skip the boundary element
608  }
609  last_edge_idx=edges.size();
610  edges.resize(last_edge_idx+1);
611  edg = &( edges.back() );
612  edg->n_sides = 0;
613  edg->side_ = new struct SideIter[ intersection_list.size() ];
614 
615  // common boundary object
616  unsigned int bdr_idx=boundary_.size();
617  boundary_.resize(bdr_idx+1);
618  BoundaryData &bdr=boundary_.back();
619  bdr.bc_ele_idx_ = i;
620  bdr.edge_idx_ = last_edge_idx;
621  bdr.mesh_=this;
622 
623  // for 1d boundaries there can be more then one 1d elements connected to the boundary element
624  // we do not detect this case later in the main search over bulk elements
625  for( vector<unsigned int>::iterator isect = intersection_list.begin(); isect!=intersection_list.end(); ++isect) {
626  ElementAccessor<3> elem = this->element_accessor(*isect);
627  for (unsigned int ecs=0; ecs<elem->n_sides(); ecs++) {
628  SideIter si = elem.side(ecs);
629  if ( same_sides( si, side_nodes) ) {
630  if (elem->edge_idx(ecs) != undef_idx) {
631  OLD_ASSERT(elem->boundary_idx_!=nullptr, "Null boundary idx array.\n");
632  int last_bc_ele_idx=this->boundary_[elem->boundary_idx_[ecs]].bc_ele_idx_;
633  int new_bc_ele_idx=i;
634  THROW( ExcDuplicateBoundary()
635  << EI_ElemLast(this->find_elem_id(last_bc_ele_idx))
636  << EI_RegLast(this->element_accessor(last_bc_ele_idx).region().label())
637  << EI_ElemNew(this->find_elem_id(new_bc_ele_idx))
638  << EI_RegNew(this->element_accessor(new_bc_ele_idx).region().label())
639  );
640  }
641  element_vec_[*isect].edge_idx_[ecs] = last_edge_idx;
642  edg->side_[ edg->n_sides++ ] = si;
643 
644  if (elem->boundary_idx_ == NULL) {
645  Element *el = &(element_vec_[*isect]);
646  el->boundary_idx_ = new unsigned int [ el->n_sides() ];
647  std::fill( el->boundary_idx_, el->boundary_idx_ + el->n_sides(), undef_idx);
648  }
649  elem->boundary_idx_[ecs] = bdr_idx;
650  break; // next element in intersection list
651  }
652  }
653  }
654 
655  }
656 
657  }
658  // Now we go through all element sides and create edges and neighbours
659  unsigned int new_bc_elem_idx = element_vec_.size(); //Mesh_idx of new boundary element generated in following block
660  for (auto e : this->elements_range()) {
661  for (unsigned int s=0; s<e->n_sides(); s++)
662  {
663  // skip sides that were already found
664  if (e->edge_idx(s) != undef_idx) continue;
665 
666 
667  // Find all elements that share this side.
668  side_nodes.resize(e.side(s)->n_nodes());
669  for (unsigned n=0; n<e.side(s)->n_nodes(); n++) side_nodes[n] = e.side(s)->node(n).idx();
670  intersect_element_lists(side_nodes, intersection_list);
671 
672  bool is_neighbour = find_lower_dim_element(intersection_list, e->dim(), ngh_element_idx);
673 
674  if (is_neighbour) { // edge connects elements of different dimensions
675  // Initialize for the neighbour case.
676  neighbour.elem_idx_ = ngh_element_idx;
677  } else { // edge connects only elements of the same dimension
678  // Initialize for the edge case.
679  last_edge_idx=edges.size();
680  edges.resize(last_edge_idx+1);
681  edg = &( edges.back() );
682  edg->n_sides = 0;
683  edg->side_ = new struct SideIter[ intersection_list.size() ];
684  if (intersection_list.size() > max_edge_sides_[e->dim()-1])
685  max_edge_sides_[e->dim()-1] = intersection_list.size();
686 
687  if (intersection_list.size() == 1) {
688  // outer edge, create boundary object as well
689  Element &elm = element_vec_[e.idx()];
690  edg->n_sides=1;
691  edg->side_[0] = e.side(s);
692  element_vec_[e.idx()].edge_idx_[s] = last_edge_idx;
693 
694  if (e->boundary_idx_ == NULL) {
695  elm.boundary_idx_ = new unsigned int [ e->n_sides() ];
696  std::fill( elm.boundary_idx_, elm.boundary_idx_ + e->n_sides(), undef_idx);
697  }
698 
699  unsigned int bdr_idx=boundary_.size()+1; // need for VTK mesh that has no boundary elements
700  // and bulk elements are indexed from 0
701  boundary_.resize(bdr_idx+1);
702  BoundaryData &bdr=boundary_.back();
703  elm.boundary_idx_[s] = bdr_idx;
704 
705  // fill boundary element
706  Element * bc_ele = add_element_to_vector(-bdr_idx);
707  bc_ele->init(e->dim()-1, region_db_.implicit_boundary_region() );
709  for(unsigned int ni = 0; ni< side_nodes.size(); ni++) bc_ele->nodes_[ni] = side_nodes[ni];
710 
711  // fill Boundary object
712  bdr.edge_idx_ = last_edge_idx;
713  bdr.bc_ele_idx_ = new_bc_elem_idx; //elem_index(-bdr_idx);
714  bdr.mesh_=this;
715  new_bc_elem_idx++;
716 
717  continue; // next side of element e
718  }
719  }
720 
721  // go through the elements connected to the edge or neighbour
722  // setup neigbour or edge
723  for( vector<unsigned int>::iterator isect = intersection_list.begin(); isect!=intersection_list.end(); ++isect) {
724  ElementAccessor<3> elem = this->element_accessor(*isect);
725  for (unsigned int ecs=0; ecs<elem->n_sides(); ecs++) {
726  if (elem->edge_idx(ecs) != undef_idx) continue; // ??? This should not happen.
727  SideIter si = elem.side(ecs);
728  if ( same_sides( si, side_nodes) ) {
729  if (is_neighbour) {
730  // create a new edge and neighbour for this side, and element to the edge
731  last_edge_idx=edges.size();
732  edges.resize(last_edge_idx+1);
733  edg = &( edges.back() );
734  edg->n_sides = 1;
735  edg->side_ = new struct SideIter[1];
736  edg->side_[0] = si;
737  element_vec_[elem.idx()].edge_idx_[ecs] = last_edge_idx;
738 
739  neighbour.edge_idx_ = last_edge_idx;
740 
741  vb_neighbours_.push_back(neighbour); // copy neighbour with this edge setting
742  } else {
743  // connect the side to the edge, and side to the edge
744  ASSERT_PTR_DBG(edg);
745  edg->side_[ edg->n_sides++ ] = si;
746  ASSERT_DBG(last_edge_idx != undef_idx);
747  element_vec_[elem.idx()].edge_idx_[ecs] = last_edge_idx;
748  }
749  break; // next element from intersection list
750  }
751  } // search for side of other connected element
752  } // connected elements
753 
754  if (! is_neighbour)
755  ASSERT_EQ( (unsigned int) edg->n_sides, intersection_list.size())(e.index())(s).error("Missing edge sides.");
756  } // for element sides
757  } // for elements
758 
759  MessageOut().fmt( "Created {} edges and {} neighbours.\n", edges.size(), vb_neighbours_.size() );
760 }
761 
762 
763 
764 //=============================================================================
765 //
766 //=============================================================================
768 {
769 
770  //MessageOut() << "Element to neighbours of vb2 type... "/*orig verb 5*/;
771 
772  for (vector<Element>::iterator ele = element_vec_.begin(); ele!= element_vec_.begin()+bulk_size_; ++ele)
773  ele->n_neighs_vb_ =0;
774 
775  // count vb neighs per element
776  for (auto & ngh : this->vb_neighbours_) ngh.element()->n_neighs_vb_++;
777 
778  // Allocation of the array per element
779  for (vector<Element>::iterator ele = element_vec_.begin(); ele!= element_vec_.begin()+bulk_size_; ++ele)
780  if( ele->n_neighs_vb() > 0 ) {
781  ele->neigh_vb = new struct Neighbour* [ele->n_neighs_vb()];
782  ele->n_neighs_vb_=0;
783  }
784 
785  // fill
786  ElementAccessor<3> ele;
787  for (auto & ngh : this->vb_neighbours_) {
788  ele = ngh.element();
789  ele->neigh_vb[ ele->n_neighs_vb_++ ] = &ngh;
790  }
791 
792  //MessageOut() << "... O.K.\n"/*orig verb 6*/;
793 }
794 
795 
796 
797 
798 
799 
801  /* Algorithm:
802  *
803  * 1) create BIH tree
804  * 2) for every 1D, find list of candidates
805  * 3) compute intersections for 1d, store it to master_elements
806  *
807  */
808  if (! intersections) {
809  intersections = std::make_shared<MixedMeshIntersections>(this);
810  intersections->compute_intersections();
811  }
812  return *intersections;
813 }
814 
815 
816 
818  return ElementAccessor<3>(this, idx);
819 }
820 
821 
822 
823 NodeAccessor<3> Mesh::node(unsigned int idx) const {
824  return NodeAccessor<3>(this, idx);
825 }
826 
827 
828 
829 void Mesh::elements_id_maps( vector<LongIdx> & bulk_elements_id, vector<LongIdx> & boundary_elements_id) const
830 {
831  if (bulk_elements_id.size() ==0) {
833  LongIdx last_id;
834 
835  bulk_elements_id.resize(n_elements());
836  map_it = bulk_elements_id.begin();
837  last_id = -1;
838  for(unsigned int idx=0; idx < n_elements(); idx++, ++map_it) {
839  LongIdx id = this->find_elem_id(idx);
840  last_id=*map_it = id;
841  }
842  std::sort(bulk_elements_id.begin(), bulk_elements_id.end());
843 
844  boundary_elements_id.resize(element_ids_.size()-bulk_size_);
845  map_it = boundary_elements_id.begin();
846  last_id = -1;
847  for(unsigned int idx=bulk_size_; idx<element_ids_.size(); idx++, ++map_it) {
848  LongIdx id = this->find_elem_id(idx);
849  // We set ID for boundary elements created by the mesh itself to "-1"
850  // this force gmsh reader to skip all remaining entries in boundary_elements_id
851  // and thus report error for any remaining data lines
852  if (id < 0) last_id=*map_it=-1;
853  else {
854  if (last_id >= id) THROW( ExcElmWrongOrder() << EI_ElemId(id) );
855  last_id=*map_it = id;
856  }
857  }
858  }
859 }
860 
861 
862 bool compare_points(const arma::vec3 &p1, const arma::vec3 &p2) {
863  static const double point_tolerance = 1E-10;
864  return fabs(p1[0]-p2[0]) < point_tolerance
865  && fabs(p1[1]-p2[1]) < point_tolerance
866  && fabs(p1[2]-p2[2]) < point_tolerance;
867 }
868 
869 
870 std::shared_ptr<EquivalentMeshMap> Mesh::check_compatible_mesh(Mesh & input_mesh)
871 {
872  // Assumptions:
873  // - target (computational) mesh is continous
874  // - source mesh can be both continous (unique nodes) and discontinous (duplicit nodes)
875  // - at least one compatible element must be found (each mesh can be only subdomain of the other one)
876 
877  std::vector<unsigned int> node_ids; // indices map: nodes from source mesh to nodes of target mesh
878  std::shared_ptr<EquivalentMeshMap> map_ptr =
879  std::make_shared<EquivalentMeshMap>(n_elements(), get_bc_mesh()->n_elements(), (LongIdx)undef_idx);
880  // indices map: nodes from source mesh to nodes of target mesh
881 
882  {
883  // create map `node_ids` from node indices of source mesh to node indices of target mesh
884  // - to each node of source mesh there must be one node in target mesh at maximum
885  // - to each node of target mesh there can be more than one node in source mesh
886  // - iterate over nodes of source mesh, use BIH tree of target mesh to find candidate nodes
887  // - check equality of nodes by their L1 distance with tolerance
888  std::vector<unsigned int> searched_elements; // for BIH tree
889  unsigned int i_node, i_elm_node;
890  const BIHTree &bih_tree=this->get_bih_tree();
891 
892  // create nodes of mesh
893  node_ids.resize( input_mesh.n_nodes(), undef_idx );
894  for (auto nod : input_mesh.node_range()) {
895  uint found_i_node = undef_idx;
896  bih_tree.find_point(*nod, searched_elements);
897 
898  for (std::vector<unsigned int>::iterator it = searched_elements.begin(); it!=searched_elements.end(); it++) {
899  ElementAccessor<3> ele = this->element_accessor( *it );
900  for (i_node=0; i_node<ele->n_nodes(); i_node++)
901  {
902  static const double point_tolerance = 1E-10;
903  if ( arma::norm(*ele.node(i_node) - *nod, 1) < point_tolerance) {
904  i_elm_node = ele.node(i_node).idx();
905  if (found_i_node == undef_idx)
906  found_i_node = i_elm_node;
907  else if (found_i_node != i_elm_node) {
908  // duplicate nodes in target mesh - not compatible
909  return std::make_shared<EquivalentMeshMap>();
910  }
911  }
912  }
913  }
914 
915  if (found_i_node!=undef_idx)
916  node_ids[nod.idx()] = found_i_node;
917 
918  searched_elements.clear();
919  }
920  }
921 
922  unsigned int n_found = 0; // number of found equivalent elements
923  // create map for bulk elements
924  n_found += check_compatible_elements(&input_mesh, this, node_ids, map_ptr->bulk);
925  // create map for boundary elements
926  n_found += check_compatible_elements(input_mesh.get_bc_mesh(), this->get_bc_mesh(), node_ids, map_ptr->boundary);
927 
928  // no equivalent element found => mesh is not compatible
929  if (n_found==0)
930  return std::make_shared<EquivalentMeshMap>();
931  else
932  return map_ptr;
933 }
934 
935 unsigned int Mesh::check_compatible_elements(Mesh* source_mesh, Mesh* target_mesh,
936  const std::vector<unsigned int>& node_ids,
938 {
939  // create map `element_ids_map` from ele indices of source mesh to ele indices of target mesh
940  // - iterate over elements of source mesh
941  // - get adjacent nodes of target mesh using `node_ids` map
942  // - find adjacent element of target mesh using the found nodes
943 
944  std::vector<unsigned int> result_list; // list of elements with same dimension as vtk element
945  std::vector<unsigned int> node_list; // auxiliary vector of node indices of a single element
946  std::vector<unsigned int> candidate_list; // auxiliary output vector for intersect_element_lists function
947  bool valid_nodes;
948 
949  unsigned int n_found = 0; // number of found equivalent elements
950 
951  for (auto elm : source_mesh->elements_range()) {
952  valid_nodes = true;
953  for (unsigned int j=0; j<elm->n_nodes(); j++) { // iterate trough all nodes of any element
954  if (node_ids[ elm->node_idx(j) ] == undef_idx)
955  valid_nodes = false;
956  node_list.push_back( node_ids[ elm->node_idx(j) ] );
957  }
958 
959  if (valid_nodes) {
960  target_mesh->intersect_element_lists(node_list, candidate_list);
961  for (auto i_elm : candidate_list) {
962  if ( target_mesh->element_accessor(i_elm)->dim() == elm.dim() )
963  result_list.push_back(i_elm);
964  }
965  }
966 
967  if (result_list.size() == 1) {
968  map[result_list[0]] = elm.idx();
969  n_found++;
970  }
971 
972  node_list.clear();
973  result_list.clear();
974  }
975  return n_found;
976 }
977 
978 
980 {
982  it != region_list.end();
983  ++it) {
984  // constructor has side effect in the mesh - create new region or set and store them to Mesh::region_db_
985  (*it).factory< RegionSetBase, const Input::Record &, Mesh * >(*it, this);
986  }
987 }
988 
990 {
992  region_db_.el_to_reg_map_.clear();
993  region_db_.close();
995 
996  if ( in_record_.val<bool>("print_regions") ) {
997  stringstream ss;
999  MessageOut() << ss.str();
1000  }
1001 }
1002 
1003 
1005  START_TIMER("Mesh::compute_element_boxes");
1007 
1008  // make element boxes
1009  unsigned int i=0;
1010  boxes.resize(this->n_elements());
1011  for (auto element : this->elements_range()) {
1012  boxes[i] = element.bounding_box();
1013  i++;
1014  }
1015 
1016  return boxes;
1017 }
1018 
1020  if (! this->bih_tree_) {
1021  bih_tree_ = std::make_shared<BIHTree>();
1022  bih_tree_->add_boxes( this->get_element_boxes() );
1023  bih_tree_->construct();
1024  }
1025  return *bih_tree_;
1026 }
1027 
1029  return in_record_.val<double>("global_snap_radius");
1030 }
1031 
1032 void Mesh::add_physical_name(unsigned int dim, unsigned int id, std::string name) {
1033  region_db_.add_region(id, name, dim, "$PhysicalNames");
1034 }
1035 
1036 
1037 void Mesh::add_node(unsigned int node_id, arma::vec3 coords) {
1038 
1039  nodes_.append(coords);
1040  node_ids_.add_item(node_id);
1041  node_permutation_.push_back(node_permutation_.size());
1042 }
1043 
1044 
1045 void Mesh::add_element(unsigned int elm_id, unsigned int dim, unsigned int region_id, unsigned int partition_id,
1046  std::vector<unsigned int> node_ids) {
1047  RegionIdx region_idx = region_db_.get_region( region_id, dim );
1048  if ( !region_idx.is_valid() ) {
1049  region_idx = region_db_.add_region( region_id, region_db_.create_label_from_id(region_id), dim, "$Element" );
1050  }
1051  region_db_.mark_used_region(region_idx.idx());
1052 
1053  if (region_idx.is_boundary()) {
1054  bc_element_tmp_.push_back( ElementTmpData(elm_id, dim, region_idx, partition_id, node_ids) );
1055  } else {
1056  if(dim == 0 ) {
1057  WarningOut().fmt("Bulk elements of zero size(dim=0) are not supported. Element ID: {}.\n", elm_id);
1058  }
1059  else {
1060  Element *ele = add_element_to_vector(elm_id);
1061  bulk_size_++;
1062  this->init_element(ele, elm_id, dim, region_idx, partition_id, node_ids);
1063  }
1064  }
1065 }
1066 
1067 
1068 void Mesh::init_element(Element *ele, unsigned int elm_id, unsigned int dim, RegionIdx region_idx, unsigned int partition_id,
1069  std::vector<unsigned int> node_ids) {
1070  ele->init(dim, region_idx);
1071  ele->pid_ = partition_id;
1072 
1073  unsigned int ni=0;
1074  for (; ni<ele->n_nodes(); ni++) {
1075  ele->nodes_[ni] = this->node_index(node_ids[ni]);
1076  }
1077  for( ;ni < 4; ni++) ele->nodes_[ni] = undef_idx;
1078 
1079  // check that tetrahedron element is numbered correctly and is not degenerated
1080  if(ele->dim() == 3)
1081  {
1082  ElementAccessor<3> ea = this->element_accessor( this->elem_index(elm_id) );
1083  double jac = ea.jacobian_S3();
1084  if( ! (jac > 0) ) {
1085  WarningOut().fmt("Tetrahedron element with id {} has wrong numbering or is degenerated (Jacobian = {}).",elm_id, jac);
1086  }
1087  }
1088 }
1089 
1090 
1092  if (node_elements_.size() == 0) {
1093  this->create_node_element_lists();
1094  }
1095  return node_elements_;
1096 }
1097 
1098 
1099 void Mesh::init_element_vector(unsigned int size) {
1100  element_vec_.clear();
1101  element_ids_.clear();
1102  elem_permutation_.clear();
1103  element_vec_.reserve(size);
1104  element_ids_.reserve(size);
1105  elem_permutation_.reserve(size);
1106  bc_element_tmp_.clear();
1107  bc_element_tmp_.reserve(size);
1108  bulk_size_ = 0;
1110 }
1111 
1112 
1113 void Mesh::init_node_vector(unsigned int size) {
1114  nodes_.reinit(size);
1115  node_ids_.clear();
1116  node_ids_.reserve(size);
1117  node_permutation_.clear();
1118  node_permutation_.reserve(size);
1119 }
1120 
1121 
1123  element_vec_.push_back( Element() );
1124  Element * elem = &element_vec_.back(); //[element_vec_.size()-1];
1125  element_ids_.add_item((unsigned int)(id));
1126  elem_permutation_.push_back(elem_permutation_.size());
1127  return elem;
1128 }
1129 
1131  auto bgn_it = make_iter<ElementAccessor<3>>( ElementAccessor<3>(this, 0) );
1132  auto end_it = make_iter<ElementAccessor<3>>( ElementAccessor<3>(this, bulk_size_) );
1133  return Range<ElementAccessor<3>>(bgn_it, end_it);
1134 }
1135 
1137  auto bgn_it = make_iter<NodeAccessor<3>>( NodeAccessor<3>(this, 0) );
1138  auto end_it = make_iter<NodeAccessor<3>>( NodeAccessor<3>(this, n_nodes()) );
1139  return Range<NodeAccessor<3>>(bgn_it, end_it);
1140 }
1141 
1143  auto bgn_it = make_iter<Edge>( Edge(this, 0) );
1144  auto end_it = make_iter<Edge>( Edge(this, edges.size()) );
1145  return Range<Edge>(bgn_it, end_it);
1146 }
1147 
1148 inline void Mesh::check_element_size(unsigned int elem_idx) const
1149 {
1150  ASSERT(elem_idx < element_vec_.size())(elem_idx)(element_vec_.size()).error("Index of element is out of bound of element vector!");
1151 }
1152 
1153 /*
1154  * Output of internal flow data.
1155  */
1157 {
1158  START_TIMER("Mesh::output_internal_ngh_data");
1159 
1160  FilePath raw_output_file_path;
1161  if (! in_record_.opt_val("raw_ngh_output", raw_output_file_path)) return;
1162 
1163  ofstream raw_ngh_output_file;
1164  int rank;
1165  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1166  if (rank == 0) {
1167  MessageOut() << "Opening raw ngh output: " << raw_output_file_path << "\n";
1168  try {
1169  raw_output_file_path.open_stream(raw_ngh_output_file);
1170  } INPUT_CATCH(FilePath::ExcFileOpen, FilePath::EI_Address_String, (in_record_))
1171  }
1172 
1173  if (! raw_ngh_output_file.is_open()) return;
1174 
1175  // header
1176  raw_ngh_output_file << "// fields:\n//ele_id n_sides ns_side_neighbors[n] neighbors[n*ns] n_vb_neighbors vb_neighbors[n_vb]\n";
1177  raw_ngh_output_file << fmt::format("{}\n" , n_elements());
1178 
1179  int cit = 0;
1180 
1181  // map from higher dim elements to its lower dim neighbors, using gmsh IDs: ele->id()
1182  unsigned int undefined_ele_id = -1;
1184  for (auto ele : this->elements_range()) {
1185  if(ele->n_neighs_vb() > 0){
1186  for (unsigned int i = 0; i < ele->n_neighs_vb(); i++){
1187  ElementAccessor<3> higher_ele = ele->neigh_vb[i]->side()->element();
1188 
1189  auto search = neigh_vb_map.find(higher_ele.idx());
1190  if(search != neigh_vb_map.end()){
1191  // if found, add id to correct local side idx
1192  search->second[ele->neigh_vb[i]->side()->side_idx()] = ele.idx();
1193  }
1194  else{
1195  // if not found, create new vector, each side can have one vb neighbour
1196  std::vector<unsigned int> higher_ele_side_ngh(higher_ele->n_sides(), undefined_ele_id);
1197  higher_ele_side_ngh[ele->neigh_vb[i]->side()->side_idx()] = ele.idx();
1198  neigh_vb_map[higher_ele.idx()] = higher_ele_side_ngh;
1199  }
1200  }
1201  }
1202  }
1203 
1204  for (auto ele : this->elements_range()) {
1205  raw_ngh_output_file << ele.idx() << " ";
1206  raw_ngh_output_file << ele->n_sides() << " ";
1207 
1208  auto search_neigh = neigh_vb_map.end();
1209  for (unsigned int i = 0; i < ele->n_sides(); i++) {
1210  unsigned int n_side_neighs = ele.side(i)->edge().n_sides()-1; //n_sides - the current one
1211  // check vb neighbors (lower dimension)
1212  if(n_side_neighs == 0){
1213  //update search
1214  if(search_neigh == neigh_vb_map.end())
1215  search_neigh = neigh_vb_map.find(ele.idx());
1216 
1217  if(search_neigh != neigh_vb_map.end())
1218  if(search_neigh->second[i] != undefined_ele_id)
1219  n_side_neighs = 1;
1220  }
1221  raw_ngh_output_file << n_side_neighs << " ";
1222  }
1223 
1224  for (unsigned int i = 0; i < ele->n_sides(); i++) {
1225  Edge edge = ele.side(i)->edge();
1226  if(edge.n_sides() > 1){
1227  for (uint j = 0; j < edge.n_sides(); j++) {
1228  if(edge.side(j) != ele.side(i))
1229  raw_ngh_output_file << edge.side(j)->element().idx() << " ";
1230  }
1231  }
1232  //check vb neighbour
1233  else if(search_neigh != neigh_vb_map.end()
1234  && search_neigh->second[i] != undefined_ele_id){
1235  raw_ngh_output_file << search_neigh->second[i] << " ";
1236  }
1237  }
1238 
1239  // list higher dim neighbours
1240  raw_ngh_output_file << ele->n_neighs_vb() << " ";
1241  for (unsigned int i = 0; i < ele->n_neighs_vb(); i++)
1242  raw_ngh_output_file << ele->neigh_vb[i]->side()->element().idx() << " ";
1243 
1244  raw_ngh_output_file << endl;
1245  cit ++;
1246  }
1247  raw_ngh_output_file << "$EndFlowField\n" << endl;
1248 }
1249 
1250 
1252  // Copy boundary elements in temporary storage to the second part of the element vector
1253  for(ElementTmpData &e_data : bc_element_tmp_) {
1254  Element *ele = add_element_to_vector(e_data.elm_id);
1255  this->init_element(ele, e_data.elm_id, e_data.dim, e_data.region_idx,
1256  e_data.partition_id, e_data.node_ids);
1257  }
1258  // release memory
1259  unsigned int bdr_size = bc_element_tmp_.size();
1261  return bdr_size;
1262 }
1263 
1264 
1266  if (bc_mesh_ == nullptr) bc_mesh_ = new BCMesh(this);
1267  return bc_mesh_;
1268 }
1269 
1270 
1272  ASSERT_PTR(el_4_loc).error("Array 'el_4_loc' is not initialized. Did you call Partitioning::id_maps?\n");
1273 
1274  unsigned int i_proc, i_node, i_ghost_node, elm_node;
1275  unsigned int my_proc = el_ds->myp();
1276  unsigned int n_proc = el_ds->np();
1277 
1278  // distribute nodes between processes, every node is assigned to minimal process of elements that own node
1279  // fill node_proc vector with same values on all processes
1280  std::vector<unsigned int> node_proc( this->n_nodes(), n_proc );
1281  std::vector<bool> local_node_flag( this->n_nodes(), false );
1282 
1283  for ( auto elm : this->elements_range() ) {
1284  i_proc = elm.proc();
1285  for (elm_node=0; elm_node<elm->n_nodes(); elm_node++) {
1286  i_node = elm->node_idx(elm_node);
1287  if (i_proc == my_proc) local_node_flag[i_node] = true;
1288  if (i_proc < node_proc[i_node]) node_proc[i_node] = i_proc;
1289  }
1290  }
1291 
1292  unsigned int n_own_nodes=0, n_local_nodes=0; // number of own and ghost nodes
1293  for(uint loc_flag : local_node_flag) if (loc_flag) n_local_nodes++;
1294  for(uint i_proc : node_proc) {
1295  if (i_proc == my_proc)
1296  n_own_nodes++;
1297  else if (i_proc == n_proc)
1298  ASSERT(0)(find_node_id(n_own_nodes)).error("A node does not belong to any element!");
1299  }
1300 
1301  //DebugOut() << print_var(n_own_nodes) << print_var(n_local_nodes) << this->n_nodes();
1302  // create and fill node_4_loc_ (mapping local to global indexes)
1303  node_4_loc_ = new LongIdx [ n_local_nodes ];
1304  i_node=0;
1305  i_ghost_node = n_own_nodes;
1306  for (unsigned int i=0; i<this->n_nodes(); ++i) {
1307  if (local_node_flag[i]) {
1308  if (node_proc[i]==my_proc)
1309  node_4_loc_[i_node++] = i;
1310  else
1311  node_4_loc_[i_ghost_node++] = i;
1312  }
1313  }
1314 
1315  // Construct node distribution object, set number of local nodes (own+ghost)
1316  node_ds_ = new Distribution(n_own_nodes, PETSC_COMM_WORLD);
1317  node_ds_->get_lsizes_array(); // need to initialize lsizes data member
1319 
1320 }
1321 
1322 //-----------------------------------------------------------------------------
1323 // vim: set cindent:
Side::edge
Edge edge() const
Returns pointer to the edge connected to the side.
Definition: accessors_impl.hh:221
MeshOptimizer::calculate_node_curve_values_as_hilbert
void calculate_node_curve_values_as_hilbert()
Definition: mesh_optimizer.hh:59
Mesh::create_boundary_elements
unsigned int create_boundary_elements()
Create boundary elements from data of temporary structure, this method MUST be call after read mesh f...
Definition: mesh.cc:1251
Input::Type::Bool
Class for declaration of the input of type Bool.
Definition: type_base.hh:452
Distribution::np
unsigned int np() const
get num of processors
Definition: distribution.hh:105
reader_to_storage.hh
Boundary
Definition: accessors.hh:363
BidirectionalMap< int >
Mesh::BCMesh
friend class BCMesh
Definition: mesh.h:636
ref_element.hh
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.
Mesh::nodes_
Armor::Array< double > nodes_
Definition: mesh.h:615
Neighbour::mesh_
Mesh * mesh_
Pointer to Mesh to which belonged.
Definition: neighbours.h:136
MeshOptimizer::sort_elements
std::vector< int > sort_elements(std::vector< unsigned int > &elem_permutation)
Definition: mesh_optimizer.hh:108
Region::label
std::string label() const
Returns label of the region (using RegionDB)
Definition: region.cc:32
Mesh::get_input_intersection_variant
static const Input::Type::Selection & get_input_intersection_variant()
The definition of input record for selection of variant of file format.
Definition: mesh.cc:65
Mesh::get_element_boxes
std::vector< BoundingBox > get_element_boxes()
Compute bounding boxes of elements contained in mesh.
Definition: mesh.cc:1004
Neighbour::edge_idx_
unsigned int edge_idx_
Index of Edge in Mesh.
Definition: neighbours.h:138
Input::ReaderToStorage
Reader for (slightly) modified input files.
Definition: reader_to_storage.hh:96
Mesh::node_4_loc_
LongIdx * node_4_loc_
Index set assigning to local node index its global index.
Definition: mesh.h:654
bih_tree.hh
Input::ReaderToStorage::read_stream
void read_stream(istream &in, const Type::TypeBase &root_type, FileFormat format)
This method actually reads the given stream in.
Definition: reader_to_storage.cc:110
mesh_optimizer.hh
Mesh::in_record_
Input::Record in_record_
Definition: mesh.h:586
RefElement
Definition: ref_element.hh:339
Mesh::init_element_vector
void init_element_vector(unsigned int size)
Initialize element_vec_, set size and reset counters of boundary and bulk elements.
Definition: mesh.cc:1099
neighbours.h
Neighbour::elem_idx_
unsigned int elem_idx_
Index of element in Mesh::element_vec_.
Definition: neighbours.h:137
Mesh::n_exsides
int n_exsides
Definition: mesh.h:325
Mesh::init_element
void init_element(Element *ele, unsigned int elm_id, unsigned int dim, RegionIdx region_idx, unsigned int partition_id, std::vector< unsigned int > node_ids)
Initialize element.
Definition: mesh.cc:1068
Distribution::myp
unsigned int myp() const
get my processor
Definition: distribution.hh:107
RegionSetBase::get_input_type
static Input::Type::Abstract & get_input_type()
Definition: region_set.cc:25
MeshOptimizer::calculate_element_curve_values_as_hilbert_of_centers
void calculate_element_curve_values_as_hilbert_of_centers()
Definition: mesh_optimizer.hh:86
ASSERT
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:347
RegionDB::check_regions
void check_regions()
Definition: region.cc:461
Element::dim
unsigned int dim() const
Definition: elements.h:127
Mesh::elements_range
virtual Range< ElementAccessor< 3 > > elements_range() const
Returns range of bulk elements.
Definition: mesh.cc:1130
Mesh::n_sides
unsigned int n_sides() const
Definition: mesh.cc:238
Mesh::check_compatible_mesh
virtual std::shared_ptr< EquivalentMeshMap > check_compatible_mesh(Mesh &input_mesh)
Definition: mesh.cc:870
distribution.hh
Support classes for parallel programing.
bc_mesh.hh
Edge::n_sides
unsigned int n_sides() const
Returns number of sides aligned with the edge.
Definition: accessors.hh:343
BoundaryData::edge_idx_
unsigned int edge_idx_
Definition: mesh_data.hh:49
MixedMeshIntersections
Main class for computation of intersection of meshes of combined dimensions.
Definition: mixed_mesh_intersections.hh:64
compare_points
bool compare_points(const arma::vec3 &p1, const arma::vec3 &p2)
Definition: mesh.cc:862
Input::Record::val
const Ret val(const string &key) const
Definition: accessors_impl.hh:31
Mesh::boundary_loaded_size_
unsigned int boundary_loaded_size_
Count of boundary elements loaded from mesh file.
Definition: mesh.h:607
Armor::Array::vec
ArmaVec< Type, nr > vec(uint mat_index) const
Definition: armor.hh:821
Input::Type::Selection::close
const Selection & close() const
Close the Selection, no more values can be added.
Definition: type_selection.cc:65
Mesh::tree
DuplicateNodes * tree
Definition: mesh.h:317
Mesh::element_vec_
vector< Element > element_vec_
Definition: mesh.h:598
Neighbour::side
SideIter side()
Definition: neighbours.h:145
Element::pid_
int pid_
Id # of mesh partition.
Definition: elements.h:104
FilePath
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
ASSERT_DBG
#define ASSERT_DBG(expr)
Definition: include_fadbad.hh:28
MeshOptimizer::calculate_sizes
void calculate_sizes()
Definition: mesh_optimizer.hh:40
RegionIdx::is_valid
bool is_valid() const
Returns false if the region has undefined/invalid value.
Definition: region.hh:78
Input::Type::Double
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
Mesh::sort_permuted_nodes_elements
void sort_permuted_nodes_elements(std::vector< int > new_node_ids, std::vector< int > new_elem_ids)
Sort elements and nodes by order stored in permutation vectors.
Definition: mesh.cc:445
RegionDB::mark_used_region
void mark_used_region(unsigned int idx)
Definition: region.cc:235
Mesh::get_input_type
static const Input::Type::Record & get_input_type()
Definition: mesh.cc:76
BIHTree
Class for O(log N) lookup for intersections with a set of bounding boxes.
Definition: bih_tree.hh:38
THROW
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
Mesh::RegionSetBase
friend class RegionSetBase
Definition: mesh.h:632
Mesh::max_edge_sides_
unsigned int max_edge_sides_[3]
Maximal number of sides per one edge in the actual mesh (set in make_neighbours_and_edges()).
Definition: mesh.h:535
std::vector< uint >
Mesh::count_side_types
void count_side_types()
Definition: mesh.cc:472
ElementAccessor< 3 >
Partitioning
Class for the mesh partitioning. This should provide:
Definition: partitioning.hh:52
Input::Type::FileName::output
static FileName output()
The factory function for declaring type FileName for input files.
Definition: type_base.cc:531
system.hh
arma::vec3
Definition: doxy_dummy_defs.hh:17
RegionDB::get_region
Region get_region(unsigned int id, unsigned int dim)
Definition: region.cc:150
BoundaryData
Definition: mesh_data.hh:40
Mesh::node_range
Range< NodeAccessor< 3 > > node_range() const
Returns range of nodes.
Definition: mesh.cc:1136
Partitioning::get_loc_part
const LongIdx * get_loc_part() const
Definition: partitioning.cc:85
Mesh::boundary
Boundary boundary(uint edge_idx) const
Definition: mesh.cc:268
EdgeData
Definition: mesh_data.hh:25
Armor::Array::resize
void resize(uint size)
Definition: armor.hh:710
uint
unsigned int uint
Definition: mh_dofhandler.hh:101
RegionIdx::idx
unsigned int idx() const
Returns a global index of the region.
Definition: region.hh:82
Mesh::n_local_nodes_
unsigned int n_local_nodes_
Hold number of local nodes (own + ghost), value is equal with size of node_4_loc array.
Definition: mesh.h:658
fmt::format
std::string format(CStringRef format_str, ArgList args)
Definition: format.h:3141
duplicate_nodes.h
Neighbour
Definition: neighbours.h:117
ElementAccessor::element
const Element * element() const
Definition: accessors.hh:195
FilePath::open_stream
void open_stream(Stream &stream) const
Definition: file_path.cc:211
Input::Array::begin
Iterator< ValueType > begin() const
Definition: accessors_impl.hh:145
Mesh::read_regions_from_input
void read_regions_from_input(Input::Array region_list)
Definition: mesh.cc:979
RegionDB::create_label_from_id
std::string create_label_from_id(unsigned int id) const
Definition: region.cc:337
index_types.hh
RegionDB::el_to_reg_map_
MapElementIDToRegionID el_to_reg_map_
Definition: region.hh:571
BoundaryData::bc_ele_idx_
unsigned int bc_ele_idx_
Definition: mesh_data.hh:53
Mesh::get_intersection_search
IntersectionSearch get_intersection_search()
Getter for input type selection for intersection search algorithm.
Definition: mesh.cc:156
Mesh::optimize_memory_locality
bool optimize_memory_locality
Definition: mesh.h:565
Mesh::n_vb_neighbours
unsigned int n_vb_neighbours() const
Definition: mesh.cc:247
exceptions.hh
Mesh::n_elements
virtual unsigned int n_elements() const
Returns count of boundary or bulk elements.
Definition: mesh.h:398
Element
Definition: elements.h:40
MPI_Comm_rank
#define MPI_Comm_rank
Definition: mpi.h:236
ASSERT_PTR_DBG
#define ASSERT_PTR_DBG(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:340
Mesh::optimize
void optimize()
Definition: mesh.cc:435
Mesh::el_ds
Distribution * el_ds
Parallel distribution of elements.
Definition: mesh.h:652
Region::id
unsigned int id() const
Returns id of the region (using RegionDB)
Definition: region.cc:37
Mesh::intersect_element_lists
void intersect_element_lists(vector< unsigned int > const &nodes_list, vector< unsigned int > &intersection_element_list)
Definition: mesh.cc:499
Mesh::n_sides_
int n_sides_
Definition: mesh.h:326
EdgeData::side_
SideIter * side_
Definition: mesh_data.hh:32
Input::Iterator
Definition: accessors.hh:143
Input::Type::Default
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
RegionDB::undefined_dim
static const unsigned int undefined_dim
Definition: region.hh:332
MeshOptimizer::sort_nodes
std::vector< int > sort_nodes(std::vector< unsigned int > &node_permutation)
Definition: mesh_optimizer.hh:104
Mesh::output_internal_ngh_data
void output_internal_ngh_data()
Output of neighboring data into raw output.
Definition: mesh.cc:1156
Mesh::get_local_part
virtual const LongIdx * get_local_part()
Definition: mesh.cc:278
Distribution
Definition: distribution.hh:50
Mesh::node_ids_
BidirectionalMap< int > node_ids_
Maps node ids to indexes into vector node_vec_.
Definition: mesh.h:618
ElementAccessor::bounding_box
BoundingBox bounding_box() const
Definition: accessors.hh:253
accessors.hh
Mesh::BBsearch
@ BBsearch
Definition: mesh.h:138
Mesh::check_mesh_on_read
void check_mesh_on_read()
Definition: mesh.cc:310
Mesh::n_corners
unsigned int n_corners()
Definition: mesh.cc:252
BIHTree::find_point
void find_point(const Space< 3 >::Point &point, std::vector< unsigned int > &result_list, bool full_list=false) const
Definition: bih_tree.cc:287
Mesh::find_lower_dim_element
bool find_lower_dim_element(vector< unsigned int > &element_list, unsigned int dim, unsigned int &element_idx)
Definition: mesh.cc:531
Input::Record
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
Mesh::get_bih_tree
const BIHTree & get_bih_tree()
Getter for BIH. Creates and compute BIH at first call.
Definition: mesh.cc:1019
Mesh::distribute_nodes
void distribute_nodes()
Fill array node_4_loc_ and create object node_ds_ according to element distribution.
Definition: mesh.cc:1271
sys_profiler.hh
Edge::side
SideIter side(const unsigned int i) const
Gets side iterator of the i -th side.
Definition: accessors_impl.hh:171
undef_idx
const unsigned int undef_idx
Definition: index_types.hh:32
Mesh::get_bc_mesh
BCMesh * get_bc_mesh()
Create boundary mesh if doesn't exist and return it.
Definition: mesh.cc:1265
RegionIdx::is_boundary
bool is_boundary() const
Returns true if it is a Boundary region and false if it is a Bulk region.
Definition: region.hh:74
mixed_mesh_intersections.hh
Input::Type::Record::allow_auto_conversion
virtual Record & allow_auto_conversion(const string &from_key)
Allows shorter input of the Record providing only value of the from_key given as the parameter.
Definition: type_record.cc:133
accessors.hh
Mesh::el_4_loc
LongIdx * el_4_loc
Index set assigning to local element index its global index.
Definition: mesh.h:650
Mesh::elements_id_maps
void elements_id_maps(vector< LongIdx > &bulk_elements_id, vector< LongIdx > &boundary_elements_id) const
Definition: mesh.cc:829
INPUT_CATCH
#define INPUT_CATCH(ExceptionType, AddressEITag, input_accessor)
Definition: accessors.hh:63
Element::n_neighs_vb_
unsigned int n_neighs_vb_
Definition: elements.h:106
Mesh::side_nodes
vector< vector< vector< unsigned int > > > side_nodes
Definition: mesh.h:339
Mesh::node_permutation_
std::vector< unsigned int > node_permutation_
Vector of node permutations of optimized mesh (see class MeshOptimizer)
Definition: mesh.h:624
Mesh::canonical_faces
void canonical_faces()
Definition: mesh.cc:390
Mesh::region_db_
RegionDB region_db_
Definition: mesh.h:571
Distribution::get_lsizes_array
const unsigned int * get_lsizes_array()
get local sizes array
Definition: distribution.cc:142
Input::AbstractRecord
Accessor to the polymorphic input data of a type given by an AbstracRecord object.
Definition: accessors.hh:458
Neighbour::element
ElementAccessor< 3 > element()
Definition: neighbours.h:161
Input::Type::Default::obligatory
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
ElementAccessor::jacobian_S3
double jacobian_S3() const
Definition: accessors.hh:131
Mesh::IntersectionSearch
IntersectionSearch
Types of search algorithm for finding intersection candidates.
Definition: mesh.h:135
RegionDB::close
void close()
Definition: region.cc:249
Mesh::find_elem_id
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:412
ASSERT_EQ
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual)
Definition: asserts.hh:328
Element::init
void init(unsigned int dim, RegionIdx reg)
Definition: elements.cc:53
Mesh::init
void init()
Definition: mesh.cc:162
RegionDB::find_id
Region find_id(unsigned int id, unsigned int dim) const
Definition: region.cc:180
Input::Record::opt_val
bool opt_val(const string &key, Ret &value) const
Definition: accessors_impl.hh:107
Input::format_JSON
@ format_JSON
Definition: reader_to_storage.hh:60
Input::Type::Record::declare_key
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:503
Input::ReaderToStorage::get_root_interface
T get_root_interface() const
Returns the root accessor.
Definition: reader_to_storage.cc:150
Side::element
ElementAccessor< 3 > element() const
Returns iterator to the element of the side.
Definition: accessors_impl.hh:212
mesh.h
Mesh::bih_tree_
std::shared_ptr< BIHTree > bih_tree_
Definition: mesh.h:580
Mesh::node_ds_
Distribution * node_ds_
Parallel distribution of nodes. Depends on elements distribution.
Definition: mesh.h:656
Input::Type::Selection
Template for classes storing finite set of named values.
Definition: type_selection.hh:65
std::map< unsigned int, unsigned int >
Side::n_nodes
unsigned int n_nodes() const
Returns number of nodes of the side.
Definition: accessors.hh:450
Element::region_idx_
RegionIdx region_idx_
Definition: elements.h:113
ElementAccessor::node
NodeAccessor< 3 > node(unsigned int ni) const
Definition: accessors.hh:245
Input::Type::FileName::input
static FileName input()
The factory function for declaring type FileName for input files.
Definition: type_base.cc:524
Mesh::same_sides
bool same_sides(const SideIter &si, vector< unsigned int > &side_nodes)
Definition: mesh.cc:553
Armor::Array::set
ArrayMatSet set(uint index)
Definition: armor.hh:838
Input::Type::Record::close
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Mesh::element_to_neigh_vb
void element_to_neigh_vb()
Definition: mesh.cc:767
Element::node_idx
unsigned int node_idx(unsigned int ni) const
Return index (in Mesh::node_vec) of ni-th node.
Definition: elements.h:73
Mesh::check_element_size
void check_element_size(unsigned int elem_idx) const
Check if given index is in element_vec_.
Definition: mesh.cc:1148
Mesh::add_element
void add_element(unsigned int elm_id, unsigned int dim, unsigned int region_id, unsigned int partition_id, std::vector< unsigned int > node_ids)
Add new element of given id to mesh.
Definition: mesh.cc:1045
Side::side_idx
unsigned int side_idx() const
Returns local index of the side on the element.
Definition: accessors.hh:458
Input::Type
Definition: balance.hh:41
partitioning.hh
Partitioning::get_input_type
static const Input::Type::Record & get_input_type()
Definition: partitioning.cc:49
RegionDB::implicit_boundary_region
Region implicit_boundary_region()
Definition: region.cc:75
Input::Type::Record
Record type proxy class.
Definition: type_record.hh:182
MPI_Comm
int MPI_Comm
Definition: mpi.h:141
Mesh::Element
friend class Element
Definition: mesh.h:633
Mesh::~Mesh
virtual ~Mesh()
Destructor.
Definition: mesh.cc:213
Mesh::Edge
friend class Edge
Definition: mesh.h:630
Element::neigh_vb
Neighbour ** neigh_vb
Definition: elements.h:86
Element::nodes_
std::array< unsigned int, 4 > nodes_
indices to element's nodes
Definition: elements.h:117
Interaction
Definition: ref_element.hh:286
Mesh::create_node_element_lists
void create_node_element_lists()
Definition: mesh.cc:486
Element::edge_idx
unsigned int edge_idx(unsigned int edg_idx) const
Return edge_idx of given index.
Definition: elements.h:142
EdgeData::n_sides
unsigned int n_sides
Definition: mesh_data.hh:29
LongIdx
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Definition: index_types.hh:24
Mesh::setup_topology
void setup_topology()
Definition: mesh.cc:399
input_type.hh
Mesh::add_element_to_vector
Element * add_element_to_vector(int id)
Adds element to mesh data structures (element_vec_, element_ids_), returns pointer to this element.
Definition: mesh.cc:1122
Mesh::modify_element_ids
void modify_element_ids(const RegionDB::MapElementIDToRegionID &map)
Definition: mesh.cc:285
Mesh
Definition: mesh.h:98
Armor::Array::reinit
void reinit(uint size)
Definition: armor.hh:698
Element::n_sides
unsigned int n_sides() const
Definition: elements.h:138
Mesh::BIHsearch
@ BIHsearch
Definition: mesh.h:136
Mesh::elem_index
int elem_index(int elem_id) const
For element of given elem_id returns index in element_vec_ or (-1) if element doesn't exist.
Definition: mesh.h:406
Range
Range helper class.
Definition: range_wrapper.hh:65
Input::Type::Array
Class for declaration of inputs sequences.
Definition: type_base.hh:339
OLD_ASSERT
#define OLD_ASSERT(...)
Definition: global_defs.h:108
Mesh::edges
std::vector< EdgeData > edges
Vector of MH edges, this should not be part of the geometrical mesh.
Definition: mesh.h:621
Mesh::make_neighbours_and_edges
void make_neighbours_and_edges()
Definition: mesh.cc:570
BidirectionalMap::clear
void clear()
Clear the content. Do not release memory.
Definition: bidirectional_map.hh:121
Mesh::init_node_vector
void init_node_vector(unsigned int size)
Initialize node_vec_, set size.
Definition: mesh.cc:1113
Mesh::elem_permutation_
std::vector< unsigned int > elem_permutation_
Vector of element permutations of optimized mesh (see class MeshOptimizer)
Definition: mesh.h:627
Edge
Definition: accessors.hh:308
Input::Array
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
Mesh::Mesh
Mesh()
Definition: mesh.cc:102
WarningOut
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
RegionIdx
Definition: region.hh:67
Input::Type::Record::finish
FinishStatus finish(FinishStatus finish_type=FinishStatus::regular_) override
Finish declaration of the Record type.
Definition: type_record.cc:243
Armor::Array::append
void append(const ArmaMat< Type, nr, nc > &item)
Definition: armor.hh:736
region_set.hh
ElementAccessor::region
Region region() const
Definition: accessors.hh:201
Mesh::node_index
int node_index(int node_id) const
For node of given node_id returns index in element_vec_ or (-1) if node doesn't exist.
Definition: mesh.h:424
BCMesh
Class represents boundary part of mesh.
Definition: bc_mesh.hh:35
Mesh::check_and_finish
void check_and_finish()
Definition: mesh.cc:989
Mesh::edge_range
Range< Edge > edge_range() const
Returns range of edges.
Definition: mesh.cc:1142
NodeAccessor
Definition: mesh.h:55
Mesh::node_elements_
vector< vector< unsigned int > > node_elements_
For each node the vector contains a list of elements that use this node.
Definition: mesh.h:403
MPI_COMM_WORLD
#define MPI_COMM_WORLD
Definition: mpi.h:123
ASSERT_LT_DBG
#define ASSERT_LT_DBG(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:300
region.hh
Mesh::element_ids_
BidirectionalMap< int > element_ids_
Maps element ids to indexes into vector element_vec_.
Definition: mesh.h:610
Armor::Array< double >
Mesh::add_node
void add_node(unsigned int node_id, arma::vec3 coords)
Add new node of given id and coordinates to mesh.
Definition: mesh.cc:1037
Mesh::bc_element_tmp_
vector< ElementTmpData > bc_element_tmp_
Hold data of boundary elements during reading mesh (allow to preserve correct order during reading of...
Definition: mesh.h:601
ElementAccessor::idx
unsigned int idx() const
Return local idx of element in boundary / bulk part of element vector.
Definition: accessors.hh:223
NDEF
#define NDEF
Definition: mesh.cc:61
RegionDB::add_region
Region add_region(unsigned int id, const std::string &label, unsigned int dim, const std::string &address="implicit")
Definition: region.cc:85
MeshOptimizer
Definition: mesh_optimizer.hh:35
Mesh::part_
std::shared_ptr< Partitioning > part_
Definition: mesh.h:575
Mesh::check_compatible_elements
unsigned int check_compatible_elements(Mesh *source_mesh, Mesh *target_mesh, const std::vector< unsigned int > &node_ids, std::vector< LongIdx > &map)
Definition: mesh.cc:935
print_var
#define print_var(var)
Definition: logger.hh:292
Mesh::bulk_size_
unsigned int bulk_size_
Count of bulk elements.
Definition: mesh.h:604
BidirectionalMap::add_item
unsigned int add_item(T val)
Add new item at the end position of map.
Definition: bidirectional_map.hh:106
RegionDB::print_region_table
void print_region_table(std::ostream &stream) const
Definition: region.cc:409
Mesh::row_4_el
LongIdx * row_4_el
Index set assigning to global element index the local index used in parallel vectors.
Definition: mesh.h:648
Mesh::node_elements
const vector< vector< unsigned int > > & node_elements()
Definition: mesh.cc:1091
Mesh::global_snap_radius
double global_snap_radius() const
Maximal distance of observe point from Mesh relative to its size.
Definition: mesh.cc:1028
Mesh::boundary_
vector< BoundaryData > boundary_
Definition: mesh.h:301
Mesh::Boundary
friend class Boundary
Definition: mesh.h:635
Mesh::n_local_nodes
unsigned int n_local_nodes() const
Definition: mesh.h:213
Mesh::get_part
virtual Partitioning * get_part()
Definition: mesh.cc:274
Mesh::n_insides
int n_insides
Definition: mesh.h:324
Side::node
NodeAccessor< 3 > node(unsigned int i) const
Returns node for given local index i on the side.
Definition: accessors_impl.hh:206
Mesh::bc_mesh_
BCMesh * bc_mesh_
Boundary mesh, object is created only if it's necessary.
Definition: mesh.h:660
Mesh::mixed_intersections
MixedMeshIntersections & mixed_intersections()
Definition: mesh.cc:800
DebugOut
#define DebugOut()
Macro defining 'debug' record of log.
Definition: logger.hh:284
Input::Type::Selection::add_value
Selection & add_value(const int value, const std::string &key, const std::string &description="", TypeBase::attribute_map attributes=TypeBase::attribute_map())
Adds one new value with name given by key to the Selection.
Definition: type_selection.cc:50
Mesh::BIHonly
@ BIHonly
Definition: mesh.h:137
ASSERT_PTR
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR)
Definition: asserts.hh:336
ElementAccessor::side
SideIter side(const unsigned int loc_index)
Definition: accessors_impl.hh:139
Mesh::node
NodeAccessor< 3 > node(unsigned int idx) const
Create and return NodeAccessor to node of given idx.
Definition: mesh.cc:823
Mesh::element_accessor
virtual ElementAccessor< 3 > element_accessor(unsigned int idx) const
Create and return ElementAccessor to element of given idx.
Definition: mesh.cc:817
BoundaryData::mesh_
Mesh * mesh_
Definition: mesh_data.hh:54
BidirectionalMap::reserve
void reserve(unsigned int init_size=0)
Reset data of map, reserve space for given size.
Definition: bidirectional_map.hh:139
Input::Type::Default::optional
static Default optional()
The factory function to make an empty default value which is optional.
Definition: type_record.hh:124
Mesh::n_nodes
virtual unsigned int n_nodes() const
Definition: mesh.h:170
START_TIMER
#define START_TIMER(tag)
Starts a timer with specified tag.
Definition: sys_profiler.hh:115
Input::Array::end
IteratorBase end() const
Definition: accessors_impl.hh:157
Mesh::find_node_id
int find_node_id(unsigned int pos) const
Return node id (in GMSH file) of node of given position in node vector.
Definition: mesh.h:430
SideIter
Definition: accessors.hh:504
Input::Record::is_empty
bool is_empty() const
Definition: accessors.hh:365
Element::boundary_idx_
unsigned int * boundary_idx_
Definition: elements.h:82
Mesh::edge
Edge edge(uint edge_idx) const
Definition: mesh.cc:262
Mesh::add_physical_name
void add_physical_name(unsigned int dim, unsigned int id, std::string name)
Add new node of given id and coordinates to mesh.
Definition: mesh.cc:1032
Mesh::vb_neighbours_
vector< Neighbour > vb_neighbours_
Definition: mesh.h:322
END_TIMER
#define END_TIMER(tag)
Ends a timer with specified tag.
Definition: sys_profiler.hh:149
Element::n_nodes
unsigned int n_nodes() const
Definition: elements.h:132
Mesh::ElementTmpData
Definition: mesh.h:463
DuplicateNodes
Definition: duplicate_nodes.h:96
BidirectionalMap::size
unsigned int size() const
Return size of map.
Definition: bidirectional_map.hh:80
range_wrapper.hh
Implementation of range helper class.
MessageOut
#define MessageOut()
Macro defining 'message' record of log.
Definition: logger.hh:275
node_accessor.hh
Mesh::intersections
std::shared_ptr< MixedMeshIntersections > intersections
Definition: mesh.h:309