sample_point_container.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
27
28
29namespace oomph
30{
31 /// ///////////////////////////////////////////////////////////////////////////
32 /// ///////////////////////////////////////////////////////////////////////////
33 /// ///////////////////////////////////////////////////////////////////////////
34 // RefineableBin class
35 /// ///////////////////////////////////////////////////////////////////////////
36 /// ///////////////////////////////////////////////////////////////////////////
37 /// ///////////////////////////////////////////////////////////////////////////
38
39 //==============================================================================
40 /// Destructor
41 //==============================================================================
43 {
44 if (Sub_bin_array_pt != 0)
45 {
46 delete Sub_bin_array_pt;
47 }
48 Sub_bin_array_pt = 0;
49
50 if (Sample_point_pt != 0)
51 {
52 unsigned n = Sample_point_pt->size();
53 for (unsigned i = 0; i < n; i++)
54 {
55 delete (*Sample_point_pt)[i];
56 }
57 delete Sample_point_pt;
58 }
59 }
60
61
62 //==============================================================================
63 /// Compute total number of sample points recursively
64 //==============================================================================
66 const
67 {
68 unsigned count = 0;
69
70 // Recurse?
71 if (Sub_bin_array_pt != 0)
72 {
73 count =
74 Sub_bin_array_pt->total_number_of_sample_points_computed_recursively();
75 }
76 else
77 {
78 if (Sample_point_pt != 0)
79 {
80 count = Sample_point_pt->size();
81 }
82 }
83 return count;
84 }
85
86
87 //==============================================================================
88 /// Function called for making a sub bin array in a given RefineableBin.
89 /// Pass the vector of min and max coordinates of the NEW bin array.
90 //==============================================================================
92 const Vector<std::pair<double, double>>& min_and_max_coordinates)
93 {
94 // Setup parameters for sub-bin
95 RefineableBinArrayParameters* ref_bin_array_parameters_pt =
96 new RefineableBinArrayParameters(Bin_array_pt->mesh_pt());
97
98 // Pass coordinates and dimensions
99 ref_bin_array_parameters_pt->min_and_max_coordinates() =
100 min_and_max_coordinates;
101
102 ref_bin_array_parameters_pt->dimensions_of_bin_array() =
103 Bin_array_pt->dimensions_of_bin_array();
104
105 // Eulerian coordinates or zeta?
106 if (Bin_array_pt->use_eulerian_coordinates_during_setup())
107 {
108 ref_bin_array_parameters_pt
109 ->enable_use_eulerian_coordinates_during_setup();
110 }
111 else
112 {
113 ref_bin_array_parameters_pt
114 ->disable_use_eulerian_coordinates_during_setup();
115 }
116
117
118#ifdef OOMPH_HAS_MPI
119
120 // How do we handle halo elements?
121 if (Bin_array_pt->ignore_halo_elements_during_locate_zeta_search())
122 {
123 ref_bin_array_parameters_pt
124 ->enable_ignore_halo_elements_during_locate_zeta_search();
125 }
126 else
127 {
128 ref_bin_array_parameters_pt
129 ->disable_ignore_halo_elements_during_locate_zeta_search();
130 }
131
132#endif
133
134 // "Measure of" number of sample points per element
135 ref_bin_array_parameters_pt->nsample_points_generated_per_element() =
136 Bin_array_pt->nsample_points_generated_per_element();
137
138
139 // Is it recursive?
140 if (Bin_array_pt->bin_array_is_recursive())
141 {
142 ref_bin_array_parameters_pt->enable_bin_array_is_recursive();
143 }
144 else
145 {
146 ref_bin_array_parameters_pt->disable_bin_array_is_recursive();
147 }
148
149 // Depth
150 ref_bin_array_parameters_pt->depth() = Bin_array_pt->depth() + 1;
151
152 // Max. depth
153 ref_bin_array_parameters_pt->max_depth() = Bin_array_pt->max_depth();
154
155 // Max. number of sample points before it's subdivided
156 ref_bin_array_parameters_pt->max_number_of_sample_point_per_bin() =
157 Bin_array_pt->max_number_of_sample_point_per_bin();
158
159 // Root bin array
160 ref_bin_array_parameters_pt->root_bin_array_pt() =
161 Bin_array_pt->root_bin_array_pt();
162
163 // We first construct a new bin array, providing the right parameters
164 BinArrayParameters* bin_array_parameters_pt = ref_bin_array_parameters_pt;
165 Sub_bin_array_pt = new RefineableBinArray(bin_array_parameters_pt);
166 delete ref_bin_array_parameters_pt;
167
168 // Fill it
169 Sub_bin_array_pt->fill_bin_array(*Sample_point_pt);
170
171 // Now deleting Sample_point_pt, we no longer need it; note that
172 // sample points themselves stay alive!
173 delete Sample_point_pt;
174 Sample_point_pt = 0;
175 }
176
177
178 //============================================================================
179 /// Output bin; x,[y,[z]],n_sample_points.
180 //============================================================================
181 void RefineableBin::output(std::ofstream& outfile, const bool& don_t_recurse)
182 {
183 // Recurse?
184 if ((Sub_bin_array_pt != 0) && (!don_t_recurse))
185 {
186 Sub_bin_array_pt->output_bins(outfile);
187 }
188 else
189 {
190 unsigned n_lagr = Bin_array_pt->ndim_zeta();
191 Vector<std::pair<double, double>> min_and_max_coordinates(n_lagr);
192 get_bin_boundaries(min_and_max_coordinates);
193
194 // How many sample points do we have in this bin?
195 unsigned n_sample_points = 0;
196 if (Sample_point_pt != 0)
197 {
198 n_sample_points = Sample_point_pt->size();
199 }
200
201 switch (n_lagr)
202 {
203 case 1:
204 outfile << "ZONE I=2\n"
205 << min_and_max_coordinates[0].first << " " << n_sample_points
206 << std::endl
207 << min_and_max_coordinates[0].second << " " << n_sample_points
208 << std::endl;
209 break;
210
211 case 2:
212
213 outfile << "ZONE I=2, J=2\n"
214 << min_and_max_coordinates[0].first << " "
215 << min_and_max_coordinates[1].first << " " << n_sample_points
216 << "\n"
217
218 << min_and_max_coordinates[0].second << " "
219 << min_and_max_coordinates[1].first << " " << n_sample_points
220 << "\n"
221
222 << min_and_max_coordinates[0].first << " "
223 << min_and_max_coordinates[1].second << " " << n_sample_points
224 << "\n"
225
226 << min_and_max_coordinates[0].second << " "
227 << min_and_max_coordinates[1].second << " " << n_sample_points
228 << "\n";
229 break;
230
231 case 3:
232
233
234 outfile << "ZONE I=2, J=2, K=2\n"
235 << min_and_max_coordinates[0].first << " "
236 << min_and_max_coordinates[1].first << " "
237 << min_and_max_coordinates[2].first << " " << n_sample_points
238 << "\n"
239
240 << min_and_max_coordinates[0].second << " "
241 << min_and_max_coordinates[1].first << " "
242 << min_and_max_coordinates[2].first << " " << n_sample_points
243 << "\n"
244
245 << min_and_max_coordinates[0].first << " "
246 << min_and_max_coordinates[1].second << " "
247 << min_and_max_coordinates[2].first << " " << n_sample_points
248 << "\n"
249
250 << min_and_max_coordinates[0].second << " "
251 << min_and_max_coordinates[1].second << " "
252 << min_and_max_coordinates[2].first << " " << n_sample_points
253 << "\n"
254
255 << min_and_max_coordinates[0].first << " "
256 << min_and_max_coordinates[1].first << " "
257 << min_and_max_coordinates[2].second << " " << n_sample_points
258 << "\n"
259
260 << min_and_max_coordinates[0].second << " "
261 << min_and_max_coordinates[1].first << " "
262 << min_and_max_coordinates[2].second << " " << n_sample_points
263 << "\n"
264
265 << min_and_max_coordinates[0].first << " "
266 << min_and_max_coordinates[1].second << " "
267 << min_and_max_coordinates[2].second << " " << n_sample_points
268 << "\n"
269
270 << min_and_max_coordinates[0].second << " "
271 << min_and_max_coordinates[1].second << " "
272 << min_and_max_coordinates[2].second << " " << n_sample_points
273 << "\n";
274
275 break;
276
277 default:
278
279 oomph_info << "n_lagr=" << n_lagr << std::endl;
280 throw OomphLibError("Wrong dimension",
281 OOMPH_CURRENT_FUNCTION,
282 OOMPH_EXCEPTION_LOCATION);
283 }
284 }
285 }
286
287
288 //============================================================================
289 /// Output bin; x,[y,[z]]
290 //============================================================================
291 void RefineableBin::output_bin_vertices(std::ofstream& outfile)
292 {
293 // Recurse?
294 if (Sub_bin_array_pt != 0)
295 {
296 Sub_bin_array_pt->output_bin_vertices(outfile);
297 }
298 else
299 {
300 unsigned n_lagr = Bin_array_pt->ndim_zeta();
301 Vector<std::pair<double, double>> min_and_max_coordinates(n_lagr);
302 get_bin_boundaries(min_and_max_coordinates);
303
304 switch (n_lagr)
305 {
306 case 1:
307 outfile << "ZONE I=2\n"
308 << min_and_max_coordinates[0].first << std::endl
309 << min_and_max_coordinates[0].second << std::endl;
310 break;
311
312 case 2:
313
314 outfile << "ZONE I=2, J=2\n"
315 << min_and_max_coordinates[0].first << " "
316 << min_and_max_coordinates[1].first << " "
317 << "\n"
318
319 << min_and_max_coordinates[0].second << " "
320 << min_and_max_coordinates[1].first << " "
321 << "\n"
322
323 << min_and_max_coordinates[0].first << " "
324 << min_and_max_coordinates[1].second << " "
325 << "\n"
326
327 << min_and_max_coordinates[0].second << " "
328 << min_and_max_coordinates[1].second << " "
329 << "\n";
330 break;
331
332 case 3:
333
334
335 outfile << "ZONE I=2, J=2, K=2\n"
336 << min_and_max_coordinates[0].first << " "
337 << min_and_max_coordinates[1].first << " "
338 << min_and_max_coordinates[2].first << " "
339 << "\n"
340
341 << min_and_max_coordinates[0].second << " "
342 << min_and_max_coordinates[1].first << " "
343 << min_and_max_coordinates[2].first << " "
344 << "\n"
345
346 << min_and_max_coordinates[0].first << " "
347 << min_and_max_coordinates[1].second << " "
348 << min_and_max_coordinates[2].first << " "
349 << "\n"
350
351 << min_and_max_coordinates[0].second << " "
352 << min_and_max_coordinates[1].second << " "
353 << min_and_max_coordinates[2].first << " "
354 << "\n"
355
356 << min_and_max_coordinates[0].first << " "
357 << min_and_max_coordinates[1].first << " "
358 << min_and_max_coordinates[2].second << " "
359 << "\n"
360
361 << min_and_max_coordinates[0].second << " "
362 << min_and_max_coordinates[1].first << " "
363 << min_and_max_coordinates[2].second << " "
364 << "\n"
365
366 << min_and_max_coordinates[0].first << " "
367 << min_and_max_coordinates[1].second << " "
368 << min_and_max_coordinates[2].second << " "
369 << "\n"
370
371 << min_and_max_coordinates[0].second << " "
372 << min_and_max_coordinates[1].second << " "
373 << min_and_max_coordinates[2].second << " "
374 << "\n";
375
376 break;
377
378 default:
379
380 oomph_info << "n_lagr=" << n_lagr << std::endl;
381 throw OomphLibError("Wrong dimension",
382 OOMPH_CURRENT_FUNCTION,
383 OOMPH_EXCEPTION_LOCATION);
384 }
385 }
386 }
387
388
389 //==============================================================================
390 /// Add a SamplePoint* to a RefineableBin object.
391 //==============================================================================
393 const Vector<double>& zeta_coordinates)
394 {
395 // If the bin is a "leaf" (ie no sub bin array)
396 if (Sub_bin_array_pt == 0)
397 {
398 // if there is no Sample_point_pt create it
399 if (Sample_point_pt == 0)
400 {
401 Sample_point_pt = new Vector<SamplePoint*>;
402 }
403 this->Sample_point_pt->push_back(new_sample_point_pt);
404
405 // If we are recursive (ie not at the maximum depth or not
406 // in fill bin by diffusion configuration) and if the number
407 // of elements there are in the RefineableBin is bigger than the maximum
408 // one
409 if ((Bin_array_pt->bin_array_is_recursive()) &&
410 (Sample_point_pt->size() >
411 Bin_array_pt->max_number_of_sample_point_per_bin()) &&
412 (Bin_array_pt->depth() < Bin_array_pt->max_depth()))
413 {
414 // Get min and max coordinates of current bin...
415 Vector<std::pair<double, double>> min_and_max_coordinates(
416 Bin_array_pt->ndim_zeta());
417 get_bin_boundaries(min_and_max_coordinates);
418
419 // ...and use them as the boundaries for new sub-bin-array
420 // (this transfers all the new points into the new sub-bin-array
421 this->make_sub_bin_array(min_and_max_coordinates);
422 }
423 }
424 else // if the bin has a sub bin array
425 {
426 // we call the corresponding method of the sub bin array
427 this->Sub_bin_array_pt->add_sample_point(new_sample_point_pt,
428 zeta_coordinates);
429 }
430 }
431
432
433 //==============================================================================
434 /// Find sub-GeomObject (finite element) and the local coordinate
435 /// s within it that contains point with global coordinate zeta.
436 /// sub_geom_object_pt=0 if point can't be found.
437 //==============================================================================
438 void RefineableBin::locate_zeta(const Vector<double>& zeta,
439 GeomObject*& sub_geom_object_pt,
440 Vector<double>& s)
441 {
442 // Haven't found zeta yet!
443 sub_geom_object_pt = 0;
444
445 // Descend?
446 if (Sub_bin_array_pt != 0)
447 {
448 Sub_bin_array_pt->locate_zeta(zeta, sub_geom_object_pt, s);
449 }
450 else
451 {
452 // Do we have to look into any sample points in this bin?
453 // NOTE: There's some slight potential for overlap/duplication
454 // because we always search through all the sample points in a bin
455 // (unless we find the required point in which case we stop).
456 bool do_it = true;
457 if (
458 Bin_array_pt->root_bin_array_pt()
459 ->total_number_of_sample_points_visited_during_locate_zeta_from_top_level() <
460 Bin_array_pt->root_bin_array_pt()
461 ->first_sample_point_to_actually_lookup_during_locate_zeta())
462 {
463 // oomph_info << "Not doing it (ref-bin) because counter less than
464 // first" << std::endl;
465 do_it = false;
466 }
467 if (
468 Bin_array_pt->root_bin_array_pt()
469 ->total_number_of_sample_points_visited_during_locate_zeta_from_top_level() >
470 Bin_array_pt->root_bin_array_pt()
471 ->last_sample_point_to_actually_lookup_during_locate_zeta())
472 {
473 // oomph_info << "Not doing it (ref-bin) because counter more than last"
474 // << std::endl;
475 do_it = false;
476 }
477 double max_search_radius =
478 Bin_array_pt->root_bin_array_pt()->max_search_radius();
479 bool dont_do_it_because_of_radius = false;
480 if (max_search_radius < DBL_MAX)
481 {
482 // Base "radius" of bin on centre of gravity
483 unsigned n = zeta.size();
484 double dist_squared = 0.0;
485 double cog = 0.0;
486 double aux = 0.0;
487 Vector<std::pair<double, double>> min_and_max_coordinates(n);
488 get_bin_boundaries(min_and_max_coordinates);
489 for (unsigned i = 0; i < n; i++)
490 {
491 cog = 0.5 * (min_and_max_coordinates[i].first +
492 min_and_max_coordinates[i].second);
493 aux = (cog - zeta[i]);
494 dist_squared += aux * aux;
495 }
496 if (dist_squared > max_search_radius * max_search_radius)
497 {
498 do_it = false;
499 dont_do_it_because_of_radius = true;
500 }
501 }
502
503 if (!do_it)
504 {
505 if (!dont_do_it_because_of_radius)
506 {
507 // Skip all the sample points in this bin
508 Bin_array_pt->root_bin_array_pt()
509 ->total_number_of_sample_points_visited_during_locate_zeta_from_top_level() +=
510 Sample_point_pt->size();
511 }
512 return;
513 }
514
515
516 // Now search through (at most) all the sample points in this bin
517 unsigned n_sample_point = Sample_point_pt->size();
518 unsigned i = 0;
519 while ((i < n_sample_point) && (sub_geom_object_pt == 0))
520 {
521 // Get the corresponding finite element
522 FiniteElement* el_pt = Bin_array_pt->mesh_pt()->finite_element_pt(
523 (*Sample_point_pt)[i]->element_index_in_mesh());
524
525#ifdef OOMPH_HAS_MPI
526 // We only look at the sample point if it isn't halo
527 // if we are set up to ignore the halo elements
528 if ((Bin_array_pt->ignore_halo_elements_during_locate_zeta_search()) &&
529 (el_pt->is_halo()))
530 {
531 // Halo
532 }
533 else
534 {
535#endif
536 // Provide initial guess for Newton search using local coordinate
537 // of sample point
538 bool use_equally_spaced_interior_sample_points =
540 unsigned j = (*Sample_point_pt)[i]->sample_point_index_in_element();
541 el_pt->get_s_plot(
542 j,
543 Bin_array_pt->nsample_points_generated_per_element(),
544 s,
545 use_equally_spaced_interior_sample_points);
546
547
548 // History of sample points visited
550 {
551 unsigned cached_dim_zeta = Bin_array_pt->ndim_zeta();
552 Vector<double> zeta_sample(cached_dim_zeta);
553 if (Bin_array_pt->use_eulerian_coordinates_during_setup())
554 {
555 el_pt->interpolated_x(s, zeta_sample);
556 }
557 else
558 {
559 el_pt->interpolated_zeta(s, zeta_sample);
560 }
561 double dist = 0.0;
562 for (unsigned ii = 0; ii < cached_dim_zeta; ii++)
563 {
564 BinArray::Visited_sample_points_file << zeta_sample[ii] << " ";
565 dist +=
566 (zeta[ii] - zeta_sample[ii]) * (zeta[ii] - zeta_sample[ii]);
567 }
569 << Bin_array_pt->root_bin_array_pt()
570 ->total_number_of_sample_points_visited_during_locate_zeta_from_top_level()
571 << " " << sqrt(dist) << std::endl;
572 }
573
574
575 // Bump counter
576 Bin_array_pt->root_bin_array_pt()
577 ->total_number_of_sample_points_visited_during_locate_zeta_from_top_level()++;
578
579 bool use_coordinate_as_initial_guess = true;
580 el_pt->locate_zeta(
581 zeta, sub_geom_object_pt, s, use_coordinate_as_initial_guess);
582
583 // Always fail? (Used for debugging, e.g. to trace out
584 // spiral path)
586 {
587 sub_geom_object_pt = 0;
588 }
589
590#ifdef OOMPH_HAS_MPI
591 }
592#endif
593 // Next one please
594 i++;
595 }
596 }
597 }
598
599 //==============================================================================
600 /// Boundaries of bin in each coordinate direction. *.first = min;
601 /// *.second = max.
602 //==============================================================================
604 Vector<std::pair<double, double>>& min_and_max_coordinates)
605 {
606 unsigned n_bin = Bin_index_in_bin_array;
607
608 // temporary storage for the eulerian dim
609 unsigned current_dim = Bin_array_pt->ndim_zeta();
610 min_and_max_coordinates.resize(current_dim);
611 for (unsigned u = 0; u < current_dim; u++)
612 {
613 // The number of bins there are according to the u-th dimension
614 double nbin_in_dir = n_bin % Bin_array_pt->dimension_of_bin_array(u);
615 n_bin /= Bin_array_pt->dimension_of_bin_array(u);
616
617 // The range between the maximum and minimum u-th coordinates of a bin
618 double range = (Bin_array_pt->min_and_max_coordinates(u).second -
619 Bin_array_pt->min_and_max_coordinates(u).first) /
620 double(Bin_array_pt->dimension_of_bin_array(u));
621
622 // Now updating the minimum and maximum u-th coordinates for this bin.
623 min_and_max_coordinates[u].first =
624 Bin_array_pt->min_and_max_coordinates(u).first + nbin_in_dir * range;
625 min_and_max_coordinates[u].second =
626 min_and_max_coordinates[u].first + range;
627 }
628 }
629
630
631 /// ///////////////////////////////////////////////////////////////////////////
632 /// ///////////////////////////////////////////////////////////////////////////
633 /// ///////////////////////////////////////////////////////////////////////////
634 /// SamplePointContainer base class
635 /// ///////////////////////////////////////////////////////////////////////////
636 /// ///////////////////////////////////////////////////////////////////////////
637 /// ///////////////////////////////////////////////////////////////////////////
638
639 /// File to record sequence of visited sample points in
641
642 /// Boolean flag to make to make locate zeta fail
644
645 /// Use equally spaced sample points? (otherwise vertices are sampled
646 /// repeatedly
648
649 /// Time setup?
651
652 /// Offset of sample point container boundaries beyond max/min coords
654
655
656 //==============================================================================
657 /// Max. bin dimension (number of bins in coordinate directions)
658 //==============================================================================
660 {
661 unsigned dim = ndim_zeta();
662 unsigned n_max_level = Dimensions_of_bin_array[0];
663 if (dim >= 2)
664 {
665 if (Dimensions_of_bin_array[1] > n_max_level)
666 {
667 n_max_level = Dimensions_of_bin_array[1];
668 }
669 }
670 if (dim == 3)
671 {
672 if (Dimensions_of_bin_array[2] > n_max_level)
673 {
674 n_max_level = Dimensions_of_bin_array[2];
675 }
676 }
677 return n_max_level;
678 }
679
680
681 //========================================================================
682 /// Setup the min and max coordinates for the mesh, in each dimension
683 //========================================================================
685 {
686 // Get the lagrangian dimension
687 int n_lagrangian = ndim_zeta();
688
689 // Storage locally (i.e. in parallel on each processor)
690 // for the minimum and maximum coordinates
691 double zeta_min_local[n_lagrangian];
692 double zeta_max_local[n_lagrangian];
693 for (int i = 0; i < n_lagrangian; i++)
694 {
695 zeta_min_local[i] = DBL_MAX;
696 zeta_max_local[i] = -DBL_MAX;
697 }
698
699 // Loop over the elements of the mesh
700 unsigned n_el = Mesh_pt->nelement();
701 for (unsigned e = 0; e < n_el; e++)
702 {
703 FiniteElement* el_pt = Mesh_pt->finite_element_pt(e);
704
705 // Get the number of vertices (nplot=2 does the trick)
706 unsigned n_plot = 2;
707 unsigned n_plot_points = el_pt->nplot_points(n_plot);
708
709 // Loop over the number of plot points
710 for (unsigned iplot = 0; iplot < n_plot_points; iplot++)
711 {
712 Vector<double> s_local(n_lagrangian);
713 Vector<double> zeta_global(n_lagrangian);
714
715 // Get the local s -- need to sample over the entire range
716 // of the elements!
717 bool use_equally_spaced_interior_sample_points = false;
718 el_pt->get_s_plot(
719 iplot, n_plot, s_local, use_equally_spaced_interior_sample_points);
720
721 // Now interpolate to global coordinates
722 if (Use_eulerian_coordinates_during_setup)
723 {
724 el_pt->interpolated_x(s_local, zeta_global);
725 }
726 else
727 {
728 el_pt->interpolated_zeta(s_local, zeta_global);
729 }
730
731 // Check the max and min in each direction
732 for (int i = 0; i < n_lagrangian; i++)
733 {
734 // Is the coordinate less than the minimum?
735 if (zeta_global[i] < zeta_min_local[i])
736 {
737 zeta_min_local[i] = zeta_global[i];
738 }
739 // Is the coordinate bigger than the maximum?
740 if (zeta_global[i] > zeta_max_local[i])
741 {
742 zeta_max_local[i] = zeta_global[i];
743 }
744 }
745 }
746 }
747
748 // Global extrema - in parallel, need to get max/min across all processors
749 double zeta_min[n_lagrangian];
750 double zeta_max[n_lagrangian];
751 for (int i = 0; i < n_lagrangian; i++)
752 {
753 zeta_min[i] = 0.0;
754 zeta_max[i] = 0.0;
755 }
756
757#ifdef OOMPH_HAS_MPI
758 // If the mesh has been distributed and we want consistent bins
759 // across all processors
760 if (Mesh_pt->is_mesh_distributed())
761 {
762 // .. we need a non-null communicator!
763 if (Mesh_pt->communicator_pt() != 0)
764 {
765 int n_proc = Mesh_pt->communicator_pt()->nproc();
766 if (n_proc > 1)
767 {
768 // Get the minima and maxima over all processors
769 MPI_Allreduce(zeta_min_local,
770 zeta_min,
771 n_lagrangian,
772 MPI_DOUBLE,
773 MPI_MIN,
774 Mesh_pt->communicator_pt()->mpi_comm());
775 MPI_Allreduce(zeta_max_local,
776 zeta_max,
777 n_lagrangian,
778 MPI_DOUBLE,
779 MPI_MAX,
780 Mesh_pt->communicator_pt()->mpi_comm());
781 }
782 }
783 else // Null communicator - throw an error
784 {
785 std::ostringstream error_message_stream;
786 error_message_stream << "Communicator not set for a Mesh\n"
787 << "that was created from a distributed Mesh";
788 throw OomphLibError(error_message_stream.str(),
789 OOMPH_CURRENT_FUNCTION,
790 OOMPH_EXCEPTION_LOCATION);
791 }
792 }
793 else // If the mesh hasn't been distributed then the
794 // max and min are the same on all processors
795 {
796 for (int i = 0; i < n_lagrangian; i++)
797 {
798 zeta_min[i] = zeta_min_local[i];
799 zeta_max[i] = zeta_max_local[i];
800 }
801 }
802#else // If we're not using MPI then the mesh can't be distributed
803 for (int i = 0; i < n_lagrangian; i++)
804 {
805 zeta_min[i] = zeta_min_local[i];
806 zeta_max[i] = zeta_max_local[i];
807 }
808#endif
809
810 // Decrease/increase min and max to allow for any overshoot in
811 // meshes that may move around
812 // There's no point in doing this for DIM_LAGRANGIAN==1
813 for (int i = 0; i < n_lagrangian; i++)
814 {
815 double length = zeta_max[i] - zeta_min[i];
816 zeta_min[i] -= ((Percentage_offset / 100.0) * length);
817 zeta_max[i] += ((Percentage_offset / 100.0) * length);
818 }
819
820 // Set the entries as the Min/Max_coords vector
821 Min_and_max_coordinates.resize(n_lagrangian);
822 for (int i = 0; i < n_lagrangian; i++)
823 {
824 Min_and_max_coordinates[i].first = zeta_min[i];
825 Min_and_max_coordinates[i].second = zeta_max[i];
826 }
827 }
828
829
830 /// ///////////////////////////////////////////////////////////////////////////
831 /// ///////////////////////////////////////////////////////////////////////////
832 /// ///////////////////////////////////////////////////////////////////////////
833 /// RefineableBin array class
834 /// ///////////////////////////////////////////////////////////////////////////
835 /// ///////////////////////////////////////////////////////////////////////////
836 /// ///////////////////////////////////////////////////////////////////////////
837
838 //==============================================================================
839 /// Constructor
840 //==============================================================================
842 SamplePointContainerParameters* sample_point_container_parameters_pt)
844 sample_point_container_parameters_pt->mesh_pt(),
845 sample_point_container_parameters_pt->min_and_max_coordinates(),
846 sample_point_container_parameters_pt
847 ->use_eulerian_coordinates_during_setup(),
848 sample_point_container_parameters_pt
849 ->ignore_halo_elements_during_locate_zeta_search(),
850 sample_point_container_parameters_pt
851 ->nsample_points_generated_per_element()),
852 BinArray(
853 sample_point_container_parameters_pt->mesh_pt(),
854 sample_point_container_parameters_pt->min_and_max_coordinates(),
855 dynamic_cast<BinArrayParameters*>(sample_point_container_parameters_pt)
856 ->dimensions_of_bin_array(),
857 sample_point_container_parameters_pt
858 ->use_eulerian_coordinates_during_setup(),
859 sample_point_container_parameters_pt
860 ->ignore_halo_elements_during_locate_zeta_search(),
861 sample_point_container_parameters_pt
862 ->nsample_points_generated_per_element())
863 {
864 RefineableBinArrayParameters* ref_bin_array_parameters_pt =
865 dynamic_cast<RefineableBinArrayParameters*>(
866 sample_point_container_parameters_pt);
867
868#ifdef PARANOID
869 if (ref_bin_array_parameters_pt == 0)
870 {
871 throw OomphLibError("Wrong sample_point_container_parameters_pt",
872 OOMPH_CURRENT_FUNCTION,
873 OOMPH_EXCEPTION_LOCATION);
874 }
875#endif
876
878 ref_bin_array_parameters_pt->bin_array_is_recursive();
879 Depth = ref_bin_array_parameters_pt->depth();
880 Max_depth = ref_bin_array_parameters_pt->max_depth();
882 ref_bin_array_parameters_pt->max_number_of_sample_point_per_bin();
883 Root_bin_array_pt = ref_bin_array_parameters_pt->root_bin_array_pt();
884
885 // Set default size of bin array (and spatial dimension!)
886 if (Dimensions_of_bin_array.size() == 0)
887 {
888 int dim = 0;
889 if (Mesh_pt->nelement() != 0)
890 {
891 dim = Mesh_pt->finite_element_pt(0)->dim();
892 }
893
894 // Need to do an Allreduce to ensure that the dimension is consistent
895 // even when no elements are assigned to a certain processor
896#ifdef OOMPH_HAS_MPI
897 // Only a problem if the mesh has been distributed
898 if (Mesh_pt->is_mesh_distributed())
899 {
900 // Need a non-null communicator
901 if (Mesh_pt->communicator_pt() != 0)
902 {
903 int n_proc = Mesh_pt->communicator_pt()->nproc();
904 if (n_proc > 1)
905 {
906 int dim_reduce;
907 MPI_Allreduce(&dim,
908 &dim_reduce,
909 1,
910 MPI_INT,
911 MPI_MAX,
912 Mesh_pt->communicator_pt()->mpi_comm());
913 dim = dim_reduce;
914 }
915 }
916 }
917#endif
918
920 }
921
922 // Have we specified max/min coordinates?
923 // If not, compute them on the fly from mesh
924 if (Min_and_max_coordinates.size() == 0)
925 {
927 }
928
929 // Get total number of bins and make space
930 unsigned dim = Dimensions_of_bin_array.size();
931 unsigned n_bin = 1;
932 for (unsigned i = 0; i < dim; i++)
933 {
934 n_bin *= Dimensions_of_bin_array[i];
935 }
936 Bin_pt.resize(n_bin, 0);
937
938 // I'm my own root bin array
939 if (Depth == 0)
940 {
941 Root_bin_array_pt = this;
942 }
943
944#ifdef PARANOID
945 if (Depth > 0)
946 {
947 if (Root_bin_array_pt == 0)
948 {
949 throw OomphLibError(
950 "Must specify root_bin_array for lower-level bin arrays\n",
951 OOMPH_CURRENT_FUNCTION,
952 OOMPH_EXCEPTION_LOCATION);
953 }
954 }
955#endif
956
957 // Initialise
962 2; // hierher tune this and create public static default
964 10; // hierher UINT MAX is temporary bypass! tune this and create public
965 // static default
966
967 // Now fill the bastard...
968 if (Depth == 0)
969 {
970 // Time it
971 double t_start = 0.0;
973 {
974 t_start = TimingHelpers::timer();
975 }
978 {
979 double t_end = TimingHelpers::timer();
981 oomph_info << "Time for setup of " << dim
982 << "-dimensional sample point container containing " << npts
983 << " sample points: " << t_end - t_start
984 << " sec (ref_bin); third party: 0 sec ( = 0 %)"
985 << std::endl;
986 }
987 }
988 }
989
990
991 //==============================================================================
992 /// Boundaries of specified bin in each coordinate direction.
993 /// *.first = min; *.second = max.
994 //==============================================================================
996 const unsigned& bin_index,
997 Vector<std::pair<double, double>>& min_and_max_coordinates_of_bin)
998 {
999 unsigned bin_index_local = bin_index;
1000
1001 // temporary storage for the eulerian dim
1002 unsigned current_dim = ndim_zeta();
1003 min_and_max_coordinates_of_bin.resize(current_dim);
1004 for (unsigned u = 0; u < current_dim; u++)
1005 {
1006 // The number of bins there are according to the u-th dimension
1007 unsigned nbin_in_dir = bin_index_local % dimension_of_bin_array(u);
1008 bin_index_local /= dimension_of_bin_array(u);
1009
1010 // The range between the maximum and minimum u-th coordinates of a bin
1011 double range =
1012 (Min_and_max_coordinates[u].second - Min_and_max_coordinates[u].first) /
1013 double(Dimensions_of_bin_array[u]);
1014
1015 // Now updating the minimum and maximum u-th coordinates for this bin.
1016 min_and_max_coordinates_of_bin[u].first =
1017 Min_and_max_coordinates[u].first + double(nbin_in_dir) * range;
1018 min_and_max_coordinates_of_bin[u].second =
1019 min_and_max_coordinates_of_bin[u].first + range;
1020 }
1021 }
1022
1023
1024 //==============================================================================
1025 /// Output neighbouring bins up to given "radius" of the specified bin
1026 //==============================================================================
1027 void RefineableBinArray::output_bin_vertices(std::ofstream& outfile)
1028 {
1029 // Loop over bins
1030 unsigned n_bin = Bin_pt.size();
1031 for (unsigned i = 0; i < n_bin; i++)
1032 {
1033 if (Bin_pt[i] != 0)
1034 {
1035 Bin_pt[i]->output_bin_vertices(outfile);
1036 }
1037 }
1038 }
1039
1040
1041 //==============================================================================
1042 /// Output neighbouring bins up to given "radius" of the specified bin
1043 //==============================================================================
1044 void RefineableBinArray::output_neighbouring_bins(const unsigned& bin_index,
1045 const unsigned& radius,
1046 std::ofstream& outfile)
1047 {
1048 unsigned n_lagr = ndim_zeta();
1049
1050 Vector<unsigned> neighbouring_bin_index;
1051 get_neighbouring_bins_helper(bin_index, radius, neighbouring_bin_index);
1052 unsigned nneigh = neighbouring_bin_index.size();
1053
1054 // Outline of bin structure
1055 switch (n_lagr)
1056 {
1057 case 1:
1058 outfile << "ZONE I=2\n"
1059 << Min_and_max_coordinates[0].first << std::endl
1060 << Min_and_max_coordinates[0].second << std::endl;
1061 break;
1062
1063 case 2:
1064
1065 outfile << "ZONE I=2, J=2\n"
1066 << Min_and_max_coordinates[0].first << " "
1067 << Min_and_max_coordinates[1].first << " "
1068 << "\n"
1069
1070 << Min_and_max_coordinates[0].second << " "
1071 << Min_and_max_coordinates[1].first << " "
1072 << "\n"
1073
1074 << Min_and_max_coordinates[0].first << " "
1075 << Min_and_max_coordinates[1].second << " "
1076 << "\n"
1077
1078 << Min_and_max_coordinates[0].second << " "
1079 << Min_and_max_coordinates[1].second << " "
1080 << "\n";
1081 break;
1082
1083 case 3:
1084
1085 outfile << "ZONE I=2, J=2, K=2 \n"
1086 << Min_and_max_coordinates[0].first << " "
1087 << Min_and_max_coordinates[1].first << " "
1088 << Min_and_max_coordinates[2].first << " "
1089 << "\n"
1090
1091 << Min_and_max_coordinates[0].second << " "
1092 << Min_and_max_coordinates[1].first << " "
1093 << Min_and_max_coordinates[2].first << " "
1094 << "\n"
1095
1096 << Min_and_max_coordinates[0].first << " "
1097 << Min_and_max_coordinates[1].second << " "
1098 << Min_and_max_coordinates[2].first << " "
1099 << "\n"
1100
1101 << Min_and_max_coordinates[0].second << " "
1102 << Min_and_max_coordinates[1].second << " "
1103 << Min_and_max_coordinates[2].first << " "
1104 << "\n"
1105
1106 << Min_and_max_coordinates[0].first << " "
1107 << Min_and_max_coordinates[1].first << " "
1108 << Min_and_max_coordinates[2].second << " "
1109 << "\n"
1110
1111 << Min_and_max_coordinates[0].second << " "
1112 << Min_and_max_coordinates[1].first << " "
1113 << Min_and_max_coordinates[2].second << " "
1114 << "\n"
1115
1116 << Min_and_max_coordinates[0].first << " "
1117 << Min_and_max_coordinates[1].second << " "
1118 << Min_and_max_coordinates[2].second << " "
1119 << "\n"
1120
1121 << Min_and_max_coordinates[0].second << " "
1122 << Min_and_max_coordinates[1].second << " "
1123 << Min_and_max_coordinates[2].second << " "
1124 << "\n";
1125 break;
1126
1127 default:
1128
1129 oomph_info << "n_lagr=" << n_lagr << std::endl;
1130 throw OomphLibError(
1131 "Wrong dimension!", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1132 }
1133
1134 // Loop over neighbours
1135 for (unsigned i = 0; i < nneigh; i++)
1136 {
1137 Vector<std::pair<double, double>> min_and_max_coordinates(n_lagr);
1138 get_bin_boundaries(neighbouring_bin_index[i], min_and_max_coordinates);
1139
1140 switch (n_lagr)
1141 {
1142 case 1:
1143 outfile << "ZONE I=2\n"
1144 << min_and_max_coordinates[0].first << std::endl
1145 << min_and_max_coordinates[0].second << std::endl;
1146 break;
1147
1148 case 2:
1149
1150 outfile << "ZONE I=2, J=2\n"
1151 << min_and_max_coordinates[0].first << " "
1152 << min_and_max_coordinates[1].first << " "
1153 << "\n"
1154
1155 << min_and_max_coordinates[0].second << " "
1156 << min_and_max_coordinates[1].first << " "
1157 << "\n"
1158
1159 << min_and_max_coordinates[0].first << " "
1160 << min_and_max_coordinates[1].second << " "
1161 << "\n"
1162
1163 << min_and_max_coordinates[0].second << " "
1164 << min_and_max_coordinates[1].second << " "
1165 << "\n";
1166 break;
1167
1168 case 3:
1169
1170 outfile << "ZONE I=2, J=2, K=2\n"
1171 << min_and_max_coordinates[0].first << " "
1172 << min_and_max_coordinates[1].first << " "
1173 << min_and_max_coordinates[2].first << " "
1174 << "\n"
1175
1176 << min_and_max_coordinates[0].second << " "
1177 << min_and_max_coordinates[1].first << " "
1178 << min_and_max_coordinates[2].first << " "
1179 << "\n"
1180
1181 << min_and_max_coordinates[0].first << " "
1182 << min_and_max_coordinates[1].second << " "
1183 << min_and_max_coordinates[2].first << " "
1184 << "\n"
1185
1186 << min_and_max_coordinates[0].second << " "
1187 << min_and_max_coordinates[1].second << " "
1188 << min_and_max_coordinates[2].first << " "
1189 << "\n"
1190
1191 << min_and_max_coordinates[0].first << " "
1192 << min_and_max_coordinates[1].first << " "
1193 << min_and_max_coordinates[2].second << " "
1194 << "\n"
1195
1196 << min_and_max_coordinates[0].second << " "
1197 << min_and_max_coordinates[1].first << " "
1198 << min_and_max_coordinates[2].second << " "
1199 << "\n"
1200
1201 << min_and_max_coordinates[0].first << " "
1202 << min_and_max_coordinates[1].second << " "
1203 << min_and_max_coordinates[2].second << " "
1204 << "\n"
1205
1206 << min_and_max_coordinates[0].second << " "
1207 << min_and_max_coordinates[1].second << " "
1208 << min_and_max_coordinates[2].second << " "
1209 << "\n";
1210 break;
1211
1212 default:
1213
1214 oomph_info << "n_lagr=" << n_lagr << std::endl;
1215 throw OomphLibError("Wrong dimension!",
1216 OOMPH_CURRENT_FUNCTION,
1217 OOMPH_EXCEPTION_LOCATION);
1218 }
1219 }
1220 }
1221
1222
1223 //========================================================================
1224 /// Profiling function to compare performance of two different
1225 /// versions of the get_neighbouring_bins_helper(...) function
1226 //========================================================================
1228 {
1229 unsigned old_faster = 0;
1230 unsigned new_faster = 0;
1231 double t_total_new = 0.0;
1232 double t_total_old = 0.0;
1233
1234 unsigned dim = ndim_zeta();
1235
1236 // Choose bin in middle of the domain
1237 Vector<double> zeta(dim);
1238 for (unsigned i = 0; i < dim; i++)
1239 {
1240 zeta[i] = 0.5 * (Min_and_max_coordinates[i].first +
1241 Min_and_max_coordinates[i].second);
1242 }
1243
1244 // Finding the bin in which the point is located
1245 int bin_index = coords_to_bin_index(zeta);
1246
1247#ifdef PARANOID
1248 if (bin_index < 0)
1249 {
1250 throw OomphLibError("Negative bin index...",
1251 OOMPH_CURRENT_FUNCTION,
1252 OOMPH_EXCEPTION_LOCATION);
1253 }
1254#endif
1255
1256 // Start at this radius (radius = 0 is the central bin)
1257 unsigned radius = 0;
1258
1259 // "coordinates" of the bin which is most likely to contain the
1260 // point
1261 Vector<unsigned> bin_index_v(dim);
1262 coords_to_vectorial_bin_index(zeta, bin_index_v);
1263
1264 // We loop over all the dimensions to find the maximum radius we have to
1265 // do the spiraling if we want to be sure there is no bin left at the end
1266 unsigned max_radius = 0;
1267 for (unsigned k = 0; k < dim; k++)
1268 {
1269 unsigned local =
1270 std::max((bin_index_v[k] + 1),
1271 (Dimensions_of_bin_array[k] - bin_index_v[k] - 1));
1272 if (local > max_radius)
1273 {
1274 max_radius = local;
1275 }
1276 }
1277
1278 // Vector which will store the indices of the neighbouring bins
1279 // at the current radius
1280 Vector<unsigned> bin_index_at_current_radius_old;
1281 Vector<unsigned> bin_index_at_current_radius_new;
1282 while (radius <= max_radius)
1283 {
1284 // Get the neighbouring bins
1285 bin_index_at_current_radius_old.clear();
1286 double t_start = TimingHelpers::timer();
1288 bin_index, radius, bin_index_at_current_radius_old, true);
1289 unsigned nbin_at_current_radius_old =
1290 bin_index_at_current_radius_old.size();
1291 double t_end = TimingHelpers::timer();
1292 double t_old = t_end - t_start;
1293
1294 bin_index_at_current_radius_new.clear();
1295 double t_start_new = TimingHelpers::timer();
1297 bin_index, radius, bin_index_at_current_radius_new, false);
1298 unsigned nbin_at_current_radius_new =
1299 bin_index_at_current_radius_new.size();
1300 double t_end_new = TimingHelpers::timer();
1301
1302 double t_new = t_end_new - t_start_new;
1303
1304 if (nbin_at_current_radius_new != nbin_at_current_radius_old)
1305 {
1306 oomph_info << "Number of bins don't match: new = "
1307 << nbin_at_current_radius_new
1308 << "old = " << nbin_at_current_radius_old
1309 << " radius = " << radius << std::endl;
1310 oomph_info << "Old: " << std::endl;
1311 for (unsigned i = 0; i < nbin_at_current_radius_old; i++)
1312 {
1313 oomph_info << bin_index_at_current_radius_old[i] << " ";
1314 }
1315 oomph_info << std::endl;
1316 oomph_info << "New: " << std::endl;
1317 for (unsigned i = 0; i < nbin_at_current_radius_new; i++)
1318 {
1319 oomph_info << bin_index_at_current_radius_new[i] << " ";
1320 }
1321 oomph_info << std::endl;
1322 }
1323
1324 t_total_new += t_new;
1325 t_total_old += t_old;
1326
1327 if (t_new < t_old)
1328 {
1329 new_faster++;
1330 }
1331 else
1332 {
1333 old_faster++;
1334 }
1335 radius++;
1336 }
1337
1338
1339 oomph_info << "Number of times old/new version was faster: " << old_faster
1340 << " " << new_faster << std::endl
1341 << "Total old/new time: " << t_total_old << " " << t_total_new
1342 << " " << std::endl;
1343 }
1344
1345
1346 //==============================================================================
1347 /// Helper function for computing the bin indices of neighbouring bins
1348 /// at a given "radius" of the specified bin
1349 //==============================================================================
1351 const unsigned& bin_index,
1352 const unsigned& radius,
1353 Vector<unsigned>& neighbouring_bin_index,
1354 const bool& use_old_version)
1355 {
1356 // OLD VERSION
1357 if (use_old_version)
1358 {
1359 neighbouring_bin_index.clear();
1360
1361 // Quick return (slightly hacky -- the machinery below adds the
1362 // "home" bin twice...)
1363 if (radius == 0)
1364 {
1365 neighbouring_bin_index.push_back(bin_index);
1366 return;
1367 }
1368
1369 // translate old code
1370 unsigned level = radius;
1371 unsigned bin = bin_index;
1372
1373 // Dimension
1374 const unsigned n_lagrangian = this->ndim_zeta();
1375
1376 // This will be different depending on the number of Lagrangian
1377 // coordinates
1378 if (n_lagrangian == 1)
1379 {
1380 // Reserve memory for the container where we return the indices
1381 // of the neighbouring bins (2 bins max, left and right)
1382 neighbouring_bin_index.reserve(2);
1383
1384 // Single "loop" in one direction - always a vector of max size 2
1385 unsigned nbr_bin_left = bin - level;
1386 if (nbr_bin_left < Dimensions_of_bin_array[0])
1387 {
1388 unsigned nbr_bin = nbr_bin_left;
1389 neighbouring_bin_index.push_back(nbr_bin);
1390 }
1391 unsigned nbr_bin_right = bin + level;
1392 if ((nbr_bin_right < Dimensions_of_bin_array[0]) &&
1393 (nbr_bin_right != nbr_bin_left))
1394 {
1395 unsigned nbr_bin = nbr_bin_right;
1396 neighbouring_bin_index.push_back(nbr_bin);
1397 }
1398 }
1399 else if (n_lagrangian == 2)
1400 {
1401 // Reserve memory for the container where we return the indices
1402 // of the neighbouring bins
1403 const unsigned n_max_neighbour_bins = 8 * level;
1404 neighbouring_bin_index.reserve(n_max_neighbour_bins);
1405
1406 const unsigned n_total_bin =
1408
1409 // Which row of the bin structure is the current bin on?
1410 // This is just given by the integer answer of dividing bin
1411 // by Nbin_x (the number of bins in a single row)
1412 // e.g. in a 6x6 grid, bins 6 through 11 would all have bin_row=1
1413 const unsigned bin_row = bin / Dimensions_of_bin_array[0];
1414
1415 // The neighbour_bin vector contains all bin numbers at the
1416 // specified "distance" (level) away from the current bin
1417
1418 // Row/column length
1419 const unsigned n_length = (level * 2) + 1;
1420
1421 {
1422 // Visit all the bins at the specified distance (level) away
1423 // from the current bin. In order to obtain the same order in
1424 // the visited bins as the previous algorithm we visit all the
1425 // bins at the specified distance (level) as follows:
1426
1427 // Suppose we want the bins at distance (level=2) of the
1428 // specified bin, then we visit them as indicated below
1429
1430 // 01 02 03 04 05 // First row
1431 // 06 07
1432 // 08 B 09
1433 // 10 11
1434 // 12 13 14 15 16 // Last row
1435 // ^--------------- First column
1436 // ^-- Last column
1437
1438 // ----------------------------------------------------------------
1439 // Visit all the bins in the first row at the specified
1440 // distance (level) away from the current bin
1441
1442 // ------------------ FIRST ROW ------------------------
1443 // Pre-compute the distance in the j-direction
1444 const unsigned j_precomputed = level * Dimensions_of_bin_array[0];
1445 // Pre-compute the row where the bin should lie on
1446 const unsigned j_initial_row = bin_row - level;
1447
1448 // Loop over the columns (of the first row)
1449 for (unsigned i = 0; i < n_length; i++)
1450 {
1451 // --------- First row ------------------------------------------
1452 const unsigned initial_neighbour_bin =
1453 bin - (level - i) - j_precomputed;
1454 // This number might fall on the wrong row of the bin
1455 // structure; this needs to be tested? Not sure why, but leave
1456 // the test!
1457
1458 // Which row is this number on? (see above)
1459 const unsigned initial_neighbour_bin_row =
1460 initial_neighbour_bin / Dimensions_of_bin_array[0];
1461 // These numbers for the rows must match; if it is then add
1462 // initial_neighbour_bin to the neighbour scheme (The bin
1463 // number must also be greater than zero and less than the
1464 // total number of bins)
1465 if ((j_initial_row == initial_neighbour_bin_row) &&
1466 (initial_neighbour_bin < n_total_bin))
1467 {
1468 neighbouring_bin_index.push_back(initial_neighbour_bin);
1469 }
1470
1471 } // for (unsigned i=0;i<n_length;i++)
1472
1473 // Then visit all the bins in the first and last column at the
1474 // specified distance (level) away from the current bin
1475
1476 // ------------------ FIRST AND LAST COLUMNS ---------------------
1477 // Loop over the rows (of the first and last column)
1478 for (unsigned j = 1; j < n_length - 1; j++)
1479 {
1480 // --------- First column ---------------------------------------
1481 const unsigned initial_neighbour_bin =
1482 bin - (level) - ((level - j) * Dimensions_of_bin_array[0]);
1483 // This number might fall on the wrong row of the bin
1484 // structure; this needs to be tested? Not sure why, but leave
1485 // the test!
1486
1487 // Which row is this number on? (see above)
1488 const unsigned initial_neighbour_bin_row =
1489 initial_neighbour_bin / Dimensions_of_bin_array[0];
1490
1491 // Which row should it be on?
1492 const unsigned initial_row = bin_row - (level - j);
1493
1494 // These numbers for the rows must match; if it is then add
1495 // initial_neighbour_bin to the neighbour scheme (The bin
1496 // number must also be greater than zero and less than the
1497 // total number of bins)
1498 if ((initial_row == initial_neighbour_bin_row) &&
1499 (initial_neighbour_bin < n_total_bin))
1500 {
1501 neighbouring_bin_index.push_back(initial_neighbour_bin);
1502 }
1503
1504 // --------- Last column -----------------------------------------
1505 const unsigned final_neighbour_bin =
1506 bin + (level) - ((level - j) * Dimensions_of_bin_array[0]);
1507 // This number might fall on the wrong row of the bin
1508 // structure; this needs to be tested? Not sure why, but leave
1509 // the test!
1510
1511 // Which row is this number on? (see above)
1512 const unsigned final_neighbour_bin_row =
1513 final_neighbour_bin / Dimensions_of_bin_array[0];
1514
1515 // Which row should it be on?
1516 const unsigned final_row = bin_row - (level - j);
1517
1518 // These numbers for the rows must match; if it is then add
1519 // initial_neighbour_bin to the neighbour scheme (The bin
1520 // number must also be greater than zero and less than the
1521 // total number of bins)
1522 if ((final_row == final_neighbour_bin_row) &&
1523 (final_neighbour_bin < n_total_bin))
1524 {
1525 neighbouring_bin_index.push_back(final_neighbour_bin);
1526 }
1527
1528 } // for (unsigned j=1;j<n_length-1;j++)
1529
1530 // ------------------ LAST ROW ------------------------
1531 // Pre-compute the row where the bin should lie on
1532 const unsigned j_final_row = bin_row + level;
1533
1534 // Loop over the columns (of the last row)
1535 for (unsigned i = 0; i < n_length; i++)
1536 {
1537 // --------- Last row ------------------------------------------
1538 const unsigned final_neighbour_bin =
1539 bin - (level - i) + j_precomputed;
1540 // This number might fall on the wrong row of the bin
1541 // structure; this needs to be tested? Not sure why, but leave
1542 // the test!
1543
1544 // Which row is this number on? (see above)
1545 const unsigned final_neighbour_bin_row =
1546 final_neighbour_bin / Dimensions_of_bin_array[0];
1547 // These numbers for the rows must match; if it is then add
1548 // initial_neighbour_bin to the neighbour scheme (The bin
1549 // number must also be greater than zero and less than the
1550 // total number of bins)
1551 if ((j_final_row == final_neighbour_bin_row) &&
1552 (final_neighbour_bin < n_total_bin))
1553 {
1554 neighbouring_bin_index.push_back(final_neighbour_bin);
1555 }
1556
1557 } // for (unsigned i=0;i<n_length;i++)
1558 }
1559 }
1560 else if (n_lagrangian == 3)
1561 {
1562 // Reserve memory for the container where we return the indices
1563 // of the neighbouring bins
1564 const unsigned n_max_neighbour_bins =
1565 8 * level * (3 + 2 * (level - 1)) +
1566 2 * (2 * (level - 1) + 1) * (2 * (level - 1) + 1);
1567 neighbouring_bin_index.reserve(n_max_neighbour_bins);
1568
1569 unsigned n_total_bin = Dimensions_of_bin_array[0] *
1572
1573 // Which layer of the bin structure is the current bin on?
1574 // This is just given by the integer answer of dividing bin
1575 // by Nbin_x*Nbin_y (the number of bins in a single layer
1576 // e.g. in a 6x6x6 grid, bins 72 through 107 would all have bin_layer=2
1577 unsigned bin_layer =
1579
1580 // Which row in this layer is the bin number on?
1581 unsigned bin_row = (bin / Dimensions_of_bin_array[0]) -
1582 (bin_layer * Dimensions_of_bin_array[1]);
1583
1584 // The neighbour_bin vector contains all bin numbers at the
1585 // specified "distance" (level) away from the current bin
1586
1587 // Row/column/layer length
1588 unsigned n_length = (level * 2) + 1;
1589
1590 // Loop over the layers
1591 for (unsigned k = 0; k < n_length; k++)
1592 {
1593 // Loop over the rows
1594 for (unsigned j = 0; j < n_length; j++)
1595 {
1596 // Loop over the columns
1597 for (unsigned i = 0; i < n_length; i++)
1598 {
1599 // Only do this for the end points of every row/layer/column
1600 if ((k == 0) || (k == n_length - 1) || (j == 0) ||
1601 (j == n_length - 1) || (i == 0) || (i == n_length - 1))
1602 {
1603 unsigned nbr_bin = bin - level + i -
1604 ((level - j) * Dimensions_of_bin_array[0]) -
1605 ((level - k) * Dimensions_of_bin_array[0] *
1607 // This number might fall on the wrong
1608 // row or layer of the bin structure; this needs to be tested
1609
1610 // Which layer is this number on?
1611 unsigned nbr_bin_layer = nbr_bin / (Dimensions_of_bin_array[0] *
1613
1614 // Which row is this number on? (see above)
1615 unsigned nbr_bin_row =
1616 (nbr_bin / Dimensions_of_bin_array[0]) -
1617 (nbr_bin_layer * Dimensions_of_bin_array[1]);
1618
1619 // Which layer and row should it be on, given level?
1620 unsigned layer = bin_layer - level + k;
1621 unsigned row = bin_row - level + j;
1622
1623 // These layers and rows must match up:
1624 // if so then add nbr_bin to the neighbour schemes
1625 // (The bin number must also be greater than zero
1626 // and less than the total number of bins)
1627 if ((row == nbr_bin_row) && (layer == nbr_bin_layer) &&
1628 (nbr_bin < n_total_bin))
1629 {
1630 neighbouring_bin_index.push_back(nbr_bin);
1631 }
1632 }
1633 }
1634 }
1635 }
1636 }
1637 }
1638 // LOUIS'S VERSION
1639 else
1640 {
1641 neighbouring_bin_index.clear();
1642
1643 // Just testing if the radius is equal to 0
1644 if (radius == 0)
1645 {
1646 // in this case we only have to push back the bin
1647 neighbouring_bin_index.push_back(bin_index);
1648 }
1649 // Else, if the radius is different from 0
1650 else
1651 {
1652 unsigned dim = ndim_zeta();
1653
1654 // The vector which will store the coordinates of the bin in the bin
1655 // arrray
1656 Vector<int> vector_of_positions(3);
1657
1658 // Will store locally the dimension of the bin array
1659 Vector<int> vector_of_dimensions(3);
1660
1661 // Will store the radiuses for each dimension
1662 // If it is 0, that means the dimension is not vector_of_active_dim
1663 Vector<int> vector_of_radiuses(3);
1664
1665 // Stores true if the dimension is "active" or false if the
1666 // dimension is inactive (for example the third dimension is inactive in
1667 // a 2D mesh).
1668 std::vector<bool> vector_of_active_dim(3);
1669
1670 // Will store the coefficients you have to multiply each bin coordinate
1671 // to have the correct unsigned corresponding to the bin index.
1672 Vector<int> vector_of_coef(3);
1673
1674
1675 // First initializing this MyArrays for the active dimensions
1676 unsigned coef = 1;
1677 for (unsigned u = 0; u < dim; u++)
1678 {
1679 vector_of_positions[u] =
1680 (((bin_index / coef) % (Dimensions_of_bin_array[u])));
1681 vector_of_coef[u] = coef;
1682 coef *= Dimensions_of_bin_array[u];
1683 vector_of_dimensions[u] = Dimensions_of_bin_array[u];
1684 vector_of_radiuses[u] = radius;
1685 vector_of_active_dim[u] = true;
1686 }
1687 // Filling the rest with default values
1688 for (unsigned u = dim; u < 3; u++)
1689 {
1690 vector_of_positions[u] = 0;
1691 vector_of_coef[u] = 0;
1692 vector_of_dimensions[u] = 1;
1693 vector_of_radiuses[u] = 0;
1694 vector_of_active_dim[u] = false;
1695 }
1696
1697 // First we fill the bins which corresponds to x+radius and x-radius in
1698 // the bin array (x being the first coordinate/dimension).
1699 // For j equal to radius or -radius
1700 for (int j = -vector_of_radiuses[0]; j <= vector_of_radiuses[0];
1701 j += 2 * vector_of_radiuses[0]) // j corresponds to x
1702 {
1703 int local_tempj =
1704 vector_of_positions[0] + j; // Corresponding bin index
1705 // if we are in the bin array
1706 if (local_tempj >= 0 && local_tempj < vector_of_dimensions[0])
1707 {
1708 // Loop over all the bins
1709 for (int i = -vector_of_radiuses[1]; i <= vector_of_radiuses[1];
1710 i++)
1711 {
1712 // i corresponds to y (2nd dim)
1713 int local_tempi = vector_of_positions[1] + i;
1714 // if we are still in the bin array
1715 if (local_tempi >= 0 && local_tempi < vector_of_dimensions[1])
1716 {
1717 for (int k = -vector_of_radiuses[2]; k <= vector_of_radiuses[2];
1718 k++)
1719 {
1720 // k corresponds to z (third dimension)
1721 int local_tempk = vector_of_positions[2] + k;
1722 // if we are still in the bin array
1723 if (local_tempk >= 0 && local_tempk < vector_of_dimensions[2])
1724 {
1725 neighbouring_bin_index.push_back(
1726 local_tempj * vector_of_coef[0] +
1727 local_tempi * vector_of_coef[1] +
1728 local_tempk * vector_of_coef[2]);
1729 }
1730 }
1731 }
1732 }
1733 }
1734 }
1735 // Secondly we get the bins corresponding to y+radius and y-radius
1736 // only if the second dimension is active
1737 if (vector_of_active_dim[1])
1738 {
1739 // For i equal to radius or -radius (i corresponds to y)
1740 for (int i = -vector_of_radiuses[1]; i <= vector_of_radiuses[1];
1741 i += 2 * vector_of_radiuses[1])
1742 {
1743 int local_tempi = vector_of_positions[1] + i;
1744 if (local_tempi >= 0 && local_tempi < vector_of_dimensions[1])
1745 {
1746 // We loop over all the surface
1747 for (int j = -vector_of_radiuses[0] + 1;
1748 j <= vector_of_radiuses[0] - 1;
1749 j++)
1750 {
1751 int local_tempj = vector_of_positions[0] + j;
1752 if (local_tempj >= 0 && local_tempj < vector_of_dimensions[0])
1753 {
1754 for (int k = -vector_of_radiuses[2];
1755 k <= vector_of_radiuses[2];
1756 k++)
1757 {
1758 int local_tempk = vector_of_positions[2] + k;
1759 if (local_tempk >= 0 &&
1760 local_tempk < vector_of_dimensions[2])
1761 {
1762 neighbouring_bin_index.push_back(
1763 local_tempj * vector_of_coef[0] +
1764 local_tempi * vector_of_coef[1] +
1765 local_tempk * vector_of_coef[2]);
1766 }
1767 }
1768 }
1769 }
1770 }
1771 }
1772 }
1773 // Thirdly we get the bins corresponding to z+radius and z-radius
1774 // if the third dimension is active.
1775 if (vector_of_active_dim[2])
1776 {
1777 // for k equal to radius or -radius (k corresponds to z)
1778 for (int k = -vector_of_radiuses[2]; k <= vector_of_radiuses[2];
1779 k += 2 * vector_of_radiuses[2])
1780 {
1781 int local_tempk = vector_of_positions[2] + k;
1782 if (local_tempk >= 0 && local_tempk < vector_of_dimensions[2])
1783 {
1784 // We loop over all the surface
1785 for (int j = -vector_of_radiuses[0] + 1;
1786 j <= vector_of_radiuses[0] - 1;
1787 j++)
1788 {
1789 int local_tempj = vector_of_positions[0] + j;
1790 if (local_tempj >= 0 && local_tempj < vector_of_dimensions[0])
1791 {
1792 for (int i = -vector_of_radiuses[1] + 1;
1793 i <= vector_of_radiuses[1] - 1;
1794 i++)
1795 {
1796 int local_tempi = vector_of_positions[1] + i;
1797 if (local_tempi >= 0 &&
1798 local_tempi < vector_of_dimensions[1])
1799 {
1800 neighbouring_bin_index.push_back(
1801 local_tempj * vector_of_coef[0] +
1802 local_tempi * vector_of_coef[1] +
1803 local_tempk * vector_of_coef[2]);
1804 }
1805 }
1806 }
1807 }
1808 }
1809 }
1810 }
1811 }
1812 } // end new version
1813 }
1814
1815
1816 //==============================================================================
1817 /// Compute total number of sample points recursively
1818 //==============================================================================
1821 {
1822 unsigned count = 0;
1823 unsigned n_bin = nbin();
1824 for (unsigned i = 0; i < n_bin; i++)
1825 {
1826 if (Bin_pt[i] != 0)
1827 {
1828 count +=
1829 Bin_pt[i]->total_number_of_sample_points_computed_recursively();
1830 }
1831 }
1832 return count;
1833 }
1834
1835
1836 //============================================================
1837 /// Get (linearly enumerated) bin index of bin that
1838 /// contains specified zeta
1839 //============================================================
1840 unsigned BinArray::coords_to_bin_index(const Vector<double>& zeta)
1841 {
1842 unsigned coef = 1;
1843 unsigned n_bin = 0;
1844 const unsigned dim = ndim_zeta();
1845
1846 // Loop over all the dimensions
1847 for (unsigned u = 0; u < dim; u++)
1848 {
1849 // for each one get the correct bin "coordinate"
1850 unsigned local_bin_number = 0;
1851
1852 if (zeta[u] < Min_and_max_coordinates[u].first)
1853 {
1854 local_bin_number = 0;
1855 }
1856 else if (zeta[u] > Min_and_max_coordinates[u].second)
1857 {
1858 local_bin_number = dimensions_of_bin_array(u) - 1;
1859 }
1860 else
1861 {
1862 local_bin_number = (int(std::min(
1864 unsigned(floor(double(zeta[u] - Min_and_max_coordinates[u].first) /
1865 double(Min_and_max_coordinates[u].second -
1866 Min_and_max_coordinates[u].first) *
1867 double(dimensions_of_bin_array(u)))))));
1868 }
1869
1870 /// Update the coef and the bin index
1871 n_bin += local_bin_number * coef;
1872 coef *= dimensions_of_bin_array(u);
1873 }
1874
1875 // return the correct bin index
1876 return (n_bin);
1877 }
1878
1879 //==========================================================================
1880 /// Fill the bin array with sample points from FiniteElements stored in mesh
1881 //==========================================================================
1883 {
1884 // Fill 'em in:
1885 unsigned nel = Mesh_pt->nelement();
1886 for (unsigned e = 0; e < nel; e++)
1887 {
1888 FiniteElement* el_pt = Mesh_pt->finite_element_pt(e);
1889
1890 // Total number of sample point we will create
1891 unsigned nplot =
1892 el_pt->nplot_points(Nsample_points_generated_per_element);
1893
1894 /// For all the sample points we have to create ...
1895 for (unsigned j = 0; j < nplot; j++)
1896 {
1897 // ... create it: Pass element index in mesh (vector
1898 // of elements and index of sample point within element
1899 SamplePoint* new_sample_point_pt = new SamplePoint(e, j);
1900
1901 // Coordinates of this point
1902 Vector<double> zeta(ndim_zeta());
1903 Vector<double> s(ndim_zeta());
1904 bool use_equally_spaced_interior_sample_points =
1906 el_pt->get_s_plot(j,
1908 s,
1909 use_equally_spaced_interior_sample_points);
1911 {
1912 el_pt->interpolated_x(s, zeta);
1913 }
1914 else
1915 {
1916 el_pt->interpolated_zeta(s, zeta);
1917 }
1918
1919#ifdef PARANOID
1920
1921 // Check if point is inside
1922 bool is_inside = true;
1923 std::ostringstream error_message;
1924 unsigned dim = ndim_zeta();
1925 for (unsigned i = 0; i < dim; i++)
1926 {
1927 if ((zeta[i] < Min_and_max_coordinates[i].first) ||
1928 (zeta[i] > Min_and_max_coordinates[i].second))
1929 {
1930 is_inside = false;
1931 error_message << "Sample point at zeta[" << i << "] = " << zeta[i]
1932 << " is outside limits of bin array: "
1933 << Min_and_max_coordinates[i].first << " and "
1934 << Min_and_max_coordinates[i].second << std::endl;
1935 }
1936 }
1937
1938 if (!is_inside)
1939 {
1940 error_message << "Please correct the limits passed to the "
1941 << "constructor." << std::endl;
1942 throw OomphLibError(error_message.str(),
1943 OOMPH_CURRENT_FUNCTION,
1944 OOMPH_EXCEPTION_LOCATION);
1945 }
1946
1947#endif
1948
1949
1950 // Finding the correct bin to put the sample point
1951 unsigned bin_index = coords_to_bin_index(zeta);
1952
1953 // if the bin is not yet created, create it
1954 if (Bin_pt[bin_index] == 0)
1955 {
1956 Bin_pt[bin_index] = new RefineableBin(this, bin_index);
1957 }
1958
1959 // ... and then fill the bin with this new sample point
1960 Bin_pt[bin_index]->add_sample_point(new_sample_point_pt, zeta);
1961 }
1962 }
1963 }
1964
1965
1966 //==================================================================
1967 /// Get "coordinates" of bin that contains specified zeta
1968 //==================================================================
1969 void BinArray::coords_to_vectorial_bin_index(const Vector<double>& zeta,
1970 Vector<unsigned>& bin_index)
1971 {
1972 unsigned dim = ndim_zeta();
1973 bin_index.resize(dim);
1974 for (unsigned u = 0; u < dim; u++)
1975 {
1976 if (zeta[u] < Min_and_max_coordinates[u].first)
1977 {
1978 bin_index[u] = 0;
1979 }
1980 else if (zeta[u] > Min_and_max_coordinates[u].second)
1981 {
1982 bin_index[u] = dimensions_of_bin_array(u) - 1;
1983 }
1984 else
1985 {
1986 bin_index[u] = (int(
1987 std::min(dimensions_of_bin_array(u) - 1,
1988 unsigned(floor((zeta[u] - Min_and_max_coordinates[u].first) /
1989 double(Min_and_max_coordinates[u].second -
1990 Min_and_max_coordinates[u].first) *
1991 double(dimensions_of_bin_array(u)))))));
1992 }
1993 }
1994 }
1995
1996
1997 //==============================================================================
1998 /// Find sub-GeomObject (finite element) and the local coordinate
1999 /// s within it that contains point with global coordinate zeta.
2000 /// sub_geom_object_pt=0 if point can't be found.
2001 //==============================================================================
2002 void RefineableBinArray::locate_zeta(const Vector<double>& zeta,
2003 GeomObject*& sub_geom_object_pt,
2004 Vector<double>& s)
2005 {
2006 // Default: we've failed miserably
2007 sub_geom_object_pt = 0;
2008
2009 unsigned dim = ndim_zeta();
2010
2011 // Top level book keeping and sanity checking
2012 if (Depth == 0)
2013 {
2014 // Reset counter for number of sample points visited.
2015 // If we can't find the point we should at least make sure that
2016 // we've visited all the sample points before giving up.
2018 0;
2019
2020 // Does the zeta coordinate lie within the current (top level!) bin
2021 // structure? Skip this test if we want to always fail because that's
2022 // usually done to trace out the spiral path
2024 {
2025 // Loop over the lagrangian dimension
2026 for (unsigned i = 0; i < dim; i++)
2027 {
2028 // If the i-th coordinate is less than the minimum
2029 if (zeta[i] < Min_and_max_coordinates[i].first)
2030 {
2031 return;
2032 }
2033 // Otherwise coordinate may be bigger than the maximum
2034 else if (zeta[i] > Min_and_max_coordinates[i].second)
2035 {
2036 return;
2037 }
2038 }
2039 }
2040 }
2041
2042 // Finding the bin in which the point is located
2043 int bin_index = coords_to_bin_index(zeta);
2044
2045#ifdef PARANOID
2046 if (bin_index < 0)
2047 {
2048 throw OomphLibError("Negative bin index...",
2049 OOMPH_CURRENT_FUNCTION,
2050 OOMPH_EXCEPTION_LOCATION);
2051 }
2052#endif
2053
2054 // Start at this radius (radius = 0 is the central bin)
2055 unsigned radius = 0;
2056
2057 // "coordinates" of the bin which is most likely to contain the
2058 // point
2059 Vector<unsigned> bin_index_v(dim);
2060 coords_to_vectorial_bin_index(zeta, bin_index_v);
2061
2062 // We loop over all the dimensions to find the maximum radius we have to
2063 // do the spiraling if we want to be sure there is no bin left at the end
2064 unsigned max_radius = 0;
2065 for (unsigned k = 0; k < dim; k++)
2066 {
2067 unsigned local =
2068 std::max((bin_index_v[k] + 1),
2069 (Dimensions_of_bin_array[k] - bin_index_v[k] - 1));
2070 if (local > max_radius)
2071 {
2072 max_radius = local;
2073 }
2074 }
2075
2076 // Vector which will store the indices of the neighbouring bins
2077 // at the current radius
2078 Vector<unsigned> bin_index_at_current_radius;
2079 while (radius <= max_radius)
2080 {
2081 // Get the neighbouring bins
2082 bin_index_at_current_radius.clear();
2084 bin_index, radius, bin_index_at_current_radius);
2085 // How many are there
2086 unsigned nbin_at_current_radius = bin_index_at_current_radius.size();
2087 unsigned n_bin = nbin();
2088
2089 // Keep looping over entries (stop if we found geom object)
2090 unsigned k = 0;
2091 while ((k < nbin_at_current_radius) && (sub_geom_object_pt == 0))
2092 {
2093 int neigh_bin_index = bin_index_at_current_radius[k];
2094#ifdef PARANOID
2095 if (neigh_bin_index < 0)
2096 {
2097 throw OomphLibError("Negative neighbour bin index...",
2098 OOMPH_CURRENT_FUNCTION,
2099 OOMPH_EXCEPTION_LOCATION);
2100 }
2101#endif
2102 if (neigh_bin_index < int(n_bin))
2103 {
2104 // If the bin exists
2105 if (Bin_pt[neigh_bin_index] != 0)
2106 {
2107 // We call the correct method to locate_zeta in the bin
2108 Bin_pt[neigh_bin_index]->locate_zeta(zeta, sub_geom_object_pt, s);
2109 }
2110 }
2111 k++;
2112 }
2113
2114 // Reached end of loop over bins at this radius (or found the point)
2115 // Which one?
2116 if (sub_geom_object_pt != 0)
2117 {
2118 return;
2119 }
2120 // Increment radius
2121 else
2122 {
2123 radius++;
2124 }
2125 }
2126
2127
2128#ifdef PARANOID
2129
2130 // If we still haven't found the point, check that we've at least visited
2131 // all the sample points
2132 if (Depth == 0)
2133 {
2134 if (
2137 {
2138 if (max_search_radius() == DBL_MAX)
2139 {
2140 std::ostringstream error_message;
2141 error_message
2142 << "Zeta not found after visiting "
2144 << " sample points out of "
2146 << "Where are the missing sample points???\n";
2147 throw OomphLibError(error_message.str(),
2148 OOMPH_CURRENT_FUNCTION,
2149 OOMPH_EXCEPTION_LOCATION);
2150 }
2151 }
2152 }
2153
2154#endif
2155 }
2156
2157 /// Default number of bins (in each coordinate direction)
2159
2160
2161 /// /////////////////////////////////////////////////////////////////////////////
2162 /// /////////////////////////////////////////////////////////////////////////////
2163 /// /////////////////////////////////////////////////////////////////////////////
2164
2165
2166 //======================================================================
2167 /// Constructor
2168 //======================================================================
2170 SamplePointContainerParameters* sample_point_container_parameters_pt)
2172 sample_point_container_parameters_pt->mesh_pt(),
2173 sample_point_container_parameters_pt->min_and_max_coordinates(),
2174 sample_point_container_parameters_pt
2175 ->use_eulerian_coordinates_during_setup(),
2176 sample_point_container_parameters_pt
2177 ->ignore_halo_elements_during_locate_zeta_search(),
2178 sample_point_container_parameters_pt
2179 ->nsample_points_generated_per_element()),
2180 BinArray(
2181 sample_point_container_parameters_pt->mesh_pt(),
2182 sample_point_container_parameters_pt->min_and_max_coordinates(),
2183 dynamic_cast<BinArrayParameters*>(sample_point_container_parameters_pt)
2184 ->dimensions_of_bin_array(),
2185 sample_point_container_parameters_pt
2186 ->use_eulerian_coordinates_during_setup(),
2187 sample_point_container_parameters_pt
2188 ->ignore_halo_elements_during_locate_zeta_search(),
2189 sample_point_container_parameters_pt
2190 ->nsample_points_generated_per_element())
2191 {
2192 // Set default size of bin array (and spatial dimension!)
2193 if (Dimensions_of_bin_array.size() == 0)
2194 {
2195 int dim = 0;
2196 if (Mesh_pt->nelement() != 0)
2197 {
2198 dim = Mesh_pt->finite_element_pt(0)->dim();
2199 }
2200
2201
2202 // Need to do an Allreduce to ensure that the dimension is consistent
2203 // even when no elements are assigned to a certain processor
2204#ifdef OOMPH_HAS_MPI
2205 // Only a problem if the mesh has been distributed
2206 if (Mesh_pt->is_mesh_distributed())
2207 {
2208 // Need a non-null communicator
2209 if (Mesh_pt->communicator_pt() != 0)
2210 {
2211 int n_proc = Mesh_pt->communicator_pt()->nproc();
2212 if (n_proc > 1)
2213 {
2214 int dim_reduce;
2215 MPI_Allreduce(&dim,
2216 &dim_reduce,
2217 1,
2218 MPI_INT,
2219 MPI_MAX,
2220 Mesh_pt->communicator_pt()->mpi_comm());
2221 dim = dim_reduce;
2222 }
2223 }
2224 }
2225#endif
2226
2228 }
2229
2230 // Have we specified max/min coordinates?
2231 // If not, compute them on the fly from mesh
2232 if (Min_and_max_coordinates.size() == 0)
2233 {
2235 }
2236
2237 // Spiraling parameters
2238 Max_spiral_level = UINT_MAX;
2240
2241 NonRefineableBinArrayParameters* non_ref_bin_array_parameters_pt =
2242 dynamic_cast<NonRefineableBinArrayParameters*>(
2243 sample_point_container_parameters_pt);
2244
2245#ifdef PARANOID
2246 if (non_ref_bin_array_parameters_pt == 0)
2247 {
2248 throw OomphLibError("Wrong sample_point_container_parameters_pt",
2249 OOMPH_CURRENT_FUNCTION,
2250 OOMPH_EXCEPTION_LOCATION);
2251 }
2252#endif
2253
2254 Nspiral_chunk = non_ref_bin_array_parameters_pt->nspiral_chunk();
2256
2257 // Time it
2258 double t_start = 0.0;
2260 {
2261 t_start = TimingHelpers::timer();
2262 }
2263
2264 // Now fill the bastard...
2266
2268 {
2269 double t_end = TimingHelpers::timer();
2271 oomph_info << "Time for setup of " << Dimensions_of_bin_array.size()
2272 << "-dimensional sample point container containing " << npts
2273 << " sample points: " << t_end - t_start
2274 << " sec (non-ref_bin); third party: 0 sec ( = 0 %)"
2275 << std::endl;
2276 }
2277 }
2278
2279
2280 //==============================================================================
2281 /// Compute total number of sample points recursively
2282 //==============================================================================
2285 {
2286 // Get pointer to map-based representation
2287 const std::map<unsigned, Vector<std::pair<FiniteElement*, Vector<double>>>>*
2288 map_pt = Bin_object_coord_pairs.map_pt();
2289
2290 // Initialise
2291 unsigned count = 0;
2292
2293 // loop...
2294 typedef std::map<
2295 unsigned,
2296 Vector<std::pair<FiniteElement*, Vector<double>>>>::const_iterator IT;
2297 for (IT it = map_pt->begin(); it != map_pt->end(); it++)
2298 {
2299 count += (*it).second.size();
2300 }
2301 return count;
2302 }
2303
2304
2305 //========================================================================
2306 /// Output bins
2307 //========================================================================
2308 void NonRefineableBinArray::output_bins(std::ofstream& outfile)
2309 {
2310 // Spatial dimension of bin
2311 const unsigned n_lagrangian = this->ndim_zeta();
2312
2313 unsigned nbin_x = Dimensions_of_bin_array[0];
2314 unsigned nbin_y = 1;
2315 if (n_lagrangian > 1) nbin_y = Dimensions_of_bin_array[1];
2316 unsigned nbin_z = 1;
2317 if (n_lagrangian > 2) nbin_z = Dimensions_of_bin_array[2];
2318
2319 unsigned b = 0;
2320 for (unsigned iz = 0; iz < nbin_z; iz++)
2321 {
2322 for (unsigned iy = 0; iy < nbin_y; iy++)
2323 {
2324 for (unsigned ix = 0; ix < nbin_x; ix++)
2325 {
2326 unsigned nentry = Bin_object_coord_pairs[b].size();
2327 for (unsigned e = 0; e < nentry; e++)
2328 {
2329 FiniteElement* el_pt = Bin_object_coord_pairs[b][e].first;
2330 Vector<double> s(Bin_object_coord_pairs[b][e].second);
2331 Vector<double> zeta(n_lagrangian);
2333 {
2334 el_pt->interpolated_x(s, zeta);
2335 }
2336 else
2337 {
2338 el_pt->interpolated_zeta(s, zeta);
2339 }
2340 for (unsigned i = 0; i < n_lagrangian; i++)
2341 {
2342 outfile << zeta[i] << " ";
2343 }
2344 outfile << ix << " " << iy << " " << iz << " " << b << " "
2345 << std::endl;
2346 }
2347 b++;
2348 }
2349 }
2350 }
2351 }
2352
2353
2354 //========================================================================
2355 /// Output bin vertices (allowing display of bins as zones).
2356 //========================================================================
2358 {
2359 // Store the lagrangian dimension
2360 const unsigned n_lagrangian = this->ndim_zeta();
2361
2362 unsigned nbin = Dimensions_of_bin_array[0];
2363 if (n_lagrangian > 1) nbin *= Dimensions_of_bin_array[1];
2364 if (n_lagrangian > 2) nbin *= Dimensions_of_bin_array[2];
2365
2366 for (unsigned i_bin = 0; i_bin < nbin; i_bin++)
2367 {
2368 // Get bin vertices
2369 Vector<Vector<double>> bin_vertex;
2370 get_bin_vertices(i_bin, bin_vertex);
2371 switch (n_lagrangian)
2372 {
2373 case 1:
2374 outfile << "ZONE I=2\n";
2375 break;
2376
2377 case 2:
2378 outfile << "ZONE I=2, J=2\n";
2379 break;
2380
2381 case 3:
2382 outfile << "ZONE I=2, J=2, K=2\n";
2383 break;
2384 }
2385
2386 unsigned nvertex = bin_vertex.size();
2387 for (unsigned i = 0; i < nvertex; i++)
2388 {
2389 for (unsigned j = 0; j < n_lagrangian; j++)
2390 {
2391 outfile << bin_vertex[i][j] << " ";
2392 }
2393 outfile << std::endl;
2394 }
2395 }
2396 }
2397
2398
2399 //========================================================================
2400 /// Fill the bin array with sample points from FiniteElements stored in mesh
2401 //========================================================================
2403 {
2404 // Store the lagrangian dimension
2405 const unsigned n_lagrangian = this->ndim_zeta();
2406
2407 // Flush all objects out of the bin structure
2409
2410 // Temporary storage in bin
2411 std::map<unsigned, Vector<std::pair<FiniteElement*, Vector<double>>>>
2412 tmp_bin_object_coord_pairs;
2413
2414 // Total number of bins
2415 unsigned ntotalbin = nbin();
2416
2417 // Initialise the structure that keep tracks of the occupided bins
2418 Bin_object_coord_pairs.initialise(ntotalbin);
2419
2420
2421 // Issue warning about small number of bins
2423 {
2424 // Calculate the (nearest integer) number of elements per bin
2425 unsigned n_mesh_element = Mesh_pt->nelement();
2426 unsigned elements_per_bin = n_mesh_element / ntotalbin;
2427 // If it is above the threshold then issue a warning
2428 if (elements_per_bin > Threshold_for_elements_per_bin_warning)
2429 {
2431 {
2433 std::ostringstream warn_message;
2434 warn_message
2435 << "The average (integer) number of elements per bin is \n\n"
2436 << elements_per_bin << ", which is more than \n\n"
2437 << " "
2438 "NonRefineableBinArray::Threshold_for_elements_per_bin_warning="
2440 << "\n\nIf the lookup seems slow (and you have the memory),"
2441 << "consider increasing\n"
2442 << "BinArray::Dimensions_of_bin_array from their current\n"
2443 << "values of { ";
2444 unsigned nn = Dimensions_of_bin_array.size();
2445 for (unsigned ii = 0; ii < nn; ii++)
2446 {
2447 warn_message << Dimensions_of_bin_array[ii] << " ";
2448 }
2449 warn_message
2450 << " }.\n"
2451 << "\nNOTE: You can suppress this warning by increasing\n\n"
2452 << " NonRefineableBinArray::"
2453 << "Threshold_for_elements_per_bin_warning\n\n"
2454 << "or by setting \n\n "
2455 << "NonRefineableBinArray::Suppress_warning_about_small_number_of_"
2456 "bins\n\n"
2457 << "to true (both are public static data).\n\n";
2458 OomphLibWarning(warn_message.str(),
2459 OOMPH_CURRENT_FUNCTION,
2460 OOMPH_EXCEPTION_LOCATION);
2461 }
2462 }
2463 }
2464
2465
2466 // Increase overall counter
2467 Total_nbin_cells_counter += ntotalbin;
2468
2469
2470 // Issue warning?
2472 {
2475 {
2477 {
2479 std::ostringstream warn_message;
2480 warn_message
2481 << "Have allocated \n\n"
2482 << " NonRefineableBinArray::Total_nbin_cells_counter="
2484 << "\n\nbin cells, which is more than \n\n"
2485 << " NonRefineableBinArray::"
2486 << "Threshold_for_total_bin_cell_number_warning="
2487 << NonRefineableBinArray::
2488 Threshold_for_total_bin_cell_number_warning
2489 << "\n\nIf you run out of memory, consider reducing\n"
2490 << "BinArray::Dimensions_of_bin_array from their current\n"
2491 << "values of { ";
2492 unsigned nn = Dimensions_of_bin_array.size();
2493 for (unsigned ii = 0; ii < nn; ii++)
2494 {
2495 warn_message << Dimensions_of_bin_array[ii] << " ";
2496 }
2497 warn_message
2498 << " }.\n"
2499 << "\nNOTE: You can suppress this warning by increasing\n\n"
2500 << " NonRefineableBinArray::"
2501 << "Threshold_for_total_bin_cell_number_warning\n\n"
2502 << "or by setting \n\n NonRefineableBinArray::"
2503 << "Suppress_warning_about_large_total_number_of_bins\n\n"
2504 << "to true (both are public static data).\n\n"
2505 << "NOTE: I'll only issue this warning once -- total number of\n"
2506 << "bins may grow yet further!\n";
2507 OomphLibWarning(warn_message.str(),
2508 OOMPH_CURRENT_FUNCTION,
2509 OOMPH_EXCEPTION_LOCATION);
2510 }
2511 }
2512 }
2513
2514 /// Loop over subobjects (elements) to decide which bin they belong in...
2515 unsigned n_sub = Mesh_pt->nelement();
2516 for (unsigned e = 0; e < n_sub; e++)
2517 {
2518 // Cast to the element (sub-object) first
2519 FiniteElement* el_pt =
2520 dynamic_cast<FiniteElement*>(Mesh_pt->finite_element_pt(e));
2521
2522 // Get specified number of points within the element
2523 unsigned n_plot_points =
2524 el_pt->nplot_points(Nsample_points_generated_per_element);
2525
2526 for (unsigned iplot = 0; iplot < n_plot_points; iplot++)
2527 {
2528 // Storage for local and global coordinates
2529 Vector<double> local_coord(n_lagrangian, 0.0);
2530 Vector<double> global_coord(n_lagrangian, 0.0);
2531
2532 // Get local coordinate and interpolate to global
2533 bool use_equally_spaced_interior_sample_points =
2535 el_pt->get_s_plot(iplot,
2537 local_coord,
2538 use_equally_spaced_interior_sample_points);
2539
2540 // Now get appropriate global coordinate
2542 {
2543 el_pt->interpolated_x(local_coord, global_coord);
2544 }
2545 else
2546 {
2547 el_pt->interpolated_zeta(local_coord, global_coord);
2548 }
2549
2550 // Which bin are the global coordinates in?
2551 unsigned bin_number = 0;
2552 unsigned multiplier = 1;
2553 // Loop over the dimension
2554 for (unsigned i = 0; i < n_lagrangian; i++)
2555 {
2556#ifdef PARANOID
2557 if ((global_coord[i] < Min_and_max_coordinates[i].first) ||
2558 (global_coord[i] > Min_and_max_coordinates[i].second))
2559 {
2560 std::ostringstream error_message;
2561 error_message
2562 << "Bin sample point " << iplot << " in element " << e << "\n"
2563 << "is outside bin limits in coordinate direction " << i << ":\n"
2564 << "Sample point coordinate: " << global_coord[i] << "\n"
2565 << "Max bin coordinate : "
2566 << Min_and_max_coordinates[i].second << "\n"
2567 << "Min bin coordinate : " << Min_and_max_coordinates[i].first
2568 << "\n"
2569 << "You should either setup the bin boundaries manually\n"
2570 << "or increase the percentage offset by which the\n"
2571 << "automatically computed bin limits are increased \n"
2572 << "beyond their sampled max/mins. This is defined in\n"
2573 << "the (public) namespace member\n\n"
2574 << "SamplePointContainer::Percentage_offset \n\n which \n"
2575 << "currently has the value: "
2577 throw OomphLibError(error_message.str(),
2578 OOMPH_CURRENT_FUNCTION,
2579 OOMPH_EXCEPTION_LOCATION);
2580 }
2581
2582#endif
2583 unsigned bin_number_i =
2585 ((global_coord[i] - Min_and_max_coordinates[i].first) /
2586 (Min_and_max_coordinates[i].second -
2587 Min_and_max_coordinates[i].first)));
2588
2589 // Buffer the case when the global coordinate is the maximum
2590 // value
2591 if (bin_number_i == Dimensions_of_bin_array[i])
2592 {
2593 bin_number_i -= 1;
2594 }
2595
2596 // Add to the bin number
2597 bin_number += multiplier * bin_number_i;
2598
2599 // Sort out the multiplier
2600 multiplier *= Dimensions_of_bin_array[i];
2601 }
2602
2603 // Add element-sample local coord pair to the calculated bin
2604 tmp_bin_object_coord_pairs[bin_number].push_back(
2605 std::make_pair(el_pt, local_coord));
2606 }
2607 }
2608
2609 // Finally copy filled vectors across -- wiping memory from temporary
2610 // map as we go along is good and wouldn't be possible if we
2611 // filled the SparseVector's internal map within that class.
2612 typedef std::map<
2613 unsigned,
2614 Vector<std::pair<FiniteElement*, Vector<double>>>>::iterator IT;
2615 for (IT it = tmp_bin_object_coord_pairs.begin();
2616 it != tmp_bin_object_coord_pairs.end();
2617 it++)
2618 {
2619 Bin_object_coord_pairs.set_value((*it).first, (*it).second);
2620 // Make space immediately
2621 (*it).second.clear();
2622 }
2623 }
2624
2625
2626 //========================================================================
2627 /// Provide some stats on the fill level of the associated bin
2628 //========================================================================
2630 unsigned& max_n_entry,
2631 unsigned& min_n_entry,
2632 unsigned& tot_n_entry,
2633 unsigned& n_empty) const
2634 {
2635 // Total number of bins
2636 n_bin = nbin();
2637 n_empty = n_bin - Bin_object_coord_pairs.nnz();
2638
2639 // Get pointer to map-based representation
2640 const std::map<unsigned, Vector<std::pair<FiniteElement*, Vector<double>>>>*
2641 map_pt = Bin_object_coord_pairs.map_pt();
2642
2643 // Initialise
2644 max_n_entry = 0;
2645 min_n_entry = UINT_MAX;
2646 tot_n_entry = 0;
2647
2648 // Do stats
2649 typedef std::map<
2650 unsigned,
2651 Vector<std::pair<FiniteElement*, Vector<double>>>>::const_iterator IT;
2652 for (IT it = map_pt->begin(); it != map_pt->end(); it++)
2653 {
2654 unsigned nentry = (*it).second.size();
2655 if (nentry > max_n_entry) max_n_entry = nentry;
2656 if (nentry < min_n_entry) min_n_entry = nentry;
2657 tot_n_entry += nentry;
2658 }
2659 }
2660
2661
2662 //========================================================================
2663 /// Fill bin by diffusion, populating each empty bin with the same content
2664 /// as the first non-empty bin found during a spiral-based search
2665 /// up to the specified "radius" (default 1)
2666 //========================================================================
2668 const unsigned& bin_diffusion_radius)
2669 {
2670 // oomph_info << "PROFILING GET_NEIGHBOURING_BINS_HELPER():" << std::endl;
2671 // profile_get_neighbouring_bins_helper();
2672
2673 // Loop over all bins to check if they're empty
2674 std::list<unsigned> empty_bins;
2675 unsigned n_bin = nbin();
2676 std::vector<bool> was_empty_until_current_iteration(n_bin, false);
2677 for (unsigned i = 0; i < n_bin; i++)
2678 {
2679 if (Bin_object_coord_pairs[i].size() == 0)
2680 {
2681 empty_bins.push_front(i);
2682 was_empty_until_current_iteration[i] = true;
2683 }
2684 }
2685
2686 // Now keep processing the empty bins until there are none left
2687 unsigned iter = 0;
2688 Vector<unsigned> newly_filled_bin;
2689 while (empty_bins.size() != 0)
2690 {
2691 newly_filled_bin.clear();
2692 newly_filled_bin.reserve(empty_bins.size());
2693 for (std::list<unsigned>::iterator it = empty_bins.begin();
2694 it != empty_bins.end();
2695 it++)
2696 {
2697 unsigned bin = (*it);
2698
2699 // Look for immediate neighbours within the specified "bin radius"
2700 unsigned level = bin_diffusion_radius;
2701 Vector<unsigned> neighbour_bin;
2702 get_neighbouring_bins_helper(bin, level, neighbour_bin);
2703 unsigned n_neigh = neighbour_bin.size();
2704
2705 // Find closest pair
2706 double min_dist = DBL_MAX;
2707 std::pair<FiniteElement*, Vector<double>> closest_pair;
2708 for (unsigned i = 0; i < n_neigh; i++)
2709 {
2710 unsigned neigh_bin = neighbour_bin[i];
2711
2712 // Only allow to re-populate from bins that were already filled at
2713 // previous iteration, otherwise things can progate too fast
2714 if (!was_empty_until_current_iteration[neigh_bin])
2715 {
2716 unsigned nbin_content = Bin_object_coord_pairs[neigh_bin].size();
2717 for (unsigned j = 0; j < nbin_content; j++)
2718 {
2719 FiniteElement* el_pt = Bin_object_coord_pairs[neigh_bin][j].first;
2720 Vector<double> s(Bin_object_coord_pairs[neigh_bin][j].second);
2721 Vector<double> x(2);
2722 el_pt->interpolated_x(s, x);
2723 // Get minimum distance of sample point from any of the vertices
2724 // of current bin
2725 // hierher Louis questions if this is the right thing to do!
2726 // Matthias agrees
2727 // with him but hasn't done anything yet...
2728 // should use the distance to the centre of gravity of
2729 // the empty bin instead!
2730 double dist = min_distance(bin, x);
2731 if (dist < min_dist)
2732 {
2733 min_dist = dist;
2734 closest_pair = Bin_object_coord_pairs[neigh_bin][j];
2735 }
2736 }
2737 }
2738 }
2739
2740 // Have we filled the bin?
2741 if (min_dist != DBL_MAX)
2742 {
2743 Vector<std::pair<FiniteElement*, Vector<double>>> new_entry;
2744 new_entry.push_back(closest_pair);
2745 Bin_object_coord_pairs.set_value(bin, new_entry);
2746
2747 // Record that we've filled it.
2748 newly_filled_bin.push_back(bin);
2749
2750 // Wipe entry without breaking the linked list (Andrew's trick --
2751 // nice!)
2752 std::list<unsigned>::iterator it_to_be_deleted = it;
2753 it--;
2754 empty_bins.erase(it_to_be_deleted);
2755 }
2756 }
2757
2758 // Get ready for next iteration on remaining empty bins
2759 iter++;
2760 // Now update the vector that records which bins were empty up to now
2761 unsigned n = newly_filled_bin.size();
2762 for (unsigned i = 0; i < n; i++)
2763 {
2764 was_empty_until_current_iteration[newly_filled_bin[i]] = false;
2765 }
2766 }
2767
2768
2769#ifdef PARANOID
2770 // Loop over all bins to check if they're empty
2771 n_bin = nbin();
2772 for (unsigned i = 0; i < n_bin; i++)
2773 {
2774 if (Bin_object_coord_pairs[i].size() == 0)
2775 {
2776 std::ostringstream error_message_stream;
2777 error_message_stream << "Bin " << i << " is still empty\n"
2778 << "after trying to fill it by diffusion\n";
2779 throw OomphLibError(error_message_stream.str(),
2780 OOMPH_CURRENT_FUNCTION,
2781 OOMPH_EXCEPTION_LOCATION);
2782 }
2783 }
2784
2785#endif
2786 }
2787
2788
2789 //=================================================================
2790 /// Get the number of the bin containing the specified coordinate.
2791 /// Bin number is negative if the coordinate is outside
2792 /// the bin structure.
2793 //=================================================================
2794 void NonRefineableBinArray::get_bin(const Vector<double>& zeta,
2795 int& bin_number)
2796 {
2797 // Default for point not found
2798 bin_number = -1;
2799
2800 // Get the lagrangian dimension
2801 const unsigned n_lagrangian = this->ndim_zeta();
2802
2803 // Does the zeta coordinate lie within the current bin structure?
2804
2805 // Loop over the lagrangian dimension
2806 for (unsigned i = 0; i < n_lagrangian; i++)
2807 {
2808 // If the i-th coordinate is less than the minimum
2809 if (zeta[i] < Min_and_max_coordinates[i].first)
2810 {
2811 return;
2812 }
2813 // Otherwise coordinate may be bigger than the maximum
2814 else if (zeta[i] > Min_and_max_coordinates[i].second)
2815 {
2816 return;
2817 }
2818 }
2819
2820 // Use the min and max coords of the bin structure, to find
2821 // the bin structure containing the current zeta cooordinate
2822
2823 // Initialise for subsequent computation
2824 bin_number = 0;
2825
2826 // Offset for rows/matrices in higher dimensions
2827 unsigned multiplier = 1;
2828
2829 // Loop over the dimension
2830 for (unsigned i = 0; i < n_lagrangian; i++)
2831 {
2832 // Find the bin number of the current coordinate
2833 unsigned bin_number_i =
2835 ((zeta[i] - Min_and_max_coordinates[i].first) /
2836 (Min_and_max_coordinates[i].second -
2837 Min_and_max_coordinates[i].first)));
2838 // Buffer the case when we are exactly on the edge
2839 if (bin_number_i == Dimensions_of_bin_array[i])
2840 {
2841 bin_number_i -= 1;
2842 }
2843 // Now add to the bin number using the multiplier
2844 bin_number += multiplier * bin_number_i;
2845 // Increase the current row/matrix multiplier for the next loop
2846 multiplier *= Dimensions_of_bin_array[i];
2847 }
2848
2849
2850#ifdef PARANOID
2851
2852 // Tolerance for "out of bin" test
2853 double tol = 1.0e-10;
2854
2855 unsigned nvertex = (int)pow(2, n_lagrangian);
2856 Vector<Vector<double>> bin_vertex(nvertex);
2857 for (unsigned j = 0; j < nvertex; j++)
2858 {
2859 bin_vertex[j].resize(n_lagrangian);
2860 }
2861 get_bin_vertices(bin_number, bin_vertex);
2862 for (unsigned i = 0; i < n_lagrangian; i++)
2863 {
2864 double min_vertex_coord = DBL_MAX;
2865 double max_vertex_coord = -DBL_MAX;
2866 for (unsigned j = 0; j < nvertex; j++)
2867 {
2868 if (bin_vertex[j][i] < min_vertex_coord)
2869 {
2870 min_vertex_coord = bin_vertex[j][i];
2871 }
2872 if (bin_vertex[j][i] > max_vertex_coord)
2873 {
2874 max_vertex_coord = bin_vertex[j][i];
2875 }
2876 }
2877 if (zeta[i] < min_vertex_coord - tol)
2878 {
2879 std::ostringstream error_message_stream;
2880 error_message_stream
2881 << "Trouble! " << i << " -th coordinate of sample point, " << zeta[i]
2882 << " , isn't actually between limits, " << min_vertex_coord << " and "
2883 << max_vertex_coord << " [it's below by more than " << tol << " !] "
2884 << std::endl;
2885 throw OomphLibError(error_message_stream.str(),
2886 OOMPH_CURRENT_FUNCTION,
2887 OOMPH_EXCEPTION_LOCATION);
2888 }
2889
2890 if (zeta[i] > max_vertex_coord + tol)
2891 {
2892 std::ostringstream error_message_stream;
2893 error_message_stream
2894 << "Trouble! " << i << " -th coordinate of sample point, " << zeta[i]
2895 << " , isn't actually between limits, " << min_vertex_coord << " and "
2896 << max_vertex_coord << " [it's above by more than " << tol << "!] "
2897 << std::endl;
2898 throw OomphLibError(error_message_stream.str(),
2899 OOMPH_CURRENT_FUNCTION,
2900 OOMPH_EXCEPTION_LOCATION);
2901 }
2902 }
2903#endif
2904 }
2905
2906
2907 //========================================================================
2908 /// Get vector of vectors containing the coordinates of the
2909 /// vertices of the i_bin-th bin: bin_vertex[j][i] contains the
2910 /// i-th coordinate of the j-th vertex.
2911 //========================================================================
2913 const unsigned& i_bin, Vector<Vector<double>>& bin_vertex)
2914 {
2915 // Spatial dimension of bin
2916 const unsigned n_lagrangian = this->ndim_zeta();
2917
2918 // How many vertices do we have?
2919 unsigned n_vertices = 1;
2920 for (unsigned i = 0; i < n_lagrangian; i++)
2921 {
2922 n_vertices *= 2;
2923 }
2924 bin_vertex.resize(n_vertices);
2925
2926 // Vectors for min [0] and max [1] coordinates of the bin in each
2927 // coordinate direction
2928 Vector<Vector<double>> zeta_vertex_bin(2);
2929 zeta_vertex_bin[0].resize(n_lagrangian);
2930 zeta_vertex_bin[1].resize(n_lagrangian);
2931
2932 Vector<double> dzeta;
2933 unsigned count = 0;
2934 Vector<unsigned> i_1d;
2935
2936 // Deal with different spatial dimensions separately
2937 switch (n_lagrangian)
2938 {
2939 case 1:
2940
2941 // Assign vertex coordinates of the bin directly
2942 dzeta.resize(1);
2943 dzeta[0] = (Min_and_max_coordinates[0].second -
2944 Min_and_max_coordinates[0].first) /
2945 double(Dimensions_of_bin_array[0]);
2946 bin_vertex[0].resize(1);
2947 bin_vertex[0][0] =
2948 Min_and_max_coordinates[0].first + double(i_bin) * dzeta[0];
2949 bin_vertex[1].resize(1);
2950 bin_vertex[1][0] =
2951 Min_and_max_coordinates[0].first + double(i_bin + 1) * dzeta[0];
2952
2953 break;
2954
2955 case 2:
2956
2957 dzeta.resize(2);
2958 dzeta[0] = (Min_and_max_coordinates[0].second -
2959 Min_and_max_coordinates[0].first) /
2960 double(Dimensions_of_bin_array[0]);
2961 dzeta[1] = (Min_and_max_coordinates[1].second -
2962 Min_and_max_coordinates[1].first) /
2963 double(Dimensions_of_bin_array[1]);
2964
2965 i_1d.resize(2);
2966 i_1d[0] = i_bin % Dimensions_of_bin_array[0];
2967 i_1d[1] = (i_bin - i_1d[0]) / Dimensions_of_bin_array[0];
2968
2969 // Max/min coordinates of the bin
2970 for (unsigned i = 0; i < n_lagrangian; i++)
2971 {
2972 zeta_vertex_bin[0][i] =
2973 Min_and_max_coordinates[i].first + double(i_1d[i]) * dzeta[i];
2974 zeta_vertex_bin[1][i] =
2975 Min_and_max_coordinates[i].first + double(i_1d[i] + 1) * dzeta[i];
2976 }
2977
2978 // Build 4 vertex coordinates
2979 count = 0;
2980 for (unsigned i_min_max = 0; i_min_max < 2; i_min_max++)
2981 {
2982 for (unsigned j_min_max = 0; j_min_max < 2; j_min_max++)
2983 {
2984 bin_vertex[count].resize(2);
2985 bin_vertex[count][0] = zeta_vertex_bin[i_min_max][0];
2986 bin_vertex[count][1] = zeta_vertex_bin[j_min_max][1];
2987 count++;
2988 }
2989 }
2990
2991 break;
2992
2993 case 3:
2994
2995 dzeta.resize(3);
2996 dzeta[0] = (Min_and_max_coordinates[0].second -
2997 Min_and_max_coordinates[0].first) /
2998 double(Dimensions_of_bin_array[0]);
2999 dzeta[1] = (Min_and_max_coordinates[1].second -
3000 Min_and_max_coordinates[1].first) /
3001 double(Dimensions_of_bin_array[1]);
3002 dzeta[2] = (Min_and_max_coordinates[2].second -
3003 Min_and_max_coordinates[2].first) /
3004 double(Dimensions_of_bin_array[2]);
3005
3006 i_1d.resize(3);
3007 i_1d[0] = i_bin % Dimensions_of_bin_array[0];
3008 i_1d[1] = ((i_bin - i_1d[0]) / Dimensions_of_bin_array[0]) %
3010 i_1d[2] = (i_bin - (i_1d[1] * Dimensions_of_bin_array[0] + i_1d[0])) /
3012
3013 // Max/min coordinates of the bin
3014 for (unsigned i = 0; i < n_lagrangian; i++)
3015 {
3016 zeta_vertex_bin[0][i] =
3017 Min_and_max_coordinates[i].first + double(i_1d[i]) * dzeta[i];
3018 zeta_vertex_bin[1][i] =
3019 Min_and_max_coordinates[i].first + double(i_1d[i] + 1) * dzeta[i];
3020 }
3021
3022 // Build 8 vertex coordinates
3023 count = 0;
3024 for (unsigned i_min_max = 0; i_min_max < 2; i_min_max++)
3025 {
3026 for (unsigned j_min_max = 0; j_min_max < 2; j_min_max++)
3027 {
3028 for (unsigned k_min_max = 0; k_min_max < 2; k_min_max++)
3029 {
3030 bin_vertex[count].resize(3);
3031 bin_vertex[count][0] = zeta_vertex_bin[i_min_max][0];
3032 bin_vertex[count][1] = zeta_vertex_bin[j_min_max][1];
3033 bin_vertex[count][2] = zeta_vertex_bin[k_min_max][2];
3034 count++;
3035 }
3036 }
3037 }
3038
3039 break;
3040
3041 default:
3042
3043 std::ostringstream error_message;
3044 error_message << "Can't deal with bins in dimension " << n_lagrangian
3045 << "\n";
3046 throw OomphLibError(error_message.str(),
3047 OOMPH_CURRENT_FUNCTION,
3048 OOMPH_EXCEPTION_LOCATION);
3049 }
3050 }
3051
3052
3053 //========================================================================
3054 /// Compute the minimum distance of any vertex in the specified bin
3055 /// from the specified Lagrangian coordinate zeta.
3056 //========================================================================
3057 double NonRefineableBinArray::min_distance(const unsigned& i_bin,
3058 const Vector<double>& zeta)
3059 {
3060 // Spatial dimension of bin
3061 const unsigned n_lagrangian = this->ndim_zeta();
3062
3063 // Get bin vertices
3064 Vector<Vector<double>> bin_vertex;
3065 get_bin_vertices(i_bin, bin_vertex);
3066 double min_dist = DBL_MAX;
3067 unsigned nvertex = bin_vertex.size();
3068 for (unsigned v = 0; v < nvertex; v++)
3069 {
3070 double dist = 0.0;
3071 for (unsigned i = 0; i < n_lagrangian; i++)
3072 {
3073 dist += pow(bin_vertex[v][i] - zeta[i], 2);
3074 }
3075 dist = sqrt(dist);
3076 if (dist < min_dist) min_dist = dist;
3077 }
3078 return min_dist;
3079 }
3080
3081
3082 //==============================================================================
3083 /// Find the sub geometric object and local coordinate therein that
3084 /// corresponds to the intrinsic coordinate zeta. If sub_geom_object_pt=0
3085 /// on return from this function, none of the constituent sub-objects
3086 /// contain the required coordinate.
3087 //==============================================================================
3088 void NonRefineableBinArray::locate_zeta(const Vector<double>& zeta,
3089 GeomObject*& sub_geom_object_pt,
3090 Vector<double>& s)
3091 {
3092 // Reset counter for number of sample points visited only if we're
3093 // starting from the beginning
3094 if (current_min_spiral_level() == 0)
3095 {
3097 0;
3098 }
3099
3100 // Initialise return to null -- if it's still null when we're
3101 // leaving we've failed!
3102 sub_geom_object_pt = 0;
3103
3104 // Get the lagrangian dimension
3105 const unsigned n_lagrangian = this->ndim_zeta();
3106
3107 // Does the zeta coordinate lie within the current bin structure?
3108 // Skip this test if we want to always fail because that's usually
3109 // done to trace out the spiral path
3111 {
3112 // Loop over the lagrangian dimension
3113 for (unsigned i = 0; i < n_lagrangian; i++)
3114 {
3115 // If the i-th coordinate is less than the minimum
3116 if (zeta[i] < Min_and_max_coordinates[i].first)
3117 {
3118 return;
3119 }
3120 // Otherwise coordinate may be bigger than the maximum
3121 else if (zeta[i] > Min_and_max_coordinates[i].second)
3122 {
3123 return;
3124 }
3125 }
3126 }
3127
3128 // Use the min and max coords of the bin structure, to find
3129 // the bin structure containing the current zeta cooordinate
3130 unsigned bin_number = 0;
3131
3132 // Offset for rows/matrices in higher dimensions
3133 unsigned multiplier = 1;
3134
3135 // Loop over the dimension
3136 for (unsigned i = 0; i < n_lagrangian; i++)
3137 {
3138 // Find the bin number of the current coordinate
3139 unsigned bin_number_i =
3141 ((zeta[i] - Min_and_max_coordinates[i].first) /
3142 (Min_and_max_coordinates[i].second -
3143 Min_and_max_coordinates[i].first)));
3144
3145 // Buffer the case when we are exactly on the edge
3146 if (bin_number_i == Dimensions_of_bin_array[i])
3147 {
3148 bin_number_i -= 1;
3149 }
3150
3151 // Now add to the bin number using the multiplier
3152 bin_number += multiplier * bin_number_i;
3153
3154 // Increase the current row/matrix multiplier for the next loop
3155 multiplier *= Dimensions_of_bin_array[i];
3156 }
3157
3158 // Loop over spirals that are to be visited in one go
3159 Vector<unsigned> neighbour_bin;
3160 unsigned min_level = current_min_spiral_level();
3161 unsigned max_level = current_max_spiral_level();
3162 for (unsigned i_level = min_level; i_level <= max_level; i_level++)
3163 {
3164 // Call helper function to find the neighbouring bins at this
3165 // level
3166 get_neighbouring_bins_helper(bin_number, i_level, neighbour_bin);
3167
3168 // Total number of bins to be visited
3169 unsigned n_nbr_bin = neighbour_bin.size();
3170
3171 // // keep around and add "false" as last argument to previous call
3172 // // to get_neighbouring_bins_helper(...) for annotation below to make
3173 // sense
3174 // {
3175 // Vector<unsigned> old_neighbour_bin;
3176 // get_neighbouring_bins_helper(bin_number,
3177 // i_level,
3178 // old_neighbour_bin,
3179 // true); // true=use old version!
3180 // unsigned nbin_new = neighbour_bin.size();
3181 // unsigned nbin_old = old_neighbour_bin.size();
3182 // unsigned n=std::min(nbin_new,nbin_old);
3183 // std::sort(old_neighbour_bin.begin(),
3184 // old_neighbour_bin.end());
3185 // std::sort(neighbour_bin.begin(),
3186 // neighbour_bin.end());
3187 // oomph_info << "# of neighbouring bins: "
3188 // << nbin_new << " " << nbin_old;
3189 // if (nbin_new!=nbin_old)
3190 // {
3191 // oomph_info << " DIFFERENT!";
3192 // }
3193 // oomph_info << std::endl;
3194 // for (unsigned i=0;i<n;i++)
3195 // {
3196 // oomph_info << neighbour_bin[i] << " "
3197 // << old_neighbour_bin[i] << " ";
3198 // if (neighbour_bin[i]!=
3199 // old_neighbour_bin[i])
3200 // {
3201 // oomph_info << "DIFFFFERENT";
3202 // }
3203 // oomph_info << std::endl;
3204 // }
3205 // if (nbin_new>nbin_old)
3206 // {
3207 // for (unsigned i=n;i<nbin_new;i++)
3208 // {
3209 // oomph_info << "NNEW:" << neighbour_bin[i]
3210 // << std::endl;
3211 // }
3212 // }
3213 // if (nbin_old>nbin_new)
3214 // {
3215 // for (unsigned i=n;i<nbin_old;i++)
3216 // {
3217 // oomph_info << "OOLD:" << old_neighbour_bin[i]
3218 // << std::endl;
3219 // }
3220 // }
3221 // }
3222
3223
3224 // Set bool for finding zeta
3225 bool found_zeta = false;
3226 for (unsigned i_nbr = 0; i_nbr < n_nbr_bin; i_nbr++)
3227 {
3228 // Only search if bin is within the max. radius but min_distance
3229 // is expensive...
3230 bool do_it = true;
3231 if (Max_search_radius < DBL_MAX)
3232 {
3233 if (min_distance(neighbour_bin[i_nbr], zeta) > Max_search_radius)
3234 {
3235 do_it = false;
3236 }
3237 }
3238 if (do_it)
3239 {
3240 // Get the number of element-sample point pairs in this bin
3241 unsigned n_sample =
3242 Bin_object_coord_pairs[neighbour_bin[i_nbr]].size();
3243
3244 // Don't do anything if this bin has no sample points
3245 if (n_sample > 0)
3246 {
3247 for (unsigned i_sam = 0; i_sam < n_sample; i_sam++)
3248 {
3249 // Get the element
3250 FiniteElement* el_pt =
3251 Bin_object_coord_pairs[neighbour_bin[i_nbr]][i_sam].first;
3252
3253 // Get the local coordinate
3254 s = Bin_object_coord_pairs[neighbour_bin[i_nbr]][i_sam].second;
3255
3256 // History of sample points visited
3258 {
3259 unsigned cached_dim_zeta = this->ndim_zeta();
3260 Vector<double> zeta_sample(cached_dim_zeta);
3262 {
3263 el_pt->interpolated_x(s, zeta_sample);
3264 }
3265 else
3266 {
3267 el_pt->interpolated_zeta(s, zeta_sample);
3268 }
3269 double dist = 0.0;
3270 for (unsigned ii = 0; ii < cached_dim_zeta; ii++)
3271 {
3272 BinArray::Visited_sample_points_file << zeta_sample[ii]
3273 << " ";
3274 dist +=
3275 (zeta[ii] - zeta_sample[ii]) * (zeta[ii] - zeta_sample[ii]);
3276 }
3279 << " " << sqrt(dist) << std::endl;
3280 }
3281
3282 // Use this coordinate as the initial guess
3283 bool use_coordinate_as_initial_guess = true;
3284
3285 // Attempt to find zeta within a sub-object
3286 el_pt->locate_zeta(
3287 zeta, sub_geom_object_pt, s, use_coordinate_as_initial_guess);
3288
3289
3291
3292 // Always fail? (Used for debugging, e.g. to trace out
3293 // spiral path)
3295 {
3296 sub_geom_object_pt = 0;
3297 }
3298
3299
3300#ifdef OOMPH_HAS_MPI
3301 // Ignore halos?
3303 {
3304 // Dynamic cast the result to a FiniteElement
3305 FiniteElement* test_el_pt =
3306 dynamic_cast<FiniteElement*>(sub_geom_object_pt);
3307 if (test_el_pt != 0)
3308 {
3309 // We only want to exit if this is a non-halo element
3310 if (test_el_pt->is_halo())
3311 {
3312 sub_geom_object_pt = 0;
3313 }
3314 }
3315 }
3316#endif
3317
3318 // If the FiniteElement is non-halo and has been located, exit
3319 if (sub_geom_object_pt != 0)
3320 {
3321 found_zeta = true;
3322 break;
3323 }
3324 } // end loop over sample points
3325 }
3326
3327
3328 if (found_zeta)
3329 {
3330 return; // break;
3331 }
3332
3333 } // end of don't search if outside search radius
3334 } // end loop over bins at this level
3335 } // End of loop over levels
3336 }
3337
3338
3339 /// Default number of bins (in each coordinate direction)
3341
3342 /// Counter for overall number of bins allocated -- used to
3343 /// issue warning if this exceeds a threshhold. (Default assignment
3344 /// of 100^DIM bins per MeshAsGeomObject can be a killer if there
3345 /// are huge numbers of sub-meshes (e.g. in unstructured FSI).
3347
3348 /// Total number of bins above which warning is issued.
3349 /// (Default assignment of 100^DIM bins per MeshAsGeomObject can
3350 /// be a killer if there are huge numbers of sub-meshes (e.g. in
3351 /// unstructured FSI).
3352 unsigned long
3354 50000000;
3355
3356 /// Boolean to supppress warnings about large number of bins
3357 bool
3359 false;
3360
3361 /// Boolean flag to make sure that warning about large number
3362 /// of bin cells only gets triggered once.
3364 false;
3365
3366 /// Fraction of elements/bin that triggers warning. Too many
3367 /// elements per bin can lead to very slow computations
3369
3370 /// Boolean to supppress warnings about small number of bins
3372 false;
3373
3374 /// Boolean flag to make sure that warning about small number
3375 /// of bin cells only gets triggered once.
3377 false;
3378
3379
3380 /// /////////////////////////////////////////////////////////////////////////////
3381 /// /////////////////////////////////////////////////////////////////////////////
3382 /// /////////////////////////////////////////////////////////////////////////////
3383
3384#ifdef OOMPH_HAS_CGAL
3385
3386
3387 //====================================================================
3388 /// Constructor
3389 //====================================================================
3391 SamplePointContainerParameters* sample_point_container_parameters_pt)
3393 sample_point_container_parameters_pt->mesh_pt(),
3394 sample_point_container_parameters_pt->min_and_max_coordinates(),
3395 sample_point_container_parameters_pt
3396 ->use_eulerian_coordinates_during_setup(),
3397 sample_point_container_parameters_pt
3398 ->ignore_halo_elements_during_locate_zeta_search(),
3399 sample_point_container_parameters_pt
3400 ->nsample_points_generated_per_element())
3401 {
3402 // Get the spatial dimension (int because of mpi below)
3403 int dim = 0;
3404 if (Mesh_pt->nelement() != 0)
3405 {
3406 dim = Mesh_pt->finite_element_pt(0)->dim();
3407 }
3408
3409 // Need to do an Allreduce to ensure that the dimension is consistent
3410 // even when no elements are assigned to a certain processor
3411#ifdef OOMPH_HAS_MPI
3412 // Only a problem if the mesh has been distributed
3413 if (Mesh_pt->is_mesh_distributed())
3414 {
3415 // Need a non-null communicator
3416 if (Mesh_pt->communicator_pt() != 0)
3417 {
3418 int n_proc = Mesh_pt->communicator_pt()->nproc();
3419 if (n_proc > 1)
3420 {
3421 int dim_reduce;
3422 MPI_Allreduce(&dim,
3423 &dim_reduce,
3424 1,
3425 MPI_INT,
3426 MPI_MAX,
3427 Mesh_pt->communicator_pt()->mpi_comm());
3428 dim = dim_reduce;
3429 }
3430 }
3431 }
3432#endif
3433
3434 Ndim_zeta = dim;
3435
3436 // Have we specified max/min coordinates?
3437 // If not, compute them on the fly from mesh
3438 if (Min_and_max_coordinates.size() == 0)
3439 {
3441 }
3442
3443
3444 // Time it
3445 double t_start = 0.0;
3447 {
3448 t_start = TimingHelpers::timer();
3449 }
3450
3451 // Fill the bastard!
3452 double CGAL_setup_time = get_sample_points();
3453
3455 {
3456 double t_end = TimingHelpers::timer();
3458 oomph_info << "Time for setup of " << dim
3459 << "-dimensional sample point container containing " << npts
3460 << " sample points: " << t_end - t_start
3461 << " sec (cgal); third party: " << CGAL_setup_time
3462 << " sec ( = " << CGAL_setup_time / (t_end - t_start) * 100.0
3463 << " %)" << std::endl;
3464 }
3465
3466 // Initialise
3471 2; // hierher tune this and create public static default
3473 10; // UINT MAX is temporary bypass! tune this and create public static
3474 // default
3475 }
3476
3477 //==============================================================================
3478 /// Get the sample points; returns time taken for setup of CGAL tree
3479 //==============================================================================
3481 {
3482 // Number of elements
3483 unsigned nel = Mesh_pt->nelement();
3484
3485 // Estimate number of sample points
3486 unsigned n_sample_estimate = 0;
3487 if (nel > 0)
3488 {
3489 FiniteElement* el_pt = Mesh_pt->finite_element_pt(0);
3490 if (el_pt != 0)
3491 {
3492 // Total number of sample point we will create
3493 n_sample_estimate =
3494 nel * el_pt->nplot_points(Nsample_points_generated_per_element);
3495 }
3496 }
3497 CGAL_sample_point_zeta_d.reserve(n_sample_estimate);
3498 Sample_point_pt.reserve(n_sample_estimate);
3499
3500 // Fill 'em in:
3501 for (unsigned e = 0; e < nel; e++)
3502 {
3503 FiniteElement* el_pt = Mesh_pt->finite_element_pt(e);
3504
3505 // Total number of sample point we will create
3506 unsigned nplot =
3507 el_pt->nplot_points(Nsample_points_generated_per_element);
3508
3509 /// For all the sample points we have to create ...
3510 for (unsigned j = 0; j < nplot; j++)
3511 {
3512 // ... create it: Pass element index in mesh (vector
3513 // of elements and index of sample point within element
3514 SamplePoint* new_sample_point_pt = new SamplePoint(e, j);
3515
3516 // Coordinates of this point
3517 Vector<double> zeta(ndim_zeta());
3518 Vector<double> s(ndim_zeta());
3519 bool use_equally_spaced_interior_sample_points =
3521 el_pt->get_s_plot(j,
3523 s,
3524 use_equally_spaced_interior_sample_points);
3526 {
3527 el_pt->interpolated_x(s, zeta);
3528 }
3529 else
3530 {
3531 el_pt->interpolated_zeta(s, zeta);
3532 }
3533
3534#ifdef PARANOID
3535
3536 // Check if point is inside
3537 bool is_inside = true;
3538 std::ostringstream error_message;
3539 unsigned dim = ndim_zeta();
3540 for (unsigned i = 0; i < dim; i++)
3541 {
3542 if ((zeta[i] < Min_and_max_coordinates[i].first) ||
3543 (zeta[i] > Min_and_max_coordinates[i].second))
3544 {
3545 is_inside = false;
3546 error_message << "Sample point at zeta[" << i << "] = " << zeta[i]
3547 << " is outside limits of sample point container: "
3548 << Min_and_max_coordinates[i].first << " and "
3549 << Min_and_max_coordinates[i].second << std::endl;
3550 }
3551 }
3552
3553 if (!is_inside)
3554 {
3555 error_message << "Please correct the limits passed to the "
3556 << "constructor." << std::endl;
3557 throw OomphLibError(error_message.str(),
3558 OOMPH_CURRENT_FUNCTION,
3559 OOMPH_EXCEPTION_LOCATION);
3560 }
3561
3562#endif
3563
3564 CGAL_sample_point_zeta_d.push_back(
3565 Kernel_d::Point_d(zeta.size(), zeta.begin(), zeta.end()));
3566 Sample_point_pt.push_back(new_sample_point_pt);
3567 }
3568 }
3569
3570 // Make tree structure
3571 double CGAL_setup_time = 0.0;
3573 {
3574 CGAL_setup_time = TimingHelpers::timer();
3575 }
3576
3577 CGAL_tree_d_pt = new K_neighbor_search_d::Tree(
3578 boost::make_zip_iterator(boost::make_tuple(
3579 CGAL_sample_point_zeta_d.begin(), Sample_point_pt.begin())),
3580 boost::make_zip_iterator(boost::make_tuple(CGAL_sample_point_zeta_d.end(),
3581 Sample_point_pt.end())));
3583 {
3584 CGAL_setup_time = TimingHelpers::timer() - CGAL_setup_time;
3585 }
3586 return CGAL_setup_time;
3587 }
3588
3589
3590 //==============================================================================
3591 /// Compute total number of sample points in sample point container
3592 //==============================================================================
3595 {
3596 return Sample_point_pt.size();
3597 }
3598
3599
3600 //==============================================================================
3601 /// Find the sub geometric object and local coordinate therein that
3602 /// corresponds to the intrinsic coordinate zeta. If sub_geom_object_pt=0
3603 /// on return from this function, none of the constituent sub-objects
3604 /// contain the required coordinate.
3605 //==============================================================================
3606 void CGALSamplePointContainer::locate_zeta(const Vector<double>& zeta,
3607 GeomObject*& sub_geom_object_pt,
3608 Vector<double>& s)
3609 {
3610 // Top level book keeping and sanity checking
3612 {
3613 // Reset counter for number of sample points visited.
3614 // If we can't find the point we should at least make sure that
3615 // we've visited all the sample points before giving up.
3617 0;
3618 }
3619
3620 // Initialise return to null -- if it's still null when we're
3621 // leaving we've failed!
3622 sub_geom_object_pt = 0;
3623
3624 // Get the lagrangian dimension
3625 const unsigned n_lagrangian = this->ndim_zeta();
3626
3627 // Does the zeta coordinate lie within the current bin structure?
3628 // Skip this test if we want to always fail because that's usually
3629 // done to trace out the spiral path
3630 if (!SamplePointContainer::Always_fail_elemental_locate_zeta)
3631 {
3632 // Loop over the lagrangian dimension
3633 for (unsigned i = 0; i < n_lagrangian; i++)
3634 {
3635 // If the i-th coordinate is less than the minimum
3636 if (zeta[i] < Min_and_max_coordinates[i].first)
3637 {
3638 return;
3639 }
3640 // Otherwise coordinate may be bigger than the maximum
3641 else if (zeta[i] > Min_and_max_coordinates[i].second)
3642 {
3643 return;
3644 }
3645 }
3646 }
3647
3648 // Number of sample points
3649 unsigned n_sample = Sample_point_pt.size();
3650
3651 // Create CGAL query -- this is the point we want!
3652 Point_d query(zeta.size(), zeta.begin(), zeta.end());
3653
3654 // Max. number of nearest neighbours
3655 const unsigned n_nearest_neighbours_max = std::min(
3657
3658 // Start with just one...
3659 unsigned n_nearest_neighbours = 1;
3660 unsigned n_neighbours_visited_last_time = 0;
3661 bool can_increase_n_nearest_neighbours = true;
3662 bool keep_going = true;
3663 while (keep_going)
3664 {
3665 // Find specified number of nearest neighbours only
3666 const unsigned n_nearest_neighbours_actual =
3667 std::min(n_nearest_neighbours, n_nearest_neighbours_max);
3668 K_neighbor_search_d search(
3669 *CGAL_tree_d_pt, query, n_nearest_neighbours_actual);
3670
3671 // Search
3672 unsigned count = 0;
3673 for (K_neighbor_search_d::iterator it = search.begin();
3674 it != search.end();
3675 it++)
3676 {
3677 count++;
3678
3679 if ((count > n_neighbours_visited_last_time) &&
3680 (count >
3682 {
3683 // Recover the sample point
3684 SamplePoint* sample_point_pt = boost::get<1>(it->first);
3685
3686 // Get the element
3687 FiniteElement* el_pt = Mesh_pt->finite_element_pt(
3688 sample_point_pt->element_index_in_mesh());
3689
3690
3691#ifdef OOMPH_HAS_MPI
3692 // We only look at the sample point if it isn't halo
3693 // if we are set up to ignore the halo elements
3695 (el_pt->is_halo()))
3696 {
3697 // Halo
3698 }
3699 else
3700 {
3701#endif
3702
3703 // Provide initial guess for Newton search using local coordinate
3704 // of sample point
3705 bool use_equally_spaced_interior_sample_points =
3707 unsigned i = sample_point_pt->sample_point_index_in_element();
3708 el_pt->get_s_plot(i,
3710 s,
3711 use_equally_spaced_interior_sample_points);
3712
3713 bool do_it = true;
3714 if (Max_search_radius < DBL_MAX)
3715 {
3716 unsigned cached_dim_zeta = ndim_zeta();
3717 Vector<double> zeta_sample(cached_dim_zeta);
3719 {
3720 el_pt->interpolated_x(s, zeta_sample);
3721 }
3722 else
3723 {
3724 el_pt->interpolated_zeta(s, zeta_sample);
3725 }
3726 double dist_sq = 0.0;
3727 for (unsigned ii = 0; ii < cached_dim_zeta; ii++)
3728 {
3729 dist_sq +=
3730 (zeta[ii] - zeta_sample[ii]) * (zeta[ii] - zeta_sample[ii]);
3731 }
3732 if (dist_sq > Max_search_radius * Max_search_radius)
3733 {
3734 do_it = false;
3735 }
3736 }
3737 if (do_it)
3738 {
3739 // History of sample points visited
3741 {
3742 unsigned cached_dim_zeta = ndim_zeta();
3743 Vector<double> zeta_sample(cached_dim_zeta);
3745 {
3746 el_pt->interpolated_x(s, zeta_sample);
3747 }
3748 else
3749 {
3750 el_pt->interpolated_zeta(s, zeta_sample);
3751 }
3752 double dist = 0.0;
3753 for (unsigned ii = 0; ii < cached_dim_zeta; ii++)
3754 {
3756 << zeta_sample[ii] << " ";
3757 dist +=
3758 (zeta[ii] - zeta_sample[ii]) * (zeta[ii] - zeta_sample[ii]);
3759 }
3762 << " " << sqrt(dist) << std::endl;
3763 }
3764
3765 // Bump counter
3767
3768 bool use_coordinate_as_initial_guess = true;
3769 el_pt->locate_zeta(
3770 zeta, sub_geom_object_pt, s, use_coordinate_as_initial_guess);
3771
3772 // Always fail? (Used for debugging, e.g. to trace out
3773 // spiral path)
3775 {
3776 sub_geom_object_pt = 0;
3777 }
3778
3779 if (sub_geom_object_pt != 0)
3780 {
3781 return;
3782 }
3783 }
3784
3785#ifdef OOMPH_HAS_MPI
3786 }
3787#endif
3788 }
3789 }
3790
3791 n_neighbours_visited_last_time = count;
3792
3793 // Can we increase the number of neighbours further?
3794 if (can_increase_n_nearest_neighbours)
3795 {
3796 unsigned factor_for_increase_in_nearest_neighbours = 10;
3797 n_nearest_neighbours *= factor_for_increase_in_nearest_neighbours;
3798
3799 // There's no point going any further (next time)
3800 if (n_nearest_neighbours > n_nearest_neighbours_max)
3801 {
3802 can_increase_n_nearest_neighbours = false;
3803 }
3804 }
3805 // Bailing out; not found but we can't increase number of search pts
3806 // further
3807 else
3808 {
3809 keep_going = false;
3810 }
3811 } // while loop to increase number of nearest neighbours
3812 }
3813
3814
3815 //==============================================================================
3816 /// Find the sub geometric object and local coordinate therein that
3817 /// corresponds to the intrinsic coordinate zeta, using up to the specified
3818 /// number of sample points as initial guess for the Newton-based search.
3819 /// If this fails, return the nearest sample point.
3820 //==============================================================================
3822 const Vector<double>& zeta,
3823 const unsigned& max_sample_points_for_newton_based_search,
3824 GeomObject*& sub_geom_object_pt,
3825 Vector<double>& s)
3826 {
3827 // Reset counter for number of sample points visited.
3829
3830 // Initialise return to null -- if it's still null when we're
3831 // leaving we've failed!
3832 sub_geom_object_pt = 0;
3833
3834 // Number of sample points
3835 unsigned n_sample = Sample_point_pt.size();
3836
3837 // Create CGAL query -- this is the point we want!
3838 Point_d query(zeta.size(), zeta.begin(), zeta.end());
3839
3840 // Max. number of nearest neighbours
3841 const unsigned n_nearest_neighbours_max =
3842 std::min(n_sample, max_sample_points_for_newton_based_search);
3843
3844 // Find 'em
3845 K_neighbor_search_d search(
3846 *CGAL_tree_d_pt, query, n_nearest_neighbours_max);
3847
3848 // Do Newton method starting from each of the nearest sample points
3849 for (K_neighbor_search_d::iterator it = search.begin(); it != search.end();
3850 it++)
3851 {
3852 // Recover the sample point
3853 SamplePoint* sample_point_pt = boost::get<1>(it->first);
3854
3855 // Get the element
3856 FiniteElement* el_pt =
3857 Mesh_pt->finite_element_pt(sample_point_pt->element_index_in_mesh());
3858
3859
3860#ifdef OOMPH_HAS_MPI
3861
3862 // We only look at the sample point if it isn't halo
3863 // if we are set up to ignore the halo elements
3865 (el_pt->is_halo()))
3866 {
3867 // Halo
3868 }
3869 else
3870 { // not halo
3871
3872#endif
3873
3874 // Provide initial guess for Newton search using local coordinate
3875 // of sample point
3876 bool use_equally_spaced_interior_sample_points =
3878 unsigned i = sample_point_pt->sample_point_index_in_element();
3879 el_pt->get_s_plot(i,
3881 s,
3882 use_equally_spaced_interior_sample_points);
3883
3884 bool do_it = true;
3885 if (Max_search_radius < DBL_MAX)
3886 {
3887 unsigned cached_dim_zeta = ndim_zeta();
3888 Vector<double> zeta_sample(cached_dim_zeta);
3890 {
3891 el_pt->interpolated_x(s, zeta_sample);
3892 }
3893 else
3894 {
3895 el_pt->interpolated_zeta(s, zeta_sample);
3896 }
3897 double dist_sq = 0.0;
3898 for (unsigned ii = 0; ii < cached_dim_zeta; ii++)
3899 {
3900 dist_sq +=
3901 (zeta[ii] - zeta_sample[ii]) * (zeta[ii] - zeta_sample[ii]);
3902 }
3903 if (dist_sq > Max_search_radius * Max_search_radius)
3904 {
3905 do_it = false;
3906 }
3907 }
3908 if (do_it)
3909 {
3910 // History of sample points visited
3912 {
3913 unsigned cached_dim_zeta = ndim_zeta();
3914 Vector<double> zeta_sample(cached_dim_zeta);
3916 {
3917 el_pt->interpolated_x(s, zeta_sample);
3918 }
3919 else
3920 {
3921 el_pt->interpolated_zeta(s, zeta_sample);
3922 }
3923 double dist = 0.0;
3924 for (unsigned ii = 0; ii < cached_dim_zeta; ii++)
3925 {
3927 << zeta_sample[ii] << " ";
3928 dist +=
3929 (zeta[ii] - zeta_sample[ii]) * (zeta[ii] - zeta_sample[ii]);
3930 }
3933 << " " << sqrt(dist) << std::endl;
3934 }
3935
3936 // Bump counter
3938
3939 bool use_coordinate_as_initial_guess = true;
3940 el_pt->locate_zeta(
3941 zeta, sub_geom_object_pt, s, use_coordinate_as_initial_guess);
3942
3943 // Always fail? (Used for debugging, e.g. to trace out
3944 // spiral path)
3946 {
3947 sub_geom_object_pt = 0;
3948 }
3949
3950 if (sub_geom_object_pt != 0)
3951 {
3952 return;
3953 }
3954 }
3955
3956#ifdef OOMPH_HAS_MPI
3957 }
3958#endif
3959 }
3960
3961 // We've searched over all the sample points but the Newton method
3962 // hasn't converged from any, so just use the nearest one
3963 K_neighbor_search_d::iterator it = search.begin();
3964
3965 // Recover the sample point
3966 SamplePoint* sample_point_pt = boost::get<1>(it->first);
3967
3968 // Get the element
3969 FiniteElement* el_pt =
3970 Mesh_pt->finite_element_pt(sample_point_pt->element_index_in_mesh());
3971
3972 // Get local coordinate of sample point in element
3973 bool use_equally_spaced_interior_sample_points =
3975 unsigned i = sample_point_pt->sample_point_index_in_element();
3976 el_pt->get_s_plot(i,
3978 s,
3979 use_equally_spaced_interior_sample_points);
3980
3981 sub_geom_object_pt = el_pt;
3982 }
3983
3984
3985#endif // cgal
3986
3987} // namespace oomph
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
unsigned dimension_of_bin_array(const unsigned &i) const
Number of bins in coordinate direction i.
unsigned max_bin_dimension() const
Max. bin dimension (number of bins in coordinate directions)
Vector< unsigned > dimensions_of_bin_array() const
Number of bins in coordinate directions. Const vector-based version.
void get_neighbouring_bins_helper(const unsigned &bin_index, const unsigned &radius, Vector< unsigned > &neighbouring_bin_index, const bool &use_old_version=true)
Helper function for computing the bin indices of neighbouring bins at a given "radius" of the specifi...
unsigned ndim_zeta() const
Dimension of the zeta ( = dim of local coordinate of elements)
Vector< unsigned > Dimensions_of_bin_array
Number of bins in each coordinate direction.
unsigned coords_to_bin_index(const Vector< double > &zeta)
Get (linearly enumerated) bin index of bin that contains specified zeta.
void profile_get_neighbouring_bins_helper()
Profiling function to compare performance of two different versions of the get_neighbouring_bins_help...
void coords_to_vectorial_bin_index(const Vector< double > &zeta, Vector< unsigned > &bin_index)
Get "coordinates" of bin that contains specified zeta.
void locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s)
Find sub-GeomObject (finite element) and the local coordinate s within it that contains point with gl...
unsigned ndim_zeta() const
Dimension of the zeta ( = dim of local coordinate of elements)
CGAL::Orthogonal_k_neighbor_search< Traits_d > K_neighbor_search_d
CGALSamplePointContainer(SamplePointContainerParameters *sample_point_container_parameters_pt)
Constructor.
void limited_locate_zeta(const Vector< double > &zeta, const unsigned &max_sample_points_for_newton_based_search, GeomObject *&sub_geom_object_pt, Vector< double > &s)
Find the sub geometric object and local coordinate therein that corresponds to the intrinsic coordina...
unsigned Multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta
Every time we've completed a "spiral", visiting a finite number of sample points in a deterministic o...
unsigned & last_sample_point_to_actually_lookup_during_locate_zeta()
When searching through sample points only actually do the locate_zeta calls when when the counter is ...
unsigned Last_sample_point_to_actually_lookup_during_locate_zeta
When searching through sample points only actually do the locate_zeta calls when when the counter is ...
Vector< SamplePoint * > Sample_point_pt
Vector storing pointers to sample point objects (which represent sample point in terms of number of e...
double get_sample_points()
Get the sample points; return time for setup of CGAL tree.
unsigned Initial_last_sample_point_to_actually_lookup_during_locate_zeta
When searching through sample points only actually do the locate_zeta calls when when the counter exc...
unsigned First_sample_point_to_actually_lookup_during_locate_zeta
When searching through sample points only actually do the locate_zeta calls when when the counter exc...
Vector< Point_d > CGAL_sample_point_zeta_d
Vector containing sample point coordinates.
K_neighbor_search_d::Tree * CGAL_tree_d_pt
Pointer to tree-based representation of sample points.
unsigned & first_sample_point_to_actually_lookup_during_locate_zeta()
When searching through sample points only actually do the locate_zeta calls when when the counter exc...
unsigned Ndim_zeta
Dimension of the zeta ( = dim of local coordinate of elements)
unsigned total_number_of_sample_points_computed_recursively() const
Compute total number of sample points in sample point container.
static bool Already_warned_about_small_number_of_bin_cells
Boolean flag to make sure that warning about small number of bin cells only gets triggered once.
void get_bin(const Vector< double > &zeta, int &bin_number)
Get the number of the bin containing the specified coordinate. Bin number is negative if the coordina...
static bool Suppress_warning_about_large_total_number_of_bins
Boolean to supppress warnings about large number of bins.
void fill_bin_array()
Fill the bin array with sample points from FiniteElements stored in mesh.
unsigned Nspiral_chunk
Number of spirals to be searched in one go.
unsigned Max_spiral_level
Max. spiralling level (for efficiency; effect similar to max_search_radius)
static unsigned long Total_nbin_cells_counter
Counter for overall number of bins allocated – used to issue warning if this exceeds a threshhold....
void locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s)
Find sub-GeomObject (finite element) and the local coordinate s within it that contains point with gl...
static unsigned Threshold_for_elements_per_bin_warning
Fraction of elements/bin that triggers warning. Too many elements per bin can lead to very slow compu...
unsigned Current_min_spiral_level
Current min. spiralling level.
unsigned Current_max_spiral_level
Current max. spiralling level.
static bool Suppress_warning_about_small_number_of_bins
Boolean to supppress warnings about small number of bins.
unsigned & current_min_spiral_level()
Access function to current min. spiral level.
void fill_bin_by_diffusion(const unsigned &bin_diffusion_radius=1)
Fill bin by diffusion, populating each empty bin with the same content as the first non-empty bin fou...
void get_bin_vertices(const unsigned &i_bin, Vector< Vector< double > > &bin_vertex)
Get vector of vectors containing the coordinates of the vertices of the i_bin-th bin: bin_vertex[j][i...
unsigned total_number_of_sample_points_computed_recursively() const
Compute total number of sample points recursively.
void output_bin_vertices(std::ofstream &outfile)
Output bin vertices (allowing display of bins as zones).
static bool Already_warned_about_large_number_of_bin_cells
Boolean flag to make sure that warning about large number of bin cells only gets triggered once.
SparseVector< Vector< std::pair< FiniteElement *, Vector< double > > > > Bin_object_coord_pairs
Storage for paired objects and coords in each bin.
void output_bins(std::ofstream &outfile)
Output bins.
static unsigned Default_n_bin_1d
Default number of bins (in each coordinate direction). (Note: don't move this into a common base clas...
void get_fill_stats(unsigned &n_bin, unsigned &max_n_entry, unsigned &min_n_entry, unsigned &tot_n_entry, unsigned &n_empty) const
Provide some stats on the fill level of the associated bin.
void flush_bins_of_objects()
Flush the storage for the binning method (and decrement counter for total number of bins in active us...
static unsigned long Threshold_for_total_bin_cell_number_warning
Total number of bins above which warning is issued. (Default assignment of 100^DIM bins per MeshAsGeo...
NonRefineableBinArray(SamplePointContainerParameters *bin_array_parameters_pt)
Constructor.
double min_distance(const unsigned &i_bin, const Vector< double > &zeta)
Compute the minimum distance of any vertex in the specified bin from the specified Lagrangian coordin...
unsigned nbin() const
Total number of bins (empty or not)
unsigned & current_max_spiral_level()
Access function to current max. spiral level.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
unsigned Multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta
Every time we've completed a "spiral", visiting a finite number of sample points in a deterministic o...
Vector< RefineableBin * > Bin_pt
Vector of pointers to constituent RefineableBins.
unsigned Max_depth
Max depth of the hierarchical bin_array.
unsigned total_number_of_sample_points_computed_recursively() const
Compute total number of sample points recursively.
unsigned First_sample_point_to_actually_lookup_during_locate_zeta
When searching through sample points recursively from the top level RefineableBinArray (in determinis...
static unsigned Default_n_bin_1d
Default number of bins (in each coordinate direction) (Note: don't move this into a common base class...
void locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s)
Find sub-GeomObject (finite element) and the local coordinate s within it that contains point with gl...
unsigned nbin() const
Number of bins (not taking recursion into account)
RefineableBinArray * Root_bin_array_pt
Pointer to root bin array.
unsigned Depth
Variable which stores the Depth value of the bin_array. Useful for debugging and for preventing "infi...
RefineableBinArray * root_bin_array_pt() const
Root bin array.
unsigned Max_number_of_sample_point_per_bin
Maximum number of sample points in bin (before it's subdivided recursively)
void fill_bin_array()
Fill the bin array with sample points from FiniteElements stored in mesh.
void output_bin_vertices(std::ofstream &outfile)
Output bin vertices (allowing display of bins as zones).
void output_neighbouring_bins(const unsigned &bin_index, const unsigned &radius, std::ofstream &outfile)
Output neighbouring bins at given "radius" of the specified bin.
unsigned Last_sample_point_to_actually_lookup_during_locate_zeta
When searching through sample points recursively from the top level RefineableBinArray (in determinis...
void get_bin_boundaries(const unsigned &bin_index, Vector< std::pair< double, double > > &min_and_max_coordinates)
Boundaries of specified bin in each coordinate direction. *.first = min; *.second = max.
unsigned Initial_last_sample_point_to_actually_lookup_during_locate_zeta
When searching through sample points recursively from the top level RefineableBinArray (in determinis...
RefineableBinArray(SamplePointContainerParameters *bin_array_parameters_pt)
Constructor.
bool Bin_array_is_recursive
Variable which stores if the RefineableBinArray is recursive or not.
RefineableBin class. Contains sample points and is embedded in a RefineableBinArray....
void make_sub_bin_array(const Vector< std::pair< double, double > > &min_and_max_coordinates)
Method for building a new subbin_array (called when the Bin size is greater than the Max_number_of_sa...
void get_bin_boundaries(Vector< std::pair< double, double > > &min_and_max_coordinates)
Boundaries of bin in each coordinate direction. *.first = min; *.second = max.
void output(std::ofstream &outfile, const bool &don_t_recurse=false)
Output bin; x,[y,[z]],n_sample_points.
void add_sample_point(SamplePoint *new_sample_point_pt, const Vector< double > &zeta_coordinates)
Add a new sample point to RefineableBin.
void output_bin_vertices(std::ofstream &outfile)
Output bin vertices (allowing display of bins as zones).
unsigned total_number_of_sample_points_computed_recursively() const
Compute total number of sample points recursively.
~RefineableBin()
Destructor.
void locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s)
Find sub-GeomObject (finite element) and the local coordinate s within it that contains point with gl...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
static bool Always_fail_elemental_locate_zeta
Boolean flag to make to make locate zeta fail. Used for debugging/ illustration of search procedures.
static double Percentage_offset
Offset of sample point container boundaries beyond max/min coords.
Vector< std::pair< double, double > > Min_and_max_coordinates
Vector of pairs of doubles for min and maximum coordinates. Call: Min_and_max_coordinates[j] gives me...
virtual unsigned & total_number_of_sample_points_visited_during_locate_zeta_from_top_level()
Counter to keep track of how many sample points we've visited during top level call to locate_zeta....
static bool Enable_timing_of_setup
Time setup?
bool ignore_halo_elements_during_locate_zeta_search() const
Ignore halo elements?
bool Ignore_halo_elements_during_locate_zeta_search
Ignore halo elements?
bool use_eulerian_coordinates_during_setup() const
Use Eulerian coordinates (i.e. interpolated_x) rather than zeta itself (i.e. interpolated_zeta) to id...
unsigned & nsample_points_generated_per_element()
"Measure of" number of sample points generated in each element
unsigned Nsample_points_generated_per_element
"Measure of" number of sample points generated in each element
unsigned Total_number_of_sample_points_visited_during_locate_zeta_from_top_level
Counter to keep track of how many sample points we've visited during top level call to locate_zeta.
bool Use_eulerian_coordinates_during_setup
Use Eulerian coordinates (i.e. interpolated_x) rather than zeta itself (i.e. interpolated_zeta) to id...
Mesh * Mesh_pt
Pointer to mesh from whose FiniteElements sample points are created.
const Vector< std::pair< double, double > > & min_and_max_coordinates() const
Vector of pair of doubles for min and maximum coordinates. min (first) and max. (second) coordinates.
double Max_search_radius
Max radius beyond which we stop searching the bin. Initialised to DBL_MAX so keep going until the poi...
double & max_search_radius()
Set maximum search radius for locate zeta. This is initialised do DBL_MAX so we brutally search throu...
static std::ofstream Visited_sample_points_file
File to record sequence of visited sample points in. Used for debugging/ illustration of search proce...
void setup_min_and_max_coordinates()
Helper function to compute the min and max coordinates for the mesh, in each dimension.
static bool Use_equally_spaced_interior_sample_points
Use equally spaced sample points? (otherwise vertices are sampled repeatedly.
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
unsigned element_index_in_mesh() const
Access function to the index of finite element in its mesh.
unsigned sample_point_index_in_element() const
Index of sample point within element.
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...