multi_domain.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 //Non-templated multi-domain functions which act on more than one mesh
31 //and set up the storage and interaction between the two
32 
33 //oomph-lib header
34 #include "multi_domain.h"
35 #include "multi_domain.template.cc"
36 #include "mesh.h"
37 #include "algebraic_elements.h"
39 #include "Qelements.h"
40 
41 namespace oomph
42 {
43 
44 //======================================================================
45 // Namespace for "global" multi-domain functions
46 //======================================================================
47 namespace Multi_domain_functions
48  {
49 
50  /// \short Output file to document the boundary coordinate
51  /// along the mesh boundary of the bulk mesh during call to
52  /// setup_bulk_elements_adjacent_to_face_mesh(...)
54 
55  // Workspace for locate zeta methods
56  //----------------------------------
57 
58  /// \short Boolean to indicate that failure in setup multi domain
59  /// functions is acceptable; defaults to false. If set to true
60  /// external element pointers are set to null for those elements
61  /// for which external elements couldn't be located.
63 
64  /// \short Dimension of zeta tuples (set by get_dim_helper) -- needed
65  /// because we store the scalar coordinates in flat-packed form.
66  unsigned Dim;
67 
68  /// \short Lookup scheme for whether a local element's integration point
69  /// has had an external element assigned to it -- essentially boolean.
70  /// External_element_located[e][ipt] = {0,1} if external element
71  /// for ipt-th integration in local element e {has not, has} been found.
72  /// Used locally to ensure that we're not searching for the same
73  /// elements over and over again when we go around the spirals.
75 
76  /// \short Vector of flat-packed zeta coordinates for which the external
77  /// element could not be found during current local search. These
78  /// will be sent to the next processor in the ring-like parallel search.
79  /// The zeta coordinates come in groups of Dim (scalar) coordinates.
81 
82  /// \short Vector of flat-packed zeta coordinates for which the external
83  /// element could not be found on another processor and for which
84  /// we're currently searching here. Whatever can't be found here,
85  /// gets written into Flat_packed_zetas_not_found_locally and then
86  /// passed on to the next processor during the ring-like parallel search.
87  /// The zeta coordinates come in groups of Dim (scalar) coordinates.
89 
90  /// \short Proc_id_plus_one_of_external_element[i] contains the
91  /// processor id (plus one) of the processor
92  /// on which the i-th zeta coordinate tuple received from elsewhere
93  /// (in the order in which these are stored in
94  /// Received_flat_packed_zetas_to_be_found) was located; it's zero if
95  /// it wasn't found during the current stage of the ring-like parallel
96  /// search.
98 
99  /// \short Vector to indicate (to another processor) whether a
100  /// located element (that will have to represented as an external
101  /// halo element on that processor) should be newly created on that
102  /// processor (2), already exists on that processor (1), or
103  /// is not on the current processor either (0).
105 
106  /// \short Vector of flat-packed local coordinates for zeta tuples
107  /// that have been located
109 
110  /// \short Vector of flat-packed doubles to be communicated with
111  /// other processors
113 
114  /// \short Counter used when processing vector of flat-packed
115  /// doubles -- this is really "private" data, declared here
116  /// to avoid having to pass it (and the associated array)
117  /// between the various helper functions
119 
120  /// \short Vector of flat-packed unsigneds to be communicated with
121  /// other processors -- this is really "private" data, declared here
122  /// to avoid having to pass the array between the various helper
123  /// functions
125 
126 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
127 
128  // Temporary vector of strings to enable full annotation of multi domain
129  // comms (but keep alive because it would be such a bloody pain to
130  // rewrite it if things ever go wrong again...)
132 
133 
134 #endif
135 
136  /// \short Counter used when processing vector of flat-packed
137  /// unsigneds -- this is really "private" data, declared here
138  /// to avoid having to pass it (and the associated array)
139  /// between the various helper functions
141 
142 
143  // Default parameters for the binning method
144  //------------------------------------------
145 
146  /// \short Bool to tell the MeshAsGeomObject whether to calculate
147  /// the extreme coordinates of the bin structure
149 
150  /// \short Number of bins in the first dimension in binning method in
151  /// setup_multi_domain_interaction().
152  unsigned Nx_bin=100;
153 
154  /// \short Number of bins in the second dimension in binning method in
155  /// setup_multi_domain_interaction().
156  unsigned Ny_bin=100;
157 
158  /// \short Number of bins in the third dimension in binning method in
159  /// setup_multi_domain_interaction().
160  unsigned Nz_bin=100;
161 
162  /// Number of spirals to be searched in one go
163  unsigned N_spiral_chunk=1;
164 
165  /// \short (Measure of) the number of sampling points within the elements
166  /// when populating the bin
167  unsigned Nsample_points=5;
168 
169  /// \short Minimum and maximum coordinates for
170  /// each dimension of the external mesh used to "create" the bins in
171  /// setup_multi_domain_interaction(). These can be set by user if they
172  /// want to (otherwise the MeshAsGeomObject calculates these values based
173  /// upon the mesh itself; see MeshAsGeomObject::get_max_and_min_coords(...))
174  /// They default to "incorrect" values initially
175 
176  /// \short Minimum coordinate in first dimension
177  double X_min=1.0;
178 
179  /// \short Maximum coordinate in first dimension
180  double X_max=-1.0;
181 
182  /// \short Minimum coordinate in second dimension
183  double Y_min=1.0;
184 
185  /// \short Maximum coordinate in second dimension
186  double Y_max=-1.0;
187 
188  /// \short Minimum coordinate in third dimension
189  double Z_min=1.0;
190 
191  /// \short Maximum coordinate in third dimension
192  double Z_max=-1.0;
193 
194  /// \short Percentage offset to add to each extreme of the bin structure.
195  /// Default value of 5 (note that this is in percent!)
196  double Percentage_offset=5.0;
197 
198 
199  // Other parameters
200  //-----------------
201 
202  /// \short Boolean to indicate when to use the bulk element as the
203  /// external element. Defaults to false, you must have set up FaceElements
204  /// properly first in order for it to work
206 
207  /// \short Boolean to indicate if we're allowed to use halo elements
208  /// as external elements. Can drastically reduce the number of
209  /// external halo elements -- currently not aware of any problems
210  /// therefore set to true by default but retention
211  /// of this flag allows easy return to previous implementation.
213 
214  /// \short The timing for sorting the elements in the bins
216 
217  /// Set up multi-domain for projection
219 
220  /// \short Indicate whether we are allowed to use halo elements as
221  /// external elements for projection, possibly only required in
222  /// parallel unstructured mesh generation during the projection
223  /// stage. Default set to true
225 
226  /// \short Boolean to indicate whether to doc timings or not.
227  bool Doc_timings=false;
228 
229  /// \short Boolean to indicate whether to output basic info during
230  /// setup_multi_domain_interaction() routines
231  bool Doc_stats=false;
232 
233  /// \short Boolean to indicate whether to output further info during
234  /// setup_multi_domain_interaction() routines
235  bool Doc_full_stats=false;
236 
237 #ifdef OOMPH_HAS_MPI
238 
239 
240  // Functions for location method in multi-domain problems
241 
242 
243 
244 //========================================================================
245 /// Send the zeta coordinates from the current process to
246 /// the next process; receive from the previous process
247 //========================================================================
249  {
250  // MPI info
251  MPI_Status status;
252  MPI_Request request;
253 
254  // Storage for number of processors, current process and communicator
255  int n_proc=problem_pt->communicator_pt()->nproc();
256  int my_rank=problem_pt->communicator_pt()->my_rank();
257  OomphCommunicator* comm_pt=problem_pt->communicator_pt();
258 
259  // Work out processors to send and receive from
260  int send_to_proc=my_rank+1;
261  int recv_from_proc=my_rank-1;
262  if (send_to_proc==n_proc) { send_to_proc=0; }
263  if (recv_from_proc<0) { recv_from_proc=n_proc-1; }
264 
265  // Send the number of flat-packed zetas that we couldn't find
266  // locally to the next processor
267  int n_missing_local_zetas=Flat_packed_zetas_not_found_locally.size();
268  MPI_Isend(&n_missing_local_zetas,1,MPI_INT,send_to_proc,4,
269  comm_pt->mpi_comm(),&request);
270 
271  // Receive the number of flat-packed zetas that couldn't be found
272  // on the "previous" processor
273  int count_zetas=0;
274  MPI_Recv(&count_zetas,1,MPI_INT,recv_from_proc,
275  4,comm_pt->mpi_comm(),&status);
276 
277  MPI_Wait(&request,MPI_STATUS_IGNORE);
278 
279  // Send the vector of flat-packed zetas that we couldn't find
280  // locally to the next processor
281  if (n_missing_local_zetas!=0)
282  {
284  n_missing_local_zetas,MPI_DOUBLE,
285  send_to_proc,5,comm_pt->mpi_comm(),&request);
286  }
287 
288  // Receive the vector of flat-packed zetas that couldn't be found
289  // on the "previous" processor
290  if (count_zetas!=0)
291  {
292  Received_flat_packed_zetas_to_be_found.resize(count_zetas);
294  count_zetas,MPI_DOUBLE,recv_from_proc,5,
295  comm_pt->mpi_comm(),&status);
296  MPI_Wait(&request,MPI_STATUS_IGNORE);
297  }
298  else
299  {
301  }
302 
303  // Now we should have the Zeta arrays set up correctly
304  // for the next round of locations
305  }
306 
307 //========start of send_and_receive_located_info==========================
308 /// Send location information from current process; Received location
309 /// information from (current process + iproc) modulo (nproc)
310 //========================================================================
312  (int& iproc, Mesh* const &external_mesh_pt, Problem* problem_pt)
313  {
314  // Set MPI info
315  MPI_Status status;
316  MPI_Request request;
317 
318  // Storage for number of processors, current process and communicator
319  OomphCommunicator* comm_pt=problem_pt->communicator_pt();
320  int n_proc=comm_pt->nproc();
321  int my_rank=comm_pt->my_rank();
322 
323  // Prepare vectors to receive information
324  Vector<double> received_double_values;
325  Vector<unsigned> received_unsigned_values;
326  Vector<double> received_located_coord;
327  Vector<int> received_proc_id_plus_one_of_external_element;
328  Vector<unsigned> received_located_element_status;
329 
330  // Communicate the located information back to the original process
331  int orig_send_proc=my_rank-iproc;
332  if (my_rank<iproc) { orig_send_proc=n_proc+orig_send_proc; }
333  int orig_recv_proc=my_rank+iproc;
334  if ((my_rank+iproc)>=n_proc) { orig_recv_proc=orig_recv_proc-n_proc; }
335 
336  // Send the double values associated with external halos
337  //------------------------------------------------------
338  unsigned send_count_double_values=Flat_packed_doubles.size();
339  MPI_Isend(&send_count_double_values,1,MPI_UNSIGNED,
340  orig_send_proc,1,comm_pt->mpi_comm(),&request);
341  int receive_count_double_values=0;
342  MPI_Recv(&receive_count_double_values,1,MPI_INT,
343  orig_recv_proc,1,comm_pt->mpi_comm(),&status);
344  MPI_Wait(&request,MPI_STATUS_IGNORE);
345 
346  if (send_count_double_values!=0)
347  {
348  MPI_Isend(&Flat_packed_doubles[0],send_count_double_values,MPI_DOUBLE,
349  orig_send_proc,2,comm_pt->mpi_comm(),&request);
350  }
351  if (receive_count_double_values!=0)
352  {
353  received_double_values.resize(receive_count_double_values);
354  MPI_Recv(&received_double_values[0],receive_count_double_values,
355  MPI_DOUBLE,orig_recv_proc,2,comm_pt->mpi_comm(),&status);
356  }
357  if (send_count_double_values!=0)
358  {
359  MPI_Wait(&request,MPI_STATUS_IGNORE);
360  }
361 
362  // Now send unsigned values associated with external halos
363  //---------------------------------------------------------
364  unsigned send_count_unsigned_values=Flat_packed_unsigneds.size();
365  MPI_Isend(&send_count_unsigned_values,1,MPI_UNSIGNED,
366  orig_send_proc,14,comm_pt->mpi_comm(),&request);
367 
368  int receive_count_unsigned_values=0;
369  MPI_Recv(&receive_count_unsigned_values,1,MPI_INT,orig_recv_proc,14,
370  comm_pt->mpi_comm(),&status);
371 
372  MPI_Wait(&request,MPI_STATUS_IGNORE);
373 
374  if (send_count_unsigned_values!=0)
375  {
376  MPI_Isend(&Flat_packed_unsigneds[0],send_count_unsigned_values,MPI_UNSIGNED,
377  orig_send_proc,15,comm_pt->mpi_comm(),&request);
378 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
379  for (unsigned i=0;i<send_count_unsigned_values;i++)
380  {
381  oomph_info << "Sent:" << i << " to orig_proc:" << orig_send_proc
382  << " " << Flat_packed_unsigneds_string[i]
383  << ": " << Flat_packed_unsigneds[i] << std::endl;
384  }
385 #endif
386  }
387  if (receive_count_unsigned_values!=0)
388  {
389  received_unsigned_values.resize(receive_count_unsigned_values);
390  MPI_Recv(&received_unsigned_values[0],receive_count_unsigned_values,
391  MPI_UNSIGNED,orig_recv_proc,15,comm_pt->mpi_comm(),&status);
392  }
393 
394  if (send_count_unsigned_values!=0)
395  {
396  MPI_Wait(&request,MPI_STATUS_IGNORE);
397  }
398 
399  // Send and receive the Located_element_status
400  //--------------------------------------------
401  int send_count=Received_flat_packed_zetas_to_be_found.size()/Dim;
402  MPI_Isend(&send_count,1,MPI_INT,orig_send_proc,
403  20,comm_pt->mpi_comm(),&request);
404  int receive_count=0;
405  MPI_Recv(&receive_count,1,MPI_INT,orig_recv_proc,
406  20,comm_pt->mpi_comm(),&status);
407  MPI_Wait(&request,MPI_STATUS_IGNORE);
408 
409  if (send_count!=0)
410  {
411  MPI_Isend(&Located_element_status[0],send_count,MPI_UNSIGNED,
412  orig_send_proc,3,comm_pt->mpi_comm(),&request);
413  }
414 
415  if (receive_count!=0)
416  {
417  received_located_element_status.resize(receive_count);
418  MPI_Recv(&received_located_element_status[0],receive_count,
419  MPI_UNSIGNED,orig_recv_proc,3,
420  comm_pt->mpi_comm(),&status);
421  }
422  if (send_count!=0)
423  {
424  MPI_Wait(&request,MPI_STATUS_IGNORE);
425  }
426 
427  // Send and receive Proc_id_plus_one_of_external_element array
428  //------------------------------------------------------------
429  if (send_count!=0)
430  {
431  MPI_Isend(&Proc_id_plus_one_of_external_element[0],send_count,MPI_INT,
432  orig_send_proc,13,comm_pt->mpi_comm(),&request);
433  }
434  if (receive_count!=0)
435  {
436  received_proc_id_plus_one_of_external_element.resize(receive_count);
437  MPI_Recv(&received_proc_id_plus_one_of_external_element[0],
438  receive_count,MPI_INT,orig_recv_proc,13,
439  comm_pt->mpi_comm(),&status);
440  }
441  if (send_count!=0)
442  {
443  MPI_Wait(&request,MPI_STATUS_IGNORE);
444  }
445 
446 
447  // And finally the Flat_packed_located_coordinates array
448  //------------------------------------------------------
449  unsigned send_count_located_coord=Flat_packed_located_coordinates.size();
450  MPI_Isend(&send_count_located_coord,1,MPI_UNSIGNED,
451  orig_send_proc,4,comm_pt->mpi_comm(),&request);
452  unsigned receive_count_located_coord=0;
453  MPI_Recv(&receive_count_located_coord,1,MPI_UNSIGNED,orig_recv_proc,4,
454  comm_pt->mpi_comm(),&status);
455  MPI_Wait(&request,MPI_STATUS_IGNORE);
456 
457  if (send_count_located_coord!=0)
458  {
459  MPI_Isend(&Flat_packed_located_coordinates[0],
460  send_count_located_coord,MPI_DOUBLE,
461  orig_send_proc,5,comm_pt->mpi_comm(),&request);
462  }
463  if (receive_count_located_coord!=0)
464  {
465  received_located_coord.resize(receive_count_located_coord);
466  MPI_Recv(&received_located_coord[0],receive_count_located_coord,MPI_DOUBLE,
467  orig_recv_proc,5,comm_pt->mpi_comm(),&status);
468  }
469  if (send_count_located_coord!=0)
470  {
471  MPI_Wait(&request,MPI_STATUS_IGNORE);
472  }
473 
474  // Copy across into original containers -- these can now
475  //------------------------------------------------------
476  // be processed by create_external_halo_elements() to generate
477  //------------------------------------------------------------
478  // external halo elements
479  //------------------------
480  Flat_packed_doubles.resize(receive_count_double_values);
481  for (int ii=0;ii<receive_count_double_values;ii++)
482  {
483  Flat_packed_doubles[ii]=received_double_values[ii];
484  }
485  Flat_packed_unsigneds.resize(receive_count_unsigned_values);
486  for (int ii=0;ii<receive_count_unsigned_values;ii++)
487  {
488  Flat_packed_unsigneds[ii]=received_unsigned_values[ii];
489  }
490  Proc_id_plus_one_of_external_element.resize(receive_count);
491  Located_element_status.resize(receive_count);
492  for (int ii=0;ii<receive_count;ii++)
493  {
495  received_proc_id_plus_one_of_external_element[ii];
496  Located_element_status[ii]=received_located_element_status[ii];
497  }
498  Flat_packed_located_coordinates.resize(receive_count_located_coord);
499  for (int ii=0;ii<int(receive_count_located_coord);ii++)
500  {
501  Flat_packed_located_coordinates[ii]=received_located_coord[ii];
502  }
503 
504  }
505 
506 
507 //========start of add_external_haloed_node_to_storage====================
508 /// Helper function to add external haloed nodes, including any masters
509 //========================================================================
510  void add_external_haloed_node_to_storage(int& iproc, Node* nod_pt,
511  Problem* problem_pt,
512  Mesh* const &external_mesh_pt,
513  int& n_cont_inter_values)
514  {
515  // Add the node if required
516  add_external_haloed_node_helper(iproc,nod_pt,problem_pt,external_mesh_pt,
517  n_cont_inter_values);
518 
519  // Recursively add any master nodes (and their master nodes etc)
521  problem_pt,
522  external_mesh_pt,
523  n_cont_inter_values);
524  }
525 
526 
527 
528  //========================================================================
529  ///Recursively add any master nodes (and their master nodes etc) of
530  /// external nodes
531  //========================================================================
533  int& iproc, Node* nod_pt,
534  Problem* problem_pt,
535  Mesh* const &external_mesh_pt,
536  int& n_cont_inter_values)
537  {
538 
539  // Loop over continuously interpolated values and add masters
540  for (int i_cont=-1;i_cont<n_cont_inter_values;i_cont++)
541  {
542  if (nod_pt->is_hanging(i_cont))
543  {
544  // Indicate that this node is a hanging node so the other
545  // process knows to create HangInfo and masters, etc.
546  Flat_packed_unsigneds.push_back(1);
547 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
548  Flat_packed_unsigneds_string.push_back("Is hanging");
549 #endif
550  // If this is a hanging node then add all its masters as
551  // external halo nodes if they have not yet been added
552  HangInfo* hang_pt=nod_pt->hanging_pt(i_cont);
553  // Loop over masters
554  unsigned n_master=hang_pt->nmaster();
555 
556  // Indicate number of master nodes to add on other process
557  Flat_packed_unsigneds.push_back(n_master);
558 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
559  Flat_packed_unsigneds_string.push_back("nmaster");
560 #endif
561  for (unsigned m=0;m<n_master;m++)
562  {
563  Node* master_nod_pt=hang_pt->master_node_pt(m);
564 
565  // Call the helper function for master nodes
566  add_external_haloed_master_node_helper(iproc,master_nod_pt,
567  problem_pt,
568  external_mesh_pt,
569  n_cont_inter_values);
570 
571  // Indicate the weight of this master
572  Flat_packed_doubles.push_back(hang_pt->master_weight(m));
573 
574  // Recursively add any master nodes (and their master nodes etc)
576  problem_pt,
577  external_mesh_pt,
578  n_cont_inter_values);
579  }
580  }
581  else
582  {
583  // Indicate that it's not a hanging node in this variable
584  Flat_packed_unsigneds.push_back(0);
585 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
586  Flat_packed_unsigneds_string.push_back("Not hanging");
587 #endif
588  }
589  } // end loop over continously interpolated values
590 
591  }
592 
593 //==========start of add_external_haloed_node_helper======================
594 /// Helper to add external haloed node that is not a master
595 //========================================================================
596  void add_external_haloed_node_helper(int& iproc, Node* nod_pt,
597  Problem* problem_pt,
598  Mesh* const &external_mesh_pt,
599  int& n_cont_inter_values)
600  {
601  // Attempt to add this node as an external haloed node
602  unsigned n_ext_haloed_nod=external_mesh_pt->nexternal_haloed_node(iproc);
603  unsigned external_haloed_node_index=
604  external_mesh_pt->add_external_haloed_node_pt(iproc,nod_pt);
605 
606  // If it was added then the new index should match the size of the storage
607  if (external_haloed_node_index==n_ext_haloed_nod)
608  {
609  Flat_packed_unsigneds.push_back(1);
610 
611 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
612  std::stringstream junk;
613  junk << "Node needs to be constructed [size="
614  << Flat_packed_unsigneds.size() << "]; last entry: "
616  Flat_packed_unsigneds_string.push_back(junk.str());
617 #endif
618 
619  // This helper function gets all the required information for the
620  // specified node and stores it into MPI-sendable information
621  // so that a halo copy can be made on the receiving process
623  problem_pt,external_mesh_pt,
624  n_cont_inter_values);
625  }
626  else // It was already added
627  {
628  Flat_packed_unsigneds.push_back(0);
629 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
630  std::stringstream junk;
631  junk << "Node was already added [size="
632  << Flat_packed_unsigneds.size() << "]; last entry: "
634 
635  Flat_packed_unsigneds_string.push_back(junk.str());
636 #endif
637 
638  // This node is already an external haloed node, so tell
639  // the other process its index in the equivalent external halo storage
640  Flat_packed_unsigneds.push_back(external_haloed_node_index);
641 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
642  Flat_packed_unsigneds_string.push_back("external haloed node index");
643 #endif
644  }
645 
646 
647 
648  }
649 
650 
651 //==========start of add_external_haloed_master_node_helper===============
652 /// Helper function to add external haloed node that is a master
653 //========================================================================
654  void add_external_haloed_master_node_helper(int& iproc, Node* master_nod_pt,
655  Problem* problem_pt,
656  Mesh* const &external_mesh_pt,
657  int& n_cont_inter_values)
658  {
659  // Attempt to add node as an external haloed node
660  unsigned n_ext_haloed_nod=external_mesh_pt->nexternal_haloed_node(iproc);
661  unsigned external_haloed_node_index;
662  external_haloed_node_index=
663  external_mesh_pt->add_external_haloed_node_pt(iproc,master_nod_pt);
664 
665  // If it was added the returned index is the same as current storage size
666  if (external_haloed_node_index==n_ext_haloed_nod)
667  {
668  // Indicate that this node needs to be constructed on
669  // the other process
670  Flat_packed_unsigneds.push_back(1);
671 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
672  Flat_packed_unsigneds_string.push_back("Node needs to be constructed[2]");
673 #endif
674 
675  // This gets all the required information for the specified
676  // master node and stores it into MPI-sendable information
677  // so that a halo copy can be made on the receiving process
679  problem_pt,
680  external_mesh_pt,
681  n_cont_inter_values);
682  }
683  else // It was already added
684  {
685  Flat_packed_unsigneds.push_back(0);
686 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
687  Flat_packed_unsigneds_string.push_back("Node was already added[2]");
688 #endif
689 
690  // This node is already an external haloed node, so tell
691  // the other process its index in the equivalent external halo storage
692  Flat_packed_unsigneds.push_back(external_haloed_node_index);
693 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
694  Flat_packed_unsigneds_string.push_back("external haloed node index[2]");
695 #endif
696  }
697  }
698 
699 
700 
701 
702 //========start of get_required_nodal_information_helper==================
703 /// Helper function to get the required nodal information from an
704 /// external haloed node so that a fully-functional external halo
705 /// node (and therefore element) can be created on the receiving process
706 //========================================================================
707  void get_required_nodal_information_helper(int& iproc, Node* nod_pt,
708  Problem* problem_pt,
709  Mesh* const &external_mesh_pt,
710  int& n_cont_inter_values)
711  {
712  // Tell the halo copy of this node how many values there are
713  // [NB this may be different for nodes within the same element, e.g.
714  // when using Lagrange multipliers]
715  unsigned n_val=nod_pt->nvalue();
716  Flat_packed_unsigneds.push_back(n_val);
717 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
718  Flat_packed_unsigneds_string.push_back("Number of values");
719 #endif
720 
721  unsigned n_dim=nod_pt->ndim();
722  TimeStepper* time_stepper_pt=nod_pt->time_stepper_pt();
723 
724  // Find the timestepper in the list of problem timesteppers
725  bool found_timestepper=false;
726  unsigned time_stepper_index;
727  unsigned n_time_steppers=problem_pt->ntime_stepper();
728  for (unsigned i=0;i<n_time_steppers;i++)
729  {
730  if (time_stepper_pt==problem_pt->time_stepper_pt(i))
731  {
732  // Indicate the timestepper's index
733  found_timestepper=true;
734  time_stepper_index=i;
735  break;
736  }
737  }
738 
739  if (found_timestepper)
740  {
741  Flat_packed_unsigneds.push_back(1);
742 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
743  Flat_packed_unsigneds_string.push_back("Found timestepper");
744 #endif
745  Flat_packed_unsigneds.push_back(time_stepper_index);
746 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
747  Flat_packed_unsigneds_string.push_back("Timestepper index");
748 #endif
749  }
750  else
751  {
752  Flat_packed_unsigneds.push_back(0);
753 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
754  Flat_packed_unsigneds_string.push_back("Not found timestepper");
755 #endif
756  }
757 
758  // Default number of previous values to 1
759  unsigned n_prev=1;
760  if (time_stepper_pt!=0)
761  {
762  // Add number of history values to n_prev
763  n_prev=time_stepper_pt->ntstorage();
764  }
765 
766  // Is the node on any boundaries?
767  if (nod_pt->is_on_boundary())
768  {
769  Flat_packed_unsigneds.push_back(1);
770 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
771  Flat_packed_unsigneds_string.push_back("Node is on boundary");
772 #endif
773 
774  // Loop over the boundaries of the external mesh
775  Vector<unsigned> boundaries;
776  unsigned n_bnd=external_mesh_pt->nboundary();
777  for (unsigned i_bnd=0;i_bnd<n_bnd;i_bnd++)
778  {
779  // Which boundaries (could be more than one) is it on?
780  if (nod_pt->is_on_boundary(i_bnd))
781  {
782  boundaries.push_back(i_bnd);
783  }
784  }
785  unsigned nb=boundaries.size();
786  Flat_packed_unsigneds.push_back(nb);
787 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
788  std::stringstream junk;
789  junk << "Node is on "<< nb << " boundaries";
790  Flat_packed_unsigneds_string.push_back(junk.str());
791 #endif
792  for (unsigned i=0;i<nb;i++)
793  {
794  Flat_packed_unsigneds.push_back(boundaries[i]);
795 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
796  std::stringstream junk;
797  junk << "Node is on boundary " << boundaries[i] << " of " << n_bnd;
798  Flat_packed_unsigneds_string.push_back(junk.str());
799 #endif
800  }
801 
802  // Get pointer to the map of indices associated with
803  // additional values created by face elements
804  BoundaryNodeBase* bnod_pt=
805  dynamic_cast<BoundaryNodeBase*>(nod_pt);
806 #ifdef PARANOID
807  if (bnod_pt==0)
808  {
809  throw OomphLibError(
810  "Failed to cast new node to boundary node\n",
811  OOMPH_CURRENT_FUNCTION,
812  OOMPH_EXCEPTION_LOCATION);
813  }
814 #endif
815  std::map<unsigned, unsigned>* map_pt=
817 
818  // No additional values created
819  if (map_pt==0)
820  {
821  Flat_packed_unsigneds.push_back(0);
822 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
823  std::stringstream junk;
824  Flat_packed_unsigneds_string.push_back("No additional values were created by face element");
825 #endif
826  }
827  // Created additional values
828  else
829  {
830  // How many?
831  Flat_packed_unsigneds.push_back(map_pt->size());
832 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
833  std::stringstream junk;
834  junk << "Map size " << map_pt->size() << n_bnd;
835  Flat_packed_unsigneds_string.push_back(junk.str());
836 #endif
837  // Loop over entries in map and add to send data
838  for (std::map<unsigned, unsigned>::iterator p=
839  map_pt->begin();
840  p!=map_pt->end();p++)
841  {
842  Flat_packed_unsigneds.push_back((*p).first);
843 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
844  std::stringstream junk;
845  Flat_packed_unsigneds_string.push_back("Key of map entry");
846 #endif
847  Flat_packed_unsigneds.push_back((*p).second);
848 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
849  Flat_packed_unsigneds_string.push_back("Value of map entry");
850 #endif
851  }
852  }
853  }
854  else
855  {
856  // Not on any boundary
857  Flat_packed_unsigneds.push_back(0);
858 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
859  Flat_packed_unsigneds_string.push_back("Node is not on any boundary");
860 #endif
861  }
862 
863  // Is the Node algebraic? If so, send its ref values and
864  // an indication of its geometric objects if they are stored
865  // in the algebraic mesh
866  AlgebraicNode* alg_nod_pt=dynamic_cast<AlgebraicNode*>(nod_pt);
867  if (alg_nod_pt!=0)
868  {
869  // The external mesh should be algebraic
870  AlgebraicMesh* alg_mesh_pt=dynamic_cast<AlgebraicMesh*>
871  (external_mesh_pt);
872 
873  // Get default node update function ID
874  unsigned update_id=alg_nod_pt->node_update_fct_id();
875  Flat_packed_unsigneds.push_back(update_id);
876 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
877  Flat_packed_unsigneds_string.push_back("Alg Node update id");
878 #endif
879 
880  // Get reference values at default...
881  unsigned n_ref_val=alg_nod_pt->nref_value();
882  Flat_packed_unsigneds.push_back(n_ref_val);
883 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
884  Flat_packed_unsigneds_string.push_back("Alg Node n ref values");
885 #endif
886  for (unsigned i_ref_val=0;i_ref_val<n_ref_val;i_ref_val++)
887  {
888  Flat_packed_doubles.push_back(alg_nod_pt->ref_value(i_ref_val));
889  }
890 
891  // Access geometric objects at default...
892  unsigned n_geom_obj=alg_nod_pt->ngeom_object();
893  Flat_packed_unsigneds.push_back(n_geom_obj);
894 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
895  Flat_packed_unsigneds_string.push_back("Alg Node n geom objects");
896 #endif
897  for (unsigned i_geom=0;i_geom<n_geom_obj;i_geom++)
898  {
899  GeomObject* geom_obj_pt=alg_nod_pt->geom_object_pt(i_geom);
900 
901  // Check this against the stored geometric objects in mesh
902  unsigned n_geom_list=alg_mesh_pt->ngeom_object_list_pt();
903 
904  // Default found index to zero
905  unsigned found_geom_object=0;
906  for (unsigned i_list=0;i_list<n_geom_list;i_list++)
907  {
908  if (geom_obj_pt==alg_mesh_pt->geom_object_list_pt(i_list))
909  {
910  found_geom_object=i_list;
911  }
912  }
913  Flat_packed_unsigneds.push_back(found_geom_object);
914 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
915  Flat_packed_unsigneds_string.push_back("Found geom object");
916 #endif
917  }
918  }
919 
920  // If it is a MacroElementNodeUpdateNode, everything has been
921  // dealt with by the new element already
922 
923  // Is it a SolidNode?
924  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
925  if (solid_nod_pt!=0)
926  {
927  unsigned n_solid_val=solid_nod_pt->variable_position_pt()->nvalue();
928  for (unsigned i_val=0;i_val<n_solid_val;i_val++)
929  {
930  for (unsigned t=0;t<n_prev;t++)
931  {
932  Flat_packed_doubles.push_back(solid_nod_pt->variable_position_pt()->
933  value(t,i_val));
934  }
935  }
936  }
937 
938  // Finally copy info required for all node types
939  for (unsigned i_val=0;i_val<n_val;i_val++)
940  {
941  for (unsigned t=0;t<n_prev;t++)
942  {
943  Flat_packed_doubles.push_back(nod_pt->value(t,i_val));
944  }
945  }
946 
947  // Now do positions
948  for (unsigned idim=0;idim<n_dim;idim++)
949  {
950  for (unsigned t=0;t<n_prev;t++)
951  {
952  Flat_packed_doubles.push_back(nod_pt->x(t,idim));
953  }
954  }
955  }
956 
957 //=========start of get_required_master_nodal_information_helper==========
958 /// Helper function to get the required master nodal information from an
959 /// external haloed master node so that a fully-functional external halo
960 /// master node (and possible element) can be created on the receiving process
961 //========================================================================
963  (int& iproc, Node* master_nod_pt, Problem* problem_pt,
964  Mesh* const &external_mesh_pt, int& n_cont_inter_values)
965  {
966  // Need to send over dimension, position type and number of values
967  Flat_packed_unsigneds.push_back(master_nod_pt->ndim());
968 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
969  Flat_packed_unsigneds_string.push_back("Master node ndim");
970 #endif
971  Flat_packed_unsigneds.push_back(master_nod_pt->nposition_type());
972 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
973  Flat_packed_unsigneds_string.push_back("Master node npos_type");
974 #endif
975  Flat_packed_unsigneds.push_back(master_nod_pt->nvalue());
976 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
977  Flat_packed_unsigneds_string.push_back("Master node nvalue");
978 #endif
979 
980  // If it's a solid node, also need to send lagrangian dim and type
981  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(master_nod_pt);
982  if (solid_nod_pt!=0)
983  {
984  Flat_packed_unsigneds.push_back(solid_nod_pt->nlagrangian());
985 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
986  Flat_packed_unsigneds_string.push_back("Master solid node nlagr");
987 #endif
988  Flat_packed_unsigneds.push_back(solid_nod_pt->nlagrangian_type());
989 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
990  Flat_packed_unsigneds_string.push_back("Master solid node nlagr_type");
991 #endif
992  }
993 
994  unsigned n_dim=master_nod_pt->ndim();
995  TimeStepper* time_stepper_pt=master_nod_pt->time_stepper_pt();
996 
997  // Find the timestepper in the list of problem timesteppers
998  bool found_timestepper=false;
999  unsigned time_stepper_index;
1000  unsigned n_time_steppers=problem_pt->ntime_stepper();
1001  for (unsigned i=0;i<n_time_steppers;i++)
1002  {
1003  if (time_stepper_pt==problem_pt->time_stepper_pt(i))
1004  {
1005  // Indicate the timestepper's index
1006  // add 1 to the index so that 0 indicates no timestepper?
1007  found_timestepper=true;
1008  time_stepper_index=i;
1009  break;
1010  }
1011  }
1012 
1013  if (found_timestepper)
1014  {
1015  Flat_packed_unsigneds.push_back(1);
1016 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1017  Flat_packed_unsigneds_string.push_back("Master node Found timestepper");
1018 #endif
1019  Flat_packed_unsigneds.push_back(time_stepper_index);
1020 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1021  Flat_packed_unsigneds_string.push_back("Master node Timestepper index");
1022 #endif
1023  }
1024  else
1025  {
1026  Flat_packed_unsigneds.push_back(0);
1027 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1028  Flat_packed_unsigneds_string.push_back("Master node Not found timestepper");
1029 #endif
1030  }
1031 
1032  // Default number of previous values to 1
1033  unsigned n_prev=1;
1034  if (time_stepper_pt!=0)
1035  {
1036  // Add number of history values to n_prev
1037  n_prev=time_stepper_pt->ntstorage();
1038  }
1039 
1040  // Is the node on any boundaries?
1041  if (master_nod_pt->is_on_boundary())
1042  {
1043  Flat_packed_unsigneds.push_back(1);
1044 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1045  Flat_packed_unsigneds_string.push_back("Master node is on boundary");
1046 #endif
1047  // Loop over the boundaries of the external mesh
1048  Vector<unsigned> boundaries;
1049  unsigned n_bnd=external_mesh_pt->nboundary();
1050  for (unsigned i_bnd=0;i_bnd<n_bnd;i_bnd++)
1051  {
1052  // Which boundaries (could be more than one) is it on?
1053  if (master_nod_pt->is_on_boundary(i_bnd))
1054  {
1055  boundaries.push_back(i_bnd);
1056  }
1057  }
1058  unsigned nb=boundaries.size();
1059  Flat_packed_unsigneds.push_back(nb);
1060 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1061  std::stringstream junk;
1062  junk << "Master node is on "<< nb << " boundaries";
1063  Flat_packed_unsigneds_string.push_back(junk.str());
1064 #endif
1065  for (unsigned i=0;i<nb;i++)
1066  {
1067  Flat_packed_unsigneds.push_back(boundaries[i]);
1068 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1069  std::stringstream junk;
1070  junk << "Master noode is on boundary "
1071  << boundaries[i] << " of " << n_bnd;
1072  Flat_packed_unsigneds_string.push_back(junk.str());
1073 #endif
1074  }
1075 
1076  // Get pointer to the map of indices associated with
1077  // additional values created by face elements
1078  BoundaryNodeBase* bnod_pt=
1079  dynamic_cast<BoundaryNodeBase*>(master_nod_pt);
1080 #ifdef PARANOID
1081  if (bnod_pt==0)
1082  {
1083  throw OomphLibError(
1084  "Failed to cast new node to boundary node\n",
1085  OOMPH_CURRENT_FUNCTION,
1086  OOMPH_EXCEPTION_LOCATION);
1087  }
1088 #endif
1089  std::map<unsigned, unsigned>* map_pt=
1091 
1092  // No additional values created
1093  if (map_pt==0)
1094  {
1095  Flat_packed_unsigneds.push_back(0);
1096 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1097  std::stringstream junk;
1098  Flat_packed_unsigneds_string.push_back("No additional values were created by face element for this master node");
1099 #endif
1100  }
1101  // Created additional values
1102  else
1103  {
1104  // How many?
1105  Flat_packed_unsigneds.push_back(map_pt->size());
1106 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1107  std::stringstream junk;
1108  junk << "Map size for master node " << map_pt->size() << n_bnd;
1109  Flat_packed_unsigneds_string.push_back(junk.str());
1110 #endif
1111  // Loop over entries in map and add to send data
1112  for (std::map<unsigned, unsigned>::iterator p=
1113  map_pt->begin();
1114  p!=map_pt->end();p++)
1115  {
1116  Flat_packed_unsigneds.push_back((*p).first);
1117 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1118  std::stringstream junk;
1119  Flat_packed_unsigneds_string.push_back(
1120  "Key of map entry for master node");
1121 #endif
1122  Flat_packed_unsigneds.push_back((*p).second);
1123 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1124  Flat_packed_unsigneds_string.push_back(
1125  "Value of map entry for master node");
1126 #endif
1127  }
1128  }
1129  }
1130  else
1131  {
1132  // Not on any boundary
1133  Flat_packed_unsigneds.push_back(0);
1134 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1135  Flat_packed_unsigneds_string.push_back("Master node is not on any boundary");
1136 #endif
1137  }
1138 
1139  // Is the Node algebraic? If so, send its ref values and
1140  // an indication of its geometric objects if they are stored
1141  // in the algebraic mesh
1142  AlgebraicNode* alg_nod_pt=dynamic_cast<AlgebraicNode*>(master_nod_pt);
1143  if (alg_nod_pt!=0)
1144  {
1145  // The external mesh should be algebraic
1146  AlgebraicMesh* alg_mesh_pt=dynamic_cast<AlgebraicMesh*>
1147  (external_mesh_pt);
1148 
1149  // Get default node update function ID
1150  unsigned update_id=alg_nod_pt->node_update_fct_id();
1151  Flat_packed_unsigneds.push_back(update_id);
1152 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1153  Flat_packed_unsigneds_string.push_back("Master Alg Node update id");
1154 #endif
1155 
1156  // Get reference values at default...
1157  unsigned n_ref_val=alg_nod_pt->nref_value();
1158  Flat_packed_unsigneds.push_back(n_ref_val);
1159 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1160  Flat_packed_unsigneds_string.push_back("Master Alg Node n ref values");
1161 #endif
1162  for (unsigned i_ref_val=0;i_ref_val<n_ref_val;i_ref_val++)
1163  {
1164  Flat_packed_doubles.push_back(alg_nod_pt->ref_value(i_ref_val));
1165  }
1166 
1167  // Access geometric objects at default...
1168  unsigned n_geom_obj=alg_nod_pt->ngeom_object();
1169  Flat_packed_unsigneds.push_back(n_geom_obj);
1170 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1171  Flat_packed_unsigneds_string.push_back("Master Alg Node n geom objects");
1172 #endif
1173  for (unsigned i_geom=0;i_geom<n_geom_obj;i_geom++)
1174  {
1175  GeomObject* geom_obj_pt=alg_nod_pt->geom_object_pt(i_geom);
1176  // Check this against the stored geometric objects in mesh
1177  unsigned n_geom_list=alg_mesh_pt->ngeom_object_list_pt();
1178  // Default found index to zero
1179  unsigned found_geom_object=0;
1180  for (unsigned i_list=0;i_list<n_geom_list;i_list++)
1181  {
1182  if (geom_obj_pt==alg_mesh_pt->geom_object_list_pt(i_list))
1183  {
1184  found_geom_object=i_list;
1185  }
1186  }
1187  Flat_packed_unsigneds.push_back(found_geom_object);
1188 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1189  Flat_packed_unsigneds_string.push_back("Master node Found geom object");
1190 #endif
1191  }
1192  } // end AlgebraicNode check
1193 
1194  // Is it a MacroElementNodeUpdateNode?
1195  MacroElementNodeUpdateNode* macro_nod_pt=
1196  dynamic_cast<MacroElementNodeUpdateNode*>(master_nod_pt);
1197  if (macro_nod_pt!=0)
1198  {
1199  // Loop over current external haloed elements - has the element which
1200  // controls the node update for this node been added yet?
1201  GeneralisedElement* macro_node_update_el_pt=
1202  macro_nod_pt->node_update_element_pt();
1203 
1204  unsigned n_ext_haloed_el=external_mesh_pt->
1205  nexternal_haloed_element(iproc);
1206  unsigned external_haloed_el_index;
1207  external_haloed_el_index=external_mesh_pt->
1208  add_external_haloed_element_pt(iproc,macro_node_update_el_pt);
1209 
1210  // If it wasn't already added, we need to create a halo copy
1211  if (external_haloed_el_index==n_ext_haloed_el)
1212  {
1213  Flat_packed_unsigneds.push_back(1);
1214 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1215  Flat_packed_unsigneds_string.push_back("Master Node needs to be constructed");
1216 #endif
1217  //Cast to a finite elemnet
1218  FiniteElement* macro_node_update_finite_el_pt
1219  = dynamic_cast<FiniteElement*>(macro_node_update_el_pt);
1220 
1221  // We're using macro elements to update...
1222  MacroElementNodeUpdateMesh* macro_mesh_pt=
1223  dynamic_cast<MacroElementNodeUpdateMesh*>(external_mesh_pt);
1224  if (macro_mesh_pt!=0)
1225  {
1226  Flat_packed_unsigneds.push_back(1);
1227 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1228  Flat_packed_unsigneds_string.push_back("Mesh is macro element mesh");
1229 #endif
1230  // Need to send the macro element number in the mesh across
1231  MacroElement* macro_el_pt= macro_node_update_finite_el_pt
1232  ->macro_elem_pt();
1233  unsigned macro_el_num=macro_el_pt->macro_element_number();
1234  Flat_packed_unsigneds.push_back(macro_el_num);
1235 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1236  Flat_packed_unsigneds_string.push_back("Number of macro element");
1237 #endif
1238  // Also need to send
1239  // the lower left and upper right coordinates of the macro element
1240  QElementBase* q_el_pt=dynamic_cast<QElementBase*>
1241  (macro_node_update_el_pt);
1242  if (q_el_pt!=0)
1243  {
1244  // The macro element needs to be set first before
1245  // its lower left and upper right coordinates can be accessed
1246  // Now send the lower left and upper right coordinates
1247  unsigned el_dim=q_el_pt->dim();
1248  for (unsigned i_dim=0;i_dim<el_dim;i_dim++)
1249  {
1250  Flat_packed_doubles.push_back(q_el_pt->s_macro_ll(i_dim));
1251  Flat_packed_doubles.push_back(q_el_pt->s_macro_ur(i_dim));
1252  }
1253  }
1254  else // Throw an error
1255  {
1256  std::ostringstream error_stream;
1257  error_stream << "You are using a MacroElement node update\n"
1258  << "in a case with non-QElements. This has not\n"
1259  << "yet been implemented.\n";
1260  throw OomphLibError
1261  (error_stream.str(),
1262  OOMPH_CURRENT_FUNCTION,
1263  OOMPH_EXCEPTION_LOCATION);
1264  }
1265 
1266  }
1267  else // Not using macro elements for node update... umm, we're
1268  // already inside a loop over macro elements, so this
1269  // should never get here... an error should be thrown I suppose
1270  {
1271  Flat_packed_unsigneds.push_back(0);
1272 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1273  Flat_packed_unsigneds_string.push_back("Mesh is not a macro element mesh");
1274 #endif
1275  }
1276 
1277  // This element needs to be fully functioning on the other
1278  // process, so send all the information required to create it
1279  unsigned n_node=macro_node_update_finite_el_pt->nnode();
1280  for (unsigned j=0;j<n_node;j++)
1281  {
1282  Node* new_nod_pt=macro_node_update_finite_el_pt->node_pt(j);
1283  add_external_haloed_node_to_storage(iproc,new_nod_pt,
1284  problem_pt,
1285  external_mesh_pt,
1286  n_cont_inter_values);
1287  }
1288  }
1289  else // The external haloed element already exists
1290  {
1291  Flat_packed_unsigneds.push_back(0);
1292 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1293  Flat_packed_unsigneds_string.push_back("External haloed element already exists");
1294 #endif
1295  Flat_packed_unsigneds.push_back(external_haloed_el_index);
1296 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1297  Flat_packed_unsigneds_string.push_back("Index of existing external haloed element");
1298 #endif
1299  }
1300 
1301  } // end of MacroElementNodeUpdateNode check
1302 
1303  // Is it a SolidNode?
1304  if (solid_nod_pt!=0)
1305  {
1306  unsigned n_val=solid_nod_pt->variable_position_pt()->nvalue();
1307  for (unsigned i_val=0;i_val<n_val;i_val++)
1308  {
1309  for (unsigned t=0;t<n_prev;t++)
1310  {
1311  Flat_packed_doubles.push_back(solid_nod_pt->variable_position_pt()->
1312  value(t,i_val));
1313  }
1314  }
1315  }
1316 
1317  // Finally copy info required for all node types
1318 
1319  // Halo copy needs to know all the history values
1320  unsigned n_val=master_nod_pt->nvalue();
1321  for (unsigned i_val=0;i_val<n_val;i_val++)
1322  {
1323  for (unsigned t=0;t<n_prev;t++)
1324  {
1325  Flat_packed_doubles.push_back(master_nod_pt->value(t,i_val));
1326  }
1327  }
1328 
1329  // Now do positions
1330  for (unsigned idim=0;idim<n_dim;idim++)
1331  {
1332  for (unsigned t=0;t<n_prev;t++)
1333  {
1334  Flat_packed_doubles.push_back(master_nod_pt->x(t,idim));
1335  }
1336  }
1337 
1338  }
1339 
1340 
1341 
1342 
1343 
1344 
1345 
1346 //=======start of add_external_halo_node_helper===========================
1347 /// Helper functiono to add external halo node that is not a master
1348 //========================================================================
1350  (Node* &new_nod_pt, Mesh* const &external_mesh_pt, unsigned& loc_p,
1351  unsigned& node_index, FiniteElement* const &new_el_pt,
1352  int& n_cont_inter_values,
1353  Problem* problem_pt)
1354  {
1355  // Given the node and the external mesh, and received information
1356  // about them from process loc_p, construct them on the current process
1357 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1359  << " Bool: New node needs to be constructed "
1361  << std::endl;
1362 #endif
1364  {
1365  // Construct a new node based upon sent information
1366  construct_new_external_halo_node_helper(new_nod_pt,loc_p,
1367  node_index,new_el_pt,
1368  external_mesh_pt,problem_pt);
1369  }
1370  else
1371  {
1372 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1374  << " Index of existing external halo node "
1376  << std::endl;
1377 #endif
1378 
1379  // Copy node from received location
1380  new_nod_pt=external_mesh_pt->external_halo_node_pt
1382 
1383  new_el_pt->node_pt(node_index)=new_nod_pt;
1384 
1385 
1386  }
1387  }
1388 
1389 
1390 
1391 
1392 //========start of construct_new_external_halo_node_helper=================
1393 /// Helper function which constructs a new external halo node (on new element)
1394 /// with the required information sent from the haloed process
1395 //========================================================================
1397  (Node* &new_nod_pt, unsigned& loc_p, unsigned& node_index,
1398  FiniteElement* const &new_el_pt,
1399  Mesh* const &external_mesh_pt, Problem* problem_pt)
1400  {
1401  // The first entry indicates the number of values at this new Node
1402  // (which may be different across the same element e.g. Lagrange multipliers)
1403 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1405  << " Number of values of external halo node "
1407  << std::endl;
1408 #endif
1410 
1411  // Null TimeStepper for now
1412  TimeStepper* time_stepper_pt=0;
1413  // Default number of previous values to 1
1414  unsigned n_prev=1;
1415 
1416  // The next entry in Flat_packed_unsigneds indicates
1417  // if a timestepper is required for this halo node
1418 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1420  << " Timestepper req'd for node "
1422  << std::endl;
1423 #endif
1425  {
1426 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1428  << " Index of timestepper "
1430  << std::endl;
1431 #endif
1432  // Index
1433  time_stepper_pt=problem_pt->time_stepper_pt
1435 
1436  // Check whether number of prev values is "sent" across
1437  n_prev=time_stepper_pt->ntstorage();
1438  }
1439 
1440  // If this node was on a boundary then it needs to
1441  // be on the same boundary here
1442 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1444  << " Is node on boundary? "
1446  << std::endl;
1447 #endif
1449  {
1450  // Construct a new boundary node
1451  if (time_stepper_pt!=0)
1452  {
1453  new_nod_pt=new_el_pt->construct_boundary_node
1454  (node_index,time_stepper_pt);
1455  }
1456  else
1457  {
1458  new_nod_pt=new_el_pt->construct_boundary_node(node_index);
1459  }
1460 
1461  // How many boundaries does the node live on?
1462 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1464  << " Number of boundaries the node is on: "
1466  << std::endl;
1467 #endif
1469  for (unsigned i=0;i<nb;i++)
1470  {
1471  // Boundary number
1472 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1474  << " Node is on boundary "
1476  << std::endl;
1477 #endif
1478  unsigned i_bnd=
1480  external_mesh_pt->add_boundary_node(i_bnd,new_nod_pt);
1481  }
1482 
1483  // Do we have additional values created by face elements?
1484 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1486  << " Number of additional values created by face element "
1488  << std::endl;
1489 #endif
1491  if (n_entry>0)
1492  {
1493  // Create storage, if it doesn't already exist, for the map
1494  // that will contain the position of the first entry of
1495  // this face element's additional values,
1496  BoundaryNodeBase* bnew_nod_pt=
1497  dynamic_cast<BoundaryNodeBase*>(new_nod_pt);
1498 #ifdef PARANOID
1499  if (bnew_nod_pt==0)
1500  {
1501  throw OomphLibError(
1502  "Failed to cast new node to boundary node\n",
1503  OOMPH_CURRENT_FUNCTION,
1504  OOMPH_EXCEPTION_LOCATION);
1505  }
1506 #endif
1507  if(bnew_nod_pt->
1508  index_of_first_value_assigned_by_face_element_pt()==0)
1509  {
1510  bnew_nod_pt->
1511  index_of_first_value_assigned_by_face_element_pt()=
1512  new std::map<unsigned, unsigned>;
1513  }
1514 
1515  // Get pointer to the map of indices associated with
1516  // additional values created by face elements
1517  std::map<unsigned, unsigned>* map_pt=
1519 
1520  // Loop over number of entries in map
1521  for (unsigned i=0;i<n_entry;i++)
1522  {
1523  // Read out pairs...
1524 
1525 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1527  << " Key of map entry"
1529  << std::endl;
1530 #endif
1532 
1533 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1535  << " Value of map entry"
1537  << std::endl;
1538 #endif
1540 
1541  // ...and assign
1542  (*map_pt)[first]=second;
1543  }
1544  }
1545  }
1546  else
1547  {
1548  // Construct an ordinary (non-boundary) node
1549  if (time_stepper_pt!=0)
1550  {
1551  new_nod_pt=new_el_pt->construct_node
1552  (node_index,time_stepper_pt);
1553  }
1554  else
1555  {
1556  new_nod_pt=new_el_pt->construct_node(node_index);
1557  }
1558  }
1559 
1560  // Node constructed: add to external halo nodes
1561  external_mesh_pt->add_external_halo_node_pt(loc_p,new_nod_pt);
1562 
1563  // Is the new constructed node Algebraic?
1564  AlgebraicNode* new_alg_nod_pt=dynamic_cast<AlgebraicNode*>
1565  (new_nod_pt);
1566 
1567  // If it is algebraic, its node update functions will
1568  // not yet have been set up properly
1569  if (new_alg_nod_pt!=0)
1570  {
1571  // The AlgebraicMesh is the external mesh
1572  AlgebraicMesh* alg_mesh_pt=dynamic_cast<AlgebraicMesh*>
1573  (external_mesh_pt);
1574 
1575  /// The first entry of All_alg_nodal_info contains
1576  /// the default node update id
1577  /// e.g. for the quarter circle there are
1578  /// "Upper_left_box", "Lower right box" etc...
1579 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1581  << " Alg node update id "
1583  << std::endl;
1584 #endif
1585 
1586  unsigned update_id=Flat_packed_unsigneds
1588 
1589  Vector<double> ref_value;
1590 
1591  // The size of this vector is in the next entry
1592  // of All_alg_nodal_info
1593 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1595  << " Alg node # of ref values "
1597  << std::endl;
1598 #endif
1599  unsigned n_ref_val=Flat_packed_unsigneds
1601 
1602  // The reference values themselves are in
1603  // All_alg_ref_value
1604  ref_value.resize(n_ref_val);
1605  for (unsigned i_ref=0;i_ref<n_ref_val;i_ref++)
1606  {
1607  ref_value[i_ref]=Flat_packed_doubles
1609  }
1610 
1611  Vector<GeomObject*> geom_object_pt;
1612  /// again we need the size of this vector as it varies
1613  /// between meshes; we also need some indication
1614  /// as to which geometric object should be used...
1615 
1616  // The size of this vector is in the next entry
1617  // of All_alg_nodal_info
1618 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1620  << " Alg node # of geom objects "
1622  << std::endl;
1623 #endif
1624  unsigned n_geom_obj=Flat_packed_unsigneds
1626 
1627  // The remaining indices are in the rest of
1628  // All_alg_nodal_info
1629  geom_object_pt.resize(n_geom_obj);
1630  for (unsigned i_geom=0;i_geom<n_geom_obj;i_geom++)
1631  {
1632 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1634  << " Alg node: geom object index "
1636  << std::endl;
1637 #endif
1638  unsigned geom_index=Flat_packed_unsigneds
1640  // This index indicates which of the AlgebraicMesh's
1641  // stored geometric objects should be used
1642  // (0 is a null pointer; everything else should have
1643  // been filled in by the specific Mesh). If it
1644  // hasn't been filled in then the update_node_update
1645  // call should fix it
1646  geom_object_pt[i_geom]=alg_mesh_pt->
1647  geom_object_list_pt(geom_index);
1648  }
1649 
1650  /// For the received update_id, ref_value, geom_object
1651  /// call add_node_update_info
1652  new_alg_nod_pt->add_node_update_info
1653  (update_id,alg_mesh_pt,geom_object_pt,ref_value);
1654 
1655  /// Now call update_node_update
1656  alg_mesh_pt->update_node_update(new_alg_nod_pt);
1657  }
1658 
1659  // Is the node a MacroElementNodeUpdateNode?
1660  MacroElementNodeUpdateNode* macro_nod_pt=
1661  dynamic_cast<MacroElementNodeUpdateNode*>(new_nod_pt);
1662 
1663  if (macro_nod_pt!=0)
1664  {
1665  // Need to call set_node_update_info; this requires
1666  // a Vector<GeomObject*> (taken from the mesh)
1667  Vector<GeomObject*> geom_object_vector_pt;
1668 
1669  // Access the required geom objects from the
1670  // MacroElementNodeUpdateMesh
1671  MacroElementNodeUpdateMesh* macro_mesh_pt=
1672  dynamic_cast<MacroElementNodeUpdateMesh*>
1673  (external_mesh_pt);
1674  geom_object_vector_pt=
1675  macro_mesh_pt->geom_object_vector_pt();
1676 
1677  // Get local coordinate of node in new element
1678  Vector<double> s_in_macro_node_update_element;
1679  new_el_pt->local_coordinate_of_node
1680  (node_index,s_in_macro_node_update_element);
1681 
1682  // Set node update info for this node
1683  macro_nod_pt->set_node_update_info
1684  (new_el_pt,s_in_macro_node_update_element,
1685  geom_object_vector_pt);
1686  }
1687 
1688  // Is the new node a SolidNode?
1689  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(new_nod_pt);
1690  if (solid_nod_pt!=0)
1691  {
1692  unsigned n_solid_val=solid_nod_pt->variable_position_pt()->nvalue();
1693  for (unsigned i_val=0;i_val<n_solid_val;i_val++)
1694  {
1695  for (unsigned t=0;t<n_prev;t++)
1696  {
1697  solid_nod_pt->variable_position_pt()->
1698  set_value(t,i_val,
1700  }
1701  }
1702  }
1703 
1704  // If there are additional values, resize the node
1705  unsigned n_new_val=new_nod_pt->nvalue();
1706  if (n_val>n_new_val)
1707  {
1708  new_nod_pt->resize(n_val);
1709  }
1710 
1711  // Get copied history values
1712  // unsigned n_val=new_nod_pt->nvalue();
1713  for (unsigned i_val=0;i_val<n_val;i_val++)
1714  {
1715  for (unsigned t=0;t<n_prev;t++)
1716  {
1717  new_nod_pt->set_value(t,i_val,Flat_packed_doubles
1719  }
1720  }
1721 
1722  // Get copied history values for positions
1723  unsigned n_dim=new_nod_pt->ndim();
1724  for (unsigned idim=0;idim<n_dim;idim++)
1725  {
1726  for (unsigned t=0;t<n_prev;t++)
1727  {
1728  // Copy to coordinate
1729  new_nod_pt->x(t,idim)=Flat_packed_doubles
1731  }
1732  }
1733  }
1734 
1735 
1736  //=====================================================================
1737  /// Locate zeta for current set of missing coordinates
1738  //=====================================================================
1740  (int& iproc, Mesh* const &external_mesh_pt, Problem* problem_pt,
1741  MeshAsGeomObject* &mesh_geom_obj_pt)
1742  {
1743  // Storage for number of processors, current process and communicator
1744  OomphCommunicator* comm_pt=problem_pt->communicator_pt();
1745  int n_proc=comm_pt->nproc();
1746  int my_rank=comm_pt->my_rank();
1747 
1748  // Clear vectors containing data to be sent
1749  Flat_packed_doubles.resize(0);
1750 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1751  Flat_packed_unsigneds_string.resize(0);
1752 #endif
1753  Flat_packed_unsigneds.resize(0);
1755 
1756  // Flush storage for zetas not found locally (when
1757  // processing the zeta coordinates received from "previous"
1758  // processor)
1760 
1761  // Number of zeta tuples to be dealt with
1762  unsigned n_zeta=Received_flat_packed_zetas_to_be_found.size()/Dim;
1763 
1764 
1765 
1766 
1767  // Create storage for the processor id (plus one) on which
1768  // the zetas stored in Flat_packed_zetas_not_found_locally[...]
1769  // were located
1770  Proc_id_plus_one_of_external_element.resize(n_zeta);
1771 
1772  // Create storage for the status of the (external halo) element associated
1773  // the zetas stored in Flat_packed_zetas_not_found_locally[...].
1774  // It either hasn't been found, already exists on the processor
1775  // that needs it, or needs to be newly created.
1776  Located_element_status.resize(n_zeta);
1777 
1778  // Counter for flat-packed array of external zeta coordinates
1779  unsigned count=0;
1780 
1781  // Loop over the zeta tuples that we received from elsewhere and
1782  // are trying to find here
1783  for (unsigned i=0;i<n_zeta;i++)
1784  {
1785  // Storage for global coordinates to be located
1786  Vector<double> x_global(Dim);
1787 
1788  // Loop to fill in coordinates
1789  for (unsigned ii=0;ii<Dim;ii++)
1790  {
1791  x_global[ii]=Received_flat_packed_zetas_to_be_found[count];
1792  count++;
1793  }
1794 
1795  // Perform locate_zeta for these coordinates
1796  GeomObject *sub_geom_obj_pt;
1797  Vector<double> ss(Dim);
1798  bool called_within_spiral=true;
1799  mesh_geom_obj_pt->spiraling_locate_zeta(x_global,sub_geom_obj_pt,ss,
1800  called_within_spiral);
1801 
1802  // Did the locate method work?
1803  if (sub_geom_obj_pt!=0)
1804  {
1805  // Get the source element - bulk or not?
1806  GeneralisedElement *source_el_pt=0;
1808  {
1809  source_el_pt=dynamic_cast<FiniteElement*>(sub_geom_obj_pt);
1810  }
1811  else
1812  {
1813  FaceElement *face_el_pt=dynamic_cast<FaceElement*>(sub_geom_obj_pt);
1814  source_el_pt=dynamic_cast<FiniteElement*>(face_el_pt->
1815  bulk_element_pt());
1816  }
1817 
1818  // Check if the returned element is halo
1819  if (!source_el_pt->is_halo()) // cannot accept halo here
1820  {
1821  // The correct non-halo element has been located; this will become
1822  // an external haloed element on the current process, and an
1823  // external halo copy needs to be created on the current process
1824  // minus wherever we are in the "ring-loop"
1825  int halo_copy_proc=my_rank-iproc;
1826 
1827  // If iproc is bigger than my_rank then we've "gone through" nproc-1
1828  if (my_rank<iproc) { halo_copy_proc=n_proc+halo_copy_proc; }
1829 
1830  // So, we found zeta on the current processor
1832 
1833  // This source element is an external halo on process halo_copy_proc
1834  // but it should only be added to the storage if it hasn't
1835  // been added already, and this information also needs to be
1836  // communicated over to the other process
1837 
1838  unsigned n_extern_haloed=external_mesh_pt->
1839  nexternal_haloed_element(halo_copy_proc);
1840  unsigned external_haloed_el_index=
1841  external_mesh_pt->add_external_haloed_element_pt(halo_copy_proc,
1842  source_el_pt);
1843 
1844  // If it was added to the storage then the returned index
1845  // will be the same as the (old) size of the storage
1846  if (external_haloed_el_index==n_extern_haloed)
1847  {
1848  // Set Located_element_status to say it
1849  // should be newly created
1851 
1852  // How many continuously interpolated values are there?
1853  int n_cont_inter_values=-1;
1854  if (dynamic_cast<RefineableElement*>(source_el_pt)!=0)
1855  {
1856  n_cont_inter_values=dynamic_cast<RefineableElement*>
1857  (source_el_pt)->ncont_interpolated_values();
1858  }
1859 
1860  // Since it is (externally) haloed from the current process,
1861  // the info required to create a new element in the equivalent
1862  // external halo layer on process halo_copy_proc needs to be
1863  // sent there
1864 
1865  // If we're using macro elements to update...
1866  MacroElementNodeUpdateMesh* macro_mesh_pt=
1867  dynamic_cast<MacroElementNodeUpdateMesh*>(external_mesh_pt);
1868  if (macro_mesh_pt!=0)
1869  {
1870  Flat_packed_unsigneds.push_back(1);
1871 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1872  Flat_packed_unsigneds_string.push_back("Mesh is macro element mesh[2]");
1873 #endif
1874  //Cast to finite element... this must work because it's
1875  //a macroelement no update mesh
1876  FiniteElement* source_finite_el_pt
1877  = dynamic_cast<FiniteElement*>(source_el_pt);
1878 
1879  MacroElement* macro_el_pt=source_finite_el_pt->macro_elem_pt();
1880  // Send the macro element number across
1881  unsigned macro_el_num=macro_el_pt->macro_element_number();
1882  Flat_packed_unsigneds.push_back(macro_el_num);
1883 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1884  Flat_packed_unsigneds_string.push_back("Number of macro element[2]");
1885 #endif
1886  // we need to send
1887  // the lower left and upper right coordinates of the macro
1888  QElementBase* q_el_pt=dynamic_cast<QElementBase*>(source_el_pt);
1889  if (q_el_pt!=0)
1890  {
1891  // The macro element needs to be set first before
1892  // its lower left and upper right coordinates can be accessed
1893  // Now send the lower left and upper right coordinates
1894  unsigned el_dim=q_el_pt->dim();
1895  for (unsigned i_dim=0;i_dim<el_dim;i_dim++)
1896  {
1897  Flat_packed_doubles.push_back(q_el_pt->s_macro_ll(i_dim));
1898  Flat_packed_doubles.push_back(q_el_pt->s_macro_ur(i_dim));
1899  }
1900  }
1901  else // Throw an error
1902  {
1903  std::ostringstream error_stream;
1904  error_stream << "You are using a MacroElement node update\n"
1905  << "in a case with non-QElements. This has not\n"
1906  << "yet been implemented.\n";
1907  throw OomphLibError
1908  (error_stream.str(),
1909  OOMPH_CURRENT_FUNCTION,
1910  OOMPH_EXCEPTION_LOCATION);
1911  }
1912 
1913  }
1914  else // Not using macro elements to update
1915  {
1916  Flat_packed_unsigneds.push_back(0);
1917 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1918  Flat_packed_unsigneds_string.push_back("Mesh is not a macro element mesh [2]");
1919 #endif
1920  }
1921 
1922 
1923  //Cast to finite element... this must work because it's
1924  //a macroelement no update mesh
1925  FiniteElement* source_finite_el_pt
1926  = dynamic_cast<FiniteElement*>(source_el_pt);
1927 #ifdef PARANOID
1928  if(source_finite_el_pt==0)
1929  {
1930  throw
1931  OomphLibError(
1932  "Unable to cast source function to finite element\n",
1933  "Multi_domain_functions::locate_zeta_for_missing_coordinates()",
1934  OOMPH_EXCEPTION_LOCATION);
1935  }
1936 #endif
1937 
1938 
1939  // Loop over the nodes of the new source element
1940  unsigned n_node=source_finite_el_pt->nnode();
1941  for (unsigned j=0;j<n_node;j++)
1942  {
1943  Node* nod_pt=source_finite_el_pt->node_pt(j);
1944 
1945  // Add the node to the storage; this routine
1946  // also takes care of any master nodes if the
1947  // node is hanging
1948  add_external_haloed_node_to_storage(halo_copy_proc,nod_pt,
1949  problem_pt,
1950  external_mesh_pt,
1951  n_cont_inter_values);
1952  }
1953  }
1954  else // it has already been added, so tell the other process
1955  {
1956  // Set Located_element_status to indicate an element has
1957  // already been added
1959  Flat_packed_unsigneds.push_back(external_haloed_el_index);
1960 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1961  Flat_packed_unsigneds_string.push_back("Index of existing external haloed element[2]");
1962 #endif
1963  }
1964 
1965  // The coordinates returned by locate_zeta are also needed
1966  // in the setup of the source elements on the other process
1968  {
1969  for (unsigned ii=0;ii<Dim;ii++)
1970  {
1971  Flat_packed_located_coordinates.push_back(ss[ii]);
1972  }
1973  }
1974  else // translate the coordinates to the bulk element
1975  {
1976  // The translation is from Lagrangian to Eulerian
1977  FaceElement *face_el_pt=
1978  dynamic_cast<FaceElement*>(sub_geom_obj_pt);
1979  //Get the dimension of the BulkElement
1980  unsigned bulk_el_dim =
1981  dynamic_cast<FiniteElement*>(source_el_pt)->dim();
1982  Vector<double> s_trans(bulk_el_dim);
1983  face_el_pt->get_local_coordinate_in_bulk(ss,s_trans);
1984  for (unsigned ii=0;ii<bulk_el_dim;ii++)
1985  {
1986  Flat_packed_located_coordinates.push_back(s_trans[ii]);
1987  }
1988  }
1989  }
1990  else // halo, so search again until non-halo equivalent is located
1991  {
1992  // Add required information to arrays (as below)
1993  for (unsigned ii=0;ii<Dim;ii++)
1994  {
1995  Flat_packed_zetas_not_found_locally.push_back(x_global[ii]);
1996  }
1997  // It wasn't found here
1999 
2000  // Set Located_element_status to indicate not found
2002  }
2003  }
2004  else // not successful this time, so prepare for next process to try
2005  {
2006  // Add this global coordinate to the LOCAL zeta array
2007  for (unsigned ii=0;ii<Dim;ii++)
2008  {
2009  Flat_packed_zetas_not_found_locally.push_back(x_global[ii]);
2010  }
2011  // It wasn't found here
2013 
2014  // Set Located_element_status to indicate not found
2016  }
2017  }
2018 
2019  }
2020 
2021  //=====================================================================
2022  /// Locate zeta for current set of missing coordinates; vector-based version
2023  //=====================================================================
2025  (int& iproc, Mesh* const &external_mesh_pt, Problem* problem_pt,
2026  Vector<MeshAsGeomObject*>& mesh_geom_obj_pt)
2027  {
2028  // How many meshes are we dealing with?
2029  unsigned n_mesh=mesh_geom_obj_pt.size();
2030 
2031  // Storage for number of processors, current process and communicator
2032  OomphCommunicator* comm_pt=problem_pt->communicator_pt();
2033  int n_proc=comm_pt->nproc();
2034  int my_rank=comm_pt->my_rank();
2035 
2036  // Clear vectors containing data to be sent
2037  Flat_packed_doubles.resize(0);
2038 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2039  Flat_packed_unsigneds_string.resize(0);
2040 #endif
2041  Flat_packed_unsigneds.resize(0);
2043 
2044  // Flush storage for zetas not found locally (when
2045  // processing the zeta coordinates received from "previous"
2046  // processor)
2048 
2049  // Number of zeta tuples to be dealt with (includes padding!)
2050  unsigned n_zeta=Received_flat_packed_zetas_to_be_found.size()/Dim;
2051 
2052  // Create storage for the processor id (plus one) on which
2053  // the zetas stored in Flat_packed_zetas_not_found_locally[...]
2054  // were located. (Remains zero for padded entries).
2055  Proc_id_plus_one_of_external_element.resize(n_zeta,0);
2056 
2057  // Create storage for the status of the (external halo) element associated
2058  // the zetas stored in Flat_packed_zetas_not_found_locally[...].
2059  // It either hasn't been found, already exists on the processor
2060  // that needs it, or needs to be newly created. (Remains Not_found
2061  // for padded entries).
2062  Located_element_status.resize(n_zeta,Not_found);
2063 
2064  // Counter for flat-packed array of external zeta coordinates
2065  unsigned count=0;
2066 
2067  // Current mesh
2068  unsigned i_mesh=0;
2069 
2070  // Loop over the zeta tuples that we received from elsewhere and
2071  // are trying to find here for current mesh
2072  for (unsigned i=0;i<n_zeta;i++)
2073  {
2074  // Storage for global coordinates to be located
2075  Vector<double> x_global(Dim);
2076 
2077  // Loop to fill in coordinates
2078  for (unsigned ii=0;ii<Dim;ii++)
2079  {
2080  x_global[ii]=Received_flat_packed_zetas_to_be_found[count];
2081  count++;
2082  }
2083 
2084  // Check if we've reached the end of the mesh
2085  bool reached_end_of_mesh=false;
2086  unsigned dbl_max_count=0;
2087  for (unsigned ii=0;ii<Dim;ii++)
2088  {
2089  if (x_global[ii]==DBL_MAX)
2090  {
2091  dbl_max_count++;
2092  reached_end_of_mesh=true;
2093  }
2094  }
2095 
2096  // Reached end of mesh
2097  if (reached_end_of_mesh)
2098  {
2099 #ifdef PARANOID
2100  // Check if all coordinates were set to DBX_MAX
2101  if (dbl_max_count!=Dim)
2102  {
2103  std::ostringstream error_stream;
2104  error_stream << "Appear to have reached end of mesh " << i_mesh
2105  << " but only " << dbl_max_count << " out of "
2106  << Dim << " zeta coordinates have been set to DBX_MAX\n";
2107  throw OomphLibError
2108  (error_stream.str(),
2109  OOMPH_CURRENT_FUNCTION,
2110  OOMPH_EXCEPTION_LOCATION);
2111  }
2112 #endif
2113  // Indicate end of mesh in flat packed data
2114  for (unsigned ii=0;ii<Dim;ii++)
2115  {
2116  Flat_packed_zetas_not_found_locally.push_back(DBL_MAX);
2117  }
2118 
2119  // Bump mesh counter
2120  i_mesh++;
2121 
2122  // Bail out if we're done
2123  if (i_mesh==n_mesh)
2124  {
2125  return;
2126  }
2127  }
2128 
2129  // Perform locate_zeta for these coordinates and current mesh
2130  GeomObject *sub_geom_obj_pt=0;
2131  Vector<double> ss(Dim);
2132  if (!reached_end_of_mesh)
2133  {
2134  bool called_within_spiral=true;
2135  mesh_geom_obj_pt[i_mesh]->spiraling_locate_zeta(x_global,
2136  sub_geom_obj_pt,ss,
2137  called_within_spiral);
2138 
2139 
2140  // Did the locate method work?
2141  if (sub_geom_obj_pt!=0)
2142  {
2143  // Get the source element - bulk or not?
2144  GeneralisedElement *source_el_pt=0;
2146  {
2147  source_el_pt=dynamic_cast<FiniteElement*>(sub_geom_obj_pt);
2148  }
2149  else
2150  {
2151  FaceElement *face_el_pt=dynamic_cast<FaceElement*>(sub_geom_obj_pt);
2152  source_el_pt=dynamic_cast<FiniteElement*>(face_el_pt->
2153  bulk_element_pt());
2154  }
2155 
2156  // Check if the returned element is halo
2157  if (!source_el_pt->is_halo()) // cannot accept halo here
2158  {
2159  // The correct non-halo element has been located; this will become
2160  // an external haloed element on the current process, and an
2161  // external halo copy needs to be created on the current process
2162  // minus wherever we are in the "ring-loop"
2163  int halo_copy_proc=my_rank-iproc;
2164 
2165  // If iproc is bigger than my_rank then we've "gone through" nproc-1
2166  if (my_rank<iproc) { halo_copy_proc=n_proc+halo_copy_proc; }
2167 
2168  // So, we found zeta on the current processor
2170 
2171  // This source element is an external halo on process halo_copy_proc
2172  // but it should only be added to the storage if it hasn't
2173  // been added already, and this information also needs to be
2174  // communicated over to the other process
2175 
2176  unsigned n_extern_haloed=external_mesh_pt->
2177  nexternal_haloed_element(halo_copy_proc);
2178  unsigned external_haloed_el_index=
2179  external_mesh_pt->add_external_haloed_element_pt(halo_copy_proc,
2180  source_el_pt);
2181 
2182  // If it was added to the storage then the returned index
2183  // will be the same as the (old) size of the storage
2184  if (external_haloed_el_index==n_extern_haloed)
2185  {
2186  // Set Located_element_status to say it
2187  // should be newly created
2189 
2190  // How many continuously interpolated values are there?
2191  int n_cont_inter_values=-1;
2192  if (dynamic_cast<RefineableElement*>(source_el_pt)!=0)
2193  {
2194  n_cont_inter_values=dynamic_cast<RefineableElement*>
2195  (source_el_pt)->ncont_interpolated_values();
2196  }
2197 
2198  // Since it is (externally) haloed from the current process,
2199  // the info required to create a new element in the equivalent
2200  // external halo layer on process halo_copy_proc needs to be
2201  // sent there
2202 
2203  // If we're using macro elements to update...
2204  MacroElementNodeUpdateMesh* macro_mesh_pt=
2205  dynamic_cast<MacroElementNodeUpdateMesh*>(external_mesh_pt);
2206  if (macro_mesh_pt!=0)
2207  {
2208  Flat_packed_unsigneds.push_back(1);
2209 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2210  Flat_packed_unsigneds_string.push_back("Mesh is macro element mesh[2]");
2211 #endif
2212  //Cast to finite element... this must work because it's
2213  //a macroelement no update mesh
2214  FiniteElement* source_finite_el_pt
2215  = dynamic_cast<FiniteElement*>(source_el_pt);
2216 
2217  MacroElement* macro_el_pt=source_finite_el_pt->macro_elem_pt();
2218  // Send the macro element number across
2219  unsigned macro_el_num=macro_el_pt->macro_element_number();
2220  Flat_packed_unsigneds.push_back(macro_el_num);
2221 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2222  Flat_packed_unsigneds_string.push_back("Number of macro element[2]");
2223 #endif
2224  // we need to send
2225  // the lower left and upper right coordinates of the macro
2226  QElementBase* q_el_pt=dynamic_cast<QElementBase*>(source_el_pt);
2227  if (q_el_pt!=0)
2228  {
2229  // The macro element needs to be set first before
2230  // its lower left and upper right coordinates can be accessed
2231  // Now send the lower left and upper right coordinates
2232  unsigned el_dim=q_el_pt->dim();
2233  for (unsigned i_dim=0;i_dim<el_dim;i_dim++)
2234  {
2235  Flat_packed_doubles.push_back(q_el_pt->s_macro_ll(i_dim));
2236  Flat_packed_doubles.push_back(q_el_pt->s_macro_ur(i_dim));
2237  }
2238  }
2239  else // Throw an error
2240  {
2241  std::ostringstream error_stream;
2242  error_stream << "You are using a MacroElement node update\n"
2243  << "in a case with non-QElements. This has not\n"
2244  << "yet been implemented.\n";
2245  throw OomphLibError
2246  (error_stream.str(),
2247  OOMPH_CURRENT_FUNCTION,
2248  OOMPH_EXCEPTION_LOCATION);
2249  }
2250 
2251  }
2252  else // Not using macro elements to update
2253  {
2254  Flat_packed_unsigneds.push_back(0);
2255 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2256  Flat_packed_unsigneds_string.push_back("Mesh is not a macro element mesh [2]");
2257 #endif
2258  }
2259 
2260 
2261  //Cast to finite element... this must work because it's
2262  //a macroelement no update mesh
2263  FiniteElement* source_finite_el_pt
2264  = dynamic_cast<FiniteElement*>(source_el_pt);
2265 #ifdef PARANOID
2266  if(source_finite_el_pt==0)
2267  {
2268  throw
2269  OomphLibError(
2270  "Unable to cast source function to finite element\n",
2271  "Multi_domain_functions::locate_zeta_for_missing_coordinates()",
2272  OOMPH_EXCEPTION_LOCATION);
2273  }
2274 #endif
2275 
2276 
2277  // Loop over the nodes of the new source element
2278  unsigned n_node=source_finite_el_pt->nnode();
2279  for (unsigned j=0;j<n_node;j++)
2280  {
2281  Node* nod_pt=source_finite_el_pt->node_pt(j);
2282 
2283  // Add the node to the storage; this routine
2284  // also takes care of any master nodes if the
2285  // node is hanging
2286  add_external_haloed_node_to_storage(halo_copy_proc,nod_pt,
2287  problem_pt,
2288  external_mesh_pt,
2289  n_cont_inter_values);
2290  }
2291  }
2292  else // it has already been added, so tell the other process
2293  {
2294  // Set Located_element_status to indicate an element has
2295  // already been added
2297  Flat_packed_unsigneds.push_back(external_haloed_el_index);
2298 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2299  Flat_packed_unsigneds_string.push_back("Index of existing external haloed element[2]");
2300 #endif
2301  }
2302 
2303  // The coordinates returned by locate_zeta are also needed
2304  // in the setup of the source elements on the other process
2306  {
2307  for (unsigned ii=0;ii<Dim;ii++)
2308  {
2309  Flat_packed_located_coordinates.push_back(ss[ii]);
2310  }
2311  }
2312  else // translate the coordinates to the bulk element
2313  {
2314  // The translation is from Lagrangian to Eulerian
2315  FaceElement *face_el_pt=
2316  dynamic_cast<FaceElement*>(sub_geom_obj_pt);
2317  //Get the dimension of the BulkElement
2318  unsigned bulk_el_dim =
2319  dynamic_cast<FiniteElement*>(source_el_pt)->dim();
2320  Vector<double> s_trans(bulk_el_dim);
2321  face_el_pt->get_local_coordinate_in_bulk(ss,s_trans);
2322  for (unsigned ii=0;ii<bulk_el_dim;ii++)
2323  {
2324  Flat_packed_located_coordinates.push_back(s_trans[ii]);
2325  }
2326  }
2327  }
2328  else // halo, so search again until non-halo equivalent is located
2329  {
2330  // Add required information to arrays (as below)
2331  for (unsigned ii=0;ii<Dim;ii++)
2332  {
2333  Flat_packed_zetas_not_found_locally.push_back(x_global[ii]);
2334  }
2335  // It wasn't found here
2337 
2338  // Set Located_element_status to indicate not found
2340  }
2341  }
2342  else // not successful this time (i.e. sub_geom_obj_pt==0), so
2343  // prepare for next process to try
2344  {
2345  // Add this global coordinate to the LOCAL zeta array
2346  for (unsigned ii=0;ii<Dim;ii++)
2347  {
2348  Flat_packed_zetas_not_found_locally.push_back(x_global[ii]);
2349  }
2350  // It wasn't found here
2352 
2353  // Set Located_element_status to indicate not found
2355  }
2356 
2357  } // end of mesh not reached
2358 
2359  } // end of loop over flat-packed zeta tuples
2360 
2361  }
2362 
2363 
2364 
2365 #endif
2366 
2367 
2368 
2369  //=====================================================================
2370  /// locate zeta for current set of "local" coordinates
2371  //=====================================================================
2373  (Mesh* const &mesh_pt, Mesh* const &external_mesh_pt,
2374  MeshAsGeomObject* &mesh_geom_obj_pt,
2375  const unsigned& interaction_index)
2376  {
2377  // Number of local elements
2378  unsigned n_element=mesh_pt->nelement();
2379 
2380  // Flush storage for zetas not found locally
2382 
2383  // Loop over this processor's elements
2384  for (unsigned e=0;e<n_element;e++)
2385  {
2387  dynamic_cast<ElementWithExternalElement*>(mesh_pt->element_pt(e));
2388 #ifdef OOMPH_HAS_MPI
2389  // Only visit non-halo elements -- we're not setting up external elements
2390  // for on-halos!
2391  if (!el_pt->is_halo())
2392 #endif
2393  {
2394  // Find number of Gauss points and element dimension
2395  unsigned n_intpt=el_pt->integral_pt()->nweight();
2396  unsigned el_dim=el_pt->dim();
2397 
2398  // Set storage for local and global coordinates
2399  Vector<double> s_local(el_dim);
2400  Vector<double> x_global(el_dim);
2401 
2402  // Loop over integration points
2403  for (unsigned ipt=0;ipt<n_intpt;ipt++)
2404  {
2405  // Has this integration point been done yet?
2406  if (External_element_located[e][ipt]==0)
2407  {
2408  // Get local coordinates
2409  for (unsigned i=0;i<el_dim;i++)
2410  {
2411  s_local[i]=el_pt->integral_pt()->knot(ipt,i);
2412  }
2413  // Interpolate to global coordinates
2414  el_pt->interpolated_zeta(s_local,x_global);
2415 
2416  // Storage for geometric object and its local coordinates
2417  GeomObject* sub_geom_obj_pt=0;
2418  Vector<double> s_ext(el_dim);
2419 
2420  // Perform locate_zeta locally for this coordinate
2421  bool called_within_spiral=true;
2422 
2423  // Which version of spiraling locate zeta should we call?
2424  // The logic is to call the specialised version
2425  // "my_spiraling_locate_zeta()" only if we are running in
2426  // serial AND doing the setup of multidomain for
2427  // projection, any other case we call the original
2428  // "spiraling_locate_zeta method()"
2429 #ifdef OOMPH_HAS_MPI
2430  // We need to check this to avoid segfault when checking
2431  // for number of processors if mpi has not been initialised
2433  {
2434  if (MPI_Helpers::communicator_pt()->nproc()>1)
2435  {
2436  mesh_geom_obj_pt->spiraling_locate_zeta(x_global,
2437  sub_geom_obj_pt,s_ext,
2438  called_within_spiral);
2439  }
2440  else
2441  {
2443  {
2444  mesh_geom_obj_pt->spiraling_locate_zeta(x_global,
2445  sub_geom_obj_pt,s_ext,
2446  called_within_spiral);
2447  }
2448  else
2449  {
2450  mesh_geom_obj_pt->
2451  my_spiraling_locate_zeta(x_global,
2452  sub_geom_obj_pt,s_ext,
2453  called_within_spiral);
2454  }
2455  }
2456  }
2457  else
2458  {
2460  {
2461  mesh_geom_obj_pt->spiraling_locate_zeta(x_global,
2462  sub_geom_obj_pt,s_ext,
2463  called_within_spiral);
2464  }
2465  else
2466  {
2467  mesh_geom_obj_pt->my_spiraling_locate_zeta(x_global,
2468  sub_geom_obj_pt,s_ext,
2469  called_within_spiral);
2470  }
2471  }
2472 #else
2474  {
2475  mesh_geom_obj_pt->spiraling_locate_zeta(x_global,
2476  sub_geom_obj_pt,s_ext,
2477  called_within_spiral);
2478  }
2479  else
2480  {
2481  mesh_geom_obj_pt->my_spiraling_locate_zeta(x_global,
2482  sub_geom_obj_pt,s_ext,
2483  called_within_spiral);
2484  }
2485 #endif // #ifdef OOMPH_HAS_MPI
2486 
2487  // Has the required element been located?
2488  if (sub_geom_obj_pt!=0)
2489  {
2490  // The required element has been located
2491  // The located coordinates have the same dimension as the bulk
2492  GeneralisedElement* source_el_pt;
2493  Vector<double> s_source(el_dim);
2494 
2495  // Is the bulk element the actual external element?
2497  {
2498  // Use the object directly (it must be a finite element)
2499  source_el_pt=dynamic_cast<FiniteElement*>(sub_geom_obj_pt);
2500  s_source=s_ext;
2501  }
2502  else
2503  {
2504  // Cast to a FaceElement and use the bulk element
2505  FaceElement* face_el_pt=
2506  dynamic_cast<FaceElement*>(sub_geom_obj_pt);
2507  source_el_pt=face_el_pt->bulk_element_pt();
2508 
2509  //Need to resize the located coordinates to have the same
2510  //dimension as the bulk element
2511  s_source.resize(dynamic_cast<FiniteElement*>
2512  (source_el_pt)->dim());
2513 
2514  // Translate the returned local coords into the bulk element
2515  face_el_pt->get_local_coordinate_in_bulk(s_ext,s_source);
2516  }
2517 
2518  // Check if it's a halo; if it is then the non-halo equivalent
2519  // needs to be located from another processor (unless we
2520  // accept halo elements as external elements
2521 #ifdef OOMPH_HAS_MPI
2523  (!source_el_pt->is_halo()))
2524 #endif
2525  {
2526  //Need to cast to a FiniteElement
2527  FiniteElement* source_finite_el_pt =
2528  dynamic_cast<FiniteElement*>(source_el_pt);
2529 
2530  // Set the external element pointer and local coordinates
2531  el_pt->external_element_pt(interaction_index,ipt)
2532  = source_finite_el_pt;
2533  el_pt->external_element_local_coord(interaction_index,ipt)
2534  =s_source;
2535 
2536  // Set the lookup array to 1/true
2537  External_element_located[e][ipt]=1;
2538  }
2539 #ifdef OOMPH_HAS_MPI
2540  else // elements can only be halo if MPI is turned on
2541  {
2542  // Add required information to arrays
2543  for (unsigned i=0;i<el_dim;i++)
2544  {
2545  Flat_packed_zetas_not_found_locally.push_back(x_global[i]);
2546  }
2547  }
2548 #endif
2549  }
2550  else
2551  {
2552  // If it has failed then add the required information to the
2553  // arrays which need to be sent to the other processors so that
2554  // they can perform the locate_zeta
2555 
2556  // Add this global coordinate to the LOCAL zeta array
2557  for (unsigned i=0;i<el_dim;i++)
2558  {
2559  Flat_packed_zetas_not_found_locally.push_back(x_global[i]);
2560  }
2561  }
2562  }
2563  } // end loop over integration points
2564  } //end for halo
2565  } // end loop over local elements
2566  }
2567 
2568 
2569 
2570  //=====================================================================
2571  /// locate zeta for current set of "local" coordinates
2572  /// vector-based version
2573  //=====================================================================
2575  (const Vector<Mesh*>& mesh_pt, Mesh* const &external_mesh_pt,
2576  Vector<MeshAsGeomObject*>& mesh_geom_obj_pt,
2577  const unsigned& interaction_index)
2578  {
2579  // Flush storage for zetas not found locally
2581 
2582  // Number of meshes
2583  unsigned n_mesh=mesh_pt.size();
2584 
2585 #ifdef PARANOID
2586  if (mesh_geom_obj_pt.size()!=n_mesh)
2587  {
2588  std::ostringstream error_stream;
2589  error_stream << "Sizes of mesh_geom_obj_pt [ "
2590  << mesh_geom_obj_pt.size() << " ] and "
2591  << "mesh_pt [ " << n_mesh << " ] don't match.\n";
2592  throw OomphLibError
2593  (error_stream.str(),
2594  OOMPH_CURRENT_FUNCTION,
2595  OOMPH_EXCEPTION_LOCATION);
2596  }
2597 #endif
2598 
2599  // Element counter
2600  unsigned e_count=0;
2601 
2602  // Loop over meshes
2603  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
2604  {
2605  // Number of local elements
2606  unsigned n_element=mesh_pt[i_mesh]->nelement();
2607 
2608  // Loop over this processor's elements
2609  for (unsigned e=0;e<n_element;e++)
2610  {
2612  dynamic_cast<ElementWithExternalElement*>(mesh_pt[i_mesh]->
2613  element_pt(e));
2614 #ifdef OOMPH_HAS_MPI
2615  // Only visit non-halo elements -- we're not setting up external elements
2616  // for on-halos!
2617  if (!el_pt->is_halo())
2618 #endif
2619  {
2620  // Find number of Gauss points and element dimension
2621  unsigned n_intpt=el_pt->integral_pt()->nweight();
2622  unsigned el_dim=el_pt->dim();
2623 
2624 
2625 #ifdef PARANOID
2626  if (el_dim!=Dim)
2627  {
2628  std::ostringstream error_stream;
2629  error_stream
2630  << "Dimension of element " << el_dim
2631  << " is not consitent with dimension assumed \n"
2632  << " in multidomain namespace, " << Dim << std::endl;
2633  throw OomphLibError
2634  (error_stream.str(),
2635  OOMPH_CURRENT_FUNCTION,
2636  OOMPH_EXCEPTION_LOCATION);
2637  }
2638 #endif
2639 
2640  // Set storage for local and global coordinates
2641  Vector<double> s_local(el_dim);
2642  Vector<double> x_global(el_dim);
2643 
2644  // Loop over integration points
2645  for (unsigned ipt=0;ipt<n_intpt;ipt++)
2646  {
2647  // Has this integration point been done yet?
2648  if (External_element_located[e_count][ipt]==0)
2649  {
2650  // Get local coordinates
2651  for (unsigned i=0;i<el_dim;i++)
2652  {
2653  s_local[i]=el_pt->integral_pt()->knot(ipt,i);
2654  }
2655  // Interpolate to global coordinates
2656  el_pt->interpolated_zeta(s_local,x_global);
2657 
2658  // Storage for geometric object and its local coordinates
2659  GeomObject* sub_geom_obj_pt=0;
2660  Vector<double> s_ext(el_dim);
2661 
2662  // Perform locate_zeta locally for this coordinate
2663  bool called_within_spiral=true;
2664 
2665  // Which version of spiraling locate zeta should we call?
2666  // The logic is to call the specialised version
2667  // "my_spiraling_locate_zeta()" only if we are running in
2668  // serial AND doing the setup of multidomain for
2669  // projection, any other case we call the original
2670  // "spiraling_locate_zeta method()"
2671 #ifdef OOMPH_HAS_MPI
2672  // We need to check this to avoid segfault when checking
2673  // for number of processors if mpi has not been
2674  // initialised
2676  {
2677  if (MPI_Helpers::communicator_pt()->nproc()>1)
2678  {
2679  mesh_geom_obj_pt[i_mesh]->
2680  spiraling_locate_zeta(x_global,
2681  sub_geom_obj_pt,s_ext,
2682  called_within_spiral);
2683  }
2684  else
2685  {
2687  {
2688  mesh_geom_obj_pt[i_mesh]->
2689  spiraling_locate_zeta(x_global,
2690  sub_geom_obj_pt,s_ext,
2691  called_within_spiral);
2692  }
2693  else
2694  {
2695  mesh_geom_obj_pt[i_mesh]->
2696  my_spiraling_locate_zeta(x_global,
2697  sub_geom_obj_pt,s_ext,
2698  called_within_spiral);
2699  }
2700  }
2701  }
2702  else
2703  {
2705  {
2706  mesh_geom_obj_pt[i_mesh]->
2707  spiraling_locate_zeta(x_global,
2708  sub_geom_obj_pt,s_ext,
2709  called_within_spiral);
2710  }
2711  else
2712  {
2713  mesh_geom_obj_pt[i_mesh]->
2714  my_spiraling_locate_zeta(x_global,
2715  sub_geom_obj_pt,s_ext,
2716  called_within_spiral);
2717  }
2718  }
2719 #else
2721  {
2722  mesh_geom_obj_pt[i_mesh]->
2723  spiraling_locate_zeta(x_global,
2724  sub_geom_obj_pt,s_ext,
2725  called_within_spiral);
2726  }
2727  else
2728  {
2729  mesh_geom_obj_pt[i_mesh]->
2730  my_spiraling_locate_zeta(x_global,
2731  sub_geom_obj_pt,s_ext,
2732  called_within_spiral);
2733  }
2734 #endif // #ifdef OOMPH_HAS_MPI
2735 
2736  // Which version of spiraling locate zeta should we call?
2737  // The logic is to call the specialised version
2738  // "my_spiraling_locate_zeta()" only if we are running in
2739  // serial AND doing the setup of multidomain for
2740  // projection, any other case we call the original
2741  // "spiraling_locate_zeta method()"
2742 #ifdef OOMPH_HAS_MPI
2743  // We need to check this to avoid segfault when checking
2744  // for number of processors if mpi has not been
2745  // initialised
2747  {
2748  if (MPI_Helpers::communicator_pt()->nproc()>1)
2749  {
2750  mesh_geom_obj_pt[i_mesh]->
2751  spiraling_locate_zeta(x_global,
2752  sub_geom_obj_pt,s_ext,
2753  called_within_spiral);
2754  }
2755  else
2756  {
2758  {
2759  mesh_geom_obj_pt[i_mesh]->
2760  spiraling_locate_zeta(x_global,
2761  sub_geom_obj_pt,s_ext,
2762  called_within_spiral);
2763  }
2764  else
2765  {
2766  mesh_geom_obj_pt[i_mesh]->
2767  my_spiraling_locate_zeta(x_global,
2768  sub_geom_obj_pt,s_ext,
2769  called_within_spiral);
2770  }
2771  }
2772  }
2773  else
2774  {
2776  {
2777  mesh_geom_obj_pt[i_mesh]->
2778  spiraling_locate_zeta(x_global,
2779  sub_geom_obj_pt,s_ext,
2780  called_within_spiral);
2781  }
2782  else
2783  {
2784  mesh_geom_obj_pt[i_mesh]->
2785  my_spiraling_locate_zeta(x_global,
2786  sub_geom_obj_pt,s_ext,
2787  called_within_spiral);
2788  }
2789  }
2790 #else
2792  {
2793  mesh_geom_obj_pt[i_mesh]->
2794  spiraling_locate_zeta(x_global,
2795  sub_geom_obj_pt,s_ext,
2796  called_within_spiral);
2797  }
2798  else
2799  {
2800  mesh_geom_obj_pt[i_mesh]->
2801  my_spiraling_locate_zeta(x_global,
2802  sub_geom_obj_pt,s_ext,
2803  called_within_spiral);
2804  }
2805 #endif // #ifdef OOMPH_HAS_MPI
2806 
2807  // Has the required element been located?
2808  if (sub_geom_obj_pt!=0)
2809  {
2810  // The required element has been located
2811  // The located coordinates have the same dimension as the bulk
2812  GeneralisedElement* source_el_pt;
2813  Vector<double> s_source(el_dim);
2814 
2815  // Is the bulk element the actual external element?
2817  {
2818  // Use the object directly (it must be a finite element)
2819  source_el_pt=dynamic_cast<FiniteElement*>(sub_geom_obj_pt);
2820  s_source=s_ext;
2821  }
2822  else
2823  {
2824  // Cast to a FaceElement and use the bulk element
2825  FaceElement* face_el_pt=
2826  dynamic_cast<FaceElement*>(sub_geom_obj_pt);
2827  source_el_pt=face_el_pt->bulk_element_pt();
2828 
2829  //Need to resize the located coordinates to have the same
2830  //dimension as the bulk element
2831  s_source.resize(dynamic_cast<FiniteElement*>
2832  (source_el_pt)->dim());
2833 
2834  // Translate the returned local coords into the bulk element
2835  face_el_pt->get_local_coordinate_in_bulk(s_ext,s_source);
2836  }
2837 
2838  // Check if it's a halo; if it is then the non-halo equivalent
2839  // needs to be located from another processor (unless we
2840  // accept halo elements as external elements)
2841 #ifdef OOMPH_HAS_MPI
2843  (!source_el_pt->is_halo()))
2844 #endif
2845  {
2846  //Need to cast to a FiniteElement
2847  FiniteElement* source_finite_el_pt =
2848  dynamic_cast<FiniteElement*>(source_el_pt);
2849 
2850  // Set the external element pointer and local coordinates
2851  el_pt->external_element_pt(interaction_index,ipt)
2852  = source_finite_el_pt;
2853  el_pt->external_element_local_coord(interaction_index,ipt)
2854  =s_source;
2855 
2856  // Set the lookup array to 1/true
2857  External_element_located[e_count][ipt]=1;
2858  }
2859 #ifdef OOMPH_HAS_MPI
2860  // located element is halo and we're not accepting haloes
2861  // obviously only makes sense in mpi mode...
2862  else
2863  {
2864  // Add required information to arrays
2865  for (unsigned i=0;i<el_dim;i++)
2866  {
2867  Flat_packed_zetas_not_found_locally.push_back(x_global[i]);
2868  }
2869  }
2870 #endif
2871  }
2872  else
2873  {
2874  // Search has failed then add the required information to the
2875  // arrays which need to be sent to the other processors so that
2876  // they can perform the locate_zeta
2877 
2878  // Add this global coordinate to the LOCAL zeta array
2879  for (unsigned i=0;i<el_dim;i++)
2880  {
2881  Flat_packed_zetas_not_found_locally.push_back(x_global[i]);
2882  }
2883  }
2884  }
2885  } // end loop over integration points
2886  } //end for halo
2887 
2888  // Bump up counter for all elements
2889  e_count++;
2890 
2891  } // end loop over local elements
2892 
2893  // Mark end of mesh data in flat packed array
2894  for (unsigned i=0;i<Dim;i++)
2895  {
2896  Flat_packed_zetas_not_found_locally.push_back(DBL_MAX);
2897  }
2898 
2899  } // end of loop over meshes
2900  }
2901 
2902 
2903 //=====================================================================
2904 /// Helper function that computes the dimension of the elements within
2905 /// each of the specified meshes (and checks they are the same).
2906 /// Stores result in Dim.
2907 //=====================================================================
2908  void get_dim_helper(Problem* problem_pt, Mesh* const &mesh_pt,
2909  Mesh* const &external_mesh_pt)
2910  {
2911 #ifdef OOMPH_HAS_MPI
2912  // Storage for number of processors, current process and communicator
2913  OomphCommunicator* comm_pt=problem_pt->communicator_pt();
2914 #endif
2915 
2916  // Extract the element dimensions from the first element of each mesh
2917  unsigned mesh_dim=0;
2918  if (mesh_pt->nelement() > 0)
2919  {
2920  mesh_dim=
2921  dynamic_cast<FiniteElement*>(mesh_pt->element_pt(0))->dim();
2922  }
2923  unsigned external_mesh_dim=0;
2924  if (external_mesh_pt->nelement() > 0)
2925  {
2926  external_mesh_dim=
2927  dynamic_cast<FiniteElement*>(external_mesh_pt->element_pt(0))->dim();
2928  }
2929 
2930  // Need to do an Allreduce
2931 #ifdef OOMPH_HAS_MPI
2932  int n_proc=comm_pt->nproc();
2933  if (n_proc > 1)
2934  {
2935  unsigned mesh_dim_reduce;
2936  MPI_Allreduce(&mesh_dim,&mesh_dim_reduce,1,MPI_UNSIGNED,
2937  MPI_MAX,comm_pt->mpi_comm());
2938  mesh_dim=mesh_dim_reduce;
2939 
2940  unsigned external_mesh_dim_reduce;
2941  MPI_Allreduce(&external_mesh_dim,&external_mesh_dim_reduce,1,MPI_UNSIGNED,
2942  MPI_MAX,comm_pt->mpi_comm());
2943  external_mesh_dim=external_mesh_dim_reduce;
2944  }
2945 #endif
2946 
2947  // Check the dimensions are the same!
2948  if (mesh_dim!=external_mesh_dim)
2949  {
2950  std::ostringstream error_stream;
2951  error_stream << "The elements within the two meshes do not\n"
2952  << "have the same dimension, so the multi-domain\n"
2953  << "method will not work.\n"
2954  << "For the mesh, dim=" << mesh_dim
2955  << ", and the external mesh, dim=" << external_mesh_dim
2956  << "\n";
2957  throw OomphLibError
2958  (error_stream.str(),
2959  OOMPH_CURRENT_FUNCTION,
2960  OOMPH_EXCEPTION_LOCATION);
2961  }
2962 
2963  // Set dimension
2964  Dim=mesh_dim;
2965  }
2966 
2967 
2968  //=================================================================
2969  /// Helper function that clears all the information used
2970  /// during the external storage creation
2971  //=================================================================
2972  void clean_up()
2973  {
2977  Located_element_status.clear();
2979  Flat_packed_doubles.clear();
2980  Flat_packed_unsigneds.clear();
2981  External_element_located.clear();
2982  }
2983 
2984 
2985  /// \short Bool to decide if to sort entries in bin during locate_zeta
2986  /// operation (default: false)
2987  bool Sort_bin_entries=false;
2988 
2989  /// Vector of zeta coordinates that we're currently trying to locate;
2990  /// used in sorting of bin entries in further_away() comparison function
2992 
2993  /// Comparison function for sorting entries in bin: Returns true if
2994  /// point identified by p1 (comprising pointer to finite element and
2995  /// vector of local coordinates within that element) is closer to
2996  /// Zeta_coords_for_further_away_comparison than p2
2998  const std::pair<FiniteElement*,Vector<double> >& p1,
2999  const std::pair<FiniteElement*,Vector<double> >& p2)
3000  {
3001  // First point
3002  FiniteElement* el_pt=p1.first;
3003  Vector<double> s(p1.second);
3004  Vector<double> zeta(Dim);
3005  el_pt->interpolated_zeta(s,zeta);
3006  double dist_squared1=0.0;
3007  for (unsigned i=0;i<Dim;i++)
3008  {
3009  dist_squared1+=
3012  }
3013 
3014  // Second point
3015  el_pt=p2.first;
3016  s=p2.second;
3017  el_pt->interpolated_zeta(s,zeta);
3018  double dist_squared2=0.0;
3019  for (unsigned i=0;i<Dim;i++)
3020  {
3021  dist_squared2+=
3024  }
3025 
3026  // Which one is further
3027  if (dist_squared1<dist_squared2)
3028  {
3029  return true;
3030  }
3031  else
3032  {
3033  return false;
3034  }
3035  }
3036 }
3037 
3038 }
A Generalised Element class.
Definition: elements.h:76
void add_external_haloed_node_helper(int &iproc, Node *nod_pt, Problem *problem_pt, Mesh *const &external_mesh_pt, int &n_cont_inter_values)
Helper to add external haloed node that is not a master.
GeomObject * geom_object_pt(const unsigned &i)
Return pointer to i-th geometric object involved in default (usually first) update function...
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation...
Definition: nodes.cc:2314
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:246
unsigned & macro_element_number()
Access function to the Macro_element_number.
void add_external_halo_node_pt(const unsigned &p, Node *&nod_pt)
Add external halo node whose non-halo (external) counterpart is held on processor p to the storage sc...
Definition: mesh.h:2055
Vector< double > Received_flat_packed_zetas_to_be_found
Vector of flat-packed zeta coordinates for which the external element could not be found on another p...
Definition: multi_domain.cc:88
bool first_closer_than_second(const std::pair< FiniteElement *, Vector< double > > &p1, const std::pair< FiniteElement *, Vector< double > > &p2)
Comparison function for sorting entries in bin: Returns true if point identified by p1 (comprising po...
virtual void update_node_update(AlgebraicNode *&node_pt)=0
Update the node update info for given node, following mesh adaptation. Must be implemented for every ...
unsigned add_external_haloed_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add external haloed element whose non-halo counterpart is held on processor p to the storage scheme f...
Definition: mesh.cc:9107
void send_and_receive_located_info(int &iproc, Mesh *const &external_mesh_pt, Problem *problem_pt)
Helper function to send back any located information.
double Total_time_for_sorting_elements_in_bins
The timing for sorting the elements in the bins.
cstr elem_len * i
Definition: cfortran.h:607
unsigned Counter_for_flat_packed_unsigneds
Counter used when processing vector of flat-packed unsigneds – this is really "private" data...
unsigned nref_value(const int &id)
Number of reference values involved in id-th update function.
void add_external_haloed_master_node_helper(int &iproc, Node *master_nod_pt, Problem *problem_pt, Mesh *const &external_mesh_pt, int &n_cont_inter_values)
Helper function to add external haloed node that is a master.
void get_dim_helper(Problem *problem_pt, Mesh *const &mesh_pt, Mesh *const &external_mesh_pt)
Helper function that computes the dimension of the elements within each of the specified meshes (and ...
The Problem class.
Definition: problem.h:152
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
bool Doc_timings
Boolean to indicate whether to doc timings or not.
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary, so return false. This will be overloaded by BoundaryNodes.
Definition: nodes.h:1293
char t
Definition: cfortran.h:572
bool Accept_failed_locate_zeta_in_setup_multi_domain_interaction
Boolean to indicate that failure in setup multi domain functions is acceptable; defaults to false...
Definition: multi_domain.cc:62
virtual Node * construct_boundary_node(const unsigned &n)
Construct the local node n as a boundary node; that is a node that MAY be placed on a mesh boundary a...
Definition: elements.h:2408
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
static OomphCommunicator * communicator_pt()
access to the global oomph-lib communicator
void add_node_update_info(const int &id, AlgebraicMesh *mesh_pt, const Vector< GeomObject * > &geom_object_pt, const Vector< double > &ref_value, const bool &called_from_constructor=false)
Add algebraic update information for node: What's the ID of the mesh update function (typically used ...
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1148
OomphInfo oomph_info
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
std::ofstream Doc_boundary_coordinate_file
Output file to document the boundary coordinate along the mesh boundary of the bulk mesh during call ...
Definition: multi_domain.cc:53
static bool mpi_has_been_initialised()
return true if MPI has been initialised
bool Compute_extreme_bin_coordinates
Bool to tell the MeshAsGeomObject whether to calculate the extreme coordinates of the bin structure...
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
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1221
Vector< unsigned > Flat_packed_unsigneds
Vector of flat-packed unsigneds to be communicated with other processors – this is really "private" d...
GeomObject * geom_object_list_pt(const unsigned &i)
Access function to the ith GeomObject.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
double X_max
Maximum coordinate in first dimension.
double Y_max
Maximum coordinate in second dimension.
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:966
int node_update_fct_id()
Default (usually first if there are multiple ones) node update fct id.
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:66
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
unsigned ntime_stepper() const
Return the number of time steppers.
Definition: problem.h:1578
void locate_zeta_for_missing_coordinates(int &iproc, Mesh *const &external_mesh_pt, Problem *problem_pt, MeshAsGeomObject *&mesh_geom_obj_pt)
Locate zeta for current set of missing coordinates.
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
bool Doc_full_stats
Boolean to indicate whether to output further info during setup_multi_domain_interaction() routines...
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
Definition: nodes.h:1654
Base class for Qelements.
Definition: Qelements.h:106
double & s_macro_ur(const unsigned &i)
Access fct to the i-th coordinate of the element's "upper right" vertex in the associated MacroElemen...
Definition: Qelements.h:240
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
unsigned add_external_haloed_node_pt(const unsigned &p, Node *&nod_pt)
Add external haloed node whose halo (external) counterpart is held on processor p to the storage sche...
Definition: mesh.cc:9148
A class that contains the information required by Nodes that are located on Mesh boundaries. A BoundaryNode of a particular type is obtained by combining a given Node with this class. By differentiating between Nodes and BoundaryNodes we avoid a lot of un-necessary storage in the bulk Nodes.
Definition: nodes.h:1854
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...
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
void add_external_haloed_node_to_storage(int &iproc, Node *nod_pt, Problem *problem_pt, Mesh *const &external_mesh_pt, int &n_cont_inter_values)
Helper function to add external haloed nodes, including any masters.
bool Allow_use_of_halo_elements_as_external_elements
Boolean to indicate if we're allowed to use halo elements as external elements. Can drastically reduc...
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
Vector< double > Flat_packed_located_coordinates
Vector of flat-packed local coordinates for zeta tuples that have been located.
Vector< std::string > Flat_packed_unsigneds_string
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object. ...
Definition: elements.h:2380
Node *& external_halo_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th external halo node in this Mesh whose non-halo external counterpart is held on...
Definition: mesh.h:2064
static char t char * s
Definition: cfortran.h:572
unsigned nexternal_haloed_node()
Total number of external haloed nodes in this Mesh.
Definition: mesh.h:2099
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
void clean_up()
Helper function that clears all the intermediate information used during the external storage creatio...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
void get_required_master_nodal_information_helper(int &iproc, Node *master_nod_pt, Problem *problem_pt, Mesh *const &external_mesh_pt, int &n_cont_inter_values)
Helper function to get the required master nodal information from an external haloed master node so t...
void recursively_add_masters_of_external_haloed_node(int &iproc, Node *nod_pt, Problem *problem_pt, Mesh *const &external_mesh_pt, int &n_cont_inter_values)
Recursively add any master nodes (and their master nodes etc) of external haloed nodes.
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...
void get_local_coordinate_in_bulk(const Vector< double > &s, Vector< double > &s_bulk) const
Calculate the vector of local coordinate in the bulk element given the local coordinates in this Face...
Definition: elements.cc:6064
TimeStepper *& time_stepper_pt()
Access function for the pointer to the first (presumably only) timestepper.
Definition: problem.h:1460
void send_and_receive_missing_zetas(Problem *problem_pt)
Helper function to send any "missing" zeta coordinates to the next process and receive any coordinate...
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
FiniteElement *& node_update_element_pt()
Pointer to finite element that performs the update by referring to its macro-element representation (...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
Definition: elements.h:1829
Vector< double > Flat_packed_doubles
Vector of flat-packed doubles to be communicated with other processors.
bool Sort_bin_entries
Bool to decide if to sort entries in bin during locate_zeta operation (default: false) ...
Class that contains data for hanging nodes.
Definition: nodes.h:684
void construct_new_external_halo_node_helper(Node *&new_nod_pt, unsigned &loc_p, unsigned &node_index, FiniteElement *const &new_el_pt, Mesh *const &external_mesh_pt, Problem *problem_pt)
Helper function which constructs a new external halo node (on an element) with the information sent f...
void get_required_nodal_information_helper(int &iproc, Node *nod_pt, Problem *problem_pt, Mesh *const &external_mesh_pt, int &n_cont_inter_values)
Helper function to get the required nodal information from an external haloed node so that a fully-fu...
bool Use_bulk_element_as_external
Boolean to indicate when to use the bulk element as the external element. Defaults to false...
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().
double & s_macro_ll(const unsigned &i)
Access fct to the i-th coordinate of the element's "lower left" vertex in the associated MacroElement...
Definition: Qelements.h:223
double ref_value(const unsigned &i)
Return i-th reference value involved in default (usually first) update function.
A Class for nodes that deform elastically (i.e. position is an unknown in the problem). The idea is that the Eulerian positions are stored in a Data object and the Lagrangian coordinates are stored in addition. The pointer that addresses the Eulerian positions is set to the pointer to Value in the Data object. Hence, SolidNode uses knowledge of the internal structure of Data and must be a friend of the Data class. In order to allow a mesh to deform via an elastic-style equation in deforming-domain problems, the positions are stored separately from the values, so that elastic problems may be combined with any other type of problem.
Definition: nodes.h:1568
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size (broken virtual) ...
Definition: elements.h:1795
bool Setup_multi_domain_for_projection
Set up multi-domain for projection.
unsigned ngeom_object_list_pt()
Return number of geometric objects associated with AlgebraicMesh.
Vector< int > Proc_id_plus_one_of_external_element
Proc_id_plus_one_of_external_element[i] contains the processor id (plus one) of the processor on whic...
Definition: multi_domain.cc:97
void resize(const unsigned &n_value)
Resize the number of equations.
Definition: nodes.cc:2089
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
Vector< Vector< unsigned > > External_element_located
Lookup scheme for whether a local element's integration point has had an external element assigned to...
Definition: multi_domain.cc:74
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point...
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
unsigned ngeom_object(const int &id)
Number of geometric objects involved in id-th update function.
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4470
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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 set_node_update_info(FiniteElement *node_update_element_pt, const Vector< double > &s_in_node_update_element, const Vector< GeomObject * > &geom_object_pt)
Set node update information for node: Pass the pointer to the element that performs the update operat...
Vector< double > Zeta_coords_for_further_away_comparison
Vector of zeta coordinates that we're currently trying to locate; used in sorting of bin entries in f...
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
unsigned nlagrangian_type() const
Number of types of Lagrangian coordinates used to interpolate the Lagrangian coordinates within the e...
Definition: nodes.h:1737
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
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:806
double Z_min
Minimum coordinate in third dimension.
std::map< unsigned, unsigned > *& index_of_first_value_assigned_by_face_element_pt()
Return pointer to the map giving the index of the first face element value.
Definition: nodes.h:1909
Vector< double > Flat_packed_zetas_not_found_locally
Vector of flat-packed zeta coordinates for which the external element could not be found during curre...
Definition: multi_domain.cc:80
unsigned N_spiral_chunk
Number of spirals to be searched in one go.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
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().
bool is_halo() const
Is this element a halo?
Definition: elements.h:1141
Vector< unsigned > Located_element_status
Vector to indicate (to another processor) whether a located element (that will have to represented as...
unsigned nlagrangian() const
Return number of lagrangian coordinates.
Definition: nodes.h:1733
Vector< GeomObject * > geom_object_vector_pt()
Access function to the vector of GeomObject.
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 add_external_halo_node_helper(Node *&new_nod_pt, Mesh *const &external_mesh_pt, unsigned &loc_p, unsigned &node_index, FiniteElement *const &new_el_pt, int &n_cont_inter_values, Problem *problem_pt)
Helper functiono to add external halo node that is not a master.
void locate_zeta_for_local_coordinates(Mesh *const &mesh_pt, Mesh *const &external_mesh_pt, MeshAsGeomObject *&mesh_geom_obj_pt, const unsigned &interaction_index)
locate zeta for current set of "local" coordinates
unsigned Counter_for_flat_packed_doubles
Counter used when processing vector of flat-packed doubles – this is really "private" data...
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:57