multi_domain.template.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1132 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-02-23 05:43:26 +0000 (Tue, 23 Feb 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Templated multi-domain functions
31 
32 //Include guards to prevent multiple inclusion of the header
33 #ifndef OOMPH_MULTI_DOMAIN_CC
34 #define OOMPH_MULTI_DOMAIN_CC
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38 #include <oomph-lib-config.h>
39 #endif
40 
41 //Oomph-lib headers
42 #include "geom_objects.h"
43 #include "problem.h"
44 #include "shape.h"
45 
46 #include "mesh.h"
48 #include "algebraic_elements.h"
50 #include "Qelements.h"
52 #include "multi_domain.h"
54 
55 //Needed to check if elements have nonuniformly spaced nodes
56 #include "refineable_elements.h"
57 #include "Qspectral_elements.h"
58 
59 namespace oomph
60 {
61 
62 
63 //// Templated helper functions for multi-domain methods using locate_zeta
64 
65 
66  //============================================================================
67  /// Identify the \c FaceElements (stored in the mesh pointed to by
68  /// \c face_mesh_pt) that are adjacent to the bulk elements next to the
69  /// \c boundary_in_bulk_mesh -th boundary of the mesh pointed to by
70  /// \c bulk_mesh_pt. The \c FaceElements must be derived
71  /// from the \c ElementWithExternalElement base class and the adjacent
72  /// bulk elements are stored as their external elements.
73  ///
74  /// This is the vector-based version which deals with multiple bulk
75  /// mesh boundaries at the same time.
76  //============================================================================
77  template<class BULK_ELEMENT, unsigned DIM>
79  Problem* problem_pt,
80  Vector<unsigned>& boundary_in_bulk_mesh,
81  Mesh* const &bulk_mesh_pt,
82  Vector<Mesh*>& face_mesh_pt,
83  const unsigned& interaction)
84  {
85 
86  unsigned n_mesh=boundary_in_bulk_mesh.size();
87 
88 #ifdef PARANOID
89  // Check sizes match
90  if (boundary_in_bulk_mesh.size()!=face_mesh_pt.size())
91  {
92  std::ostringstream error_message;
93  error_message
94  << "Sizes of vector of boundary ids in bulk mesh ("
95  << boundary_in_bulk_mesh.size() << ") and vector of pointers\n"
96  << "to FaceElements (" << face_mesh_pt.size() << " doesn't match.\n";
97  throw OomphLibError(
98  error_message.str(),
99  OOMPH_CURRENT_FUNCTION,
100  OOMPH_EXCEPTION_LOCATION);
101  }
102 #endif
103 
104  // Create face meshes adjacent to the bulk mesh's b-th boundary.
105  // Each face mesh consists of FaceElements that may also be
106  // interpreted as GeomObjects
107  Vector<Mesh*> bulk_face_mesh_pt(n_mesh);
108 
109  // Loop over all meshes
110  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
111  {
112  bulk_face_mesh_pt[i_mesh] = new Mesh;
113  bulk_mesh_pt->template build_face_mesh
114  <BULK_ELEMENT,FaceElementAsGeomObject>(boundary_in_bulk_mesh[i_mesh],
115  bulk_face_mesh_pt[i_mesh]);
116 
117  // Loop over these new face elements and tell them the boundary number
118  // from the bulk mesh -- this is required to they can
119  // get access to the boundary coordinates!
120  unsigned n_face_element = bulk_face_mesh_pt[i_mesh]->nelement();
121  for(unsigned e=0;e<n_face_element;e++)
122  {
123  //Cast the element pointer to the correct thing!
126  (bulk_face_mesh_pt[i_mesh]->element_pt(e));
127 
128  // Set bulk boundary number
129  el_pt->set_boundary_number_in_bulk_mesh(boundary_in_bulk_mesh[i_mesh]);
130 
131  // Doc?
132  if (Doc_boundary_coordinate_file.is_open())
133  {
134  Vector<double> s(DIM-1);
135  Vector<double> zeta(DIM-1);
136  Vector<double> x(DIM);
137  unsigned n_plot=5;
139 
140  // Loop over plot points
141  unsigned num_plot_points=el_pt->nplot_points(n_plot);
142  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
143  {
144  // Get local coordinates of plot point
145  el_pt->get_s_plot(iplot,n_plot,s);
146  el_pt->interpolated_zeta(s,zeta);
147  el_pt->interpolated_x(s,x);
148  for (unsigned i=0;i<DIM;i++)
149  {
150  Doc_boundary_coordinate_file << x[i] << " ";
151  }
152  for (unsigned i=0;i<DIM-1;i++)
153  {
154  Doc_boundary_coordinate_file << zeta[i] << " ";
155  }
156  Doc_boundary_coordinate_file << std::endl;
157  }
159  }
160  }
161 
162  // Now sort the elements based on the boundary coordinates.
163  // This may allow a faster implementation of the locate_zeta
164  // function for the MeshAsGeomObject representation of this
165  // mesh, but also creates a unique ordering of the elements
166  std::sort(bulk_face_mesh_pt[i_mesh]->element_pt().begin(),
167  bulk_face_mesh_pt[i_mesh]->element_pt().end(),
169  } // end of loop over meshes
170 
171 
172 
173  // Setup the interactions for this problem using the surface mesh
174  // on the fluid mesh. The interaction parameter in this instance is
175  // given by the "face" parameter.
178  (problem_pt,face_mesh_pt,bulk_mesh_pt,
179  bulk_face_mesh_pt,interaction);
180 
181 
182  // Loop over all meshes to clean up
183  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
184  {
185  unsigned n_face_element = bulk_face_mesh_pt[i_mesh]->nelement();
186 
187  //The MeshAsGeomObject has already been deleted (in set_external_storage)
188 
189  //Must be careful with the FaceMesh, because we cannot delete the nodes
190  //Loop over the elements backwards (paranoid!) and delete them
191  for(unsigned e=n_face_element;e>0;e--)
192  {
193  delete bulk_face_mesh_pt[i_mesh]->element_pt(e-1);
194  bulk_face_mesh_pt[i_mesh]->element_pt(e-1) = 0;
195  }
196  //Now clear all element and node storage (should maybe fine-grain this)
197  bulk_face_mesh_pt[i_mesh]->flush_element_and_node_storage();
198 
199  //Finally delete the mesh
200  delete bulk_face_mesh_pt[i_mesh];
201 
202  } // end of loop over meshes
203 
204  }
205 
206 
207 //========================================================================
208  /// Identify the \c FaceElements (stored in the mesh pointed to by
209  /// \c face_mesh_pt) that are adjacent to the bulk elements next to the
210  /// \c boundary_in_bulk_mesh -th boundary of the mesh pointed to by
211  /// \c bulk_mesh_pt. The \c FaceElements must be derived
212  /// from the \c ElementWithExternalElement base class and the adjacent
213  /// bulk elements are stored as their external elements.
214 //========================================================================
215  template<class BULK_ELEMENT,unsigned DIM>
217  Problem* problem_pt,
218  const unsigned& boundary_in_bulk_mesh,
219  Mesh* const& bulk_mesh_pt,
220  Mesh* const& face_mesh_pt,
221  const unsigned& interaction)
222  {
223 
224 #ifdef USE_VECTOR_BASED_MD
225 
226  // Convert to vector-based storage
227  Vector<unsigned> boundary_in_bulk_mesh_vect(1);
228  boundary_in_bulk_mesh_vect[0]=boundary_in_bulk_mesh;
229  Vector<Mesh*> face_mesh_pt_vect(1);
230  face_mesh_pt_vect[0]=face_mesh_pt;
231 
232  // Call vector-based version
233  setup_bulk_elements_adjacent_to_face_mesh<BULK_ELEMENT,DIM>(
234  problem_pt, boundary_in_bulk_mesh_vect, bulk_mesh_pt,
235  face_mesh_pt_vect,interaction);
236 
237 #else
238 
239  // Create a face mesh adjacent to the bulk mesh's b-th boundary.
240  // The face mesh consists of FaceElements that may also be
241  // interpreted as GeomObjects
242  Mesh* bulk_face_mesh_pt = new Mesh;
243  bulk_mesh_pt->template build_face_mesh
244  <BULK_ELEMENT,FaceElementAsGeomObject>(boundary_in_bulk_mesh,
245  bulk_face_mesh_pt);
246 
247  // Loop over these new face elements and tell them the boundary number
248  // from the bulk mesh -- this is required to they can
249  // get access to the boundary coordinates!
250  unsigned n_face_element = bulk_face_mesh_pt->nelement();
251  for(unsigned e=0;e<n_face_element;e++)
252  {
253  //Cast the element pointer to the correct thing!
256  (bulk_face_mesh_pt->element_pt(e));
257 
258  // Set bulk boundary number
259  el_pt->set_boundary_number_in_bulk_mesh(boundary_in_bulk_mesh);
260 
261  // Doc?
262  if (Doc_boundary_coordinate_file.is_open())
263  {
264  Vector<double> s(DIM-1);
265  Vector<double> zeta(DIM-1);
266  Vector<double> x(DIM);
267  unsigned n_plot=5;
269 
270  // Loop over plot points
271  unsigned num_plot_points=el_pt->nplot_points(n_plot);
272  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
273  {
274  // Get local coordinates of plot point
275  el_pt->get_s_plot(iplot,n_plot,s);
276  el_pt->interpolated_zeta(s,zeta);
277  el_pt->interpolated_x(s,x);
278  for (unsigned i=0;i<DIM;i++)
279  {
280  Doc_boundary_coordinate_file << x[i] << " ";
281  }
282  for (unsigned i=0;i<DIM-1;i++)
283  {
284  Doc_boundary_coordinate_file << zeta[i] << " ";
285  }
286  Doc_boundary_coordinate_file << std::endl;
287  }
289  }
290  }
291 
292  // Now sort the elements based on the boundary coordinates.
293  // This may allow a faster implementation of the locate_zeta
294  // function for the MeshAsGeomObject representation of this
295  // mesh, but also creates a unique ordering of the elements
296  std::sort(bulk_face_mesh_pt->element_pt().begin(),
297  bulk_face_mesh_pt->element_pt().end(),
299 
300  // Setup the interactions for this problem using the surface mesh
301  // on the bulk mesh.
304  (problem_pt,face_mesh_pt,bulk_mesh_pt,bulk_face_mesh_pt,interaction);
305 
306  // The source elements and coordinates have now all been set
307 
308  //Clean up the memory allocated:
309 
310  //The MeshAsGeomObject has already been deleted (in set_external_storage)
311 
312  //Must be careful with the FaceMesh, because we cannot delete the nodes
313  //Loop over the elements backwards (paranoid!) and delete them
314  for(unsigned e=n_face_element;e>0;e--)
315  {
316  delete bulk_face_mesh_pt->element_pt(e-1);
317  bulk_face_mesh_pt->element_pt(e-1) = 0;
318  }
319  //Now clear all element and node storage (should maybe fine-grain this)
320  bulk_face_mesh_pt->flush_element_and_node_storage();
321 
322  //Finally delete the mesh
323  delete bulk_face_mesh_pt;
324 
325 #endif
326 
327  }
328 
329 
330 //========================================================================
331 /// Set up the two-way multi-domain interactions for the
332 /// problem pointed to by \c problem_pt.
333 /// Use this for cases where first_mesh_pt and second_mesh_pt
334 /// occupy the same physical space and are populated by
335 /// ELEMENT_0 and ELEMENT_1 respectively, and are combined to solve
336 /// a single problem. The elements in two meshes interact both ways
337 /// the elements in each mesh act as "external elements" for the
338 /// elements in the "other" mesh. The interaction indices allow the
339 /// specification of which interaction we're setting up in the two
340 /// meshes. They default to zero, which is appropriate if there's
341 /// only a single interaction.
342 //========================================================================
343  template<class ELEMENT_0,class ELEMENT_1>
345  (Problem* problem_pt,Mesh* const &first_mesh_pt, Mesh* const &second_mesh_pt,
346  const unsigned& first_interaction, const unsigned& second_interaction)
347  {
348  // Delete all the current external halo(ed) element and node storage
349  first_mesh_pt->delete_all_external_storage();
350  second_mesh_pt->delete_all_external_storage();
351 
352  // Call setup_multi_domain_interaction in both directions
353  setup_multi_domain_interaction<ELEMENT_1>
354  (problem_pt,first_mesh_pt,second_mesh_pt,first_interaction);
355 
356  setup_multi_domain_interaction<ELEMENT_0>
357  (problem_pt,second_mesh_pt,first_mesh_pt,second_interaction);
358 
359  }
360 
361 
362 //========================================================================
363  /// Function to set up the one-way multi-domain interaction for
364  /// problems where the meshes pointed to by \c mesh_pt and \c external_mesh_pt
365  /// occupy the same physical space, and the elements in \c external_mesh_pt
366  /// act as "external elements" for the \c ElementWithExternalElements
367  /// in \c mesh_pt (but not vice versa):
368  /// - \c mesh_pt points to the mesh of ElemenWithExternalElements for which
369  /// the interaction is set up.
370  /// - \c external_mesh_pt points to the mesh that contains the elements
371  /// of type EXT_ELEMENT that act as "external elements" for the
372  /// \c ElementWithExternalElements in \c mesh_pt.
373  /// - The interaction_index parameter defaults to zero and must be otherwise
374  /// set by the user if there is more than one mesh that provides sources
375  /// for the Mesh pointed to by mesh_pt.
376 //========================================================================
377  template<class EXT_ELEMENT>
379  (Problem* problem_pt, Mesh* const &mesh_pt, Mesh* const &external_mesh_pt,
380  const unsigned& interaction_index)
381  {
382  // Bulk elements must not be external elements in this case
384 
385  // Call the auxiliary function with GEOM_OBJECT=EXT_ELEMENT
386  // and EL_DIM_EUL=EL_DIM_LAG=dimension returned from helper function
387  get_dim_helper(problem_pt,mesh_pt,external_mesh_pt);
388 
389  if(Dim > 3)
390  {
391  std::ostringstream error_stream;
392  error_stream << "The elements within the two interacting meshes have a\n"
393  << "dimension not equal to 1, 2 or 3.\n"
394  << "The multi-domain method will not work in this case.\n"
395  << "The dimension is: " << Dim << "\n";
396  throw OomphLibError
397  (error_stream.str(),
398  OOMPH_CURRENT_FUNCTION,
399  OOMPH_EXCEPTION_LOCATION);
400  }
401 
402  // Wrapper for each dimension (template parameter)
403  aux_setup_multi_domain_interaction<EXT_ELEMENT,EXT_ELEMENT>
404  (problem_pt,mesh_pt,external_mesh_pt,interaction_index);
405 
406  }
407 
408 
409 //========================================================================
410  /// Function to set up the one-way multi-domain interaction for
411  /// FSI-like problems.
412  /// - \c mesh_pt points to the mesh of \c ElemenWithExternalElements for which
413  /// the interaction is set up. In an FSI example, this mesh would contain
414  /// the \c FSIWallElements (either beam/shell elements or the
415  /// \c FSISolidTractionElements that apply the traction to
416  /// a "bulk" solid mesh that is loaded by the fluid.)
417  /// - \c external_mesh_pt points to the mesh that contains the elements
418  /// of type EXT_ELEMENT that provide the "source" for the
419  /// \c ElementWithExternalElements. In an FSI example, this
420  /// mesh would contain the "bulk" fluid elements.
421  /// - \c external_face_mesh_pt points to the mesh of \c FaceElements
422  /// attached to the \c external_mesh_pt. The mesh pointed to by
423  /// \c external_face_mesh_pt has the same dimension as \c mesh_pt.
424  /// The elements contained in \c external_face_mesh_pt are of type
425  /// FACE_ELEMENT_GEOM_OBJECT. In an FSI example, these elements
426  /// are usually the \c FaceElementAsGeomObjects (templated by the
427  /// type of the "bulk" fluid elements to which they are attached)
428  /// that define the FSI boundary of the fluid domain.
429  /// - The interaction_index parameter defaults to zero and must otherwise be
430  /// set by the user if there is more than one mesh that provides "external
431  /// elements" for the Mesh pointed to by mesh_pt (e.g. in the case
432  /// when a beam or shell structure is loaded by fluid from both sides.)
433 //========================================================================
434  template<class EXT_ELEMENT,class FACE_ELEMENT_GEOM_OBJECT>
436  (Problem* problem_pt, Mesh* const &mesh_pt, Mesh* const &external_mesh_pt,
437  Mesh* const &external_face_mesh_pt, const unsigned& interaction_index)
438  {
439  // Bulk elements must be external elements in this case
441 
442  // Call the auxiliary routine with GEOM_OBJECT=FACE_ELEMENT_GEOM_OBJECT
443  // and EL_DIM_LAG=Dim, EL_DIM_EUL=Dim+1
444  get_dim_helper(problem_pt,mesh_pt,external_face_mesh_pt);
445 
446  if(Dim > 2)
447  {
448  std::ostringstream error_stream;
449  error_stream << "The elements within the two interacting meshes have a\n"
450  << "dimension not equal to 1 or 2.\n"
451  << "The multi-domain method will not work in this case.\n"
452  << "The dimension is: " << Dim << "\n";
453  throw OomphLibError
454  (error_stream.str(),
455  OOMPH_CURRENT_FUNCTION,
456  OOMPH_EXCEPTION_LOCATION);
457  }
458 
459  // Call the function that actually does the work
461  <EXT_ELEMENT,FACE_ELEMENT_GEOM_OBJECT>
462  (problem_pt,mesh_pt,external_mesh_pt,
463  interaction_index,external_face_mesh_pt);
464  }
465 
466 
467 //========================================================================
468  /// Function to set up the one-way multi-domain interaction for
469  /// FSI-like problems.
470  /// - \c mesh_pt points to the mesh of \c ElemenWithExternalElements for which
471  /// the interaction is set up. In an FSI example, this mesh would contain
472  /// the \c FSIWallElements (either beam/shell elements or the
473  /// \c FSISolidTractionElements that apply the traction to
474  /// a "bulk" solid mesh that is loaded by the fluid.)
475  /// - \c external_mesh_pt points to the mesh that contains the elements
476  /// of type EXT_ELEMENT that provide the "source" for the
477  /// \c ElementWithExternalElements. In an FSI example, this
478  /// mesh would contain the "bulk" fluid elements.
479  /// - \c external_face_mesh_pt points to the mesh of \c FaceElements
480  /// attached to the \c external_mesh_pt. The mesh pointed to by
481  /// \c external_face_mesh_pt has the same dimension as \c mesh_pt.
482  /// The elements contained in \c external_face_mesh_pt are of type
483  /// FACE_ELEMENT_GEOM_OBJECT. In an FSI example, these elements
484  /// are usually the \c FaceElementAsGeomObjects (templated by the
485  /// type of the "bulk" fluid elements to which they are attached)
486  /// that define the FSI boundary of the fluid domain.
487  /// - The interaction_index parameter defaults to zero and must otherwise be
488  /// set by the user if there is more than one mesh that provides "external
489  /// elements" for the Mesh pointed to by mesh_pt (e.g. in the case
490  /// when a beam or shell structure is loaded by fluid from both sides.)
491  /// .
492  /// Vector-based version operates simultaneously on the meshes contained
493  /// in the vectors.
494 //========================================================================
495  template<class EXT_ELEMENT,class FACE_ELEMENT_GEOM_OBJECT>
497  (Problem* problem_pt, const Vector<Mesh*>& mesh_pt,
498  Mesh* const &external_mesh_pt, const Vector<Mesh*>& external_face_mesh_pt,
499  const unsigned& interaction_index)
500  {
501 
502  // How many meshes do we have?
503  unsigned n_mesh=mesh_pt.size();
504 
505 #ifdef PARANOID
506  if (external_face_mesh_pt.size()!=n_mesh)
507  {
508  std::ostringstream error_stream;
509  error_stream << "Sizes of external_face_mesh_pt [ "
510  << external_face_mesh_pt.size() << " ] and "
511  << "mesh_pt [ " << n_mesh << " ] don't match.\n";
512  throw OomphLibError
513  (error_stream.str(),
514  OOMPH_CURRENT_FUNCTION,
515  OOMPH_EXCEPTION_LOCATION);
516  }
517 #endif
518 
519  // Bail out?
520  if (n_mesh==0) return;
521 
522  // Bulk elements must be external elements in this case
524 
525  // Call the auxiliary routine with GEOM_OBJECT=FACE_ELEMENT_GEOM_OBJECT
526  // and EL_DIM_LAG=Dim, EL_DIM_EUL=Dim+1. Use first mesh only.
527  get_dim_helper(problem_pt,mesh_pt[0],external_face_mesh_pt[0]);
528 
529 
530 #ifdef PARANOID
531  // Check consitency
532  unsigned old_dim=Dim;
533  for (unsigned i=1;i<n_mesh;i++)
534  {
535  // Set up Dim again
536  get_dim_helper(problem_pt,mesh_pt[i],external_face_mesh_pt[i]);
537 
538  if (Dim!=old_dim)
539  {
540  std::ostringstream error_stream;
541  error_stream
542  << "Inconsistency: Mesh " << i << " has Dim=" << Dim
543  << "while mesh 0 has Dim="<< old_dim << std::endl;
544  throw OomphLibError
545  (error_stream.str(),
546  OOMPH_CURRENT_FUNCTION,
547  OOMPH_EXCEPTION_LOCATION);
548  }
549  }
550 #endif
551 
552  if(Dim > 2)
553  {
554  std::ostringstream error_stream;
555  error_stream << "The elements within the two interacting meshes have a\n"
556  << "dimension not equal to 1 or 2.\n"
557  << "The multi-domain method will not work in this case.\n"
558  << "The dimension is: " << Dim << "\n";
559  throw OomphLibError
560  (error_stream.str(),
561  OOMPH_CURRENT_FUNCTION,
562  OOMPH_EXCEPTION_LOCATION);
563  }
564 
565  // Now do the actual work for all meshes simultaneously
567  <EXT_ELEMENT,FACE_ELEMENT_GEOM_OBJECT>
568  (problem_pt,mesh_pt,external_mesh_pt,
569  interaction_index,external_face_mesh_pt);
570 
571  }
572 
573 
574 
575 
576 //========================================================================
577 /// This routine calls the locate_zeta routine (simultaneously on each
578 /// processor for each individual processor's element set if necessary)
579 /// and sets up the external (halo) element and node storage as
580 /// necessary. The locate_zeta procedure here works for all multi-domain
581 /// problems where either two meshes occupy the same physical space but have
582 /// differing element types (e.g. a Boussinesq convection problem where
583 /// AdvectionDiffusion elements interact with Navier-Stokes type elements)
584 /// or two meshes interact along some boundary of the external mesh,
585 /// represented by a "face mesh", such as an FSI problem.
586 //========================================================================
587  template<class EXT_ELEMENT,class GEOM_OBJECT>
589  (Problem* problem_pt, Mesh* const &mesh_pt, Mesh* const &external_mesh_pt,
590  const unsigned& interaction_index, Mesh* const &external_face_mesh_pt)
591  {
592 
593 #ifdef USE_VECTOR_BASED_MD
594 
595  // Convert to vector-based storage
596  Vector<Mesh*> mesh_pt_vector(1);
597  mesh_pt_vector[0]=mesh_pt;
598  Vector<Mesh*> external_face_mesh_pt_vector(1);
599  external_face_mesh_pt_vector[0]=external_face_mesh_pt;
600 
601  // Call vector-based version
602  aux_setup_multi_domain_interaction<EXT_ELEMENT,GEOM_OBJECT>
603  (problem_pt, mesh_pt_vector,
604  external_mesh_pt, interaction_index,
605  external_face_mesh_pt_vector);
606 
607 #else
608 
609  //oomph_info << "OLD aux_setup_multi_domain_interaction\n";
610 
611  //Multi-domain setup will not work for elements with
612  //nonuniformly spaced nodes because the same storage is used for
613  //the missing master nodes as for the external halo/haloed nodes
614 #ifdef PARANOID
615  //Must check type of elements in the mesh and in the external mesh
616  //(assume element type is the same for all elements in each mesh)
617 
618  //Get pointer to first element in each mesh
619  GeneralisedElement* el_pt_0=0;
620  if (mesh_pt->nelement()!=0)
621  {
622  el_pt_0 = mesh_pt->element_pt(0);
623  }
624  GeneralisedElement* ext_el_pt_0=0;
625  if (external_mesh_pt->nelement()!=0)
626  {
627  ext_el_pt_0 = external_mesh_pt->element_pt(0);
628  }
629 
630 
631  //Check they are not spectral elements
632  if(dynamic_cast<SpectralElement*>(el_pt_0)!=0
633  || dynamic_cast<SpectralElement*>(ext_el_pt_0)!=0)
634  {
635  throw OomphLibError(
636  "Multi-domain setup does not work with spectral elements.",
637  OOMPH_CURRENT_FUNCTION,
638  OOMPH_EXCEPTION_LOCATION);
639  }
640 
641  //Check they are not hp-refineable elements
642  if(dynamic_cast<PRefineableElement*>(el_pt_0)!=0
643  || dynamic_cast<PRefineableElement*>(ext_el_pt_0)!=0)
644  {
645  throw OomphLibError(
646  "Multi-domain setup does not work with hp-refineable elements.",
647  OOMPH_CURRENT_FUNCTION,
648  OOMPH_EXCEPTION_LOCATION);
649  }
650 #endif
651 
652 #ifdef OOMPH_HAS_MPI
653  // Storage for number of processors and my rank
654  int n_proc=problem_pt->communicator_pt()->nproc();
655  int my_rank=problem_pt->communicator_pt()->my_rank();
656 #endif
657 
658  // Timing
659  double t_start=0.0; double t_end=0.0; double t_local=0.0;
660  double t_set=0.0; double t_locate=0.0; double t_spiral_start=0.0;
661 #ifdef OOMPH_HAS_MPI
662  double t_loop_start=0.0;
663  double t_sendrecv=0.0;
664  double t_missing=0.0;
665  double t_send_info=0.0; double t_create_halo=0.0;
666 #endif
667 
668  if (Doc_timings)
669  {
670  t_start=TimingHelpers::timer();
671  }
672 
673  // Initialize number of zeta coordinates not found yet
674  unsigned n_zeta_not_found=0;
675 
676  // Geometric object used to represent the external (face) mesh
677  MeshAsGeomObject* mesh_geom_obj_pt=0;
678 
679  // Are bulk elements used as external elements?
681  {
682  // Set the geometric object from the external mesh
683  mesh_geom_obj_pt=new MeshAsGeomObject
684  (external_mesh_pt,problem_pt->communicator_pt(),
686  }
687  else
688  {
689  // Set the geometric object from the external face mesh argument
690  mesh_geom_obj_pt=new MeshAsGeomObject
691  (external_face_mesh_pt,problem_pt->communicator_pt(),
693  }
694 
695  unsigned EL_DIM_LAG = mesh_geom_obj_pt->nlagrangian();
696 
697  double t_setup_lookups=0.0;
698  if (Doc_timings)
699  {
700  t_set=TimingHelpers::timer();
701  oomph_info << "CPU for creation of MeshAsGeomObject and bin structure: "
702  << t_set-t_start << std::endl;
703  t_setup_lookups=TimingHelpers::timer();
704  }
705 
706  // Total number of integration points
707  unsigned tot_int=0;
708 
709  // Loop over (this processor's) elements and set lookup array
710  unsigned n_element=mesh_pt->nelement();
711  External_element_located.resize(n_element);
712  for (unsigned e=0;e<n_element;e++)
713  {
714  // Zero-sized vector means its a halo
715  External_element_located[e].resize(0);
716 
718  dynamic_cast<ElementWithExternalElement*>(mesh_pt->element_pt(e));
719 
720 #ifdef OOMPH_HAS_MPI
721  // We're not setting up external elements for halo elements
722  if (!el_pt->is_halo())
723 #endif
724  {
725  //We need to allocate storage for the external elements
726  //within the element. Memory will actually only be
727  //allocated the first time this function is called for
728  //each element, or if the number of interactions or integration
729  //points within the element has changed.
731 
732  unsigned n_intpt=el_pt->integral_pt()->nweight();
733  External_element_located[e].resize(n_intpt);
734  for (unsigned ipt=0;ipt<n_intpt;ipt++)
735  {
736  External_element_located[e][ipt]=0;
737  tot_int++;
738  }
739  }
740  }
741 
742  if (Doc_timings)
743  {
744  double t=TimingHelpers::timer();
745  oomph_info
746  << "CPU for setup of lookup schemes for located elements/coords: "
747  << t-t_setup_lookups << std::endl;
748  }
749 
750  // Find maximum spiral level within the cartesian bin structure
751  unsigned n_max_level=Nx_bin;
752  if (EL_DIM_LAG>=2)
753  {
754  if (Ny_bin > n_max_level)
755  {
756  n_max_level=Ny_bin;
757  }
758  }
759  if (EL_DIM_LAG==3)
760  {
761  if (Nz_bin > n_max_level)
762  {
763  n_max_level=Nz_bin;
764  }
765  }
766 
767  // Storage for info about coordinate location -- initialise to 100%.
768  // Gets overwritten below if we've actually had any location to do here.
769  Vector<double> percentage_coords_located_locally(n_max_level,100.0);
770  Vector<double> percentage_coords_located_elsewhere(n_max_level,100.0);
771 #ifdef OOMPH_HAS_MPI
772  unsigned max_level_reached;
773 #endif
774  // Initialise spiral levels
775  mesh_geom_obj_pt->current_min_spiral_level()=0;
776  mesh_geom_obj_pt->current_max_spiral_level()=N_spiral_chunk-1;
777 
778  // Limit it so that we search at least once
779  if (mesh_geom_obj_pt->current_max_spiral_level()>n_max_level)
780  {
781  mesh_geom_obj_pt->current_max_spiral_level()=n_max_level-1;
782  }
783 
784  // Loop over "spirals/levels" away from the current position
785  unsigned i_level=0;
786  while (mesh_geom_obj_pt->current_max_spiral_level()<n_max_level)
787  {
788  // Record time at start of spiral loop
789  if (Doc_timings)
790  {
791  t_spiral_start=TimingHelpers::timer();
792  // Initialize the total time for sorting elements in bins
793  if (Sort_bin_entries)
794  {
796  }
797  }
798 
799  // Perform locate_zeta locally first! This looks locally for
800  // all not-yet-located zetas within the current spiral range
801  locate_zeta_for_local_coordinates(mesh_pt,external_mesh_pt,
802  mesh_geom_obj_pt,interaction_index);
803 
804  // Store stats about successful locates for reporting later
805  if (Doc_stats)
806  {
807  unsigned count_locates=0;
808  unsigned n_ext_loc=External_element_located.size();
809  for (unsigned e=0;e<n_ext_loc;e++)
810  {
811  unsigned n_intpt=External_element_located[e].size();
812  for (unsigned ipt=0;ipt<n_intpt;ipt++)
813  {
814  count_locates+=External_element_located[e][ipt];
815  }
816  }
817 
818  // Store percentage of integration points successfully located.
819  // Only assign if we had anything to allocte, otherwise 100%
820  // (default assignment; see above) is correct
821  if (tot_int!=0)
822  {
823  percentage_coords_located_locally[i_level]=
824  100.0*double(count_locates)/double(tot_int);
825  }
826  }
827 
828  if (Doc_timings)
829  {
830  t_local=TimingHelpers::timer();
831  oomph_info
832  << "CPU for local location of zeta coordinate [spiral level "
833  << i_level << "]: "
834  << t_local-t_spiral_start << std::endl;
835  }
836 
837 
838  // Now test whether anything needs to be broadcast elsewhere
839  // (i.e. were there any failures in the locate method above?)
840  // If there are, then the zetas for these failures need to be
841  // broadcast...
842 
843  // How many zetas have we failed to find?
844  n_zeta_not_found=Flat_packed_zetas_not_found_locally.size();
845 
846 #ifdef OOMPH_HAS_MPI
847  // Only perform the reduction operation if there's more than one process
848  if (problem_pt->communicator_pt()->nproc() > 1)
849  {
850  unsigned count_local_zetas=n_zeta_not_found;
851  MPI_Allreduce(&count_local_zetas,&n_zeta_not_found,1,MPI_UNSIGNED,MPI_SUM,
852  problem_pt->communicator_pt()->mpi_comm());
853  }
854 
855  // If we have missing zetas on any process
856  // and the problem is distributed, we need to locate elsewhere
857  if ((n_zeta_not_found!=0) && (problem_pt->problem_has_been_distributed()))
858  {
859  // Timings
860  double t_sendrecv_min= DBL_MAX;
861  double t_sendrecv_max=-DBL_MAX;
862  double t_sendrecv_tot=0.0;
863 
864  double t_missing_min= DBL_MAX;
865  double t_missing_max=-DBL_MAX;
866  double t_missing_tot=0.0;
867 
868  double t_send_info_min= DBL_MAX;
869  double t_send_info_max=-DBL_MAX;
870  double t_send_info_tot=0.0;
871 
872  double t_create_halo_min= DBL_MAX;
873  double t_create_halo_max=-DBL_MAX;
874  double t_create_halo_tot=0.0;
875 
876  // Start ring communication: Loop (number of processes - 1)
877  // starting from 1. The variable iproc represents the "distance" from
878  // the current process to the process for which it is attempting
879  // to locate an element for the current set of not-yet-located
880  // zeta coordinates
881  unsigned ring_count=0;
882  for (int iproc=1;iproc<n_proc;iproc++)
883  {
884  // Record time at start of loop
885  if (Doc_timings)
886  {
887  t_loop_start=TimingHelpers::timer();
888  }
889 
890  // Send the zeta values you haven't found to the
891  // next process, receive from the previous process
892  send_and_receive_missing_zetas(problem_pt);
893 
894  if (Doc_timings)
895  {
896  ring_count++;
897  t_sendrecv=TimingHelpers::timer();
898  t_sendrecv_max=std::max(t_sendrecv_max,t_sendrecv-t_loop_start);
899  t_sendrecv_min=std::min(t_sendrecv_min,t_sendrecv-t_loop_start);
900  t_sendrecv_tot+=(t_sendrecv-t_loop_start);
901  }
902 
903  // Perform the locate_zeta for the new set of zetas on this process
905  (iproc,external_mesh_pt,problem_pt,mesh_geom_obj_pt);
906 
907  if (Doc_timings)
908  {
909  t_missing=TimingHelpers::timer();
910  t_missing_max=std::max(t_missing_max,t_missing-t_sendrecv);
911  t_missing_min=std::min(t_missing_min,t_missing-t_sendrecv);
912  t_missing_tot+=(t_missing-t_sendrecv);
913  }
914 
915  // Send any located coordinates back to the correct process, and
916  // prepare to send on to the next process if necessary
917  send_and_receive_located_info(iproc,external_mesh_pt,problem_pt);
918 
919  if (Doc_timings)
920  {
921  t_send_info=TimingHelpers::timer();
922  t_send_info_max=std::max(t_send_info_max,t_send_info-t_missing);
923  t_send_info_min=std::min(t_send_info_min,t_send_info-t_missing);
924  t_send_info_tot+=(t_send_info-t_missing);
925  }
926 
927  // Create any located external halo elements on the current process
928  create_external_halo_elements<EXT_ELEMENT>
929  (iproc,mesh_pt,external_mesh_pt,problem_pt,interaction_index);
930 
931  if (Doc_timings)
932  {
933  t_create_halo=TimingHelpers::timer();
934  t_create_halo_max=std::max(t_create_halo_max,
935  t_create_halo-t_send_info);
936  t_create_halo_min=std::min(t_create_halo_min,
937  t_create_halo-t_send_info);
938  t_create_halo_tot+=(t_create_halo-t_send_info);
939  }
940 
941 
942  // Do we have any further locating to do or have we found
943  // everything at this level of the ring communication?
944  // Only perform the reduction operation if there's more than
945  // one process
946  n_zeta_not_found=Flat_packed_zetas_not_found_locally.size();
947 
948 #ifdef OOMPH_HAS_MPI
949  if (problem_pt->communicator_pt()->nproc() > 1)
950  {
951  unsigned count_local_zetas=n_zeta_not_found;
952  MPI_Allreduce(&count_local_zetas,&n_zeta_not_found,1,
953  MPI_UNSIGNED,MPI_SUM,
954  problem_pt->communicator_pt()->mpi_comm());
955  }
956 #endif
957 
958  // If its is now zero then break out of the ring comms loop
959  if (n_zeta_not_found==0)
960  {
961  if (Doc_timings)
962  {
963  t_local=TimingHelpers::timer();
964  oomph_info
965  << "BREAK N-1: CPU for entrire spiral [spiral level "
966  << i_level << "]: "
967  << t_local-t_spiral_start << std::endl;
968  }
969  break;
970  }
971  }
972 
973 
974  // Doc timings
975  if (Doc_timings)
976  {
977  oomph_info
978  << "Ring-based search continued until iteration "
979  << ring_count << " out of a maximum of "
980  << problem_pt->communicator_pt()->nproc()-1 << "\n";
981  oomph_info
982  << "Total, av, max, min CPU for send/recv of remaining zeta coordinates: "
983  << t_sendrecv_tot << " "
984  << t_sendrecv_tot/double(ring_count) << " "
985  << t_sendrecv_max << " "
986  << t_sendrecv_min << "\n";
987  oomph_info
988  << "Total, av, max, min CPU for location of missing zeta coordinates : "
989  << t_missing_tot << " "
990  << t_missing_tot/double(ring_count) << " "
991  << t_missing_max << " "
992  << t_missing_min << "\n";
993  oomph_info
994  << "Total, av, max, min CPU for send/recv of new element info : "
995  << t_send_info_tot << " "
996  << t_send_info_tot/double(ring_count) << " "
997  << t_send_info_max << " "
998  << t_send_info_min << "\n";
999  oomph_info
1000  << "Total, av, max, min CPU for local creation of external halo objects: "
1001  << t_create_halo_tot << " "
1002  << t_create_halo_tot/double(ring_count) << " "
1003  << t_create_halo_max << " "
1004  << t_create_halo_min << "\n";
1005  }
1006 
1007  } // end if for missing zetas on any process
1008 #endif
1009 
1010  // Store information about location of elements for integration points
1011  if (Doc_stats)
1012  {
1013  unsigned count_locates=0;
1014  unsigned n_ext_loc=External_element_located.size();
1015  for (unsigned e=0;e<n_ext_loc;e++)
1016  {
1017  unsigned n_intpt=External_element_located[e].size();
1018  for (unsigned ipt=0;ipt<n_intpt;ipt++)
1019  {
1020  count_locates+=External_element_located[e][ipt];
1021  }
1022  }
1023 
1024  // Store total percentage of locates so far.
1025  // Only assign if we had anything to allocte, otherwise 100%
1026  // (default assignment; see above) is correct
1027  if (tot_int!=0)
1028  {
1029  percentage_coords_located_elsewhere[i_level]=
1030  100.0*double(count_locates)/double(tot_int);
1031  }
1032  }
1033 
1034  // Do we have any further locating to do? If so, the remaining
1035  // zetas will (hopefully) be found at the next spiral level.
1036  // Only perform the reduction operation if there's more than one process
1037  n_zeta_not_found=Flat_packed_zetas_not_found_locally.size();
1038 
1039 #ifdef OOMPH_HAS_MPI
1040  if (problem_pt->communicator_pt()->nproc() > 1)
1041  {
1042  unsigned count_local_zetas=n_zeta_not_found;
1043  MPI_Allreduce(&count_local_zetas,&n_zeta_not_found,1,MPI_UNSIGNED,MPI_SUM,
1044  problem_pt->communicator_pt()->mpi_comm());
1045  }
1046 
1047  // Specify max level reached for later loop
1048  max_level_reached=i_level+1;
1049 #endif
1050 
1051  /// If it's is now zero then break out of the spirals loop
1052  if (n_zeta_not_found==0)
1053  {
1054  if (Doc_timings)
1055  {
1056  t_local=TimingHelpers::timer();
1057  oomph_info
1058  << "BREAK N: CPU for entrire spiral [spiral level "
1059  << i_level << "]: "
1060  << t_local-t_spiral_start << std::endl;
1061  }
1062  break;
1063  }
1064 
1065  if (Doc_timings)
1066  {
1067  t_local=TimingHelpers::timer();
1068  oomph_info
1069  << "CPU for entrire spiral [spiral level "
1070  << i_level << "]: "
1071  << t_local-t_spiral_start << std::endl;
1072 
1073  // Document the time for sorting the bins
1074  if (Sort_bin_entries)
1075  {
1076  oomph_info << "CPU for sorting elements in bins [spiral level "
1077  << i_level << "]: "
1079  << std::endl;
1080  }
1081 
1082  }
1083 
1084  // Bump up spiral levels
1085  i_level++;
1086  mesh_geom_obj_pt->current_min_spiral_level()+=N_spiral_chunk;
1087  mesh_geom_obj_pt->current_max_spiral_level()+=N_spiral_chunk;
1088 
1089 
1090  } // end of "spirals" loop
1091 
1092  // If we haven't found all zetas we're dead now...
1093  if (n_zeta_not_found!=0)
1094  {
1095  std::ostringstream error_stream;
1096  error_stream
1097  << "Multi_domain_functions::locate_zeta_for_local_coordinates()"
1098  << "\nhas failed ";
1099 
1100 #ifdef OOMPH_HAS_MPI
1101  if (problem_pt->communicator_pt()->nproc() > 1)
1102  {
1103  error_stream << " on proc: " << problem_pt->communicator_pt()->my_rank()
1104  << std::endl;
1105  }
1106 #endif
1107  error_stream
1108  << "\n\n\nThis is most likely to arise because the two meshes\n"
1109  << "that are to be matched don't overlap perfectly or\n"
1110  << "because the elements are distorted and too small a \n"
1111  << "number of sampling points has been used when setting\n"
1112  << "up the bin structure.\n\n"
1113  << "You should try to increase the value of \n"
1114  << "Multi_domain_functions::Nsample_points from \n"
1115  << "its current value of "
1117 
1118  std::ostringstream modifier;
1119 #ifdef OOMPH_HAS_MPI
1120  if (problem_pt->communicator_pt()->nproc() > 1)
1121  {
1122  modifier << "_proc" << problem_pt->communicator_pt()->my_rank();
1123  }
1124 #endif
1125 
1126  std::ofstream outfile;
1127  char filename[100];
1128  sprintf(filename,"missing_coords_mesh%s.dat",modifier.str().c_str());
1129  outfile.open(filename);
1130  unsigned nel=mesh_pt->nelement();
1131  for (unsigned e=0;e<nel;e++)
1132  {
1133  mesh_pt->finite_element_pt(e)->FiniteElement::output(outfile);
1134  }
1135  outfile.close();
1136 
1137  sprintf(filename,"missing_coords_ext_mesh%s.dat",modifier.str().c_str());
1138  outfile.open(filename);
1139  nel=external_mesh_pt->nelement();
1140  for (unsigned e=0;e<nel;e++)
1141  {
1142  external_mesh_pt->finite_element_pt(e)->FiniteElement::output(outfile);
1143  }
1144  outfile.close();
1145 
1146  sprintf(filename,"missing_coords_bin%s.dat",modifier.str().c_str());
1147  outfile.open(filename);
1148  mesh_geom_obj_pt->output_bins(outfile);
1149  outfile.close();
1150 
1151  sprintf(filename,"missing_coords%s.dat",modifier.str().c_str());
1152  outfile.open(filename);
1153  unsigned n=External_element_located.size();
1154  error_stream << "Number of unlocated elements " << n << std::endl;
1155  for (unsigned e=0;e<n;e++)
1156  {
1157  unsigned n_intpt=External_element_located[e].size();
1158  if (n_intpt==0)
1159  {
1160  error_stream
1161  << "Failure to locate in halo element! Why are we searching there?"
1162  << std::endl;
1163  }
1164  for (unsigned ipt=0;ipt<n_intpt;ipt++)
1165  {
1166  if (External_element_located[e][ipt]==0)
1167  {
1168  error_stream << "Failure at element/intpt: "
1169  << e << " " << ipt << std::endl;
1170 
1171  // Cast
1173  dynamic_cast<ElementWithExternalElement*>(mesh_pt->element_pt(e));
1174 
1175  unsigned n_dim_el=el_pt->dim();
1176  Vector<double> s(n_dim_el);
1177  for (unsigned i=0;i<n_dim_el;i++)
1178  {
1179  s[i]=el_pt->integral_pt()->knot(ipt,i);
1180  }
1181  unsigned n_dim=el_pt->node_pt(0)->ndim();
1182  Vector<double> x(n_dim);
1183  el_pt->interpolated_x(s,x);
1184  for (unsigned i=0;i<n_dim;i++)
1185  {
1186  error_stream << x[i] << " ";
1187  outfile<< x[i] << " ";
1188  }
1189  error_stream << std::endl;
1190  outfile << std::endl;
1191  }
1192  }
1193  }
1194  outfile.close();
1195 
1196  error_stream
1197  << "Mesh and external mesh documented in missing_coords_mesh*.dat\n"
1198  << "and missing_coords_ext_mesh*.dat, respectively. Missing \n"
1199  << "coordinates in missing_coords*.dat\n";
1200  throw OomphLibError
1201  (error_stream.str(),
1202  OOMPH_CURRENT_FUNCTION,
1203  OOMPH_EXCEPTION_LOCATION);
1204  }
1205 
1206 
1207  // Doc timings if required
1208  if (Doc_timings)
1209  {
1210  t_locate=TimingHelpers::timer();
1211  oomph_info
1212  << "Total CPU for location and creation of all external elements: "
1213  << t_locate-t_start << std::endl;
1214  }
1215 
1216  // Delete the geometric object representing the mesh
1217  delete mesh_geom_obj_pt;
1218 
1219  // Clean up all the (extern) Vectors associated with creating the
1220  // external storage information
1221  clean_up();
1222 
1223 #ifdef OOMPH_HAS_MPI
1224  // Output information about external storage if required
1225  if (Doc_stats)
1226  {
1227  // Report stats regarding location method
1228  bool comm_was_required=false;
1229  oomph_info << "-------------------------------------------" << std::endl;
1230  oomph_info << "- Cumulative percentage of locate success -" << std::endl;
1231  oomph_info << "--- Spiral -- Found local -- Found else ---" << std::endl;
1232  for (unsigned level=0; level<max_level_reached; level++)
1233  {
1234  oomph_info << "--- " << level << " -- "
1235  << percentage_coords_located_locally[level] << " -- "
1236  << percentage_coords_located_elsewhere[level] << " ---"
1237  << std::endl;
1238  // Has communication with other processors at this level actually
1239  // produced any results?
1240  if (percentage_coords_located_elsewhere[level]>
1241  percentage_coords_located_locally[level])
1242  {
1243  comm_was_required=true;
1244  }
1245  }
1246  oomph_info << "-------------------------------------------" << std::endl;
1247 
1248  // Initialise to indicate that none of the zetas required
1249  // on this processor were located through parallel ring search,
1250  // i.e. comm was not required and we could have done some
1251  // more local searching first
1252  oomph_info << std::endl;
1253  oomph_info <<"ASSESSMENT OF NEED FOR PARALLEL SEARCH: \n";
1254  oomph_info <<"=======================================\n";
1255  unsigned status=0;
1256  if (comm_was_required)
1257  {
1258  oomph_info
1259  <<"- Ring-based parallel search did successfully locate zetas on proc "
1260  << my_rank << std::endl;
1261  status=1;
1262  }
1263  else
1264  {
1265  if (max_level_reached>1)
1266  {
1267  oomph_info
1268  << "- Ring-based parallel search did NOT locate zetas on proc"
1269  << my_rank << std::endl;
1270  }
1271  else
1272  {
1273  oomph_info
1274  << "- No ring-based parallel search was performed on proc"
1275  << my_rank << std::endl;
1276  }
1277  }
1278 
1279  // Allreduce to check if anyone has benefitted from parallel ring
1280  // search
1281  unsigned overall_status=0;
1282  // Only perform the reduction operation if there's more than one process
1283  if (problem_pt->communicator_pt()->nproc() > 1)
1284  {
1285  MPI_Allreduce(&status,&overall_status,1,
1286  MPI_UNSIGNED,MPI_MAX,
1287  problem_pt->communicator_pt()->mpi_comm());
1288  }
1289  else
1290  {
1291  overall_status=status;
1292  }
1293 
1294  // Report of mpi was useful to anyone
1295  if (overall_status==1)
1296  {
1297  oomph_info << "- Ring-based, parallel search did succesfully\n";
1298  oomph_info << " locate zetas on at least one other proc, so it\n";
1299  oomph_info << " was worthwhile.\n";
1300  oomph_info << std::endl;
1301  }
1302  else
1303  {
1304  if (max_level_reached>1)
1305  {
1306  oomph_info << "- Ring-based, parallel search did NOT locate zetas\n";
1307  oomph_info << " on ANY other procs, i.e it was useless.\n";
1308  oomph_info << " --> We should really have done more local search\n";
1309  oomph_info << " by reducing number of bins, or doing more spirals\n";
1310  oomph_info << " in one go before initiating parallel search.\n";
1311  oomph_info << std::endl;
1312  }
1313  else
1314  {
1315  oomph_info << "- No ring-based, parallel search was performed\n";
1316  oomph_info << " or necessary. Perfect!\n";
1317  oomph_info << std::endl;
1318  }
1319  }
1320 
1321 
1322  // How many external elements does the external mesh have now?
1323  oomph_info << "------------------------------------------" << std::endl;
1324  oomph_info << "External mesh: I have " << external_mesh_pt->nelement()
1325  << " elements, and " << std::endl
1326  << external_mesh_pt->nexternal_halo_element()
1327  << " external halo elements, "
1328  << external_mesh_pt->nexternal_haloed_element()
1329  << " external haloed elements"
1330  << std::endl;
1331 
1332  // How many external nodes does each submesh have now?
1333  oomph_info << "------------------------------------------" << std::endl;
1334  oomph_info << "External mesh: I have " << external_mesh_pt->nnode()
1335  << " nodes, and " << std::endl
1336  << external_mesh_pt->nexternal_halo_node()
1337  << " external halo nodes, "
1338  << external_mesh_pt->nexternal_haloed_node()
1339  << " external haloed nodes"
1340  << std::endl;
1341  oomph_info << "------------------------------------------" << std::endl;
1342  }
1343 
1344  // Output further information about (external) halo(ed)
1345  // elements and nodes if required
1346  if (Doc_full_stats)
1347  {
1348  // How many elements does this submesh have for each of the processors?
1349  for (int iproc=0;iproc<n_proc;iproc++)
1350  {
1351  oomph_info << "----------------------------------------" << std::endl;
1352  oomph_info << "With process " << iproc << " there are "
1353  << external_mesh_pt->nroot_halo_element(iproc)
1354  << " root halo elements, and "
1355  << external_mesh_pt->nroot_haloed_element(iproc)
1356  << " root haloed elements" << std::endl
1357  << "and there are "
1358  << external_mesh_pt->nexternal_halo_element(iproc)
1359  << " external halo elements, and "
1360  << external_mesh_pt->nexternal_haloed_element(iproc)
1361  << " external haloed elements." << std::endl;
1362 
1363  oomph_info << "----------------------------------------" << std::endl;
1364  oomph_info << "With process " << iproc << " there are "
1365  << external_mesh_pt->nhalo_node(iproc)
1366  << " halo nodes, and "
1367  << external_mesh_pt->nhaloed_node(iproc)
1368  << " haloed nodes" << std::endl
1369  << "and there are "
1370  << external_mesh_pt->nexternal_halo_node(iproc)
1371  << " external halo nodes, and "
1372  << external_mesh_pt->nexternal_haloed_node(iproc)
1373  << " external haloed nodes." << std::endl;
1374  }
1375  oomph_info << "-----------------------------------------" << std::endl
1376  << std::endl;
1377  }
1378 
1379  #endif
1380 
1381  // Doc timings if required
1382  if (Doc_timings)
1383  {
1384  t_end=TimingHelpers::timer();
1385  oomph_info << "CPU for (one way) aux_setup_multi_domain_interaction: "
1386  << t_end-t_start << std::endl;
1387  }
1388 
1389 #endif
1390 
1391  } // end of aux_setup_multi_domain_interaction
1392 
1393 
1394 
1395 
1396 
1397 //========================================================================
1398 /// This routine calls the locate_zeta routine (simultaneously on each
1399 /// processor for each individual processor's element set if necessary)
1400 /// and sets up the external (halo) element and node storage as
1401 /// necessary. The locate_zeta procedure here works for all multi-domain
1402 /// problems where either two meshes occupy the same physical space but have
1403 /// differing element types (e.g. a Boussinesq convection problem where
1404 /// AdvectionDiffusion elements interact with Navier-Stokes type elements)
1405 /// or two meshes interact along some boundary of the external mesh,
1406 /// represented by a "face mesh", such as an FSI problem.
1407 ///
1408 /// Vector-based version operates simultaneously on the meshes contained
1409 /// in the vectors.
1410 //========================================================================
1411  template<class EXT_ELEMENT,class GEOM_OBJECT>
1413  (Problem* problem_pt, const Vector<Mesh*>& mesh_pt,
1414  Mesh* const &external_mesh_pt, const unsigned& interaction_index,
1415  const Vector<Mesh*>& external_face_mesh_pt)
1416  {
1417 
1418  // oomph_info << "NEW aux_setup_multi_domain_interaction\n";
1419 
1420  // How many meshes do we have?
1421  unsigned n_mesh=mesh_pt.size();
1422 
1423 #ifdef PARANOID
1424  if (external_face_mesh_pt.size()!=n_mesh)
1425  {
1426  std::ostringstream error_stream;
1427  error_stream << "Sizes of external_face_mesh_pt [ "
1428  << external_face_mesh_pt.size() << " ] and "
1429  << "mesh_pt [ " << n_mesh << " ] don't match.\n";
1430  throw OomphLibError
1431  (error_stream.str(),
1432  OOMPH_CURRENT_FUNCTION,
1433  OOMPH_EXCEPTION_LOCATION);
1434  }
1435 #endif
1436 
1437  // Bail out?
1438  if (n_mesh==0) return;
1439 
1440  //Multi-domain setup will not work for elements with
1441  //nonuniformly spaced nodes
1442  //Must check type of elements in the mesh and in the external mesh
1443  //(assume element type is the same for all elements in each mesh)
1444 
1445 #ifdef PARANOID
1446 
1447  // Pointer to first element in external mesh
1448  GeneralisedElement* ext_el_pt_0=0;
1449  if (external_mesh_pt->nelement()!=0)
1450  {
1451  ext_el_pt_0 = external_mesh_pt->element_pt(0);
1452  }
1453 
1454  // Loop over all meshes
1455  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
1456  {
1457  //Get pointer to first element in each mesh
1458  GeneralisedElement* el_pt_0=0;
1459  if (mesh_pt[i_mesh]->nelement()!=0)
1460  {
1461  el_pt_0 = mesh_pt[i_mesh]->element_pt(0);
1462  }
1463 
1464  //Check they are not spectral elements
1465  if(dynamic_cast<SpectralElement*>(el_pt_0)!=0
1466  || dynamic_cast<SpectralElement*>(ext_el_pt_0)!=0)
1467  {
1468  throw OomphLibError(
1469  "Multi-domain setup does not work with spectral elements.",
1470  OOMPH_CURRENT_FUNCTION,
1471  OOMPH_EXCEPTION_LOCATION);
1472  }
1473 
1474  //Check they are not hp-refineable elements
1475  if(dynamic_cast<PRefineableElement*>(el_pt_0)!=0
1476  || dynamic_cast<PRefineableElement*>(ext_el_pt_0)!=0)
1477  {
1478  throw OomphLibError(
1479  "Multi-domain setup does not work with hp-refineable elements.",
1480  OOMPH_CURRENT_FUNCTION,
1481  OOMPH_EXCEPTION_LOCATION);
1482  }
1483  } // end over initial loop over meshes
1484 
1485 #endif
1486 
1487 
1488 
1489 #ifdef OOMPH_HAS_MPI
1490  // Storage for number of processors and my rank
1491  int n_proc=problem_pt->communicator_pt()->nproc();
1492  int my_rank=problem_pt->communicator_pt()->my_rank();
1493 #endif
1494 
1495  // Timing
1496  double t_start=0.0; double t_end=0.0; double t_local=0.0;
1497  double t_set=0.0; double t_locate=0.0; double t_spiral_start=0.0;
1498 #ifdef OOMPH_HAS_MPI
1499  double t_loop_start=0.0;
1500  double t_sendrecv=0.0;
1501  double t_missing=0.0;
1502  double t_send_info=0.0; double t_create_halo=0.0;
1503 #endif
1504 
1505  if (Doc_timings)
1506  {
1507  t_start=TimingHelpers::timer();
1508  }
1509 
1510  // Initialize number of zeta coordinates not found yet
1511  unsigned n_zeta_not_found=0;
1512 
1513  // Geometric objects used to represent the external (face) meshes
1514  Vector<MeshAsGeomObject*> mesh_geom_obj_pt(n_mesh,0);
1515 
1516  // Initialise lagrangian dimension of element
1517  unsigned EL_DIM_LAG = 0;
1518 
1519  // Create mesh as geom objects for all meshes
1520  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
1521  {
1522  // Are bulk elements used as external elements?
1524  {
1525  // Fix this when required
1526  if (n_mesh!=1)
1527  {
1528  std::ostringstream error_stream;
1529  error_stream
1530  << "Sorry I currently can't deal with non-bulk external elements\n"
1531  << "in multi-domain setup for multiple meshes.\n"
1532  << "The functionality should be easy to implement now that you\n"
1533  << "have a test case. If you're not willinig to do this, call\n"
1534  << "the multi-domain setup mesh-by-mesh instead (though this can\n"
1535  << "be costly in parallel because of global comms. \n";
1536  throw OomphLibError
1537  (error_stream.str(),
1538  OOMPH_CURRENT_FUNCTION,
1539  OOMPH_EXCEPTION_LOCATION);
1540  }
1541 
1542  // Set the geometric object from the external mesh
1543  mesh_geom_obj_pt[0]=new MeshAsGeomObject
1544  (external_mesh_pt,problem_pt->communicator_pt(),
1546  }
1547  else
1548  {
1549  // Set the geometric object from the external face mesh argument
1550  mesh_geom_obj_pt[i_mesh]=new MeshAsGeomObject
1551  (external_face_mesh_pt[i_mesh],problem_pt->communicator_pt(),
1553  }
1554 
1555 #ifdef PARANOID
1556  unsigned old_el_dim_lag=EL_DIM_LAG;
1557 #endif
1558 
1559  // Set lagrangian dimension of element
1560  EL_DIM_LAG = mesh_geom_obj_pt[i_mesh]->nlagrangian();
1561 
1562 #ifdef PARANOID
1563  // Check consistency
1564  if (i_mesh>0)
1565  {
1566  if (EL_DIM_LAG!=old_el_dim_lag)
1567  {
1568  std::ostringstream error_stream;
1569  error_stream << "Lagrangian dimensions of elements don't match \n "
1570  << "between meshes: " << EL_DIM_LAG << " "
1571  << old_el_dim_lag << "\n";
1572  throw OomphLibError
1573  (error_stream.str(),
1574  OOMPH_CURRENT_FUNCTION,
1575  OOMPH_EXCEPTION_LOCATION);
1576  }
1577  }
1578 #endif
1579 
1580 
1581  }// end of loop over meshes
1582 
1583  double t_setup_lookups=0.0;
1584  if (Doc_timings)
1585  {
1586  t_set=TimingHelpers::timer();
1587  oomph_info << "CPU for creation of MeshAsGeomObjects and bin structure: "
1588  << t_set-t_start << std::endl;
1589  t_setup_lookups=TimingHelpers::timer();
1590  }
1591 
1592  // Total number of integration points
1593  unsigned tot_int=0;
1594 
1595  // Counter for total number of elements (in flat packed order)
1596  unsigned e_count=0;
1597 
1598  // Loop over all meshes to get total number of elements
1599  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
1600  {
1601  e_count+=mesh_pt[i_mesh]->nelement();
1602  }
1603  External_element_located.resize(e_count);
1604 
1605  // Reset counter for elements in flat packed storage
1606  e_count=0;
1607 
1608  // Loop over all meshes
1609  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
1610  {
1611  // Loop over (this processor's) elements and set lookup array
1612  unsigned n_element=mesh_pt[i_mesh]->nelement();
1613  for (unsigned e=0;e<n_element;e++)
1614  {
1615  // Zero-sized vector means its a halo
1616  External_element_located[e_count].resize(0);
1618  dynamic_cast<ElementWithExternalElement*>(
1619  mesh_pt[i_mesh]->element_pt(e));
1620 
1621 #ifdef OOMPH_HAS_MPI
1622  // We're not setting up external elements for halo elements
1623  if (!el_pt->is_halo())
1624 #endif
1625  {
1626  //We need to allocate storage for the external elements
1627  //within the element. Memory will actually only be
1628  //allocated the first time this function is called for
1629  //each element, or if the number of interactions or integration
1630  //points within the element has changed.
1632 
1633  // Clear any previous allocation
1634  unsigned n_intpt=el_pt->integral_pt()->nweight();
1635  for (unsigned ipt=0;ipt<n_intpt;ipt++)
1636  {
1637  el_pt->external_element_pt(interaction_index,ipt)=0;
1638  }
1639 
1640  External_element_located[e_count].resize(n_intpt);
1641  for (unsigned ipt=0;ipt<n_intpt;ipt++)
1642  {
1643  External_element_located[e_count][ipt]=0;
1644  tot_int++;
1645  }
1646  }
1647  // next element
1648  e_count++;
1649  }
1650  } // end of loop over meshes
1651 
1652  if (Doc_timings)
1653  {
1654  double t=TimingHelpers::timer();
1655  oomph_info
1656  << "CPU for setup of lookup schemes for located elements/coords: "
1657  << t-t_setup_lookups << std::endl;
1658  }
1659 
1660  // Find maximum spiral level within the cartesian bin structure
1661  unsigned n_max_level=Nx_bin;
1662  if (EL_DIM_LAG>=2)
1663  {
1664  if (Ny_bin > n_max_level)
1665  {
1666  n_max_level=Ny_bin;
1667  }
1668  }
1669  if (EL_DIM_LAG==3)
1670  {
1671  if (Nz_bin > n_max_level)
1672  {
1673  n_max_level=Nz_bin;
1674  }
1675  }
1676 
1677  // Storage for info about coordinate location -- initialise to 100%.
1678  // Gets overwritten below if we've actually had any location to do here.
1679  Vector<double> percentage_coords_located_locally(n_max_level,100.0);
1680  Vector<double> percentage_coords_located_elsewhere(n_max_level,100.0);
1681 #ifdef OOMPH_HAS_MPI
1682  unsigned max_level_reached=1;
1683 #endif
1684 
1685  // Loop over all meshes
1686  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
1687  {
1688  // Initialise spiral levels
1689  mesh_geom_obj_pt[i_mesh]->current_min_spiral_level()=0;
1690  mesh_geom_obj_pt[i_mesh]->current_max_spiral_level()=N_spiral_chunk-1;
1691 
1692  // Limit it so that we search at least once
1693  if (mesh_geom_obj_pt[i_mesh]->current_max_spiral_level()>n_max_level)
1694  {
1695  mesh_geom_obj_pt[i_mesh]->current_max_spiral_level()=n_max_level-1;
1696  }
1697  }
1698 
1699  // Loop over "spirals/levels" away from the current position
1700  // Note: All meshes go through their spirals simultaneously;
1701  // read out spiral level from first one
1702  unsigned i_level=0;
1703  while (mesh_geom_obj_pt[0]->current_max_spiral_level()<n_max_level)
1704  {
1705 
1706  // Record time at start of spiral loop
1707  if (Doc_timings)
1708  {
1709  t_spiral_start=TimingHelpers::timer();
1710  // Initialize the total time for sorting elements in bins
1711  if (Sort_bin_entries)
1712  {
1714  }
1715 
1716  }
1717 
1718  // Perform locate_zeta locally first! This looks locally for
1719  // all not-yet-located zetas within the current spiral range.
1720  locate_zeta_for_local_coordinates(mesh_pt,external_mesh_pt,
1721  mesh_geom_obj_pt,
1722  interaction_index);
1723 
1724  // Store stats about successful locates for reporting later
1725  if (Doc_stats)
1726  {
1727  unsigned count_locates=0;
1728  unsigned n_ext_loc=External_element_located.size();
1729  for (unsigned e=0;e<n_ext_loc;e++)
1730  {
1731  unsigned n_intpt=External_element_located[e].size();
1732  for (unsigned ipt=0;ipt<n_intpt;ipt++)
1733  {
1734  count_locates+=External_element_located[e][ipt];
1735  }
1736  }
1737 
1738  // Store percentage of integration points successfully located.
1739  // Only assign if we had anything to allocte, otherwise 100%
1740  // (default assignment; see above) is correct
1741  if (tot_int!=0)
1742  {
1743  percentage_coords_located_locally[i_level]=
1744  100.0*double(count_locates)/double(tot_int);
1745  }
1746  }
1747 
1748  if (Doc_timings)
1749  {
1750  t_local=TimingHelpers::timer();
1751  oomph_info
1752  << "CPU for local location of zeta coordinate [spiral level "
1753  << i_level << "]: "
1754  << t_local-t_spiral_start << std::endl;
1755  }
1756 
1757 
1758  // Now test whether anything needs to be broadcast elsewhere
1759  // (i.e. were there any failures in the locate method above?)
1760  // If there are, then the zetas for these failures need to be
1761  // broadcast...
1762 
1763  // How many zetas have we failed to find? [Note: Array is padded
1764  // by Dim padded entries (DBL_MAX) for each mesh]
1765  n_zeta_not_found=Flat_packed_zetas_not_found_locally.size()-
1766  Dim*n_mesh;
1767 
1768 #ifdef OOMPH_HAS_MPI
1769  // Only perform the reduction operation if there's more than one process
1770  if (problem_pt->communicator_pt()->nproc() > 1)
1771  {
1772  unsigned count_local_zetas=n_zeta_not_found;
1773  MPI_Allreduce(&count_local_zetas,&n_zeta_not_found,1,
1774  MPI_UNSIGNED,MPI_SUM,
1775  problem_pt->communicator_pt()->mpi_comm());
1776  }
1777 
1778  // If we have missing zetas on any process
1779  // and the problem is distributed, we need to locate elsewhere
1780  if ((n_zeta_not_found!=0) && (problem_pt->problem_has_been_distributed()))
1781  {
1782  // Timings
1783  double t_sendrecv_min= DBL_MAX;
1784  double t_sendrecv_max=-DBL_MAX;
1785  double t_sendrecv_tot=0.0;
1786 
1787  double t_missing_min= DBL_MAX;
1788  double t_missing_max=-DBL_MAX;
1789  double t_missing_tot=0.0;
1790 
1791  double t_send_info_min= DBL_MAX;
1792  double t_send_info_max=-DBL_MAX;
1793  double t_send_info_tot=0.0;
1794 
1795  double t_create_halo_min= DBL_MAX;
1796  double t_create_halo_max=-DBL_MAX;
1797  double t_create_halo_tot=0.0;
1798 
1799  // Start ring communication: Loop (number of processes - 1)
1800  // starting from 1. The variable iproc represents the "distance" from
1801  // the current process to the process for which it is attempting
1802  // to locate an element for the current set of not-yet-located
1803  // zeta coordinates
1804  unsigned ring_count=0;
1805  for (int iproc=1;iproc<n_proc;iproc++)
1806  {
1807  // Record time at start of loop
1808  if (Doc_timings)
1809  {
1810  t_loop_start=TimingHelpers::timer();
1811  }
1812 
1813  // Send the zeta values you haven't found to the
1814  // next process, receive from the previous process:
1815  // (Padded) Flat_packed_zetas_not_found_locally are sent
1816  // to next processor where they are received as
1817  // (padded) Received_flat_packed_zetas_to_be_found.
1818  send_and_receive_missing_zetas(problem_pt);
1819 
1820  if (Doc_timings)
1821  {
1822  ring_count++;
1823  t_sendrecv=TimingHelpers::timer();
1824  t_sendrecv_max=std::max(t_sendrecv_max,t_sendrecv-t_loop_start);
1825  t_sendrecv_min=std::min(t_sendrecv_min,t_sendrecv-t_loop_start);
1826  t_sendrecv_tot+=(t_sendrecv-t_loop_start);
1827  }
1828 
1829  // Perform the locate_zeta for the new set of zetas on this process
1831  (iproc,external_mesh_pt,problem_pt,mesh_geom_obj_pt);
1832 
1833  if (Doc_timings)
1834  {
1835  t_missing=TimingHelpers::timer();
1836  t_missing_max=std::max(t_missing_max,t_missing-t_sendrecv);
1837  t_missing_min=std::min(t_missing_min,t_missing-t_sendrecv);
1838  t_missing_tot+=(t_missing-t_sendrecv);
1839  }
1840 
1841  // Send any located coordinates back to the correct process, and
1842  // prepare to send on to the next process if necessary
1843  send_and_receive_located_info(iproc,external_mesh_pt,problem_pt);
1844 
1845  if (Doc_timings)
1846  {
1847  t_send_info=TimingHelpers::timer();
1848  t_send_info_max=std::max(t_send_info_max,t_send_info-t_missing);
1849  t_send_info_min=std::min(t_send_info_min,t_send_info-t_missing);
1850  t_send_info_tot+=(t_send_info-t_missing);
1851  }
1852 
1853  // Create any located external halo elements on the current process
1854  create_external_halo_elements<EXT_ELEMENT>
1855  (iproc,mesh_pt,external_mesh_pt,problem_pt,interaction_index);
1856 
1857  if (Doc_timings)
1858  {
1859  t_create_halo=TimingHelpers::timer();
1860  t_create_halo_max=std::max(t_create_halo_max,
1861  t_create_halo-t_send_info);
1862  t_create_halo_min=std::min(t_create_halo_min,
1863  t_create_halo-t_send_info);
1864  t_create_halo_tot+=(t_create_halo-t_send_info);
1865  }
1866 
1867  // Do we have any further locating to do or have we found
1868  // everything at this level of the ring communication?
1869  // Only perform the reduction operation if there's more than
1870  // one process [Note: Array is padded
1871  // by DIM times DBL_MAX entries for each mesh]
1872  n_zeta_not_found=Flat_packed_zetas_not_found_locally.size()-
1873  Dim*n_mesh;
1874 
1875 
1876 #ifdef OOMPH_HAS_MPI
1877  if (problem_pt->communicator_pt()->nproc() > 1)
1878  {
1879  unsigned count_local_zetas=n_zeta_not_found;
1880  MPI_Allreduce(&count_local_zetas,&n_zeta_not_found,1,
1881  MPI_UNSIGNED,MPI_SUM,
1882  problem_pt->communicator_pt()->mpi_comm());
1883  }
1884 #endif
1885 
1886  // If its is now zero then break out of the ring comms loop
1887  if (n_zeta_not_found==0)
1888  {
1889  if (Doc_timings)
1890  {
1891  t_local=TimingHelpers::timer();
1892  oomph_info
1893  << "BREAK N-1: CPU for entrire spiral [spiral level "
1894  << i_level << "]: "
1895  << t_local-t_spiral_start << std::endl;
1896  }
1897  break;
1898  }
1899  }
1900 
1901 
1902  // Doc timings
1903  if (Doc_timings)
1904  {
1905  oomph_info
1906  << "Ring-based search continued until iteration "
1907  << ring_count << " out of a maximum of "
1908  << problem_pt->communicator_pt()->nproc()-1 << "\n";
1909  oomph_info
1910  << "Total, av, max, min CPU for send/recv of remaining zeta coordinates: "
1911  << t_sendrecv_tot << " "
1912  << t_sendrecv_tot/double(ring_count) << " "
1913  << t_sendrecv_max << " "
1914  << t_sendrecv_min << "\n";
1915  oomph_info
1916  << "Total, av, max, min CPU for location of missing zeta coordinates : "
1917  << t_missing_tot << " "
1918  << t_missing_tot/double(ring_count) << " "
1919  << t_missing_max << " "
1920  << t_missing_min << "\n";
1921  oomph_info
1922  << "Total, av, max, min CPU for send/recv of new element info : "
1923  << t_send_info_tot << " "
1924  << t_send_info_tot/double(ring_count) << " "
1925  << t_send_info_max << " "
1926  << t_send_info_min << "\n";
1927  oomph_info
1928  << "Total, av, max, min CPU for local creation of external halo objects: "
1929  << t_create_halo_tot << " "
1930  << t_create_halo_tot/double(ring_count) << " "
1931  << t_create_halo_max << " "
1932  << t_create_halo_min << "\n";
1933  }
1934 
1935  } // end if for missing zetas on any process
1936 #endif
1937 
1938 
1939  // Store information about location of elements for integration points
1940  if (Doc_stats)
1941  {
1942  unsigned count_locates=0;
1943  unsigned n_ext_loc=External_element_located.size();
1944  for (unsigned e=0;e<n_ext_loc;e++)
1945  {
1946  unsigned n_intpt=External_element_located[e].size();
1947  for (unsigned ipt=0;ipt<n_intpt;ipt++)
1948  {
1949  count_locates+=External_element_located[e][ipt];
1950  }
1951  }
1952 
1953  // Store total percentage of locates so far.
1954  // Only assign if we had anything to allocte, otherwise 100%
1955  // (default assignment; see above) is correct
1956  if (tot_int!=0)
1957  {
1958  percentage_coords_located_elsewhere[i_level]=
1959  100.0*double(count_locates)/double(tot_int);
1960  }
1961  }
1962 
1963  // Do we have any further locating to do? If so, the remaining
1964  // zetas will (hopefully) be found at the next spiral level.
1965  // Only perform the reduction operation if there's more than one process
1966  // [Note: Array is padded
1967  // by DIM times DBL_MAX entries for each mesh]
1968  n_zeta_not_found=Flat_packed_zetas_not_found_locally.size()-
1969  Dim*n_mesh;
1970 
1971 
1972 #ifdef OOMPH_HAS_MPI
1973  if (problem_pt->communicator_pt()->nproc() > 1)
1974  {
1975  unsigned count_local_zetas=n_zeta_not_found;
1976  MPI_Allreduce(&count_local_zetas,&n_zeta_not_found,1,
1977  MPI_UNSIGNED,MPI_SUM,
1978  problem_pt->communicator_pt()->mpi_comm());
1979  }
1980 
1981  // Specify max level reached for later loop
1982  max_level_reached=i_level+1;
1983 #endif
1984 
1985  /// If it's is now zero then break out of the spirals loop
1986  if (n_zeta_not_found==0)
1987  {
1988  if (Doc_timings)
1989  {
1990  t_local=TimingHelpers::timer();
1991  oomph_info
1992  << "BREAK N: CPU for entrire spiral [spiral level "
1993  << i_level << "]: "
1994  << t_local-t_spiral_start << std::endl;
1995  }
1996  break;
1997  }
1998 
1999  if (Doc_timings)
2000  {
2001  t_local=TimingHelpers::timer();
2002  oomph_info
2003  << "CPU for entrire spiral [spiral level "
2004  << i_level << "]: "
2005  << t_local-t_spiral_start << std::endl;
2006 
2007  // Document the time for sorting the bins
2008  if (Sort_bin_entries)
2009  {
2010  oomph_info << "CPU for sorting elements in bins [spiral level "
2011  << i_level << "]: "
2013  << std::endl;
2014  }
2015 
2016  }
2017 
2018  // Bump up spiral levels for all meshes
2019  i_level++;
2020  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
2021  {
2022  mesh_geom_obj_pt[i_mesh]->current_min_spiral_level()+=N_spiral_chunk;
2023  mesh_geom_obj_pt[i_mesh]->current_max_spiral_level()+=N_spiral_chunk;
2024  }
2025 
2026  } // end of "spirals" loop
2027 
2028 
2029 
2030  // If we haven't found all zetas we're dead now...
2031  //-------------------------------------------------
2032  if (n_zeta_not_found!=0)
2033  {
2034  // Shout?
2036  {
2037 
2038  std::ostringstream error_stream;
2039  error_stream
2040  << "Multi_domain_functions::locate_zeta_for_local_coordinates()"
2041  << "\nhas failed ";
2042 
2043 #ifdef OOMPH_HAS_MPI
2044  if (problem_pt->communicator_pt()->nproc() > 1)
2045  {
2046  error_stream << " on proc: "
2047  << problem_pt->communicator_pt()->my_rank()
2048  << std::endl;
2049  }
2050 #endif
2051  error_stream
2052  << "\n\n\nThis is most likely to arise because the two meshes\n"
2053  << "that are to be matched don't overlap perfectly or\n"
2054  << "because the elements are distorted and too small a \n"
2055  << "number of sampling points has been used when setting\n"
2056  << "up the bin structure.\n\n"
2057  << "You should try to increase the value of \n"
2058  << "Multi_domain_functions::Nsample_points from \n"
2059  << "its current value of "
2061  << "\n\n"
2062  << "NOTE: You can suppress this error and \"accept failure\""
2063  << " by setting the public boolean \n\n"
2064  << " Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_domain_interaction\n\n"
2065  << " to true. In this case, the pointers to external elements\n"
2066  << " that couldn't be located will remain null\n";
2067 
2068  std::ostringstream modifier;
2069 #ifdef OOMPH_HAS_MPI
2070  if (problem_pt->communicator_pt()->nproc() > 1)
2071  {
2072  modifier << "_proc" << problem_pt->communicator_pt()->my_rank();
2073  }
2074 #endif
2075 
2076  // Loop over all meshes
2077  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
2078  {
2079 
2080  // Add yet another modifier to distinguish meshes if reqd
2081  if (n_mesh>1)
2082  {
2083  modifier << "_mesh" << i_mesh;
2084  }
2085 
2086  std::ofstream outfile;
2087  char filename[100];
2088  sprintf(filename,"missing_coords_mesh%s.dat",modifier.str().c_str());
2089  outfile.open(filename);
2090  unsigned nel=mesh_pt[i_mesh]->nelement();
2091  for (unsigned e=0;e<nel;e++)
2092  {
2093  mesh_pt[i_mesh]->finite_element_pt(e)->
2094  FiniteElement::output(outfile);
2095  }
2096  outfile.close();
2097 
2098  sprintf(filename,"missing_coords_ext_mesh%s.dat",
2099  modifier.str().c_str());
2100  outfile.open(filename);
2101  nel=external_mesh_pt->nelement();
2102  for (unsigned e=0;e<nel;e++)
2103  {
2104  external_mesh_pt->finite_element_pt(e)->
2105  FiniteElement::output(outfile);
2106  }
2107  outfile.close();
2108 
2109  sprintf(filename,"missing_coords_bin%s.dat",modifier.str().c_str());
2110  outfile.open(filename);
2111  mesh_geom_obj_pt[i_mesh]->output_bins(outfile);
2112  outfile.close();
2113 
2114  sprintf(filename,"missing_coords%s.dat",modifier.str().c_str());
2115  outfile.open(filename);
2116  unsigned n=External_element_located.size();
2117  error_stream << "Number of unlocated elements " << n << std::endl;
2118  for (unsigned e=0;e<n;e++)
2119  {
2120  unsigned n_intpt=External_element_located[e].size();
2121  if (n_intpt==0)
2122  {
2123  error_stream
2124  << "Failure to locate in halo element! "
2125  << "Why are we searching there?"
2126  << std::endl;
2127  }
2128  for (unsigned ipt=0;ipt<n_intpt;ipt++)
2129  {
2130  if (External_element_located[e][ipt]==0)
2131  {
2132  error_stream << "Failure at element/intpt: "
2133  << e << " " << ipt << std::endl;
2134 
2135  // Cast
2137  dynamic_cast<ElementWithExternalElement*>(
2138  mesh_pt[i_mesh]->element_pt(e));
2139 
2140  unsigned n_dim_el=el_pt->dim();
2141  Vector<double> s(n_dim_el);
2142  for (unsigned i=0;i<n_dim_el;i++)
2143  {
2144  s[i]=el_pt->integral_pt()->knot(ipt,i);
2145  }
2146  unsigned n_dim=el_pt->node_pt(0)->ndim();
2147  Vector<double> x(n_dim);
2148  el_pt->interpolated_x(s,x);
2149  for (unsigned i=0;i<n_dim;i++)
2150  {
2151  error_stream << x[i] << " ";
2152  outfile<< x[i] << " ";
2153  }
2154  error_stream << std::endl;
2155  outfile << std::endl;
2156  }
2157  }
2158  }
2159  outfile.close();
2160  }
2161 
2162  error_stream
2163  << "Mesh and external mesh documented in missing_coords_mesh*.dat\n"
2164  << "and missing_coords_ext_mesh*.dat, respectively. Missing \n"
2165  << "coordinates in missing_coords*.dat\n";
2166  throw OomphLibError
2167  (error_stream.str(),
2168  OOMPH_CURRENT_FUNCTION,
2169  OOMPH_EXCEPTION_LOCATION);
2170  }
2171  // Failure is deeemed to be acceptable
2172  else
2173  {
2174  oomph_info
2175  << "NOTE: Haven't found " << n_zeta_not_found
2176  << " external elements in \n"
2177  << "Multi_domain_functions::aux_setup_multi_domain_interaction(...)\n"
2178  << "but this deemed to be acceptable because \n"
2179  << "Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_domain_interaction\n"
2180  << "is true.\n";
2181  }
2182  }
2183 
2184 
2185  // Doc timings if required
2186  if (Doc_timings)
2187  {
2188  t_locate=TimingHelpers::timer();
2189  oomph_info
2190  << "Total CPU for location and creation of all external elements: "
2191  << t_locate-t_start << std::endl;
2192  }
2193 
2194  // Delete the geometric object representing the mesh
2195  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
2196  {
2197  delete mesh_geom_obj_pt[i_mesh];
2198  }
2199 
2200  // Clean up all the (extern) Vectors associated with creating the
2201  // external storage information
2202  clean_up();
2203 
2204 #ifdef OOMPH_HAS_MPI
2205  // Output information about external storage if required
2206  if (Doc_stats)
2207  {
2208  // Report stats regarding location method
2209  bool comm_was_required=false;
2210  oomph_info << "-------------------------------------------" << std::endl;
2211  oomph_info << "- Cumulative percentage of locate success -" << std::endl;
2212  oomph_info << "--- Spiral -- Found local -- Found else ---" << std::endl;
2213  for (unsigned level=0; level<max_level_reached; level++)
2214  {
2215  oomph_info << "--- " << level << " -- "
2216  << percentage_coords_located_locally[level] << " -- "
2217  << percentage_coords_located_elsewhere[level] << " ---"
2218  << std::endl;
2219  // Has communication with other processors at this level actually
2220  // produced any results?
2221  if (percentage_coords_located_elsewhere[level]>
2222  percentage_coords_located_locally[level])
2223  {
2224  comm_was_required=true;
2225  }
2226  }
2227  oomph_info << "-------------------------------------------" << std::endl;
2228 
2229  // Initialise to indicate that none of the zetas required
2230  // on this processor were located through parallel ring search,
2231  // i.e. comm was not required and we could have done some
2232  // more local searching first
2233  oomph_info << std::endl;
2234  oomph_info <<"ASSESSMENT OF NEED FOR PARALLEL SEARCH: \n";
2235  oomph_info <<"=======================================\n";
2236  unsigned status=0;
2237  if (comm_was_required)
2238  {
2239  oomph_info
2240  <<"- Ring-based parallel search did successfully locate zetas on proc "
2241  << my_rank << std::endl;
2242  status=1;
2243  }
2244  else
2245  {
2246  if (max_level_reached>1)
2247  {
2248  oomph_info
2249  << "- Ring-based parallel search did NOT locate zetas on proc"
2250  << my_rank << std::endl;
2251  }
2252  else
2253  {
2254  oomph_info
2255  << "- No ring-based parallel search was performed on proc"
2256  << my_rank << std::endl;
2257  }
2258  }
2259 
2260  // Allreduce to check if anyone has benefitted from parallel ring
2261  // search
2262  unsigned overall_status=0;
2263  MPI_Allreduce(&status,&overall_status,1,
2264  MPI_UNSIGNED,MPI_MAX,
2265  problem_pt->communicator_pt()->mpi_comm());
2266 
2267  // Report of mpi was useful to anyone
2268  if (overall_status==1)
2269  {
2270  oomph_info << "- Ring-based, parallel search did succesfully\n";
2271  oomph_info << " locate zetas on at least one other proc, so it\n";
2272  oomph_info << " was worthwhile.\n";
2273  oomph_info << std::endl;
2274  }
2275  else
2276  {
2277  if (max_level_reached>1)
2278  {
2279  oomph_info << "- Ring-based, parallel search did NOT locate zetas\n";
2280  oomph_info << " on ANY other procs, i.e it was useless.\n";
2281  oomph_info << " --> We should really have done more local search\n";
2282  oomph_info << " by reducing number of bins, or doing more spirals\n";
2283  oomph_info << " in one go before initiating parallel search.\n";
2284  oomph_info << std::endl;
2285  }
2286  else
2287  {
2288  oomph_info << "- No ring-based, parallel search was performed\n";
2289  oomph_info << " or necessary. Perfect!\n";
2290  oomph_info << std::endl;
2291  }
2292  }
2293 
2294 
2295  // How many external elements does the external mesh have now?
2296  oomph_info << "------------------------------------------" << std::endl;
2297  oomph_info << "External mesh: I have " << external_mesh_pt->nelement()
2298  << " elements, and " << std::endl
2299  << external_mesh_pt->nexternal_halo_element()
2300  << " external halo elements, "
2301  << external_mesh_pt->nexternal_haloed_element()
2302  << " external haloed elements"
2303  << std::endl;
2304 
2305  // How many external nodes does each submesh have now?
2306  oomph_info << "------------------------------------------" << std::endl;
2307  oomph_info << "External mesh: I have " << external_mesh_pt->nnode()
2308  << " nodes, and " << std::endl
2309  << external_mesh_pt->nexternal_halo_node()
2310  << " external halo nodes, "
2311  << external_mesh_pt->nexternal_haloed_node()
2312  << " external haloed nodes"
2313  << std::endl;
2314  oomph_info << "------------------------------------------" << std::endl;
2315  }
2316 
2317  // Output further information about (external) halo(ed)
2318  // elements and nodes if required
2319  if (Doc_full_stats)
2320  {
2321  // How many elements does this submesh have for each of the processors?
2322  for (int iproc=0;iproc<n_proc;iproc++)
2323  {
2324  oomph_info << "----------------------------------------" << std::endl;
2325  oomph_info << "With process " << iproc << " there are "
2326  << external_mesh_pt->nroot_halo_element(iproc)
2327  << " root halo elements, and "
2328  << external_mesh_pt->nroot_haloed_element(iproc)
2329  << " root haloed elements" << std::endl
2330  << "and there are "
2331  << external_mesh_pt->nexternal_halo_element(iproc)
2332  << " external halo elements, and "
2333  << external_mesh_pt->nexternal_haloed_element(iproc)
2334  << " external haloed elements." << std::endl;
2335 
2336  oomph_info << "----------------------------------------" << std::endl;
2337  oomph_info << "With process " << iproc << " there are "
2338  << external_mesh_pt->nhalo_node(iproc)
2339  << " halo nodes, and "
2340  << external_mesh_pt->nhaloed_node(iproc)
2341  << " haloed nodes" << std::endl
2342  << "and there are "
2343  << external_mesh_pt->nexternal_halo_node(iproc)
2344  << " external halo nodes, and "
2345  << external_mesh_pt->nexternal_haloed_node(iproc)
2346  << " external haloed nodes." << std::endl;
2347  }
2348  oomph_info << "-----------------------------------------" << std::endl
2349  << std::endl;
2350  }
2351 
2352  #endif
2353 
2354  // Doc timings if required
2355  if (Doc_timings)
2356  {
2357  t_end=TimingHelpers::timer();
2358  oomph_info << "CPU for (one way) aux_setup_multi_domain_interaction: "
2359  << t_end-t_start << std::endl;
2360  }
2361 
2362  } // end of aux_setup_multi_domain_interaction
2363 
2364 #ifdef OOMPH_HAS_MPI
2365 
2366 
2367 //=====================================================================
2368 /// Creates external (halo) elements on the loop process based on the
2369 /// information received from each locate_zeta call on other processes
2370 //=====================================================================
2371  template<class EXT_ELEMENT>
2373  (int& iproc, Mesh* const &mesh_pt, Mesh* const &external_mesh_pt,
2374  Problem* problem_pt, const unsigned& interaction_index)
2375  {
2376  OomphCommunicator* comm_pt=problem_pt->communicator_pt();
2377  int my_rank=comm_pt->my_rank();
2378 
2379  // Reset counters for flat packed unsigneds (namespace data because
2380  // it's also accessed by helper functions)
2383 
2384  // Initialise counter for stepping through zetas
2385  unsigned zeta_counter=0;
2386 
2387  // Initialise counter for stepping through flat-packed located
2388  // coordinates
2389  unsigned counter_for_located_coord=0;
2390 
2391  // The creation all happens on the current processor
2392  // Loop over this processors elements
2393  unsigned n_element=mesh_pt->nelement();
2394  for (unsigned e=0;e<n_element;e++)
2395  {
2396  // Cast to ElementWithExternalElement to set external element (if located)
2398  dynamic_cast<ElementWithExternalElement*>(mesh_pt->element_pt(e));
2399 
2400  // We're not setting up external elements for halo elements
2401  if (!el_pt->is_halo())
2402  {
2403  // Loop over integration points
2404  unsigned n_intpt=el_pt->integral_pt()->nweight();
2405  for (unsigned ipt=0;ipt<n_intpt;ipt++)
2406  {
2407  // Has an external element been assigned to this integration point?
2408  if (External_element_located[e][ipt]==0)
2409  {
2410  // Was a (non-halo) element located for this integration point
2411  if (((Proc_id_plus_one_of_external_element[zeta_counter]-1)==
2412  my_rank) ||
2413  (Proc_id_plus_one_of_external_element[zeta_counter]==0))
2414  {
2415  // Either it was already found, or not found on the current proc.
2416  // In either case, we don't need to do anything for this
2417  // integration point
2418  }
2419  else
2420  {
2421  // Get the process number on which the element was located
2422  unsigned loc_p=
2423  Proc_id_plus_one_of_external_element[zeta_counter]-1;
2424 
2425  // Is it a new external halo element or not?
2426  // If so, then create it, populate it, and add it as a
2427  // source; if not, then find the right one which
2428  // has already been created and use it as the source
2429  // element.
2430 
2431  // FiniteElement stored at this integration point
2432  FiniteElement* f_el_pt=0;
2433 
2434  // Is it a new element?
2435  if (Located_element_status[zeta_counter]==New)
2436  {
2437  // Create a new element from the communicated values
2438  // and coords from the process that located zeta
2439  GeneralisedElement *new_el_pt= new EXT_ELEMENT;
2440 
2441  // Add external halo element to this mesh
2442  external_mesh_pt->
2443  add_external_halo_element_pt(loc_p, new_el_pt);
2444 
2445  // Cast to the FE pointer
2446  f_el_pt=dynamic_cast<FiniteElement*>(new_el_pt);
2447 
2448  // We need the number of interpolated values if Refineable
2449  int n_cont_inter_values=-1;
2450  if (dynamic_cast<RefineableElement*>(new_el_pt)!=0)
2451  {
2452  n_cont_inter_values=dynamic_cast<RefineableElement*>
2453  (new_el_pt)->ncont_interpolated_values();
2454  }
2455 
2456  // If we're using macro elements to update
2457 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2458  oomph_info
2459  << "Rec:" << Counter_for_flat_packed_unsigneds
2460  << " Using macro element node update "
2462  << std::endl;
2463 #endif
2465  ==1)
2466  {
2467  // Set the macro element
2468  MacroElementNodeUpdateMesh* macro_mesh_pt=
2469  dynamic_cast<MacroElementNodeUpdateMesh*>
2470  (external_mesh_pt);
2471 
2472 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2473  oomph_info
2474  << "Rec:" << Counter_for_flat_packed_unsigneds
2475  << " Number of macro element "
2477  << std::endl;
2478 #endif
2479  unsigned macro_el_num=
2481  f_el_pt->set_macro_elem_pt
2482  (macro_mesh_pt->macro_domain_pt()->
2483  macro_element_pt(macro_el_num));
2484 
2485 
2486  // We need to receive the lower left
2487  // and upper right coordinates of the macro element
2488  QElementBase* q_el_pt=
2489  dynamic_cast<QElementBase*>(new_el_pt);
2490  if (q_el_pt!=0)
2491  {
2492  unsigned el_dim=q_el_pt->dim();
2493  for (unsigned i_dim=0;i_dim<el_dim;i_dim++)
2494  {
2495  q_el_pt->s_macro_ll(i_dim)=
2497  q_el_pt->s_macro_ur(i_dim)=
2499  }
2500  }
2501  else // Throw an error, since this is only implemented for Q
2502  {
2503  std::ostringstream error_stream;
2504  error_stream << "Using MacroElement node update\n"
2505  << "in a case with non-QElements\n"
2506  << "has not yet been implemented.\n";
2507  throw OomphLibError
2508  (error_stream.str(),
2509  OOMPH_CURRENT_FUNCTION,
2510  OOMPH_EXCEPTION_LOCATION);
2511 
2512  }
2513  }
2514 
2515  // Now we add nodes to the new element
2516  unsigned n_node=f_el_pt->nnode();
2517  for (unsigned j=0;j<n_node;j++)
2518  {
2519  Node* new_nod_pt=0;
2520 
2521  // Call the add external halo node helper function
2522  add_external_halo_node_to_storage<EXT_ELEMENT>
2523  (new_nod_pt,external_mesh_pt,loc_p,j,f_el_pt,
2524  n_cont_inter_values,problem_pt);
2525  }
2526  }
2527  else // the element already exists as an external_halo
2528  {
2529 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2530  oomph_info
2531  << "Rec:" << Counter_for_flat_packed_unsigneds
2532  << " Index of existing external halo element "
2534  << std::endl;
2535 #endif
2536  // The index itself is in Flat_packed_unsigneds[...]
2537  unsigned external_halo_el_index=
2539 
2540  // Use this index to get the element
2541  f_el_pt=dynamic_cast<FiniteElement*>(external_mesh_pt->
2542  external_halo_element_pt
2543  (loc_p,external_halo_el_index));
2544 
2545  //If it's not a finite element die
2546  if(f_el_pt==0)
2547  {
2548  throw OomphLibError(
2549  "External halo element is not a FiniteElement\n",
2550  OOMPH_CURRENT_FUNCTION,
2551  OOMPH_EXCEPTION_LOCATION);
2552  }
2553  }
2554 
2555  // The source element storage was initialised but
2556  // not filled earlier, so do it now
2557  // The located coordinates are required
2558  // (which could be a different dimension to zeta, e.g. in FSI)
2559  unsigned el_dim=f_el_pt->dim();
2560  Vector<double> s_located(el_dim);
2561  for (unsigned i=0;i<el_dim;i++)
2562  {
2563  s_located[i]=
2564  Flat_packed_located_coordinates[counter_for_located_coord];
2565  counter_for_located_coord++;
2566  }
2567 
2568  // Set the element for this integration point
2569  el_pt->external_element_pt(interaction_index,ipt)=f_el_pt;
2570  el_pt->
2571  external_element_local_coord(interaction_index,ipt)=s_located;
2572 
2573  // Set the lookup array to true
2574  External_element_located[e][ipt]=1;
2575  }
2576 
2577  // Increment the integration point counter
2578  zeta_counter++;
2579  }
2580  } // end loop over integration points
2581  }
2582  } // end loop over local processor's elements
2583 
2584  }
2585 
2586 
2587 
2588 //=====================================================================
2589 // vector based version
2590 
2591 /// Creates external (halo) elements on the loop process based on the
2592 /// information received from each locate_zeta call on other processes
2593 //=====================================================================
2594  template<class EXT_ELEMENT>
2596  (int& iproc, const Vector<Mesh*>& mesh_pt, Mesh* const &external_mesh_pt,
2597  Problem* problem_pt, const unsigned& interaction_index)
2598  {
2599  OomphCommunicator* comm_pt=problem_pt->communicator_pt();
2600  int my_rank=comm_pt->my_rank();
2601 
2602  // Reset counters for flat packed unsigneds (namespace data because
2603  // it's also accessed by helper functions)
2606 
2607  // Initialise counter for stepping through zetas
2608  unsigned zeta_counter=0;
2609 
2610  // Initialise counter for stepping through flat-packed located
2611  // coordinates
2612  unsigned counter_for_located_coord=0;
2613 
2614  // Counter for elements in flag packed storage
2615  unsigned e_count=0;
2616 
2617 
2618  // Loop over all meshes
2619  unsigned n_mesh=mesh_pt.size();
2620  for (unsigned i_mesh=0;i_mesh<n_mesh;i_mesh++)
2621  {
2622 
2623  // The creation all happens on the current processor
2624  // Loop over this processors elements
2625  unsigned n_element=mesh_pt[i_mesh]->nelement();
2626  for (unsigned e=0;e<n_element;e++)
2627  {
2628  // Cast to ElementWithExternalElement to set external element
2629  // (if located)
2631  dynamic_cast<ElementWithExternalElement*>(mesh_pt[i_mesh]
2632  ->element_pt(e));
2633 
2634  // We're not setting up external elements for halo elements
2635  // (Note: this is consistent with padding introduced when
2636  // External_element_located was first assigned)
2637  if (!el_pt->is_halo())
2638  {
2639  // Loop over integration points
2640  unsigned n_intpt=el_pt->integral_pt()->nweight();
2641  for (unsigned ipt=0;ipt<n_intpt;ipt++)
2642  {
2643  // Has an external element been assigned to this integration point?
2644  if (External_element_located[e_count][ipt]==0)
2645  {
2646  // Was a (non-halo) element located for this integration point
2647  if (((Proc_id_plus_one_of_external_element[zeta_counter]-1)==
2648  my_rank) ||
2649  (Proc_id_plus_one_of_external_element[zeta_counter]==0))
2650  {
2651  // Either it was already found, or not found on the current proc.
2652  // In either case, we don't need to do anything for this
2653  // integration point
2654  }
2655  else
2656  {
2657  // Get the process number on which the element was located
2658  unsigned loc_p=
2659  Proc_id_plus_one_of_external_element[zeta_counter]-1;
2660 
2661  // Is it a new external halo element or not?
2662  // If so, then create it, populate it, and add it as a
2663  // source; if not, then find the right one which
2664  // has already been created and use it as the source
2665  // element.
2666 
2667  // FiniteElement stored at this integration point
2668  FiniteElement* f_el_pt=0;
2669 
2670  // Is it a new element?
2671  if (Located_element_status[zeta_counter]==New)
2672  {
2673  // Create a new element from the communicated values
2674  // and coords from the process that located zeta
2675  GeneralisedElement *new_el_pt= new EXT_ELEMENT;
2676 
2677  // Add external halo element to this mesh
2678  external_mesh_pt->
2679  add_external_halo_element_pt(loc_p, new_el_pt);
2680 
2681  // Cast to the FE pointer
2682  f_el_pt=dynamic_cast<FiniteElement*>(new_el_pt);
2683 
2684  // We need the number of interpolated values if Refineable
2685  int n_cont_inter_values=-1;
2686  if (dynamic_cast<RefineableElement*>(new_el_pt)!=0)
2687  {
2688  n_cont_inter_values=dynamic_cast<RefineableElement*>
2689  (new_el_pt)->ncont_interpolated_values();
2690  }
2691 
2692  // If we're using macro elements to update
2693 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2694  oomph_info
2695  << "Rec:" << Counter_for_flat_packed_unsigneds
2696  << " Using macro element node update "
2698  << std::endl;
2699 #endif
2701  ==1)
2702  {
2703  // Set the macro element
2704  MacroElementNodeUpdateMesh* macro_mesh_pt=
2705  dynamic_cast<MacroElementNodeUpdateMesh*>
2706  (external_mesh_pt);
2707 
2708 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2709  oomph_info
2710  << "Rec:" << Counter_for_flat_packed_unsigneds
2711  << " Number of macro element "
2713  << std::endl;
2714 #endif
2715  unsigned macro_el_num=
2717  f_el_pt->set_macro_elem_pt
2718  (macro_mesh_pt->macro_domain_pt()->
2719  macro_element_pt(macro_el_num));
2720 
2721 
2722  // We need to receive the lower left
2723  // and upper right coordinates of the macro element
2724  QElementBase* q_el_pt=
2725  dynamic_cast<QElementBase*>(new_el_pt);
2726  if (q_el_pt!=0)
2727  {
2728  unsigned el_dim=q_el_pt->dim();
2729  for (unsigned i_dim=0;i_dim<el_dim;i_dim++)
2730  {
2731  q_el_pt->s_macro_ll(i_dim)=
2733  q_el_pt->s_macro_ur(i_dim)=
2735  }
2736  }
2737  else // Throw an error, since this is only implemented for Q
2738  {
2739  std::ostringstream error_stream;
2740  error_stream << "Using MacroElement node update\n"
2741  << "in a case with non-QElements\n"
2742  << "has not yet been implemented.\n";
2743  throw OomphLibError
2744  (error_stream.str(),
2745  OOMPH_CURRENT_FUNCTION,
2746  OOMPH_EXCEPTION_LOCATION);
2747 
2748  }
2749  }
2750 
2751  // Now we add nodes to the new element
2752  unsigned n_node=f_el_pt->nnode();
2753  for (unsigned j=0;j<n_node;j++)
2754  {
2755  Node* new_nod_pt=0;
2756 
2757  // Call the add external halo node helper function
2758  add_external_halo_node_to_storage<EXT_ELEMENT>
2759  (new_nod_pt,external_mesh_pt,loc_p,j,f_el_pt,
2760  n_cont_inter_values,problem_pt);
2761  }
2762  }
2763  else // the element already exists as an external_halo
2764  {
2765 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2766  oomph_info
2767  << "Rec:" << Counter_for_flat_packed_unsigneds
2768  << " Index of existing external halo element "
2770  << std::endl;
2771 #endif
2772  // The index itself is in Flat_packed_unsigneds[...]
2773  unsigned external_halo_el_index=
2775 
2776  // Use this index to get the element
2777  f_el_pt=dynamic_cast<FiniteElement*>(external_mesh_pt->
2778  external_halo_element_pt
2779  (loc_p,
2780  external_halo_el_index));
2781 
2782  //If it's not a finite element die
2783  if(f_el_pt==0)
2784  {
2785  throw OomphLibError(
2786  "External halo element is not a FiniteElement\n",
2787  OOMPH_CURRENT_FUNCTION,
2788  OOMPH_EXCEPTION_LOCATION);
2789  }
2790  }
2791 
2792  // The source element storage was initialised but
2793  // not filled earlier, so do it now
2794  // The located coordinates are required
2795  // (which could be a different dimension to zeta, e.g. in FSI)
2796  unsigned el_dim=f_el_pt->dim();
2797  Vector<double> s_located(el_dim);
2798  for (unsigned i=0;i<el_dim;i++)
2799  {
2800  s_located[i]=
2801  Flat_packed_located_coordinates[counter_for_located_coord];
2802  counter_for_located_coord++;
2803  }
2804 
2805  // Set the element for this integration point
2806  el_pt->external_element_pt(interaction_index,ipt)=f_el_pt;
2807  el_pt->
2808  external_element_local_coord(interaction_index,ipt)=s_located;
2809 
2810  // Set the lookup array to true
2811  External_element_located[e_count][ipt]=1;
2812  }
2813 
2814  // Increment the integration point counter
2815  zeta_counter++;
2816  }
2817  } // end loop over integration points
2818  } // end of is halo
2819 
2820  // Bump flat-packed element counter
2821  e_count++;
2822 
2823  } // end of loop over elements
2824 
2825  // Bump up zeta counter to skip over padding entry at end of
2826  // mesh
2827  zeta_counter++;
2828 
2829  } // end loop over meshes
2830  }
2831 
2832 
2833 
2834 
2835 
2836 
2837 
2838 
2839 
2840 //============start of add_external_halo_node_to_storage===============
2841 /// Helper function to add external halo nodes, including any masters,
2842 /// based on information received from the haloed process
2843 //=========================================================================
2844  template<class EXT_ELEMENT>
2846  (Node* &new_nod_pt, Mesh* const &external_mesh_pt, unsigned& loc_p,
2847  unsigned& node_index, FiniteElement* const &new_el_pt,
2848  int& n_cont_inter_values,
2849  Problem* problem_pt)
2850  {
2851  // Add the external halo node if required
2852  add_external_halo_node_helper(new_nod_pt,external_mesh_pt,loc_p,
2853  node_index,new_el_pt,
2854  n_cont_inter_values,problem_pt);
2855 
2856  // Recursively add masters
2857  recursively_add_masters_of_external_halo_node_to_storage<EXT_ELEMENT>
2858  (new_nod_pt, external_mesh_pt, loc_p,
2859  node_index, new_el_pt,
2860  n_cont_inter_values,
2861  problem_pt);
2862 
2863 
2864  }
2865 
2866 
2867 //========================================================================
2868 /// Recursively add masters of external halo nodes (and their masters, etc)
2869 /// based on information received from the haloed process
2870 //=========================================================================
2871  template<class EXT_ELEMENT>
2874  (Node* &new_nod_pt, Mesh* const &external_mesh_pt, unsigned& loc_p,
2875  unsigned& node_index, FiniteElement* const &new_el_pt,
2876  int& n_cont_inter_values,
2877  Problem* problem_pt)
2878  {
2879 
2880  for (int i_cont=-1;i_cont<n_cont_inter_values;i_cont++)
2881  {
2882 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2883  oomph_info
2884  << "Rec:" << Counter_for_flat_packed_unsigneds
2885  << " Boolean to indicate that continuously interpolated variable i_cont "
2886  << i_cont << " is hanging "
2888  << std::endl;
2889 #endif
2891  {
2892 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2893  oomph_info
2894  << "Rec:" << Counter_for_flat_packed_unsigneds
2895  << " Number of master nodes "
2897  << std::endl;
2898 #endif
2899  unsigned n_master=Flat_packed_unsigneds
2901 
2902  // Setup new HangInfo
2903  HangInfo* hang_pt=new HangInfo(n_master);
2904  for (unsigned m=0;m<n_master;m++)
2905  {
2906  Node* master_nod_pt=0;
2907  // Get the master node (creating and adding it if required)
2908  add_external_halo_master_node_helper<EXT_ELEMENT>
2909  (master_nod_pt,new_nod_pt,external_mesh_pt,loc_p,
2910  n_cont_inter_values,problem_pt);
2911 
2912  // Get the weight and set the HangInfo
2913  double master_weight=Flat_packed_doubles
2915  hang_pt->set_master_node_pt(m,master_nod_pt,master_weight);
2916 
2917  // Recursively add masters of master
2918  recursively_add_masters_of_external_halo_node_to_storage<EXT_ELEMENT>
2919  (master_nod_pt, external_mesh_pt, loc_p,
2920  node_index, new_el_pt,
2921  n_cont_inter_values,
2922  problem_pt);
2923  }
2924  new_nod_pt->set_hanging_pt(hang_pt,i_cont);
2925  }
2926  } // end loop over continous interpolated values
2927 
2928  }
2929 
2930 //========================================================================
2931 /// Helper function to add external halo node that is a master
2932 //========================================================================
2933 template<class EXT_ELEMENT>
2935  (Node* &new_master_nod_pt, Node* &new_nod_pt, Mesh* const &external_mesh_pt,
2936  unsigned& loc_p, int& ncont_inter_values,Problem* problem_pt)
2937  {
2938  // Given the node and the external mesh, and received information
2939  // about them from process loc_p, construct them on the current process
2940 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2941  oomph_info
2942  << "Rec:" << Counter_for_flat_packed_unsigneds
2943  << " Boolean to trigger construction of new external halo master node "
2945  << std::endl;
2946 #endif
2948  {
2949  // Construct a new node based upon sent information
2950  construct_new_external_halo_master_node_helper<EXT_ELEMENT>
2951  (new_master_nod_pt,new_nod_pt,loc_p,external_mesh_pt,problem_pt);
2952  }
2953  else
2954  {
2955 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2956  oomph_info
2957  << "Rec:" << Counter_for_flat_packed_unsigneds
2958  << " index of existing external halo master node "
2960  << std::endl;
2961 #endif
2962  // Copy node from received location
2963  new_master_nod_pt=external_mesh_pt->external_halo_node_pt
2965  }
2966  }
2967 
2968 //======start of construct_new_external_halo_master_node_helper===========
2969 /// Helper function which constructs a new external halo master node
2970 /// with the required information sent from the haloed process
2971 //========================================================================
2972 template<class EXT_ELEMENT>
2974  (Node* &new_master_nod_pt, Node* &nod_pt, unsigned& loc_p,
2975  Mesh* const &external_mesh_pt, Problem* problem_pt)
2976  {
2977  // First three sent numbers are dimension, position type and nvalue
2978  // (to be used in Node constructors)
2979 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2980  oomph_info
2981  << "Rec:" << Counter_for_flat_packed_unsigneds
2982  << " ndim for external halo master node "
2984  << std::endl;
2985 #endif
2987 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2988  oomph_info
2989  << "Rec:" << Counter_for_flat_packed_unsigneds
2990  << " nposition type for external halo master node "
2992  << std::endl;
2993 #endif
2994  unsigned n_position_type=Flat_packed_unsigneds
2996 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
2997  oomph_info
2998  << "Rec:" << Counter_for_flat_packed_unsigneds
2999  << " nvalue for external halo master node "
3001  << std::endl;
3002 #endif
3003  unsigned n_value=Flat_packed_unsigneds
3005 
3006  // If it's a solid node also receive the lagrangian dimension and pos type
3007  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
3008  unsigned n_lag_dim;
3009  unsigned n_lag_type;
3010  if (solid_nod_pt!=0)
3011  {
3012 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3013  oomph_info
3014  << "Rec:" << Counter_for_flat_packed_unsigneds
3015  << " nlagrdim for external halo master solid node "
3017  << std::endl;
3018 #endif
3020 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3021  oomph_info
3022  << "Rec:" << Counter_for_flat_packed_unsigneds
3023  << " nlagrtype for external halo master solid node "
3025  << std::endl;
3026 #endif
3028  }
3029 
3030  // Null TimeStepper for now
3031  TimeStepper* time_stepper_pt=0;
3032  // Default number of previous values to 1
3033  unsigned n_prev=1;
3034 
3035  // The first entry in nodal_info indicates
3036  // the timestepper required for this halo node
3037 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3038  oomph_info
3039  << "Rec:" << Counter_for_flat_packed_unsigneds
3040  << " Boolean: need timestepper "
3042  << std::endl;
3043 #endif
3045  {
3046 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3047  oomph_info
3048  << "Rec:" << Counter_for_flat_packed_unsigneds
3049  << " Index minus one of timestepper "
3051  << std::endl;
3052 #endif
3053  // Index minus one!
3054  time_stepper_pt=problem_pt->time_stepper_pt
3056 
3057  // Check whether number of prev values is "sent" across
3058  n_prev=time_stepper_pt->ntstorage();
3059  }
3060 
3061  // Is the node for which the master is required Algebraic, Macro or Solid?
3062  AlgebraicNode* alg_nod_pt=dynamic_cast<AlgebraicNode*>(nod_pt);
3063  MacroElementNodeUpdateNode* macro_nod_pt=
3064  dynamic_cast<MacroElementNodeUpdateNode*>(nod_pt);
3065 
3066  // What type of node was the node for which we are constructing a master?
3067  if (alg_nod_pt!=0)
3068  {
3069  // The master node should also be algebraic
3070 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3071  oomph_info
3072  << "Rec:" << Counter_for_flat_packed_unsigneds
3073  << " Boolean for algebraic boundary node "
3075  << std::endl;
3076 #endif
3077  // If this master node's haloed copy is on a boundary then
3078  // it needs to be on the same boundary here
3080  {
3081  // Create a new BoundaryNode (not attached to an element)
3082  if (time_stepper_pt!=0)
3083  {
3084  new_master_nod_pt = new BoundaryNode<AlgebraicNode>
3085  (time_stepper_pt,n_dim,n_position_type,n_value);
3086  }
3087  else
3088  {
3089  new_master_nod_pt = new BoundaryNode<AlgebraicNode>
3090  (n_dim,n_position_type,n_value);
3091  }
3092 
3093  // How many boundaries does the algebraic master node live on?
3094 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3096  << " Number of boundaries the algebraic master node is on: "
3098  << std::endl;
3099 #endif
3101  for (unsigned i=0;i<nb;i++)
3102  {
3103  // Boundary number
3104 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3106  << " Algebraic master node is on boundary "
3108  << std::endl;
3109 #endif
3110  unsigned i_bnd=
3112  external_mesh_pt->add_boundary_node(i_bnd,new_master_nod_pt);
3113  }
3114 
3115 
3116  // Do we have additional values created by face elements?
3117 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3118  oomph_info
3119  << "Rec:" << Counter_for_flat_packed_unsigneds << " "
3120  << "Number of additional values created by face element "
3121  << "for master node "
3123  << std::endl;
3124 #endif
3125  unsigned n_entry=Flat_packed_unsigneds[
3127  if (n_entry>0)
3128  {
3129  // Create storage, if it doesn't already exist, for the map
3130  // that will contain the position of the first entry of
3131  // this face element's additional values,
3132  BoundaryNodeBase* bnew_master_nod_pt=
3133  dynamic_cast<BoundaryNodeBase*>(new_master_nod_pt);
3134 #ifdef PARANOID
3135  if (bnew_master_nod_pt==0)
3136  {
3137  throw OomphLibError(
3138  "Failed to cast new node to boundary node\n",
3139  OOMPH_CURRENT_FUNCTION,
3140  OOMPH_EXCEPTION_LOCATION);
3141  }
3142 #endif
3143  if(bnew_master_nod_pt->
3144  index_of_first_value_assigned_by_face_element_pt()==0)
3145  {
3146  bnew_master_nod_pt->
3147  index_of_first_value_assigned_by_face_element_pt()=
3148  new std::map<unsigned, unsigned>;
3149  }
3150 
3151  // Get pointer to the map of indices associated with
3152  // additional values created by face elements
3153  std::map<unsigned, unsigned>* map_pt=
3155 
3156  // Loop over number of entries in map
3157  for (unsigned i=0;i<n_entry;i++)
3158  {
3159  // Read out pairs...
3160 
3161 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3163  << " Key of map entry for master node"
3165  << std::endl;
3166 #endif
3167  unsigned first=Flat_packed_unsigneds[
3169 
3170 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3172  << " Value of map entry for master node"
3174  << std::endl;
3175 #endif
3176  unsigned second=Flat_packed_unsigneds[
3178 
3179  // ...and assign
3180  (*map_pt)[first]=second;
3181  }
3182  }
3183 
3184  }
3185  else
3186  {
3187  // Create node (not attached to any element)
3188  if (time_stepper_pt!=0)
3189  {
3190  new_master_nod_pt = new AlgebraicNode
3191  (time_stepper_pt,n_dim,n_position_type,n_value);
3192  }
3193  else
3194  {
3195  new_master_nod_pt = new AlgebraicNode
3196  (n_dim,n_position_type,n_value);
3197  }
3198  }
3199 
3200  // Add this as an external halo node BEFORE considering node update!
3201  external_mesh_pt->add_external_halo_node_pt(loc_p,new_master_nod_pt);
3202 
3203  // The external mesh is itself Algebraic...
3204  AlgebraicMesh* alg_mesh_pt=dynamic_cast<AlgebraicMesh*>
3205  (external_mesh_pt);
3206 
3207 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3208  oomph_info
3209  << "Rec:" << Counter_for_flat_packed_unsigneds
3210  << " algebraic node update id for master node "
3212  << std::endl;
3213 #endif
3214  /// The first entry of All_unsigned_values is the default node update id
3215  unsigned update_id=Flat_packed_unsigneds
3217 
3218  // Setup algebraic node update info for this new node
3219  Vector<double> ref_value;
3220 
3221 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3222  oomph_info
3223  << "Rec:" << Counter_for_flat_packed_unsigneds
3224  << " algebraic node number of ref values for master node "
3226  << std::endl;
3227 #endif
3228  // The size of this vector is in the next entry
3229  unsigned n_ref_val=Flat_packed_unsigneds
3231 
3232  // The reference values are in the subsequent entries of All_double_values
3233  ref_value.resize(n_ref_val);
3234  for (unsigned i_ref=0;i_ref<n_ref_val;i_ref++)
3235  {
3236  ref_value[i_ref]=Flat_packed_doubles
3238  }
3239 
3240  // Also require a Vector of geometric objects
3241  Vector<GeomObject*> geom_object_pt;
3242 
3243 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3244  oomph_info
3245  << "Rec:" << Counter_for_flat_packed_unsigneds
3246  << " algebraic node number of geom objects for master node "
3248  << std::endl;
3249 #endif
3250 
3251  // The size of this vector is in the next entry of All_unsigned_values
3252  unsigned n_geom_obj=Flat_packed_unsigneds
3254 
3255  // The remaining indices are in the rest of
3256  // All_alg_nodal_info
3257  geom_object_pt.resize(n_geom_obj);
3258  for (unsigned i_geom=0;i_geom<n_geom_obj;i_geom++)
3259  {
3260 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3261  oomph_info
3262  << "Rec:" << Counter_for_flat_packed_unsigneds
3263  << " algebraic node: " << i_geom << "-th out of "
3264  << n_geom_obj << "-th geom index "
3266  << std::endl;
3267 #endif
3268  unsigned geom_index=Flat_packed_unsigneds
3270 
3271  // This index indicates which (if any) of the AlgebraicMesh's
3272  // stored geometric objects should be used
3273  geom_object_pt[i_geom]=alg_mesh_pt->geom_object_list_pt(geom_index);
3274  }
3275 
3276  AlgebraicNode* alg_master_nod_pt=
3277  dynamic_cast<AlgebraicNode*>(new_master_nod_pt);
3278 
3279  /// ... so for the specified update_id, call
3280  /// add_node_update_info
3281  alg_master_nod_pt->add_node_update_info
3282  (update_id,alg_mesh_pt,geom_object_pt,ref_value);
3283 
3284  /// Now call update_node_update
3285  alg_mesh_pt->update_node_update(alg_master_nod_pt);
3286  }
3287  else if (macro_nod_pt!=0)
3288  {
3289  // The master node should also be a macro node
3290  // If this master node's haloed copy is on a boundary then
3291  // it needs to be on the same boundary here
3292 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3294  << " Boolean for master algebraic node is boundary node "
3296  << std::endl;
3297 #endif
3299  {
3300  // Create a new BoundaryNode (not attached to an element)
3301  if (time_stepper_pt!=0)
3302  {
3303  new_master_nod_pt = new BoundaryNode<MacroElementNodeUpdateNode>
3304  (time_stepper_pt,n_dim,n_position_type,n_value);
3305  }
3306  else
3307  {
3308  new_master_nod_pt = new BoundaryNode<MacroElementNodeUpdateNode>
3309  (n_dim,n_position_type,n_value);
3310  }
3311 
3312 
3313  // How many boundaries does the macro element master node live on?
3314 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3315  oomph_info
3316  << "Rec:" << Counter_for_flat_packed_unsigneds
3317  << " Number of boundaries the macro element master node is on: "
3319  << std::endl;
3320 #endif
3322  for (unsigned i=0;i<nb;i++)
3323  {
3324  // Boundary number
3325 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3327  << " Macro element master node is on boundary "
3329  << std::endl;
3330 #endif
3331  unsigned i_bnd=
3333  external_mesh_pt->add_boundary_node(i_bnd,new_master_nod_pt);
3334  }
3335 
3336  // Do we have additional values created by face elements?
3337 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3338  oomph_info
3339  << "Rec:" << Counter_for_flat_packed_unsigneds
3340  << " Number of additional values created by face element "
3341  << "for macro element master node "
3343  << std::endl;
3344 #endif
3345  unsigned n_entry=Flat_packed_unsigneds[
3347  if (n_entry>0)
3348  {
3349  // Create storage, if it doesn't already exist, for the map
3350  // that will contain the position of the first entry of
3351  // this face element's additional values,
3352  BoundaryNodeBase* bnew_master_nod_pt=
3353  dynamic_cast<BoundaryNodeBase*>(new_master_nod_pt);
3354 #ifdef PARANOID
3355  if (bnew_master_nod_pt==0)
3356  {
3357  throw OomphLibError(
3358  "Failed to cast new node to boundary node\n",
3359  OOMPH_CURRENT_FUNCTION,
3360  OOMPH_EXCEPTION_LOCATION);
3361  }
3362 #endif
3363  if(bnew_master_nod_pt->
3364  index_of_first_value_assigned_by_face_element_pt()==0)
3365  {
3366  bnew_master_nod_pt->
3367  index_of_first_value_assigned_by_face_element_pt()=
3368  new std::map<unsigned, unsigned>;
3369  }
3370 
3371  // Get pointer to the map of indices associated with
3372  // additional values created by face elements
3373  std::map<unsigned, unsigned>* map_pt=
3375 
3376  // Loop over number of entries in map
3377  for (unsigned i=0;i<n_entry;i++)
3378  {
3379  // Read out pairs...
3380 
3381 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3383  << " Key of map entry for macro element master node"
3385  << std::endl;
3386 #endif
3387  unsigned first=Flat_packed_unsigneds[
3389 
3390 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3392  << " Value of map entry for macro element master node"
3394  << std::endl;
3395 #endif
3396  unsigned second=Flat_packed_unsigneds[
3398 
3399  // ...and assign
3400  (*map_pt)[first]=second;
3401  }
3402  }
3403 
3404  }
3405  else
3406  {
3407  // Create node (not attached to any element)
3408  if (time_stepper_pt!=0)
3409  {
3410  new_master_nod_pt = new MacroElementNodeUpdateNode
3411  (time_stepper_pt,n_dim,n_position_type,n_value);
3412  }
3413  else
3414  {
3415  new_master_nod_pt = new MacroElementNodeUpdateNode
3416  (n_dim,n_position_type,n_value);
3417  }
3418  }
3419 
3420  // Add this as an external halo node
3421  external_mesh_pt->add_external_halo_node_pt(loc_p,new_master_nod_pt);
3422 
3423  // Create a new node update element for this master node if required
3424  FiniteElement *new_node_update_f_el_pt=0;
3425 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3427  << " Bool: need new external halo element "
3429  << std::endl;
3430 #endif
3432  {
3433  GeneralisedElement* new_node_update_el_pt = new EXT_ELEMENT;
3434 
3435  //Add external hal element to this mesh
3436  external_mesh_pt->add_external_halo_element_pt(
3437  loc_p,new_node_update_el_pt);
3438 
3439  //Cast to finite element
3440  new_node_update_f_el_pt =
3441  dynamic_cast<FiniteElement*>(new_node_update_el_pt);
3442 
3443  // Need number of interpolated values if Refineable
3444  int n_cont_inter_values;
3445  if (dynamic_cast<RefineableElement*>(new_node_update_f_el_pt)!=0)
3446  {
3447  n_cont_inter_values=dynamic_cast<RefineableElement*>
3448  (new_node_update_f_el_pt)->ncont_interpolated_values();
3449  }
3450  else
3451  {
3452  n_cont_inter_values=-1;
3453  }
3454 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3456  << " Bool: we have a macro element mesh "
3458  << std::endl;
3459 #endif
3460  // If we're using macro elements to update,
3462  {
3463  // Set the macro element
3464  MacroElementNodeUpdateMesh* macro_mesh_pt=
3465  dynamic_cast<MacroElementNodeUpdateMesh*>
3466  (external_mesh_pt);
3467 
3468 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3470  << " Number of macro element "
3472  << std::endl;
3473 #endif
3474  unsigned macro_el_num=
3476  new_node_update_f_el_pt->set_macro_elem_pt
3477  (macro_mesh_pt->macro_domain_pt()->macro_element_pt(macro_el_num));
3478 
3479  // we need to receive
3480  // the lower left and upper right coordinates of the macro
3481  QElementBase* q_el_pt=
3482  dynamic_cast<QElementBase*>(new_node_update_f_el_pt);
3483  if (q_el_pt!=0)
3484  {
3485  unsigned el_dim=q_el_pt->dim();
3486  for (unsigned i_dim=0;i_dim<el_dim;i_dim++)
3487  {
3488  q_el_pt->s_macro_ll(i_dim)=Flat_packed_doubles
3490  q_el_pt->s_macro_ur(i_dim)=Flat_packed_doubles
3492  }
3493  }
3494  else // Throw an error
3495  {
3496  std::ostringstream error_stream;
3497  error_stream << "You are using a MacroElement node update\n"
3498  << "in a case with non-QElements. This has not\n"
3499  << "yet been implemented.\n";
3500  throw OomphLibError
3501  (error_stream.str(),
3502  OOMPH_CURRENT_FUNCTION,
3503  OOMPH_EXCEPTION_LOCATION);
3504  }
3505  }
3506 
3507  unsigned n_node=new_node_update_f_el_pt->nnode();
3508  for (unsigned j=0;j<n_node;j++)
3509  {
3510  Node* new_nod_pt=0;
3511  add_external_halo_node_to_storage<EXT_ELEMENT>
3512  (new_nod_pt,external_mesh_pt,loc_p,j,new_node_update_f_el_pt,
3513  n_cont_inter_values,problem_pt);
3514  }
3515 
3516  }
3517  else // The node update element exists already
3518  {
3519 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3521  << " Number of already existing external halo element "
3523  << std::endl;
3524 #endif
3525  new_node_update_f_el_pt=dynamic_cast<FiniteElement*>(
3526  external_mesh_pt->external_halo_element_pt
3528  }
3529 
3530  // Remaining required information to create functioning
3531  // MacroElementNodeUpdateNode...
3532 
3533  // Get the required geom objects for the node update
3534  // from the mesh
3535  Vector<GeomObject*> geom_object_vector_pt;
3536  MacroElementNodeUpdateMesh* macro_mesh_pt=
3537  dynamic_cast<MacroElementNodeUpdateMesh*>
3538  (external_mesh_pt);
3539  geom_object_vector_pt=macro_mesh_pt->geom_object_vector_pt();
3540 
3541  // Cast to MacroElementNodeUpdateNode
3542  MacroElementNodeUpdateNode* macro_master_nod_pt=
3543  dynamic_cast<MacroElementNodeUpdateNode*>(new_master_nod_pt);
3544 
3545  // Set all required information - node update element,
3546  // local coordinate in this element, and then set node update info
3547  macro_master_nod_pt->node_update_element_pt()=
3548  new_node_update_f_el_pt;
3549 
3550  // Need to get the local node index of the macro_master_nod_pt
3551  unsigned local_node_index;
3552  unsigned n_node=new_node_update_f_el_pt->nnode();
3553  for (unsigned j=0;j<n_node;j++)
3554  {
3555  if (macro_master_nod_pt==new_node_update_f_el_pt->node_pt(j))
3556  {
3557  local_node_index=j;
3558  break;
3559  }
3560  }
3561 
3562  Vector<double> s_in_macro_node_update_element;
3563  new_node_update_f_el_pt->local_coordinate_of_node
3564  (local_node_index,s_in_macro_node_update_element);
3565 
3566  macro_master_nod_pt->set_node_update_info
3567  (new_node_update_f_el_pt,s_in_macro_node_update_element,
3568  geom_object_vector_pt);
3569  }
3570  else if (solid_nod_pt!=0)
3571  {
3572  // The master node should also be a SolidNode
3573  // If this node was on a boundary then it needs to
3574  // be on the same boundary here
3575 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3577  << " Bool master is a boundary (solid) node "
3579  << std::endl;
3580 #endif
3582  {
3583  // Construct a new boundary node
3584  if (time_stepper_pt!=0)
3585  {
3586  new_master_nod_pt=new BoundaryNode<SolidNode>
3587  (time_stepper_pt,n_lag_dim,n_lag_type,n_dim,n_position_type,n_value);
3588  }
3589  else
3590  {
3591  new_master_nod_pt=new BoundaryNode<SolidNode>
3592  (n_lag_dim,n_lag_type,n_dim,n_position_type,n_value);
3593  }
3594 
3595 
3596  // How many boundaries does the macro element master node live on?
3597 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3598  oomph_info
3599  << "Rec:" << Counter_for_flat_packed_unsigneds
3600  << " Number of boundaries the solid master node is on: "
3602  << std::endl;
3603 #endif
3605  for (unsigned i=0;i<nb;i++)
3606  {
3607  // Boundary number
3608 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3610  << " Solid master node is on boundary "
3612  << std::endl;
3613 #endif
3614  unsigned i_bnd=
3616  external_mesh_pt->add_boundary_node(i_bnd,new_master_nod_pt);
3617  }
3618 
3619  // Do we have additional values created by face elements?
3620 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3621  oomph_info
3622  << "Rec:" << Counter_for_flat_packed_unsigneds
3623  << " Number of additional values created by face element "
3624  << "for solid master node "
3626  << std::endl;
3627 #endif
3628  unsigned n_entry=Flat_packed_unsigneds[
3630  if (n_entry>0)
3631  {
3632  // Create storage, if it doesn't already exist, for the map
3633  // that will contain the position of the first entry of
3634  // this face element's additional values,
3635  BoundaryNodeBase* bnew_master_nod_pt=
3636  dynamic_cast<BoundaryNodeBase*>(new_master_nod_pt);
3637 #ifdef PARANOID
3638  if (bnew_master_nod_pt==0)
3639  {
3640  throw OomphLibError(
3641  "Failed to cast new node to boundary node\n",
3642  OOMPH_CURRENT_FUNCTION,
3643  OOMPH_EXCEPTION_LOCATION);
3644  }
3645 #endif
3646  if(bnew_master_nod_pt->
3647  index_of_first_value_assigned_by_face_element_pt()==0)
3648  {
3649  bnew_master_nod_pt->
3650  index_of_first_value_assigned_by_face_element_pt()=
3651  new std::map<unsigned, unsigned>;
3652  }
3653 
3654  // Get pointer to the map of indices associated with
3655  // additional values created by face elements
3656  std::map<unsigned, unsigned>* map_pt=
3658 
3659  // Loop over number of entries in map
3660  for (unsigned i=0;i<n_entry;i++)
3661  {
3662  // Read out pairs...
3663 
3664 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3666  << " Key of map entry for solid master node"
3668  << std::endl;
3669 #endif
3670  unsigned first=Flat_packed_unsigneds[
3672 
3673 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3675  << " Value of map entry for solid master node"
3677  << std::endl;
3678 #endif
3679  unsigned second=Flat_packed_unsigneds[
3681 
3682  // ...and assign
3683  (*map_pt)[first]=second;
3684  }
3685  }
3686 
3687  }
3688  else
3689  {
3690  // Construct an ordinary (non-boundary) node
3691  if (time_stepper_pt!=0)
3692  {
3693  new_master_nod_pt=new SolidNode
3694  (time_stepper_pt,n_lag_dim,n_lag_type,n_dim,n_position_type,n_value);
3695  }
3696  else
3697  {
3698  new_master_nod_pt=new SolidNode
3699  (n_lag_dim,n_lag_type,n_dim,n_position_type,n_value);
3700  }
3701  }
3702 
3703  // Add this as an external halo node
3704  external_mesh_pt->add_external_halo_node_pt(loc_p,new_master_nod_pt);
3705 
3706  // Copy across particular info required for SolidNode
3707  // NOTE: Are there any problems with additional values for SolidNodes?
3708  SolidNode* solid_master_nod_pt=dynamic_cast<SolidNode*>(new_master_nod_pt);
3709  unsigned n_solid_val=solid_master_nod_pt->variable_position_pt()->nvalue();
3710  for (unsigned i_val=0;i_val<n_solid_val;i_val++)
3711  {
3712  for (unsigned t=0;t<n_prev;t++)
3713  {
3714  solid_master_nod_pt->variable_position_pt()->
3715  set_value(t,i_val,
3717  }
3718  }
3719  }
3720  else // Just an ordinary node!
3721  {
3722 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3724  << " Bool node is on boundary "
3726  << std::endl;
3727 #endif
3728 
3729  // If this node was on a boundary then it needs to
3730  // be on the same boundary here
3732  {
3733  // Construct a new boundary node
3734  if (time_stepper_pt!=0)
3735  {
3736  new_master_nod_pt=new BoundaryNode<Node>
3737  (time_stepper_pt,n_dim,n_position_type,n_value);
3738  }
3739  else
3740  {
3741  new_master_nod_pt=new BoundaryNode<Node>
3742  (n_dim,n_position_type,n_value);
3743  }
3744 
3745  // How many boundaries does the master node live on?
3746 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3748  << " Number of boundaries the master node is on: "
3750  << std::endl;
3751 #endif
3753  for (unsigned i=0;i<nb;i++)
3754  {
3755  // Boundary number
3756 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3758  << " Master node is on boundary "
3760  << std::endl;
3761 #endif
3762  unsigned i_bnd=
3764  external_mesh_pt->add_boundary_node(i_bnd,new_master_nod_pt);
3765  }
3766 
3767 
3768  // Do we have additional values created by face elements?
3769 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3770  oomph_info
3771  << "Rec:" << Counter_for_flat_packed_unsigneds
3772  << " Number of additional values created by face element "
3773  << "for master node "
3775  << std::endl;
3776 #endif
3777  unsigned n_entry=Flat_packed_unsigneds[
3779  if (n_entry>0)
3780  {
3781  // Create storage, if it doesn't already exist, for the map
3782  // that will contain the position of the first entry of
3783  // this face element's additional values,
3784  BoundaryNodeBase* bnew_master_nod_pt=
3785  dynamic_cast<BoundaryNodeBase*>(new_master_nod_pt);
3786 #ifdef PARANOID
3787  if (bnew_master_nod_pt==0)
3788  {
3789  throw OomphLibError(
3790  "Failed to cast new node to boundary node\n",
3791  OOMPH_CURRENT_FUNCTION,
3792  OOMPH_EXCEPTION_LOCATION);
3793  }
3794 #endif
3795  if(bnew_master_nod_pt->
3796  index_of_first_value_assigned_by_face_element_pt()==0)
3797  {
3798  bnew_master_nod_pt->
3799  index_of_first_value_assigned_by_face_element_pt()=
3800  new std::map<unsigned, unsigned>;
3801  }
3802 
3803  // Get pointer to the map of indices associated with
3804  // additional values created by face elements
3805  std::map<unsigned, unsigned>* map_pt=
3807 
3808  // Loop over number of entries in map
3809  for (unsigned i=0;i<n_entry;i++)
3810  {
3811  // Read out pairs...
3812 
3813 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3815  << " Key of map entry for master node"
3817  << std::endl;
3818 #endif
3819  unsigned first=Flat_packed_unsigneds[
3821 
3822 #ifdef ANNOTATE_MULTI_DOMAIN_COMMUNICATION
3824  << " Value of map entry for master node"
3826  << std::endl;
3827 #endif
3828  unsigned second=Flat_packed_unsigneds[
3830 
3831  // ...and assign
3832  (*map_pt)[first]=second;
3833  }
3834  }
3835  }
3836  else
3837  {
3838  // Construct an ordinary (non-boundary) node
3839  if (time_stepper_pt!=0)
3840  {
3841  new_master_nod_pt=new Node
3842  (time_stepper_pt,n_dim,n_position_type,n_value);
3843  }
3844  else
3845  {
3846  new_master_nod_pt=new Node(n_dim,n_position_type,n_value);
3847  }
3848  }
3849 
3850  // Add this as an external halo node
3851  external_mesh_pt->add_external_halo_node_pt(loc_p,new_master_nod_pt);
3852  }
3853 
3854  // Remaining info received for all node types
3855  // Get copied history values
3856  // unsigned n_val=new_master_nod_pt->nvalue();
3857  for (unsigned i_val=0;i_val<n_value;i_val++)
3858  {
3859  for (unsigned t=0;t<n_prev;t++)
3860  {
3861  new_master_nod_pt->set_value(t,i_val,Flat_packed_doubles
3863  }
3864  }
3865 
3866  // Get copied history values for positions
3867  unsigned n_nod_dim=new_master_nod_pt->ndim();
3868  for (unsigned idim=0;idim<n_nod_dim;idim++)
3869  {
3870  for (unsigned t=0;t<n_prev;t++)
3871  {
3872  // Copy to coordinate
3873  new_master_nod_pt->x(t,idim)=Flat_packed_doubles
3875  }
3876  }
3877 
3878  }
3879 
3880 
3881 #endif
3882 
3883 
3884 }
3885 
3886 #endif
3887 
3888 
3889 
3890 
3891 
void initialise_external_element_storage()
Initialise storage for pointers to external elements and their local coordinates. This must be called...
A Generalised Element class.
Definition: elements.h:76
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
unsigned nexternal_halo_element()
Total number of external halo elements in this Mesh.
Definition: mesh.h:1913
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:246
void setup_multi_domain_interaction(Problem *problem_pt, Mesh *const &mesh_pt, Mesh *const &external_mesh_pt, const unsigned &interaction_index=0)
Function to set up the one-way multi-domain interaction for problems where the meshes pointed to by m...
void add_external_halo_node_pt(const unsigned &p, Node *&nod_pt)
Add external halo node whose non-halo (external) counterpart is held on processor p to the storage sc...
Definition: mesh.h:2055
void aux_setup_multi_domain_interaction(Problem *problem_pt, Mesh *const &mesh_pt, Mesh *const &external_mesh_pt, const unsigned &interaction_index, Mesh *const &external_face_mesh_pt=0)
Auxiliary helper function.
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:590
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
virtual void update_node_update(AlgebraicNode *&node_pt)=0
Update the node update info for given node, following mesh adaptation. Must be implemented for every ...
unsigned & current_min_spiral_level()
Access function to current min. spiral level.
void send_and_receive_located_info(int &iproc, Mesh *const &external_mesh_pt, Problem *problem_pt)
Helper function to send back any located information.
double Total_time_for_sorting_elements_in_bins
The timing for sorting the elements in the bins.
cstr elem_len * i
Definition: cfortran.h:607
unsigned Counter_for_flat_packed_unsigneds
Counter used when processing vector of flat-packed unsigneds – this is really "private" data...
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:2962
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
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 ...
GeneralisedElement *& external_halo_element_pt(const unsigned &p, const unsigned &e)
Access fct to the e-th external halo element in this Mesh whose non-halo counterpart is held on proce...
Definition: mesh.h:1941
The Problem class.
Definition: problem.h:152
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
void recursively_add_masters_of_external_halo_node_to_storage(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)
Recursively add masters of external halo nodes (and their masters, etc) based on information received...
unsigned nexternal_halo_node()
Total number of external halo nodes in this Mesh.
Definition: mesh.h:1997
unsigned Nz_bin
Number of bins in the third dimension in binning method in setup_multi_domain_interaction().
A general Finite Element class.
Definition: elements.h:1271
char t
Definition: cfortran.h:572
bool Accept_failed_locate_zeta_in_setup_multi_domain_interaction
Boolean to indicate that failure in setup multi domain functions is acceptable; defaults to false...
Definition: multi_domain.cc:62
bool Doc_full_stats
Boolean to indicate whether to output further info during setup_multi_domain_interaction() routines...
MacroElement * macro_element_pt(const unsigned &i)
Access to i-th macro element.
Definition: domain.h:100
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
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 ...
OomphInfo oomph_info
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
std::ofstream Doc_boundary_coordinate_file
Output file to document the boundary coordinate along the mesh boundary of the bulk mesh during call ...
Definition: multi_domain.cc:53
unsigned nroot_haloed_element()
Total number of root haloed elements in this Mesh.
Definition: mesh.h:1623
void set_master_node_pt(const unsigned &i, Node *const &master_node_pt, const double &weight)
Set the pointer to the i-th master node and its weight.
Definition: nodes.cc:1417
bool Compute_extreme_bin_coordinates
Bool to tell the MeshAsGeomObject whether to calculate the extreme coordinates of the bin structure...
e
Definition: cfortran.h:575
void set_boundary_number_in_bulk_mesh(const unsigned &b)
Set function for the boundary number in bulk mesh.
Definition: elements.h:4236
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1221
Vector< unsigned > Flat_packed_unsigneds
Vector of flat-packed unsigneds to be communicated with other processors – this is really "private" d...
GeomObject * geom_object_list_pt(const unsigned &i)
Access function to the ith GeomObject.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
void output_bins(std::ofstream &outfile)
Output bins.
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:66
void locate_zeta_for_missing_coordinates(int &iproc, Mesh *const &external_mesh_pt, Problem *problem_pt, MeshAsGeomObject *&mesh_geom_obj_pt)
Locate zeta for current set of missing coordinates.
void setup_bulk_elements_adjacent_to_face_mesh(Problem *problem_pt, Vector< unsigned > &boundary_in_bulk_mesh, Mesh *const &bulk_mesh_pt, Vector< Mesh * > &face_mesh_pt, const unsigned &interaction=0)
Identify the FaceElements (stored in the mesh pointed to by face_mesh_pt) that are adjacent to the bu...
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
bool Doc_timings
Boolean to indicate whether to doc timings or not.
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
Definition: nodes.h:1654
Base class for Qelements.
Definition: Qelements.h:106
virtual void set_macro_elem_pt(MacroElement *macro_elem_pt)
Set pointer to macro element – can be overloaded in derived elements to perform additional tasks...
Definition: elements.h:1825
double & s_macro_ur(const unsigned &i)
Access fct to the i-th coordinate of the element's "upper right" vertex in the associated MacroElemen...
Definition: Qelements.h:240
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
OomphCommunicator * communicator_pt() const
Definition: mesh.h:1317
A class that contains the information required by Nodes that are located on Mesh boundaries. A BoundaryNode of a particular type is obtained by combining a given Node with this class. By differentiating between Nodes and BoundaryNodes we avoid a lot of un-necessary storage in the bulk Nodes.
Definition: nodes.h:1854
void add_external_halo_node_helper(Node *&new_nod_pt, Mesh *const &mesh_pt, unsigned &loc_p, unsigned &node_index, FiniteElement *const &new_el_pt, int &n_cont_inter_values, unsigned &counter_for_recv_unsigneds, Vector< unsigned > &recv_unsigneds, unsigned &counter_for_recv_doubles, Vector< double > &recv_doubles)
Helper functiono to add external halo node that is not a master.
unsigned nhalo_node()
Total number of halo nodes in this Mesh.
Definition: mesh.h:1581
unsigned nlagrangian() const
Access function to # of Lagrangian coordinates.
Definition: geom_objects.h:182
void create_external_halo_elements(int &iproc, Mesh *const &mesh_pt, Mesh *const &external_mesh_pt, Problem *problem_pt, const unsigned &interaction_index)
Helper function to create external (halo) elements on the loop process based on the info received in ...
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
Vector< double > Flat_packed_located_coordinates
Vector of flat-packed local coordinates for zeta tuples that have been located.
void construct_new_external_halo_master_node_helper(Node *&new_master_nod_pt, Node *&nod_pt, unsigned &loc_p, Mesh *const &external_mesh_pt, Problem *problem_pt)
Helper function which constructs a new external halo master node with the information sent from the h...
A template Class for BoundaryNodes; that is Nodes that MAY live on the boundary of a Mesh...
Definition: nodes.h:67
Node *& external_halo_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th external halo node in this Mesh whose non-halo external counterpart is held on...
Definition: mesh.h:2064
unsigned & current_max_spiral_level()
Access function to current max. spiral level.
unsigned nexternal_haloed_element()
Total number of external haloed elements in this Mesh.
Definition: mesh.h:1957
static char t char * s
Definition: cfortran.h:572
unsigned nexternal_haloed_node()
Total number of external haloed nodes in this Mesh.
Definition: mesh.h:2099
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
Domain *& macro_domain_pt()
Broken assignment operator.
void clean_up()
Helper function that clears all the intermediate information used during the external storage creatio...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
TimeStepper *& time_stepper_pt()
Access function for the pointer to the first (presumably only) timestepper.
Definition: problem.h:1460
void send_and_receive_missing_zetas(Problem *problem_pt)
Helper function to send any "missing" zeta coordinates to the next process and receive any coordinate...
FiniteElement *& node_update_element_pt()
Pointer to finite element that performs the update by referring to its macro-element representation (...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double timer()
returns the time in seconds after some point in past
unsigned nhaloed_node()
Total number of haloed nodes in this Mesh.
Definition: mesh.h:1684
Vector< double > Flat_packed_doubles
Vector of flat-packed doubles to be communicated with other processors.
bool Sort_bin_entries
Bool to decide if to sort entries in bin during locate_zeta operation (default: false) ...
Class that contains data for hanging nodes.
Definition: nodes.h:684
bool Use_bulk_element_as_external
Boolean to indicate when to use the bulk element as the external element. Defaults to false...
unsigned Nx_bin
Number of bins in the first dimension in binning method in setup_multi_domain_interaction().
void add_external_halo_master_node_helper(Node *&new_master_nod_pt, Node *&new_nod_pt, Mesh *const &external_mesh_pt, unsigned &loc_p, int &n_cont_inter_values, Problem *problem_pt)
Helper function to add external halo node that is a master.
void set_hanging_pt(HangInfo *const &hang_pt, const int &i)
Set the hanging data for the i-th value. (hang_pt=0 to make non-hanging)
Definition: nodes.cc:1997
double & s_macro_ll(const unsigned &i)
Access fct to the i-th coordinate of the element's "lower left" vertex in the associated MacroElement...
Definition: Qelements.h:223
A Class for nodes that deform elastically (i.e. position is an unknown in the problem). The idea is that the Eulerian positions are stored in a Data object and the Lagrangian coordinates are stored in addition. The pointer that addresses the Eulerian positions is set to the pointer to Value in the Data object. Hence, SolidNode uses knowledge of the internal structure of Data and must be a friend of the Data class. In order to allow a mesh to deform via an elastic-style equation in deforming-domain problems, the positions are stored separately from the values, so that elastic problems may be combined with any other type of problem.
Definition: nodes.h:1568
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size (broken virtual) ...
Definition: elements.h:1795
unsigned nroot_halo_element()
Total number of root halo elements in this Mesh.
Definition: mesh.h:1535
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4277
Vector< int > Proc_id_plus_one_of_external_element
Proc_id_plus_one_of_external_element[i] contains the processor id (plus one) of the processor on whic...
Definition: multi_domain.cc:97
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
Vector< Vector< unsigned > > External_element_located
Lookup scheme for whether a local element's integration point has had an external element assigned to...
Definition: multi_domain.cc:74
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
Definition: elements.cc:4480
void set_node_update_info(FiniteElement *node_update_element_pt, const Vector< double > &s_in_node_update_element, const Vector< GeomObject * > &geom_object_pt)
Set node update information for node: Pass the pointer to the element that performs the update operat...
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2949
bool Doc_stats
Boolean to indicate whether to output basic info during setup_multi_domain_interaction() routines...
unsigned Nsample_points
(Measure of) the number of sampling points within the elements when populating the bin ...
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
A class to do comparison of the elements by lexicographic ordering, based on the boundary coordinates...
bool problem_has_been_distributed()
Access to Problem_has_been_distributed flag.
Definition: problem.h:2056
std::map< unsigned, unsigned > *& index_of_first_value_assigned_by_face_element_pt()
Return pointer to the map giving the index of the first face element value.
Definition: nodes.h:1909
void add_external_halo_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add external halo element whose non-halo counterpart is held on processor p to this Mesh...
Definition: mesh.h:1949
Vector< double > Flat_packed_zetas_not_found_locally
Vector of flat-packed zeta coordinates for which the external element could not be found during curre...
Definition: multi_domain.cc:80
unsigned N_spiral_chunk
Number of spirals to be searched in one go.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
void delete_all_external_storage()
Wipe the storage for all externally-based elements.
Definition: mesh.cc:8836
unsigned Ny_bin
Number of bins in the second dimension in binning method in setup_multi_domain_interaction().
void add_external_halo_node_to_storage(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 function to add external halo nodes, including any masters, based on information received from...
bool is_halo() const
Is this element a halo?
Definition: elements.h:1141
Vector< unsigned > Located_element_status
Vector to indicate (to another processor) whether a located element (that will have to represented as...
Vector< GeomObject * > geom_object_vector_pt()
Access function to the vector of GeomObject.
A general mesh class.
Definition: mesh.h:74
void locate_zeta_for_local_coordinates(Mesh *const &mesh_pt, Mesh *const &external_mesh_pt, MeshAsGeomObject *&mesh_geom_obj_pt, const unsigned &interaction_index)
locate zeta for current set of "local" coordinates
unsigned Counter_for_flat_packed_doubles
Counter used when processing vector of flat-packed doubles – this is really "private" data...
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:57
void setup_multi_domain_interactions(Problem *problem_pt, Mesh *const &first_mesh_pt, Mesh *const &second_mesh_pt, const unsigned &first_interaction=0, const unsigned &second_interaction=0)
Set up the two-way multi-domain interactions for the problem pointed to by problem_pt. Use this for cases where first_mesh_pt and second_mesh_pt occupy the same physical space and are populated by ELEMENT_0 and ELEMENT_1 respectively, and are combined to solve a single problem. The elements in two meshes interact both ways the elements in each mesh act as "external elements" for the elements in the "other" mesh. The interaction indices allow the specification of which interaction we're setting up in the two meshes. They default to zero, which is appropriate if there's only a single interaction.