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$
7 //LIC//
8 //LIC// $LastChangedDate$
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"
38 #include "macro_element_node_update_element.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(...)
53  std::ofstream Doc_boundary_coordinate_file;
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.
62  bool Accept_failed_locate_zeta_in_setup_multi_domain_interaction=false;
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.
74  Vector<Vector<unsigned> > External_element_located;
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.
80  Vector<double> Flat_packed_zetas_not_found_locally;
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.
88  Vector<double> Received_flat_packed_zetas_to_be_found;
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.
97  Vector<int> Proc_id_plus_one_of_external_element;
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).
104  Vector<unsigned> Located_element_status;
105 
106  /// \short Vector of flat-packed local coordinates for zeta tuples
107  /// that have been located
108  Vector<double> Flat_packed_located_coordinates;
109 
110  /// \short Vector of flat-packed doubles to be communicated with
111  /// other processors
112  Vector<double> Flat_packed_doubles;
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
118  unsigned Counter_for_flat_packed_doubles;
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
124  Vector<unsigned> Flat_packed_unsigneds;
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...)
131  Vector<std::string> Flat_packed_unsigneds_string;
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
140  unsigned Counter_for_flat_packed_unsigneds;
141 
142  // Other parameters
143  //-----------------
144 
145  /// \short Boolean to indicate when to use the bulk element as the
146  /// external element. Defaults to false, you must have set up FaceElements
147  /// properly first in order for it to work
148  bool Use_bulk_element_as_external=false;
149 
150  /// \short Boolean to indicate if we're allowed to use halo elements
151  /// as external elements. Can drastically reduce the number of
152  /// external halo elements -- currently not aware of any problems
153  /// therefore set to true by default but retention
154  /// of this flag allows easy return to previous implementation.
155  bool Allow_use_of_halo_elements_as_external_elements=true;
156 
157  /// \short Indicate whether we are allowed to use halo elements as
158  /// external elements for projection, possibly only required in
159  /// parallel unstructured mesh generation during the projection
160  /// stage. Default set to true
161  bool Allow_use_of_halo_elements_as_external_elements_for_projection=true;
162 
163  /// \short Boolean to indicate whether to doc timings or not.
164  bool Doc_timings=false;
165 
166  /// \short Boolean to indicate whether to output basic info during
167  /// setup_multi_domain_interaction() routines
168  bool Doc_stats=false;
169 
170  /// \short Boolean to indicate whether to output further info during
171  /// setup_multi_domain_interaction() routines
172  bool Doc_full_stats=false;
173 
174 #ifdef OOMPH_HAS_MPI
175 
176  // Functions for location method in multi-domain problems
177 
178 
179 
180 //========================================================================
181 /// Send the zeta coordinates from the current process to
182 /// the next process; receive from the previous process
183 //========================================================================
184  void send_and_receive_missing_zetas(Problem* problem_pt)
185  {
186  // MPI info
187  MPI_Status status;
188  MPI_Request request;
189 
190  // Storage for number of processors, current process and communicator
191  int n_proc=problem_pt->communicator_pt()->nproc();
192  int my_rank=problem_pt->communicator_pt()->my_rank();
193  OomphCommunicator* comm_pt=problem_pt->communicator_pt();
194 
195  // Work out processors to send and receive from
196  int send_to_proc=my_rank+1;
197  int recv_from_proc=my_rank-1;
198  if (send_to_proc==n_proc) { send_to_proc=0; }
199  if (recv_from_proc<0) { recv_from_proc=n_proc-1; }
200 
201  // Send the number of flat-packed zetas that we couldn't find
202  // locally to the next processor
203  int n_missing_local_zetas=Flat_packed_zetas_not_found_locally.size();
204  MPI_Isend(&n_missing_local_zetas,1,MPI_INT,send_to_proc,4,
205  comm_pt->mpi_comm(),&request);
206 
207  // Receive the number of flat-packed zetas that couldn't be found
208  // on the "previous" processor
209  int count_zetas=0;
210  MPI_Recv(&count_zetas,1,MPI_INT,recv_from_proc,
211  4,comm_pt->mpi_comm(),&status);
212 
213  MPI_Wait(&request,MPI_STATUS_IGNORE);
214 
215  // Send the vector of flat-packed zetas that we couldn't find
216  // locally to the next processor
217  if (n_missing_local_zetas!=0)
218  {
219  MPI_Isend(&Flat_packed_zetas_not_found_locally[0],
220  n_missing_local_zetas,MPI_DOUBLE,
221  send_to_proc,5,comm_pt->mpi_comm(),&request);
222  }
223 
224  // Receive the vector of flat-packed zetas that couldn't be found
225  // on the "previous" processor
226  if (count_zetas!=0)
227  {
228  Received_flat_packed_zetas_to_be_found.resize(count_zetas);
229  MPI_Recv(&Received_flat_packed_zetas_to_be_found[0],
230  count_zetas,MPI_DOUBLE,recv_from_proc,5,
231  comm_pt->mpi_comm(),&status);
232  MPI_Wait(&request,MPI_STATUS_IGNORE);
233  }
234  else
235  {
236  Received_flat_packed_zetas_to_be_found.resize(0);
237  }
238 
239  // Now we should have the Zeta arrays set up correctly
240  // for the next round of locations
241  }
242 
243 //========start of send_and_receive_located_info==========================
244 /// Send location information from current process; Received location
245 /// information from (current process + iproc) modulo (nproc)
246 //========================================================================
247  void send_and_receive_located_info
248  (int& iproc, Mesh* const &external_mesh_pt, Problem* problem_pt)
249  {
250  // Set MPI info
251  MPI_Status status;
252  MPI_Request request;
253 
254  // Storage for number of processors, current process and communicator
255  OomphCommunicator* comm_pt=problem_pt->communicator_pt();
256  int n_proc=comm_pt->nproc();
257  int my_rank=comm_pt->my_rank();
258 
259  // Prepare vectors to receive information
260  Vector<double> received_double_values;
261  Vector<unsigned> received_unsigned_values;
262  Vector<double> received_located_coord;
263  Vector<int> received_proc_id_plus_one_of_external_element;
264  Vector<unsigned> received_located_element_status;
265 
266  // Communicate the located information back to the original process
267  int orig_send_proc=my_rank-iproc;
268  if (my_rank<iproc) { orig_send_proc=n_proc+orig_send_proc; }
269  int orig_recv_proc=my_rank+iproc;
270  if ((my_rank+iproc)>=n_proc) { orig_recv_proc=orig_recv_proc-n_proc; }
271 
272  // Send the double values associated with external halos
273  //------------------------------------------------------
274  unsigned send_count_double_values=Flat_packed_doubles.size();
275  MPI_Isend(&send_count_double_values,1,MPI_UNSIGNED,
276  orig_send_proc,1,comm_pt->mpi_comm(),&request);
277  int receive_count_double_values=0;
278  MPI_Recv(&receive_count_double_values,1,MPI_INT,
279  orig_recv_proc,1,comm_pt->mpi_comm(),&status);
280  MPI_Wait(&request,MPI_STATUS_IGNORE);
281 
282  if (send_count_double_values!=0)
283  {
284  MPI_Isend(&Flat_packed_doubles[0],send_count_double_values,MPI_DOUBLE,
285  orig_send_proc,2,comm_pt->mpi_comm(),&request);
286  }
287  if (receive_count_double_values!=0)
288  {
289  received_double_values.resize(receive_count_double_values);
290  MPI_Recv(&received_double_values[0],receive_count_double_values,
291  MPI_DOUBLE,orig_recv_proc,2,comm_pt->mpi_comm(),&status);
292  }
293  if (send_count_double_values!=0)
294  {
295  MPI_Wait(&request,MPI_STATUS_IGNORE);
296  }
297 
298  // Now send unsigned values associated with external halos
299  //---------------------------------------------------------
300  unsigned send_count_unsigned_values=Flat_packed_unsigneds.size();
301  MPI_Isend(&send_count_unsigned_values,1,MPI_UNSIGNED,
302  orig_send_proc,14,comm_pt->mpi_comm(),&request);
303 
304  int receive_count_unsigned_values=0;
305  MPI_Recv(&receive_count_unsigned_values,1,MPI_INT,orig_recv_proc,14,
306  comm_pt->mpi_comm(),&status);
307 
308  MPI_Wait(&request,MPI_STATUS_IGNORE);
309 
310  if (send_count_unsigned_values!=0)
311  {
312  MPI_Isend(&Flat_packed_unsigneds[0],send_count_unsigned_values,MPI_UNSIGNED,
313  orig_send_proc,15,comm_pt->mpi_comm(),&request);
314 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
315  for (unsigned i=0;i<send_count_unsigned_values;i++)
316  {
317  oomph_info << "Sent:" << i << " to orig_proc:" << orig_send_proc
318  << " " << Flat_packed_unsigneds_string[i]
319  << ": " << Flat_packed_unsigneds[i] << std::endl;
320  }
321 #endif
322  }
323  if (receive_count_unsigned_values!=0)
324  {
325  received_unsigned_values.resize(receive_count_unsigned_values);
326  MPI_Recv(&received_unsigned_values[0],receive_count_unsigned_values,
327  MPI_UNSIGNED,orig_recv_proc,15,comm_pt->mpi_comm(),&status);
328  }
329 
330  if (send_count_unsigned_values!=0)
331  {
332  MPI_Wait(&request,MPI_STATUS_IGNORE);
333  }
334 
335  // Send and receive the Located_element_status
336  //--------------------------------------------
337  int send_count=Received_flat_packed_zetas_to_be_found.size()/Dim;
338  MPI_Isend(&send_count,1,MPI_INT,orig_send_proc,
339  20,comm_pt->mpi_comm(),&request);
340  int receive_count=0;
341  MPI_Recv(&receive_count,1,MPI_INT,orig_recv_proc,
342  20,comm_pt->mpi_comm(),&status);
343  MPI_Wait(&request,MPI_STATUS_IGNORE);
344 
345  if (send_count!=0)
346  {
347  MPI_Isend(&Located_element_status[0],send_count,MPI_UNSIGNED,
348  orig_send_proc,3,comm_pt->mpi_comm(),&request);
349  }
350 
351  if (receive_count!=0)
352  {
353  received_located_element_status.resize(receive_count);
354  MPI_Recv(&received_located_element_status[0],receive_count,
355  MPI_UNSIGNED,orig_recv_proc,3,
356  comm_pt->mpi_comm(),&status);
357  }
358  if (send_count!=0)
359  {
360  MPI_Wait(&request,MPI_STATUS_IGNORE);
361  }
362 
363  // Send and receive Proc_id_plus_one_of_external_element array
364  //------------------------------------------------------------
365  if (send_count!=0)
366  {
367  MPI_Isend(&Proc_id_plus_one_of_external_element[0],send_count,MPI_INT,
368  orig_send_proc,13,comm_pt->mpi_comm(),&request);
369  }
370  if (receive_count!=0)
371  {
372  received_proc_id_plus_one_of_external_element.resize(receive_count);
373  MPI_Recv(&received_proc_id_plus_one_of_external_element[0],
374  receive_count,MPI_INT,orig_recv_proc,13,
375  comm_pt->mpi_comm(),&status);
376  }
377  if (send_count!=0)
378  {
379  MPI_Wait(&request,MPI_STATUS_IGNORE);
380  }
381 
382 
383  // And finally the Flat_packed_located_coordinates array
384  //------------------------------------------------------
385  unsigned send_count_located_coord=Flat_packed_located_coordinates.size();
386  MPI_Isend(&send_count_located_coord,1,MPI_UNSIGNED,
387  orig_send_proc,4,comm_pt->mpi_comm(),&request);
388  unsigned receive_count_located_coord=0;
389  MPI_Recv(&receive_count_located_coord,1,MPI_UNSIGNED,orig_recv_proc,4,
390  comm_pt->mpi_comm(),&status);
391  MPI_Wait(&request,MPI_STATUS_IGNORE);
392 
393  if (send_count_located_coord!=0)
394  {
395  MPI_Isend(&Flat_packed_located_coordinates[0],
396  send_count_located_coord,MPI_DOUBLE,
397  orig_send_proc,5,comm_pt->mpi_comm(),&request);
398  }
399  if (receive_count_located_coord!=0)
400  {
401  received_located_coord.resize(receive_count_located_coord);
402  MPI_Recv(&received_located_coord[0],receive_count_located_coord,MPI_DOUBLE,
403  orig_recv_proc,5,comm_pt->mpi_comm(),&status);
404  }
405  if (send_count_located_coord!=0)
406  {
407  MPI_Wait(&request,MPI_STATUS_IGNORE);
408  }
409 
410  // Copy across into original containers -- these can now
411  //------------------------------------------------------
412  // be processed by create_external_halo_elements() to generate
413  //------------------------------------------------------------
414  // external halo elements
415  //------------------------
416  Flat_packed_doubles.resize(receive_count_double_values);
417  for (int ii=0;ii<receive_count_double_values;ii++)
418  {
419  Flat_packed_doubles[ii]=received_double_values[ii];
420  }
421  Flat_packed_unsigneds.resize(receive_count_unsigned_values);
422  for (int ii=0;ii<receive_count_unsigned_values;ii++)
423  {
424  Flat_packed_unsigneds[ii]=received_unsigned_values[ii];
425  }
426  Proc_id_plus_one_of_external_element.resize(receive_count);
427  Located_element_status.resize(receive_count);
428  for (int ii=0;ii<receive_count;ii++)
429  {
430  Proc_id_plus_one_of_external_element[ii]=
431  received_proc_id_plus_one_of_external_element[ii];
432  Located_element_status[ii]=received_located_element_status[ii];
433  }
434  Flat_packed_located_coordinates.resize(receive_count_located_coord);
435  for (int ii=0;ii<int(receive_count_located_coord);ii++)
436  {
437  Flat_packed_located_coordinates[ii]=received_located_coord[ii];
438  }
439 
440  }
441 
442 
443 //========start of add_external_haloed_node_to_storage====================
444 /// Helper function to add external haloed nodes, including any masters
445 //========================================================================
446  void add_external_haloed_node_to_storage(int& iproc, Node* nod_pt,
447  Problem* problem_pt,
448  Mesh* const &external_mesh_pt,
449  int& n_cont_inter_values)
450  {
451  // Add the node if required
452  add_external_haloed_node_helper(iproc,nod_pt,problem_pt,external_mesh_pt,
453  n_cont_inter_values);
454 
455  // Recursively add any master nodes (and their master nodes etc)
456  recursively_add_masters_of_external_haloed_node(iproc, nod_pt,
457  problem_pt,
458  external_mesh_pt,
459  n_cont_inter_values);
460  }
461 
462 
463 
464  //========================================================================
465  ///Recursively add any master nodes (and their master nodes etc) of
466  /// external nodes
467  //========================================================================
468  void recursively_add_masters_of_external_haloed_node(
469  int& iproc, Node* nod_pt,
470  Problem* problem_pt,
471  Mesh* const &external_mesh_pt,
472  int& n_cont_inter_values)
473  {
474 
475  // Loop over continuously interpolated values and add masters
476  for (int i_cont=-1;i_cont<n_cont_inter_values;i_cont++)
477  {
478  if (nod_pt->is_hanging(i_cont))
479  {
480  // Indicate that this node is a hanging node so the other
481  // process knows to create HangInfo and masters, etc.
482  Flat_packed_unsigneds.push_back(1);
483 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
484  Flat_packed_unsigneds_string.push_back("Is hanging");
485 #endif
486  // If this is a hanging node then add all its masters as
487  // external halo nodes if they have not yet been added
488  HangInfo* hang_pt=nod_pt->hanging_pt(i_cont);
489  // Loop over masters
490  unsigned n_master=hang_pt->nmaster();
491 
492  // Indicate number of master nodes to add on other process
493  Flat_packed_unsigneds.push_back(n_master);
494 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
495  Flat_packed_unsigneds_string.push_back("nmaster");
496 #endif
497  for (unsigned m=0;m<n_master;m++)
498  {
499  Node* master_nod_pt=hang_pt->master_node_pt(m);
500 
501  // Call the helper function for master nodes
502  add_external_haloed_master_node_helper(iproc,master_nod_pt,
503  problem_pt,
504  external_mesh_pt,
505  n_cont_inter_values);
506 
507  // Indicate the weight of this master
508  Flat_packed_doubles.push_back(hang_pt->master_weight(m));
509 
510  // Recursively add any master nodes (and their master nodes etc)
511  recursively_add_masters_of_external_haloed_node(iproc, master_nod_pt,
512  problem_pt,
513  external_mesh_pt,
514  n_cont_inter_values);
515  }
516  }
517  else
518  {
519  // Indicate that it's not a hanging node in this variable
520  Flat_packed_unsigneds.push_back(0);
521 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
522  Flat_packed_unsigneds_string.push_back("Not hanging");
523 #endif
524  }
525  } // end loop over continously interpolated values
526 
527  }
528 
529 //==========start of add_external_haloed_node_helper======================
530 /// Helper to add external haloed node that is not a master
531 //========================================================================
532  void add_external_haloed_node_helper(int& iproc, Node* nod_pt,
533  Problem* problem_pt,
534  Mesh* const &external_mesh_pt,
535  int& n_cont_inter_values)
536  {
537  // Attempt to add this node as an external haloed node
538  unsigned n_ext_haloed_nod=external_mesh_pt->nexternal_haloed_node(iproc);
539  unsigned external_haloed_node_index=
540  external_mesh_pt->add_external_haloed_node_pt(iproc,nod_pt);
541 
542  // If it was added then the new index should match the size of the storage
543  if (external_haloed_node_index==n_ext_haloed_nod)
544  {
545  Flat_packed_unsigneds.push_back(1);
546 
547 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
548  std::stringstream junk;
549  junk << "Node needs to be constructed [size="
550  << Flat_packed_unsigneds.size() << "]; last entry: "
551  << Flat_packed_unsigneds[Flat_packed_unsigneds.size()-1];
552  Flat_packed_unsigneds_string.push_back(junk.str());
553 #endif
554 
555  // This helper function gets all the required information for the
556  // specified node and stores it into MPI-sendable information
557  // so that a halo copy can be made on the receiving process
558  get_required_nodal_information_helper(iproc,nod_pt,
559  problem_pt,external_mesh_pt,
560  n_cont_inter_values);
561  }
562  else // It was already added
563  {
564  Flat_packed_unsigneds.push_back(0);
565 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
566  std::stringstream junk;
567  junk << "Node was already added [size="
568  << Flat_packed_unsigneds.size() << "]; last entry: "
569  << Flat_packed_unsigneds[Flat_packed_unsigneds.size()-1];
570 
571  Flat_packed_unsigneds_string.push_back(junk.str());
572 #endif
573 
574  // This node is already an external haloed node, so tell
575  // the other process its index in the equivalent external halo storage
576  Flat_packed_unsigneds.push_back(external_haloed_node_index);
577 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
578  Flat_packed_unsigneds_string.push_back("external haloed node index");
579 #endif
580  }
581 
582 
583 
584  }
585 
586 
587 //==========start of add_external_haloed_master_node_helper===============
588 /// Helper function to add external haloed node that is a master
589 //========================================================================
590  void add_external_haloed_master_node_helper(int& iproc, Node* master_nod_pt,
591  Problem* problem_pt,
592  Mesh* const &external_mesh_pt,
593  int& n_cont_inter_values)
594  {
595  // Attempt to add node as an external haloed node
596  unsigned n_ext_haloed_nod=external_mesh_pt->nexternal_haloed_node(iproc);
597  unsigned external_haloed_node_index;
598  external_haloed_node_index=
599  external_mesh_pt->add_external_haloed_node_pt(iproc,master_nod_pt);
600 
601  // If it was added the returned index is the same as current storage size
602  if (external_haloed_node_index==n_ext_haloed_nod)
603  {
604  // Indicate that this node needs to be constructed on
605  // the other process
606  Flat_packed_unsigneds.push_back(1);
607 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
608  Flat_packed_unsigneds_string.push_back("Node needs to be constructed[2]");
609 #endif
610 
611  // This gets all the required information for the specified
612  // master node and stores it into MPI-sendable information
613  // so that a halo copy can be made on the receiving process
614  get_required_master_nodal_information_helper(iproc,master_nod_pt,
615  problem_pt,
616  external_mesh_pt,
617  n_cont_inter_values);
618  }
619  else // It was already added
620  {
621  Flat_packed_unsigneds.push_back(0);
622 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
623  Flat_packed_unsigneds_string.push_back("Node was already added[2]");
624 #endif
625 
626  // This node is already an external haloed node, so tell
627  // the other process its index in the equivalent external halo storage
628  Flat_packed_unsigneds.push_back(external_haloed_node_index);
629 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
630  Flat_packed_unsigneds_string.push_back("external haloed node index[2]");
631 #endif
632  }
633  }
634 
635 
636 
637 
638 //========start of get_required_nodal_information_helper==================
639 /// Helper function to get the required nodal information from an
640 /// external haloed node so that a fully-functional external halo
641 /// node (and therefore element) can be created on the receiving process
642 //========================================================================
643  void get_required_nodal_information_helper(int& iproc, Node* nod_pt,
644  Problem* problem_pt,
645  Mesh* const &external_mesh_pt,
646  int& n_cont_inter_values)
647  {
648  // Tell the halo copy of this node how many values there are
649  // [NB this may be different for nodes within the same element, e.g.
650  // when using Lagrange multipliers]
651  unsigned n_val=nod_pt->nvalue();
652  Flat_packed_unsigneds.push_back(n_val);
653 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
654  Flat_packed_unsigneds_string.push_back("Number of values");
655 #endif
656 
657  unsigned n_dim=nod_pt->ndim();
658  TimeStepper* time_stepper_pt=nod_pt->time_stepper_pt();
659 
660  // Find the timestepper in the list of problem timesteppers
661  bool found_timestepper=false;
662  unsigned time_stepper_index;
663  unsigned n_time_steppers=problem_pt->ntime_stepper();
664  for (unsigned i=0;i<n_time_steppers;i++)
665  {
666  if (time_stepper_pt==problem_pt->time_stepper_pt(i))
667  {
668  // Indicate the timestepper's index
669  found_timestepper=true;
670  time_stepper_index=i;
671  break;
672  }
673  }
674 
675  if (found_timestepper)
676  {
677  Flat_packed_unsigneds.push_back(1);
678 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
679  Flat_packed_unsigneds_string.push_back("Found timestepper");
680 #endif
681  Flat_packed_unsigneds.push_back(time_stepper_index);
682 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
683  Flat_packed_unsigneds_string.push_back("Timestepper index");
684 #endif
685  }
686  else
687  {
688  Flat_packed_unsigneds.push_back(0);
689 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
690  Flat_packed_unsigneds_string.push_back("Not found timestepper");
691 #endif
692  }
693 
694  // Default number of previous values to 1
695  unsigned n_prev=1;
696  if (time_stepper_pt!=0)
697  {
698  // Add number of history values to n_prev
699  n_prev=time_stepper_pt->ntstorage();
700  }
701 
702  // Is the node on any boundaries?
703  if (nod_pt->is_on_boundary())
704  {
705  Flat_packed_unsigneds.push_back(1);
706 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
707  Flat_packed_unsigneds_string.push_back("Node is on boundary");
708 #endif
709 
710  // Loop over the boundaries of the external mesh
711  Vector<unsigned> boundaries;
712  unsigned n_bnd=external_mesh_pt->nboundary();
713  for (unsigned i_bnd=0;i_bnd<n_bnd;i_bnd++)
714  {
715  // Which boundaries (could be more than one) is it on?
716  if (nod_pt->is_on_boundary(i_bnd))
717  {
718  boundaries.push_back(i_bnd);
719  }
720  }
721  unsigned nb=boundaries.size();
722  Flat_packed_unsigneds.push_back(nb);
723 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
724  std::stringstream junk;
725  junk << "Node is on "<< nb << " boundaries";
726  Flat_packed_unsigneds_string.push_back(junk.str());
727 #endif
728  for (unsigned i=0;i<nb;i++)
729  {
730  Flat_packed_unsigneds.push_back(boundaries[i]);
731 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
732  std::stringstream junk;
733  junk << "Node is on boundary " << boundaries[i] << " of " << n_bnd;
734  Flat_packed_unsigneds_string.push_back(junk.str());
735 #endif
736  }
737 
738  // Get pointer to the map of indices associated with
739  // additional values created by face elements
740  BoundaryNodeBase* bnod_pt=
741  dynamic_cast<BoundaryNodeBase*>(nod_pt);
742 #ifdef PARANOID
743  if (bnod_pt==0)
744  {
745  throw OomphLibError(
746  "Failed to cast new node to boundary node\n",
747  OOMPH_CURRENT_FUNCTION,
748  OOMPH_EXCEPTION_LOCATION);
749  }
750 #endif
751  std::map<unsigned, unsigned>* map_pt=
752  bnod_pt->index_of_first_value_assigned_by_face_element_pt();
753 
754  // No additional values created
755  if (map_pt==0)
756  {
757  Flat_packed_unsigneds.push_back(0);
758 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
759  std::stringstream junk;
760  Flat_packed_unsigneds_string.push_back("No additional values were created by face element");
761 #endif
762  }
763  // Created additional values
764  else
765  {
766  // How many?
767  Flat_packed_unsigneds.push_back(map_pt->size());
768 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
769  std::stringstream junk;
770  junk << "Map size " << map_pt->size() << n_bnd;
771  Flat_packed_unsigneds_string.push_back(junk.str());
772 #endif
773  // Loop over entries in map and add to send data
774  for (std::map<unsigned, unsigned>::iterator p=
775  map_pt->begin();
776  p!=map_pt->end();p++)
777  {
778  Flat_packed_unsigneds.push_back((*p).first);
779 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
780  std::stringstream junk;
781  Flat_packed_unsigneds_string.push_back("Key of map entry");
782 #endif
783  Flat_packed_unsigneds.push_back((*p).second);
784 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
785  Flat_packed_unsigneds_string.push_back("Value of map entry");
786 #endif
787  }
788  }
789  }
790  else
791  {
792  // Not on any boundary
793  Flat_packed_unsigneds.push_back(0);
794 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
795  Flat_packed_unsigneds_string.push_back("Node is not on any boundary");
796 #endif
797  }
798 
799  // Is the Node algebraic? If so, send its ref values and
800  // an indication of its geometric objects if they are stored
801  // in the algebraic mesh
802  AlgebraicNode* alg_nod_pt=dynamic_cast<AlgebraicNode*>(nod_pt);
803  if (alg_nod_pt!=0)
804  {
805  // The external mesh should be algebraic
806  AlgebraicMesh* alg_mesh_pt=dynamic_cast<AlgebraicMesh*>
807  (external_mesh_pt);
808 
809  // Get default node update function ID
810  unsigned update_id=alg_nod_pt->node_update_fct_id();
811  Flat_packed_unsigneds.push_back(update_id);
812 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
813  Flat_packed_unsigneds_string.push_back("Alg Node update id");
814 #endif
815 
816  // Get reference values at default...
817  unsigned n_ref_val=alg_nod_pt->nref_value();
818  Flat_packed_unsigneds.push_back(n_ref_val);
819 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
820  Flat_packed_unsigneds_string.push_back("Alg Node n ref values");
821 #endif
822  for (unsigned i_ref_val=0;i_ref_val<n_ref_val;i_ref_val++)
823  {
824  Flat_packed_doubles.push_back(alg_nod_pt->ref_value(i_ref_val));
825  }
826 
827  // Access geometric objects at default...
828  unsigned n_geom_obj=alg_nod_pt->ngeom_object();
829  Flat_packed_unsigneds.push_back(n_geom_obj);
830 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
831  Flat_packed_unsigneds_string.push_back("Alg Node n geom objects");
832 #endif
833  for (unsigned i_geom=0;i_geom<n_geom_obj;i_geom++)
834  {
835  GeomObject* geom_obj_pt=alg_nod_pt->geom_object_pt(i_geom);
836 
837  // Check this against the stored geometric objects in mesh
838  unsigned n_geom_list=alg_mesh_pt->ngeom_object_list_pt();
839 
840  // Default found index to zero
841  unsigned found_geom_object=0;
842  for (unsigned i_list=0;i_list<n_geom_list;i_list++)
843  {
844  if (geom_obj_pt==alg_mesh_pt->geom_object_list_pt(i_list))
845  {
846  found_geom_object=i_list;
847  }
848  }
849  Flat_packed_unsigneds.push_back(found_geom_object);
850 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
851  Flat_packed_unsigneds_string.push_back("Found geom object");
852 #endif
853  }
854  }
855 
856  // If it is a MacroElementNodeUpdateNode, everything has been
857  // dealt with by the new element already
858 
859  // Is it a SolidNode?
860  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
861  if (solid_nod_pt!=0)
862  {
863  unsigned n_solid_val=solid_nod_pt->variable_position_pt()->nvalue();
864  for (unsigned i_val=0;i_val<n_solid_val;i_val++)
865  {
866  for (unsigned t=0;t<n_prev;t++)
867  {
868  Flat_packed_doubles.push_back(solid_nod_pt->variable_position_pt()->
869  value(t,i_val));
870  }
871  }
872  }
873 
874  // Finally copy info required for all node types
875  for (unsigned i_val=0;i_val<n_val;i_val++)
876  {
877  for (unsigned t=0;t<n_prev;t++)
878  {
879  Flat_packed_doubles.push_back(nod_pt->value(t,i_val));
880  }
881  }
882 
883  // Now do positions
884  for (unsigned idim=0;idim<n_dim;idim++)
885  {
886  for (unsigned t=0;t<n_prev;t++)
887  {
888  Flat_packed_doubles.push_back(nod_pt->x(t,idim));
889  }
890  }
891  }
892 
893 //=========start of get_required_master_nodal_information_helper==========
894 /// Helper function to get the required master nodal information from an
895 /// external haloed master node so that a fully-functional external halo
896 /// master node (and possible element) can be created on the receiving process
897 //========================================================================
898  void get_required_master_nodal_information_helper
899  (int& iproc, Node* master_nod_pt, Problem* problem_pt,
900  Mesh* const &external_mesh_pt, int& n_cont_inter_values)
901  {
902  // Need to send over dimension, position type and number of values
903  Flat_packed_unsigneds.push_back(master_nod_pt->ndim());
904 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
905  Flat_packed_unsigneds_string.push_back("Master node ndim");
906 #endif
907  Flat_packed_unsigneds.push_back(master_nod_pt->nposition_type());
908 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
909  Flat_packed_unsigneds_string.push_back("Master node npos_type");
910 #endif
911  Flat_packed_unsigneds.push_back(master_nod_pt->nvalue());
912 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
913  Flat_packed_unsigneds_string.push_back("Master node nvalue");
914 #endif
915 
916  // If it's a solid node, also need to send lagrangian dim and type
917  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(master_nod_pt);
918  if (solid_nod_pt!=0)
919  {
920  Flat_packed_unsigneds.push_back(solid_nod_pt->nlagrangian());
921 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
922  Flat_packed_unsigneds_string.push_back("Master solid node nlagr");
923 #endif
924  Flat_packed_unsigneds.push_back(solid_nod_pt->nlagrangian_type());
925 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
926  Flat_packed_unsigneds_string.push_back("Master solid node nlagr_type");
927 #endif
928  }
929 
930  unsigned n_dim=master_nod_pt->ndim();
931  TimeStepper* time_stepper_pt=master_nod_pt->time_stepper_pt();
932 
933  // Find the timestepper in the list of problem timesteppers
934  bool found_timestepper=false;
935  unsigned time_stepper_index;
936  unsigned n_time_steppers=problem_pt->ntime_stepper();
937  for (unsigned i=0;i<n_time_steppers;i++)
938  {
939  if (time_stepper_pt==problem_pt->time_stepper_pt(i))
940  {
941  // Indicate the timestepper's index
942  // add 1 to the index so that 0 indicates no timestepper?
943  found_timestepper=true;
944  time_stepper_index=i;
945  break;
946  }
947  }
948 
949  if (found_timestepper)
950  {
951  Flat_packed_unsigneds.push_back(1);
952 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
953  Flat_packed_unsigneds_string.push_back("Master node Found timestepper");
954 #endif
955  Flat_packed_unsigneds.push_back(time_stepper_index);
956 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
957  Flat_packed_unsigneds_string.push_back("Master node Timestepper index");
958 #endif
959  }
960  else
961  {
962  Flat_packed_unsigneds.push_back(0);
963 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
964  Flat_packed_unsigneds_string.push_back("Master node Not found timestepper");
965 #endif
966  }
967 
968  // Default number of previous values to 1
969  unsigned n_prev=1;
970  if (time_stepper_pt!=0)
971  {
972  // Add number of history values to n_prev
973  n_prev=time_stepper_pt->ntstorage();
974  }
975 
976  // Is the node on any boundaries?
977  if (master_nod_pt->is_on_boundary())
978  {
979  Flat_packed_unsigneds.push_back(1);
980 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
981  Flat_packed_unsigneds_string.push_back("Master node is on boundary");
982 #endif
983  // Loop over the boundaries of the external mesh
984  Vector<unsigned> boundaries;
985  unsigned n_bnd=external_mesh_pt->nboundary();
986  for (unsigned i_bnd=0;i_bnd<n_bnd;i_bnd++)
987  {
988  // Which boundaries (could be more than one) is it on?
989  if (master_nod_pt->is_on_boundary(i_bnd))
990  {
991  boundaries.push_back(i_bnd);
992  }
993  }
994  unsigned nb=boundaries.size();
995  Flat_packed_unsigneds.push_back(nb);
996 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
997  std::stringstream junk;
998  junk << "Master node is on "<< nb << " boundaries";
999  Flat_packed_unsigneds_string.push_back(junk.str());
1000 #endif
1001  for (unsigned i=0;i<nb;i++)
1002  {
1003  Flat_packed_unsigneds.push_back(boundaries[i]);
1004 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1005  std::stringstream junk;
1006  junk << "Master noode is on boundary "
1007  << boundaries[i] << " of " << n_bnd;
1008  Flat_packed_unsigneds_string.push_back(junk.str());
1009 #endif
1010  }
1011 
1012  // Get pointer to the map of indices associated with
1013  // additional values created by face elements
1014  BoundaryNodeBase* bnod_pt=
1015  dynamic_cast<BoundaryNodeBase*>(master_nod_pt);
1016 #ifdef PARANOID
1017  if (bnod_pt==0)
1018  {
1019  throw OomphLibError(
1020  "Failed to cast new node to boundary node\n",
1021  OOMPH_CURRENT_FUNCTION,
1022  OOMPH_EXCEPTION_LOCATION);
1023  }
1024 #endif
1025  std::map<unsigned, unsigned>* map_pt=
1026  bnod_pt->index_of_first_value_assigned_by_face_element_pt();
1027 
1028  // No additional values created
1029  if (map_pt==0)
1030  {
1031  Flat_packed_unsigneds.push_back(0);
1032 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1033  std::stringstream junk;
1034  Flat_packed_unsigneds_string.push_back("No additional values were created by face element for this master node");
1035 #endif
1036  }
1037  // Created additional values
1038  else
1039  {
1040  // How many?
1041  Flat_packed_unsigneds.push_back(map_pt->size());
1042 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1043  std::stringstream junk;
1044  junk << "Map size for master node " << map_pt->size() << n_bnd;
1045  Flat_packed_unsigneds_string.push_back(junk.str());
1046 #endif
1047  // Loop over entries in map and add to send data
1048  for (std::map<unsigned, unsigned>::iterator p=
1049  map_pt->begin();
1050  p!=map_pt->end();p++)
1051  {
1052  Flat_packed_unsigneds.push_back((*p).first);
1053 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1054  std::stringstream junk;
1055  Flat_packed_unsigneds_string.push_back(
1056  "Key of map entry for master node");
1057 #endif
1058  Flat_packed_unsigneds.push_back((*p).second);
1059 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1060  Flat_packed_unsigneds_string.push_back(
1061  "Value of map entry for master node");
1062 #endif
1063  }
1064  }
1065  }
1066  else
1067  {
1068  // Not on any boundary
1069  Flat_packed_unsigneds.push_back(0);
1070 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1071  Flat_packed_unsigneds_string.push_back("Master node is not on any boundary");
1072 #endif
1073  }
1074 
1075  // Is the Node algebraic? If so, send its ref values and
1076  // an indication of its geometric objects if they are stored
1077  // in the algebraic mesh
1078  AlgebraicNode* alg_nod_pt=dynamic_cast<AlgebraicNode*>(master_nod_pt);
1079  if (alg_nod_pt!=0)
1080  {
1081  // The external mesh should be algebraic
1082  AlgebraicMesh* alg_mesh_pt=dynamic_cast<AlgebraicMesh*>
1083  (external_mesh_pt);
1084 
1085  // Get default node update function ID
1086  unsigned update_id=alg_nod_pt->node_update_fct_id();
1087  Flat_packed_unsigneds.push_back(update_id);
1088 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1089  Flat_packed_unsigneds_string.push_back("Master Alg Node update id");
1090 #endif
1091 
1092  // Get reference values at default...
1093  unsigned n_ref_val=alg_nod_pt->nref_value();
1094  Flat_packed_unsigneds.push_back(n_ref_val);
1095 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1096  Flat_packed_unsigneds_string.push_back("Master Alg Node n ref values");
1097 #endif
1098  for (unsigned i_ref_val=0;i_ref_val<n_ref_val;i_ref_val++)
1099  {
1100  Flat_packed_doubles.push_back(alg_nod_pt->ref_value(i_ref_val));
1101  }
1102 
1103  // Access geometric objects at default...
1104  unsigned n_geom_obj=alg_nod_pt->ngeom_object();
1105  Flat_packed_unsigneds.push_back(n_geom_obj);
1106 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1107  Flat_packed_unsigneds_string.push_back("Master Alg Node n geom objects");
1108 #endif
1109  for (unsigned i_geom=0;i_geom<n_geom_obj;i_geom++)
1110  {
1111  GeomObject* geom_obj_pt=alg_nod_pt->geom_object_pt(i_geom);
1112  // Check this against the stored geometric objects in mesh
1113  unsigned n_geom_list=alg_mesh_pt->ngeom_object_list_pt();
1114  // Default found index to zero
1115  unsigned found_geom_object=0;
1116  for (unsigned i_list=0;i_list<n_geom_list;i_list++)
1117  {
1118  if (geom_obj_pt==alg_mesh_pt->geom_object_list_pt(i_list))
1119  {
1120  found_geom_object=i_list;
1121  }
1122  }
1123  Flat_packed_unsigneds.push_back(found_geom_object);
1124 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1125  Flat_packed_unsigneds_string.push_back("Master node Found geom object");
1126 #endif
1127  }
1128  } // end AlgebraicNode check
1129 
1130  // Is it a MacroElementNodeUpdateNode?
1131  MacroElementNodeUpdateNode* macro_nod_pt=
1132  dynamic_cast<MacroElementNodeUpdateNode*>(master_nod_pt);
1133  if (macro_nod_pt!=0)
1134  {
1135  // Loop over current external haloed elements - has the element which
1136  // controls the node update for this node been added yet?
1137  GeneralisedElement* macro_node_update_el_pt=
1138  macro_nod_pt->node_update_element_pt();
1139 
1140  unsigned n_ext_haloed_el=external_mesh_pt->
1141  nexternal_haloed_element(iproc);
1142  unsigned external_haloed_el_index;
1143  external_haloed_el_index=external_mesh_pt->
1144  add_external_haloed_element_pt(iproc,macro_node_update_el_pt);
1145 
1146  // If it wasn't already added, we need to create a halo copy
1147  if (external_haloed_el_index==n_ext_haloed_el)
1148  {
1149  Flat_packed_unsigneds.push_back(1);
1150 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1151  Flat_packed_unsigneds_string.push_back("Master Node needs to be constructed");
1152 #endif
1153  //Cast to a finite elemnet
1154  FiniteElement* macro_node_update_finite_el_pt
1155  = dynamic_cast<FiniteElement*>(macro_node_update_el_pt);
1156 
1157  // We're using macro elements to update...
1158  MacroElementNodeUpdateMesh* macro_mesh_pt=
1159  dynamic_cast<MacroElementNodeUpdateMesh*>(external_mesh_pt);
1160  if (macro_mesh_pt!=0)
1161  {
1162  Flat_packed_unsigneds.push_back(1);
1163 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1164  Flat_packed_unsigneds_string.push_back("Mesh is macro element mesh");
1165 #endif
1166  // Need to send the macro element number in the mesh across
1167  MacroElement* macro_el_pt= macro_node_update_finite_el_pt
1168  ->macro_elem_pt();
1169  unsigned macro_el_num=macro_el_pt->macro_element_number();
1170  Flat_packed_unsigneds.push_back(macro_el_num);
1171 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1172  Flat_packed_unsigneds_string.push_back("Number of macro element");
1173 #endif
1174  // Also need to send
1175  // the lower left and upper right coordinates of the macro element
1176  QElementBase* q_el_pt=dynamic_cast<QElementBase*>
1177  (macro_node_update_el_pt);
1178  if (q_el_pt!=0)
1179  {
1180  // The macro element needs to be set first before
1181  // its lower left and upper right coordinates can be accessed
1182  // Now send the lower left and upper right coordinates
1183  unsigned el_dim=q_el_pt->dim();
1184  for (unsigned i_dim=0;i_dim<el_dim;i_dim++)
1185  {
1186  Flat_packed_doubles.push_back(q_el_pt->s_macro_ll(i_dim));
1187  Flat_packed_doubles.push_back(q_el_pt->s_macro_ur(i_dim));
1188  }
1189  }
1190  else // Throw an error
1191  {
1192  std::ostringstream error_stream;
1193  error_stream << "You are using a MacroElement node update\n"
1194  << "in a case with non-QElements. This has not\n"
1195  << "yet been implemented.\n";
1196  throw OomphLibError
1197  (error_stream.str(),
1198  OOMPH_CURRENT_FUNCTION,
1199  OOMPH_EXCEPTION_LOCATION);
1200  }
1201 
1202  }
1203  else // Not using macro elements for node update... umm, we're
1204  // already inside a loop over macro elements, so this
1205  // should never get here... an error should be thrown I suppose
1206  {
1207  Flat_packed_unsigneds.push_back(0);
1208 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1209  Flat_packed_unsigneds_string.push_back("Mesh is not a macro element mesh");
1210 #endif
1211  }
1212 
1213  // This element needs to be fully functioning on the other
1214  // process, so send all the information required to create it
1215  unsigned n_node=macro_node_update_finite_el_pt->nnode();
1216  for (unsigned j=0;j<n_node;j++)
1217  {
1218  Node* new_nod_pt=macro_node_update_finite_el_pt->node_pt(j);
1219  add_external_haloed_node_to_storage(iproc,new_nod_pt,
1220  problem_pt,
1221  external_mesh_pt,
1222  n_cont_inter_values);
1223  }
1224  }
1225  else // The external haloed element already exists
1226  {
1227  Flat_packed_unsigneds.push_back(0);
1228 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1229  Flat_packed_unsigneds_string.push_back("External haloed element already exists");
1230 #endif
1231  Flat_packed_unsigneds.push_back(external_haloed_el_index);
1232 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1233  Flat_packed_unsigneds_string.push_back("Index of existing external haloed element");
1234 #endif
1235  }
1236 
1237  } // end of MacroElementNodeUpdateNode check
1238 
1239  // Is it a SolidNode?
1240  if (solid_nod_pt!=0)
1241  {
1242  unsigned n_val=solid_nod_pt->variable_position_pt()->nvalue();
1243  for (unsigned i_val=0;i_val<n_val;i_val++)
1244  {
1245  for (unsigned t=0;t<n_prev;t++)
1246  {
1247  Flat_packed_doubles.push_back(solid_nod_pt->variable_position_pt()->
1248  value(t,i_val));
1249  }
1250  }
1251  }
1252 
1253  // Finally copy info required for all node types
1254 
1255  // Halo copy needs to know all the history values
1256  unsigned n_val=master_nod_pt->nvalue();
1257  for (unsigned i_val=0;i_val<n_val;i_val++)
1258  {
1259  for (unsigned t=0;t<n_prev;t++)
1260  {
1261  Flat_packed_doubles.push_back(master_nod_pt->value(t,i_val));
1262  }
1263  }
1264 
1265  // Now do positions
1266  for (unsigned idim=0;idim<n_dim;idim++)
1267  {
1268  for (unsigned t=0;t<n_prev;t++)
1269  {
1270  Flat_packed_doubles.push_back(master_nod_pt->x(t,idim));
1271  }
1272  }
1273 
1274  }
1275 
1276 
1277 
1278 
1279 
1280 
1281 
1282 //=======start of add_external_halo_node_helper===========================
1283 /// Helper functiono to add external halo node that is not a master
1284 //========================================================================
1285  void add_external_halo_node_helper
1286  (Node* &new_nod_pt, Mesh* const &external_mesh_pt, unsigned& loc_p,
1287  unsigned& node_index, FiniteElement* const &new_el_pt,
1288  int& n_cont_inter_values,
1289  Problem* problem_pt)
1290  {
1291  // Given the node and the external mesh, and received information
1292  // about them from process loc_p, construct them on the current process
1293 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1294  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1295  << " Bool: New node needs to be constructed "
1296  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1297  << std::endl;
1298 #endif
1299  if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]==1)
1300  {
1301  // Construct a new node based upon sent information
1302  construct_new_external_halo_node_helper(new_nod_pt,loc_p,
1303  node_index,new_el_pt,
1304  external_mesh_pt,problem_pt);
1305  }
1306  else
1307  {
1308 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1309  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1310  << " Index of existing external halo node "
1311  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1312  << std::endl;
1313 #endif
1314 
1315  // Copy node from received location
1316  new_nod_pt=external_mesh_pt->external_halo_node_pt
1317  (loc_p,Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]);
1318 
1319  new_el_pt->node_pt(node_index)=new_nod_pt;
1320 
1321 
1322  }
1323  }
1324 
1325 
1326 
1327 
1328 //========start of construct_new_external_halo_node_helper=================
1329 /// Helper function which constructs a new external halo node (on new element)
1330 /// with the required information sent from the haloed process
1331 //========================================================================
1332  void construct_new_external_halo_node_helper
1333  (Node* &new_nod_pt, unsigned& loc_p, unsigned& node_index,
1334  FiniteElement* const &new_el_pt,
1335  Mesh* const &external_mesh_pt, Problem* problem_pt)
1336  {
1337  // The first entry indicates the number of values at this new Node
1338  // (which may be different across the same element e.g. Lagrange multipliers)
1339 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1340  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1341  << " Number of values of external halo node "
1342  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1343  << std::endl;
1344 #endif
1345  unsigned n_val=Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
1346 
1347  // Null TimeStepper for now
1348  TimeStepper* time_stepper_pt=0;
1349  // Default number of previous values to 1
1350  unsigned n_prev=1;
1351 
1352  // The next entry in Flat_packed_unsigneds indicates
1353  // if a timestepper is required for this halo node
1354 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1355  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1356  << " Timestepper req'd for node "
1357  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1358  << std::endl;
1359 #endif
1360  if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]==1)
1361  {
1362 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1363  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1364  << " Index of timestepper "
1365  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1366  << std::endl;
1367 #endif
1368  // Index
1369  time_stepper_pt=problem_pt->time_stepper_pt
1370  (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]);
1371 
1372  // Check whether number of prev values is "sent" across
1373  n_prev=time_stepper_pt->ntstorage();
1374  }
1375 
1376  // If this node was on a boundary then it needs to
1377  // be on the same boundary here
1378 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1379  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1380  << " Is node on boundary? "
1381  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1382  << std::endl;
1383 #endif
1384  if (Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++]==1)
1385  {
1386  // Construct a new boundary node
1387  if (time_stepper_pt!=0)
1388  {
1389  new_nod_pt=new_el_pt->construct_boundary_node
1390  (node_index,time_stepper_pt);
1391  }
1392  else
1393  {
1394  new_nod_pt=new_el_pt->construct_boundary_node(node_index);
1395  }
1396 
1397  // How many boundaries does the node live on?
1398 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1399  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1400  << " Number of boundaries the node is on: "
1401  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1402  << std::endl;
1403 #endif
1404  unsigned nb=Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
1405  for (unsigned i=0;i<nb;i++)
1406  {
1407  // Boundary number
1408 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1409  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1410  << " Node is on boundary "
1411  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1412  << std::endl;
1413 #endif
1414  unsigned i_bnd=
1415  Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
1416  external_mesh_pt->add_boundary_node(i_bnd,new_nod_pt);
1417  }
1418 
1419  // Do we have additional values created by face elements?
1420 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1421  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1422  << " Number of additional values created by face element "
1423  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1424  << std::endl;
1425 #endif
1426  unsigned n_entry=Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
1427  if (n_entry>0)
1428  {
1429  // Create storage, if it doesn't already exist, for the map
1430  // that will contain the position of the first entry of
1431  // this face element's additional values,
1432  BoundaryNodeBase* bnew_nod_pt=
1433  dynamic_cast<BoundaryNodeBase*>(new_nod_pt);
1434 #ifdef PARANOID
1435  if (bnew_nod_pt==0)
1436  {
1437  throw OomphLibError(
1438  "Failed to cast new node to boundary node\n",
1439  OOMPH_CURRENT_FUNCTION,
1440  OOMPH_EXCEPTION_LOCATION);
1441  }
1442 #endif
1443  if(bnew_nod_pt->
1444  index_of_first_value_assigned_by_face_element_pt()==0)
1445  {
1446  bnew_nod_pt->
1447  index_of_first_value_assigned_by_face_element_pt()=
1448  new std::map<unsigned, unsigned>;
1449  }
1450 
1451  // Get pointer to the map of indices associated with
1452  // additional values created by face elements
1453  std::map<unsigned, unsigned>* map_pt=
1454  bnew_nod_pt->index_of_first_value_assigned_by_face_element_pt();
1455 
1456  // Loop over number of entries in map
1457  for (unsigned i=0;i<n_entry;i++)
1458  {
1459  // Read out pairs...
1460 
1461 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1462  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1463  << " Key of map entry"
1464  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1465  << std::endl;
1466 #endif
1467  unsigned first=Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
1468 
1469 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1470  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1471  << " Value of map entry"
1472  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1473  << std::endl;
1474 #endif
1475  unsigned second=Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds++];
1476 
1477  // ...and assign
1478  (*map_pt)[first]=second;
1479  }
1480  }
1481  }
1482  else
1483  {
1484  // Construct an ordinary (non-boundary) node
1485  if (time_stepper_pt!=0)
1486  {
1487  new_nod_pt=new_el_pt->construct_node
1488  (node_index,time_stepper_pt);
1489  }
1490  else
1491  {
1492  new_nod_pt=new_el_pt->construct_node(node_index);
1493  }
1494  }
1495 
1496  // Node constructed: add to external halo nodes
1497  external_mesh_pt->add_external_halo_node_pt(loc_p,new_nod_pt);
1498 
1499  // Is the new constructed node Algebraic?
1500  AlgebraicNode* new_alg_nod_pt=dynamic_cast<AlgebraicNode*>
1501  (new_nod_pt);
1502 
1503  // If it is algebraic, its node update functions will
1504  // not yet have been set up properly
1505  if (new_alg_nod_pt!=0)
1506  {
1507  // The AlgebraicMesh is the external mesh
1508  AlgebraicMesh* alg_mesh_pt=dynamic_cast<AlgebraicMesh*>
1509  (external_mesh_pt);
1510 
1511  /// The first entry of All_alg_nodal_info contains
1512  /// the default node update id
1513  /// e.g. for the quarter circle there are
1514  /// "Upper_left_box", "Lower right box" etc...
1515 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1516  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1517  << " Alg node update id "
1518  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1519  << std::endl;
1520 #endif
1521 
1522  unsigned update_id=Flat_packed_unsigneds
1523  [Counter_for_flat_packed_unsigneds++];
1524 
1525  Vector<double> ref_value;
1526 
1527  // The size of this vector is in the next entry
1528  // of All_alg_nodal_info
1529 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1530  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1531  << " Alg node # of ref values "
1532  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1533  << std::endl;
1534 #endif
1535  unsigned n_ref_val=Flat_packed_unsigneds
1536  [Counter_for_flat_packed_unsigneds++];
1537 
1538  // The reference values themselves are in
1539  // All_alg_ref_value
1540  ref_value.resize(n_ref_val);
1541  for (unsigned i_ref=0;i_ref<n_ref_val;i_ref++)
1542  {
1543  ref_value[i_ref]=Flat_packed_doubles
1544  [Counter_for_flat_packed_doubles++];
1545  }
1546 
1547  Vector<GeomObject*> geom_object_pt;
1548  /// again we need the size of this vector as it varies
1549  /// between meshes; we also need some indication
1550  /// as to which geometric object should be used...
1551 
1552  // The size of this vector is in the next entry
1553  // of All_alg_nodal_info
1554 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1555  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1556  << " Alg node # of geom objects "
1557  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1558  << std::endl;
1559 #endif
1560  unsigned n_geom_obj=Flat_packed_unsigneds
1561  [Counter_for_flat_packed_unsigneds++];
1562 
1563  // The remaining indices are in the rest of
1564  // All_alg_nodal_info
1565  geom_object_pt.resize(n_geom_obj);
1566  for (unsigned i_geom=0;i_geom<n_geom_obj;i_geom++)
1567  {
1568 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1569  oomph_info << "Rec:" << Counter_for_flat_packed_unsigneds
1570  << " Alg node: geom object index "
1571  << Flat_packed_unsigneds[Counter_for_flat_packed_unsigneds]
1572  << std::endl;
1573 #endif
1574  unsigned geom_index=Flat_packed_unsigneds
1575  [Counter_for_flat_packed_unsigneds++];
1576  // This index indicates which of the AlgebraicMesh's
1577  // stored geometric objects should be used
1578  // (0 is a null pointer; everything else should have
1579  // been filled in by the specific Mesh). If it
1580  // hasn't been filled in then the update_node_update
1581  // call should fix it
1582  geom_object_pt[i_geom]=alg_mesh_pt->
1583  geom_object_list_pt(geom_index);
1584  }
1585 
1586  /// For the received update_id, ref_value, geom_object
1587  /// call add_node_update_info
1588  new_alg_nod_pt->add_node_update_info
1589  (update_id,alg_mesh_pt,geom_object_pt,ref_value);
1590 
1591  /// Now call update_node_update
1592  alg_mesh_pt->update_node_update(new_alg_nod_pt);
1593  }
1594 
1595  // Is the node a MacroElementNodeUpdateNode?
1596  MacroElementNodeUpdateNode* macro_nod_pt=
1597  dynamic_cast<MacroElementNodeUpdateNode*>(new_nod_pt);
1598 
1599  if (macro_nod_pt!=0)
1600  {
1601  // Need to call set_node_update_info; this requires
1602  // a Vector<GeomObject*> (taken from the mesh)
1603  Vector<GeomObject*> geom_object_vector_pt;
1604 
1605  // Access the required geom objects from the
1606  // MacroElementNodeUpdateMesh
1607  MacroElementNodeUpdateMesh* macro_mesh_pt=
1608  dynamic_cast<MacroElementNodeUpdateMesh*>
1609  (external_mesh_pt);
1610  geom_object_vector_pt=
1611  macro_mesh_pt->geom_object_vector_pt();
1612 
1613  // Get local coordinate of node in new element
1614  Vector<double> s_in_macro_node_update_element;
1615  new_el_pt->local_coordinate_of_node
1616  (node_index,s_in_macro_node_update_element);
1617 
1618  // Set node update info for this node
1619  macro_nod_pt->set_node_update_info
1620  (new_el_pt,s_in_macro_node_update_element,
1621  geom_object_vector_pt);
1622  }
1623 
1624  // Is the new node a SolidNode?
1625  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(new_nod_pt);
1626  if (solid_nod_pt!=0)
1627  {
1628  unsigned n_solid_val=solid_nod_pt->variable_position_pt()->nvalue();
1629  for (unsigned i_val=0;i_val<n_solid_val;i_val++)
1630  {
1631  for (unsigned t=0;t<n_prev;t++)
1632  {
1633  solid_nod_pt->variable_position_pt()->
1634  set_value(t,i_val,
1635  Flat_packed_doubles[Counter_for_flat_packed_doubles++]);
1636  }
1637  }
1638  }
1639 
1640  // If there are additional values, resize the node
1641  unsigned n_new_val=new_nod_pt->nvalue();
1642  if (n_val>n_new_val)
1643  {
1644  new_nod_pt->resize(n_val);
1645  }
1646 
1647  // Get copied history values
1648  // unsigned n_val=new_nod_pt->nvalue();
1649  for (unsigned i_val=0;i_val<n_val;i_val++)
1650  {
1651  for (unsigned t=0;t<n_prev;t++)
1652  {
1653  new_nod_pt->set_value(t,i_val,Flat_packed_doubles
1654  [Counter_for_flat_packed_doubles++]);
1655  }
1656  }
1657 
1658  // Get copied history values for positions
1659  unsigned n_dim=new_nod_pt->ndim();
1660  for (unsigned idim=0;idim<n_dim;idim++)
1661  {
1662  for (unsigned t=0;t<n_prev;t++)
1663  {
1664  // Copy to coordinate
1665  new_nod_pt->x(t,idim)=Flat_packed_doubles
1666  [Counter_for_flat_packed_doubles++];
1667  }
1668  }
1669  }
1670 
1671 
1672 
1673  //=====================================================================
1674  /// Locate zeta for current set of missing coordinates; vector-based version
1675  //=====================================================================
1676  void locate_zeta_for_missing_coordinates
1677  (int& iproc, Mesh* const &external_mesh_pt, Problem* problem_pt,
1678  Vector<MeshAsGeomObject*>& mesh_geom_obj_pt)
1679  {
1680  // How many meshes are we dealing with?
1681  unsigned n_mesh=mesh_geom_obj_pt.size();
1682 
1683  // Storage for number of processors, current process and communicator
1684  OomphCommunicator* comm_pt=problem_pt->communicator_pt();
1685  int n_proc=comm_pt->nproc();
1686  int my_rank=comm_pt->my_rank();
1687 
1688  // Clear vectors containing data to be sent
1689  Flat_packed_doubles.resize(0);
1690 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1691  Flat_packed_unsigneds_string.resize(0);
1692 #endif
1693  Flat_packed_unsigneds.resize(0);
1694  Flat_packed_located_coordinates.resize(0);
1695 
1696  // Flush storage for zetas not found locally (when
1697  // processing the zeta coordinates received from "previous"
1698  // processor)
1699  Flat_packed_zetas_not_found_locally.resize(0);
1700 
1701  // Number of zeta tuples to be dealt with (includes padding!)
1702  unsigned n_zeta=Received_flat_packed_zetas_to_be_found.size()/Dim;
1703 
1704  // Create storage for the processor id (plus one) on which
1705  // the zetas stored in Flat_packed_zetas_not_found_locally[...]
1706  // were located. (Remains zero for padded entries).
1707  Proc_id_plus_one_of_external_element.resize(n_zeta,0);
1708 
1709  // Create storage for the status of the (external halo) element associated
1710  // the zetas stored in Flat_packed_zetas_not_found_locally[...].
1711  // It either hasn't been found, already exists on the processor
1712  // that needs it, or needs to be newly created. (Remains Not_found
1713  // for padded entries).
1714  Located_element_status.resize(n_zeta,Not_found);
1715 
1716  // Counter for flat-packed array of external zeta coordinates
1717  unsigned count=0;
1718 
1719  // Current mesh
1720  unsigned i_mesh=0;
1721 
1722  // Loop over the zeta tuples that we received from elsewhere and
1723  // are trying to find here for current mesh
1724  for (unsigned i=0;i<n_zeta;i++)
1725  {
1726  // Storage for global coordinates to be located
1727  Vector<double> x_global(Dim);
1728 
1729  // Loop to fill in coordinates
1730  for (unsigned ii=0;ii<Dim;ii++)
1731  {
1732  x_global[ii]=Received_flat_packed_zetas_to_be_found[count];
1733  count++;
1734  }
1735 
1736  // Check if we've reached the end of the mesh
1737  bool reached_end_of_mesh=false;
1738  unsigned dbl_max_count=0;
1739  for (unsigned ii=0;ii<Dim;ii++)
1740  {
1741  if (x_global[ii]==DBL_MAX)
1742  {
1743  dbl_max_count++;
1744  reached_end_of_mesh=true;
1745  }
1746  }
1747 
1748  // Reached end of mesh
1749  if (reached_end_of_mesh)
1750  {
1751 #ifdef PARANOID
1752  // Check if all coordinates were set to DBX_MAX
1753  if (dbl_max_count!=Dim)
1754  {
1755  std::ostringstream error_stream;
1756  error_stream << "Appear to have reached end of mesh " << i_mesh
1757  << " but only " << dbl_max_count << " out of "
1758  << Dim << " zeta coordinates have been set to DBX_MAX\n";
1759  throw OomphLibError
1760  (error_stream.str(),
1761  OOMPH_CURRENT_FUNCTION,
1762  OOMPH_EXCEPTION_LOCATION);
1763  }
1764 #endif
1765  // Indicate end of mesh in flat packed data
1766  for (unsigned ii=0;ii<Dim;ii++)
1767  {
1768  Flat_packed_zetas_not_found_locally.push_back(DBL_MAX);
1769  }
1770 
1771  // Bump mesh counter
1772  i_mesh++;
1773 
1774  // Bail out if we're done
1775  if (i_mesh==n_mesh)
1776  {
1777  return;
1778  }
1779  }
1780 
1781  // Perform locate_zeta for these coordinates and current mesh
1782  GeomObject *sub_geom_obj_pt=0;
1783  Vector<double> ss(Dim);
1784  if (!reached_end_of_mesh)
1785  {
1786  mesh_geom_obj_pt[i_mesh]->locate_zeta(x_global,
1787  sub_geom_obj_pt,ss);
1788 
1789  // Did the locate method work?
1790  if (sub_geom_obj_pt!=0)
1791  {
1792  // Get the source element - bulk or not?
1793  GeneralisedElement *source_el_pt=0;
1794  if (!Use_bulk_element_as_external)
1795  {
1796  source_el_pt=dynamic_cast<FiniteElement*>(sub_geom_obj_pt);
1797  }
1798  else
1799  {
1800  FaceElement *face_el_pt=dynamic_cast<FaceElement*>(sub_geom_obj_pt);
1801  source_el_pt=dynamic_cast<FiniteElement*>(face_el_pt->
1802  bulk_element_pt());
1803  }
1804 
1805  // Check if the returned element is halo
1806  if (!source_el_pt->is_halo()) // cannot accept halo here
1807  {
1808  // The correct non-halo element has been located; this will become
1809  // an external haloed element on the current process, and an
1810  // external halo copy needs to be created on the current process
1811  // minus wherever we are in the "ring-loop"
1812  int halo_copy_proc=my_rank-iproc;
1813 
1814  // If iproc is bigger than my_rank then we've "gone through" nproc-1
1815  if (my_rank<iproc) { halo_copy_proc=n_proc+halo_copy_proc; }
1816 
1817  // So, we found zeta on the current processor
1818  Proc_id_plus_one_of_external_element[i]=my_rank+1;
1819 
1820  // This source element is an external halo on process halo_copy_proc
1821  // but it should only be added to the storage if it hasn't
1822  // been added already, and this information also needs to be
1823  // communicated over to the other process
1824 
1825  unsigned n_extern_haloed=external_mesh_pt->
1826  nexternal_haloed_element(halo_copy_proc);
1827  unsigned external_haloed_el_index=
1828  external_mesh_pt->add_external_haloed_element_pt(halo_copy_proc,
1829  source_el_pt);
1830 
1831  // If it was added to the storage then the returned index
1832  // will be the same as the (old) size of the storage
1833  if (external_haloed_el_index==n_extern_haloed)
1834  {
1835  // Set Located_element_status to say it
1836  // should be newly created
1837  Located_element_status[i]=New;
1838 
1839  // How many continuously interpolated values are there?
1840  int n_cont_inter_values=-1;
1841  if (dynamic_cast<RefineableElement*>(source_el_pt)!=0)
1842  {
1843  n_cont_inter_values=dynamic_cast<RefineableElement*>
1844  (source_el_pt)->ncont_interpolated_values();
1845  }
1846 
1847  // Since it is (externally) haloed from the current process,
1848  // the info required to create a new element in the equivalent
1849  // external halo layer on process halo_copy_proc needs to be
1850  // sent there
1851 
1852  // If we're using macro elements to update...
1853  MacroElementNodeUpdateMesh* macro_mesh_pt=
1854  dynamic_cast<MacroElementNodeUpdateMesh*>(external_mesh_pt);
1855  if (macro_mesh_pt!=0)
1856  {
1857  Flat_packed_unsigneds.push_back(1);
1858 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1859  Flat_packed_unsigneds_string.push_back("Mesh is macro element mesh[2]");
1860 #endif
1861  //Cast to finite element... this must work because it's
1862  //a macroelement no update mesh
1863  FiniteElement* source_finite_el_pt
1864  = dynamic_cast<FiniteElement*>(source_el_pt);
1865 
1866  MacroElement* macro_el_pt=source_finite_el_pt->macro_elem_pt();
1867  // Send the macro element number across
1868  unsigned macro_el_num=macro_el_pt->macro_element_number();
1869  Flat_packed_unsigneds.push_back(macro_el_num);
1870 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1871  Flat_packed_unsigneds_string.push_back("Number of macro element[2]");
1872 #endif
1873  // we need to send
1874  // the lower left and upper right coordinates of the macro
1875  QElementBase* q_el_pt=dynamic_cast<QElementBase*>(source_el_pt);
1876  if (q_el_pt!=0)
1877  {
1878  // The macro element needs to be set first before
1879  // its lower left and upper right coordinates can be accessed
1880  // Now send the lower left and upper right coordinates
1881  unsigned el_dim=q_el_pt->dim();
1882  for (unsigned i_dim=0;i_dim<el_dim;i_dim++)
1883  {
1884  Flat_packed_doubles.push_back(q_el_pt->s_macro_ll(i_dim));
1885  Flat_packed_doubles.push_back(q_el_pt->s_macro_ur(i_dim));
1886  }
1887  }
1888  else // Throw an error
1889  {
1890  std::ostringstream error_stream;
1891  error_stream << "You are using a MacroElement node update\n"
1892  << "in a case with non-QElements. This has not\n"
1893  << "yet been implemented.\n";
1894  throw OomphLibError
1895  (error_stream.str(),
1896  OOMPH_CURRENT_FUNCTION,
1897  OOMPH_EXCEPTION_LOCATION);
1898  }
1899 
1900  }
1901  else // Not using macro elements to update
1902  {
1903  Flat_packed_unsigneds.push_back(0);
1904 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1905  Flat_packed_unsigneds_string.push_back("Mesh is not a macro element mesh [2]");
1906 #endif
1907  }
1908 
1909 
1910  //Cast to finite element... this must work because it's
1911  //a macroelement no update mesh
1912  FiniteElement* source_finite_el_pt
1913  = dynamic_cast<FiniteElement*>(source_el_pt);
1914 #ifdef PARANOID
1915  if(source_finite_el_pt==0)
1916  {
1917  throw
1918  OomphLibError(
1919  "Unable to cast source function to finite element\n",
1920  "Multi_domain_functions::locate_zeta_for_missing_coordinates()",
1921  OOMPH_EXCEPTION_LOCATION);
1922  }
1923 #endif
1924 
1925 
1926  // Loop over the nodes of the new source element
1927  unsigned n_node=source_finite_el_pt->nnode();
1928  for (unsigned j=0;j<n_node;j++)
1929  {
1930  Node* nod_pt=source_finite_el_pt->node_pt(j);
1931 
1932  // Add the node to the storage; this routine
1933  // also takes care of any master nodes if the
1934  // node is hanging
1935  add_external_haloed_node_to_storage(halo_copy_proc,nod_pt,
1936  problem_pt,
1937  external_mesh_pt,
1938  n_cont_inter_values);
1939  }
1940  }
1941  else // it has already been added, so tell the other process
1942  {
1943  // Set Located_element_status to indicate an element has
1944  // already been added
1945  Located_element_status[i]=Exists;
1946  Flat_packed_unsigneds.push_back(external_haloed_el_index);
1947 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
1948  Flat_packed_unsigneds_string.push_back("Index of existing external haloed element[2]");
1949 #endif
1950  }
1951 
1952  // The coordinates returned by locate_zeta are also needed
1953  // in the setup of the source elements on the other process
1954  if (!Use_bulk_element_as_external)
1955  {
1956  for (unsigned ii=0;ii<Dim;ii++)
1957  {
1958  Flat_packed_located_coordinates.push_back(ss[ii]);
1959  }
1960  }
1961  else // translate the coordinates to the bulk element
1962  {
1963  // The translation is from Lagrangian to Eulerian
1964  FaceElement *face_el_pt=
1965  dynamic_cast<FaceElement*>(sub_geom_obj_pt);
1966  //Get the dimension of the BulkElement
1967  unsigned bulk_el_dim =
1968  dynamic_cast<FiniteElement*>(source_el_pt)->dim();
1969  Vector<double> s_trans(bulk_el_dim);
1970  face_el_pt->get_local_coordinate_in_bulk(ss,s_trans);
1971  for (unsigned ii=0;ii<bulk_el_dim;ii++)
1972  {
1973  Flat_packed_located_coordinates.push_back(s_trans[ii]);
1974  }
1975  }
1976  }
1977  else // halo, so search again until non-halo equivalent is located
1978  {
1979  // Add required information to arrays (as below)
1980  for (unsigned ii=0;ii<Dim;ii++)
1981  {
1982  Flat_packed_zetas_not_found_locally.push_back(x_global[ii]);
1983  }
1984  // It wasn't found here
1985  Proc_id_plus_one_of_external_element[i]=0;
1986 
1987  // Set Located_element_status to indicate not found
1988  Located_element_status[i]=Not_found;
1989  }
1990  }
1991  else // not successful this time (i.e. sub_geom_obj_pt==0), so
1992  // prepare for next process to try
1993  {
1994  // Add this global coordinate to the LOCAL zeta array
1995  for (unsigned ii=0;ii<Dim;ii++)
1996  {
1997  Flat_packed_zetas_not_found_locally.push_back(x_global[ii]);
1998  }
1999  // It wasn't found here
2000  Proc_id_plus_one_of_external_element[i]=0;
2001 
2002  // Set Located_element_status to indicate not found
2003  Located_element_status[i]=Not_found;
2004  }
2005 
2006  } // end of mesh not reached
2007 
2008  } // end of loop over flat-packed zeta tuples
2009 
2010  }
2011 
2012 
2013 
2014 #endif
2015 
2016 
2017  //=====================================================================
2018  /// locate zeta for current set of "local" coordinates
2019  /// vector-based version
2020  //=====================================================================
2021  void locate_zeta_for_local_coordinates
2022  (const Vector<Mesh*>& mesh_pt, Mesh* const &external_mesh_pt,
2023  Vector<MeshAsGeomObject*>& mesh_geom_obj_pt,
2024  const unsigned& interaction_index)
2025  {
2026  // Flush storage for zetas not found locally
2027  Flat_packed_zetas_not_found_locally.resize(0);
2028 
2029  // Number of meshes
2030  unsigned n_mesh=mesh_pt.size();
2031 
2032 #ifdef PARANOID
2033  if (mesh_geom_obj_pt.size()!=n_mesh)
2034  {
2035  std::ostringstream error_stream;
2036  error_stream << "Sizes of mesh_geom_obj_pt [ "
2037  << mesh_geom_obj_pt.size() << " ] and "
2038  << "mesh_pt [ " << n_mesh << " ] don't match.\n";
2039  throw OomphLibError
2040  (error_stream.str(),
2041  OOMPH_CURRENT_FUNCTION,
2042  OOMPH_EXCEPTION_LOCATION);
2043  }
2044 #endif
2045 
2046  // Element counter
2047  unsigned e_count=0;
2048 
2049  // Loop over meshes
2050  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
2051  {
2052  // Number of local elements
2053  unsigned n_element=mesh_pt[i_mesh]->nelement();
2054 
2055  // Loop over this processor's elements
2056  for (unsigned e=0;e<n_element;e++)
2057  {
2058  ElementWithExternalElement *el_pt=
2059  dynamic_cast<ElementWithExternalElement*>(mesh_pt[i_mesh]->
2060  element_pt(e));
2061 #ifdef OOMPH_HAS_MPI
2062  // Only visit non-halo elements -- we're not setting up external elements
2063  // for on-halos!
2064  if (!el_pt->is_halo())
2065 #endif
2066  {
2067  // Find number of Gauss points and element dimension
2068  unsigned n_intpt=el_pt->integral_pt()->nweight();
2069  unsigned el_dim=el_pt->dim();
2070 
2071 
2072 #ifdef PARANOID
2073  if (el_dim!=Dim)
2074  {
2075  std::ostringstream error_stream;
2076  error_stream
2077  << "Dimension of element " << el_dim
2078  << " is not consitent with dimension assumed \n"
2079  << " in multidomain namespace, " << Dim << std::endl;
2080  throw OomphLibError
2081  (error_stream.str(),
2082  OOMPH_CURRENT_FUNCTION,
2083  OOMPH_EXCEPTION_LOCATION);
2084  }
2085 #endif
2086 
2087  // Set storage for local and global coordinates
2088  Vector<double> s_local(el_dim);
2089  Vector<double> x_global(el_dim);
2090 
2091  // Loop over integration points
2092  for (unsigned ipt=0;ipt<n_intpt;ipt++)
2093  {
2094  // Has this integration point been done yet?
2095  if (External_element_located[e_count][ipt]==0)
2096  {
2097  // Get local coordinates
2098  for (unsigned i=0;i<el_dim;i++)
2099  {
2100  s_local[i]=el_pt->integral_pt()->knot(ipt,i);
2101  }
2102  // Interpolate to global coordinates
2103  el_pt->interpolated_zeta(s_local,x_global);
2104 
2105  // Storage for geometric object and its local coordinates
2106  GeomObject* sub_geom_obj_pt=0;
2107  Vector<double> s_ext(el_dim);
2108  mesh_geom_obj_pt[i_mesh]->locate_zeta(x_global,
2109  sub_geom_obj_pt,s_ext);
2110 
2111  // Has the required element been located?
2112  if (sub_geom_obj_pt!=0)
2113  {
2114  // The required element has been located
2115  // The located coordinates have the same dimension as the bulk
2116  GeneralisedElement* source_el_pt;
2117  Vector<double> s_source(el_dim);
2118 
2119  // Is the bulk element the actual external element?
2120  if (!Use_bulk_element_as_external)
2121  {
2122  // Use the object directly (it must be a finite element)
2123  source_el_pt=dynamic_cast<FiniteElement*>(sub_geom_obj_pt);
2124  s_source=s_ext;
2125  }
2126  else
2127  {
2128  // Cast to a FaceElement and use the bulk element
2129  FaceElement* face_el_pt=
2130  dynamic_cast<FaceElement*>(sub_geom_obj_pt);
2131  source_el_pt=face_el_pt->bulk_element_pt();
2132 
2133  //Need to resize the located coordinates to have the same
2134  //dimension as the bulk element
2135  s_source.resize(dynamic_cast<FiniteElement*>
2136  (source_el_pt)->dim());
2137 
2138  // Translate the returned local coords into the bulk element
2139  face_el_pt->get_local_coordinate_in_bulk(s_ext,s_source);
2140  }
2141 
2142  // Check if it's a halo; if it is then the non-halo equivalent
2143  // needs to be located from another processor (unless we
2144  // accept halo elements as external elements)
2145 #ifdef OOMPH_HAS_MPI
2146  if (Allow_use_of_halo_elements_as_external_elements||
2147  (!source_el_pt->is_halo()))
2148 #endif
2149  {
2150  //Need to cast to a FiniteElement
2151  FiniteElement* source_finite_el_pt =
2152  dynamic_cast<FiniteElement*>(source_el_pt);
2153 
2154  // Set the external element pointer and local coordinates
2155  el_pt->external_element_pt(interaction_index,ipt)
2156  = source_finite_el_pt;
2157  el_pt->external_element_local_coord(interaction_index,ipt)
2158  =s_source;
2159 
2160  // Set the lookup array to 1/true
2161  External_element_located[e_count][ipt]=1;
2162  }
2163 #ifdef OOMPH_HAS_MPI
2164  // located element is halo and we're not accepting haloes
2165  // obviously only makes sense in mpi mode...
2166  else
2167  {
2168  // Add required information to arrays
2169  for (unsigned i=0;i<el_dim;i++)
2170  {
2171  Flat_packed_zetas_not_found_locally.push_back(x_global[i]);
2172  }
2173  }
2174 #endif
2175  }
2176  else
2177  {
2178  // Search has failed then add the required information to the
2179  // arrays which need to be sent to the other processors so that
2180  // they can perform the locate_zeta
2181 
2182  // Add this global coordinate to the LOCAL zeta array
2183  for (unsigned i=0;i<el_dim;i++)
2184  {
2185  Flat_packed_zetas_not_found_locally.push_back(x_global[i]);
2186  }
2187  }
2188  }
2189  } // end loop over integration points
2190  } //end for halo
2191 
2192  // Bump up counter for all elements
2193  e_count++;
2194 
2195  } // end loop over local elements
2196 
2197  // Mark end of mesh data in flat packed array
2198  for (unsigned i=0;i<Dim;i++)
2199  {
2200  Flat_packed_zetas_not_found_locally.push_back(DBL_MAX);
2201  }
2202 
2203  } // end of loop over meshes
2204  }
2205 
2206 
2207 //=====================================================================
2208 /// Helper function that computes the dimension of the elements within
2209 /// each of the specified meshes (and checks they are the same).
2210 /// Stores result in Dim.
2211 //=====================================================================
2212  void get_dim_helper(Problem* problem_pt, Mesh* const &mesh_pt,
2213  Mesh* const &external_mesh_pt)
2214  {
2215 #ifdef OOMPH_HAS_MPI
2216  // Storage for number of processors, current process and communicator
2217  OomphCommunicator* comm_pt=problem_pt->communicator_pt();
2218 #endif
2219 
2220  // Extract the element dimensions from the first element of each mesh
2221  unsigned mesh_dim=0;
2222  if (mesh_pt->nelement() > 0)
2223  {
2224  mesh_dim=
2225  dynamic_cast<FiniteElement*>(mesh_pt->element_pt(0))->dim();
2226  }
2227  unsigned external_mesh_dim=0;
2228  if (external_mesh_pt->nelement() > 0)
2229  {
2230  external_mesh_dim=
2231  dynamic_cast<FiniteElement*>(external_mesh_pt->element_pt(0))->dim();
2232  }
2233 
2234  // Need to do an Allreduce
2235 #ifdef OOMPH_HAS_MPI
2236  int n_proc=comm_pt->nproc();
2237  if (n_proc > 1)
2238  {
2239  unsigned mesh_dim_reduce;
2240  MPI_Allreduce(&mesh_dim,&mesh_dim_reduce,1,MPI_UNSIGNED,
2241  MPI_MAX,comm_pt->mpi_comm());
2242  mesh_dim=mesh_dim_reduce;
2243 
2244  unsigned external_mesh_dim_reduce;
2245  MPI_Allreduce(&external_mesh_dim,&external_mesh_dim_reduce,1,MPI_UNSIGNED,
2246  MPI_MAX,comm_pt->mpi_comm());
2247  external_mesh_dim=external_mesh_dim_reduce;
2248  }
2249 #endif
2250 
2251  // Check the dimensions are the same!
2252  if (mesh_dim!=external_mesh_dim)
2253  {
2254  std::ostringstream error_stream;
2255  error_stream << "The elements within the two meshes do not\n"
2256  << "have the same dimension, so the multi-domain\n"
2257  << "method will not work.\n"
2258  << "For the mesh, dim=" << mesh_dim
2259  << ", and the external mesh, dim=" << external_mesh_dim
2260  << "\n";
2261  throw OomphLibError
2262  (error_stream.str(),
2263  OOMPH_CURRENT_FUNCTION,
2264  OOMPH_EXCEPTION_LOCATION);
2265  }
2266 
2267  // Set dimension
2268  Dim=mesh_dim;
2269  }
2270 
2271 
2272  //=================================================================
2273  /// Helper function that clears all the information used
2274  /// during the external storage creation
2275  //=================================================================
2276  void clean_up()
2277  {
2278  Flat_packed_zetas_not_found_locally.clear();
2279  Received_flat_packed_zetas_to_be_found.clear();
2280  Proc_id_plus_one_of_external_element.clear();
2281  Located_element_status.clear();
2282  Flat_packed_located_coordinates.clear();
2283  Flat_packed_doubles.clear();
2284  Flat_packed_unsigneds.clear();
2285  External_element_located.clear();
2286  }
2287 
2288  /// Vector of zeta coordinates that we're currently trying to locate;
2289  /// used in sorting of bin entries in further_away() comparison function
2290  Vector<double> Zeta_coords_for_further_away_comparison;
2291 
2292  /// Comparison function for sorting entries in bin: Returns true if
2293  /// point identified by p1 (comprising pointer to finite element and
2294  /// vector of local coordinates within that element) is closer to
2295  /// Zeta_coords_for_further_away_comparison than p2
2296  bool first_closer_than_second(
2297  const std::pair<FiniteElement*,Vector<double> >& p1,
2298  const std::pair<FiniteElement*,Vector<double> >& p2)
2299  {
2300  // First point
2301  FiniteElement* el_pt=p1.first;
2302  Vector<double> s(p1.second);
2303  Vector<double> zeta(Dim);
2304  el_pt->interpolated_zeta(s,zeta);
2305  double dist_squared1=0.0;
2306  for (unsigned i=0;i<Dim;i++)
2307  {
2308  dist_squared1+=
2309  (zeta[i]-Zeta_coords_for_further_away_comparison[i])*
2310  (zeta[i]-Zeta_coords_for_further_away_comparison[i]);
2311  }
2312 
2313  // Second point
2314  el_pt=p2.first;
2315  s=p2.second;
2316  el_pt->interpolated_zeta(s,zeta);
2317  double dist_squared2=0.0;
2318  for (unsigned i=0;i<Dim;i++)
2319  {
2320  dist_squared2+=
2321  (zeta[i]-Zeta_coords_for_further_away_comparison[i])*
2322  (zeta[i]-Zeta_coords_for_further_away_comparison[i]);
2323  }
2324 
2325  // Which one is further
2326  if (dist_squared1<dist_squared2)
2327  {
2328  return true;
2329  }
2330  else
2331  {
2332  return false;
2333  }
2334  }
2335 }
2336 
2337 }
oomph::GeneralisedElement
A Generalised Element class.
Definition: elements.h:76
oomph::FiniteElement::local_coordinate_of_node
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:1803
oomph::Multi_domain_functions::add_external_haloed_node_helper
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.
Definition: multi_domain.cc:532
oomph::AlgebraicNode::geom_object_pt
GeomObject * geom_object_pt(const unsigned &i)
Return pointer to i-th geometric object involved in default (usually first) update function...
Definition: algebraic_elements.h:238
oomph::ElementWithExternalElement::external_element_local_coord
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element&#39;s local coords for specified interaction index at specified int...
Definition: element_with_external_element.h:146
oomph::Mesh::add_boundary_node
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
oomph::MacroElement::macro_element_number
unsigned & macro_element_number()
Access function to the Macro_element_number.
Definition: macro_element.h:214
oomph::Mesh::add_external_halo_node_pt
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:2016
oomph::Multi_domain_functions::Received_flat_packed_zetas_to_be_found
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
oomph::Multi_domain_functions::first_closer_than_second
bool first_closer_than_second(const std::pair< FiniteElement *, Vector< double > > &p1, const std::pair< FiniteElement *, Vector< double > > &p2)
Definition: multi_domain.cc:2296
oomph::AlgebraicMesh::update_node_update
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 ...
oomph::Node::hanging_pt
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
oomph::Mesh::add_external_haloed_element_pt
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:9087
oomph::Multi_domain_functions::send_and_receive_located_info
void send_and_receive_located_info(int &iproc, Mesh *const &external_mesh_pt, Problem *problem_pt)
Helper function to send back any located information.
Definition: multi_domain.cc:248
i
cstr elem_len * i
Definition: cfortran.h:607
oomph::Multi_domain_functions::Counter_for_flat_packed_unsigneds
unsigned Counter_for_flat_packed_unsigneds
Counter used when processing vector of flat-packed unsigneds – this is really "private" data...
Definition: multi_domain.cc:140
oomph::Problem::ntime_stepper
unsigned ntime_stepper() const
Return the number of time steppers.
Definition: problem.h:1578
oomph::AlgebraicNode::add_node_update_info
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&#39;s the ID of the mesh update function (typically used ...
Definition: algebraic_elements.h:289
oomph::AlgebraicNode::nref_value
unsigned nref_value(const int &id)
Number of reference values involved in id-th update function.
Definition: algebraic_elements.h:244
oomph::Multi_domain_functions::add_external_haloed_master_node_helper
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.
Definition: multi_domain.cc:590
oomph::Multi_domain_functions::get_dim_helper
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 ...
Definition: multi_domain.cc:2212
oomph::Problem
The Problem class.
Definition: problem.h:152
oomph::FiniteElement
A general Finite Element class.
Definition: elements.h:1274
oomph::Multi_domain_functions::Doc_timings
bool Doc_timings
Boolean to indicate whether to doc timings or not.
Definition: multi_domain.cc:164
t
char t
Definition: cfortran.h:572
oomph::Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_domain_interaction
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
oomph::HangInfo::master_weight
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
oomph::GeneralisedElement::is_halo
bool is_halo() const
Is this element a halo?
Definition: elements.h:1141
oomph::FiniteElement::construct_boundary_node
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:2420
oomph::Node
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
oomph::HangInfo::nmaster
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:311
oomph::Node::is_hanging
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
oomph::Data::nvalue
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
oomph::Multi_domain_functions::Doc_boundary_coordinate_file
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
e
e
Definition: cfortran.h:575
oomph::Problem::communicator_pt
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1221
oomph::Multi_domain_functions::Flat_packed_unsigneds
Vector< unsigned > Flat_packed_unsigneds
Vector of flat-packed unsigneds to be communicated with other processors – this is really "private" ...
Definition: multi_domain.cc:124
oomph::AlgebraicMesh::geom_object_list_pt
GeomObject * geom_object_list_pt(const unsigned &i)
Access function to the ith GeomObject.
Definition: algebraic_elements.h:879
oomph::FaceElement::get_local_coordinate_in_bulk
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:6242
multi_domain.template.cc
oomph::AlgebraicNode::node_update_fct_id
int node_update_fct_id()
Default (usually first if there are multiple ones) node update fct id.
Definition: algebraic_elements.h:161
oomph::Multi_domain_functions::Dim
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in...
Definition: multi_domain.cc:66
oomph::Mesh::nelement
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
oomph::Multi_domain_functions::Doc_full_stats
bool Doc_full_stats
Boolean to indicate whether to output further info during setup_multi_domain_interaction() routines...
Definition: multi_domain.cc:172
oomph::Node::ndim
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
oomph::QElementBase
Base class for Qelements.
Definition: Qelements.h:106
oomph::QElementBase::s_macro_ur
double & s_macro_ur(const unsigned &i)
Access fct to the i-th coordinate of the element&#39;s "upper right" vertex in the associated MacroElemen...
Definition: Qelements.h:228
oomph::Node::x
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
oomph::Mesh::add_external_haloed_node_pt
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:9128
oomph::BoundaryNodeBase
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:1855
multi_domain.h
oomph::Multi_domain_functions::add_external_haloed_node_to_storage
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.
Definition: multi_domain.cc:446
oomph::Multi_domain_functions::Allow_use_of_halo_elements_as_external_elements
bool Allow_use_of_halo_elements_as_external_elements
Boolean to indicate if we&#39;re allowed to use halo elements as external elements. Can drastically reduc...
Definition: multi_domain.cc:155
oomph::Mesh::element_pt
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
oomph::Multi_domain_functions::Flat_packed_located_coordinates
Vector< double > Flat_packed_located_coordinates
Vector of flat-packed local coordinates for zeta tuples that have been located.
Definition: multi_domain.cc:108
oomph::Multi_domain_functions::Flat_packed_unsigneds_string
Vector< std::string > Flat_packed_unsigneds_string
Definition: multi_domain.cc:131
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1908
oomph::FiniteElement::construct_node
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:2392
oomph::FiniteElement::interpolated_zeta
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:4564
oomph::Mesh::external_halo_node_pt
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:2025
oomph::Mesh::nboundary
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:806
s
static char t char * s
Definition: cfortran.h:572
oomph::Mesh::nexternal_haloed_node
unsigned nexternal_haloed_node()
Total number of external haloed nodes in this Mesh.
Definition: mesh.h:2060
oomph::Integral::knot
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
oomph::Data::set_value
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
oomph::Multi_domain_functions::clean_up
void clean_up()
Helper function that clears all the intermediate information used during the external storage creatio...
Definition: multi_domain.cc:2276
oomph::Node::nposition_type
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:966
oomph::Multi_domain_functions::get_required_master_nodal_information_helper
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...
Definition: multi_domain.cc:899
oomph::Multi_domain_functions::recursively_add_masters_of_external_haloed_node
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.
Definition: multi_domain.cc:468
oomph::SolidNode::variable_position_pt
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
Definition: nodes.h:1655
oomph::Multi_domain_functions::Not_found
Definition: multi_domain.h:153
oomph::Problem::time_stepper_pt
TimeStepper *& time_stepper_pt()
Access function for the pointer to the first (presumably only) timestepper.
Definition: problem.h:1460
oomph::Multi_domain_functions::send_and_receive_missing_zetas
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...
Definition: multi_domain.cc:184
oomph::HangInfo::master_node_pt
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
oomph::MacroElementNodeUpdateNode::node_update_element_pt
FiniteElement *& node_update_element_pt()
Pointer to finite element that performs the update by referring to its macro-element representation (...
Definition: macro_element_node_update_element.h:116
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2109
oomph::SolidNode::nlagrangian_type
unsigned nlagrangian_type() const
Number of types of Lagrangian coordinates used to interpolate the Lagrangian coordinates within the e...
Definition: nodes.h:1738
oomph::FiniteElement::macro_elem_pt
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
Definition: elements.h:1837
oomph::Multi_domain_functions::Flat_packed_doubles
Vector< double > Flat_packed_doubles
Vector of flat-packed doubles to be communicated with other processors.
Definition: multi_domain.cc:112
oomph::HangInfo
Class that contains data for hanging nodes.
Definition: nodes.h:684
oomph::TimeStepper::ntstorage
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
oomph::Multi_domain_functions::construct_new_external_halo_node_helper
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...
Definition: multi_domain.cc:1333
oomph::Multi_domain_functions::get_required_nodal_information_helper
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...
Definition: multi_domain.cc:643
oomph::Multi_domain_functions::Use_bulk_element_as_external
bool Use_bulk_element_as_external
Boolean to indicate when to use the bulk element as the external element. Defaults to false...
Definition: multi_domain.cc:148
oomph::Multi_domain_functions::Allow_use_of_halo_elements_as_external_elements_for_projection
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.
Definition: multi_domain.cc:161
oomph::MacroElementNodeUpdateNode::set_node_update_info
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...
Definition: macro_element_node_update_element.h:170
oomph::FiniteElement::dim
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2482
oomph::QElementBase::s_macro_ll
double & s_macro_ll(const unsigned &i)
Access fct to the i-th coordinate of the element&#39;s "lower left" vertex in the associated MacroElement...
Definition: Qelements.h:211
oomph::AlgebraicNode::ref_value
double ref_value(const unsigned &i)
Return i-th reference value involved in default (usually first) update function.
Definition: algebraic_elements.h:276
oomph::SolidNode
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:1569
oomph::AlgebraicMesh::ngeom_object_list_pt
unsigned ngeom_object_list_pt()
Return number of geometric objects associated with AlgebraicMesh.
Definition: algebraic_elements.h:873
oomph::Multi_domain_functions::Proc_id_plus_one_of_external_element
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
oomph::Node::resize
void resize(const unsigned &n_value)
Resize the number of equations.
Definition: nodes.cc:2089
oomph::Integral::nweight
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
oomph::Multi_domain_functions::External_element_located
Vector< Vector< unsigned > > External_element_located
Lookup scheme for whether a local element&#39;s integration point has had an external element assigned to...
Definition: multi_domain.cc:74
oomph::ElementWithExternalElement::external_element_pt
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point...
Definition: element_with_external_element.h:117
oomph::Data::time_stepper_pt
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
oomph::AlgebraicNode::ngeom_object
unsigned ngeom_object(const int &id)
Number of geometric objects involved in id-th update function.
Definition: algebraic_elements.h:199
oomph::FaceElement::bulk_element_pt
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4515
algebraic_elements.h
macro_element_node_update_element.h
oomph::Node::is_on_boundary
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:1290
oomph::SolidNode::nlagrangian
unsigned nlagrangian() const
Return number of lagrangian coordinates.
Definition: nodes.h:1734
oomph::BoundaryNodeBase::index_of_first_value_assigned_by_face_element_pt
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:1910
oomph::Multi_domain_functions::Flat_packed_zetas_not_found_locally
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
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2146
oomph::TimeStepper
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
oomph::Multi_domain_functions::Doc_stats
bool Doc_stats
Boolean to indicate whether to output basic info during setup_multi_domain_interaction() routines...
Definition: multi_domain.cc:168
oomph::Multi_domain_functions::Located_element_status
Vector< unsigned > Located_element_status
Vector to indicate (to another processor) whether a located element (that will have to represented as...
Definition: multi_domain.cc:104
oomph::Multi_domain_functions::New
Definition: multi_domain.h:153
oomph::MacroElementNodeUpdateMesh::geom_object_vector_pt
Vector< GeomObject * > geom_object_vector_pt()
Access function to the vector of GeomObject.
Definition: macro_element_node_update_element.h:539
oomph::Mesh
A general mesh class.
Definition: mesh.h:74
oomph::Multi_domain_functions::add_external_halo_node_helper
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.
Definition: multi_domain.cc:1286
oomph::Node::value
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:2328
oomph::Multi_domain_functions::locate_zeta_for_missing_coordinates
void locate_zeta_for_missing_coordinates(int &iproc, Mesh *const &external_mesh_pt, Problem *problem_pt, Vector< MeshAsGeomObject *> &mesh_geom_obj_pt)
Locate zeta for current set of missing coordinates; vector-based version.
Definition: multi_domain.cc:1677
oomph::Multi_domain_functions::Zeta_coords_for_further_away_comparison
Vector< double > Zeta_coords_for_further_away_comparison
Definition: multi_domain.cc:2290
oomph::Multi_domain_functions::Counter_for_flat_packed_doubles
unsigned Counter_for_flat_packed_doubles
Counter used when processing vector of flat-packed doubles – this is really "private" data...
Definition: multi_domain.cc:118
oomph::Multi_domain_functions::Exists
Definition: multi_domain.h:153
oomph::Multi_domain_functions::locate_zeta_for_local_coordinates
void locate_zeta_for_local_coordinates(const Vector< Mesh *> &mesh_pt, Mesh *const &external_mesh_pt, Vector< MeshAsGeomObject *> &mesh_geom_obj_pt, const unsigned &interaction_index)
Helper function to locate "local" zeta coordinates This is the vector-based version which operates si...
Definition: multi_domain.cc:2022
oomph::OomphCommunicator
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:57