mesh_as_geometric_object.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// Version 1.0; svn revision $LastChangedRevision: 1132 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-02-23 05:43:26 +0000 (Tue, 23 Feb 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 // Templated functions for MeshAsGeomObject
31 // Config header generated by autoconfig
32 #ifdef HAVE_CONFIG_H
33  #include <oomph-lib-config.h>
34 #endif
35 
36 //oomph-lib headers
37 #include "mesh.h"
39 #include "multi_domain.h"
40 
41 #include <cstdio>
42 
43 namespace oomph
44 {
45 
46  //==========================================================================
47  /// \short Globally set-able (sorry!) flag to indicate that MeshAsGeomObject
48  /// is to use Eulerian coordinates when setting up bin.
49  //==========================================================================
51 
52 
53 //========================================================================
54  /// \short Counter for overall number of bins allocated -- used to
55  /// issue warning if this exceeds a threshhold. (Default assignment
56  /// of 100^DIM bins per MeshAsGeomObject can be a killer if there
57  /// are huge numbers of sub-meshes (e.g. in unstructured FSI).
58 //========================================================================
60 
61 
62 //========================================================================
63  /// \short Total number of bins above which warning is issued.
64  /// (Default assignment of 100^DIM bins per MeshAsGeomObject can
65  /// be a killer if there are huge numbers of sub-meshes (e.g. in
66  /// unstructured FSI).
67 //========================================================================
69  50000000;
70 
71 //========================================================================
72  /// \short Boolean to supppress warnings about large number of bins
73 //========================================================================
75 
76 
77 //========================================================================
78  /// \short Boolean flag to make sure that warning about large number
79 /// of bin cells only gets triggered once.
80 //========================================================================
82 
83 //========================================================================
84  /// \short Number elements/bin that triggers a warning when it gets
85  /// too large, i.e. too many elements per bin, which can lead to very slow
86  /// computations. Sart at 100
87 //========================================================================
89 
90 //========================================================================
91  /// \short Boolean to supppress warnings about small number of bins
92 //========================================================================
94 
95 //========================================================================
96  /// \short Boolean flag to make sure that warning about small number
97 /// of bin cells only gets triggered once.
98 //========================================================================
100 
101 
102 //========================================================================
103 /// Helper function for constructor: Pass the pointer to the mesh,
104 /// communicator and boolean
105 /// to specify whether to calculate coordinate extrema or not
106 /// The dimensions for the GeomObject are read out from the elements and
107 /// nodes of the mesh. Bin dimensions are specified in
108 /// Multi_domain_functions:N{x,y,z}_bin.
109 //========================================================================
111  (Mesh* const &mesh_pt,
112  const bool& compute_extreme_bin_coords)
113  {
114 #ifdef OOMPH_HAS_MPI
115  // Set communicator
116  Communicator_pt=mesh_pt->communicator_pt();
117 #endif
118 
119  // Remember how bin was (or rather: will have been) set up
120  Have_used_eulerian_coordinates_during_setup=
122 
123  //Storage for the Lagrangian and Eulerian dimension
124  int dim[2]={0,0};
125 
126  //Set the Lagrangian dimension from the dimension of the first element
127  //if it exists (if not the Lagrangian dimension will be zero)
128  if(mesh_pt->nelement()!=0)
129  {dim[0] = mesh_pt->finite_element_pt(0)->dim();}
130 
131  //Read out the Eulerian dimension from the first node, if it exists.
132  //(if not the Eulerian dimension will be zero);
133  if(mesh_pt->nnode()!=0)
134  {dim[1] = mesh_pt->node_pt(0)->ndim();}
135 
136  // Need to do an Allreduce to ensure that the dimension is consistent
137  // even when no elements are assigned to a certain processor
138 #ifdef OOMPH_HAS_MPI
139  //Only a problem if the mesh has been distributed
140  if(mesh_pt->is_mesh_distributed())
141  {
142  //Need a non-null communicator
143  if(Communicator_pt!=0)
144  {
145  int n_proc=Communicator_pt->nproc();
146  if (n_proc > 1)
147  {
148  int dim_reduce[2];
149  MPI_Allreduce(&dim,&dim_reduce,2,MPI_INT,
150  MPI_MAX,Communicator_pt->mpi_comm());
151 
152  dim[0] = dim_reduce[0];
153  dim[1] = dim_reduce[1];
154  }
155  }
156  }
157 #endif
158 
159  //Set the Lagrangian and Eulerian dimensions within this geometric object
160  this->set_nlagrangian_and_ndim(static_cast<unsigned>(dim[0]),
161  static_cast<unsigned>(dim[1]));
162 
163  // Create temporary storage for geometric Data (don't count
164  // Data twice!
165  std::set<Data*> tmp_geom_data;
166 
167  double t_cast_finite_element = TimingHelpers::timer();
168 
169  //Copy all the elements in the mesh into local storage
170  //N.B. elements must be able to have a geometric object representation.
171  unsigned n_sub_object = mesh_pt->nelement();
172  Sub_geom_object_pt.resize(n_sub_object);
173  for(unsigned e=0;e<n_sub_object;e++)
174  {
175 
176  // (Try to) cast to a finite element:
177  Sub_geom_object_pt[e]=
178  dynamic_cast<FiniteElement*>(mesh_pt->element_pt(e));
179 
180 #ifdef PARANOID
181  if (Sub_geom_object_pt[e]==0)
182  {
183  std::ostringstream error_message;
184  error_message
185  << "Unable to dynamic cast element: " << std::endl
186  << "into a FiniteElement: GeomObject representation is not possible\n";
187  throw OomphLibError(
188  error_message.str(),
189  OOMPH_CURRENT_FUNCTION,
190  OOMPH_EXCEPTION_LOCATION);
191  }
192 #endif
193 
194  // Add the geometric Data of each element into set
195  unsigned ngeom_data=Sub_geom_object_pt[e]->ngeom_data();
196  for (unsigned i=0;i<ngeom_data;i++)
197  {
198  tmp_geom_data.insert(Sub_geom_object_pt[e]->geom_data_pt(i));
199  }
200  }
201 
202  // Now copy unique geom Data values across into vector
203  unsigned ngeom=tmp_geom_data.size();
204  Geom_data_pt.resize(ngeom);
205  typedef std::set<Data*>::iterator IT;
206  unsigned count=0;
207  for (IT it=tmp_geom_data.begin();it!=tmp_geom_data.end();it++)
208  {
209  Geom_data_pt[count]=*it;
210  count++;
211  }
212 
214  {
215  oomph_info << "CPU for cast_FiniteElement*: "
216  << TimingHelpers::timer()-t_cast_finite_element
217  << std::endl;
218  }
219 
220  // Set storage for minimum and maximum coordinates
221  const unsigned dim_lagrangian = this->nlagrangian();
222  Min_coords.resize(dim_lagrangian);
223  Max_coords.resize(dim_lagrangian);
224 
225  // Get the default parameters for the number of bins in each
226  // dimension from the Multi_domain_functions namespace
230 
231  // Are we computing the extreme bin coordinates here?
232  if (compute_extreme_bin_coords)
233  {
234  double t_compute_extreme_coordinates = TimingHelpers::timer();
235  // Find the maximum and minimum coordinates for the mesh
236  get_min_and_max_coordinates(mesh_pt);
238  {
239  oomph_info << "CPU for computing extreme bin coordinates: "
240  << TimingHelpers::timer()-t_compute_extreme_coordinates
241  << std::endl;
242  }
243  }
244  else
245  {
246 #ifdef PARANOID
247  // Check X_min is less than X_max
249  {
250  std::ostringstream error_stream;
251  error_stream << "Minimum coordinate in the X direction of the bin\n"
252  << "to be used in the multi-domain method is greater\n"
253  << "than or equal to the maximum coordinate.\n"
254  << " X_min=" << Multi_domain_functions::X_min
255  << ", X_max=" << Multi_domain_functions::X_max << "\n";
256  throw OomphLibError
257  (error_stream.str(),
258  OOMPH_CURRENT_FUNCTION,
259  OOMPH_EXCEPTION_LOCATION);
260  }
261  if (dim_lagrangian>=2)
262  {
263  // Check Y_min is less than Y_max
265  {
266  std::ostringstream error_stream;
267  error_stream << "Minimum coordinate in the Y direction of the bin\n"
268  << "to be used in the multi-domain method is greater\n"
269  << "than or equal to the maximum coordinate.\n"
270  << " Y_min=" << Multi_domain_functions::Y_min
271  << ", Y_max=" << Multi_domain_functions::Y_max << "\n";
272  throw OomphLibError
273  (error_stream.str(),
274  OOMPH_CURRENT_FUNCTION,
275  OOMPH_EXCEPTION_LOCATION);
276  }
277  }
278  if (dim_lagrangian==3)
279  {
280  // Check Z_min is less than Z_max
282  {
283  std::ostringstream error_stream;
284  error_stream << "Minimum coordinate in the Z direction of the bin\n"
285  << "to be used in the multi-domain method is greater\n"
286  << "than or equal to the maximum coordinate.\n"
287  << " Z_min=" << Multi_domain_functions::Z_min
288  << ", Z_max=" << Multi_domain_functions::Z_max << "\n";
289  throw OomphLibError
290  (error_stream.str(),
291  OOMPH_CURRENT_FUNCTION,
292  OOMPH_EXCEPTION_LOCATION);
293  }
294  }
295 #endif
296 
297  Min_coords[0] = Multi_domain_functions::X_min-
300  Max_coords[0] = Multi_domain_functions::X_max+
303  if (nlagrangian()>1)
304  {
305  Min_coords[1] = Multi_domain_functions::Y_min-
308  Max_coords[1] = Multi_domain_functions::Y_max+
311  if (nlagrangian()>2)
312  {
313  Min_coords[2] = Multi_domain_functions::Z_min-
316  Max_coords[2] = Multi_domain_functions::Z_max+
319  }
320  }
321  }
322 
323  double t_create_bins_of_object = TimingHelpers::timer();
324 
325  // Create the bin structure, passing in the number of elements in
326  // the mesh
327  create_bins_of_objects(mesh_pt->nelement());
328 
330  {
331 
332  oomph_info << "CPU to create bins of objects [TOTAL]: "
333  << TimingHelpers::timer()-t_create_bins_of_object
334  << std::endl;
335  }
336 
337  // Do some stats
339  {
340  unsigned nbin=0;
341  unsigned max_entry=0;
342  unsigned min_entry=UINT_MAX;
343  unsigned tot_entry=0;
344  unsigned nempty=0;
345  get_fill_stats(nbin,max_entry,min_entry,
346  tot_entry,nempty);
347 
348  double average_entries = double(tot_entry)/double(nbin);
349  oomph_info
350  << "After creation of bins: nbin:("<<nbin<<")"
351  << " nempty:("<<nempty<<")"
352  << " min:("<<min_entry<<")"
353  << " max:("<<max_entry<<")"
354  << " av entries:("<<average_entries<<")"
355  << std::endl;
356 
357 #ifdef OOMPH_HAS_MPI
359  mesh_pt->is_mesh_distributed())
360  {
361  // Get the communicator
362  OomphCommunicator* comm_pt = Communicator_pt;
363  unsigned global_min_entry = 0;
364  unsigned global_max_entry = 0;
365  unsigned global_max_n_empty = 0;
366  double global_max_average_entries = 0;
367 
368  // Get the minimum and maximum time for projection
369  MPI_Reduce(&min_entry, &global_min_entry,
370  1, MPI_UNSIGNED, MPI_MIN, 0, comm_pt->mpi_comm());
371  MPI_Reduce(&max_entry, &global_max_entry,
372  1, MPI_UNSIGNED, MPI_MAX, 0, comm_pt->mpi_comm());
373  MPI_Reduce(&nempty, &global_max_n_empty,
374  1, MPI_UNSIGNED, MPI_MAX, 0, comm_pt->mpi_comm());
375  MPI_Reduce(&average_entries, &global_max_average_entries,
376  1, MPI_DOUBLE, MPI_MAX, 0, comm_pt->mpi_comm());
377 
378  if (comm_pt->my_rank() == 0)
379  {
380  oomph_info
381  << "After creation of bins: (GLOBAL)"
382  << " n_empty:("<<global_max_n_empty<<")"
383  << " min:("<<global_min_entry<<")"
384  << " max:("<<global_max_entry<<")"
385  << " av.entries:("<<global_max_average_entries<<")"
386  << std::endl;
387 
388  }
389 
390  }
391 #endif // #ifdef OOMPH_HAS_MPI
392 
393  }
394 
395  }
396 
397 
398 
399  //=================================================================
400  /// Sort the sampling points in the specified bin by distance from
401  /// sampling point
402  //=================================================================
404  Vector<std::pair<FiniteElement*,
405  Vector<double> > >& sample_point_pairs)
406  {
407 
408  OomphLibWarning("Multi_domain_functions::Sort_bin_entries is currently disabled\n",
409  OOMPH_CURRENT_FUNCTION,
410  OOMPH_EXCEPTION_LOCATION);
411 
412  // // Set current zeta coordinate
413  // Multi_domain_functions::Zeta_coords_for_further_away_comparison=zeta;
414 
415  // // Sort the bin
416  // std::sort(sample_point_pairs.begin(),
417  // sample_point_pairs.end(),
418  // Multi_domain_functions::first_closer_than_second);
419 
420  }
421 
422 
423 
424 
425  //=================================================================
426  /// Get the number of the bin containing the specified coordinate; also
427  /// return the contents of that bin. Bin number is negative if
428  /// the coordinate is outside the bin structure.
429  //=================================================================
430  void MeshAsGeomObject::get_bin(const Vector<double>& zeta, int& bin_number,
431  Vector<std::pair<FiniteElement*,Vector<double> > >& sample_point_pairs)
432  {
433 
434  // Default for point not found
435  sample_point_pairs.clear();
436  bin_number=-1;
437 
438  get_bin(zeta,bin_number);
439  if (bin_number>=0)
440  {
441  sample_point_pairs=Bin_object_coord_pairs[unsigned(bin_number)];
442  }
443  }
444 
445 
446  //=================================================================
447  /// Get the number of the bin containing the specified coordinate.
448  /// Bin number is negative if the coordinate is outside
449  /// the bin structure.
450  //=================================================================
451  void MeshAsGeomObject::get_bin(const Vector<double>& zeta, int& bin_number)
452  {
453 
454  // Default for point not found
455  bin_number=-1;
456 
457  //Get the lagrangian dimension
458  const unsigned n_lagrangian = this->nlagrangian();
459 
460  // Does the zeta coordinate lie within the current bin structure?
461 
462  //Loop over the lagrangian dimension
463  for(unsigned i=0;i<n_lagrangian;i++)
464  {
465  //If the i-th coordinate is less than the minimum
466  if(zeta[i] < Min_coords[i])
467  {
468  return;
469  }
470  //Otherwise coordinate may be bigger than the maximum
471  else if(zeta[i] > Max_coords[i])
472  {
473  return;
474  }
475  }
476 
477  // Use the min and max coords of the bin structure, to find
478  // the bin structure containing the current zeta cooordinate
479 
480  // Initialise for subsequent computation
481  bin_number=0;
482 
483  //Offset for rows/matrices in higher dimensions
484  unsigned multiplier=1;
485  unsigned Nbin[3]={Nbin_x,Nbin_y,Nbin_z};
486 
487  // Loop over the dimension
488  for(unsigned i=0;i<n_lagrangian;i++)
489  {
490  //Find the bin number of the current coordinate
491  unsigned bin_number_i =
492  int(Nbin[i]*((zeta[i]-Min_coords[i])/
493  (Max_coords[i] - Min_coords[i])));
494  //Buffer the case when we are exactly on the edge
495  if(bin_number_i==Nbin[i]) {bin_number_i -= 1;}
496  //Now add to the bin number using the multiplier
497  bin_number += multiplier*bin_number_i;
498  //Increase the current row/matrix multiplier for the next loop
499  multiplier *= Nbin[i];
500  }
501 
502 
503 #ifdef PARANOID
504 
505  // Tolerance for "out of bin" test
506  double tol=1.0e-10;
507 
508  unsigned nvertex=(int)pow(2,n_lagrangian);
509  Vector<Vector<double> > bin_vertex(nvertex);
510  for (unsigned j=0;j<nvertex;j++)
511  {
512  bin_vertex[j].resize(n_lagrangian);
513  }
514  get_bin_vertices(bin_number, bin_vertex);
515  for (unsigned i=0;i<n_lagrangian;i++)
516  {
517  double min_vertex_coord=DBL_MAX;
518  double max_vertex_coord=-DBL_MAX;
519  for (unsigned j=0;j<nvertex;j++)
520  {
521  if (bin_vertex[j][i]<min_vertex_coord) min_vertex_coord=bin_vertex[j][i];
522  if (bin_vertex[j][i]>max_vertex_coord) max_vertex_coord=bin_vertex[j][i];
523  }
524  if (zeta[i]<min_vertex_coord-tol)
525  {
526  std::ostringstream error_message_stream;
527  error_message_stream
528  << "Trouble! " << i << " -th coordinate of sample point, "
529  << zeta[i] << " , isn't actually between limits, "
530  << min_vertex_coord << " and " << max_vertex_coord
531  << " [it's below by more than " << tol << " !] " << std::endl;
532  throw OomphLibError(error_message_stream.str(),
533  OOMPH_CURRENT_FUNCTION,
534  OOMPH_EXCEPTION_LOCATION);
535  }
536 
537  if (zeta[i]>max_vertex_coord+tol)
538  {
539  std::ostringstream error_message_stream;
540  error_message_stream
541  << "Trouble! " << i << " -th coordinate of sample point, "
542  << zeta[i] << " , isn't actually between limits, "
543  << min_vertex_coord << " and " << max_vertex_coord
544  << " [it's above by more than " << tol << "!] " << std::endl;
545  throw OomphLibError(error_message_stream.str(),
546  OOMPH_CURRENT_FUNCTION,
547  OOMPH_EXCEPTION_LOCATION);
548  }
549  }
550 
551 #endif
552 
553  }
554 
555 
556 
557 //========================================================================
558 /// Fill bin by diffusion, populating each empty bin with the same content
559 /// as the first non-empty bin found during a spiral-based search
560 /// up to the specified "radius" (default 1)
561 //========================================================================
563  bin_diffusion_radius)
564  {
565  // Loop over all bins to check if they're empty
566  std::list<unsigned> empty_bins;
567  unsigned nbin=total_nbin();
568  std::vector<bool> was_empty_until_current_iteration(nbin,false);
569  for (unsigned i=0;i<nbin;i++)
570  {
571  if (Bin_object_coord_pairs[i].size()==0)
572  {
573  empty_bins.push_front(i);
574  was_empty_until_current_iteration[i]=true;
575  }
576  }
577 
578  // Now keep processing the empty bins until there are none left
579  unsigned iter=0;
580  Vector<unsigned> newly_filled_bin;
581  while (empty_bins.size()!=0)
582  {
583  newly_filled_bin.clear();
584  newly_filled_bin.reserve(empty_bins.size());
585  for (std::list<unsigned>::iterator it=empty_bins.begin();
586  it!=empty_bins.end();it++)
587  {
588  unsigned bin=(*it);
589 
590  // Look for immediate neighbours within the specified "bin radius"
591  unsigned level=bin_diffusion_radius;
592  Vector<unsigned> neighbour_bin;
593  get_neighbouring_bins_helper(bin,level,neighbour_bin);
594  unsigned n_neigh=neighbour_bin.size();
595 
596 
597  // Find closest pair
598  double min_dist=DBL_MAX;
599  std::pair<FiniteElement*,Vector<double> > closest_pair;
600  for (unsigned i=0;i<n_neigh;i++)
601  {
602  unsigned neigh_bin=neighbour_bin[i];
603 
604  // Only allow to re-populate from bins that were already filled at
605  // previous iteration, otherwise things can progate too fast
606  if (!was_empty_until_current_iteration[neigh_bin])
607  {
608  unsigned nbin_content=Bin_object_coord_pairs[neigh_bin].size();
609  for (unsigned j=0;j<nbin_content;j++)
610  {
611  FiniteElement* el_pt=Bin_object_coord_pairs[neigh_bin][j].first;
612  Vector<double> s(Bin_object_coord_pairs[neigh_bin][j].second);
613  Vector<double> x(2);
614  el_pt->interpolated_x(s,x);
615  // Get minimum distance of sample point from any of the vertices
616  // of current bin
617  double dist=min_distance(bin,x);
618  if (dist<min_dist)
619  {
620  min_dist=dist;
621  closest_pair=Bin_object_coord_pairs[neigh_bin][j];
622  }
623  }
624  }
625  }
626 
627 
628  // Have we filled the bin?
629  if (min_dist!=DBL_MAX)
630  {
632  new_entry.push_back(closest_pair);
633  Bin_object_coord_pairs.set_value(bin,new_entry);
634 
635  // Record that we've filled it.
636  newly_filled_bin.push_back(bin);
637 
638  // Wipe entry without breaking the linked list (Andrew's trick -- nice!)
639  std::list<unsigned>::iterator it_to_be_deleted=it;
640  it--;
641  empty_bins.erase(it_to_be_deleted);
642 
643  }
644  }
645 
646  // Get ready for next iteration on remaining empty bins
647  iter++;
648 
649  // Now update the vector that records which bins were empty up to now
650  unsigned n=newly_filled_bin.size();
651  for (unsigned i=0;i<n;i++)
652  {
653  was_empty_until_current_iteration[newly_filled_bin[i]]=false;
654  }
655  }
656 
657 
658 #ifdef PARANOID
659  // Loop over all bins to check if they're empty
660  nbin=total_nbin();
661  for (unsigned i=0;i<nbin;i++)
662  {
663  if (Bin_object_coord_pairs[i].size()==0)
664  {
665  std::ostringstream error_message_stream;
666  error_message_stream
667  << "Bin " << i << " is still empty\n"
668  << "after trying to fill it by diffusion\n";
669  throw OomphLibError(error_message_stream.str(),
670  OOMPH_CURRENT_FUNCTION,
671  OOMPH_EXCEPTION_LOCATION);
672  }
673  }
674 #endif
675 
676 
677  }
678 
679 //========================================================================
680 /// \short Find the sub geometric object and local coordinate therein that
681 /// corresponds to the intrinsic coordinate zeta. If sub_geom_object_pt=0
682 /// on return from this function, none of the constituent sub-objects
683 /// contain the required coordinate.
684 /// Setting the final bool argument to true means that we only search
685 /// for matching element within a a certain number of "spirals" within
686 /// the bin structure.
687 /// Only search in nonempty bins
688 //========================================================================
690  (const Vector<double>& zeta,GeomObject*& sub_geom_object_pt,
691  Vector<double>& s, const bool &called_within_spiral)
692  {
693  // Initialise return to null -- if it's still null when we're
694  // leaving we've failed!
695  sub_geom_object_pt=0;
696 
697  //Get the lagrangian dimension
698  const unsigned n_lagrangian = this->nlagrangian();
699 
700  // Does the zeta coordinate lie within the current bin structure?
701 
702  //Loop over the lagrangian dimension
703  for(unsigned i=0;i<n_lagrangian;i++)
704  {
705  //If the i-th coordinate is less than the minimum
706  if(zeta[i] < Min_coords[i])
707  {
708  return;
709  }
710  //Otherwise coordinate may be bigger than the maximum
711  else if(zeta[i] > Max_coords[i])
712  {
713  return;
714  }
715  }
716 
717  // Use the min and max coords of the bin structure, to find
718  // the bin structure containing the current zeta cooordinate
719  unsigned bin_number=0;
720  //Offset for rows/matrices in higher dimensions
721  unsigned multiplier=1;
722  unsigned Nbin[3]={Nbin_x,Nbin_y,Nbin_z};
723 
724  // Loop over the dimension
725  for(unsigned i=0;i<n_lagrangian;i++)
726  {
727  //Find the bin number of the current coordinate
728  unsigned bin_number_i =
729  int(Nbin[i]*((zeta[i]-Min_coords[i])/
730  (Max_coords[i] - Min_coords[i])));
731  //Buffer the case when we are exactly on the edge
732  if(bin_number_i==Nbin[i]) {bin_number_i -= 1;}
733  //Now add to the bin number using the multiplier
734  bin_number += multiplier*bin_number_i;
735  //Increase the current row/matrix multiplier for the next loop
736  multiplier *= Nbin[i];
737  }
738 
739  if (called_within_spiral)
740  {
741  // Loop over spirals that are to be visited in one go
742  Vector<unsigned> neighbour_bin;
743  for (unsigned i_level=current_min_spiral_level();
744  i_level<=current_max_spiral_level();i_level++)
745  {
746  // Only add bins that have elements if finding zeta values for
747  // projection
748  bool only_use_filled_bins=false;
750  {
751  only_use_filled_bins=true;
752  }
753  // Call helper function to find the neighbouring bins at this
754  // level
755  get_neighbouring_bins_helper(bin_number,i_level,neighbour_bin,
756  only_use_filled_bins);
757  }
758 
759  // Total number of bins to be visited
760  unsigned n_nbr_bin=neighbour_bin.size();
761 
762  // Set bool for finding zeta
763  bool found_zeta=false;
764  for (unsigned i_nbr=0;i_nbr<n_nbr_bin;i_nbr++)
765  {
766  // Only search if bin is within the max. radius
767  if (min_distance(neighbour_bin[i_nbr],zeta)<Max_search_radius)
768  {
769  // Get the number of element-sample point pairs in this bin
770  unsigned n_sample=
771  Bin_object_coord_pairs[neighbour_bin[i_nbr]].size();
772 
773  // Don't do anything if this bin has no sample points
774  if (n_sample>0)
775  {
776 
777  // Sort the bin if required
779  {
780  OomphLibWarning("Multi_domain_functions::Sort_bin_entries is currently disabled\n",
781  OOMPH_CURRENT_FUNCTION,
782  OOMPH_EXCEPTION_LOCATION);
783  //sort_the_bin(zeta,Bin_object_coord_pairs[neighbour_bin[i_nbr]]);
784  }
785 
786  for (unsigned i_sam=0;i_sam<n_sample;i_sam++)
787  {
788  // Get the element
789  FiniteElement* el_pt=Bin_object_coord_pairs
790  [neighbour_bin[i_nbr]][i_sam].first;
791 
792  // Get the local coordinate
793  s=Bin_object_coord_pairs[neighbour_bin[i_nbr]][i_sam].second;
794 
795  // Use this coordinate as the initial guess
796  bool use_coordinate_as_initial_guess=true;
797 
798  // Attempt to find zeta within a sub-object
799  el_pt->locate_zeta(zeta,sub_geom_object_pt,s,
800  use_coordinate_as_initial_guess);
801 
802 #ifdef OOMPH_HAS_MPI
804  {
805  // Dynamic cast the result to a FiniteElement
806  FiniteElement* test_el_pt=
807  dynamic_cast<FiniteElement*>(sub_geom_object_pt);
808  if (test_el_pt!=0)
809  {
810  // We only want to exit if this is a non-halo element
811  if (test_el_pt->is_halo()) {sub_geom_object_pt=0;}
812  }
813  }
814  else
815  {
817  {
818  // Dynamic cast the result to a FiniteElement
819  FiniteElement* test_el_pt=
820  dynamic_cast<FiniteElement*>(sub_geom_object_pt);
821  if (test_el_pt!=0)
822  {
823  // We only want to exit if this is a non-halo element
824  if (test_el_pt->is_halo()) {sub_geom_object_pt=0;}
825  }
826  }
827  }
828 #endif
829 
830  // If the FiniteElement is non-halo and has been located, exit
831  if (sub_geom_object_pt!=0)
832  {
833  found_zeta=true;
834  break;
835  }
836  } // end loop over sample points
837  }
838 
839  if (found_zeta)
840  {
841  break;
842  }
843 
844  } //end of don't search if outside search radius
845  else
846  {
847  //oomph_info << "Terminating: Outside search radius [1]\n";
848  }
849  } // end loop over bins at this level
850 
851  }
852  else
853  {
854  // Not called from within a spiral procedure
855  // (i.e. the loop in multi_domain.h), so do the spiralling here
856 
857  // Loop over all levels... maximum of N*_bin
858  unsigned n_level=Nbin[0];
859  for(unsigned i=1;i<n_lagrangian;i++)
860  {
861  if(n_level < Nbin[i]) {n_level = Nbin[i];}
862  }
863 
864  // Also limit to max. spiral level
865  n_level=std::min(n_level,max_spiral_level());
866 
867  // Set bool for finding zeta
868  bool found_zeta=false;
869  for (unsigned i_level=0;i_level<n_level;i_level++)
870  {
871  // Call helper function to find the neighbouring bins at this
872  // level
873  Vector<unsigned> neighbour_bin;
874  // Only add bins that have elements if finding zeta values for
875  // projection
876  bool only_use_filled_bins=false;
878  {
879  only_use_filled_bins=true;
880  }
881  get_neighbouring_bins_helper(bin_number,i_level,neighbour_bin,
882  only_use_filled_bins);
883  unsigned n_nbr_bin=neighbour_bin.size();
884 
885  // Loop over neighbouring bins
886  for (unsigned i_nbr=0;i_nbr<n_nbr_bin;i_nbr++)
887  {
888 
889  // Only search if bin is within the max. radius
890  if (min_distance(neighbour_bin[i_nbr],zeta)<Max_search_radius)
891  {
892  // Get the number of element-sample point pairs in this bin
893  unsigned n_sample=
894  Bin_object_coord_pairs[neighbour_bin[i_nbr]].size();
895 
896  // Don't do anything if this bin has no sample points
897  if (n_sample>0)
898  {
899 
900  // Sort the bin if required
902  {
903  OomphLibWarning("Multi_domain_functions::Sort_bin_entries is currently disabled\n",
904  OOMPH_CURRENT_FUNCTION,
905  OOMPH_EXCEPTION_LOCATION);
906  //sort_the_bin(zeta,Bin_object_coord_pairs[neighbour_bin[i_nbr]]);
907  }
908 
909  for (unsigned i_sam=0;i_sam<n_sample;i_sam++)
910  {
911  // Get the element
912  FiniteElement* el_pt=Bin_object_coord_pairs
913  [neighbour_bin[i_nbr]][i_sam].first;
914 
915  // Get the local coordinate
916  s=Bin_object_coord_pairs[neighbour_bin[i_nbr]][i_sam].second;
917 
918  // Use this coordinate as the initial guess in locate_zeta
919  bool use_coordinate_as_initial_guess=true;
920 
921  // Attempt to loacte the correct sub-object
922  el_pt->locate_zeta(zeta,sub_geom_object_pt,s,
923  use_coordinate_as_initial_guess);
924 
925  // If it was found then break
926  if (sub_geom_object_pt!=0)
927  {
928  found_zeta=true;
929  break;
930  }
931  } // end loop over sample points
932  }
933 
934  // Break out of the bin loop if locate was successful
935  if (found_zeta)
936  {
937  break;
938  }
939  } //end of don't search if outside search radius
940  else
941  {
942  //oomph_info << "Terminating: Outside search radius [2]\n";
943  }
944 
945  } // end loop over bins at this level
946 
947  // Break out of the spiral loop if locate was successful
948  if (found_zeta)
949  {
950  break;
951  }
952 
953  } // end loop over levels
954 
955  } // end if (called_within_spiral)
956 
957  }
958 
959 //========================================================================
960 /// \short Find the sub geometric object and local coordinate therein that
961 /// corresponds to the intrinsic coordinate zeta. If sub_geom_object_pt=0
962 /// on return from this function, none of the constituent sub-objects
963 /// contain the required coordinate.
964 /// Setting the final bool argument to true means that we only search
965 /// for matching element within a a certain number of "spirals" within
966 /// the bin structure.
967 //========================================================================
969  (const Vector<double>& zeta,GeomObject*& sub_geom_object_pt,
970  Vector<double>& s, const bool &called_within_spiral)
971  {
972  // Initialise return to null -- if it's still null when we're
973  // leaving we've failed!
974  sub_geom_object_pt=0;
975 
976  //Get the lagrangian dimension
977  const unsigned n_lagrangian = this->nlagrangian();
978 
979  // Does the zeta coordinate lie within the current bin structure?
980 
981  //Loop over the lagrangian dimension
982  for(unsigned i=0;i<n_lagrangian;i++)
983  {
984  //If the i-th coordinate is less than the minimum
985  if(zeta[i] < Min_coords[i])
986  {
987  return;
988  }
989  //Otherwise coordinate may be bigger than the maximum
990  else if(zeta[i] > Max_coords[i])
991  {
992  return;
993  }
994  }
995 
996  // Use the min and max coords of the bin structure, to find
997  // the bin structure containing the current zeta cooordinate
998  unsigned bin_number=0;
999  //Offset for rows/matrices in higher dimensions
1000  unsigned multiplier=1;
1001  unsigned Nbin[3]={Nbin_x,Nbin_y,Nbin_z};
1002 
1003  // Loop over the dimension
1004  for(unsigned i=0;i<n_lagrangian;i++)
1005  {
1006  //Find the bin number of the current coordinate
1007  unsigned bin_number_i =
1008  int(Nbin[i]*((zeta[i]-Min_coords[i])/
1009  (Max_coords[i] - Min_coords[i])));
1010  //Buffer the case when we are exactly on the edge
1011  if(bin_number_i==Nbin[i]) {bin_number_i -= 1;}
1012  //Now add to the bin number using the multiplier
1013  bin_number += multiplier*bin_number_i;
1014  //Increase the current row/matrix multiplier for the next loop
1015  multiplier *= Nbin[i];
1016  }
1017 
1018  // Visit all levels
1019  const unsigned my_max_spiral_level = UINT_MAX;
1020  if (called_within_spiral)
1021  {
1022  for (unsigned i_level=0;i_level<=my_max_spiral_level;i_level++)
1023  //for (unsigned i_level=current_min_spiral_level();
1024  // i_level<=current_max_spiral_level();i_level++)
1025  {
1026  // Loop over spirals that are to be visited in one go
1027  Vector<unsigned> neighbour_bin;
1028  // Only add bins that have elements
1029  const bool only_use_filled_bins=true;
1030  // Call helper function to find the neighbouring bins at this level
1031  get_neighbouring_bins_helper(bin_number,i_level,neighbour_bin,
1032  only_use_filled_bins);
1033 
1034  // Total number of bins to be visited
1035  unsigned n_nbr_bin=neighbour_bin.size();
1036 
1037  // Set bool for finding zeta
1038  bool found_zeta=false;
1039  for (unsigned i_nbr=0;i_nbr<n_nbr_bin;i_nbr++)
1040  {
1041  // Only search if bin is within the max. radius
1042  if (min_distance(neighbour_bin[i_nbr],zeta)<Max_search_radius)
1043  {
1044  // Get the number of element-sample point pairs in this bin
1045  unsigned n_sample=
1046  Bin_object_coord_pairs[neighbour_bin[i_nbr]].size();
1047 
1048  // Don't do anything if this bin has no sample points
1049  if (n_sample>0)
1050  {
1051  // Sort the bin if required
1053  {
1054 // if (Multi_domain_functions::Doc_timings)
1055 // {
1056 // double t_sort_bin=TimingHelpers::timer();
1057 // sort_the_bin(zeta,
1058 // Bin_object_coord_pairs[neighbour_bin[i_nbr]]);
1059 // Multi_domain_functions::
1060 // Total_time_for_sorting_elements_in_bins+=
1061 // TimingHelpers::timer()-t_sort_bin;
1062 // }
1063 // else
1064 // {
1065 // sort_the_bin(zeta,
1066 // Bin_object_coord_pairs[neighbour_bin[i_nbr]]);
1067 // }
1068  OomphLibWarning("Multi_domain_functions::Sort_bin_entries is currently disabled\n",
1069  OOMPH_CURRENT_FUNCTION,
1070  OOMPH_EXCEPTION_LOCATION);
1071  }
1072 
1073  for (unsigned i_sam=0;i_sam<n_sample;i_sam++)
1074  {
1075  // Get the element
1076  FiniteElement* el_pt=Bin_object_coord_pairs
1077  [neighbour_bin[i_nbr]][i_sam].first;
1078 
1079  // Get the local coordinate
1080  s=Bin_object_coord_pairs[neighbour_bin[i_nbr]][i_sam].second;
1081 
1082  // Use this coordinate as the initial guess
1083  bool use_coordinate_as_initial_guess=true;
1084 
1085  // Attempt to find zeta within a sub-object
1086  el_pt->locate_zeta(zeta,sub_geom_object_pt,s,
1087  use_coordinate_as_initial_guess);
1088 
1089 #ifdef OOMPH_HAS_MPI
1090  // Allow halo elements as external elements when doing
1091  // the projection
1093  {
1094  // Dynamic cast the result to a FiniteElement
1095  FiniteElement* test_el_pt=
1096  dynamic_cast<FiniteElement*>(sub_geom_object_pt);
1097  if (test_el_pt!=0)
1098  {
1099  // We only want to exit if this is a non-halo element
1100  if (test_el_pt->is_halo()) {sub_geom_object_pt=0;}
1101  }
1102  }
1103 #endif
1104 
1105  // If the FiniteElement is non-halo and has been located, exit
1106  if (sub_geom_object_pt!=0)
1107  {
1108  found_zeta=true;
1109  break;
1110  }
1111  } // end loop over sample points
1112  }
1113 
1114  if (found_zeta)
1115  {
1116  break;
1117  }
1118 
1119  } //end of don't search if outside search radius
1120  else
1121  {
1122  //oomph_info << "Terminating: Outside search radius [1]\n";
1123  }
1124 
1125  } // end loop over bins at this level
1126 
1127  if (found_zeta)
1128  {
1129  break;
1130  }
1131 
1132  } // end loop over levels in this chunk
1133 
1134  }
1135  else
1136  {
1137  // Not called from within a spiral procedure
1138  // (i.e. the loop in multi_domain.h), so do the spiralling here
1139 
1140  // Loop over all levels... maximum of N*_bin
1141  unsigned n_level=Nbin[0];
1142  for(unsigned i=1;i<n_lagrangian;i++)
1143  {
1144  if(n_level < Nbin[i]) {n_level = Nbin[i];}
1145  }
1146 
1147  // Also limit to max. spiral level
1148  n_level=std::min(n_level,max_spiral_level());
1149 
1150  // Set bool for finding zeta
1151  bool found_zeta=false;
1152  for (unsigned i_level=0;i_level<n_level;i_level++)
1153  {
1154  // Loop over spirals that are to be visited in one go
1155  Vector<unsigned> neighbour_bin;
1156  // Only add bins that have elements
1157  const bool only_use_filled_bins=true;
1158  // Call helper function to find the neighbouring bins at this level
1159  get_neighbouring_bins_helper(bin_number,i_level,neighbour_bin,
1160  only_use_filled_bins);
1161  unsigned n_nbr_bin=neighbour_bin.size();
1162 
1163  // Loop over neighbouring bins
1164  for (unsigned i_nbr=0;i_nbr<n_nbr_bin;i_nbr++)
1165  {
1166  // Only search if bin is within the max. radius
1167  if (min_distance(neighbour_bin[i_nbr],zeta)<Max_search_radius)
1168  {
1169  // Get the number of element-sample point pairs in this bin
1170  unsigned n_sample=
1171  Bin_object_coord_pairs[neighbour_bin[i_nbr]].size();
1172 
1173  // Don't do anything if this bin has no sample points
1174  if (n_sample>0)
1175  {
1176  // Sort the bin if required
1178  {
1179 // sort_the_bin(zeta,Bin_object_coord_pairs[neighbour_bin[i_nbr]]);
1180  OomphLibWarning("Multi_domain_functions::Sort_bin_entries is currently disabled\n",
1181  OOMPH_CURRENT_FUNCTION,
1182  OOMPH_EXCEPTION_LOCATION);
1183  }
1184 
1185  for (unsigned i_sam=0;i_sam<n_sample;i_sam++)
1186  {
1187  // Get the element
1188  FiniteElement* el_pt=Bin_object_coord_pairs
1189  [neighbour_bin[i_nbr]][i_sam].first;
1190 
1191  // Get the local coordinate
1192  s=Bin_object_coord_pairs[neighbour_bin[i_nbr]][i_sam].second;
1193 
1194  // Use this coordinate as the initial guess in locate_zeta
1195  bool use_coordinate_as_initial_guess=true;
1196 
1197  // Attempt to loacte the correct sub-object
1198  el_pt->locate_zeta(zeta,sub_geom_object_pt,s,
1199  use_coordinate_as_initial_guess);
1200 
1201  // If it was found then break
1202  if (sub_geom_object_pt!=0)
1203  {
1204  found_zeta=true;
1205  break;
1206  }
1207  } // end loop over sample points
1208  }
1209 
1210  // Break out of the bin loop if locate was successful
1211  if (found_zeta)
1212  {
1213  break;
1214  }
1215  } //end of don't search if outside search radius
1216  else
1217  {
1218  //oomph_info << "Terminating: Outside search radius [2]\n";
1219  }
1220 
1221  } // end loop over bins at this level
1222 
1223  // Break out of the spiral loop if locate was successful
1224  if (found_zeta)
1225  {
1226  break;
1227  }
1228 
1229  } // end loop over levels
1230 
1231  } // end if (called_within_spiral)
1232 
1233  }
1234 
1235 //========================================================================
1236 ///Get the min and max coordinates for the mesh, in each dimension
1237 //========================================================================
1238  void MeshAsGeomObject::
1240  {
1241  //Get the lagrangian dimension
1242  int n_lagrangian = this->nlagrangian();
1243 
1244  // Storage locally (i.e. in parallel on each processor)
1245  // for the minimum and maximum coordinates
1246  double zeta_min_local[n_lagrangian]; double zeta_max_local[n_lagrangian];
1247  for(int i=0;i<n_lagrangian;i++)
1248  {zeta_min_local[i] = DBL_MAX; zeta_max_local[i] = -DBL_MAX;}
1249 
1250  // Loop over the elements of the mesh
1251  unsigned n_el=mesh_pt->nelement();
1252  for (unsigned e=0;e<n_el;e++)
1253  {
1254  FiniteElement* el_pt = mesh_pt->finite_element_pt(e);
1255 
1256  // Get the number of vertices (nplot=2 does the trick)
1257  unsigned n_plot=2;
1258  unsigned n_plot_points=el_pt->nplot_points(n_plot);
1259 
1260  // Loop over the number of plot points
1261  for (unsigned iplot=0;iplot<n_plot_points;iplot++)
1262  {
1263  Vector<double> s_local(n_lagrangian);
1264  Vector<double> zeta_global(n_lagrangian);
1265 
1266  // Get the local s
1267  el_pt->get_s_plot(iplot,n_plot,s_local);
1268 
1269  // Now interpolate to global coordinates
1271  {
1272  el_pt->interpolated_x(s_local,zeta_global);
1273  }
1274  else
1275  {
1276  el_pt->interpolated_zeta(s_local,zeta_global);
1277  }
1278 
1279  // Check the max and min in each direction
1280  for(int i=0;i<n_lagrangian;i++)
1281  {
1282  //Is the coordinate less than the minimum?
1283  if(zeta_global[i] < zeta_min_local[i])
1284  {zeta_min_local[i] = zeta_global[i];}
1285  //Is the coordinate bigger than the maximum?
1286  if(zeta_global[i] > zeta_max_local[i])
1287  {zeta_max_local[i] = zeta_global[i];}
1288  }
1289  }
1290  }
1291 
1292  // Global extrema - in parallel, need to get max/min across all processors
1293  double zeta_min[n_lagrangian]; double zeta_max[n_lagrangian];
1294  for(int i=0;i<n_lagrangian;i++) {zeta_min[i] = 0.0; zeta_max[i] = 0.0;}
1295 
1296 #ifdef OOMPH_HAS_MPI
1297  // If the mesh has been distributed and we want consistent bins
1298  // across all processors
1299  if ( mesh_pt->is_mesh_distributed() &&
1301  {
1302  // .. we need a non-null communicator!
1303  if (Communicator_pt!=0)
1304  {
1305  int n_proc=Communicator_pt->nproc();
1306  if (n_proc>1)
1307  {
1308  //Get the minima and maxima over all processors
1309  MPI_Allreduce(zeta_min_local,zeta_min,n_lagrangian,MPI_DOUBLE,MPI_MIN,
1310  Communicator_pt->mpi_comm());
1311  MPI_Allreduce(zeta_max_local,zeta_max,n_lagrangian,MPI_DOUBLE,MPI_MAX,
1312  Communicator_pt->mpi_comm());
1313  }
1314  }
1315  else // Null communicator - throw an error
1316  {
1317  std::ostringstream error_message_stream;
1318  error_message_stream
1319  << "Communicator not set for a MeshAsGeomObject\n"
1320  << "that was created from a distributed Mesh";
1321  throw OomphLibError(error_message_stream.str(),
1322  OOMPH_CURRENT_FUNCTION,
1323  OOMPH_EXCEPTION_LOCATION);
1324  }
1325  }
1326  else // If the mesh hasn't been distributed then the
1327  // max and min are the same on all processors
1328  {
1329  for(int i=0;i<n_lagrangian;i++)
1330  {
1331  zeta_min[i] = zeta_min_local[i];
1332  zeta_max[i] = zeta_max_local[i];
1333  }
1334  }
1335 #else // If we're not using MPI then the mesh can't be distributed
1336  for(int i=0;i<n_lagrangian;i++)
1337  {
1338  zeta_min[i] = zeta_min_local[i];
1339  zeta_max[i] = zeta_max_local[i];
1340  }
1341 #endif
1342 
1343  // Decrease/increase min and max to allow for any overshoot in
1344  // meshes that may move around
1345  // There's no point in doing this for DIM_LAGRANGIAN==1
1346  double percentage_offset=Multi_domain_functions::Percentage_offset;
1347  for(int i=0;i<n_lagrangian;i++)
1348  {
1349  double length = zeta_max[i] - zeta_min[i];
1350  zeta_min[i] -= ((percentage_offset/100.0)*length);
1351  zeta_max[i] += ((percentage_offset/100.0)*length);
1352  }
1353 
1354  // Set the entries as the Min/Max_coords vector
1355  for(int i=0;i<n_lagrangian;i++)
1356  {
1357  Min_coords[i] = zeta_min[i];
1358  Max_coords[i] = zeta_max[i];
1359  }
1360  }
1361 
1362 
1363 
1364 //========================================================================
1365  /// Provide some stats on the fill level of the associated bin
1366 //========================================================================
1367  void MeshAsGeomObject::get_fill_stats(unsigned& n_bin,
1368  unsigned& max_n_entry,
1369  unsigned& min_n_entry,
1370  unsigned& tot_n_entry,
1371  unsigned& n_empty) const
1372  {
1373  // Total number of bins
1374  n_bin=total_nbin();
1375  n_empty=n_bin-Bin_object_coord_pairs.nnz();
1376 
1377  // Get pointer to map-based representation
1378  const std::map<unsigned,Vector<std::pair<FiniteElement*,Vector<double> > > >*
1379  map_pt=Bin_object_coord_pairs.map_pt();
1380 
1381  // Initialise
1382  max_n_entry=0;
1383  min_n_entry=UINT_MAX;
1384  tot_n_entry=0;
1385 
1386  // Do stats
1387  typedef std::map<unsigned,Vector<
1388  std::pair<FiniteElement*,Vector<double> > > >::const_iterator IT;
1389  for (IT it=map_pt->begin();it!=map_pt->end();it++)
1390  {
1391  unsigned nentry=(*it).second.size();
1392  if (nentry>max_n_entry) max_n_entry=nentry;
1393  if (nentry<min_n_entry) min_n_entry=nentry;
1394  tot_n_entry+=nentry;
1395  }
1396 
1397  }
1398 
1399 
1400 
1401 //========================================================================
1402 ///Initialise the "bin" structure for locating coordinates
1403 ///The number of elements in the mesh is passed in only for an error check
1404 ///about the number of elements / bin.
1405 //========================================================================
1407  const unsigned long &n_mesh_element)
1408  {
1409  //Store the lagrangian dimension
1410  const unsigned n_lagrangian = this->nlagrangian();
1411 
1412  // Output message regarding bin structure setup if required
1414  {
1415  oomph_info << "============================================" << std::endl;
1416  oomph_info << " MeshAsGeomObject: set up bin search with:" << std::endl;
1417  oomph_info << " Nbin_x=" << Nbin_x << " ";
1418  if (n_lagrangian>=2)
1419  {
1420  oomph_info << "Nbin_y=" << Nbin_y << " ";
1421  }
1422  if (n_lagrangian==3)
1423  {
1424  oomph_info << "Nbin_z=" << Nbin_z;
1425  }
1426  oomph_info << std::endl;
1427  oomph_info << " Xminmax=" << Min_coords[0] << " " << Max_coords[0]
1428  << " ";
1429  if (n_lagrangian>=2)
1430  {
1431  oomph_info << "Yminmax=" << Min_coords[1] << " " << Max_coords[1]
1432  << " ";
1433  }
1434  if (n_lagrangian==3)
1435  {
1436  oomph_info << "Zminmax=" << Min_coords[2] << " " << Max_coords[2]
1437  << " ";
1438  }
1439  oomph_info << std::endl;
1440  oomph_info << "============================================" << std::endl;
1441  }
1442 
1443  /// Flush all objects out of the bin structure
1445 
1446  // Temporary storage in bin
1447  std::map<unsigned,Vector<std::pair<FiniteElement*,Vector<double> > > >
1448  tmp_bin_object_coord_pairs;
1449 
1450  unsigned Nbin[3] ={Nbin_x, Nbin_y, Nbin_z};
1451 
1452  // Total number of bins
1453  unsigned ntotalbin=total_nbin();
1454 
1455  // Initialise the structure that keep tracks of the occupided bins
1456  Bin_object_coord_pairs.initialise(ntotalbin);
1457 
1458  // Issue warning about small number of bins
1460  {
1461  //Calculate the (nearest integer) number of elements per bin
1462  unsigned elements_per_bin = n_mesh_element / ntotalbin;
1463  //If it is above the threshold then issue a warning
1464  if (elements_per_bin > Threshold_for_elements_per_bin_warning)
1465  {
1467  {
1469  std::ostringstream warn_message;
1470  warn_message
1471  << "The average (integer) number of elements per bin is \n\n"
1472  << elements_per_bin
1473  << ", which is more than \n\n"
1474  << " MeshAsGeomObject::Threshold_for_elements_per_bin_warning="
1476  << "\n\nIf the lookup seems slow (and you have the memory), consider increasing\n"
1477  << "Multi_domain_functions::N{x,y,z}_bin from their current\n"
1478  << "values of { "
1481  << Multi_domain_functions::Nz_bin << " }.\n"
1482  << "\nNOTE: You can suppress this warning by increasing\n\n"
1483  << " MeshAsGeomObject::Threshold_for_elements_per_bin_warning\n\n"
1484  << "or by setting \n\n MeshAsGeomObject::Suppress_warning_about_small_number_of_bins\n\n"
1485  << "to true (both are public static data).\n\n";
1487  warn_message.str(),
1488  "MeshAsGeomObject::create_bins_of_objects()",
1489  OOMPH_EXCEPTION_LOCATION);
1490  }
1491  }
1492  }
1493 
1494 
1495 
1496  // Increase overall counter
1497  Total_nbin_cells_counter+=ntotalbin;
1498 
1499  // Issue warning?
1501  {
1503  {
1505  {
1507  std::ostringstream warn_message;
1508  warn_message
1509  << "Have allocated \n\n"
1510  << " MeshAsGeomObject::Total_nbin_cells_counter="
1512  << "\n\nbin cells, which is more than \n\n"
1513  << " MeshAsGeomObject::Threshold_for_total_bin_cell_number_warning="
1515  << "\n\nIf you run out of memory, consider reducing\n"
1516  << "Multi_domain_functions::N{x,y,z}_bin from their current\n"
1517  << "values of { "
1520  << Multi_domain_functions::Nz_bin << " }.\n"
1521  << "\nNOTE: You can suppress this warning by increasing\n\n"
1522  << " MeshAsGeomObject::Threshold_for_total_bin_cell_number_warning\n\n"
1523  << "or by setting \n\n MeshAsGeomObject::Suppress_warning_about_large_total_number_of_bins\n\n"
1524  << "to true (both are public static data).\n\n"
1525  << "NOTE: I'll only issues this warning once -- total number of\n"
1526  << "bins may grow yet further!\n";
1528  warn_message.str(),
1529  "MeshAsGeomObject::create_bins_of_objects()",
1530  OOMPH_EXCEPTION_LOCATION);
1531  }
1532  }
1533  }
1534 
1535  ///Loop over subobjects (elements) to decide which bin they belong in...
1536  unsigned n_sub=Sub_geom_object_pt.size();
1537 
1538  // Some stats
1540  {
1541  oomph_info << "There are " << n_sub << " element[s] to be put into bins"
1542  << std::endl << std::endl;
1543  }
1544 
1545  double t_put_elements_in_bins = TimingHelpers::timer();
1546 
1547  for (unsigned e=0;e<n_sub;e++)
1548  {
1549  // Cast to the element (sub-object) first
1550  FiniteElement* el_pt=dynamic_cast<FiniteElement*>(Sub_geom_object_pt[e]);
1551 
1552  // Get specified number of points within the element
1553  unsigned n_plot_points=
1555 
1556  for (unsigned iplot=0;iplot<n_plot_points;iplot++)
1557  {
1558  // Storage for local and global coordinates
1559  Vector<double> local_coord(n_lagrangian,0.0);
1560  Vector<double> global_coord(n_lagrangian,0.0);
1561 
1562  // Get local coordinate and interpolate to global
1563  el_pt->get_s_plot(iplot,
1565 
1566  // Now get appropriate global coordinate
1568  {
1569  el_pt->interpolated_x(local_coord,global_coord);
1570  }
1571  else
1572  {
1573  el_pt->interpolated_zeta(local_coord,global_coord);
1574  }
1575 
1576  //Which bin are the global coordinates in?
1577  unsigned bin_number=0;
1578  unsigned multiplier=1;
1579  // Loop over the dimension
1580  for(unsigned i=0;i<n_lagrangian;i++)
1581  {
1582 #ifdef PARANOID
1583  if ((global_coord[i]<Min_coords[i])||
1584  (global_coord[i]>Max_coords[i]))
1585  {
1586  std::ostringstream error_message;
1587  error_message
1588  << "Bin sample point " << iplot << " in element " << e << "\n"
1589  << "is outside bin limits in coordinate direction " << i << ":\n"
1590  << "Sample point coordinate: " << global_coord[i] << "\n"
1591  << "Max bin coordinate : " << Max_coords[i] << "\n"
1592  << "Min bin coordinate : " << Min_coords[i] << "\n"
1593  << "You should either setup the bin boundaries manually\n"
1594  << "or increase the percentage offset by which the\n"
1595  << "automatically computed bin limits are increased \n"
1596  << "beyond their sampled max/mins. This is defined in\n"
1597  << "the (public) namespace member\n\n"
1598  << "Multi_domain_functions::Percentage_offset \n\n which \n"
1599  << "currently has the value: "
1601  throw OomphLibError(
1602  error_message.str(),
1603  OOMPH_CURRENT_FUNCTION,
1604  OOMPH_EXCEPTION_LOCATION);
1605  }
1606 
1607 #endif
1608  unsigned bin_number_i =
1609  int(Nbin[i]*(
1610  (global_coord[i] - Min_coords[i])/
1611  (Max_coords[i] - Min_coords[i])));
1612 
1613  //Buffer the case when the global coordinate is the maximum
1614  //value
1615  if(bin_number_i==Nbin[i]) {bin_number_i -= 1;}
1616 
1617  //Add to the bin number
1618  bin_number += multiplier*bin_number_i;
1619 
1620  //Sort out the multiplier
1621  multiplier *= Nbin[i];
1622  }
1623 
1624  // //Add element-sample local coord pair to the calculated bin
1625  // Bin_object_coord_pairs[bin_number].push_back
1626  // (std::make_pair(el_pt,local_coord));
1627 
1628 
1629  //Add element-sample local coord pair to the calculated bin
1630  tmp_bin_object_coord_pairs[bin_number].push_back
1631  (std::make_pair(el_pt,local_coord));
1632 
1633  }
1634  }
1635 
1637  {
1638  oomph_info << "CPU for adding elements in bins "
1639  << TimingHelpers::timer()-t_put_elements_in_bins
1640  << std::endl;
1641  }
1642 
1643  // Finally copy filled vectors across -- wiping memory from temporary
1644  // map as we go along is good and wouldn't be possible if we
1645  // filled the SparseVector's internal map within that class.
1646  typedef std::map<unsigned,Vector<
1647  std::pair<FiniteElement*,Vector<double> > > >::iterator IT;
1648  for (IT it=tmp_bin_object_coord_pairs.begin();
1649  it!=tmp_bin_object_coord_pairs.end();it++)
1650  {
1651  Bin_object_coord_pairs.set_value((*it).first,
1652  (*it).second);
1653  // Make space immediately
1654  (*it).second.clear();
1655  }
1656 
1657  }
1658 
1659 //========================================================================
1660 /// Output bins
1661 //========================================================================
1662  void MeshAsGeomObject::output_bins(std::ofstream& outfile)
1663  {
1664  // Spatial dimension of bin
1665  const unsigned n_lagrangian = this->nlagrangian();
1666 
1667  unsigned nbin_x=Nbin_x;
1668  unsigned nbin_y=1;
1669  if (n_lagrangian>1) nbin_y=Nbin_y;
1670  unsigned nbin_z=1;
1671  if (n_lagrangian>2) nbin_z=Nbin_z;
1672 
1673  unsigned b=0;
1674  for (unsigned iz=0;iz<nbin_z;iz++)
1675  {
1676  for (unsigned iy=0;iy<nbin_y;iy++)
1677  {
1678  for (unsigned ix=0;ix<nbin_x;ix++)
1679  {
1680  unsigned nentry=Bin_object_coord_pairs[b].size();
1681  for (unsigned e=0;e<nentry;e++)
1682  {
1683  FiniteElement* el_pt=Bin_object_coord_pairs[b][e].first;
1685  unsigned dim=this->nlagrangian();
1686  Vector<double> zeta(dim);
1688  {
1689  el_pt->interpolated_x(s,zeta);
1690  }
1691  else
1692  {
1693  el_pt->interpolated_zeta(s,zeta);
1694  }
1695  for (unsigned i=0;i<dim;i++)
1696  {
1697  outfile << zeta[i] << " ";
1698  }
1699  outfile << ix << " "
1700  << iy << " "
1701  << iz << " "
1702  << b << " "
1703  << std::endl;
1704  }
1705  b++;
1706  }
1707  }
1708  }
1709  }
1710 
1711 
1712 //========================================================================
1713 /// Compute the minimum distance of any vertex in the specified bin
1714 /// from the specified Lagrangian coordinate zeta.
1715 //========================================================================
1716  double MeshAsGeomObject::min_distance(const unsigned& i_bin,
1717  const Vector<double>& zeta)
1718  {
1719  // Spatial dimension of bin
1720  const unsigned n_lagrangian = this->nlagrangian();
1721 
1722  // Get bin vertices
1723  Vector<Vector<double> > bin_vertex;
1724  get_bin_vertices(i_bin, bin_vertex);
1725 
1726  double min_dist=DBL_MAX;
1727  unsigned nvertex=bin_vertex.size();
1728  for (unsigned v=0;v<nvertex;v++)
1729  {
1730  double dist=0.0;
1731  for (unsigned i=0;i<n_lagrangian;i++)
1732  {
1733  dist+=pow(bin_vertex[v][i]-zeta[i],2);
1734  }
1735  dist=sqrt(dist);
1736  if (dist<min_dist) min_dist=dist;
1737  }
1738  return min_dist;
1739  }
1740 
1741 
1742 
1743 //========================================================================
1744 /// Output bin vertices (allowing display of bins as zones).
1745 /// Final argument specifies the coordinates of a point
1746 /// and output includes the minimum distance of any of
1747 /// the bin vertices to this point.
1748 //========================================================================
1749  void MeshAsGeomObject::output_bin_vertices(std::ostream& outfile,
1750  const Vector<double>& zeta)
1751  {
1752 
1753  // Spatial dimension of bin
1754  const unsigned n_lagrangian = this->nlagrangian();
1755 
1756  unsigned nbin=Nbin_x;
1757  if (n_lagrangian>1) nbin*=Nbin_y;
1758  if (n_lagrangian>2) nbin*=Nbin_z;
1759 
1760  for (unsigned i_bin=0;i_bin<nbin;i_bin++)
1761  {
1762  // Get bin vertices
1763  Vector<Vector<double> > bin_vertex;
1764  get_bin_vertices(i_bin, bin_vertex);
1765  switch(n_lagrangian)
1766  {
1767  case 1:
1768  outfile << "ZONE I=2\n";
1769  break;
1770 
1771  case 2:
1772  outfile << "ZONE I=2, J=2\n";
1773  break;
1774 
1775  case 3:
1776  outfile << "ZONE I=2, J=2, K=2\n";
1777  break;
1778  }
1779 
1780  unsigned nvertex=bin_vertex.size();
1781  for (unsigned i=0;i<nvertex;i++)
1782  {
1783  for (unsigned j=0;j<n_lagrangian;j++)
1784  {
1785  outfile
1786  << bin_vertex[i][j] << " ";
1787  }
1788  outfile << i_bin << " "
1789  << min_distance(i_bin,zeta) << "\n";
1790  }
1791  }
1792  }
1793 
1794 
1795 //========================================================================
1796 /// Output bin vertices of specified bin (allowing display of bins as zones).
1797 /// Final argument specifies the coordinates of a point
1798 /// and output includes the minimum distance of any of
1799 /// the bin vertices to this point.
1800 //========================================================================
1801  void MeshAsGeomObject::output_bin_vertices(std::ostream& outfile,
1802  const unsigned& i_bin,
1803  const Vector<double>& zeta)
1804  {
1805 
1806  // Spatial dimension of bin
1807  const unsigned n_lagrangian = this->nlagrangian();
1808 
1809  // Get bin vertices
1810  Vector<Vector<double> > bin_vertex;
1811  get_bin_vertices(i_bin, bin_vertex);
1812  switch(n_lagrangian)
1813  {
1814  case 1:
1815  outfile << "ZONE I=2\n";
1816  break;
1817 
1818  case 2:
1819  outfile << "ZONE I=2, J=2\n";
1820  break;
1821 
1822  case 3:
1823  outfile << "ZONE I=2, J=2, K=2\n";
1824  break;
1825  }
1826 
1827  unsigned nvertex=bin_vertex.size();
1828  for (unsigned i=0;i<nvertex;i++)
1829  {
1830  for (unsigned j=0;j<n_lagrangian;j++)
1831  {
1832  outfile
1833  << bin_vertex[i][j] << " ";
1834  }
1835  outfile << i_bin << " "
1836  << min_distance(i_bin,zeta) << "\n";
1837  }
1838  }
1839 
1840 
1841 
1842 //========================================================================
1843 /// Get vector of vectors containing the coordinates of the
1844 /// vertices of the i_bin-th bin: bin_vertex[j][i] contains the
1845 /// i-th coordinate of the j-th vertex.
1846 //========================================================================
1847  void MeshAsGeomObject::get_bin_vertices(const unsigned& i_bin,
1848  Vector<Vector<double> >& bin_vertex)
1849  {
1850 
1851  // Spatial dimension of bin
1852  const unsigned n_lagrangian = this->nlagrangian();
1853 
1854  // How many vertices do we have?
1855  unsigned n_vertices=1;
1856  for (unsigned i=0;i<n_lagrangian;i++)
1857  {
1858  n_vertices*=2;
1859  }
1860  bin_vertex.resize(n_vertices);
1861 
1862  // Vectors for min [0] and max [1] coordinates of the bin in each
1863  // coordinate direction
1864  Vector<Vector<double> > zeta_vertex_bin(2);
1865  zeta_vertex_bin[0].resize(n_lagrangian);
1866  zeta_vertex_bin[1].resize(n_lagrangian);
1867 
1868  unsigned Nbin[3]={Nbin_x,Nbin_y,Nbin_z};
1869  Vector<double> dzeta;
1870  unsigned count=0;
1871  Vector<unsigned> i_1d;
1872 
1873  // Deal with different spatial dimensions separately
1874  switch (n_lagrangian)
1875  {
1876 
1877  case 1:
1878 
1879  // Assign vertex coordinates of the bin directly
1880  dzeta.resize(1);
1881  dzeta[0]=(Max_coords[0]-Min_coords[0])/double(Nbin[0]);
1882  bin_vertex[0].resize(1);
1883  bin_vertex[0][0]=Min_coords[0]+double(i_bin)*dzeta[0];
1884  bin_vertex[1].resize(1);
1885  bin_vertex[1][0]=Min_coords[0]+double(i_bin+1)*dzeta[0];
1886 
1887  break;
1888 
1889  case 2:
1890 
1891  dzeta.resize(2);
1892  dzeta[0]=(Max_coords[0]-Min_coords[0])/double(Nbin[0]);
1893  dzeta[1]=(Max_coords[1]-Min_coords[1])/double(Nbin[1]);
1894 
1895  i_1d.resize(2);
1896  i_1d[0]=i_bin%Nbin[0];
1897  i_1d[1]=(i_bin-i_1d[0])/Nbin[0];
1898 
1899  // Max/min coordinates of the bin
1900  for (unsigned i=0;i<n_lagrangian;i++)
1901  {
1902  zeta_vertex_bin[0][i]=Min_coords[i]+double(i_1d[i])*dzeta[i];
1903  zeta_vertex_bin[1][i]=Min_coords[i]+double(i_1d[i]+1)*dzeta[i];
1904  }
1905 
1906  // Build 4 vertex coordinates
1907  count=0;
1908  for (unsigned i_min_max=0;i_min_max<2;i_min_max++)
1909  {
1910  for (unsigned j_min_max=0;j_min_max<2;j_min_max++)
1911  {
1912  bin_vertex[count].resize(2);
1913  bin_vertex[count][0]=zeta_vertex_bin[i_min_max][0];
1914  bin_vertex[count][1]=zeta_vertex_bin[j_min_max][1];
1915  count++;
1916  }
1917  }
1918 
1919  break;
1920 
1921  case 3:
1922 
1923  dzeta.resize(3);
1924  dzeta[0]=(Max_coords[0]-Min_coords[0])/double(Nbin[0]);
1925  dzeta[1]=(Max_coords[1]-Min_coords[1])/double(Nbin[1]);
1926  dzeta[2]=(Max_coords[2]-Min_coords[2])/double(Nbin[2]);
1927 
1928  i_1d.resize(3);
1929  i_1d[0]=i_bin%Nbin[0];
1930  i_1d[1]=((i_bin-i_1d[0])/Nbin[0])%Nbin[1];
1931  i_1d[2]=(i_bin-(i_1d[1]*Nbin[0]+i_1d[0]))/(Nbin[0]*Nbin[1]);
1932 
1933  // Max/min coordinates of the bin
1934  for (unsigned i=0;i<n_lagrangian;i++)
1935  {
1936  zeta_vertex_bin[0][i]=Min_coords[i]+double(i_1d[i])*dzeta[i];
1937  zeta_vertex_bin[1][i]=Min_coords[i]+double(i_1d[i]+1)*dzeta[i];
1938  }
1939 
1940  // Build 8 vertex coordinates
1941  count=0;
1942  for (unsigned i_min_max=0;i_min_max<2;i_min_max++)
1943  {
1944  for (unsigned j_min_max=0;j_min_max<2;j_min_max++)
1945  {
1946  for (unsigned k_min_max=0;k_min_max<2;k_min_max++)
1947  {
1948  bin_vertex[count].resize(3);
1949  bin_vertex[count][0]=zeta_vertex_bin[i_min_max][0];
1950  bin_vertex[count][1]=zeta_vertex_bin[j_min_max][1];
1951  bin_vertex[count][2]=zeta_vertex_bin[k_min_max][2];
1952  count++;
1953  }
1954  }
1955  }
1956 
1957  break;
1958 
1959  default:
1960 
1961  std::ostringstream error_message;
1962  error_message
1963  << "Can't deal with bins in dimension " << n_lagrangian << "\n";
1964  throw OomphLibError(
1965  error_message.str(),
1966  OOMPH_CURRENT_FUNCTION,
1967  OOMPH_EXCEPTION_LOCATION);
1968  }
1969 
1970  }
1971 
1972 //========================================================================
1973 ///Calculate the bin numbers of all the neighbours to "bin" given the level
1974 //========================================================================
1976  const unsigned& bin, const unsigned& level,
1977  Vector<unsigned>& neighbour_bin,
1978  const bool &only_use_filled_bins)
1979  {
1980  const unsigned n_lagrangian = this->nlagrangian();
1981  // This will be different depending on the number of Lagrangian
1982  // coordinates
1983  if (n_lagrangian==1)
1984  {
1985  // Reserve memory for the container where we return the indices
1986  // of the neighbouring bins (2 bins max, left and right)
1987  neighbour_bin.reserve(2);
1988 
1989  // Single "loop" in one direction - always a vector of max size 2
1990  unsigned nbr_bin_left=bin-level;
1991  if (nbr_bin_left<Nbin_x)
1992  {
1993  unsigned nbr_bin=nbr_bin_left;
1994  if (only_use_filled_bins)
1995  {
1996  if (Bin_object_coord_pairs.has_entry(nbr_bin))
1997  {
1998  neighbour_bin.push_back(nbr_bin);
1999  }
2000  }
2001  else
2002  {
2003  neighbour_bin.push_back(nbr_bin);
2004  }
2005  }
2006  unsigned nbr_bin_right=bin+level;
2007  if ((nbr_bin_right<Nbin_x) &&
2008  (nbr_bin_right!=nbr_bin_left))
2009  {
2010  unsigned nbr_bin=nbr_bin_right;
2011  if (only_use_filled_bins)
2012  {
2013  if (Bin_object_coord_pairs.has_entry(nbr_bin))
2014  {
2015  neighbour_bin.push_back(nbr_bin);
2016  }
2017  }
2018  else
2019  {
2020  neighbour_bin.push_back(nbr_bin);
2021  }
2022  }
2023  }
2024  else if (n_lagrangian==2)
2025  {
2026  // Reserve memory for the container where we return the indices
2027  // of the neighbouring bins
2028  const unsigned n_max_neighbour_bins = 8*level;
2029  neighbour_bin.reserve(n_max_neighbour_bins);
2030 
2031  const unsigned n_total_bin=Nbin_x*Nbin_y;
2032 
2033  // Which row of the bin structure is the current bin on?
2034  // This is just given by the integer answer of dividing bin
2035  // by Nbin_x (the number of bins in a single row)
2036  // e.g. in a 6x6 grid, bins 6 through 11 would all have bin_row=1
2037  const unsigned bin_row=bin/Nbin_x;
2038 
2039  // The neighbour_bin vector contains all bin numbers at the
2040  // specified "distance" (level) away from the current bin
2041 
2042  // Row/column length
2043  const unsigned n_length=(level*2)+1;
2044 
2045  {
2046  // Visit all the bins at the specified distance (level) away
2047  // from the current bin. In order to obtain the same order in
2048  // the visited bins as the previous algorithm we visit all the
2049  // bins at the specified distance (level) as follows:
2050 
2051  // Suppose we want the bins at distance (level=2) of the
2052  // specified bin, then we visit them as indicated below
2053 
2054  // 01 02 03 04 05 // First row
2055  // 06 07
2056  // 08 B 09
2057  // 10 11
2058  // 12 13 14 15 16 // Last row
2059  // ^--------------- First column
2060  // ^-- Last column
2061 
2062  // ----------------------------------------------------------------
2063  // Visit all the bins in the first row at the specified
2064  // distance (level) away from the current bin
2065 
2066  // ------------------ FIRST ROW ------------------------
2067  // Pre-compute the distance in the j-direction
2068  const unsigned j_precomputed = level*Nbin_x;
2069  // Pre-compute the row where the bin should lie on
2070  const unsigned j_initial_row=bin_row-level;
2071 
2072  // Loop over the columns (of the first row)
2073  for (unsigned i=0;i<n_length;i++)
2074  {
2075  // --------- First row ------------------------------------------
2076  const unsigned initial_neighbour_bin=bin-(level-i)-j_precomputed;
2077  // This number might fall on the wrong row of the bin
2078  // structure; this needs to be tested? Not sure why, but leave
2079  // the test!
2080 
2081  // Which row is this number on? (see above)
2082  const unsigned initial_neighbour_bin_row=initial_neighbour_bin/Nbin_x;
2083  // These numbers for the rows must match; if it is then add
2084  // initial_neighbour_bin to the neighbour scheme (The bin
2085  // number must also be greater than zero and less than the
2086  // total number of bins)
2087  if ((j_initial_row==initial_neighbour_bin_row) &&
2088  (initial_neighbour_bin<n_total_bin))
2089  {
2090  if (only_use_filled_bins)
2091  {
2092  if (Bin_object_coord_pairs.has_entry(initial_neighbour_bin))
2093  {
2094  neighbour_bin.push_back(initial_neighbour_bin);
2095  }
2096  }
2097  else
2098  {
2099  neighbour_bin.push_back(initial_neighbour_bin);
2100  }
2101  }
2102 
2103  } // for (unsigned i=0;i<n_length;i++)
2104 
2105  // Then visit all the bins in the first and last column at the
2106  // specified distance (level) away from the current bin
2107 
2108  // ------------------ FIRST AND LAST COLUMNS ---------------------
2109  // Loop over the rows (of the first and last column)
2110  for (unsigned j=1;j<n_length-1;j++)
2111  {
2112  // --------- First column ---------------------------------------
2113  const unsigned initial_neighbour_bin=bin-(level)-((level-j)*Nbin_x);
2114  // This number might fall on the wrong row of the bin
2115  // structure; this needs to be tested? Not sure why, but leave
2116  // the test!
2117 
2118  // Which row is this number on? (see above)
2119  const unsigned initial_neighbour_bin_row=initial_neighbour_bin/Nbin_x;
2120 
2121  // Which row should it be on?
2122  const unsigned initial_row=bin_row-(level-j);
2123 
2124  // These numbers for the rows must match; if it is then add
2125  // initial_neighbour_bin to the neighbour scheme (The bin
2126  // number must also be greater than zero and less than the
2127  // total number of bins)
2128  if ((initial_row==initial_neighbour_bin_row) &&
2129  (initial_neighbour_bin<n_total_bin))
2130  {
2131  if (only_use_filled_bins)
2132  {
2133  if (Bin_object_coord_pairs.has_entry(initial_neighbour_bin))
2134  {
2135  neighbour_bin.push_back(initial_neighbour_bin);
2136  }
2137  }
2138  else
2139  {
2140  neighbour_bin.push_back(initial_neighbour_bin);
2141  }
2142  }
2143 
2144  // --------- Last column -----------------------------------------
2145  const unsigned final_neighbour_bin=bin+(level)-((level-j)*Nbin_x);
2146  // This number might fall on the wrong row of the bin
2147  // structure; this needs to be tested? Not sure why, but leave
2148  // the test!
2149 
2150  // Which row is this number on? (see above)
2151  const unsigned final_neighbour_bin_row=final_neighbour_bin/Nbin_x;
2152 
2153  // Which row should it be on?
2154  const unsigned final_row=bin_row-(level-j);
2155 
2156  // These numbers for the rows must match; if it is then add
2157  // initial_neighbour_bin to the neighbour scheme (The bin
2158  // number must also be greater than zero and less than the
2159  // total number of bins)
2160  if ((final_row==final_neighbour_bin_row) &&
2161  (final_neighbour_bin<n_total_bin))
2162  {
2163  if (only_use_filled_bins)
2164  {
2165  if (Bin_object_coord_pairs.has_entry(final_neighbour_bin))
2166  {
2167  neighbour_bin.push_back(final_neighbour_bin);
2168  }
2169  }
2170  else
2171  {
2172  neighbour_bin.push_back(final_neighbour_bin);
2173  }
2174  }
2175 
2176  } // for (unsigned j=1;j<n_length-1;j++)
2177 
2178  // ------------------ LAST ROW ------------------------
2179  // Pre-compute the row where the bin should lie on
2180  const unsigned j_final_row=bin_row+level;
2181 
2182  // Loop over the columns (of the last row)
2183  for (unsigned i=0;i<n_length;i++)
2184  {
2185  // --------- Last row ------------------------------------------
2186  const unsigned final_neighbour_bin=bin-(level-i)+j_precomputed;
2187  // This number might fall on the wrong row of the bin
2188  // structure; this needs to be tested? Not sure why, but leave
2189  // the test!
2190 
2191  // Which row is this number on? (see above)
2192  const unsigned final_neighbour_bin_row=final_neighbour_bin/Nbin_x;
2193  // These numbers for the rows must match; if it is then add
2194  // initial_neighbour_bin to the neighbour scheme (The bin
2195  // number must also be greater than zero and less than the
2196  // total number of bins)
2197  if ((j_final_row==final_neighbour_bin_row) &&
2198  (final_neighbour_bin<n_total_bin))
2199  {
2200  if (only_use_filled_bins)
2201  {
2202  if (Bin_object_coord_pairs.has_entry(final_neighbour_bin))
2203  {
2204  neighbour_bin.push_back(final_neighbour_bin);
2205  }
2206  }
2207  else
2208  {
2209  neighbour_bin.push_back(final_neighbour_bin);
2210  }
2211  }
2212 
2213  } // for (unsigned i=0;i<n_length;i++)
2214 
2215  }
2216 
2217  }
2218  else if (n_lagrangian==3)
2219  {
2220  // Reserve memory for the container where we return the indices
2221  // of the neighbouring bins
2222  const unsigned n_max_neighbour_bins =
2223  8*level*(3+2*(level-1))+2*(2*(level-1)+1)*(2*(level-1)+1);
2224  neighbour_bin.reserve(n_max_neighbour_bins);
2225 
2226  unsigned n_total_bin=Nbin_x*Nbin_y*Nbin_z;
2227 
2228  // Which layer of the bin structure is the current bin on?
2229  // This is just given by the integer answer of dividing bin
2230  // by Nbin_x*Nbin_y (the number of bins in a single layer
2231  // e.g. in a 6x6x6 grid, bins 72 through 107 would all have bin_layer=2
2232  unsigned bin_layer=bin/(Nbin_x*Nbin_y);
2233 
2234  // Which row in this layer is the bin number on?
2235  unsigned bin_row=(bin/Nbin_x)-(bin_layer*Nbin_y);
2236 
2237  // The neighbour_bin vector contains all bin numbers at the
2238  // specified "distance" (level) away from the current bin
2239 
2240  // Row/column/layer length
2241  unsigned n_length=(level*2)+1;
2242 
2243  // Loop over the layers
2244  for (unsigned k=0;k<n_length;k++)
2245  {
2246  // Loop over the rows
2247  for (unsigned j=0;j<n_length;j++)
2248  {
2249  // Loop over the columns
2250  for (unsigned i=0;i<n_length;i++)
2251  {
2252  // Only do this for the end points of every row/layer/column
2253  if ((k==0) || (k==n_length-1) || (j==0) ||
2254  (j==n_length-1) || (i==0) || (i==n_length-1))
2255  {
2256  unsigned nbr_bin=bin-level+i-((level-j)*Nbin_x)-
2257  ((level-k)*Nbin_x*Nbin_y);
2258  // This number might fall on the wrong
2259  // row or layer of the bin structure; this needs to be tested
2260 
2261  // Which layer is this number on?
2262  unsigned nbr_bin_layer=nbr_bin/(Nbin_x*Nbin_y);
2263 
2264  // Which row is this number on? (see above)
2265  unsigned nbr_bin_row=(nbr_bin/Nbin_x)-(nbr_bin_layer*Nbin_y);
2266 
2267  // Which layer and row should it be on, given level?
2268  unsigned layer=bin_layer-level+k;
2269  unsigned row=bin_row-level+j;
2270 
2271  // These layers and rows must match up:
2272  // if so then add nbr_bin to the neighbour schemes
2273  // (The bin number must also be greater than zero
2274  // and less than the total number of bins)
2275  if ((row==nbr_bin_row) && (layer==nbr_bin_layer)
2276  && (nbr_bin<n_total_bin))
2277  {
2278  if (only_use_filled_bins)
2279  {
2280  if (Bin_object_coord_pairs.has_entry(nbr_bin))
2281  {
2282  neighbour_bin.push_back(nbr_bin);
2283  }
2284  }
2285  else
2286  {
2287  neighbour_bin.push_back(nbr_bin);
2288  }
2289  }
2290  }
2291 
2292  }
2293  }
2294  }
2295 
2296  }
2297  }
2298 
2299 }
unsigned nbin_z()
Number of bins in z direction.
unsigned total_nbin() const
Total number of bins (empty or not)
unsigned nbin_y()
Number of bins in y direction.
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
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 long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:590
cstr elem_len * i
Definition: cfortran.h:607
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
unsigned nbin_x()
Number of bins in x direction.
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1305
Vector< double > Max_coords
Storage for max coordinates in the mesh.
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
unsigned Nz_bin
Number of bins in the third dimension in binning method in setup_multi_domain_interaction().
A general Finite Element class.
Definition: elements.h:1271
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.
OomphInfo oomph_info
static bool mpi_has_been_initialised()
return true if MPI has been initialised
double Percentage_offset
Percentage offset to add to each extreme of the bin structure. Default value of 5 (note that this is ...
e
Definition: cfortran.h:575
static bool Use_eulerian_coordinates_during_setup
Globally set-able (sorry!) flag to indicate that MeshAsGeomObject is to use Eulerian coordinates when...
void output_bins(std::ofstream &outfile)
Output bins.
double X_max
Maximum coordinate in first dimension.
bool Have_used_eulerian_coordinates_during_setup
Flag to indicate that MeshAsGeomObject has used Eulerian coordinates when setting up bin...
double Y_max
Maximum coordinate in second dimension.
double Y_min
Minimum coordinate in second dimension.
double Z_max
Maximum coordinate in third dimension.
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
OomphCommunicator * communicator_pt() const
Definition: mesh.h:1317
void spiraling_locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s, const bool &called_within_spiral)
Find the sub geometric object and local coordinate therein that corresponds to the intrinsic coordina...
OomphCommunicator * Communicator_pt
Communicator.
unsigned nlagrangian() const
Access function to # of Lagrangian coordinates.
Definition: geom_objects.h:182
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
unsigned Nbin_y
Number of bins in y direction.
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 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 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_min_and_max_coordinates(Mesh *const &mesh_pt)
Get the min and max coordinates for the mesh, in each dimension.
Vector< double > Min_coords
Storage for min coordinates in the mesh.
static char t char * s
Definition: cfortran.h:572
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...
void my_spiraling_locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s, const bool &called_within_spiral)
Version of spiraling locate zeta used for the projection during the unstructured mesh adaptation...
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...
void create_bins_of_objects(const unsigned long &n_mesh_element)
static unsigned long Total_nbin_cells_counter
Counter for overall number of bins allocated – used to issue warning if this exceeds a threshhold...
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...
double timer()
returns the time in seconds after some point in past
static bool Suppress_warning_about_large_total_number_of_bins
Boolean to supppress warnings about large number of bins.
void get_neighbouring_bins_helper(const unsigned &bin_number, const unsigned &level, Vector< unsigned > &neighbour_bin, const bool &only_use_filled_bins=false)
Calculate the bin numbers of all the neighbours to "bin" given the level.
bool Sort_bin_entries
Bool to decide if to sort entries in bin during locate_zeta operation (default: false) ...
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:456
bool Allow_use_of_halo_elements_as_external_elements_for_projection
Indicate whether we are allowed to use halo elements as external elements for projection, possibly only required in parallel unstructured mesh generation during the projection stage. Default set to true.
unsigned Nx_bin
Number of bins in the first dimension in binning method in setup_multi_domain_interaction().
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...
static bool Suppress_warning_about_small_number_of_bins
Boolean to supppress warnings about small number of bins.
bool Setup_multi_domain_for_projection
Set up multi-domain for projection.
Vector< FiniteElement * > Sub_geom_object_pt
Internal storage for the elements that constitute the object.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
unsigned Nbin_z
Number of bins in z direction.
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
Definition: elements.cc:4480
void locate_zeta(const Vector< double > &zeta, GeomObject *&geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
For a given value of zeta, the "global" intrinsic coordinate of a mesh of FiniteElements represente...
Definition: elements.cc:4533
unsigned Nsample_points
(Measure of) the number of sampling points within the elements when populating the bin ...
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
void sort_the_bin(const Vector< double > &zeta, Vector< std::pair< FiniteElement *, Vector< double > > > &sample_point_pairs)
Sort the sampling points in the specified bin by distance from sampling point.
unsigned Nbin_x
Number of bins in x direction.
double Z_min
Minimum coordinate in third dimension.
bool Suppress_synchronisation_of_bins
In parallel computation, suppress synchronisation of bins Default is false. If set to true...
void flush_bins_of_objects()
Flush the storage for the binning method (and decrement counter for total number of bins in active us...
bool Doc_stats
Boolean to indicate whether to output basic info during setup_multi_domain_interaction() routines...
unsigned Ny_bin
Number of bins in the second dimension in binning method in setup_multi_domain_interaction().
void construct_it(Mesh *const &mesh_pt, const bool &compute_extreme_bin_coords)
Helper function for constructor: pass the pointer to the mesh, and boolean to specify whether to calc...
bool is_halo() const
Is this element a halo?
Definition: elements.h:1141
A general mesh class.
Definition: mesh.h:74
double X_min
Minimum and maximum coordinates for each dimension of the external mesh used to "create" the bins in ...
void output_bin_vertices(std::ostream &outfile, const Vector< double > &zeta)
Output bin vertices (allowing display of bins as zones). Final argument specifies the coordinates of ...
SparseVector< Vector< std::pair< FiniteElement *, Vector< double > > > > Bin_object_coord_pairs
Storage for paired objects and coords in each bin.
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:57