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