Flow123d  release_2.2.0-914-gf1a3a4f
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 #include <unistd.h>
20 #include <set>
21 
22 
23 #include "system/system.hh"
24 #include "system/exceptions.hh"
26 #include "input/input_type.hh"
27 #include "input/accessors.hh"
28 #include "system/sys_profiler.hh"
29 #include "la/distribution.hh"
30 
31 #include "mesh/mesh.h"
32 #include "mesh/ref_element.hh"
33 
34 // think about following dependencies
35 #include "mesh/boundaries.h"
36 #include "mesh/accessors.hh"
37 #include "mesh/partitioning.hh"
38 
39 
40 #include "mesh/bih_tree.hh"
41 
45 
47 
48 
49 //TODO: sources, concentrations, initial condition and similarly boundary conditions should be
50 // instances of a Element valued field
51 // concentrations is in fact reimplemented in transport REMOVE it HERE
52 
53 // After removing non-geometrical things from mesh, this should be part of mash initializing.
54 #include "mesh/region.hh"
55 
56 #define NDEF -1
57 
58 namespace IT = Input::Type;
59 
61  return Input::Type::Selection("Types of search algorithm for finding intersection candidates.")
62  .add_value(Mesh::BIHsearch, "BIHsearch",
63  "Use BIH for finding initial candidates, then continue by prolongation.")
64  .add_value(Mesh::BIHonly, "BIHonly",
65  "Use BIH for finding all candidates.")
66  .add_value(Mesh::BBsearch, "BBsearch",
67  "Use bounding boxes for finding initial candidates, then continue by prolongation.")
68  .close();
69 }
70 
72  return IT::Record("Mesh","Record with mesh related data." )
73  .allow_auto_conversion("mesh_file")
75  "Input file with mesh description.")
77  "List of additional region and region set definitions not contained in the mesh. "
78  "There are three region sets implicitly defined:\n\n"
79  "- ALL (all regions of the mesh)\n"
80  "- .BOUNDARY (all boundary regions)\n"
81  "- BULK (all bulk regions)")
82  .declare_key("partitioning", Partitioning::get_input_type(), IT::Default("\"any_neighboring\""), "Parameters of mesh partitioning algorithms.\n" )
83  .declare_key("print_regions", IT::Bool(), IT::Default("true"), "If true, print table of all used regions.")
84  .declare_key("intersection_search", Mesh::get_input_intersection_variant(),
85  IT::Default("\"BIHsearch\""), "Search algorithm for element intersections.")
86  .declare_key("global_observe_search_radius", IT::Double(0.0), IT::Default("1E-3"),
87  "Maximal distance of observe point from Mesh relative to its size (bounding box). "
88  "Value is global and it can be rewrite at arbitrary ObservePoint by setting the key search_radius.")
89  .close();
90 }
91 
92 const unsigned int Mesh::undef_idx;
93 
95 : row_4_el(nullptr),
96  el_4_loc(nullptr),
97  el_ds(nullptr)
98 {}
99 
100 
101 
103 : in_record_(in_record),
104  comm_(com),
105  row_4_el(nullptr),
106  el_4_loc(nullptr),
107  el_ds(nullptr)
108 {
109  // set in_record_, if input accessor is empty
110  if (in_record_.is_empty()) {
111  istringstream is("{mesh_file=\"\"}");
112  Input::ReaderToStorage reader;
113  IT::Record &in_rec = const_cast<IT::Record &>(Mesh::get_input_type());
114  in_rec.finish();
115  reader.read_stream(is, in_rec, Input::FileFormat::format_JSON);
117  }
118 
120 }
121 
123 {
124  return in_record_.val<Mesh::IntersectionSearch>("intersection_search");
125 }
126 
127 
129 {
130 
131  n_insides = NDEF;
132  n_exsides = NDEF;
133  n_sides_ = NDEF;
134 
135  // number of element of particular dimension
136  n_lines = 0;
137  n_triangles = 0;
138  n_tetrahedras = 0;
139 
140  for (int d=0; d<3; d++) max_edge_sides_[d] = 0;
141 
142  // Initialize numbering of nodes on sides.
143  // This is temporary solution, until class Element is templated
144  // by dimension. Then we can replace Mesh::side_nodes by
145  // RefElement<dim>::side_nodes.
146 
147  // indices of side nodes in element node array
148  // Currently this is made ad libitum
149  // with some ordering here we can get sides with correct orientation.
150  // This speedup normal calculation.
151 
152  side_nodes.resize(3); // three side dimensions
153  for(int i=0; i < 3; i++) {
154  side_nodes[i].resize(i+2); // number of sides
155  for(int j=0; j < i+2; j++)
156  side_nodes[i][j].resize(i+1);
157  }
158 
159  for (unsigned int sid=0; sid<RefElement<1>::n_sides; sid++)
160  for (unsigned int nid=0; nid<RefElement<1>::n_nodes_per_side; nid++)
161  side_nodes[0][sid][nid] = RefElement<1>::interact(Interaction<0,0>(sid))[nid];
162 
163  for (unsigned int sid=0; sid<RefElement<2>::n_sides; sid++)
164  for (unsigned int nid=0; nid<RefElement<2>::n_nodes_per_side; nid++)
165  side_nodes[1][sid][nid] = RefElement<2>::interact(Interaction<0,1>(sid))[nid];
166 
167  for (unsigned int sid=0; sid<RefElement<3>::n_sides; sid++)
168  for (unsigned int nid=0; nid<RefElement<3>::n_nodes_per_side; nid++)
169  side_nodes[2][sid][nid] = RefElement<3>::interact(Interaction<0,2>(sid))[nid];
170 }
171 
172 
174  for(Edge &edg : this->edges)
175  if (edg.side_) delete[] edg.side_;
176 
177  FOR_ELEMENTS( this, ele ) {
178  if (ele->node) delete[] ele->node;
179  if (ele->edge_idx_) delete[] ele->edge_idx_;
180  if (ele->permutation_idx_) delete[] ele->permutation_idx_;
181  if (ele->boundary_idx_) delete[] ele->boundary_idx_;
182  if (ele->neigh_vb) delete[] ele->neigh_vb;
183  }
184 
185  for(unsigned int idx=0; idx < this->bc_elements.size(); idx++) {
186  Element *ele=&(bc_elements[idx]);
187  if (ele->node) delete[] ele->node;
188  if (ele->edge_idx_) delete[] ele->edge_idx_;
189  if (ele->permutation_idx_) delete[] ele->permutation_idx_;
190  if (ele->boundary_idx_) delete[] ele->boundary_idx_;
191  }
192 
193  if (row_4_el != nullptr) delete[] row_4_el;
194  if (el_4_loc != nullptr) delete[] el_4_loc;
195  if (el_ds != nullptr) delete el_ds;
196 }
197 
198 
199 unsigned int Mesh::n_sides()
200 {
201  if (n_sides_ == NDEF) {
202  n_sides_=0;
203  FOR_ELEMENTS(this, ele) n_sides_ += ele->n_sides();
204  }
205  return n_sides_;
206 }
207 
208 unsigned int Mesh::n_corners() {
209  unsigned int li, count = 0;
210  FOR_ELEMENTS(this, ele) {
211  FOR_ELEMENT_NODES(ele, li) {
212  count++;
213  }
214  }
215  return count;
216 }
217 
219  return part_.get();
220 }
221 
222 
223 //=============================================================================
224 // COUNT ELEMENT TYPES
225 //=============================================================================
226 
228  FOR_ELEMENTS(this, elm)
229  switch (elm->dim()) {
230  case 1:
231  n_lines++;
232  break;
233  case 2:
234  n_triangles++;
235  break;
236  case 3:
237  n_tetrahedras++;
238  break;
239  }
240 }
241 
242 
243 
245  for (auto elem_to_region : map) {
246  ElementIter ele = this->element.find_id(elem_to_region.first);
247  ele->region_idx_ = region_db_.get_region( elem_to_region.second, ele->dim() );
249  }
250 }
251 
252 
253 
255  START_TIMER("MESH - setup topology");
256 
258 
259  // check mesh quality
260  FOR_ELEMENTS(this, ele)
261  if (ele->quality_measure_smooth() < 0.001) WarningOut().fmt("Bad quality (<0.001) of the element {}.\n", ele.id());
262 
267 
268  part_ = std::make_shared<Partitioning>(this, in_record_.val<Input::Record>("partitioning") );
269 
270  // create parallel distribution and numbering of elements
271  IdxInt *id_4_old = new IdxInt[element.size()];
272  int i = 0;
273  FOR_ELEMENTS(this, ele)
274  id_4_old[i++] = ele.index();
275  part_->id_maps(element.size(), id_4_old, el_ds, el_4_loc, row_4_el);
276 
277  delete[] id_4_old;
278 }
279 
280 
281 //
283 {
284 
285  n_insides = 0;
286  n_exsides = 0;
287  FOR_SIDES(this, sde ) {
288  if (sde->is_external()) n_exsides++;
289  else n_insides++;
290  }
291 }
292 
293 
294 
296  // for each node we make a list of elements that use this node
297  node_elements_.resize(node_vector.size());
298 
299  FOR_ELEMENTS( this, e )
300  for (unsigned int n=0; n<e->n_nodes(); n++)
301  node_elements_[node_vector.index(e->node[n])].push_back(e->index());
302 
303  for (vector<vector<unsigned int> >::iterator n=node_elements_.begin(); n!=node_elements_.end(); n++)
304  stable_sort(n->begin(), n->end());
305 }
306 
307 
308 void Mesh::intersect_element_lists(vector<unsigned int> const &nodes_list, vector<unsigned int> &intersection_element_list)
309 {
310  if (nodes_list.size() == 0) {
311  intersection_element_list.clear();
312  } else if (nodes_list.size() == 1) {
313  intersection_element_list = node_elements_[ nodes_list[0] ];
314  } else {
315  vector<unsigned int>::const_iterator it1=nodes_list.begin();
317  intersection_element_list.resize( node_elements_[*it1].size() ); // make enough space
318 
319  it1=set_intersection(
320  node_elements_[*it1].begin(), node_elements_[*it1].end(),
321  node_elements_[*it2].begin(), node_elements_[*it2].end(),
322  intersection_element_list.begin());
323  intersection_element_list.resize(it1-intersection_element_list.begin()); // resize to true size
324 
325  for(;it2<nodes_list.end();++it2) {
326  it1=set_intersection(
327  intersection_element_list.begin(), intersection_element_list.end(),
328  node_elements_[*it2].begin(), node_elements_[*it2].end(),
329  intersection_element_list.begin());
330  intersection_element_list.resize(it1-intersection_element_list.begin()); // resize to true size
331  }
332  }
333 }
334 
335 
337  unsigned int dim, unsigned int &element_idx) {
338  bool is_neighbour = false;
339 
340  vector<unsigned int>::iterator e_dest=element_list.begin();
341  for( vector<unsigned int>::iterator ele = element_list.begin(); ele!=element_list.end(); ++ele)
342  if (elements[*ele].dim() == dim) { // keep only indexes of elements of same dimension
343  *e_dest=*ele;
344  ++e_dest;
345  } else if (elements[*ele].dim() == dim-1) { // get only first element of lower dimension
346  if (is_neighbour) xprintf(UsrErr, "Too matching elements id: %d and id: %d in the same mesh.\n",
347  elements(*ele).id(), elements(element_idx).id() );
348 
349  is_neighbour = true;
350  element_idx = *ele;
351  }
352  element_list.resize( e_dest - element_list.begin());
353  return is_neighbour;
354 }
355 
357  // 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)
358  unsigned int ni=0;
359  while ( ni < si->n_nodes()
360  && find(side_nodes.begin(), side_nodes.end(), node_vector.index( si->node(ni) ) ) != side_nodes.end() ) ni++;
361  return ( ni == si->n_nodes() );
362 }
363 
364 /**
365  * TODO:
366  * - use std::is_any for setting is_neigbour
367  * - possibly make appropriate constructors for Edge and Neighbour
368  * - check side!=-1 when searching neigbouring element
369  * - process bc_elements first, there should be no Neigh, but check it
370  * set Edge and boundary there
371  */
372 
374 {
375  Neighbour neighbour;
376  Edge *edg;
377  unsigned int ngh_element_idx, last_edge_idx;
378 
380 
381  // pointers to created edges
382  //vector<Edge *> tmp_edges;
383  edges.resize(0); // be sure that edges are empty
384 
386  vector<unsigned int> intersection_list; // list of elements in intersection of node element lists
387 
388  for( ElementFullIter bc_ele = bc_elements.begin(); bc_ele != bc_elements.end(); ++bc_ele) {
389  // Find all elements that share this side.
390  side_nodes.resize(bc_ele->n_nodes());
391  for (unsigned n=0; n<bc_ele->n_nodes(); n++) side_nodes[n] = node_vector.index(bc_ele->node[n]);
392  intersect_element_lists(side_nodes, intersection_list);
393  bool is_neighbour = find_lower_dim_element(element, intersection_list, bc_ele->dim() +1, ngh_element_idx);
394  if (is_neighbour) {
395  xprintf(UsrErr, "Boundary element (id: %d) match a regular element (id: %d) of lower dimension.\n",
396  bc_ele.id(), element(ngh_element_idx).id());
397  } else {
398  if (intersection_list.size() == 0) {
399  // no matching dim+1 element found
400  WarningOut().fmt("Lonely boundary element, id: {}, region: {}, dimension {}.\n",
401  bc_ele.id(), bc_ele->region().id(), bc_ele->dim());
402  continue; // skip the boundary element
403  }
404  last_edge_idx=edges.size();
405  edges.resize(last_edge_idx+1);
406  edg = &( edges.back() );
407  edg->n_sides = 0;
408  edg->side_ = new struct SideIter[ intersection_list.size() ];
409 
410  // common boundary object
411  unsigned int bdr_idx=boundary_.size();
412  boundary_.resize(bdr_idx+1);
413  Boundary &bdr=boundary_.back();
414  bdr.bc_ele_idx_ = bc_ele.index();
415  bdr.edge_idx_ = last_edge_idx;
416  bdr.mesh_=this;
417 
418  // for 1d boundaries there can be more then one 1d elements connected to the boundary element
419  // we do not detect this case later in the main search over bulk elements
420  for( vector<unsigned int>::iterator isect = intersection_list.begin(); isect!=intersection_list.end(); ++isect) {
421  Element *elem = &(element[*isect]);
422  for (unsigned int ecs=0; ecs<elem->n_sides(); ecs++) {
423  SideIter si = elem->side(ecs);
424  if ( same_sides( si, side_nodes) ) {
425  if (elem->edge_idx_[ecs] != Mesh::undef_idx) {
426  OLD_ASSERT(elem->boundary_idx_!=nullptr, "Null boundary idx array.\n");
427  int last_bc_ele_idx=this->boundary_[elem->boundary_idx_[ecs]].bc_ele_idx_;
428  int new_bc_ele_idx=bc_ele.index();
429  THROW( ExcDuplicateBoundary()
430  << EI_ElemLast(this->bc_elements.get_id(last_bc_ele_idx))
431  << EI_RegLast(this->bc_elements[last_bc_ele_idx].region().label())
432  << EI_ElemNew(this->bc_elements.get_id(new_bc_ele_idx))
433  << EI_RegNew(this->bc_elements[new_bc_ele_idx].region().label())
434  );
435  }
436  elem->edge_idx_[ecs] = last_edge_idx;
437  edg->side_[ edg->n_sides++ ] = si;
438 
439  if (elem->boundary_idx_ == NULL) {
440  elem->boundary_idx_ = new unsigned int [ elem->n_sides() ];
441  std::fill( elem->boundary_idx_, elem->boundary_idx_ + elem->n_sides(), Mesh::undef_idx);
442  }
443  elem->boundary_idx_[ecs] = bdr_idx;
444  break; // next element in intersection list
445  }
446  }
447  }
448 
449  }
450 
451  }
452  // Now we go through all element sides and create edges and neighbours
453  FOR_ELEMENTS( this, e )
454  {
455  for (unsigned int s=0; s<e->n_sides(); s++)
456  {
457  // skip sides that were already found
458  if (e->edge_idx_[s] != Mesh::undef_idx) continue;
459 
460 
461  // Find all elements that share this side.
462  side_nodes.resize(e->side(s)->n_nodes());
463  for (unsigned n=0; n<e->side(s)->n_nodes(); n++) side_nodes[n] = node_vector.index(e->side(s)->node(n));
464  intersect_element_lists(side_nodes, intersection_list);
465 
466  bool is_neighbour = find_lower_dim_element(element, intersection_list, e->dim(), ngh_element_idx);
467 
468  if (is_neighbour) { // edge connects elements of different dimensions
469  neighbour.element_ = &(element[ngh_element_idx]);
470  } else { // edge connects only elements of the same dimension
471  // Allocate the array of sides.
472  last_edge_idx=edges.size();
473  edges.resize(last_edge_idx+1);
474  edg = &( edges.back() );
475  edg->n_sides = 0;
476  edg->side_ = new struct SideIter[ intersection_list.size() ];
477  if (intersection_list.size() > max_edge_sides_[e->dim()-1])
478  max_edge_sides_[e->dim()-1] = intersection_list.size();
479 
480  if (intersection_list.size() == 1) { // outer edge, create boundary object as well
481  edg->n_sides=1;
482  edg->side_[0] = e->side(s);
483  e->edge_idx_[s] = last_edge_idx;
484 
485  if (e->boundary_idx_ == NULL) {
486  e->boundary_idx_ = new unsigned int [ e->n_sides() ];
487  std::fill( e->boundary_idx_, e->boundary_idx_ + e->n_sides(), Mesh::undef_idx);
488  }
489 
490  unsigned int bdr_idx=boundary_.size();
491  boundary_.resize(bdr_idx+1);
492  Boundary &bdr=boundary_.back();
493  e->boundary_idx_[s] = bdr_idx;
494 
495  // fill boundary element
496  ElementFullIter bc_ele = bc_elements.add_item( -bdr_idx ); // use negative bcd index as ID,
497  bc_ele->init(e->dim()-1, this, region_db_.implicit_boundary_region() );
498  region_db_.mark_used_region( bc_ele->region_idx_.idx() );
499  for(unsigned int ni = 0; ni< side_nodes.size(); ni++) bc_ele->node[ni] = &( node_vector[side_nodes[ni]] );
500 
501  // fill Boundary object
502  bdr.edge_idx_ = last_edge_idx;
503  bdr.bc_ele_idx_ = bc_ele.index();
504  bdr.mesh_=this;
505 
506  continue; // next side of element e
507  }
508  }
509 
510  // go through the elements connected to the edge or neighbour
511  for( vector<unsigned int>::iterator isect = intersection_list.begin(); isect!=intersection_list.end(); ++isect) {
512  Element *elem = &(element[*isect]);
513  for (unsigned int ecs=0; ecs<elem->n_sides(); ecs++) {
514  if (elem->edge_idx_[ecs] != Mesh::undef_idx) continue;
515  SideIter si = elem->side(ecs);
516  if ( same_sides( si, side_nodes) ) {
517  if (is_neighbour) {
518  // create a new edge and neighbour for this side, and element to the edge
519  last_edge_idx=edges.size();
520  edges.resize(last_edge_idx+1);
521  edg = &( edges.back() );
522  edg->n_sides = 1;
523  edg->side_ = new struct SideIter[1];
524  edg->side_[0] = si;
525  elem->edge_idx_[ecs] = last_edge_idx;
526 
527  neighbour.edge_idx_ = last_edge_idx;
528 
529  vb_neighbours_.push_back(neighbour); // copy neighbour with this edge setting
530  } else {
531  // connect the side to the edge, and side to the edge
532  edg->side_[ edg->n_sides++ ] = si;
533  elem->edge_idx_[ecs] = last_edge_idx;
534  }
535  break; // next element from intersection list
536  }
537  } // search for side of other connected element
538  } // connected elements
539  OLD_ASSERT( is_neighbour || ( (unsigned int) edg->n_sides ) == intersection_list.size(), "Some connected sides were not found.\n");
540  } // for element sides
541  } // for elements
542 
543  MessageOut().fmt( "Created {} edges and {} neighbours.\n", edges.size(), vb_neighbours_.size() );
544 }
545 
546 
547 
549 {
550  for (EdgeVector::iterator edg=edges.begin(); edg!=edges.end(); edg++)
551  {
552  // side 0 is reference, so its permutation is 0
553  edg->side(0)->element()->permutation_idx_[edg->side(0)->el_idx()] = 0;
554 
555  if (edg->n_sides > 1)
556  {
557  map<const Node*,unsigned int> node_numbers;
558  unsigned int permutation[edg->side(0)->n_nodes()];
559 
560  for (unsigned int i=0; i<edg->side(0)->n_nodes(); i++)
561  node_numbers[edg->side(0)->node(i)] = i;
562 
563  for (int sid=1; sid<edg->n_sides; sid++)
564  {
565  for (unsigned int i=0; i<edg->side(0)->n_nodes(); i++)
566  permutation[node_numbers[edg->side(sid)->node(i)]] = i;
567 
568  switch (edg->side(0)->dim())
569  {
570  case 0:
571  edg->side(sid)->element()->permutation_idx_[edg->side(sid)->el_idx()] = RefElement<1>::permutation_index(permutation);
572  break;
573  case 1:
574  edg->side(sid)->element()->permutation_idx_[edg->side(sid)->el_idx()] = RefElement<2>::permutation_index(permutation);
575  break;
576  case 2:
577  edg->side(sid)->element()->permutation_idx_[edg->side(sid)->el_idx()] = RefElement<3>::permutation_index(permutation);
578  break;
579  }
580  }
581  }
582  }
583 
584  for (vector<Neighbour>::iterator nb=vb_neighbours_.begin(); nb!=vb_neighbours_.end(); nb++)
585  {
586  map<const Node*,unsigned int> node_numbers;
587  unsigned int permutation[nb->element()->n_nodes()];
588 
589  // element of lower dimension is reference, so
590  // we calculate permutation for the adjacent side
591  for (unsigned int i=0; i<nb->element()->n_nodes(); i++)
592  node_numbers[nb->element()->node[i]] = i;
593 
594  for (unsigned int i=0; i<nb->side()->n_nodes(); i++)
595  permutation[node_numbers[nb->side()->node(i)]] = i;
596 
597  switch (nb->side()->dim())
598  {
599  case 0:
600  nb->side()->element()->permutation_idx_[nb->side()->el_idx()] = RefElement<1>::permutation_index(permutation);
601  break;
602  case 1:
603  nb->side()->element()->permutation_idx_[nb->side()->el_idx()] = RefElement<2>::permutation_index(permutation);
604  break;
605  case 2:
606  nb->side()->element()->permutation_idx_[nb->side()->el_idx()] = RefElement<3>::permutation_index(permutation);
607  break;
608  }
609  }
610 }
611 
612 
613 
614 
615 
616 //=============================================================================
617 //
618 //=============================================================================
620 {
621 
622  //MessageOut() << "Element to neighbours of vb2 type... "/*orig verb 5*/;
623 
624  FOR_ELEMENTS(this,ele) ele->n_neighs_vb =0;
625 
626  // count vb neighs per element
627  FOR_NEIGHBOURS(this, ngh ) ngh->element_->n_neighs_vb++;
628 
629  // Allocation of the array per element
630  FOR_ELEMENTS(this, ele )
631  if( ele->n_neighs_vb > 0 ) {
632  ele->neigh_vb = new struct Neighbour* [ele->n_neighs_vb];
633  ele->n_neighs_vb=0;
634  }
635 
636  // fill
637  ElementIter ele;
638  FOR_NEIGHBOURS(this, ngh ) {
639  ele = ngh->element();
640  ele->neigh_vb[ ele->n_neighs_vb++ ] = &( *ngh );
641  }
642 
643  //MessageOut() << "... O.K.\n"/*orig verb 6*/;
644 }
645 
646 
647 
648 
649 
650 
652  /* Algorithm:
653  *
654  * 1) create BIH tree
655  * 2) for every 1D, find list of candidates
656  * 3) compute intersections for 1d, store it to master_elements
657  *
658  */
659  if (! intersections) {
660  intersections = std::make_shared<MixedMeshIntersections>(this);
661  intersections->compute_intersections();
662  }
663  return *intersections;
664 }
665 
666 
667 
668 ElementAccessor<3> Mesh::element_accessor(unsigned int idx, bool boundary) {
669  return ElementAccessor<3>(this, idx, boundary);
670 }
671 
672 
673 
674 void Mesh::elements_id_maps( vector<IdxInt> & bulk_elements_id, vector<IdxInt> & boundary_elements_id) const
675 {
676  if (bulk_elements_id.size() ==0) {
678  IdxInt last_id;
679 
680  bulk_elements_id.resize(n_elements());
681  map_it = bulk_elements_id.begin();
682  last_id = -1;
683  for(unsigned int idx=0; idx < element.size(); idx++, ++map_it) {
684  IdxInt id = element.get_id(idx);
685  if (last_id >= id) xprintf(UsrErr, "Element IDs in non-increasing order, ID: %d\n", id);
686  last_id=*map_it = id;
687  }
688 
689  boundary_elements_id.resize(bc_elements.size());
690  map_it = boundary_elements_id.begin();
691  last_id = -1;
692  for(unsigned int idx=0; idx < bc_elements.size(); idx++, ++map_it) {
693  IdxInt id = bc_elements.get_id(idx);
694  // We set ID for boundary elements created by the mesh itself to "-1"
695  // this force gmsh reader to skip all remaining entries in boundary_elements_id
696  // and thus report error for any remaining data lines
697  if (id < 0) last_id=*map_it=-1;
698  else {
699  if (last_id >= id) xprintf(UsrErr, "Element IDs in non-increasing order, ID: %d\n", id);
700  last_id=*map_it = id;
701  }
702  }
703  }
704 }
705 
707 {
709  it != region_list.end();
710  ++it) {
711  // constructor has side effect in the mesh - create new region or set and store them to Mesh::region_db_
712  (*it).factory< RegionSetBase, const Input::Record &, Mesh * >(*it, this);
713  }
714 }
715 
717 {
719  region_db_.el_to_reg_map_.clear();
720  region_db_.close();
722 
723  if ( in_record_.val<bool>("print_regions") ) {
724  stringstream ss;
726  MessageOut() << ss.str();
727  }
728 }
729 
730 
732  START_TIMER("Mesh::compute_element_boxes");
733  if (element_box_.size() > 0) return;
734 
735  // make element boxes
736  element_box_.resize(this->element.size());
737  unsigned int i=0;
738  FOR_ELEMENTS(this, element) {
739  element_box_[i] = element->bounding_box();
740  i++;
741  }
742 
743  // make mesh box
744  Node* node = this->node_vector.begin();
745  mesh_box_ = BoundingBox(node->point(), node->point());
746  FOR_NODES(this, node ) {
747  mesh_box_.expand( node->point() );
748  }
749 
750 }
751 
753  if (! this->bih_tree_)
754  bih_tree_ = std::make_shared<BIHTree>(this);
755  return *bih_tree_;
756 }
757 
759  return in_record_.val<double>("global_observe_search_radius");
760 }
761 
762 void Mesh::add_physical_name(unsigned int dim, unsigned int id, std::string name) {
763  region_db_.add_region(id, name, dim, "$PhysicalNames");
764 }
765 
766 
767 void Mesh::add_node(unsigned int node_id, arma::vec3 coords) {
768  NodeFullIter node = node_vector.add_item(node_id);
769  node->point() = coords;
770 }
771 
772 
773 void Mesh::add_element(unsigned int elm_id, unsigned int dim, unsigned int region_id, unsigned int partition_id,
774  std::vector<unsigned int> node_ids) {
775  Element *ele=nullptr;
776  RegionIdx region_idx = region_db_.get_region( region_id, dim );
777  if ( !region_idx.is_valid() ) {
778  region_idx = region_db_.add_region( region_id, region_db_.create_label_from_id(region_id), dim, "$Element" );
779  }
780  region_db_.mark_used_region(region_idx.idx());
781 
782  if (region_idx.is_boundary()) {
783  ele = bc_elements.add_item(elm_id);
784  } else {
785  if(dim == 0 ) {
786  WarningOut().fmt("Bulk elements of zero size(dim=0) are not supported. Element ID: {}.\n", elm_id);
787  return;
788  }
789  else
790  ele = element.add_item(elm_id);
791  }
792  ele->init(dim, this, region_idx);
793  ele->pid = partition_id;
794 
795  unsigned int ni;
796  FOR_ELEMENT_NODES(ele, ni) {
797  unsigned int node_id = node_ids[ni];
798  NodeIter node = node_vector.find_id( node_id );
799  INPUT_CHECK( node != node_vector.end(),
800  "Unknown node id %d in specification of element with id=%d.\n", node_id, elm_id);
801  ele->node[ni] = node;
802  }
803 
804  // check that tetrahedron element is numbered correctly and is not degenerated
805  if(ele->dim() == 3)
806  {
807  double jac = ele->tetrahedron_jacobian();
808  if( ! (jac > 0) )
809  WarningOut().fmt("Tetrahedron element with id {} has wrong numbering or is degenerated (Jacobian = {}).",ele->index(),jac);
810  }
811 }
812 
813 
815  if (node_elements_.size() == 0) {
817  }
818  return node_elements_;
819 }
820 
821 //-----------------------------------------------------------------------------
822 // vim: set cindent:
int n_triangles
Definition: mesh.h:297
Distribution * el_ds
Parallel distribution of elements.
Definition: mesh.h:467
Class for the mesh partitioning. This should provide:
Definition: partitioning.hh:43
Iterator< ValueType > begin() const
void read_stream(istream &in, const Type::TypeBase &root_type, FileFormat format)
This method actually reads the given stream in.
int IdxInt
Define integers that are indices into large arrays (elements, nodes, dofs etc.)
Definition: mesh.h:85
Bounding box in 3d ambient space.
Definition: bounding_box.hh:45
vector< vector< unsigned int > > const & node_elements()
Definition: mesh.cc:814
#define FOR_ELEMENT_NODES(i, j)
Definition: elements.h:156
Accessor to input data conforming to declared Array.
Definition: accessors.hh:567
unsigned int * boundary_idx_
Definition: elements.h:94
int get_id(const T *it) const
Definition: sys_vector.hh:458
bool is_boundary() const
Returns true if it is a Boundary region and false if it is a Bulk region.
Definition: region.hh:73
void count_side_types()
Definition: mesh.cc:282
int MPI_Comm
Definition: mpi.h:141
Definition: nodes.hh:32
void check_and_finish()
Definition: mesh.cc:716
MapElementIDToRegionID el_to_reg_map_
Definition: region.hh:570
Reader for (slightly) modified input files.
int pid
Definition: elements.h:87
MixedMeshIntersections & mixed_intersections()
Definition: mesh.cc:651
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
#define FOR_ELEMENTS(_mesh_, __i)
Definition: mesh.h:480
Class for declaration of the input of type Bool.
Definition: type_base.hh:458
void add_node(unsigned int node_id, arma::vec3 coords)
Add new node of given id and coordinates to mesh.
Definition: mesh.cc:767
#define MessageOut()
Macro defining &#39;message&#39; record of log.
Definition: logger.hh:243
unsigned int * permutation_idx_
Definition: elements.h:103
Definition: abscissa.h:25
int n_lines
Definition: mesh.h:296
static const Input::Type::Record & get_input_type()
Definition: partitioning.cc:45
static const unsigned int undef_idx
Definition: mesh.h:125
void create_node_element_lists()
Definition: mesh.cc:295
static unsigned int permutation_index(unsigned int p[n_nodes_per_side])
Definition: ref_element.cc:533
static const Input::Type::Record & get_input_type()
Definition: mesh.cc:71
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
std::shared_ptr< MixedMeshIntersections > intersections
Definition: mesh.h:279
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:416
BoundingBox mesh_box_
Bounding box of whole mesh.
Definition: mesh.h:432
unsigned int index() const
void expand(const Point &point)
Definition: mesh.h:99
double tetrahedron_jacobian() const
Definition: elements.cc:116
int index() const
Definition: sys_vector.hh:78
Input::Record in_record_
Definition: mesh.h:443
string create_label_from_id(unsigned int id) const
Definition: region.cc:337
Node ** node
Definition: elements.h:90
FullIter add_item(int id)
Definition: sys_vector.hh:359
int n_sides
Definition: edges.h:36
vector< vector< vector< unsigned int > > > side_nodes
Definition: mesh.h:309
FullIter find_id(const int id)
Definition: sys_vector.hh:434
Definition: edges.h:26
int n_tetrahedras
Definition: mesh.h:298
ElementAccessor< 3 > element_accessor(unsigned int idx, bool boundary=false)
Definition: mesh.cc:668
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:762
vector< Boundary > boundary_
Definition: mesh.h:264
Partitioning * get_part()
Definition: mesh.cc:218
unsigned int n_sides()
Definition: mesh.cc:199
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:303
Class for declaration of inputs sequences.
Definition: type_base.hh:345
void compute_element_boxes()
Precompute element bounding boxes if it is not done yet.
Definition: mesh.cc:731
int n_exsides
Definition: mesh.h:293
IteratorBase end() const
Region get_region(unsigned int id, unsigned int dim)
Definition: region.cc:150
unsigned int dim() const
unsigned int size() const
Returns size of the container. This is independent of the allocated space.
Definition: sys_vector.hh:391
static Default optional()
The factory function to make an empty default value which is optional.
Definition: type_record.hh:124
#define OLD_ASSERT(...)
Definition: global_defs.h:131
void elements_id_maps(vector< IdxInt > &bulk_elements_id, vector< IdxInt > &boundary_elements_id) const
Definition: mesh.cc:674
unsigned int edge_idx_
Definition: neighbours.h:135
unsigned int * edge_idx_
Definition: elements.h:93
void read_regions_from_input(Input::Array region_list)
Definition: mesh.cc:706
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:132
unsigned int n_elements() const
Definition: mesh.h:156
bool same_sides(const SideIter &si, vector< unsigned int > &side_nodes)
Definition: mesh.cc:356
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:540
std::vector< T >::iterator iterator
Definition: sys_vector.hh:215
void setup_topology()
Definition: mesh.cc:254
IntersectionSearch get_intersection_search()
Getter for input type selection for intersection search algorithm.
Definition: mesh.cc:122
Neighbour ** neigh_vb
Definition: elements.h:135
T get_root_interface() const
Returns the root accessor.
friend class RegionSetBase
Definition: mesh.h:454
static const Input::Type::Selection & get_input_intersection_variant()
The definition of input record for selection of variant of file format.
Definition: mesh.cc:60
static FileName input()
The factory function for declaring type FileName for input files.
Definition: type_base.cc:526
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
#define NDEF
Definition: mesh.cc:56
const Ret val(const string &key) const
unsigned int n_sides() const
#define xprintf(...)
Definition: system.hh:92
#define START_TIMER(tag)
Starts a timer with specified tag.
Class for O(log N) lookup for intersections with a set of bounding boxes.
Definition: bih_tree.hh:36
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.
unsigned int n_corners()
Definition: mesh.cc:208
unsigned int index(const T *pointer) const
Definition: sys_vector.hh:373
#define FOR_NEIGHBOURS(_mesh_, it)
Definition: mesh.h:531
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:490
ElementVector bc_elements
Definition: mesh.h:268
void close()
Definition: region.cc:249
SideIter side(const unsigned int loc_index)
FullIter begin()
Definition: sys_vector.hh:383
bool is_empty() const
Definition: accessors.hh:366
Region implicit_boundary_region()
Definition: region.cc:75
std::shared_ptr< BIHTree > bih_tree_
Definition: mesh.h:437
#define FOR_NODES(_mesh_, i)
Definition: mesh.h:56
vector< vector< unsigned int > > node_elements_
Definition: mesh.h:451
unsigned int edge_idx_
Definition: boundaries.h:75
#define INPUT_CHECK(i,...)
Debugging macros.
Definition: global_defs.h:51
Accessor to the polymorphic input data of a type given by an AbstracRecord object.
Definition: accessors.hh:459
unsigned int bc_ele_idx_
Definition: boundaries.h:76
void init(unsigned int dim, Mesh *mesh_in, RegionIdx reg)
Definition: elements.cc:61
void count_element_types()
Definition: mesh.cc:227
const BIHTree & get_bih_tree()
Getter for BIH. Creates and compute BIH at first call.
Definition: mesh.cc:752
void check_regions()
Definition: region.cc:461
Definition: system.hh:64
void mark_used_region(unsigned int idx)
Definition: region.cc:235
void print_region_table(ostream &stream) const
Definition: region.cc:409
Support classes for parallel programing.
Region add_region(unsigned int id, const std::string &label, unsigned int dim, const std::string &address="implicit")
Definition: region.cc:85
vector< Neighbour > vb_neighbours_
Definition: mesh.h:290
FinishStatus finish(FinishStatus finish_type=FinishStatus::regular_) override
Finish declaration of the Record type.
Definition: type_record.cc:242
int n_sides_
Definition: mesh.h:294
RegionIdx region_idx_
Definition: elements.h:143
const Selection & close() const
Close the Selection, no more values can be added.
IdxInt * row_4_el
Index set assigning to global element index the local index used in parallel vectors.
Definition: mesh.h:463
RegionDB region_db_
Definition: mesh.h:422
std::vector< Edge > edges
Vector of MH edges, this should not be part of the geometrical mesh.
Definition: mesh.h:271
MPI_Comm comm_
Definition: mesh.h:448
unsigned int n_neighs_vb
Definition: elements.h:133
int n_insides
Definition: mesh.h:292
void intersect_element_lists(vector< unsigned int > const &nodes_list, vector< unsigned int > &intersection_element_list)
Definition: mesh.cc:308
void modify_element_ids(const RegionDB::MapElementIDToRegionID &map)
Definition: mesh.cc:244
#define WarningOut()
Macro defining &#39;warning&#39; record of log.
Definition: logger.hh:246
void make_edge_permutations()
Definition: mesh.cc:548
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc...
static Input::Type::Abstract & get_input_type()
Definition: region_set.cc:25
Mesh * mesh_
Definition: boundaries.h:77
unsigned int n_nodes() const
Definition: mesh.h:152
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:773
Record type proxy class.
Definition: type_record.hh:182
ElementIter element_
Definition: neighbours.h:136
bool is_valid() const
Returns false if the region has undefined/invalid value.
Definition: region.hh:77
IntersectionSearch
Types of search algorithm for finding intersection candidates.
Definition: mesh.h:114
SideIter * side_
Definition: edges.h:37
void element_to_neigh_vb()
Definition: mesh.cc:619
const Node * node(unsigned int i) const
Definition: side_impl.hh:46
Mesh()
Definition: mesh.cc:94
arma::vec3 & point()
Definition: nodes.hh:68
void reinit(Input::Record in_record)
Definition: mesh.cc:128
FullIter end()
Returns FullFullIterer of the fictions past the end element.
Definition: sys_vector.hh:387
#define FOR_SIDES(_mesh_, it)
Definition: mesh.h:524
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
bool find_lower_dim_element(ElementVector &elements, vector< unsigned int > &element_list, unsigned int dim, unsigned int &element_idx)
Definition: mesh.cc:336
Template for classes storing finite set of named values.
std::shared_ptr< Partitioning > part_
Definition: mesh.h:426
void make_neighbours_and_edges()
Definition: mesh.cc:373
NodeVector node_vector
Vector of nodes of the mesh.
Definition: mesh.h:258
double global_observe_radius() const
Maximal distance of observe point from Mesh relative to its size.
Definition: mesh.cc:758
Main class for computation of intersection of meshes of combined dimensions.
ElementVector element
Vector of elements of the mesh.
Definition: mesh.h:260
unsigned int n_nodes() const
Definition: side_impl.hh:33
~Mesh()
Destructor.
Definition: mesh.cc:173
unsigned int idx() const
Returns a global index of the region.
Definition: region.hh:81
std::vector< BoundingBox > element_box_
Auxiliary vector of mesh elements bounding boxes.
Definition: mesh.h:429
IdxInt * el_4_loc
Index set assigning to local element index its global index.
Definition: mesh.h:465