Flow123d  3.9.0-127d9ed9b
mixed_mesh_intersections.cc
Go to the documentation of this file.
1 /*
2  * MixedMeshIntersections.cpp
3  *
4  * Created on: 13.4.2014
5  * Author: viktor, pe, jb
6  */
7 
8 #include "inspect_elements_algorithm.hh"
9 #include "intersection_point_aux.hh"
10 #include "intersection_aux.hh"
11 #include "intersection_local.hh"
12 
13 #include "system/global_defs.h"
14 #include "system/sys_profiler.hh"
15 
16 #include "mesh/mesh.h"
17 #include "mesh/ref_element.hh"
18 #include "mixed_mesh_intersections.hh"
19 #include "mesh/bih_tree.hh"
20 #include "mesh/accessors.hh"
21 #include "mesh/node_accessor.hh"
22 #include "mesh/range_wrapper.hh"
23 
24 
25 MixedMeshIntersections::MixedMeshIntersections(Mesh* mesh)
26 : mesh(mesh), algorithm13_(mesh), algorithm23_(mesh), algorithm22_(mesh), algorithm12_(mesh)
27 {}
28 
29 MixedMeshIntersections::~MixedMeshIntersections()
30 {}
31 
32 
33 unsigned int MixedMeshIntersections::number_of_components(unsigned int dim)
34 {
35  ASSERT_EQ(dim, 2).error("Not implemented for given dim.\n");
36 // if(dim == 1) return algorithm13_.component_counter_;
37  return algorithm22_.component_counter_; // if(dim == 2) return ...
38 }
39 
40 
41 double MixedMeshIntersections::measure_13()
42 {
43  double subtotal = 0.0;
44 
45  for(unsigned int i = 0; i < intersection_storage13_.size(); i++){
46  ElementAccessor<3> ele = mesh->element_accessor( intersection_storage13_[i].component_ele_idx() );
47  double t1d_length = ele.measure();
48  double local_length = intersection_storage13_[i].compute_measure();
49 
50  if(intersection_storage13_[i].size() == 2)
51  {
52  arma::vec3 from = intersection_storage13_[i][0].coords(ele);
53  arma::vec3 to = intersection_storage13_[i][1].coords(ele);
54 // DebugOut().fmt("sublength from [{} {} {}] to [{} {} {}] = %f\n",
55 // from[0], from[1], from[2],
56 // to[0], to[1], to[2],
57 // local_length*t1d_length);
58  }
59  subtotal += local_length*t1d_length;
60  }
61  return subtotal;
62 }
63 
64 
65 double MixedMeshIntersections::measure_23()
66 {
67  double subtotal = 0.0;
68 
69  for(unsigned int i = 0; i < intersection_storage23_.size(); i++){
70  double t2dArea = mesh->element_accessor( intersection_storage23_[i].component_ele_idx() ).measure();
71  double localArea = intersection_storage23_[i].compute_measure();
72  subtotal += 2*localArea*t2dArea;
73  }
74  return subtotal;
75 }
76 
77 
78 double MixedMeshIntersections::measure_22()
79 {
80  double subtotal = 0.0;
81  double val;
82 
83  for(unsigned int i = 0; i < intersection_storage22_.size(); i++){
84  if(intersection_storage22_[i].size() > 1)
85  {
86  ElementAccessor<3> eleA = mesh->element_accessor( intersection_storage22_[i].component_ele_idx() );
87 // ElementAccessor<3> eleB = mesh->element_accessor( intersection_storage22_[i].bulk_ele_idx() );
88 
89  arma::vec3 from = intersection_storage22_[i][0].coords(eleA);
90  arma::vec3 to = intersection_storage22_[i][1].coords(eleA);
91  val = arma::norm(from - to, 2);
92 
93 // DebugOut().fmt("{}--{}:sublength from [{} {} {}] to [{} {} {}] = {}\n",
94 // eleA.idx(), eleB.idx(),
95 // from[0], from[1], from[2],
96 // to[0], to[1], to[2],
97 // val);
98  subtotal += val;
99  }
100  }
101  return subtotal;
102 }
103 
104 
105 template<uint dim_A, uint dim_B>
106 void MixedMeshIntersections::store_intersection(std::vector<IntersectionLocal<dim_A, dim_B>> &storage, IntersectionAux<dim_A, dim_B> &isec_aux) {
107  //unsigned int ele_a_idx = isec_aux.component_ele_idx();
108  //unsigned int ele_b_idx = isec_aux.bulk_ele_idx();
109 
110  //WARNING:
111  // - not all algorithms uses this function (e.g. 2d2d pushes directly into storage)
112  // - we cannot throw away isec of zero measure generaly (e.g. 1d2d, 1d3d)
113  // - is it not better to test number of IPs according to dimensions?
114 
115  IntersectionLocal<dim_A, dim_B> isec(isec_aux);
116  //if ( (! (dim_A==1 && dim_B==2)) && isec.compute_measure() < 1e-14) return;
117 
118  storage.push_back(isec);
119  /*
120  element_intersections_[ele_a_idx].push_back(
121  std::make_pair(ele_b_idx, &(storage.back())) );
122 
123  element_intersections_[ele_b_idx].push_back(
124  std::make_pair(ele_a_idx, &(storage.back())) );
125  */
126 
127 }
128 
129 template<uint dim_A, uint dim_B>
130 void MixedMeshIntersections::append_to_index( std::vector<IntersectionLocal<dim_A, dim_B>> &storage)
131 {
132  for(auto &isec : storage) {
133 
134  unsigned int ele_a_idx = isec.component_ele_idx();
135  unsigned int ele_b_idx = isec.bulk_ele_idx();
136  ASSERT_EQ(mesh->element_accessor(ele_a_idx)->dim(), dim_A)(ele_a_idx);
137  ASSERT_EQ(mesh->element_accessor(ele_b_idx)->dim(), dim_B)(ele_b_idx);
138  element_intersections_[ele_a_idx].push_back(
139  std::make_pair(ele_b_idx, &(isec)) );
140 
141 
142  if (dim_B==3) {
143  // necessary for 2d-2d intersections
144  element_intersections_[ele_b_idx].push_back(
145  std::make_pair(ele_a_idx, &(isec)) );
146 
147  }
148 /*
149  element_intersections_[ele_b_idx].push_back(
150  std::make_pair(ele_a_idx, &(isec)) );*/
151  }
152 }
153 
154 
155 template<unsigned int dim>
156 void MixedMeshIntersections::compute_intersections(InspectElementsAlgorithm< dim >& iea,
157  std::vector< IntersectionLocal<dim,3>>& storage)
158 {
159  START_TIMER("Intersection algorithm");
160 
161  Mesh::IntersectionSearch is = mesh->get_intersection_search();
162  switch(is){
163  case Mesh::BIHsearch: iea.compute_intersections(mesh->get_bih_tree()); break;
164  case Mesh::BIHonly: iea.compute_intersections_BIHtree(mesh->get_bih_tree()); break;
165  case Mesh::BBsearch: iea.compute_intersections_BB(); break;
166  default: ASSERT_PERMANENT(0).error("Unsupported search algorithm.");
167  }
168 
169  END_TIMER("Intersection algorithm");
170 
171  START_TIMER("Intersection into storage");
172  storage.reserve(iea.n_intersections_);
173 
174  for (auto elm : mesh->elements_range()) {
175  unsigned int idx = elm.idx();
176 
177  if(elm->dim() == dim)
178  {
179 // intersection_map_[idx].resize(iea.intersection_list_[idx].size());
180  element_intersections_[idx].reserve(iea.intersection_list_[idx].size());
181  for(unsigned int j = 0; j < iea.intersection_list_[idx].size(); j++){
182 
183  // skip zero intersections (are made in iea.prolongate())
184  if(iea.intersection_list_[idx][j].size() == 0) continue;
185  store_intersection(storage, iea.intersection_list_[idx][j]); }
186  }
187  }
188  END_TIMER("Intersection into storage");
189 
190 // for(IntersectionLocal<2,3> &is : intersection_storage23_) {
191 // DebugOut().fmt("comp-bulk: {} {}\n", is.component_ele_idx(), is.bulk_ele_idx());
192 // }
193 }
194 
195 void MixedMeshIntersections::compute_intersections_22(vector< IntersectionLocal< 2, 2 > >& storage)
196 {
197  START_TIMER("Intersection algorithm");
198  algorithm22_.compute_intersections(element_intersections_, storage);
199  END_TIMER("Intersection algorithm");
200 
201 // START_TIMER("Intersection into storage");
202 //
203 // storage.reserve(algorithm22_.intersectionaux_storage22_.size());
204 //
205 // for(IntersectionAux<2,2> &is : algorithm22_.intersectionaux_storage22_) {
206 // unsigned int triaA_idx = is.component_ele_idx();
207 // unsigned int triaB_idx = is.bulk_ele_idx();
208 //
209 // //HACK: 'skip flag' move this check into algorithm12_.compute_intersections()
210 // bool skip = false;
211 // for(unsigned int i=0; i<element_intersections_[triaA_idx].size(); i++)
212 // {
213 // if(element_intersections_[triaA_idx][i].first == triaB_idx)
214 // skip = true;
215 // }
216 // if(! skip) {
217 // storage.push_back(IntersectionLocal<2,2>(is));
218 // element_intersections_[triaA_idx].push_back(std::make_pair(
219 // triaB_idx,
220 // &(storage.back())
221 // ));
222 // element_intersections_[triaB_idx].push_back(std::make_pair(
223 // triaA_idx,
224 // &(storage.back())
225 // ));
226 //
227 // // DebugOut().fmt("2D-2D intersection [{} - {}]:\n",
228 // // mesh->element_accessor(is.component_ele_idx()).idx(),
229 // // mesh->element_accessor(is.bulk_ele_idx()).idx());
230 // // for(const IntersectionPointAux<2,2>& ip : is.points()) {
231 // // DebugOut() << ip;
232 // // auto p = ip.coords(mesh->element(is.component_ele_idx()));
233 // // DebugOut() << "[" << p[0] << " " << p[1] << " " << p[2] << "]\n";
234 // // }
235 // }
236 // }
237 // DBGVAR(algorithm22_.intersectionaux_storage22_.size());
238 //
239 // END_TIMER("Intersection into storage");
240 }
241 
242 void MixedMeshIntersections::compute_intersections_12_3(vector< IntersectionLocal< 1, 2 > >& storage)
243 {
244  storage.reserve(intersection_storage13_.size());
245  algorithm12_.compute_intersections_3(element_intersections_, storage);
246  storage.shrink_to_fit();
247 
248 // START_TIMER("Intersection into storage");
249 // storage.reserve(algorithm12_.intersectionaux_storage12_.size());
250 
251 // for(IntersectionAux<1,2> &is : algorithm12_.intersectionaux_storage12_) {
252 // unsigned int abscissa_idx = is.component_ele_idx();
253 // unsigned int triangle_idx = is.bulk_ele_idx();
254 //
255 // //HACK: 'skip flag' move this check into algorithm12_.compute_intersections()
256 // bool skip = false;
257 // for(unsigned int i=0; i<intersection_map_[abscissa_idx].size(); i++)
258 // {
259 // if(intersection_map_[abscissa_idx][i].first == triangle_idx)
260 // skip = true;
261 // }
262 // if(! skip) {
263 // storage.push_back(IntersectionLocal<1,2>(is));
264 // intersection_map_[abscissa_idx].push_back(std::make_pair(
265 // triangle_idx,
266 // &(storage.back())
267 // ));
268 // intersection_map_[triangle_idx].push_back(std::make_pair(
269 // abscissa_idx,
270 // &(storage.back())
271 // ));
272 // DebugOut().fmt("1D-2D intersection [{} - {}]:\n",is.component_ele_idx(), is.bulk_ele_idx());
273 // for(const IntersectionPointAux<1,2>& ip : is.points()) {
274 // //DebugOut() << ip;
275 // auto p = ip.coords(mesh->element(is.component_ele_idx()));
276 // DebugOut() << "[" << p[0] << " " << p[1] << " " << p[2] << "]\n";
277 // }
278 // }
279 // }
280 // END_TIMER("Intersection into storage");
281 }
282 
283 void MixedMeshIntersections::compute_intersections_12_2(vector< IntersectionLocal< 1, 2 > >& storage)
284 {
285  algorithm12_.compute_intersections_2(mesh->get_bih_tree());
286 // DBGVAR(algorithm12_.intersectionaux_storage12_.size());
287 
288  START_TIMER("Intersection into storage");
289  storage.reserve(algorithm12_.intersectionaux_storage12_.size());
290 
291  for(IntersectionAux<1,2> &is : algorithm12_.intersectionaux_storage12_) {
292  store_intersection(storage, is);
293 // DebugOut().fmt("1D-2D intersection [{} - {}]:\n",is.component_ele_idx(), is.bulk_ele_idx());
294 // for(const IntersectionPointAux<1,2>& ip : is.points()) {
295 // //DebugOut() << ip;
296 // auto p = ip.coords(mesh->element(is.component_ele_idx()));
297 // DebugOut() << "[" << p[0] << " " << p[1] << " " << p[2] << "]\n";
298 // }
299  }
300  END_TIMER("Intersection into storage");
301 }
302 
303 void MixedMeshIntersections::compute_intersections_12_1(vector< IntersectionLocal< 1, 2 > >& storage)
304 {
305  algorithm12_.compute_intersections_1(mesh->get_bih_tree());
306 // DBGVAR(algorithm12_.intersectionaux_storage12_.size());
307 
308  START_TIMER("Intersection into storage");
309  storage.reserve(algorithm12_.intersectionaux_storage12_.size());
310 
311  for(IntersectionAux<1,2> &is : algorithm12_.intersectionaux_storage12_) {
312  store_intersection(storage, is);
313 // DebugOut().fmt("1D-2D intersection [{} - {}]:\n",is.component_ele_idx(), is.bulk_ele_idx());
314 // for(const IntersectionPointAux<1,2>& ip : is.points()) {
315 // //DebugOut() << ip;
316 // auto p = ip.coords(mesh->element(is.component_ele_idx()));
317 // DebugOut() << "[" << p[0] << " " << p[1] << " " << p[2] << "]\n";
318 // }
319  }
320  END_TIMER("Intersection into storage");
321 }
322 
323 void MixedMeshIntersections::compute_intersections(IntersectionType d)
324 {
325  element_intersections_.resize(mesh->n_elements());
326 
327  // check whether the mesh is in plane only
328  bool mesh_in_2d_only = false;
329  auto bb = mesh->get_bih_tree().tree_box();
330  for(uint axis = 0; axis < bb.dimension; axis++)
331  if(bb.size(axis) < geometry_epsilon) mesh_in_2d_only = true;
332 
333  if(d & (IntersectionType::d13 | IntersectionType::d12_3)){
334  START_TIMER("Intersections 1D-3D");
335 // DebugOut() << "Intersection Algorithm d13\n";
336  compute_intersections<1>(algorithm13_,intersection_storage13_);
337  END_TIMER("Intersections 1D-3D");
338  }
339  append_to_index(intersection_storage13_);
340 
341 
342  if(d & (IntersectionType::d23 | IntersectionType::d22 | IntersectionType::d12_3)){
343  START_TIMER("Intersections 2D-3D");
344 // DebugOut() << "Intersection Algorithm d23\n";
345  compute_intersections<2>(algorithm23_,intersection_storage23_);
346  END_TIMER("Intersections 2D-3D");
347  }
348  append_to_index(intersection_storage23_);
349 
350 
351  if(d & IntersectionType::d22){
352  START_TIMER("Intersections 2D-2D");
353 // DebugOut() << "Intersection Algorithm d22\n";
354  compute_intersections_22(intersection_storage22_);
355  END_TIMER("Intersections 2D-2D");
356  }
357 
358  if( mesh_in_2d_only){
359  START_TIMER("Intersections 1D-2D (1)");
360  if(d & IntersectionType::d12_1) compute_intersections_12_1(intersection_storage12_);
361  END_TIMER("Intersections 1D-2D (1)");
362  }
363  // make sence only if some intersections in 3D are computed
364  // TODO: this does NOT compute 1d-2d outside 3d bulk
365  // NOTE: create input record in mesh to decide, whether compute also outside (means to call alg. 2)
366  else if( ! intersection_storage13_.empty() &&
367  ! intersection_storage23_.empty() &&
368  (d & IntersectionType::d12_3)){
369  START_TIMER("Intersections 1D-2D (3)");
370  DebugOut() << "Intersection Algorithm d12_3\n";
371  compute_intersections_12_3(intersection_storage12_);
372  END_TIMER("Intersections 1D-2D (3)");
373  }
374  // otherwise compute 1d-2d in the most general case
375  else if(d & IntersectionType::d12_2){
376  START_TIMER("Intersections 1D-2D (2)");
377  DebugOut() << "Intersection Algorithm d12_2\n";
378  compute_intersections_12_2(intersection_storage12_);
379  END_TIMER("Intersections 1D-2D (2)");
380  }
381 
382  //ASSERT_PERMANENT_EQ(intersection_storage13_.size(), 0);
383  //ASSERT_PERMANENT_EQ(intersection_storage23_.size(), 0);
384  //ASSERT_PERMANENT_EQ(intersection_storage22_.size(), 0);
385  // compose master
386  append_to_index(intersection_storage12_);
387 
388  // release temporary links from 3d elements
389  for (auto elm : mesh->elements_range()) {
390  if(elm->dim() == 3) element_intersections_[elm.idx()].clear();
391  }
392 
393 
394 
395 }
396 
397 
398 
399 
400 void MixedMeshIntersections::print_mesh_to_file_13(string name)
401 {
402  string t_name = name;
403 
404  unsigned int number_of_intersection_points = 0;
405  unsigned int number_of_nodes = mesh->n_nodes();
406 
407  for(unsigned int j = 0; j < intersection_storage13_.size();j++){
408  number_of_intersection_points += intersection_storage13_[j].size();
409  }
410 
411  FILE * file;
412  file = fopen((t_name.append(".msh")).c_str(),"w");
413 
414  fprintf(file, "$MeshFormat\n");
415  fprintf(file, "2.2 0 8\n");
416  fprintf(file, "$EndMeshFormat\n");
417  fprintf(file, "$Nodes\n");
418  fprintf(file, "%d\n", (mesh->n_nodes() + number_of_intersection_points));
419 
420  unsigned int idx = 1;
421  for (auto nod : mesh->node_range()) {
422  arma::vec3 _nod = *nod;
423  fprintf(file,"%d %.16f %.16f %.16f\n", idx, _nod[0], _nod[1], _nod[2]);
424  idx++;
425  }
426 
427  for(unsigned int j = 0; j < intersection_storage13_.size();j++){
428  IntersectionLocal<1,3> il = intersection_storage13_[j];
429  ElementAccessor<3> el1D = mesh->element_accessor( il.component_ele_idx() );
430 // ElementAccessor<3> el3D = mesh->element_accessor( il.bulk_ele_idx() );
431 
432  for(unsigned int k = 0; k < il.size();k++){
433  number_of_nodes++;
434  IntersectionPoint<1,3> IP13 = il[k];
435  arma::vec3 global = IP13.coords(el1D);
436 
437 // if(i == 0){
438 // _global = (IP13.local_bcoords_A())[0] * el1D.node(0)->point()
439 // +(IP13.local_bcoords_A())[1] * el1D.node(1)->point();
440 // }else{
441 // _global = (IP13.local_bcoords_B())[0] * el3D.node(0)->point()
442 // +(IP13.local_bcoords_B())[1] * el3D.node(1)->point()
443 // +(IP13.local_bcoords_B())[2] * el3D.node(2)->point()
444 // +(IP13.local_bcoords_B())[3] * el3D.node(3)->point();
445 // }
446 
447  fprintf(file,"%d %.16f %.16f %.16f\n", number_of_nodes, global[0], global[1], global[2]);
448  }
449  }
450 
451  fprintf(file,"$EndNodes\n");
452  fprintf(file,"$Elements\n");
453  fprintf(file,"%d\n", ((unsigned int)intersection_storage13_.size() + mesh->n_elements()) );
454 
455  for (auto elee : mesh->elements_range()) {
456  if(elee->dim() == 3){
457  int id1 = elee.node(0).idx() + 1;
458  int id2 = elee.node(1).idx() + 1;
459  int id3 = elee.node(2).idx() + 1;
460  int id4 = elee.node(3).idx() + 1;
461 
462  fprintf(file,"%d 4 2 %d %d %d %d %d %d\n", elee.idx(), elee.region().id(), elee->pid(), id1, id2, id3, id4);
463  }else if(elee->dim() == 2){
464  int id1 = elee.node(0).idx() + 1;
465  int id2 = elee.node(1).idx() + 1;
466  int id3 = elee.node(2).idx() + 1;
467  fprintf(file,"%d 2 2 %d %d %d %d %d\n", elee.idx(), elee.region().id(), elee->pid(), id1, id2, id3);
468 
469  }else if(elee->dim() == 1){
470  int id1 = elee.node(0).idx() + 1;
471  int id2 = elee.node(1).idx() + 1;
472  fprintf(file,"%d 1 2 %d %d %d %d\n",elee.idx(), elee.region().id(), elee->pid(), id1, id2);
473  }
474  }
475 
476  unsigned int number_of_elements = mesh->n_elements();
477  unsigned int nodes = mesh->n_nodes();
478 
479  for(unsigned int j = 0; j < intersection_storage13_.size();j++){
480  IntersectionLocal<1,3> il = intersection_storage13_[j];
481  number_of_elements++;
482  nodes++;
483  if(il.size() == 1){
484  fprintf(file,"%d 1 2 1001 0 %d %d\n", number_of_elements, nodes, nodes);
485  }else if(il.size() == 2){
486  fprintf(file,"%d 1 2 1001 0 %d %d\n", number_of_elements, nodes, nodes+1);
487  nodes++;
488  }
489  }
490 
491  fprintf(file,"$EndElements\n");
492  fclose(file);
493 }
494 
495 void MixedMeshIntersections::print_mesh_to_file_23(string name)
496 {
497  //for(unsigned int i = 0; i < 2;i++){
498  string t_name = name;
499 
500  unsigned int number_of_intersection_points = 0;
501  unsigned int number_of_nodes = mesh->n_nodes();
502 
503  for(unsigned int j = 0; j < intersection_storage23_.size();j++){
504  number_of_intersection_points += intersection_storage23_[j].size();
505  }
506 
507  FILE * file;
508  file = fopen((t_name.append(".msh")).c_str(),"w");
509 
510  fprintf(file, "$MeshFormat\n");
511  fprintf(file, "2.2 0 8\n");
512  fprintf(file, "$EndMeshFormat\n");
513  fprintf(file, "$Nodes\n");
514  fprintf(file, "%d\n", (mesh->n_nodes() + number_of_intersection_points));
515 
516  unsigned int idx = 1;
517  for (auto nod : mesh->node_range()) {
518  arma::vec3 _nod = *nod;
519  fprintf(file,"%d %.16f %.16f %.16f\n", idx, _nod[0], _nod[1], _nod[2]);
520  idx++;
521  }
522 
523  for(unsigned int j = 0; j < intersection_storage23_.size();j++){
524 
525  IntersectionLocal<2,3> il = intersection_storage23_[j];
526  ElementAccessor<3> el2D = mesh->element_accessor( il.component_ele_idx() );
527 // ElementAccessor<3> el3D = mesh->element_accessor( il.bulk_ele_idx() );
528 
529  for(unsigned int k = 0; k < intersection_storage23_[j].size();k++){
530 
531  number_of_nodes++;
532  IntersectionPoint<2,3> IP23 = il[k];
533  arma::vec3 global = IP23.coords(el2D);
534 // if(i == 0){
535 // _global = (IP23.local_bcoords_A())[0] * el2D.node(0)->point()
536 // +(IP23.local_bcoords_A())[1] * el2D.node(1)->point()
537 // +(IP23.local_bcoords_A())[2] * el2D.node(2)->point();
538 // }else{
539 // _global = (IP23.local_bcoords_B())[0] * el3D.node(0)->point()
540 // +(IP23.local_bcoords_B())[1] * el3D.node(1)->point()
541 // +(IP23.local_bcoords_B())[2] * el3D.node(2)->point()
542 // +(IP23.local_bcoords_B())[3] * el3D.node(3)->point();
543 // }
544  fprintf(file,"%d %.16f %.16f %.16f\n", number_of_nodes, global[0], global[1], global[2]);
545  }
546  }
547 
548  fprintf(file,"$EndNodes\n");
549  fprintf(file,"$Elements\n");
550  fprintf(file,"%d\n", (number_of_intersection_points + mesh->n_elements()) );
551 
552  for (auto elee : mesh->elements_range()) {
553  if(elee->dim() == 3){
554  int id1 = elee.node(0).idx() + 1;
555  int id2 = elee.node(1).idx() + 1;
556  int id3 = elee.node(2).idx() + 1;
557  int id4 = elee.node(3).idx() + 1;
558 
559  fprintf(file,"%d 4 2 %d %d %d %d %d %d\n", elee.idx(), elee.region().id(), elee->pid(), id1, id2, id3, id4);
560  }else if(elee->dim() == 2){
561  int id1 = elee.node(0).idx() + 1;
562  int id2 = elee.node(1).idx() + 1;
563  int id3 = elee.node(2).idx() + 1;
564  fprintf(file,"%d 2 2 %d %d %d %d %d\n", elee.idx(), elee.region().id(), elee->pid(), id1, id2, id3);
565 
566  }else{
567  int id1 = elee.node(0).idx() + 1;
568  int id2 = elee.node(1).idx() + 1;
569  fprintf(file,"%d 1 2 %d %d %d %d\n",elee.idx(), elee.region().id(), elee->pid(), id1, id2);
570  }
571  }
572 
573  unsigned int number_of_elements = mesh->n_elements();
574  unsigned int last = 0;
575  unsigned int nodes = mesh->n_nodes();
576 
577  for(unsigned int j = 0; j < intersection_storage23_.size();j++){
578  IntersectionLocal<2,3> il = intersection_storage23_[j];
579 
580  for(unsigned int k = 0; k < il.size();k++){
581 
582  number_of_elements++;
583  nodes++;
584  if(k == 0){
585  last = nodes;
586  }
587 
588  if((k+1) == il.size() && il.size()){
589  fprintf(file,"%d 1 2 1002 0 %d %d\n", number_of_elements, nodes, last);
590  }else{
591  fprintf(file,"%d 1 2 1002 0 %d %d\n", number_of_elements, nodes, nodes+1);
592  }
593  //if(il.size() < 3) break;
594  }
595  }
596 
597  fprintf(file,"$EndElements\n");
598  fclose(file);
599  //}
600 }
601 
602 
MixedMeshIntersections::store_intersection
void store_intersection(std::vector< IntersectionLocal< dim_A, dim_B >> &storage, IntersectionAux< dim_A, dim_B > &isec_aux)
Definition: mixed_mesh_intersections.cc:106
d12_3
@ d12_3
Definition: mixed_mesh_intersections.hh:44
ref_element.hh
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.
MixedMeshIntersections::measure_23
double measure_23()
Definition: mixed_mesh_intersections.cc:65
bih_tree.hh
MixedMeshIntersections::algorithm22_
InspectElementsAlgorithm22 algorithm22_
Definition: mixed_mesh_intersections.hh:111
Element::dim
unsigned int dim() const
Definition: elements.h:120
MixedMeshIntersections::compute_intersections_12_3
void compute_intersections_12_3(std::vector< IntersectionLocal< 1, 2 >> &storage)
Definition: mixed_mesh_intersections.cc:242
d12_1
@ d12_1
Definition: mixed_mesh_intersections.hh:42
d23
@ d23
Definition: mixed_mesh_intersections.hh:40
MixedMeshIntersections::compute_intersections_12_2
void compute_intersections_12_2(std::vector< IntersectionLocal< 1, 2 >> &storage)
Definition: mixed_mesh_intersections.cc:283
InspectElementsAlgorithm::n_intersections_
unsigned int n_intersections_
Counter for intersection among elements.
Definition: inspect_elements_algorithm.hh:134
InspectElementsAlgorithm::compute_intersections_BIHtree
void compute_intersections_BIHtree(const BIHTree &bih)
Uses only BIHtree to find intersection candidates. (No prolongation).
Definition: inspect_elements_algorithm.cc:266
fmt::fprintf
FMT_FUNC int fprintf(std::ostream &os, CStringRef format, ArgList args)
Definition: ostream.cc:56
geometry_epsilon
static const double geometry_epsilon
Definition: intersection_point_aux.hh:28
d13
@ d13
Definition: mixed_mesh_intersections.hh:39
MixedMeshIntersections::algorithm23_
InspectElementsAlgorithm< 2 > algorithm23_
Definition: mixed_mesh_intersections.hh:110
MixedMeshIntersections::compute_intersections_12_1
void compute_intersections_12_1(std::vector< IntersectionLocal< 1, 2 >> &storage)
Definition: mixed_mesh_intersections.cc:303
std::vector
Definition: doxy_dummy_defs.hh:7
ElementAccessor< 3 >
arma::vec3
Definition: doxy_dummy_defs.hh:17
InspectElementsAlgorithm::compute_intersections
void compute_intersections(const BIHTree &bih)
Uses BIHtree to find the initial candidate of a component and then prolongates the component interset...
Definition: inspect_elements_algorithm.cc:142
uint
unsigned int uint
Definition: mh_dofhandler.hh:101
MeshBase::elements_range
Range< ElementAccessor< 3 > > elements_range() const
Returns range of mesh elements.
Definition: mesh.cc:1174
MixedMeshIntersections::intersection_storage13_
std::vector< IntersectionLocal< 1, 3 > > intersection_storage13_
Stores 1D-3D intersections.
Definition: mixed_mesh_intersections.hh:68
Mesh::get_intersection_search
IntersectionSearch get_intersection_search()
Getter for input type selection for intersection search algorithm.
Definition: mesh.cc:275
InspectElementsAlgorithm
Class implements algorithm for dim-dimensional intersections with 3D elements.
Definition: inspect_elements_algorithm.hh:98
MixedMeshIntersections::MixedMeshIntersections
MixedMeshIntersections(Mesh *mesh)
Definition: mixed_mesh_intersections.cc:25
InspectElementsAlgorithm22::component_counter_
unsigned int component_counter_
Definition: inspect_elements_algorithm.hh:208
MixedMeshIntersections::~MixedMeshIntersections
~MixedMeshIntersections()
Definition: mixed_mesh_intersections.cc:29
MixedMeshIntersections::number_of_components
unsigned int number_of_components(unsigned int dim)
Definition: mixed_mesh_intersections.cc:33
IntersectionType
IntersectionType
Selection of intersections of different dimensions.
Definition: mixed_mesh_intersections.hh:36
accessors.hh
Mesh::BBsearch
@ BBsearch
Definition: mesh.h:397
Mesh::get_bih_tree
const BIHTree & get_bih_tree()
Getter for BIH. Creates and compute BIH at first call.
Definition: mesh.cc:1068
inspect_elements_algorithm.hh
sys_profiler.hh
MixedMeshIntersections::algorithm12_
InspectElementsAlgorithm12 algorithm12_
Definition: mixed_mesh_intersections.hh:112
ASSERT_PERMANENT
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
MixedMeshIntersections::compute_intersections_22
void compute_intersections_22(std::vector< IntersectionLocal< 2, 2 >> &storage)
Definition: mixed_mesh_intersections.cc:195
BIHTree::tree_box
const BoundingBox & tree_box() const
Definition: bih_tree.cc:229
MixedMeshIntersections::mesh
Mesh * mesh
Mesh pointer.
Definition: mixed_mesh_intersections.hh:107
InspectElementsAlgorithm22::compute_intersections
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.
Definition: inspect_elements_algorithm.cc:708
mixed_mesh_intersections.hh
MeshBase::node_range
Range< NodeAccessor< 3 > > node_range() const
Returns range of nodes.
Definition: mesh.cc:1180
MixedMeshIntersections::compute_intersections
void compute_intersections(IntersectionType d=IntersectionType::all)
Definition: mixed_mesh_intersections.cc:323
MixedMeshIntersections::print_mesh_to_file_13
void print_mesh_to_file_13(std::string name)
Generates a mesh file of the given name, including the intersection.
Definition: mixed_mesh_intersections.cc:400
MeshBase::element_accessor
ElementAccessor< 3 > element_accessor(unsigned int idx) const
Create and return ElementAccessor to element of given idx.
Definition: mesh.cc:866
IntersectionPoint::coords
arma::vec3 coords(ElementAccessor< 3 > comp_ele) const
Definition: intersection_local.cc:107
Mesh::IntersectionSearch
IntersectionSearch
Types of search algorithm for finding intersection candidates.
Definition: mesh.h:394
ASSERT_EQ
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
MixedMeshIntersections::intersection_storage22_
std::vector< IntersectionLocal< 2, 2 > > intersection_storage22_
Stores 2D-2D intersections.
Definition: mixed_mesh_intersections.hh:72
InspectElementsAlgorithm12::compute_intersections_1
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...
Definition: inspect_elements_algorithm.cc:1035
mesh.h
MixedMeshIntersections::algorithm13_
InspectElementsAlgorithm< 1 > algorithm13_
Definition: mixed_mesh_intersections.hh:109
MixedMeshIntersections::element_intersections_
std::vector< std::vector< ILpair > > element_intersections_
Definition: mixed_mesh_intersections.hh:83
intersection_aux.hh
Internal class representing intersection object.
d22
@ d22
Definition: mixed_mesh_intersections.hh:41
Mesh
Definition: mesh.h:361
MixedMeshIntersections::intersection_storage23_
std::vector< IntersectionLocal< 2, 3 > > intersection_storage23_
Stores 2D-3D intersections.
Definition: mixed_mesh_intersections.hh:70
MixedMeshIntersections::intersection_storage12_
std::vector< IntersectionLocal< 1, 2 > > intersection_storage12_
Stores 1D-2D intersections.
Definition: mixed_mesh_intersections.hh:74
Mesh::BIHsearch
@ BIHsearch
Definition: mesh.h:395
MixedMeshIntersections::print_mesh_to_file_23
void print_mesh_to_file_23(std::string name)
Definition: mixed_mesh_intersections.cc:495
InspectElementsAlgorithm::compute_intersections_BB
void compute_intersections_BB()
Definition: inspect_elements_algorithm.cc:321
global_defs.h
Global macros to enhance readability and debugging, general constants.
IntersectionAux
Internal auxiliary class representing intersection object of simplex<dimA> and simplex<dimB>.
Definition: compute_intersection.hh:47
IntersectionLocal
Class represents intersection of two elements.
Definition: inspect_elements_algorithm.hh:34
IntersectionPoint
Class represents an intersection point of simplex<N> and simplex<M>. It contains barycentric coordina...
Definition: intersection_local.hh:32
InspectElementsAlgorithm12::intersectionaux_storage12_
std::vector< IntersectionAux< 1, 2 > > intersectionaux_storage12_
Stores temporarily 1D-2D intersections.
Definition: inspect_elements_algorithm.hh:277
InspectElementsAlgorithm12::compute_intersections_2
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...
Definition: inspect_elements_algorithm.cc:989
IntersectionLocal::size
unsigned int size() const
Returns number of intersection points.
Definition: intersection_local.hh:143
InspectElementsAlgorithm12::compute_intersections_3
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...
Definition: inspect_elements_algorithm.cc:879
MeshBase::n_elements
unsigned int n_elements() const
Definition: mesh.h:111
MixedMeshIntersections::measure_22
double measure_22()
Definition: mixed_mesh_intersections.cc:78
d12_2
@ d12_2
Definition: mixed_mesh_intersections.hh:43
DebugOut
#define DebugOut()
Macro defining 'debug' record of log.
Definition: logger.hh:284
Mesh::BIHonly
@ BIHonly
Definition: mesh.h:396
START_TIMER
#define START_TIMER(tag)
Starts a timer with specified tag.
Definition: sys_profiler.hh:115
MixedMeshIntersections::measure_13
double measure_13()
Computes the size of the intersection in dim dimenstions.
Definition: mixed_mesh_intersections.cc:41
MeshBase::n_nodes
unsigned int n_nodes() const
Definition: mesh.h:171
MixedMeshIntersections::append_to_index
void append_to_index(std::vector< IntersectionLocal< dim_A, dim_B >> &storage)
Definition: mixed_mesh_intersections.cc:130
InspectElementsAlgorithm::intersection_list_
std::vector< std::vector< IntersectionAux< dim, 3 > > > intersection_list_
Resulting vector of intersections.
Definition: inspect_elements_algorithm.hh:151
ElementAccessor::measure
double measure() const
Computes the measure of the element.
Definition: accessors_impl.hh:59
END_TIMER
#define END_TIMER(tag)
Ends a timer with specified tag.
Definition: sys_profiler.hh:149
intersection_point_aux.hh
Internal class representing intersection point.
range_wrapper.hh
Implementation of range helper class.
intersection_local.hh
Classes with algorithms for computation of intersections of meshes.
node_accessor.hh
Generated by   doxygen 1.8.17