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