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