Flow123d  release_3.0.0-1190-g0f314ad
inspect_elements_algorithm.cc
Go to the documentation of this file.
1 /*
2  *
3  *
4  * Created on: 29.4.2016
5  * Author: pe
6  */
7 
8 #include <unordered_set>
9 #include <boost/functional/hash.hpp>
10 
13 #include "intersection_aux.hh"
14 #include "intersection_local.hh"
15 #include "compute_intersection.hh"
16 
17 #include "system/sys_profiler.hh"
18 
19 #include "mesh/side_impl.hh"
20 #include "mesh/mesh.h"
21 #include "mesh/ref_element.hh"
22 #include "mesh/bih_tree.hh"
23 #include "mesh/accessors.hh"
24 #include "mesh/range_wrapper.hh"
25 
26 
27 
28 template<unsigned int dimA, unsigned int dimB>
30 : mesh(mesh)
31 {}
32 
33 // template<unsigned int dimA, unsigned int dimB>
34 // template<unsigned int simplex_dim>
35 // void IntersectionAlgorithmBase<dimA,dimB>::update_simplex(const ElementAccessor<3>& element, Simplex< simplex_dim >& simplex)
36 // {
37 // ASSERT(simplex_dim == element.dim());
38 // arma::vec3 *field_of_points[simplex_dim+1];
39 // for(unsigned int i=0; i < simplex_dim+1; i++)
40 // field_of_points[i]= &(element.node(i)->point());
41 // simplex.set_simplices(field_of_points);
42 // }
43 
44 
45 template<unsigned int dim>
47 : IntersectionAlgorithmBase<dim,3>(input_mesh)
48 {
49 }
50 
51 template<unsigned int dim>
53 {}
54 
55 
56 
57 template<unsigned int dim>
59 {
60  START_TIMER("Intersection initialization");
62  closed_elements.assign(mesh->n_elements(), false);
64  n_intersections_ = 0;
65  END_TIMER("Intersection initialization");
66 }
67 
68 template<unsigned int dim>
70 {
71  START_TIMER("Compute bounding boxes");
72  if(elements_bb.size() == 0){
73  elements_bb.resize(mesh->n_elements());
74  bool first_3d_element = true;
75  for (auto elm : mesh->elements_range()) {
76 
77  elements_bb[elm.idx()] = elm.bounding_box();
78 
79  if (elm->dim() == 3){
80  if(first_3d_element){
81  first_3d_element = false;
82  mesh_3D_bb = elements_bb[elm.idx()];
83  }else{
84  mesh_3D_bb.expand(elements_bb[elm.idx()].min());
85  mesh_3D_bb.expand(elements_bb[elm.idx()].max());
86  }
87  }
88  }
89  }
90  END_TIMER("Compute bounding boxes");
91 }
92 
93 template<unsigned int dim>
95  const ElementAccessor<3> &bulk_ele)
96 {
97  unsigned int component_ele_idx = comp_ele.idx(),
98  bulk_ele_idx = bulk_ele.idx();
99 
100  IntersectionAux<dim,3> is(component_ele_idx, bulk_ele_idx);
101  START_TIMER("Compute intersection");
102  ComputeIntersection<dim,3> CI(comp_ele, bulk_ele, mesh);
103  CI.init();
104  CI.compute(is);
105  END_TIMER("Compute intersection");
106 
107  last_slave_for_3D_elements[bulk_ele_idx] = component_ele_idx;
108 
109  if(is.points().size() > 0) {
110  intersection_list_[component_ele_idx].push_back(is);
112  return true;
113  }
114  else return false;
115 }
116 
117 
118 
119 
120 template<unsigned int dim>
121 bool InspectElementsAlgorithm<dim>::intersection_exists(unsigned int component_ele_idx, unsigned int bulk_ele_idx)
122 {
123  START_TIMER("intersection_exists");
124  bool found = false;
125 
126  START_TIMER("cycle");
127  unsigned int i = 0;
128  for(; i < intersection_list_[component_ele_idx].size();i++){
129 
130  if(intersection_list_[component_ele_idx][i].bulk_ele_idx() == bulk_ele_idx){
131  found = true;
132  break;
133  }
134  }
135  ADD_CALLS(i);
136  END_TIMER("cycle");
137  END_TIMER("intersection_exists");
138  return found;
139 }
140 
141 
142 template<unsigned int dim>
144 {
145  //DebugOut() << "######### ALGORITHM: compute_intersections #########\n";
146 
147  init();
148 
149  START_TIMER("Element iteration");
150 
151  for (auto elm : mesh->elements_range()) {
152  unsigned int component_ele_idx = elm.idx();
153 
154  if (elm->dim() == dim && // is component element
155  !closed_elements[component_ele_idx] && // is not closed yet
156  bih.ele_bounding_box(component_ele_idx).intersect(bih.tree_box())) // its bounding box intersects 3D mesh bounding box
157  {
158  std::vector<unsigned int> searchedElements;
159 
160  START_TIMER("BIHtree find");
161  bih.find_bounding_box(bih.ele_bounding_box(component_ele_idx), searchedElements);
162  END_TIMER("BIHtree find");
163 
164  START_TIMER("Bounding box element iteration");
165 
166  // Go through all element which bounding box intersects the component element bounding box
167  for (std::vector<unsigned int>::iterator it = searchedElements.begin(); it!=searchedElements.end(); it++)
168  {
169  unsigned int bulk_ele_idx = *it;
170  ElementAccessor<3> ele_3D = mesh->element_accessor( bulk_ele_idx );
171 
172  // if:
173  // check 3D only
174  // check with the last component element computed for the current 3D element
175  // intersection has not been computed already
176  if (ele_3D->dim() == 3 &&
177  (last_slave_for_3D_elements[bulk_ele_idx] != component_ele_idx &&
178  !intersection_exists(component_ele_idx,bulk_ele_idx) )
179  ) {
180  // check that tetrahedron element is numbered correctly and is not degenerated
181  ASSERT_DBG(ele_3D.tetrahedron_jacobian() > 0).add_value(ele_3D.index(),"element index").error(
182  "Tetrahedron element (%d) has wrong numbering or is degenerated (negative Jacobian).");
183 
184  // - find first intersection
185  // - if found, prolongate and possibly fill both prolongation queues
186  // do-while loop:
187  // - empty prolongation queues:
188  // - empty bulk queue:
189  // - get a candidate from queue and compute CI
190  // - prolongate and possibly push new candidates into queues
191  // - repeat until bulk queue is empty
192  // - the component element is still the same whole time in here
193  //
194  // - the component element might get fully covered by bulk elements
195  // and only then it can be closed
196  //
197  // - pop next candidate from component queue:
198  // - the component element is now changed
199  // - compute CI
200  // - prolongate and possibly push new candidates into queues
201  //
202  // - repeat until both queues are empty
203 
204  bool found = compute_initial_CI(elm, ele_3D);
205 
206  // keep the index of the current component element that is being investigated
207  unsigned int current_component_element_idx = component_ele_idx;
208 
209  if(found){
210 
211  prolongation_decide(elm, ele_3D, intersection_list_[component_ele_idx].back());
212 
213  START_TIMER("Prolongation algorithm");
214  do{
215  // flag is set false if the component element is not fully covered with tetrahedrons
216  bool element_covered = true;
217 
218  while(!bulk_queue_.empty()){
219  Prolongation pr = bulk_queue_.front();
220  //DebugOut().fmt("Bulk queue: ele_idx {}.\n",pr.elm_3D_idx);
221 
222  if( pr.elm_3D_idx == undefined_elm_idx_)
223  {
224  //DebugOut().fmt("Open intersection component element: {}\n",current_component_element_idx);
225  element_covered = false;
226  }
227  else prolongate(pr);
228 
229  bulk_queue_.pop();
230  }
231 
232  if(! closed_elements[current_component_element_idx])
233  closed_elements[current_component_element_idx] = element_covered;
234 
235 
236  if(!component_queue_.empty()){
237  Prolongation pr = component_queue_.front();
238 
239  // note the component element index
240  current_component_element_idx = pr.component_elm_idx;
241  //DebugOut().fmt("Component queue: ele_idx {}.\n",current_component_element_idx);
242 
243  prolongate(pr);
244  component_queue_.pop();
245  }
246  }
247  while( !(component_queue_.empty() && bulk_queue_.empty()) );
248  END_TIMER("Prolongation algorithm");
249 
250  // if component element is closed, do not check other bounding boxes
251  if(closed_elements[component_ele_idx])
252  break;
253  }
254  }
255  }
256  END_TIMER("Bounding box element iteration");
257  }
258  }
259 
260  END_TIMER("Element iteration");
261 
262  MessageOut().fmt("{}D-3D: number of intersections = {}\n", dim, n_intersections_);
263  // DBG write which elements are closed
264 // for (auto ele : mesh->elements_range()) {
265 // DebugOut().fmt("Element[{}] closed: {}\n",ele.index(),(closed_elements[ele.index()] ? 1 : 0));
266 // }
267 }
268 
269 template<unsigned int dim>
271 {
272  DebugOut() << "######### ALGORITHM: compute_intersections_BIHtree #########\n";
273 
274  init();
275 
276  START_TIMER("Element iteration");
277 
278  for (auto elm : mesh->elements_range()) {
279  unsigned int component_ele_idx = elm.idx();
280 
281  if (elm.dim() == dim && // is component element
282  bih.ele_bounding_box(component_ele_idx).intersect(bih.tree_box())) // its bounding box intersects 3D mesh bounding box
283  {
284  std::vector<unsigned int> searchedElements;
285 
286  START_TIMER("BIHtree find");
287  bih.find_bounding_box(bih.ele_bounding_box(component_ele_idx), searchedElements);
288  END_TIMER("BIHtree find");
289 
290  START_TIMER("Bounding box element iteration");
291 
292  // Go through all element which bounding box intersects the component element bounding box
293  for (std::vector<unsigned int>::iterator it = searchedElements.begin(); it!=searchedElements.end(); it++)
294  {
295  unsigned int bulk_ele_idx = *it;
296  ElementAccessor<3> ele_3D = mesh->element_accessor( bulk_ele_idx );
297 
298  if (ele_3D.dim() == 3
299  ) {
300  // check that tetrahedron element is numbered correctly and is not degenerated
301  ASSERT_DBG(ele_3D.tetrahedron_jacobian() > 0).add_value(ele_3D.idx(),"element index").error(
302  "Tetrahedron element (%d) has wrong numbering or is degenerated (negative Jacobian).");
303 
304  IntersectionAux<dim,3> is(component_ele_idx, bulk_ele_idx);
305  START_TIMER("Compute intersection");
306  ComputeIntersection<dim,3> CI(elm, ele_3D, mesh);
307  CI.init();
308  CI.compute(is);
309  END_TIMER("Compute intersection");
310 
311  if(is.points().size() > 0) {
312 
313  intersection_list_[component_ele_idx].push_back(is);
315  // if component element is closed, do not check other bounding boxes
316  closed_elements[component_ele_idx] = true;
317  }
318  }
319  }
320  END_TIMER("Bounding box element iteration");
321  }
322  }
323 
324  END_TIMER("Element iteration");
325 }
326 
327 template<unsigned int dim>
329 {
330  //DebugOut() << "######### ALGORITHM: compute_intersections_BB #########\n";
331  init();
333 
334  START_TIMER("Element iteration");
335 
336 
337  for (auto elm : mesh->elements_range()) {
338  unsigned int component_ele_idx = elm.idx();
339 
340  if (elm.dim() == dim && // is component element
341  !closed_elements[component_ele_idx] && // is not closed yet
342  elements_bb[component_ele_idx].intersect(mesh_3D_bb)) // its bounding box intersects 3D mesh bounding box
343  {
344 
345  START_TIMER("Bounding box element iteration");
346 
347  // Go through all element which bounding box intersects the component element bounding box
348  for (auto ele_3D : mesh->elements_range()) {
349  unsigned int bulk_ele_idx = ele_3D.idx();
350 
351  // if:
352  // check 3D only
353  // check with the last component element computed for the current 3D element
354  // check that the bounding boxes intersect
355  // intersection has not been computed already
356  if (ele_3D.dim() == 3 &&
357  (last_slave_for_3D_elements[bulk_ele_idx] != component_ele_idx &&
358  elements_bb[component_ele_idx].intersect(elements_bb[bulk_ele_idx]) &&
359  !intersection_exists(component_ele_idx,bulk_ele_idx) )
360  ){
361  // check that tetrahedron element is numbered correctly and is not degenerated
362  ASSERT_DBG(ele_3D.tetrahedron_jacobian() > 0).add_value(ele_3D.index(),"element index").error(
363  "Tetrahedron element (%d) has wrong numbering or is degenerated (negative Jacobian).");
364 
365  // - find first intersection
366  // - if found, prolongate and possibly fill both prolongation queues
367  // do-while loop:
368  // - empty prolongation queues:
369  // - empty bulk queue:
370  // - get a candidate from queue and compute CI
371  // - prolongate and possibly push new candidates into queues
372  // - repeat until bulk queue is empty
373  // - the component element is still the same whole time in here
374  //
375  // - the component element might get fully covered by bulk elements
376  // and only then it can be closed
377  //
378  // - pop next candidate from component queue:
379  // - the component element is now changed
380  // - compute CI
381  // - prolongate and possibly push new candidates into queues
382  //
383  // - repeat until both queues are empty
384 
385  bool found = compute_initial_CI(elm, ele_3D);
386 
387  // keep the index of the current component element that is being investigated
388  unsigned int current_component_element_idx = component_ele_idx;
389 
390  if(found){
391  //DebugOut().fmt("start component with elements {} {}\n",component_ele_idx, bulk_ele_idx);
392 
393  prolongation_decide(elm, ele_3D, intersection_list_[component_ele_idx].back());
394 
395  START_TIMER("Prolongation algorithm");
396  do{
397  // flag is set false if the component element is not fully covered with tetrahedrons
398  bool element_covered = true;
399 
400  while(!bulk_queue_.empty()){
401  Prolongation pr = bulk_queue_.front();
402  //DebugOut().fmt("Bulk queue: ele_idx {}.\n",pr.elm_3D_idx);
403 
404  if( pr.elm_3D_idx == undefined_elm_idx_)
405  {
406  //DebugOut().fmt("Open intersection component element: {}\n",current_component_element_idx);
407  element_covered = false;
408  }
409  else prolongate(pr);
410 
411  bulk_queue_.pop();
412  }
413 
414  if(! closed_elements[current_component_element_idx])
415  closed_elements[current_component_element_idx] = element_covered;
416 
417 
418  if(!component_queue_.empty()){
419  Prolongation pr = component_queue_.front();
420 
421  // note the component element index
422  current_component_element_idx = pr.component_elm_idx;
423  //DebugOut().fmt("Component queue: ele_idx {}.\n",current_component_element_idx);
424 
425  prolongate(pr);
426  component_queue_.pop();
427  }
428  }
429  while( !(component_queue_.empty() && bulk_queue_.empty()) );
430  END_TIMER("Prolongation algorithm");
431 
432  // if component element is closed, do not check other bounding boxes
433  if(closed_elements[component_ele_idx])
434  break;
435  }
436 
437  }
438  }
439  END_TIMER("Bounding box element iteration");
440  }
441  }
442 
443  END_TIMER("Element iteration");
444 
445  // DBG write which elements are closed
446 // for (auto ele : mesh->elements_range()) {
447 // DebugOut().fmt("Element[{}] closed: {}\n",ele.index(),(closed_elements[ele.index()] ? 1 : 0));
448 // }
449 }
450 
451 
452 
453 template<unsigned int dim>
454 template<unsigned int ele_dim>
456  unsigned int ip_dim,
457  unsigned int ip_obj_idx)
458 {
459  std::vector<Edge*> edges;
460  edges.reserve(ele_dim - ip_dim); // reserve number of possible edges
461 
462  //DebugOut() << "dim " << ele_dim << ": ";
463  switch (ip_dim)
464  {
465  // IP is at a node of tetrahedron; possible edges are from all connected sides (3)
466  case 0: if(ele_dim == 1) {
467  edges.push_back(&(mesh->edges[ele->edge_idx(ip_obj_idx)]));
468  break;
469  }
470 
471  for(unsigned int j=0; j < RefElement<ele_dim>::n_sides_per_node; j++){
472  unsigned int local_edge = RefElement<ele_dim>::interact(Interaction<ele_dim-1,0>(ip_obj_idx))[j];
473  edges.push_back(&(mesh->edges[ele->edge_idx(local_edge)]));
474  }
475  //DebugOut() << "prolong (node)\n";
476  break;
477 
478  // IP is on a line of tetrahedron; possible edges are from all connected sides (2)
479  case 1: if(ele_dim == 2) {
480  edges.push_back(&(mesh->edges[ele->edge_idx(ip_obj_idx)]));
481  break;
482  }
483 
484  ASSERT_DBG(ele_dim == 3);
485  for(unsigned int j=0; j < RefElement<ele_dim>::n_sides_per_line; j++){
486  unsigned int local_edge = RefElement<ele_dim>::interact(Interaction<2,1>(ip_obj_idx))[j];
487  edges.push_back(&(mesh->edges[ele->edge_idx(local_edge)]));
488  }
489  //DebugOut() << "prolong (edge)\n";
490  break;
491 
492  // IP is on a side of tetrahedron; only possible edge is from the given side (1)
493  case 2: ASSERT_DBG(ele_dim == 3);
494  edges.push_back(&(mesh->edges[ele->edge_idx(ip_obj_idx)]));
495  //DebugOut() << "prolong (side)\n";
496  break;
497  default: ASSERT_DBG(0);
498  }
499 
500  // get indices of neighboring bulk elements
501  std::vector<unsigned int> elements_idx;
502  elements_idx.reserve(2*edges.size()); // twice the number of edges
503  for(Edge* edg : edges)
504  for(int j=0; j < edg->n_sides;j++) {
505  if ( edg->side(j)->element().idx() != ele.idx() )
506  elements_idx.push_back(edg->side(j)->element().idx());
507  }
508 
509  return elements_idx;
510 }
511 
512 
513 template<unsigned int dim>
514 unsigned int InspectElementsAlgorithm<dim>::create_prolongation(unsigned int bulk_ele_idx,
515  unsigned int component_ele_idx,
516  std::queue< Prolongation >& queue)
517 {
518 // if(last_slave_for_3D_elements[bulk_ele_idx] == undefined_elm_idx_ ||
519 // (last_slave_for_3D_elements[bulk_ele_idx] != component_ele_idx && !intersection_exists(component_ele_idx,bulk_ele_idx)))
520 // {
521  last_slave_for_3D_elements[bulk_ele_idx] = component_ele_idx;
522 
523  //DebugOut().fmt("prolongation: c {} in b {}\n",component_ele_idx,bulk_ele_idx);
524 
525  // prepare empty intersection object
526  IntersectionAux<dim,3> il_other(component_ele_idx, bulk_ele_idx);
527  intersection_list_[component_ele_idx].push_back(il_other);
528 
529  Prolongation pr = {component_ele_idx, bulk_ele_idx, (unsigned int)intersection_list_[component_ele_idx].size() - 1};
530  queue.push(pr);
531 
532  return 1;
533 // }
534 // return 0;
535 }
536 
537 
538 
539 template<unsigned int dim>
541  const ElementAccessor<3>& bulk_ele,
543 {
544  //DebugOut() << "DECIDE\n";
545  // number of IPs that are at vertices of component element (counter used for closing element)
546  unsigned int n_ip_vertices = 0;
547 
548  if(dim == 2){
549  // check whether all IPs lie in the face
550  unsigned int sid = is.ips_in_face();
551  if(sid < RefElement<3>::count<2>()){
552  //NOTE: we are unable to combine cases A and B in 1d-2d
553  // CASE A: compatible vb neighboring
554 // DBGVAR(comp_ele.idx());
555 // DBGVAR(bulk_ele.idx());
556  if(comp_ele->n_neighs_vb() > 1){
557  //set n_neighs_vb duplicities
558  is.set_duplicities(comp_ele->n_neighs_vb());
559 // DBGCOUT(<< "vb neigh: N = " << is.duplicities() << "\n");
560  //possibly copy intersection object
561  }
562  // CASE B: incompatible neighboring with lower dim element
563  else {
564  const Edge* edg = bulk_ele.side(sid)->edge();
565  if(edg->n_sides > 1){
566  //set n_sides duplicities
567  is.set_duplicities(edg->n_sides);
568 // DBGCOUT(<< "incomp. neigh: N = " << is.duplicities() << "\n");
569  //possibly copy intersection object
570  }
571  }
572  }
573  }
574 
575  for(const IntersectionPointAux<dim,3> &IP : is.points()) {
576 
577  // 1) prolong over component, if IP is at its boundary
578  // 2) prolong over bulk, if IP is at its boundary
579  // both cases are possible, so both prolongations might happen at once..
580 
581  if(IP.dim_A() < dim) { // if IP on the boundary of component element
582  if(IP.dim_A() == 0) n_ip_vertices++;
583  //DebugOut() << "on " << dim << "D boundary, dim = " << IP.dim_A() << "\n";
584 
585  // search for indices of neighboring component elements (including the current one)
586  std::vector<unsigned int> comp_neighbors = get_element_neighbors<dim>(comp_ele,IP.dim_A(), IP.idx_A());
587 
588 // DBGCOUT( << comp_ele.idx() << "--" << bulk_ele.idx() << ": ");
589 // for(unsigned int& comp_neighbor_idx : comp_neighbors)
590 // cout << mesh->element_accessor(comp_neighbor_idx).idx() << " ";
591 // cout << "\n";
592 
593  unsigned int bulk_current = bulk_ele.idx();
594 
595  // add all component neighbors with current bulk element into component queue
596  for(unsigned int& comp_neighbor_idx : comp_neighbors) {
597  if(!intersection_exists(comp_neighbor_idx,bulk_current))
598  create_prolongation(bulk_current, comp_neighbor_idx, component_queue_);
599  }
600  }
601 
602  if(IP.dim_B() < 3)
603  {
604  //DebugOut() << "on 3D boundary, dim = " << IP.dim_B() << "\n";
605 
606  // search for indices of neighboring bulk elements (including the current one)
607  std::vector<unsigned int> bulk_neighbors = get_element_neighbors<3>(bulk_ele,IP.dim_B(),IP.idx_B());
608 
609 // DBGCOUT( << comp_ele.idx() << "--" << bulk_ele.idx() << ": ");
610 // for(unsigned int& bulk_neighbor_idx : bulk_neighbors)
611 // cout << mesh->element_accessor(bulk_neighbor_idx).idx() << " ";
612 // cout << "\n";
613 
614  unsigned int comp_current = comp_ele.idx();
615  unsigned int n_prolongations = 0;
616  // prolong over current comp element to other bulk elements (into bulk queue) (covering comp ele)
617  for(unsigned int& bulk_neighbor_idx : bulk_neighbors)
618  {
619  if(last_slave_for_3D_elements[bulk_neighbor_idx] == undefined_elm_idx_ ||
620  (last_slave_for_3D_elements[bulk_neighbor_idx] != comp_current &&
621  !intersection_exists(comp_current,bulk_neighbor_idx)))
622  n_prolongations += create_prolongation(bulk_neighbor_idx,
623  comp_current,
624  bulk_queue_);
625  }
626 
627  // if there are no sides of any edge that we can continue to prolongate over,
628  // it means we are at the boundary and cannot prolongate further
629  if(n_prolongations == 0)
630  {
632  bulk_queue_.push(pr);
633  }
634  }
635  }
636 
637  // close component element if it has all vertices inside bulk element
638  if(n_ip_vertices == is.size()) closed_elements[comp_ele.idx()] = true;
639 }
640 
641 
642 
643 template<>
644 void InspectElementsAlgorithm<2>::assert_same_intersection(unsigned int comp_ele_idx, unsigned int bulk_ele_idx)
645 {
646  for(unsigned int i=0; i < intersection_list_[comp_ele_idx].size(); i++)
647  {
648  if(intersection_list_[comp_ele_idx][i].bulk_ele_idx() == bulk_ele_idx)
649  {
650  //DebugOut().fmt("intersection comp-bulk: {} {}\n", comp_ele_idx, bulk_ele_idx);
651  ASSERT_DBG(0).add_value(bulk_ele_idx,"bulk_ele_idx").error("Want to add the same intersection!");
652  }
653  }
654 }
655 
656 
657 template<unsigned int dim>
659 {
662 
663 // DebugOut().fmt("Prolongate: {} in {}.\n", elm.idx(), ele_3D.idx());
664 
665  //TODO: optimization: this might be called before and not every time
666  //(component element is not changing when emptying bulk queue)
667 // this->update_simplex(elm, simplexA);
668 // this->update_simplex(ele_3D, simplexB);
669 
671 
672  START_TIMER("Compute intersection");
673  ComputeIntersection<dim,3> CI(elm, ele_3D, mesh);
674  CI.init();
675  CI.compute(is);
676  END_TIMER("Compute intersection");
677 
679 
680  if(is.size() > 0){
681 // for(unsigned int j=0; j < is.size(); j++)
682 // DebugOut() << is[j];
683 // DebugOut().fmt("intersection of elements {} {} [{}--{}] size {}\n",
684 // elm.idx(), ele_3D.idx(),
685 // elm.region().label(), ele_3D.region().label(),
686 // is.size()
687 // );
688 
689  prolongation_decide(elm, ele_3D,is);
691 // DBGVAR(n_intersections_);
692  }
693  else{
694  // NOTE: we get here, when create_prolongation creates an empty intersection
695  // - it can happen, that the CI will be empty
696  // - currently, we remove it when storing the final intersection objects
697  // - or we can erase it at this point
698 // WarningOut() << "zero is: c " << elm.index() << " b " << ele_3D.index();
699 // auto & v = intersection_list_[pr.component_elm_idx];
700 // v.erase( v.next(v.being(),pr.dictionary_idx) );
701  }
702 }
703 
704 
705 
706 
708 : IntersectionAlgorithmBase<2,2>(input_mesh)
709 {}
710 
711 
714 {
715 // DebugOut() << "Intersections 2d-2d\n";
716  ASSERT(storage.size() == 0);
718 
719  unsigned int ele_idx, eleA_idx, eleB_idx,
720  componentA_idx, componentB_idx,
721  temp_eleA_idx;
722 
723  typedef std::pair<unsigned int, unsigned int> ipair;
724  std::unordered_set<ipair, boost::hash<ipair>> computed_pairs;
725 
726  for (auto ele : mesh->elements_range()) {
727  if (ele->dim() == 3)
728  {
729  ele_idx = ele.idx();
730  // if there are not at least 2 2D elements intersecting 3D element; continue
731  if(intersection_map[ele_idx].size() < 2) continue;
732 
733  const std::vector<ILpair> &local_map = intersection_map[ele_idx];
734 
735 // DebugOut() << print_var(local_map.size());
736  for(unsigned int i=0; i < local_map.size(); i++)
737  {
738  //TODO: 1] compute all plucker coords at once
739  //TODO: 2] pass plucker coords from 2d-3d
740 
741  eleA_idx = local_map[i].first;
742  ElementAccessor<3> eleA = mesh->element_accessor( eleA_idx );
743  if(eleA->dim() !=2 ) continue; //skip other dimension intersection
744  componentA_idx = component_idx_[eleA_idx];
745 
746 // IntersectionLocalBase * ilb = local_map[i].second;
747 // DebugOut().fmt("2d-2d ILB: {} {} {}\n", ilb->bulk_ele_idx(), ilb->component_ele_idx(), componentA_idx);
748 
749  for(unsigned int j=i+1; j < local_map.size(); j++)
750  {
751  eleB_idx = local_map[j].first;
752  componentB_idx = component_idx_[eleB_idx];
753 
754  if(componentA_idx == componentB_idx) continue; //skip elements of the same component
755  // this also skips the compatible connections (it is still a single component in this case)
756 
757  ElementAccessor<3> eleB = mesh->element_accessor( eleB_idx );
758  if(eleB->dim() !=2 ) continue; //skip other dimension intersection
759 
760  // set master -- slave order
761  // do not overwrite the original eleA
762  temp_eleA_idx = eleA_idx;
763  ElementAccessor<3> temp_eleA = eleA;
764  if (componentA_idx < componentB_idx){
765  std::swap(temp_eleA_idx, eleB_idx);
766  std::swap(temp_eleA, eleB);
767  }
768 
769  //skip candidates already computed
770  ipair ip = std::make_pair(temp_eleA_idx, eleB_idx);
771  if(computed_pairs.count(ip) == 1){
772 // DBGCOUT(<< "skip: " << eleA_idx << " " << eleB_idx << "\n");
773  continue;
774  }
775  else{
776  compute_single_intersection(temp_eleA, eleB, storage);
777  computed_pairs.emplace(ip);
778  }
779 
780 // bool skip = false;
781 // for(unsigned int k=0; k<intersection_map[eleA_idx].size(); k++)
782 // {
783 // if(intersection_map[eleA_idx][k].first == eleB_idx) {
784 // skip = true;
785 // break;
786 // }
787 // }
788 // if(skip) continue;
789 
790 // DebugOut().fmt("compute intersection 2d-2d: e_{} e_{} c_{} c_{}\n",
791 // eleA.index(), eleB.index(), componentA_idx, componentB_idx);
792 // DebugOut().fmt("compute intersection 2d-2d: e_{} e_{} c_{} c_{}\n",
793 // eleA.idx(), eleB.idx(), componentA_idx, componentB_idx);
794 
795 // IntersectionAux<2,2> is;
796 // if (componentA_idx < componentB_idx)
797 // compute_single_intersection(eleA, eleB, storage);
798 // else
799 // compute_single_intersection(eleB, eleA, storage);
800  }
801  }
802  }
803  }
804  MessageOut() << "2D-2D: number of intersections = " << storage.size() << "\n";
805 }
806 
807 
809  const ElementAccessor<3>& eleB,
811 {
812  ASSERT_DBG(eleA.dim() == 2);
813  ASSERT_DBG(eleB.dim() == 2);
814  ASSERT_DBG(eleA.idx() != eleB.idx());
815 
816  IntersectionAux<2,2> is(eleA.idx(), eleB.idx());
817 
818  ComputeIntersection<2,2> CI(eleA, eleB, mesh);
819  CI.init();
820  unsigned int n_local_intersection = CI.compute(is);
821 
822  // do not store point intersections
823  if(n_local_intersection > 1){
824  storage.push_back(IntersectionLocal<2,2>(is));
825  }
826 
827 }
828 
830 {
832  component_counter_ = 0;
833 
834  // prolongation queue in the component mesh.
835  std::queue<unsigned int> queue;
836 
837  for (auto ele : mesh->elements_range()) {
838  if (ele->dim() == 2 &&
839  component_idx_[ele.idx()] == (unsigned int)-1)
840  {
841  // start component
842  queue.push(ele.idx());
843 
844  while(!queue.empty()){
845  unsigned int ele_idx = queue.front();
846  queue.pop();
847  const ElementAccessor<3>& elm = mesh->element_accessor( ele_idx );
848  for(unsigned int sid=0; sid < elm->n_sides(); sid++) {
849  const Edge* edg = elm.side(sid)->edge();
850 
851  for(int j=0; j < edg->n_sides;j++) {
852  uint neigh_idx = edg->side(j)->element().idx();
853  if (component_idx_[neigh_idx] == (unsigned int)-1) {
854  component_idx_[neigh_idx] = component_counter_;
855  queue.push(neigh_idx);
856  }
857  }
858  }
859  }
861  }
862  }
863 
864  MessageOut() << "2D-2D: number of components = " << component_counter_ << "\n";
865 
866 // DBGCOUT(<< "Component numbering: \n");
867 // for (auto ele : mesh->elements_range()) {
868 // if (ele->dim() == 2){
869 // cout << "2d ele " << ele.index() << ": " << component_idx_[ele.index()] << endl;
870 // }
871 // }
872 }
873 
874 
875 
876 
877 
879 : IntersectionAlgorithmBase<1,2>(input_mesh)
880 {}
881 
882 
885 {
886  //DebugOut() << "Intersections 1d-2d\n";
888  ASSERT(storage.size() == 0);
889 
890  for (auto ele : mesh->elements_range()) {
891  if (ele->dim() == 3)
892  {
893  unsigned int ele_idx = ele.idx();
894  // if there are not at least 2 elements intersecting 3D element; continue
895  if(intersection_map[ele_idx].size() < 2) continue;
896 
897  const std::vector<ILpair> &local_map = intersection_map[ele_idx];
898 
899  //DebugOut() << "more than 2 intersections in tetrahedron found\n";
900  for(unsigned int i=0; i < local_map.size(); i++)
901  {
902  //TODO: 1] compute all plucker coords at once
903  //TODO: 2] pass plucker coords from 1d-3d
904 
905  unsigned int eleA_idx = local_map[i].first;
906  ElementAccessor<3> eleA = mesh->element_accessor( eleA_idx );
907 
908  if(eleA.dim() !=1 ) continue; //skip other dimension intersection
909 
910  for(unsigned int j=0; j < local_map.size(); j++)
911  {
912  unsigned int eleB_idx = local_map[j].first;
913  ElementAccessor<3> eleB = mesh->element_accessor( local_map[j].first );
914  if(eleB.dim() !=2 ) continue; //skip other dimension intersection
915 
916  //skip candidates already computed
917  bool skip = false;
918  for(unsigned int i=0; i<intersection_map[eleA_idx].size(); i++)
919  {
920  if(intersection_map[eleA_idx][i].first == eleB_idx) {
921  skip = true;
922  break;
923  }
924  }
925 
926  if(skip) continue;
927 
928  //DebugOut().fmt("compute intersection 1d-2d: {} {}\n",eleA.index(), eleB.index());
929 // compute_single_intersection(eleA,
930 // eleB);
931 
932  IntersectionAux<1,2> is(eleA_idx, eleB_idx);
933 
934  ComputeIntersection<1,2> CI(eleA, eleB, mesh);
935  unsigned int n_local_intersection = CI.compute_final(is.points());
936 
937  if(n_local_intersection > 0)
938  {
939  storage.push_back(IntersectionLocal<1,2>(is));
940  intersection_map[eleA_idx].push_back(std::make_pair(
941  eleB_idx,
942  &(storage.back())
943  ));
944  intersection_map[eleB_idx].push_back(std::make_pair(
945  eleA_idx,
946  &(storage.back())
947  ));
948 
949  //DebugOut().fmt("1D-2D intersection [{} - {}]:\n",is.component_ele_idx(), is.bulk_ele_idx());
950  }
951  }
952  }
953  }
954  }
955 
956  MessageOut() << "1D-2D [3]: number of intersections = " << storage.size() << "\n";
957  // just dbg output
958 // for(IntersectionLocal<1,2> &is : storage)
959 // {
960 // DebugOut().fmt("1D-2D intersection [{} - {}]:\n",is.component_ele_idx(), is.bulk_ele_idx());
961 // for(const IntersectionPoint<1,2>& ip : is.points()) {
962 // DebugOut() << ip;
963 // auto p = ip.coords(mesh->element(is.component_ele_idx()));
964 // DebugOut() << "[" << p[0] << " " << p[1] << " " << p[2] << "]\n";
965 // }
966 // }
967 }
968 
969 // void InspectElementsAlgorithm12::compute_single_intersection(const ElementAccessor<3>& eleA,
970 // const ElementAccessor<3>& eleB)
971 // {
972 // ASSERT_DBG(eleA.dim() == 1);
973 // ASSERT_DBG(eleB.dim() == 2);
974 //
975 // this->update_simplex(eleA, simplexA);
976 // this->update_simplex(eleB, simplexB);
977 //
978 // IntersectionAux<1,2> is(eleA.index(), eleB.index(), 0);
979 // // std::vector<unsigned int> prolongation_table;
980 //
981 // ComputeIntersection< Simplex<1>, Simplex<2>> CI(simplexA, simplexB);
982 // unsigned int n_local_intersection = CI.compute_final(is.points());
983 //
984 // if(n_local_intersection > 0)
985 // {
986 // DebugOut().fmt("found: {}\n",n_local_intersection);
987 // intersectionaux_storage12_.push_back(is);
988 // }
989 // }
990 
991 
992 
994 {
995  //DebugOut() << "Intersections 1d-2d (2-bihtree)\n";
997  START_TIMER("Element iteration");
998 
999  for (auto elm : mesh->elements_range()) {
1000  unsigned int component_ele_idx = elm.idx();
1001 
1002  if (elm.dim() == 1) // is component element
1003  //&& elements_bb[component_ele_idx].intersect(mesh_3D_bb)) // its bounding box intersects 3D mesh bounding box
1004  {
1005  std::vector<unsigned int> searchedElements;
1006 
1007  START_TIMER("BIHtree find");
1008  bih.find_bounding_box(bih.ele_bounding_box(component_ele_idx), searchedElements);
1009  END_TIMER("BIHtree find");
1010 
1011  START_TIMER("Bounding box element iteration");
1012 
1013  // Go through all element which bounding box intersects the component element bounding box
1014  for (std::vector<unsigned int>::iterator it = searchedElements.begin(); it!=searchedElements.end(); it++)
1015  {
1016  unsigned int bulk_ele_idx = *it;
1017  ElementAccessor<3> ele_2D = mesh->element_accessor( bulk_ele_idx );
1018 
1019  if (ele_2D.dim() == 2) {
1020 
1021  IntersectionAux<1,2> is(component_ele_idx, bulk_ele_idx);
1022  START_TIMER("Compute intersection");
1023  ComputeIntersection<1,2> CI(elm, ele_2D, mesh);
1024  CI.compute_final(is.points());
1025  END_TIMER("Compute intersection");
1026 
1027  if(is.points().size() > 0) {
1028  intersectionaux_storage12_.push_back(is);
1029  }
1030  }
1031  }
1032  END_TIMER("Bounding box element iteration");
1033  }
1034  }
1035 
1036  END_TIMER("Element iteration");
1037 }
1038 
1040 {
1041  //DebugOut() << "Intersections 1d-2d (2-bihtree) in 2D plane.\n";
1043  START_TIMER("Element iteration");
1044 
1045  for (auto elm : mesh->elements_range()) {
1046  unsigned int component_ele_idx = elm.idx();
1047 
1048  if (elm->dim() == 1) // is component element
1049  //&& elements_bb[component_ele_idx].intersect(mesh_3D_bb)) // its bounding box intersects 3D mesh bounding box
1050  {
1051  std::vector<unsigned int> candidate_list;
1052  bih.find_bounding_box(bih.ele_bounding_box(component_ele_idx), candidate_list);
1053 
1054  START_TIMER("Bounding box element iteration");
1055 
1056  // Go through all element which bounding box intersects the component element bounding box
1057  for(unsigned int bulk_ele_idx : candidate_list) {
1058  ElementAccessor<3> ele_2D = mesh->element_accessor(bulk_ele_idx);
1059 
1060  if (ele_2D->dim() == 2) {
1061 
1062  IntersectionAux<1,2> is(component_ele_idx, bulk_ele_idx);
1063  START_TIMER("Compute intersection");
1064  ComputeIntersection<1,2> CI(elm, ele_2D, mesh);
1065  CI.compute_final_in_plane(is.points());
1066  END_TIMER("Compute intersection");
1067 
1068  if(is.points().size() > 0) {
1069  intersectionaux_storage12_.push_back(is);
1070  }
1071  }
1072  }
1073  END_TIMER("Bounding box element iteration");
1074  }
1075  }
1076  END_TIMER("Element iteration");
1077  MessageOut() << "1D-2D [1]: number of intersections = " << intersectionaux_storage12_.size() << "\n";
1078 }
1079 
1080 // Declaration of specializations implemented in cpp:
1081 template class IntersectionAlgorithmBase<1,3>;
1082 template class IntersectionAlgorithmBase<2,3>;
1083 template class IntersectionAlgorithmBase<1,2>;
1084 template class IntersectionAlgorithmBase<2,2>;
1085 
1086 template class InspectElementsAlgorithm<1>;
1087 template class InspectElementsAlgorithm<2>;
1088 
1089 
const Edge * edge() const
Definition: side_impl.hh:66
unsigned int uint
Classes with algorithms for computation of intersections of meshes.
std::vector< BoundingBox > elements_bb
Elements bounding boxes.
std::vector< unsigned int > get_element_neighbors(const ElementAccessor< 3 > &ele, unsigned int ip_dim, unsigned int ip_obj_idx)
void prolongate(const Prolongation &pr)
Computes the intersection for a candidate in a queue and calls prolongation_decide again...
#define MessageOut()
Macro defining &#39;message&#39; record of log.
Definition: logger.hh:243
Mesh * mesh
Auxiliary function that translates ElementAccessor<3> to Simplex<simplex_dim>.
Fundamental simplicial intersections.
bool intersection_exists(unsigned int component_ele_idx, unsigned int bulk_ele_idx)
A hard way to find whether the intersection of two elements has already been computed, or not.
std::vector< unsigned int > last_slave_for_3D_elements
void expand(const Point &point)
Definition: mesh.h:76
Auxiliary structure for prolongation process.
unsigned int size() const
Returns number of intersection points.
void prolongation_decide(const ElementAccessor< 3 > &comp_ele, const ElementAccessor< 3 > &bulk_ele, IntersectionAux< dim, 3 > &is)
int n_sides
Definition: edges.h:36
Definition: edges.h:26
SideIter side(const unsigned int loc_index)
Definition: accessors.hh:137
bool intersect(const BoundingBox &b2) const
unsigned int n_intersections_
Counter for intersection among elements.
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library...
Definition: asserts.hh:346
unsigned int compute_final_in_plane(std::vector< IPAux12 > &IP12s)
#define ADD_CALLS(n_calls)
Increase number of calls in actual timer.
std::queue< Prolongation > component_queue_
Prolongation queue in the component mesh.
ElementAccessor< 3 > element() const
Definition: side_impl.hh:53
void compute_intersections_BIHtree(const BIHTree &bih)
Uses only BIHtree to find intersection candidates. (No prolongation).
std::vector< unsigned int > component_idx_
virtual ElementAccessor< 3 > element_accessor(unsigned int idx) const
Create and return ElementAccessor to element of given idx.
Definition: mesh.cc:727
unsigned int dim() const
Definition: elements.h:124
void assert_same_intersection(unsigned int comp_ele_idx, unsigned int bulk_ele_idx)
void compute_intersections(std::vector< std::vector< ILpair >> &intersection_map, std::vector< IntersectionLocal< 2, 2 >> &storage)
Runs the core algorithm for computing 2D-2D intersection in 3D.
void set_duplicities(unsigned int n_duplicities)
std::vector< IntersectionAux< 1, 2 > > intersectionaux_storage12_
Stores temporarily 1D-2D intersections.
void find_bounding_box(const BoundingBox &boundingBox, std::vector< unsigned int > &result_list, bool full_list=false) const
Definition: bih_tree.cc:234
void compute_intersections_3(std::vector< std::vector< ILpair >> &intersection_map, std::vector< IntersectionLocal< 1, 2 >> &storage)
Runs the algorithm (3): compute 1D-2D intersection inside 3D mesh. It directly fills the intersection...
void compute_single_intersection(const ElementAccessor< 3 > &eleA, const ElementAccessor< 3 > &eleB, std::vector< IntersectionLocal< 2, 2 >> &storage)
Computes fundamental intersection of two 2D elements.
std::vector< std::vector< IntersectionAux< dim, 3 > > > intersection_list_
Resulting vector of intersections.
void create_component_numbering()
Creates numbering of the 2D components and fills component_idx_ vector.
unsigned int index() const
Definition: accessors.hh:121
unsigned int edge_idx(unsigned int edg_idx) const
Return edge_idx of given index.
Definition: elements.h:139
const BoundingBox & tree_box() const
Definition: bih_tree.cc:229
unsigned int n_sides() const
Definition: elements.h:135
static const IdxVector< (InDim >OutDim?InDim+1:dim-InDim) > interact(TInteraction< OutDim, InDim > interaction)
#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:38
std::vector< IntersectionPointAux< dimA, dimB > > & points()
Returns intersection points by a reference.
Class for 1D-2D intersections.
virtual Range< ElementAccessor< 3 > > elements_range() const
Returns range of bulk elements.
Definition: mesh.cc:1041
void swap(nlohmann::json &j1, nlohmann::json &j2) noexcept(is_nothrow_move_constructible< nlohmann::json >::value andis_nothrow_move_assignable< nlohmann::json >::value)
exchanges the values of two JSON objects
Definition: json.hpp:8688
Class for 2D-2D intersections.
unsigned int compute_final(std::vector< IPAux12 > &IP12s)
virtual unsigned int n_elements(bool boundary=false) const
Returns count of boundary or bulk elements.
Definition: mesh.h:346
bool compute_initial_CI(const ElementAccessor< 3 > &comp_ele, const ElementAccessor< 3 > &bulk_ele)
Computes the first intersection, from which we then prolongate.
void compute_intersections(const BIHTree &bih)
Uses BIHtree to find the initial candidate of a component and then prolongates the component interset...
std::vector< Edge > edges
Vector of MH edges, this should not be part of the geometrical mesh.
Definition: mesh.h:252
#define ASSERT_DBG(expr)
Definition: asserts.hh:349
BoundingBox mesh_3D_bb
Bounding box of all 3D elements.
#define END_TIMER(tag)
Ends a timer with specified tag.
void compute_bounding_boxes()
Computes bounding boxes of all elements. Fills elements_bb and mesh_3D_bb.
double tetrahedron_jacobian() const
Definition: accessors.hh:179
Internal auxiliary class represents an intersection point of simplex<N> and simplex<M>.
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc...
void compute_intersections_2(const BIHTree &bih)
Runs the algorithm (2): compute 1D-2D intersection in 3D ambient space BIH is used to find intersecti...
void compute_intersections_1(const BIHTree &bih)
Runs the algorithm (1): compute 1D-2D intersection in 2D plane. BIH is used to find intersection cand...
unsigned int dim() const
Definition: accessors.hh:87
Internal auxiliary class representing intersection object of simplex<dimA> and simplex<dimB>.
unsigned int create_prolongation(unsigned int bulk_ele_idx, unsigned int component_ele_idx, std::queue< Prolongation > &queue)
unsigned int n_neighs_vb() const
Return number of neighbours.
Definition: elements.h:70
#define DebugOut()
Macro defining &#39;debug&#39; record of log.
Definition: logger.hh:252
unsigned int idx() const
Return local idx of element in boundary / bulk part of element vector.
Definition: accessors.hh:111
Internal class representing intersection point.
std::queue< Prolongation > bulk_queue_
Prolongation queue in the bulk mesh.
SideIter side(const unsigned int i) const
Definition: edges.h:31
Implementation of range helper class.
const BoundingBox & ele_bounding_box(unsigned int ele_idx) const
Gets bounding box of element of given index ele_index.
Definition: bih_tree.cc:72
unsigned int ips_in_face() const
Returns idx of face when all IPs lie on it; -1 otherwise.
Internal class representing intersection object.