Functions
locate_zeta_tester.h File Reference

Go to the source code of this file.

Functions

void check_locate_zeta (Mesh *mesh_pt)
 

Function Documentation

◆ check_locate_zeta()

void check_locate_zeta ( Mesh *  mesh_pt)
29 {
30 
31  // Spatial dimension
32  unsigned n_dim=mesh_pt->finite_element_pt(0)->dim();
33 
34 
35  //Check elemental locate zeta:
36  {
37  FiniteElement* el_pt=mesh_pt->finite_element_pt(0);
38  Vector<double> cog(n_dim);
39  double max_radius=0.0;
40  el_pt->get_centre_of_gravity_and_max_radius_in_terms_of_zeta(cog,max_radius);
41  std::ofstream outfile;
42  outfile.open("RESLT/single_element.dat");
43  el_pt->output(outfile,el_pt->nnode_1d());
44  outfile.close();
45  if (n_dim==2)
46  {
47  outfile.open("RESLT/enclosing_circle_for_single_element.dat");
48  unsigned nplot=100;
49  outfile << "ZONE I=" << nplot << std::endl;
50  for (unsigned i=0;i<nplot;i++)
51  {
52  double phi=2.0*MathematicalConstants::Pi*double(i)/double(nplot-1);
53  outfile << cog[0]+max_radius*cos(phi) << " "
54  << cog[1]+max_radius*sin(phi) << std::endl;
55  }
56  outfile.close();
57  }
58  else if (n_dim==3)
59  {
60  outfile.open("RESLT/enclosing_sphere_for_single_element.dat");
61  unsigned nplot=100;
62  outfile << "ZONE I=" << nplot << ", J=" << nplot << std::endl;
63  for (unsigned i=0;i<nplot;i++)
64  {
65  double phi=2.0*MathematicalConstants::Pi*double(i)/double(nplot-1);
66  for (unsigned j=0;j<nplot;j++)
67  {
68  double theta=MathematicalConstants::Pi*(-0.5+double(j)/double(nplot-1));
69  outfile << cog[0]+max_radius*cos(phi)*cos(theta) << " "
70  << cog[1]+max_radius*sin(phi)*cos(theta) << " "
71  << cog[2]+max_radius*sin(theta) << " "
72  << std::endl;
73  }
74  }
75  outfile.close();
76  }
77  else
78  {
79  oomph_info << "wtf?" << std::endl;
80  abort();
81  }
82 
83  // Check elemental locate zeta
84  Vector<double> zeta(cog);
85  GeomObject* geom_object_pt=0;
86  Vector<double> s(n_dim);
87  bool use_coordinate_as_initial_guess=false;
88  el_pt->locate_zeta(zeta,
89  geom_object_pt,
90  s,
91  use_coordinate_as_initial_guess);
92 
93  if (geom_object_pt==0)
94  {
95  oomph_info << "elemental locate zeta failed\n";
96  abort();
97  }
98  else
99  {
100  oomph_info << "elemental locate zeta succeeded\n";
101 
102  }
103  }
104 
105 
106  // Create Mesh as geometric object
107  RefineableBinArrayParameters ref_bin_array_parameters(mesh_pt);
108  NonRefineableBinArrayParameters non_ref_bin_array_parameters(mesh_pt);
109 #ifdef OOMPH_HAS_CGAL
110  CGALSamplePointContainerParameters sample_point_container_parameters(mesh_pt);
111 #endif
112 
113 
114  // Loop over the various sample point containers
115  unsigned n_sample_point_container_type=2;
116 #ifdef OOMPH_HAS_CGAL
117  n_sample_point_container_type++;
118 #endif
119  for (unsigned ref_flag=0;ref_flag<n_sample_point_container_type;ref_flag++)
120  {
121 
122  SamplePointContainerParameters* sample_point_container_params_pt=0;
123  if (ref_flag==0)
124  {
125  sample_point_container_params_pt=&non_ref_bin_array_parameters;
126  oomph_info << "\n\n\nDOING NON-REFINEABLE BIN ARRAY" << std::endl;
127  oomph_info << "==============================" << std::endl;
128  }
129  else if (ref_flag==1)
130  {
131  sample_point_container_params_pt=&ref_bin_array_parameters;
132  oomph_info << "\n\n\nDOING REFINEABLE BIN ARRAY" << std::endl;
133  oomph_info << "==========================" << std::endl;
134  }
135 #ifdef OOMPH_HAS_CGAL
136  else if (ref_flag==2)
137  {
138  sample_point_container_params_pt=&sample_point_container_parameters;
139  oomph_info << "\n\n\nDOING CGAL BASED CONTAINER" << std::endl;
140  oomph_info << "==========================" << std::endl;
141  }
142 #endif
143  else
144  {
145  oomph_info << "wrong container\n";
146  abort();
147  }
148 
149  MeshAsGeomObject* mesh_geom_obj_pt=new MeshAsGeomObject(sample_point_container_params_pt);
150 
151 
152 // hierher kill
153 /* if (ref_flag==0) */
154 /* { */
155 /* sample_point_container_params_pt=&non_ref_bin_array_parameters; */
156 /* oomph_info << "\n\n\nDOING NON-REFINEABLE BIN ARRAY" << std::endl; */
157 /* oomph_info << "==============================" << std::endl; */
158 /* } */
159 /* else if (ref_flag==1) */
160 /* { */
161 /* sample_point_container_params_pt=&ref_bin_array_parameters; */
162 /* oomph_info << "\n\n\nDOING REFINEABLE BIN ARRAY" << std::endl; */
163 /* oomph_info << "==========================" << std::endl; */
164 /* } */
165 /* #ifdef OOMPH_HAS_CGAL */
166 /* else if (ref_flag==2) */
167 /* { */
168 /* sample_point_container_params_pt=&sample_point_container_parameters; */
169 /* oomph_info << "\n\n\nDOING CGAL BASED CONTAINER" << std::endl; */
170 /* oomph_info << "==========================" << std::endl; */
171 /* } */
172 /* #endif */
173 /* else */
174 /* { */
175 /* oomph_info << "wrong container\n"; */
176 /* abort(); */
177 /* } */
178 
179 
180  // Output bin structure
181  {
182 
183  if (ref_flag==0)
184  {
185  std::string filename="RESLT/non_ref_bin.dat";
186  ofstream outfile;
187  outfile.open(filename.c_str());
188  dynamic_cast<BinArray*>(mesh_geom_obj_pt->sample_point_container_pt())->
189  output_bin_vertices(outfile);
190  outfile.close();
191  }
192  else if (ref_flag==1)
193  {
194  std::string filename="RESLT/ref_bin.dat";
195  ofstream outfile;
196  outfile.open(filename.c_str());
197  dynamic_cast<BinArray*>(mesh_geom_obj_pt->sample_point_container_pt())->
198  output_bin_vertices(outfile);
199  outfile.close();
200  }
201  }
202 
203 
204  // CHECK LOTS OF POINTS INSIDE THE MESH
205  //-------------------------------------
206  {
207  oomph_info << "\n\n\nLots of zetas:" << std::endl;
208  oomph_info << "--------------\n\n" << std::endl;
209 
210 
211  // Do lots of locate zeta calls, check the outcome and doc how
212  // long it took.
213  unsigned count=0;
214  Vector<double> s(n_dim);
215  Vector<double> zeta(n_dim);
216  GeomObject* geom_obj_pt=0;
217  Vector<double> s_from_locate_zeta(n_dim);
218  Vector<double> zeta_test(n_dim);
219  double max_error=0.0;
220  double t_locate_zeta=0.0;
221  unsigned nelem=mesh_pt->nelement();
222  for (unsigned e=0;e<nelem;e++)
223  {
224  FiniteElement* el_pt=mesh_pt->finite_element_pt(e);
225  Integral* int_pt=el_pt->integral_pt();
226  unsigned nint=int_pt->nweight();
227  for (unsigned j=0;j<nint;j++)
228  {
229  // Find integration point in element
230  Vector<double> s(n_dim);
231  for (unsigned i=0;i<n_dim;i++)
232  {
233  s[i]=int_pt->knot(j,i);
234  }
235  el_pt->interpolated_zeta(s,zeta);
236 
237  // Now reverse mapping
238  double start_t = TimingHelpers::timer();
239  mesh_geom_obj_pt->locate_zeta(zeta,geom_obj_pt,s_from_locate_zeta);
240  double end_t = TimingHelpers::timer();
241  t_locate_zeta+=end_t-start_t;
242  count++;
243 
244  if (geom_obj_pt==0)
245  {
246  oomph_info << "Search for zeta = ( "
247  << zeta[0] << " "
248  << zeta[0] << " "
249  << " failed "
250  << std::endl;
251  abort();
252  }
253 
254  // Test
255  geom_obj_pt->interpolated_zeta(s_from_locate_zeta,zeta_test);
256 
257  double error=0.0;
258  for (unsigned i=0;i<n_dim;i++)
259  {
260  error+=pow(zeta[i]-zeta_test[i],2);
261  }
262  error=sqrt(error);
263  if (error>max_error) max_error=error;
264  }
265  }
266  oomph_info << "done/tested " << count << " located zetas. Max. error = "
267  << max_error << std::endl
268  << "Total time for locate zeta: " << t_locate_zeta
269  << " sec " << std::endl;
270 
271  }
272 
273  // NOW DO ONE POINT BUT MAKE LOCATE ZETA FAIL AND DOC SPIRAL
274  //----------------------------------------------------------
275  {
276  oomph_info << "\n\n\nEnforced failure for one zeta:" << std::endl;
277  oomph_info << "------------------------------\n\n" << std::endl;
278 
279  BinArray::Always_fail_elemental_locate_zeta=true;
280  std::stringstream filename;
281  filename << "RESLT/sample_points_" << n_dim << "d_";
282  if (ref_flag==0)
283  {
284  filename << "nonref_bin";
285  }
286  else if (ref_flag==1)
287  {
288  filename << "ref_bin";
289  }
290  else if (ref_flag==2)
291  {
292  filename << "cgal";
293  }
294  else
295  {
296  oomph_info << "wrong container\n";
297  abort();
298  }
299  filename << ".dat";
300  BinArray::Visited_sample_points_file.open(filename.str().c_str());
301  Vector<double> zeta(n_dim,0.5);
302  GeomObject* geom_obj_pt=0;
303  Vector<double> s_from_locate_zeta(n_dim);
304  Vector<double> zeta_test(n_dim,0.5);
305  double t_locate_zeta=0.0;
306  double start_t = TimingHelpers::timer();
307  mesh_geom_obj_pt->locate_zeta(zeta,geom_obj_pt,s_from_locate_zeta);
308  double end_t = TimingHelpers::timer();
309  t_locate_zeta+=end_t-start_t;
310 
311  // Test
312  oomph_info << "It took " << t_locate_zeta << " secs to ";
313  if (geom_obj_pt!=0)
314  {
315  geom_obj_pt->interpolated_zeta(s_from_locate_zeta,zeta_test);
316  oomph_info << "find: "
317  << zeta_test[0] << " "
318  << zeta_test[1] << " "
319  << std::endl;
320 
321  unsigned nvisit=mesh_geom_obj_pt->sample_point_container_pt()->
322  total_number_of_sample_points_visited_during_locate_zeta_from_top_level();
323  unsigned ntotal=mesh_geom_obj_pt->sample_point_container_pt()->
324  total_number_of_sample_points_computed_recursively();
325 
326  oomph_info << "Visited " << nvisit << " sample points "
327  << " out of a total of " << ntotal << std::endl;
328  }
329  else
330  {
331  oomph_info << "not find: "
332  << zeta_test[0] << " "
333  << zeta_test[1] << " "
334  << std::endl;
335 
336  unsigned nvisit=mesh_geom_obj_pt->sample_point_container_pt()->
337  total_number_of_sample_points_visited_during_locate_zeta_from_top_level();
338  unsigned ntotal=mesh_geom_obj_pt->sample_point_container_pt()->
339  total_number_of_sample_points_computed_recursively();
340 
341  oomph_info << "Visited " << nvisit << " sample points "
342  << " out of a total of " << ntotal << std::endl;
343 
344  if (nvisit!=ntotal)
345  {
346  oomph_info << "ERROR\n";
347  abort();
348  }
349 
350  }
351 
352  // Reset and close
353  BinArray::Always_fail_elemental_locate_zeta=false;
354  BinArray::Visited_sample_points_file.close();
355  }
356 
357 
358  // NOW DO ONE POINT BUT MAKE LOCATE ZETA FAIL BECAUSE OF
359  //------------------------------------------------------
360  // LIMITED RADIUS; DOC SPIRAL
361  // --------------------------
362  {
363  oomph_info
364  << "\n\n\nEnforced failure for one zeta because of limited radius:"
365  << std::endl;
366  oomph_info
367  << "--------------------------------------------------------\n\n"
368  << std::endl;
369 
370  BinArray::Always_fail_elemental_locate_zeta=true;
371  std::stringstream filename;
372  filename << "RESLT/sample_points_limited_radius" << n_dim << "d_";
373  if (ref_flag==0)
374  {
375  filename << "nonref_bin";
376  }
377  else if (ref_flag==1)
378  {
379  filename << "ref_bin";
380  }
381  else if (ref_flag==2)
382  {
383  filename << "cgal";
384  }
385  else
386  {
387  oomph_info << "wrong container\n";
388  abort();
389  }
390  filename << ".dat";
391  BinArray::Visited_sample_points_file.open(filename.str().c_str());
392  Vector<double> zeta(n_dim);
393  GeomObject* geom_obj_pt=0;
394  Vector<double> s_from_locate_zeta(n_dim);
395  Vector<double> zeta_test(n_dim);
396  Vector<std::pair<double, double> > min_and_max_coordinates=
397  mesh_geom_obj_pt->sample_point_container_pt()->
398  min_and_max_coordinates();
399  double min_dim=DBL_MAX;
400  for (unsigned i=0;i<n_dim;i++)
401  {
402  zeta[i]=0.5*(min_and_max_coordinates[i].first+
403  min_and_max_coordinates[i].second);
404  double dist=fabs(min_and_max_coordinates[i].first-
405  min_and_max_coordinates[i].second);
406  if (dist<min_dim) min_dim=dist;
407  }
408  mesh_geom_obj_pt->sample_point_container_pt()->max_search_radius()=
409  0.1*min_dim;
410 
411  oomph_info << "Search radius: "
412  << mesh_geom_obj_pt->sample_point_container_pt()->
413  max_search_radius() << std::endl;
414 
415  double t_locate_zeta=0.0;
416  double start_t = TimingHelpers::timer();
417  mesh_geom_obj_pt->locate_zeta(zeta,geom_obj_pt,s_from_locate_zeta);
418  double end_t = TimingHelpers::timer();
419  t_locate_zeta+=end_t-start_t;
420 
421  // Test
422  oomph_info << "It took " << t_locate_zeta << " secs to ";
423  if (geom_obj_pt!=0)
424  {
425  geom_obj_pt->interpolated_zeta(s_from_locate_zeta,zeta_test);
426  oomph_info << "find: "
427  << zeta_test[0] << " "
428  << zeta_test[1] << " "
429  << std::endl;
430 
431  unsigned nvisit=mesh_geom_obj_pt->sample_point_container_pt()->
432  total_number_of_sample_points_visited_during_locate_zeta_from_top_level();
433  unsigned ntotal=mesh_geom_obj_pt->sample_point_container_pt()->
434  total_number_of_sample_points_computed_recursively();
435 
436  oomph_info << "Visited " << nvisit << " sample points "
437  << " out of a total of " << ntotal << std::endl;
438  }
439  else
440  {
441  oomph_info << "not find: "
442  << zeta_test[0] << " "
443  << zeta_test[1] << " "
444  << std::endl;
445 
446  unsigned nvisit=mesh_geom_obj_pt->sample_point_container_pt()->
447  total_number_of_sample_points_visited_during_locate_zeta_from_top_level();
448  unsigned ntotal=mesh_geom_obj_pt->sample_point_container_pt()->
449  total_number_of_sample_points_computed_recursively();
450 
451  oomph_info << "Visited " << nvisit << " sample points "
452  << " out of a total of " << ntotal << std::endl;
453 
454  /* // Here we should only visit fewer than the lot... */
455  /* if (nvisit>=ntotal) */
456  /* { */
457  /* oomph_info << "ERROR\n"; */
458  /* abort(); */
459  /* } */
460 
461  }
462 
463  // Reset and close
464  BinArray::Always_fail_elemental_locate_zeta=false;
465  BinArray::Visited_sample_points_file.close();
466  }
467 
468  // Reset
469  mesh_geom_obj_pt->sample_point_container_pt()->max_search_radius()=
470  DBL_MAX;
471 
472 
473  // NOW DO ONE POINT BUT MAKE LOCATE ZETA FAIL AND DOC INCREMENTAL SPIRAL
474  //---------------------------------------------------------------------
475  {
476  oomph_info << "\n\n\nEnforced failure for one zeta (spiraling):" << std::endl;
477  oomph_info << "------------------------------------------\n\n" << std::endl;
478 
479  Vector<double> zeta(n_dim,0.5);
480  GeomObject* geom_obj_pt=0;
481 
482  // Get number of sample points
483  unsigned nsample_points_in_bin_array=mesh_geom_obj_pt->
484  sample_point_container_pt()->total_number_of_sample_points_computed_recursively();
485 
486  // Maximum spiral level within the cartesian bin structure
487  // (to be assigned below)
488  unsigned n_max_level=0;
489 
490 
491  // Initialise:
492  double t_locate_zeta=0.0;
493 
494  // Refineable
495  if (mesh_geom_obj_pt->sample_point_container_version()==
496  UseRefineableBinArray)
497  {
498  RefineableBinArray* sample_point_container_pt=
499  dynamic_cast<RefineableBinArray*>(mesh_geom_obj_pt->
500  sample_point_container_pt());
501 
502  sample_point_container_pt->
503  last_sample_point_to_actually_lookup_during_locate_zeta() =
504  sample_point_container_pt->
505  initial_last_sample_point_to_actually_lookup_during_locate_zeta();
506  sample_point_container_pt->
507  first_sample_point_to_actually_lookup_during_locate_zeta() = 0;
508  }
509  // Non-refineable
510  else if (mesh_geom_obj_pt->sample_point_container_version()==
511  UseNonRefineableBinArray)
512  {
513  NonRefineableBinArray* sample_point_container_pt=
514  dynamic_cast<NonRefineableBinArray*>(mesh_geom_obj_pt->
515  sample_point_container_pt());
516 
517  // Initialise spiral levels
518  sample_point_container_pt->current_min_spiral_level()=0;
519  sample_point_container_pt->current_max_spiral_level()=
520  sample_point_container_pt->n_spiral_chunk()-1;
521 
522  // Find maximum spiral level within the cartesian bin structure
523  n_max_level=sample_point_container_pt->max_bin_dimension();
524 
525  // Limit it
526  if (sample_point_container_pt->current_max_spiral_level()>n_max_level)
527  {
528  sample_point_container_pt->current_max_spiral_level()=n_max_level-1;
529  }
530  }
531 #ifdef OOMPH_HAS_CGAL
532  // CGAL
533  else if (mesh_geom_obj_pt->sample_point_container_version()==
534  UseCGALSamplePointContainer)
535  {
536  CGALSamplePointContainer* sample_point_container_pt=
537  dynamic_cast<CGALSamplePointContainer*>(mesh_geom_obj_pt->
538  sample_point_container_pt());
539 
540  sample_point_container_pt->
541  last_sample_point_to_actually_lookup_during_locate_zeta() =
542  sample_point_container_pt->
543  initial_last_sample_point_to_actually_lookup_during_locate_zeta();
544  sample_point_container_pt->
545  first_sample_point_to_actually_lookup_during_locate_zeta() = 0;
546  }
547 #endif
548  else
549  {
550  oomph_info << "wrong sample point container\n";
551  abort();
552  }
553 
554  // Loop over "spirals/levels" away from the current position
555  unsigned i_level = 0;
556  bool has_not_reached_max_level_of_search=true;
557  while ((geom_obj_pt==0)&&(has_not_reached_max_level_of_search))
558  {
559  oomph_info << "\n\nSearch level: " << i_level << std::endl;
560 
561  // Refineable
562  if (mesh_geom_obj_pt->sample_point_container_version()==
563  UseRefineableBinArray)
564  {
565 
566  RefineableBinArray* sample_point_container_pt=
567  dynamic_cast<RefineableBinArray*>(mesh_geom_obj_pt->
568  sample_point_container_pt());
569 
570  oomph_info
571  << "First/last sample point: "
572  << sample_point_container_pt->
573  first_sample_point_to_actually_lookup_during_locate_zeta() << " "
574  << sample_point_container_pt->
575  last_sample_point_to_actually_lookup_during_locate_zeta()
576  << std::endl;
577 
578  }
579  // Non-refineable
580  else if (mesh_geom_obj_pt->sample_point_container_version()==
581  UseNonRefineableBinArray)
582  {
583  NonRefineableBinArray* sample_point_container_pt=
584  dynamic_cast<NonRefineableBinArray*>(mesh_geom_obj_pt->
585  sample_point_container_pt());
586 
587  oomph_info << "First/last spiral level to be visited: "
588  << sample_point_container_pt->current_min_spiral_level() << " "
589  << sample_point_container_pt->current_max_spiral_level() << " "
590  << std::endl;
591  }
592 #ifdef OOMPH_HAS_CGAL
593  // CGAL
594  else if (mesh_geom_obj_pt->sample_point_container_version()==
595  UseCGALSamplePointContainer)
596  {
597  CGALSamplePointContainer* sample_point_container_pt=
598  dynamic_cast<CGALSamplePointContainer*>(mesh_geom_obj_pt->
599  sample_point_container_pt());
600  oomph_info
601  << "First/last sample point: "
602  << sample_point_container_pt->
603  first_sample_point_to_actually_lookup_during_locate_zeta() << " "
604  << sample_point_container_pt->
605  last_sample_point_to_actually_lookup_during_locate_zeta()
606  << std::endl;
607  }
608 #endif
609  else
610  {
611  oomph_info << "wrong sample point container\n";
612  abort();
613  }
614 
615  BinArray::Always_fail_elemental_locate_zeta=true;
616  std::stringstream filename;
617  filename << "RESLT/sample_points_" << n_dim << "d";
618  if (ref_flag==0)
619  {
620  filename << "_nonref_bin";
621  }
622  else if (ref_flag==1)
623  {
624  filename << "_ref_bin";
625  }
626  else if (ref_flag==2)
627  {
628  filename << "cgal";
629  }
630  else
631  {
632  oomph_info << "wrong container\n";
633  abort();
634  }
635  filename << "_spiral_level"+StringConversion::to_string(i_level)+".dat";
636  BinArray::Visited_sample_points_file.open(filename.str().c_str());
637  Vector<double> s_from_locate_zeta(n_dim);
638  Vector<double> zeta_test(n_dim,0.5);
639  double start_t = TimingHelpers::timer();
640  mesh_geom_obj_pt->locate_zeta(zeta,geom_obj_pt,s_from_locate_zeta);
641  double end_t = TimingHelpers::timer();
642  t_locate_zeta+=end_t-start_t;
643 
644 
645  // Reset and close
646  BinArray::Always_fail_elemental_locate_zeta=false;
647  BinArray::Visited_sample_points_file.close();
648 
649 
650 
651  // Bail
652  if (geom_obj_pt!=0)
653  {
654  oomph_info << "Found zeta in " << i_level << " search levels.\n";
655  break;
656  }
657 
658  // Bump
659 
660  // Refineable
661  if (mesh_geom_obj_pt->sample_point_container_version()==
662  UseRefineableBinArray)
663  {
664 
665  RefineableBinArray* sample_point_container_pt=
666  dynamic_cast<RefineableBinArray*>(mesh_geom_obj_pt->
667  sample_point_container_pt());
668  sample_point_container_pt->
669  first_sample_point_to_actually_lookup_during_locate_zeta() =
670  sample_point_container_pt->
671  last_sample_point_to_actually_lookup_during_locate_zeta();
672  sample_point_container_pt->
673  last_sample_point_to_actually_lookup_during_locate_zeta() *=
674  sample_point_container_pt->
675  multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta();
676  }
677  // Non refineable
678  else if (mesh_geom_obj_pt->sample_point_container_version()==
679  UseNonRefineableBinArray)
680  {
681  NonRefineableBinArray* sample_point_container_pt=
682  dynamic_cast<NonRefineableBinArray*>(mesh_geom_obj_pt->
683  sample_point_container_pt());
684 
685  sample_point_container_pt->current_min_spiral_level()+=sample_point_container_pt->n_spiral_chunk();
686  sample_point_container_pt->current_max_spiral_level()+=sample_point_container_pt->n_spiral_chunk();
687  }
688 #ifdef OOMPH_HAS_CGAL
689  // CGAL
690  else if (mesh_geom_obj_pt->sample_point_container_version()==
691  UseCGALSamplePointContainer)
692  {
693 
694  CGALSamplePointContainer* sample_point_container_pt=
695  dynamic_cast<CGALSamplePointContainer*>(mesh_geom_obj_pt->
696  sample_point_container_pt());
697  sample_point_container_pt->
698  first_sample_point_to_actually_lookup_during_locate_zeta() =
699  sample_point_container_pt->
700  last_sample_point_to_actually_lookup_during_locate_zeta();
701  sample_point_container_pt->
702  last_sample_point_to_actually_lookup_during_locate_zeta() *=
703  sample_point_container_pt->
704  multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta();
705  }
706 #endif
707  else
708  {
709  oomph_info << "wrong sample point container\n";
710  abort();
711  }
712 
713 
714  // Stop it because we're out of the bin array?
715 
716  // Refineable
717  if (mesh_geom_obj_pt->sample_point_container_version()==
718  UseRefineableBinArray)
719  {
720  RefineableBinArray* sample_point_container_pt=
721  dynamic_cast<RefineableBinArray*>(mesh_geom_obj_pt->
722  sample_point_container_pt());
723  if (sample_point_container_pt->
724  first_sample_point_to_actually_lookup_during_locate_zeta()
725  <= nsample_points_in_bin_array)
726  {
727  has_not_reached_max_level_of_search=true;
728  }
729  else
730  {
731  has_not_reached_max_level_of_search=false;
732  }
733  }
734  // Non-refineable
735  else if (mesh_geom_obj_pt->sample_point_container_version()==
736  UseNonRefineableBinArray)
737  {
738  NonRefineableBinArray* sample_point_container_pt=
739  dynamic_cast<NonRefineableBinArray*>(mesh_geom_obj_pt->
740  sample_point_container_pt());
741 
742  if (sample_point_container_pt->current_max_spiral_level() < n_max_level)
743  {
744  has_not_reached_max_level_of_search=true;
745  }
746  else
747  {
748  has_not_reached_max_level_of_search=false;
749  }
750 
751  }
752 #ifdef OOMPH_HAS_CGAL
753  // CGAL
754  else if (mesh_geom_obj_pt->sample_point_container_version()==
755  UseCGALSamplePointContainer)
756  {
757  CGALSamplePointContainer* sample_point_container_pt=
758  dynamic_cast<CGALSamplePointContainer*>(mesh_geom_obj_pt->
759  sample_point_container_pt());
760  if (sample_point_container_pt->
761  first_sample_point_to_actually_lookup_during_locate_zeta()
762  <= nsample_points_in_bin_array)
763  {
764  has_not_reached_max_level_of_search=true;
765  }
766  else
767  {
768  has_not_reached_max_level_of_search=false;
769  }
770  }
771 #endif
772 
773  // Bump
774  i_level++;
775  }
776 
777 
778  // Test
779  oomph_info << "\n\n\nIt took " << t_locate_zeta << " secs to ";
780  if (geom_obj_pt!=0)
781  {
782  oomph_info << "find: "
783  << zeta[0] << " "
784  << zeta[1] << " "
785  << std::endl;
786  }
787  else
788  {
789  oomph_info << "not find: "
790  << zeta[0] << " "
791  << zeta[1] << " "
792  << std::endl;
793  }
794 
795  unsigned nvisit=mesh_geom_obj_pt->sample_point_container_pt()->
796  total_number_of_sample_points_visited_during_locate_zeta_from_top_level();
797  unsigned ntotal=mesh_geom_obj_pt->sample_point_container_pt()->
798  total_number_of_sample_points_computed_recursively();
799 
800  oomph_info << "Visited " << nvisit << " sample points "
801  << " out of a total of " << ntotal << std::endl;
802 
803  if (nvisit!=ntotal)
804  {
805  oomph_info << "ERROR\n";
806  abort();
807  }
808 
809 
810  }
811 
812 
813  // Clean up
814  delete mesh_geom_obj_pt;
815  mesh_geom_obj_pt=0;
816 
817  // If we got to here without aborting; we're done
818 
819  if (ref_flag==0)
820  {
821  std::string filename="RESLT/success_non_ref_bin.dat";
822  ofstream outfile;
823  outfile.open(filename.c_str());
824  outfile << "1" << std::endl;
825  outfile.close();
826  }
827  else if (ref_flag==1)
828  {
829  std::string filename="RESLT/success_ref_bin.dat";
830  ofstream outfile;
831  outfile.open(filename.c_str());
832  outfile << "1" << std::endl;
833  outfile.close();
834  }
835  else if (ref_flag==2)
836  {
837  std::string filename="RESLT/success_cgal.dat";
838  ofstream outfile;
839  outfile.open(filename.c_str());
840  outfile << "1" << std::endl;
841  outfile.close();
842  }
843  else
844  {
845  oomph_info << "wrong container\n";
846  abort();
847  }
848 
849  }
850 
851 
852 #ifndef OOMPH_HAS_CGAL
853  // Dummy output
854  {
855  std::string filename="RESLT/success_cgal.dat";
856  ofstream outfile;
857  outfile.open(filename.c_str());
858  outfile << "1" << std::endl;
859  outfile.close();
860  }
861 #endif
862 
863 }
cos
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
sin
AnnoyingScalar sin(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:137
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
BinArray
Base class for all bin arrays.
Definition: sample_point_container.h:402
BinArray::max_bin_dimension
unsigned max_bin_dimension() const
Max. bin dimension (number of bins in coordinate directions)
Definition: sample_point_container.cc:659
Eigen::ArrayBase::pow
EIGEN_DEVICE_FUNC const GlobalUnaryPowReturnType< Derived, ScalarExponent > pow(const Eigen::ArrayBase< Derived > &x, const ScalarExponent &exponent)
Definition: GlobalFunctions.h:137
NonRefineableBinArray
NonRefineableBinArray class.
Definition: sample_point_container.h:819
NonRefineableBinArray::n_spiral_chunk
unsigned n_spiral_chunk() const
Number of spirals to be searched in one go const version.
Definition: sample_point_container.h:866
NonRefineableBinArray::current_max_spiral_level
unsigned & current_max_spiral_level()
Access function to current max. spiral level.
Definition: sample_point_container.h:891
NonRefineableBinArray::current_min_spiral_level
unsigned & current_min_spiral_level()
Access function to current min. spiral level.
Definition: sample_point_container.h:885
RefineableBinArray
RefineableBinArray class.
Definition: sample_point_container.h:521
SamplePointContainer::Always_fail_elemental_locate_zeta
static bool Always_fail_elemental_locate_zeta
Boolean flag to make to make locate zeta fail.
Definition: sample_point_container.h:339
SamplePointContainer::Visited_sample_points_file
static std::ofstream Visited_sample_points_file
File to record sequence of visited sample points in.
Definition: sample_point_container.h:335
s
RealScalar s
Definition: level1_cplx_impl.h:130
BiharmonicTestFunctions2::Pi
double Pi
Definition: two_d_biharmonic.cc:235
BiharmonicTestFunctions2::theta
double theta
Definition: two_d_biharmonic.cc:236
Eigen::zeta
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:152
MergeRestartFiles.filename
string filename
Definition: MergeRestartFiles.py:39
MeshRefinement::max_error
double max_error
Definition: MortaringCantileverCompareToNonMortaring.cpp:188
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
calibrate.error
int error
Definition: calibrate.py:297
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
oomph::StringConversion::to_string
std::string to_string(T object, unsigned float_precision=8)
Definition: oomph_utilities.h:189
oomph::UseRefineableBinArray
@ UseRefineableBinArray
Definition: sample_point_parameters.h:41
oomph::UseNonRefineableBinArray
@ UseNonRefineableBinArray
Definition: sample_point_parameters.h:42
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319
timer
double timer
Definition: oomph_metis_from_parmetis_3.1.1/struct.h:210
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References SamplePointContainer::Always_fail_elemental_locate_zeta, cos(), NonRefineableBinArray::current_max_spiral_level(), NonRefineableBinArray::current_min_spiral_level(), e(), calibrate::error, boost::multiprecision::fabs(), MergeRestartFiles::filename, i, j, BinArray::max_bin_dimension(), MeshRefinement::max_error, NonRefineableBinArray::n_spiral_chunk(), oomph::oomph_info, BiharmonicTestFunctions2::Pi, Eigen::ArrayBase< Derived >::pow(), s, sin(), sqrt(), oomph::Global_string_for_annotation::string(), BiharmonicTestFunctions2::theta, oomph::StringConversion::to_string(), oomph::UseNonRefineableBinArray, oomph::UseRefineableBinArray, SamplePointContainer::Visited_sample_points_file, and Eigen::zeta().

Referenced by main().