Enumerations
enum	{ New , Exists , Not_found }
	Enumerators for element status in location procedure. More...

Functions
void	locate_zeta_for_local_coordinates (const Vector< Mesh * > &mesh_pt, Mesh const &external_mesh_pt, Vector< MeshAsGeomObject > &mesh_geom_obj_pt, const unsigned &interaction_index)

void	get_dim_helper (Problem problem_pt, Mesh const &mesh_pt, Mesh *const &external_mesh_pt)

void	clean_up ()

bool	first_closer_than_second (const std::pair< FiniteElement , Vector< double >> &p1, const std::pair< FiniteElement , Vector< double >> &p2)

template<class BULK_ELEMENT , unsigned DIM>
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)
	/ Templated helper functions for multi-domain methods using locate_zeta More...

template<class BULK_ELEMENT , unsigned DIM>
void	setup_bulk_elements_adjacent_to_face_mesh (Problem problem_pt, const unsigned &boundary_in_bulk_mesh, Mesh const &bulk_mesh_pt, Mesh *const &face_mesh_pt, const unsigned &interaction=0)

template<class ELEMENT_0 , class ELEMENT_1 >
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)

template<class EXT_ELEMENT >
void	setup_multi_domain_interaction (Problem problem_pt, Mesh const &mesh_pt, Mesh *const &external_mesh_pt, const unsigned &interaction_index=0)

template<class EXT_ELEMENT , class FACE_ELEMENT_GEOM_OBJECT >
void	setup_multi_domain_interaction (Problem problem_pt, Mesh const &mesh_pt, Mesh const &external_mesh_pt, Mesh const &external_face_mesh_pt, const unsigned &interaction_index=0)

template<class EXT_ELEMENT , class FACE_ELEMENT_GEOM_OBJECT >
void	setup_multi_domain_interaction (Problem problem_pt, const Vector< Mesh > &mesh_pt, Mesh const &external_mesh_pt, const Vector< Mesh > &external_face_mesh_pt, const unsigned &interaction_index=0)

template<class EXT_ELEMENT , class GEOM_OBJECT >
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. More...

template<class EXT_ELEMENT , class GEOM_OBJECT >
void	aux_setup_multi_domain_interaction (Problem problem_pt, const Vector< Mesh > &mesh_pt, Mesh const &external_mesh_pt, const unsigned &interaction_index, const Vector< Mesh > &external_face_mesh_pt)
	Auxiliary helper function. More...

Variables
std::ofstream	Doc_boundary_coordinate_file

bool	Accept_failed_locate_zeta_in_setup_multi_domain_interaction = false

unsigned	Dim

Vector< Vector< unsigned > >	External_element_located

Vector< double >	Flat_packed_zetas_not_found_locally

Vector< double >	Received_flat_packed_zetas_to_be_found

Vector< int >	Proc_id_plus_one_of_external_element

Vector< unsigned >	Located_element_status

Vector< double >	Flat_packed_located_coordinates

Vector< double >	Flat_packed_doubles

unsigned	Counter_for_flat_packed_doubles

Vector< unsigned >	Flat_packed_unsigneds

unsigned	Counter_for_flat_packed_unsigneds

bool	Use_bulk_element_as_external = false

bool	Allow_use_of_halo_elements_as_external_elements = true

bool	Allow_use_of_halo_elements_as_external_elements_for_projection = true

bool	Doc_timings = false
	Boolean to indicate whether to doc timings or not. More...

bool	Doc_stats = false

bool	Doc_full_stats = false

Vector< double >	Zeta_coords_for_further_away_comparison

Enumeration Type Documentation

◆ anonymous enum

anonymous enum

Enumerators for element status in location procedure.

Enumerator
New
Exists
Not_found

     {
       New,
       Exists,
       Not_found
     };

Function Documentation

◆ aux_setup_multi_domain_interaction() [1/2]

template<class EXT_ELEMENT , class GEOM_OBJECT >

void oomph::Multi_domain_functions::aux_setup_multi_domain_interaction	(	Problem *	problem_pt,
		const Vector< Mesh * > &	mesh_pt,
		Mesh *const &	external_mesh_pt,
		const unsigned &	interaction_index,
		const Vector< Mesh * > &	external_face_mesh_pt
	)

Auxiliary helper function.

This routine calls the locate_zeta routine (simultaneously on each processor for each individual processor's element set if necessary) and sets up the external (halo) element and node storage as necessary. The locate_zeta procedure here works for all multi-domain problems where either two meshes occupy the same physical space but have differing element types (e.g. a Boussinesq convection problem where AdvectionDiffusion elements interact with Navier-Stokes type elements) or two meshes interact along some boundary of the external mesh, represented by a "face mesh", such as an FSI problem.

Vector-based version operates simultaneously on the meshes contained in the vectors.

If it's is now zero then break out of the spirals loop

   {
     // How many meshes do we have?
     unsigned n_mesh = mesh_pt.size();
  
 #ifdef PARANOID
     if (external_face_mesh_pt.size() != n_mesh)
     {
       std::ostringstream error_stream;
       error_stream << "Sizes of external_face_mesh_pt [ "
                    << external_face_mesh_pt.size() << " ] and "
                    << "mesh_pt [ " << n_mesh << " ] don't match.\n";
       throw OomphLibError(
         error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Bail out?
     if (n_mesh == 0) return;
  
       // Multi-domain setup will not work for elements with
       // nonuniformly spaced nodes
       // Must check type of elements in the mesh and in the external mesh
       //(assume element type is the same for all elements in each mesh)
  
 #ifdef PARANOID
  
     // Pointer to first element in external mesh
     GeneralisedElement* ext_el_pt_0 = 0;
     if (external_mesh_pt->nelement() != 0)
     {
       ext_el_pt_0 = external_mesh_pt->element_pt(0);
     }
  
     // Loop over all meshes
     for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
     {
       // Get pointer to first element in each mesh
       GeneralisedElement* el_pt_0 = 0;
       if (mesh_pt[i_mesh]->nelement() != 0)
       {
         el_pt_0 = mesh_pt[i_mesh]->element_pt(0);
       }
  
       // Check they are not spectral elements
       if (dynamic_cast<SpectralElement*>(el_pt_0) != 0 ||
           dynamic_cast<SpectralElement*>(ext_el_pt_0) != 0)
       {
         throw OomphLibError(
           "Multi-domain setup does not work with spectral elements.",
           OOMPH_CURRENT_FUNCTION,
           OOMPH_EXCEPTION_LOCATION);
       }
  
       // Check they are not hp-refineable elements
       if (dynamic_cast<PRefineableElement*>(el_pt_0) != 0 ||
           dynamic_cast<PRefineableElement*>(ext_el_pt_0) != 0)
       {
         throw OomphLibError(
           "Multi-domain setup does not work with hp-refineable elements.",
           OOMPH_CURRENT_FUNCTION,
           OOMPH_EXCEPTION_LOCATION);
       }
     } // end over initial loop over meshes
  
 #endif
  
  
 #ifdef OOMPH_HAS_MPI
     // Storage for number of processors and my rank
     int n_proc = problem_pt->communicator_pt()->nproc();
     int my_rank = problem_pt->communicator_pt()->my_rank();
 #endif
  
     // Timing
     double t_start = 0.0;
     double t_end = 0.0;
     double t_local = 0.0;
     double t_set = 0.0;
     double t_locate = 0.0;
     double t_spiral_start = 0.0;
 #ifdef OOMPH_HAS_MPI
     double t_loop_start = 0.0;
     double t_sendrecv = 0.0;
     double t_missing = 0.0;
     double t_send_info = 0.0;
     double t_create_halo = 0.0;
 #endif
  
     if (Doc_timings)
     {
       t_start = TimingHelpers::timer();
     }
  
     // Initialize number of zeta coordinates not found yet
     unsigned n_zeta_not_found = 0;
  
     // Geometric objects used to represent the external (face) meshes
     Vector<MeshAsGeomObject*> mesh_geom_obj_pt(n_mesh, 0);
  
 #ifdef PARANOID
  
     // Initialise lagrangian dimension of element (test only)
     unsigned el_dim_lag = 0;
  
 #endif
  
     // Create mesh as geom objects for all meshes
     for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
     {
       // Are bulk elements used as external elements?
       if (!Use_bulk_element_as_external)
       {
         // Fix this when required
         if (n_mesh != 1)
         {
           std::ostringstream error_stream;
           error_stream
             << "Sorry I currently can't deal with non-bulk external elements\n"
             << "in multi-domain setup for multiple meshes.\n"
             << "The functionality should be easy to implement now that you\n"
             << "have a test case. If you're not willinig to do this, call\n"
             << "the multi-domain setup mesh-by-mesh instead (though this can\n"
             << "be costly in parallel because of global comms. \n";
           throw OomphLibError(error_stream.str(),
                               OOMPH_CURRENT_FUNCTION,
                               OOMPH_EXCEPTION_LOCATION);
         }
  
         // Set the geometric object from the external mesh
         mesh_geom_obj_pt[0] = new MeshAsGeomObject(external_mesh_pt);
       }
       else
       {
         // Set the geometric object from the external face mesh argument
         mesh_geom_obj_pt[i_mesh] =
           new MeshAsGeomObject(external_face_mesh_pt[i_mesh]);
       }
  
 #ifdef PARANOID
       unsigned old_el_dim_lag = el_dim_lag;
  
       // Set lagrangian dimension of element
       el_dim_lag = mesh_geom_obj_pt[i_mesh]->nlagrangian();
  
       // Check consistency
       if (i_mesh > 0)
       {
         if (el_dim_lag != old_el_dim_lag)
         {
           std::ostringstream error_stream;
           error_stream << "Lagrangian dimensions of elements don't match \n "
                        << "between meshes: " << el_dim_lag << " "
                        << old_el_dim_lag << "\n";
           throw OomphLibError(error_stream.str(),
                               OOMPH_CURRENT_FUNCTION,
                               OOMPH_EXCEPTION_LOCATION);
         }
       }
 #endif
  
  
     } // end of loop over meshes
  
     double t_setup_lookups = 0.0;
     if (Doc_timings)
     {
       t_set = TimingHelpers::timer();
       oomph_info << "CPU for creation of MeshAsGeomObjects and bin structure: "
                  << t_set - t_start << std::endl;
       t_setup_lookups = TimingHelpers::timer();
     }
  
     // Total number of integration points
     unsigned tot_int = 0;
  
     // Counter for total number of elements (in flat packed order)
     unsigned e_count = 0;
  
     // Loop over all meshes to get total number of elements
     for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
     {
       e_count += mesh_pt[i_mesh]->nelement();
     }
     External_element_located.resize(e_count);
  
     // Reset counter for elements in flat packed storage
     e_count = 0;
  
     // Loop over all meshes
     for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
     {
       // Loop over (this processor's) elements and set lookup array
       unsigned n_element = mesh_pt[i_mesh]->nelement();
       for (unsigned e = 0; e < n_element; e++)
       {
         // Zero-sized vector means its a halo
         External_element_located[e_count].resize(0);
         ElementWithExternalElement* el_pt =
           dynamic_cast<ElementWithExternalElement*>(
             mesh_pt[i_mesh]->element_pt(e));
  
 #ifdef OOMPH_HAS_MPI
         // We're not setting up external elements for halo elements
         if (!el_pt->is_halo())
 #endif
         {
           // We need to allocate storage for the external elements
           // within the element. Memory will actually only be
           // allocated the first time this function is called for
           // each element, or if the number of interactions or integration
           // points within the element has changed.
           el_pt->initialise_external_element_storage();
  
           // Clear any previous allocation
           unsigned n_intpt = el_pt->integral_pt()->nweight();
           for (unsigned ipt = 0; ipt < n_intpt; ipt++)
           {
             el_pt->external_element_pt(interaction_index, ipt) = 0;
           }
  
           External_element_located[e_count].resize(n_intpt);
           for (unsigned ipt = 0; ipt < n_intpt; ipt++)
           {
             External_element_located[e_count][ipt] = 0;
             tot_int++;
           }
         }
         // next element
         e_count++;
       }
     } // end of loop over meshes
  
     if (Doc_timings)
     {
       double t = TimingHelpers::timer();
       oomph_info
         << "CPU for setup of lookup schemes for located elements/coords: "
         << t - t_setup_lookups << std::endl;
     }
  
     // Initialise maximum spiral level within the cartesian bin structure
     // Used to terminate spiraling for non-refineable bin
     unsigned n_max_level = 0;
  
 #ifdef OOMPH_HAS_MPI
     unsigned max_level_reached = 1;
 #endif
  
     // Max. number of sample points -- used to decide on termination of
     // "spiraling"
     unsigned max_n_sample_points_of_sample_point_containers = 0;
  
     // Loop over all meshes
     for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
     {
       if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version() ==
           UseRefineableBinArray)
       {
         RefineableBinArray* bin_array_pt = dynamic_cast<RefineableBinArray*>(
           mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
  
         bin_array_pt
           ->last_sample_point_to_actually_lookup_during_locate_zeta() =
           bin_array_pt
             ->initial_last_sample_point_to_actually_lookup_during_locate_zeta();
         bin_array_pt
           ->first_sample_point_to_actually_lookup_during_locate_zeta() = 0;
  
         unsigned nsp =
           bin_array_pt->total_number_of_sample_points_computed_recursively();
         if (nsp > max_n_sample_points_of_sample_point_containers)
         {
           max_n_sample_points_of_sample_point_containers = nsp;
         }
  
  
 #ifdef OOMPH_HAS_MPI
         // If the mesh has been distributed we want the max. number
         // of sample points across all processors
         if (problem_pt->communicator_pt()->nproc() > 1)
         {
           unsigned local_max_n_sample_points_of_sample_point_containers =
             max_n_sample_points_of_sample_point_containers;
  
           // Get  maximum over all processors
           MPI_Allreduce(&local_max_n_sample_points_of_sample_point_containers,
                         &max_n_sample_points_of_sample_point_containers,
                         1,
                         MPI_UNSIGNED,
                         MPI_MAX,
                         problem_pt->communicator_pt()->mpi_comm());
         }
 #endif
       }
       else if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version() ==
                UseNonRefineableBinArray)
       {
         NonRefineableBinArray* bin_array_pt =
           dynamic_cast<NonRefineableBinArray*>(
             mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
  
         // Initialise spiral levels
         bin_array_pt->current_min_spiral_level() = 0;
         bin_array_pt->current_max_spiral_level() =
           bin_array_pt->n_spiral_chunk() - 1;
  
         // Find maximum spiral level within the cartesian bin structure
         n_max_level = bin_array_pt->max_bin_dimension();
  
         // Limit it
         if (bin_array_pt->current_max_spiral_level() > n_max_level)
         {
           bin_array_pt->current_max_spiral_level() = n_max_level - 1;
         }
       }
 #ifdef OOMPH_HAS_CGAL
       // CGAL
       else if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version() ==
                UseCGALSamplePointContainer)
       {
         CGALSamplePointContainer* bin_array_pt =
           dynamic_cast<CGALSamplePointContainer*>(
             mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
         bin_array_pt
           ->last_sample_point_to_actually_lookup_during_locate_zeta() =
           bin_array_pt
             ->initial_last_sample_point_to_actually_lookup_during_locate_zeta();
         bin_array_pt
           ->first_sample_point_to_actually_lookup_during_locate_zeta() = 0;
  
         unsigned nsp =
           bin_array_pt->total_number_of_sample_points_computed_recursively();
         if (nsp > max_n_sample_points_of_sample_point_containers)
         {
           max_n_sample_points_of_sample_point_containers = nsp;
         }
  
  
 #ifdef OOMPH_HAS_MPI
         // If the mesh has been distributed we want the max. number
         // of sample points across all processors
         if (problem_pt->communicator_pt()->nproc() > 1)
         {
           unsigned local_max_n_sample_points_of_sample_point_containers =
             max_n_sample_points_of_sample_point_containers;
  
           // Get  maximum over all processors
           MPI_Allreduce(&local_max_n_sample_points_of_sample_point_containers,
                         &max_n_sample_points_of_sample_point_containers,
                         1,
                         MPI_UNSIGNED,
                         MPI_MAX,
                         problem_pt->communicator_pt()->mpi_comm());
         }
 #endif
       }
 #endif // cgal
     }
  
  
     // Storage for info about coordinate location
     Vector<double> percentage_coords_located_locally;
     Vector<double> percentage_coords_located_elsewhere;
  
     // Loop over "spirals/levels" away from the current position
     // Note: All meshes go through their spirals simultaneously;
     // read out spiral level from first one
     unsigned i_level = 0;
     bool has_not_reached_max_level_of_search = true;
     while (has_not_reached_max_level_of_search)
     {
       // Record time at start of spiral loop
       if (Doc_timings)
       {
         t_spiral_start = TimingHelpers::timer();
       }
  
       // Perform locate_zeta locally first! This looks locally for
       // all not-yet-located zetas within the current spiral range.
       locate_zeta_for_local_coordinates(
         mesh_pt, external_mesh_pt, mesh_geom_obj_pt, interaction_index);
  
       // Store stats about successful locates for reporting later
       if (Doc_stats)
       {
         unsigned count_locates = 0;
         unsigned n_ext_loc = External_element_located.size();
         for (unsigned e = 0; e < n_ext_loc; e++)
         {
           unsigned n_intpt = External_element_located[e].size();
           for (unsigned ipt = 0; ipt < n_intpt; ipt++)
           {
             count_locates += External_element_located[e][ipt];
           }
         }
  
         // Store percentage of integration points successfully located.
         // Only assign if we had anything to allocte, otherwise 100%
         // (default assignment; see above) is correct
         if (tot_int != 0)
         {
           percentage_coords_located_locally.push_back(
             100.0 * double(count_locates) / double(tot_int));
         }
         else
         {
           // Had none to find so we found them all!
           percentage_coords_located_locally.push_back(100.0);
         }
       }
  
  
       // Now test whether anything needs to be broadcast elsewhere
       // (i.e. were there any failures in the locate method above?)
       // If there are, then the zetas for these failures need to be
       // broadcast...
  
       // How many zetas have we failed to find? [Note: Array is padded
       // by Dim padded entries (DBL_MAX) for each mesh]
       n_zeta_not_found =
         Flat_packed_zetas_not_found_locally.size() - Dim * n_mesh;
  
       if (Doc_timings)
       {
         t_local = TimingHelpers::timer();
         oomph_info << "CPU for local location of zeta coordinate [spiral level "
                    << i_level << "]: " << t_local - t_spiral_start << std::endl
                    << "Number of missing zetas: " << n_zeta_not_found
                    << std::endl;
       }
  
  
 #ifdef OOMPH_HAS_MPI
       // Only perform the reduction operation if there's more than one process
       if (problem_pt->communicator_pt()->nproc() > 1)
       {
         unsigned count_local_zetas = n_zeta_not_found;
         MPI_Allreduce(&count_local_zetas,
                       &n_zeta_not_found,
                       1,
                       MPI_UNSIGNED,
                       MPI_SUM,
                       problem_pt->communicator_pt()->mpi_comm());
       }
  
       // If we have missing zetas on any process
       // and the problem is distributed, we need to locate elsewhere
       if ((n_zeta_not_found != 0) &&
           (problem_pt->problem_has_been_distributed()))
       {
         // Timings
         double t_sendrecv_min = DBL_MAX;
         double t_sendrecv_max = -DBL_MAX;
         double t_sendrecv_tot = 0.0;
  
         double t_missing_min = DBL_MAX;
         double t_missing_max = -DBL_MAX;
         double t_missing_tot = 0.0;
  
         double t_send_info_min = DBL_MAX;
         double t_send_info_max = -DBL_MAX;
         double t_send_info_tot = 0.0;
  
         double t_create_halo_min = DBL_MAX;
         double t_create_halo_max = -DBL_MAX;
         double t_create_halo_tot = 0.0;
  
         // Start ring communication: Loop (number of processes - 1)
         // starting from 1. The variable iproc represents the "distance" from
         // the current process to the process for which it is attempting
         // to locate an element for the current set of not-yet-located
         // zeta coordinates
         unsigned ring_count = 0;
         for (int iproc = 1; iproc < n_proc; iproc++)
         {
           // Record time at start of loop
           if (Doc_timings)
           {
             t_loop_start = TimingHelpers::timer();
           }
  
           // Send the zeta values you haven't found to the
           // next process, receive from the previous process:
           // (Padded) Flat_packed_zetas_not_found_locally are sent
           // to next processor where they are received as
           // (padded) Received_flat_packed_zetas_to_be_found.
           send_and_receive_missing_zetas(problem_pt);
  
           if (Doc_timings)
           {
             ring_count++;
             t_sendrecv = TimingHelpers::timer();
             t_sendrecv_max =
               std::max(t_sendrecv_max, t_sendrecv - t_loop_start);
             t_sendrecv_min =
               std::min(t_sendrecv_min, t_sendrecv - t_loop_start);
             t_sendrecv_tot += (t_sendrecv - t_loop_start);
           }
  
           // Perform the locate_zeta for the new set of zetas on this process
           locate_zeta_for_missing_coordinates(
             iproc, external_mesh_pt, problem_pt, mesh_geom_obj_pt);
  
           if (Doc_timings)
           {
             t_missing = TimingHelpers::timer();
             t_missing_max = std::max(t_missing_max, t_missing - t_sendrecv);
             t_missing_min = std::min(t_missing_min, t_missing - t_sendrecv);
             t_missing_tot += (t_missing - t_sendrecv);
           }
  
           // Send any located coordinates back to the correct process, and
           // prepare to send on to the next process if necessary
           send_and_receive_located_info(iproc, external_mesh_pt, problem_pt);
  
           if (Doc_timings)
           {
             t_send_info = TimingHelpers::timer();
             t_send_info_max =
               std::max(t_send_info_max, t_send_info - t_missing);
             t_send_info_min =
               std::min(t_send_info_min, t_send_info - t_missing);
             t_send_info_tot += (t_send_info - t_missing);
           }
  
           // Create any located external halo elements on the current process
           create_external_halo_elements<EXT_ELEMENT>(
             iproc, mesh_pt, external_mesh_pt, problem_pt, interaction_index);
  
           if (Doc_timings)
           {
             t_create_halo = TimingHelpers::timer();
             t_create_halo_max =
               std::max(t_create_halo_max, t_create_halo - t_send_info);
             t_create_halo_min =
               std::min(t_create_halo_min, t_create_halo - t_send_info);
             t_create_halo_tot += (t_create_halo - t_send_info);
           }
  
           // Do we have any further locating to do or have we found
           // everything at this level of the ring communication?
           // Only perform the reduction operation if there's more than
           // one process [Note: Array is padded
           // by DIM times DBL_MAX entries for each mesh]
           n_zeta_not_found =
             Flat_packed_zetas_not_found_locally.size() - Dim * n_mesh;
  
  
 #ifdef OOMPH_HAS_MPI
           if (problem_pt->communicator_pt()->nproc() > 1)
           {
             unsigned count_local_zetas = n_zeta_not_found;
             MPI_Allreduce(&count_local_zetas,
                           &n_zeta_not_found,
                           1,
                           MPI_UNSIGNED,
                           MPI_SUM,
                           problem_pt->communicator_pt()->mpi_comm());
           }
 #endif
  
           // If  its is now zero then break out of the ring comms loop
           if (n_zeta_not_found == 0)
           {
             if (Doc_timings)
             {
               t_local = TimingHelpers::timer();
               oomph_info << "BREAK N-1: CPU for entrire spiral [spiral level "
                          << i_level << "]: " << t_local - t_spiral_start
                          << std::endl;
             }
             break;
           }
         }
  
  
         // Doc timings
         if (Doc_timings)
         {
           oomph_info << "Ring-based search continued until iteration "
                      << ring_count << " out of a maximum of "
                      << problem_pt->communicator_pt()->nproc() - 1 << "\n";
           oomph_info << "Total, av, max, min CPU for send/recv of remaining "
                         "zeta coordinates: "
                      << t_sendrecv_tot << " "
                      << t_sendrecv_tot / double(ring_count) << " "
                      << t_sendrecv_max << " " << t_sendrecv_min << "\n";
           oomph_info << "Total, av, max, min CPU for location of missing zeta "
                         "coordinates   : "
                      << t_missing_tot << " "
                      << t_missing_tot / double(ring_count) << " "
                      << t_missing_max << " " << t_missing_min << "\n";
           oomph_info << "Total, av, max, min CPU for send/recv of new element "
                         "info          : "
                      << t_send_info_tot << " "
                      << t_send_info_tot / double(ring_count) << " "
                      << t_send_info_max << " " << t_send_info_min << "\n";
           oomph_info << "Total, av, max, min CPU for local creation of "
                         "external halo objects: "
                      << t_create_halo_tot << " "
                      << t_create_halo_tot / double(ring_count) << " "
                      << t_create_halo_max << " " << t_create_halo_min << "\n";
         }
  
       } // end if for missing zetas on any process
 #endif
  
  
       // Store information about location of elements for integration points
       if (Doc_stats)
       {
         unsigned count_locates = 0;
         unsigned n_ext_loc = External_element_located.size();
         for (unsigned e = 0; e < n_ext_loc; e++)
         {
           unsigned n_intpt = External_element_located[e].size();
           for (unsigned ipt = 0; ipt < n_intpt; ipt++)
           {
             count_locates += External_element_located[e][ipt];
           }
         }
  
  
         // Store total percentage of locates so far.
         // Only assign if we had anything to allocte, otherwise 100%
         // (default assignment) is correct
         if (tot_int != 0)
         {
           percentage_coords_located_elsewhere.push_back(
             100.0 * double(count_locates) / double(tot_int));
         }
         else
         {
           // Had none to find so we found them all!
           percentage_coords_located_locally.push_back(100.0);
         }
       }
  
       // Do we have any further locating to do? If so, the remaining
       // zetas will (hopefully) be found at the next spiral level.
       // Only perform the reduction operation if there's more than one process
       // [Note: Array is padded
       // by DIM times DBL_MAX entries for each mesh]
       n_zeta_not_found =
         Flat_packed_zetas_not_found_locally.size() - Dim * n_mesh;
  
  
 #ifdef OOMPH_HAS_MPI
       if (problem_pt->communicator_pt()->nproc() > 1)
       {
         unsigned count_local_zetas = n_zeta_not_found;
         MPI_Allreduce(&count_local_zetas,
                       &n_zeta_not_found,
                       1,
                       MPI_UNSIGNED,
                       MPI_SUM,
                       problem_pt->communicator_pt()->mpi_comm());
       }
  
       // Specify max level reached for later loop
       max_level_reached = i_level + 1;
 #endif
  
       if (n_zeta_not_found == 0)
       {
         if (Doc_timings)
         {
           t_local = TimingHelpers::timer();
           oomph_info << "BREAK N: CPU for entrire spiral [spiral level "
                      << i_level << "]: " << t_local - t_spiral_start
                      << std::endl;
         }
         break;
       }
  
       if (Doc_timings)
       {
         t_local = TimingHelpers::timer();
         oomph_info << "CPU for entrire spiral [spiral level " << i_level
                    << "]: " << t_local - t_spiral_start << std::endl;
       }
  
       // Bump up spiral levels for all meshes
       i_level++;
       for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
       {
         if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version() ==
             UseRefineableBinArray)
         {
           RefineableBinArray* bin_array_pt = dynamic_cast<RefineableBinArray*>(
             mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
           bin_array_pt
             ->first_sample_point_to_actually_lookup_during_locate_zeta() =
             bin_array_pt
               ->last_sample_point_to_actually_lookup_during_locate_zeta();
           bin_array_pt
             ->last_sample_point_to_actually_lookup_during_locate_zeta() *=
             bin_array_pt
               ->multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta();
         }
         else if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version() ==
                  UseNonRefineableBinArray)
         {
           NonRefineableBinArray* bin_array_pt =
             dynamic_cast<NonRefineableBinArray*>(
               mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
  
           bin_array_pt->current_min_spiral_level() +=
             bin_array_pt->n_spiral_chunk();
           bin_array_pt->current_max_spiral_level() +=
             bin_array_pt->n_spiral_chunk();
         }
 #ifdef OOMPH_HAS_CGAL
         else if (mesh_geom_obj_pt[i_mesh]->sample_point_container_version() ==
                  UseCGALSamplePointContainer)
         {
           CGALSamplePointContainer* bin_array_pt =
             dynamic_cast<CGALSamplePointContainer*>(
               mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
           bin_array_pt
             ->first_sample_point_to_actually_lookup_during_locate_zeta() =
             bin_array_pt
               ->last_sample_point_to_actually_lookup_during_locate_zeta();
           bin_array_pt
             ->last_sample_point_to_actually_lookup_during_locate_zeta() *=
             bin_array_pt
               ->multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta();
         }
 #endif // cgal
       }
  
       // Check termination criterion for while loop
       if (mesh_geom_obj_pt[0]->sample_point_container_version() ==
           UseRefineableBinArray)
       {
         RefineableBinArray* bin_array_pt = dynamic_cast<RefineableBinArray*>(
           mesh_geom_obj_pt[0]->sample_point_container_pt());
  
         if (bin_array_pt
               ->first_sample_point_to_actually_lookup_during_locate_zeta() <=
             max_n_sample_points_of_sample_point_containers)
         {
           has_not_reached_max_level_of_search = true;
         }
         else
         {
           has_not_reached_max_level_of_search = false;
         }
       }
       else if (mesh_geom_obj_pt[0]->sample_point_container_version() ==
                UseNonRefineableBinArray)
       {
         NonRefineableBinArray* bin_array_pt =
           dynamic_cast<NonRefineableBinArray*>(
             mesh_geom_obj_pt[0]->sample_point_container_pt());
  
         if (bin_array_pt->current_max_spiral_level() < n_max_level)
         {
           has_not_reached_max_level_of_search = true;
         }
         else
         {
           has_not_reached_max_level_of_search = false;
         }
       }
 #ifdef OOMPH_HAS_CGAL
       else if (mesh_geom_obj_pt[0]->sample_point_container_version() ==
                UseCGALSamplePointContainer)
       {
         CGALSamplePointContainer* bin_array_pt =
           dynamic_cast<CGALSamplePointContainer*>(
             mesh_geom_obj_pt[0]->sample_point_container_pt());
  
         if (bin_array_pt
               ->first_sample_point_to_actually_lookup_during_locate_zeta() <=
             max_n_sample_points_of_sample_point_containers)
         {
           has_not_reached_max_level_of_search = true;
         }
         else
         {
           has_not_reached_max_level_of_search = false;
         }
       }
 #endif // cgal
     } // end of "spirals" loop
  
  
     // If we haven't found all zetas we're dead now...
     //-------------------------------------------------
     if (n_zeta_not_found != 0)
     {
       // Shout?
       if (!Accept_failed_locate_zeta_in_setup_multi_domain_interaction)
       {
         std::ostringstream error_stream;
         error_stream
           << "Multi_domain_functions::locate_zeta_for_local_coordinates()"
           << "\nhas failed ";
  
 #ifdef OOMPH_HAS_MPI
         if (problem_pt->communicator_pt()->nproc() > 1)
         {
           error_stream << " on proc: "
                        << problem_pt->communicator_pt()->my_rank() << std::endl;
         }
 #endif
         error_stream
           << "\n\n\nThis is most likely to arise because the two meshes\n"
           << "that are to be matched don't overlap perfectly or\n"
           << "because the elements are distorted and too small a \n"
           << "number of sampling points has been used when setting\n"
           << "up the bin structure.\n\n"
           << "You should try to increase the value of \n"
           << "the number of sample points defined in \n\n"
           << "  "
              "SamplePointContainerParameters::Default_nsample_points_generated_"
              "per_element"
           << "\n\n from its current value of "
           << SamplePointContainerParameters::
                Default_nsample_points_generated_per_element
           << "\n"
           << "\n\n"
           << "NOTE: You can suppress this error and \"accept failure\""
           << "      by setting the public boolean \n\n"
           << "        "
              "Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_"
              "domain_interaction\n\n"
           << "      to true. In this case, the pointers to external elements\n"
           << "      that couldn't be located will remain null\n";
  
         std::ostringstream modifier;
 #ifdef OOMPH_HAS_MPI
         if (problem_pt->communicator_pt()->nproc() > 1)
         {
           modifier << "_proc" << problem_pt->communicator_pt()->my_rank();
         }
 #endif
  
         // Loop over all meshes
         for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
         {
           // Add yet another modifier to distinguish meshes if reqd
           if (n_mesh > 1)
           {
             modifier << "_mesh" << i_mesh;
           }
  
           std::ofstream outfile;
           char filename[100];
           sprintf(
             filename, "missing_coords_mesh%s.dat", modifier.str().c_str());
           outfile.open(filename);
           unsigned nel = mesh_pt[i_mesh]->nelement();
           for (unsigned e = 0; e < nel; e++)
           {
             mesh_pt[i_mesh]->finite_element_pt(e)->output(outfile);
             // FiniteElement::output(outfile);
           }
           outfile.close();
  
           sprintf(
             filename, "missing_coords_ext_mesh%s.dat", modifier.str().c_str());
           outfile.open(filename);
           nel = external_mesh_pt->nelement();
           for (unsigned e = 0; e < nel; e++)
           {
             external_mesh_pt->finite_element_pt(e)->output(outfile);
             // FiniteElement::output(outfile);
           }
           outfile.close();
  
           BinArray* bin_array_pt = dynamic_cast<BinArray*>(
             mesh_geom_obj_pt[i_mesh]->sample_point_container_pt());
           if (bin_array_pt != 0)
           {
             sprintf(
               filename, "missing_coords_bin%s.dat", modifier.str().c_str());
             outfile.open(filename);
             bin_array_pt->output_bins(outfile);
             outfile.close();
           }
  
           sprintf(filename, "missing_coords%s.dat", modifier.str().c_str());
           outfile.open(filename);
           unsigned n = External_element_located.size();
           error_stream << "Number of unlocated elements " << n << std::endl;
           for (unsigned e = 0; e < n; e++)
           {
             unsigned n_intpt = External_element_located[e].size();
             if (n_intpt == 0)
             {
               error_stream << "Failure to locate in halo element! "
                            << "Why are we searching there?" << std::endl;
             }
             for (unsigned ipt = 0; ipt < n_intpt; ipt++)
             {
               if (External_element_located[e][ipt] == 0)
               {
                 error_stream << "Failure at element/intpt: " << e << " " << ipt
                              << std::endl;
  
                 // Cast
                 ElementWithExternalElement* el_pt =
                   dynamic_cast<ElementWithExternalElement*>(
                     mesh_pt[i_mesh]->element_pt(e));
  
                 unsigned n_dim_el = el_pt->dim();
                 Vector<double> s(n_dim_el);
                 for (unsigned i = 0; i < n_dim_el; i++)
                 {
                   s[i] = el_pt->integral_pt()->knot(ipt, i);
                 }
                 unsigned n_dim = el_pt->node_pt(0)->ndim();
                 Vector<double> x(n_dim);
                 el_pt->interpolated_x(s, x);
                 for (unsigned i = 0; i < n_dim; i++)
                 {
                   error_stream << x[i] << " ";
                   outfile << x[i] << " ";
                 }
                 error_stream << std::endl;
                 outfile << std::endl;
               }
             }
           }
           outfile.close();
         }
  
         error_stream
           << "Mesh and external mesh documented in missing_coords_mesh*.dat\n"
           << "and missing_coords_ext_mesh*.dat, respectively. Missing \n"
           << "coordinates in missing_coords*.dat\n";
         throw OomphLibError(
           error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
       }
       // Failure is deeemed to be acceptable
       else
       {
         oomph_info
           << "NOTE: Haven't found " << n_zeta_not_found
           << " external elements in \n"
           << "Multi_domain_functions::aux_setup_multi_domain_interaction(...)\n"
           << "but this deemed to be acceptable because \n"
           << "Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_"
              "domain_interaction\n"
           << "is true.\n";
       }
     }
  
  
     // Doc timings if required
     if (Doc_timings)
     {
       t_locate = TimingHelpers::timer();
       oomph_info
         << "Total CPU for location and creation of all external elements: "
         << t_locate - t_start << std::endl;
     }
  
     // Delete the geometric object representing the mesh
     for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
     {
       delete mesh_geom_obj_pt[i_mesh];
     }
  
     // Clean up all the (extern) Vectors associated with creating the
     // external storage information
     clean_up();
  
 #ifdef OOMPH_HAS_MPI
  
     // Output information about external storage if required
     if (Doc_stats)
     {
       // Report stats regarding location method
       bool comm_was_required = false;
       oomph_info << "-------------------------------------------" << std::endl;
       oomph_info << "- Cumulative percentage of locate success -" << std::endl;
       oomph_info << "--- Spiral -- Found local -- Found else ---" << std::endl;
       for (unsigned level = 0; level < max_level_reached; level++)
       {
         oomph_info << "---   " << level << "   -- "
                    << percentage_coords_located_locally[level] << " -- "
                    << percentage_coords_located_elsewhere[level] << " ---"
                    << std::endl;
         // Has communication with other processors at this level actually
         // produced any results?
         if (percentage_coords_located_elsewhere[level] >
             percentage_coords_located_locally[level])
         {
           comm_was_required = true;
         }
       }
       oomph_info << "-------------------------------------------" << std::endl;
  
  
       // No need for any of this malaki if we're not running in parallel
       if (problem_pt->communicator_pt()->nproc() > 1)
       {
         // Initialise to indicate that none of the zetas required
         // on this processor were located through parallel ring search,
         // i.e. comm was not required and we could have done some
         // more local searching first
         oomph_info << std::endl;
         oomph_info << "ASSESSMENT OF NEED FOR PARALLEL SEARCH: \n";
         oomph_info << "=======================================\n";
         unsigned status = 0;
         if (comm_was_required)
         {
           oomph_info << "- Ring-based parallel search did successfully locate "
                         "zetas on proc "
                      << my_rank << std::endl;
           status = 1;
         }
         else
         {
           if (max_level_reached > 1)
           {
             oomph_info
               << "- Ring-based parallel search did NOT locate zetas on proc"
               << my_rank << std::endl;
           }
           else
           {
             oomph_info
               << "- No ring-based parallel search was performed on proc"
               << my_rank << std::endl;
           }
         }
  
         // Allreduce to check if anyone has benefitted from parallel ring
         // search
         unsigned overall_status = 0;
         MPI_Allreduce(&status,
                       &overall_status,
                       1,
                       MPI_UNSIGNED,
                       MPI_MAX,
                       problem_pt->communicator_pt()->mpi_comm());
  
         // Report of mpi was useful to anyone
         if (overall_status == 1)
         {
           oomph_info << "- Ring-based, parallel search did succesfully\n";
           oomph_info << "  locate zetas on at least one other proc, so it\n";
           oomph_info << "  was worthwhile.\n";
           oomph_info << std::endl;
         }
         else
         {
           if (max_level_reached > 1)
           {
             oomph_info
               << "- Ring-based, parallel search did NOT locate zetas\n";
             oomph_info << "  on ANY other procs, i.e it was useless.\n";
             oomph_info
               << "  --> We should really have done more local search\n";
             oomph_info
               << "   by reducing number of bins, or doing more spirals\n";
             oomph_info << "   in one go before initiating parallel search.\n";
             oomph_info << std::endl;
           }
           else
           {
             oomph_info << "- No ring-based, parallel search was performed\n";
             oomph_info << "  or necessary. Perfect!\n";
             oomph_info << std::endl;
           }
         }
       }
  
       // How many external elements does the external mesh have now?
       oomph_info << "------------------------------------------" << std::endl;
       oomph_info << "External mesh: I have " << external_mesh_pt->nelement()
                  << " elements, and " << std::endl
                  << external_mesh_pt->nexternal_halo_element()
                  << " external halo elements, "
                  << external_mesh_pt->nexternal_haloed_element()
                  << " external haloed elements" << std::endl;
  
       // How many external nodes does each submesh have now?
       oomph_info << "------------------------------------------" << std::endl;
       oomph_info << "External mesh: I have " << external_mesh_pt->nnode()
                  << " nodes, and " << std::endl
                  << external_mesh_pt->nexternal_halo_node()
                  << " external halo nodes, "
                  << external_mesh_pt->nexternal_haloed_node()
                  << " external haloed nodes" << std::endl;
       oomph_info << "------------------------------------------" << std::endl;
     }
  
     // Output further information about (external) halo(ed)
     // elements and nodes if required
     if (Doc_full_stats)
     {
       // How many elements does this submesh have for each of the processors?
       for (int iproc = 0; iproc < n_proc; iproc++)
       {
         oomph_info << "----------------------------------------" << std::endl;
         oomph_info << "With process " << iproc << " there are "
                    << external_mesh_pt->nroot_halo_element(iproc)
                    << " root halo elements, and "
                    << external_mesh_pt->nroot_haloed_element(iproc)
                    << " root haloed elements" << std::endl
                    << "and there are "
                    << external_mesh_pt->nexternal_halo_element(iproc)
                    << " external halo elements, and "
                    << external_mesh_pt->nexternal_haloed_element(iproc)
                    << " external haloed elements." << std::endl;
  
         oomph_info << "----------------------------------------" << std::endl;
         oomph_info << "With process " << iproc << " there are "
                    << external_mesh_pt->nhalo_node(iproc) << " halo nodes, and "
                    << external_mesh_pt->nhaloed_node(iproc) << " haloed nodes"
                    << std::endl
                    << "and there are "
                    << external_mesh_pt->nexternal_halo_node(iproc)
                    << " external halo nodes, and "
                    << external_mesh_pt->nexternal_haloed_node(iproc)
                    << " external haloed nodes." << std::endl;
       }
       oomph_info << "-----------------------------------------" << std::endl
                  << std::endl;
     }
  
 #endif
  
     // Doc timings if required
     if (Doc_timings)
     {
       t_end = TimingHelpers::timer();
       oomph_info << "CPU for (one way) aux_setup_multi_domain_interaction: "
                  << t_end - t_start << std::endl;
     }
  
   } // end of aux_setup_multi_domain_interaction

References Accept_failed_locate_zeta_in_setup_multi_domain_interaction, clean_up(), oomph::Problem::communicator_pt(), NonRefineableBinArray::current_max_spiral_level(), NonRefineableBinArray::current_min_spiral_level(), Global_Variables::Dim, oomph::FiniteElement::dim(), oomph::Missing_masters_functions::Doc_full_stats, oomph::Missing_masters_functions::Doc_stats, oomph::Missing_masters_functions::Doc_timings, e(), oomph::Mesh::element_pt(), External_element_located, oomph::ElementWithExternalElement::external_element_pt(), MergeRestartFiles::filename, oomph::Mesh::finite_element_pt(), RefineableBinArray::first_sample_point_to_actually_lookup_during_locate_zeta(), Flat_packed_zetas_not_found_locally, i, RefineableBinArray::initial_last_sample_point_to_actually_lookup_during_locate_zeta(), oomph::ElementWithExternalElement::initialise_external_element_storage(), oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_x(), oomph::Integral::knot(), RefineableBinArray::last_sample_point_to_actually_lookup_during_locate_zeta(), locate_zeta_for_local_coordinates(), max, BinArray::max_bin_dimension(), min, RefineableBinArray::multiplier_for_max_sample_point_to_actually_lookup_during_locate_zeta(), oomph::OomphCommunicator::my_rank(), n, NonRefineableBinArray::n_spiral_chunk(), oomph::Node::ndim(), oomph::Mesh::nelement(), oomph::Mesh::nnode(), oomph::FiniteElement::node_pt(), oomph::OomphCommunicator::nproc(), oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::FiniteElement::output(), BinArray::output_bins(), s, plotPSD::t, oomph::TimingHelpers::timer(), RefineableBinArray::total_number_of_sample_points_computed_recursively(), Use_bulk_element_as_external, oomph::UseNonRefineableBinArray, oomph::UseRefineableBinArray, and plotDoE::x.

◆ aux_setup_multi_domain_interaction() [2/2]

template<class EXT_ELEMENT , class GEOM_OBJECT >

void oomph::Multi_domain_functions::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.

This routine calls the locate_zeta routine (simultaneously on each processor for each individual processor's element set if necessary) and sets up the external (halo) element and node storage as necessary. The locate_zeta procedure here works for all multi-domain problems where either two meshes occupy the same physical space but have differing element types (e.g. a Boussinesq convection problem where AdvectionDiffusion elements interact with Navier-Stokes type elements) or two meshes interact along some boundary of the external mesh, represented by a "face mesh", such as an FSI problem.

   {
     // Convert to vector-based storage
     Vector<Mesh*> mesh_pt_vector(1);
     mesh_pt_vector[0] = mesh_pt;
     Vector<Mesh*> external_face_mesh_pt_vector(1);
     external_face_mesh_pt_vector[0] = external_face_mesh_pt;
  
     // Call vector-based version
     aux_setup_multi_domain_interaction<EXT_ELEMENT, GEOM_OBJECT>(
       problem_pt,
       mesh_pt_vector,
       external_mesh_pt,
       interaction_index,
       external_face_mesh_pt_vector);
  
   } // end of aux_setup_multi_domain_interaction

◆ clean_up()

void oomph::Multi_domain_functions::clean_up ( )

Helper function that clears all the information used during the external storage creation

Helper function that clears all the intermediate information used during the external storage creation at the end of the procedure

     {
       Flat_packed_zetas_not_found_locally.clear();
       Received_flat_packed_zetas_to_be_found.clear();
       Proc_id_plus_one_of_external_element.clear();
       Located_element_status.clear();
       Flat_packed_located_coordinates.clear();
       Flat_packed_doubles.clear();
       Flat_packed_unsigneds.clear();
       External_element_located.clear();
     }

References External_element_located, Flat_packed_doubles, Flat_packed_located_coordinates, Flat_packed_unsigneds, Flat_packed_zetas_not_found_locally, Located_element_status, Proc_id_plus_one_of_external_element, and Received_flat_packed_zetas_to_be_found.

Referenced by aux_setup_multi_domain_interaction().

◆ first_closer_than_second()

bool oomph::Multi_domain_functions::first_closer_than_second	(	const std::pair< FiniteElement *, Vector< double >> &	p1,
		const std::pair< FiniteElement *, Vector< double >> &	p2
	)

Comparison function for sorting entries in bin: Returns true if point identified by p1 (comprising pointer to finite element and vector of local coordinates within that element) is closer to Zeta_coords_for_further_away_comparison than p2

     {
       // First point
       FiniteElement* el_pt = p1.first;
       Vector<double> s(p1.second);
       Vector<double> zeta(Dim);
       el_pt->interpolated_zeta(s, zeta);
       double dist_squared1 = 0.0;
       for (unsigned i = 0; i < Dim; i++)
       {
         dist_squared1 +=
           (zeta[i] - Zeta_coords_for_further_away_comparison[i]) *
           (zeta[i] - Zeta_coords_for_further_away_comparison[i]);
       }
  
       // Second point
       el_pt = p2.first;
       s = p2.second;
       el_pt->interpolated_zeta(s, zeta);
       double dist_squared2 = 0.0;
       for (unsigned i = 0; i < Dim; i++)
       {
         dist_squared2 +=
           (zeta[i] - Zeta_coords_for_further_away_comparison[i]) *
           (zeta[i] - Zeta_coords_for_further_away_comparison[i]);
       }
  
       // Which one is further
       if (dist_squared1 < dist_squared2)
       {
         return true;
       }
       else
       {
         return false;
       }
     }

References Dim, i, oomph::FiniteElement::interpolated_zeta(), p1, s, Eigen::zeta(), and Zeta_coords_for_further_away_comparison.

◆ get_dim_helper()

void oomph::Multi_domain_functions::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 checks they are the same). Stores result in Dim.

Helper function that computes the dimension of the elements within each of the specified meshes (and checks they are the same) Stores result in Dim.

     {
 #ifdef OOMPH_HAS_MPI
       // Storage for number of processors, current process and communicator
       OomphCommunicator* comm_pt = problem_pt->communicator_pt();
 #endif
  
       // Extract the element dimensions from the first element of each mesh
       unsigned mesh_dim = 0;
       if (mesh_pt->nelement() > 0)
       {
         mesh_dim = dynamic_cast<FiniteElement*>(mesh_pt->element_pt(0))->dim();
       }
       unsigned external_mesh_dim = 0;
       if (external_mesh_pt->nelement() > 0)
       {
         external_mesh_dim =
           dynamic_cast<FiniteElement*>(external_mesh_pt->element_pt(0))->dim();
       }
  
       // Need to do an Allreduce
 #ifdef OOMPH_HAS_MPI
       int n_proc = comm_pt->nproc();
       if (n_proc > 1)
       {
         unsigned mesh_dim_reduce;
         MPI_Allreduce(&mesh_dim,
                       &mesh_dim_reduce,
                       1,
                       MPI_UNSIGNED,
                       MPI_MAX,
                       comm_pt->mpi_comm());
         mesh_dim = mesh_dim_reduce;
  
         unsigned external_mesh_dim_reduce;
         MPI_Allreduce(&external_mesh_dim,
                       &external_mesh_dim_reduce,
                       1,
                       MPI_UNSIGNED,
                       MPI_MAX,
                       comm_pt->mpi_comm());
         external_mesh_dim = external_mesh_dim_reduce;
       }
 #endif
  
       // Check the dimensions are the same!
       if (mesh_dim != external_mesh_dim)
       {
         std::ostringstream error_stream;
         error_stream << "The elements within the two meshes do not\n"
                      << "have the same dimension, so the multi-domain\n"
                      << "method will not work.\n"
                      << "For the mesh, dim=" << mesh_dim
                      << ", and the external mesh, dim=" << external_mesh_dim
                      << "\n";
         throw OomphLibError(
           error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
       }
  
       // Set dimension
       Dim = mesh_dim;
     }

References oomph::Problem::communicator_pt(), Dim, oomph::Mesh::element_pt(), oomph::Mesh::nelement(), oomph::OomphCommunicator::nproc(), OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

Referenced by setup_multi_domain_interaction().

◆ locate_zeta_for_local_coordinates()

void oomph::Multi_domain_functions::locate_zeta_for_local_coordinates	(	const Vector< Mesh * > &	mesh_pt,
		Mesh *const &	external_mesh_pt,
		Vector< MeshAsGeomObject * > &	mesh_geom_obj_pt,
		const unsigned &	interaction_index
	)

locate zeta for current set of "local" coordinates vector-based version

Helper function to locate "local" zeta coordinates This is the vector-based version which operates simultaenously on the meshes contained in the Vectors.

     {
       // Flush storage for zetas not found locally
       Flat_packed_zetas_not_found_locally.resize(0);
  
       // Number of meshes
       unsigned n_mesh = mesh_pt.size();
  
 #ifdef PARANOID
       if (mesh_geom_obj_pt.size() != n_mesh)
       {
         std::ostringstream error_stream;
         error_stream << "Sizes of mesh_geom_obj_pt [ "
                      << mesh_geom_obj_pt.size() << " ] and "
                      << "mesh_pt [ " << n_mesh << " ] don't match.\n";
         throw OomphLibError(
           error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
       }
 #endif
  
       // Element counter
       unsigned e_count = 0;
  
       // Loop over meshes
       for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
       {
         // Number of local elements
         unsigned n_element = mesh_pt[i_mesh]->nelement();
  
         // Loop over this processor's elements
         for (unsigned e = 0; e < n_element; e++)
         {
           ElementWithExternalElement* el_pt =
             dynamic_cast<ElementWithExternalElement*>(
               mesh_pt[i_mesh]->element_pt(e));
 #ifdef OOMPH_HAS_MPI
           // Only visit non-halo elements -- we're not setting up external
           // elements for on-halos!
           if (!el_pt->is_halo())
 #endif
           {
             // Find number of Gauss points and element dimension
             unsigned n_intpt = el_pt->integral_pt()->nweight();
             unsigned el_dim = el_pt->dim();
  
  
 #ifdef PARANOID
             if (el_dim != Dim)
             {
               std::ostringstream error_stream;
               error_stream << "Dimension of element " << el_dim
                            << " is not consitent with dimension assumed \n"
                            << " in multidomain namespace, " << Dim << std::endl;
               throw OomphLibError(error_stream.str(),
                                   OOMPH_CURRENT_FUNCTION,
                                   OOMPH_EXCEPTION_LOCATION);
             }
 #endif
  
             // Set storage for local and global coordinates
             Vector<double> s_local(el_dim);
             Vector<double> x_global(el_dim);
  
             // Loop over integration points
             for (unsigned ipt = 0; ipt < n_intpt; ipt++)
             {
               // Has this integration point been done yet?
               if (External_element_located[e_count][ipt] == 0)
               {
                 // Get local coordinates
                 for (unsigned i = 0; i < el_dim; i++)
                 {
                   s_local[i] = el_pt->integral_pt()->knot(ipt, i);
                 }
                 // Interpolate to global coordinates
                 el_pt->interpolated_zeta(s_local, x_global);
  
                 // Storage for geometric object and its local coordinates
                 GeomObject* sub_geom_obj_pt = 0;
                 Vector<double> s_ext(el_dim);
                 mesh_geom_obj_pt[i_mesh]->locate_zeta(
                   x_global, sub_geom_obj_pt, s_ext);
  
                 // Has the required element been located?
                 if (sub_geom_obj_pt != 0)
                 {
                   // The required element has been located
                   // The located coordinates have the same dimension as the bulk
                   GeneralisedElement* source_el_pt;
                   Vector<double> s_source(el_dim);
  
                   // Is the bulk element the actual external element?
                   if (!Use_bulk_element_as_external)
                   {
                     // Use the object directly (it must be a finite element)
                     source_el_pt =
                       dynamic_cast<FiniteElement*>(sub_geom_obj_pt);
                     s_source = s_ext;
                   }
                   else
                   {
                     // Cast to a FaceElement and use the bulk element
                     FaceElement* face_el_pt =
                       dynamic_cast<FaceElement*>(sub_geom_obj_pt);
                     source_el_pt = face_el_pt->bulk_element_pt();
  
                     // Need to resize the located coordinates to have the same
                     // dimension as the bulk element
                     s_source.resize(
                       dynamic_cast<FiniteElement*>(source_el_pt)->dim());
  
                     // Translate the returned local coords into the bulk element
                     face_el_pt->get_local_coordinate_in_bulk(s_ext, s_source);
                   }
  
                   // Check if it's a halo; if it is then the non-halo equivalent
                   // needs to be located from another processor (unless we
                   // accept halo elements as external elements)
 #ifdef OOMPH_HAS_MPI
                   if (Allow_use_of_halo_elements_as_external_elements ||
                       (!source_el_pt->is_halo()))
 #endif
                   {
                     // Need to cast to a FiniteElement
                     FiniteElement* source_finite_el_pt =
                       dynamic_cast<FiniteElement*>(source_el_pt);
  
                     // Set the external element pointer and local coordinates
                     el_pt->external_element_pt(interaction_index, ipt) =
                       source_finite_el_pt;
                     el_pt->external_element_local_coord(interaction_index,
                                                         ipt) = s_source;
  
                     // Set the lookup array to 1/true
                     External_element_located[e_count][ipt] = 1;
                   }
 #ifdef OOMPH_HAS_MPI
                   // located element is halo and we're not accepting haloes
                   // obviously only makes sense in mpi mode...
                   else
                   {
                     // Add required information to arrays
                     for (unsigned i = 0; i < el_dim; i++)
                     {
                       Flat_packed_zetas_not_found_locally.push_back(
                         x_global[i]);
                     }
                   }
 #endif
                 }
                 else
                 {
                   // Search has failed then add the required information to the
                   // arrays which need to be sent to the other processors so
                   // that they can perform the locate_zeta
  
                   // Add this global coordinate to the LOCAL zeta array
                   for (unsigned i = 0; i < el_dim; i++)
                   {
                     Flat_packed_zetas_not_found_locally.push_back(x_global[i]);
                   }
                 }
               }
             } // end loop over integration points
           } // end for halo
  
           // Bump up counter for all elements
           e_count++;
  
         } // end loop over local elements
  
         // Mark end of mesh data in flat packed array
         for (unsigned i = 0; i < Dim; i++)
         {
           Flat_packed_zetas_not_found_locally.push_back(DBL_MAX);
         }
  
       } // end of loop over meshes
     }

References Allow_use_of_halo_elements_as_external_elements, oomph::FaceElement::bulk_element_pt(), oomph::FiniteElement::dim(), Dim, e(), el_dim, oomph::ElementWithExternalElement::external_element_local_coord(), External_element_located, oomph::ElementWithExternalElement::external_element_pt(), Flat_packed_zetas_not_found_locally, oomph::FaceElement::get_local_coordinate_in_bulk(), i, oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_zeta(), oomph::Integral::knot(), oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and Use_bulk_element_as_external.

Referenced by aux_setup_multi_domain_interaction().

◆ setup_bulk_elements_adjacent_to_face_mesh() [1/2]

template<class BULK_ELEMENT , unsigned DIM>

void oomph::Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh	(	Problem *	problem_pt,
		const unsigned &	boundary_in_bulk_mesh,
		Mesh *const &	bulk_mesh_pt,
		Mesh *const &	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 bulk elements next to the boundary_in_bulk_mesh -th boundary of the mesh pointed to by bulk_mesh_pt. The FaceElements must be derived from the ElementWithExternalElement base class and the adjacent bulk elements are stored as their external elements.

   {
     // Convert to vector-based storage
     Vector<unsigned> boundary_in_bulk_mesh_vect(1);
     boundary_in_bulk_mesh_vect[0] = boundary_in_bulk_mesh;
     Vector<Mesh*> face_mesh_pt_vect(1);
     face_mesh_pt_vect[0] = face_mesh_pt;
  
     // Call vector-based version
     setup_bulk_elements_adjacent_to_face_mesh<BULK_ELEMENT, DIM>(
       problem_pt,
       boundary_in_bulk_mesh_vect,
       bulk_mesh_pt,
       face_mesh_pt_vect,
       interaction);
   }

◆ setup_bulk_elements_adjacent_to_face_mesh() [2/2]

template<class BULK_ELEMENT , unsigned DIM>

void oomph::Multi_domain_functions::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`
	)

/ Templated helper functions for multi-domain methods using locate_zeta

Identify the FaceElements (stored in the mesh pointed to by face_mesh_pt) that are adjacent to the bulk elements next to the boundary_in_bulk_mesh -th boundary of the mesh pointed to by bulk_mesh_pt. The FaceElements must be derived from the ElementWithExternalElement base class and the adjacent bulk elements are stored as their external elements.

This is the vector-based version which deals with multiple bulk mesh boundaries at the same time.

   {
     unsigned n_mesh = boundary_in_bulk_mesh.size();
  
 #ifdef PARANOID
     // Check sizes match
     if (boundary_in_bulk_mesh.size() != face_mesh_pt.size())
     {
       std::ostringstream error_message;
       error_message << "Sizes of vector of boundary ids in bulk mesh ("
                     << boundary_in_bulk_mesh.size()
                     << ") and vector of pointers\n"
                     << "to FaceElements (" << face_mesh_pt.size()
                     << " doesn't match.\n";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Create face meshes adjacent to the bulk mesh's b-th boundary.
     // Each face mesh consists of FaceElements that may also be
     // interpreted as GeomObjects
     Vector<Mesh*> bulk_face_mesh_pt(n_mesh);
  
     // Loop over all meshes
     for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
     {
       bulk_face_mesh_pt[i_mesh] = new Mesh;
       bulk_mesh_pt
         ->template build_face_mesh<BULK_ELEMENT, FaceElementAsGeomObject>(
           boundary_in_bulk_mesh[i_mesh], bulk_face_mesh_pt[i_mesh]);
  
       // Loop over these new face elements and tell them the boundary number
       // from the bulk mesh -- this is required to they can
       // get access to the boundary coordinates!
       unsigned n_face_element = bulk_face_mesh_pt[i_mesh]->nelement();
       for (unsigned e = 0; e < n_face_element; e++)
       {
         // Cast the element pointer to the correct thing!
         FaceElementAsGeomObject<BULK_ELEMENT>* el_pt =
           dynamic_cast<FaceElementAsGeomObject<BULK_ELEMENT>*>(
             bulk_face_mesh_pt[i_mesh]->element_pt(e));
  
         // Set bulk boundary number
         el_pt->set_boundary_number_in_bulk_mesh(boundary_in_bulk_mesh[i_mesh]);
  
         // Doc?
         if (Doc_boundary_coordinate_file.is_open())
         {
           Vector<double> s(DIM - 1);
           Vector<double> zeta(DIM - 1);
           Vector<double> x(DIM);
           unsigned n_plot = 5;
           Doc_boundary_coordinate_file << el_pt->tecplot_zone_string(n_plot);
  
           // Loop over plot points
           unsigned num_plot_points = el_pt->nplot_points(n_plot);
           for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
           {
             // Get local coordinates of plot point
             el_pt->get_s_plot(iplot, n_plot, s);
             el_pt->interpolated_zeta(s, zeta);
             el_pt->interpolated_x(s, x);
             for (unsigned i = 0; i < DIM; i++)
             {
               Doc_boundary_coordinate_file << x[i] << " ";
             }
             for (unsigned i = 0; i < DIM - 1; i++)
             {
               Doc_boundary_coordinate_file << zeta[i] << " ";
             }
             Doc_boundary_coordinate_file << std::endl;
           }
           el_pt->write_tecplot_zone_footer(Doc_boundary_coordinate_file,
                                            n_plot);
         }
       }
  
       // Now sort the elements based on the boundary coordinates.
       // This may allow a faster implementation of the locate_zeta
       // function for the MeshAsGeomObject representation of this
       // mesh, but also creates a unique ordering of the elements
       std::sort(bulk_face_mesh_pt[i_mesh]->element_pt().begin(),
                 bulk_face_mesh_pt[i_mesh]->element_pt().end(),
                 CompareBoundaryCoordinate<BULK_ELEMENT>());
     } // end of loop over meshes
  
  
     // Setup the interactions for this problem using the surface mesh
     // on the fluid mesh.  The interaction parameter in this instance is
     // given by the "face" parameter.
     Multi_domain_functions::setup_multi_domain_interaction<
       BULK_ELEMENT,
       FaceElementAsGeomObject<BULK_ELEMENT>>(
       problem_pt, face_mesh_pt, bulk_mesh_pt, bulk_face_mesh_pt, interaction);
  
  
     // Loop over all meshes to clean up
     for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
     {
       unsigned n_face_element = bulk_face_mesh_pt[i_mesh]->nelement();
  
       // The MeshAsGeomObject has already been deleted (in set_external_storage)
  
       // Must be careful with the FaceMesh, because we cannot delete the nodes
       // Loop over the elements backwards (paranoid!) and delete them
       for (unsigned e = n_face_element; e > 0; e--)
       {
         delete bulk_face_mesh_pt[i_mesh]->element_pt(e - 1);
         bulk_face_mesh_pt[i_mesh]->element_pt(e - 1) = 0;
       }
       // Now clear all element and node storage (should maybe fine-grain this)
       bulk_face_mesh_pt[i_mesh]->flush_element_and_node_storage();
  
       // Finally delete the mesh
       delete bulk_face_mesh_pt[i_mesh];
  
     } // end of loop over meshes
   }

References DIM, Doc_boundary_coordinate_file, e(), Eigen::placeholders::end, oomph::FiniteElement::get_s_plot(), i, oomph::FaceElement::interpolated_x(), oomph::FiniteElement::interpolated_zeta(), oomph::FiniteElement::nplot_points(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, oomph::FaceElement::set_boundary_number_in_bulk_mesh(), setup_multi_domain_interaction(), oomph::FiniteElement::tecplot_zone_string(), oomph::FiniteElement::write_tecplot_zone_footer(), plotDoE::x, and Eigen::zeta().

Referenced by oomph::FSI_functions::setup_fluid_load_info_for_solid_elements(), PressureWaveFSIProblem< FLUID_ELEMENT, SOLID_ELEMENT >::setup_fsi(), CoatedDiskProblem< ELASTICITY_ELEMENT, HELMHOLTZ_ELEMENT >::setup_interaction(), CoatedSphereProblem< ELASTICITY_ELEMENT, HELMHOLTZ_ELEMENT >::setup_interaction(), and oomph::FSI_functions::setup_solid_elements_for_displacement_bc().

◆ setup_multi_domain_interaction() [1/3]

template<class EXT_ELEMENT , class FACE_ELEMENT_GEOM_OBJECT >

void oomph::Multi_domain_functions::setup_multi_domain_interaction	(	Problem *	problem_pt,
		const Vector< Mesh * > &	mesh_pt,
		Mesh *const &	external_mesh_pt,
		const Vector< Mesh * > &	external_face_mesh_pt,
		const unsigned &	interaction_index = `0`
	)

Function to set up the one-way multi-domain interaction for FSI-like problems.

mesh_pt points to the mesh of ElemenWithExternalElements for which the interaction is set up. In an FSI example, this mesh would contain the FSIWallElements (either beam/shell elements or the FSISolidTractionElements that apply the traction to a "bulk" solid mesh that is loaded by the fluid.)
external_mesh_pt points to the mesh that contains the elements of type EXT_ELEMENT that provide the "source" for the ElementWithExternalElements. In an FSI example, this mesh would contain the "bulk" fluid elements.
external_face_mesh_pt points to the mesh of FaceElements attached to the external_mesh_pt. The mesh pointed to by external_face_mesh_pt has the same dimension as mesh_pt. The elements contained in external_face_mesh_pt are of type FACE_ELEMENT_GEOM_OBJECT. In an FSI example, these elements are usually the FaceElementAsGeomObjects (templated by the type of the "bulk" fluid elements to which they are attached) that define the FSI boundary of the fluid domain.
The interaction_index parameter defaults to zero and must otherwise be set by the user if there is more than one mesh that provides "external elements" for the Mesh pointed to by mesh_pt (e.g. in the case when a beam or shell structure is loaded by fluid from both sides.)

This is the vector-based version which operates simultaneously on the meshes contained in the Vector arguments.

Function to set up the one-way multi-domain interaction for FSI-like problems.

mesh_pt points to the mesh of ElemenWithExternalElements for which the interaction is set up. In an FSI example, this mesh would contain the FSIWallElements (either beam/shell elements or the FSISolidTractionElements that apply the traction to a "bulk" solid mesh that is loaded by the fluid.)
external_mesh_pt points to the mesh that contains the elements of type EXT_ELEMENT that provide the "source" for the ElementWithExternalElements. In an FSI example, this mesh would contain the "bulk" fluid elements.
external_face_mesh_pt points to the mesh of FaceElements attached to the external_mesh_pt. The mesh pointed to by external_face_mesh_pt has the same dimension as mesh_pt. The elements contained in external_face_mesh_pt are of type FACE_ELEMENT_GEOM_OBJECT. In an FSI example, these elements are usually the FaceElementAsGeomObjects (templated by the type of the "bulk" fluid elements to which they are attached) that define the FSI boundary of the fluid domain.
The interaction_index parameter defaults to zero and must otherwise be set by the user if there is more than one mesh that provides "external elements" for the Mesh pointed to by mesh_pt (e.g. in the case when a beam or shell structure is loaded by fluid from both sides.)

Vector-based version operates simultaneously on the meshes contained in the vectors.

   {
     // How many meshes do we have?
     unsigned n_mesh = mesh_pt.size();
  
 #ifdef PARANOID
     if (external_face_mesh_pt.size() != n_mesh)
     {
       std::ostringstream error_stream;
       error_stream << "Sizes of external_face_mesh_pt [ "
                    << external_face_mesh_pt.size() << " ] and "
                    << "mesh_pt [ " << n_mesh << " ] don't match.\n";
       throw OomphLibError(
         error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Bail out?
     if (n_mesh == 0) return;
  
     // Bulk elements must be external elements in this case
     Use_bulk_element_as_external = true;
  
     // Call the auxiliary routine with GEOM_OBJECT=FACE_ELEMENT_GEOM_OBJECT
     // and EL_DIM_LAG=Dim, EL_DIM_EUL=Dim+1. Use first mesh only.
     get_dim_helper(problem_pt, mesh_pt[0], external_face_mesh_pt[0]);
  
  
 #ifdef PARANOID
     // Check consitency
     unsigned old_dim = Dim;
     for (unsigned i = 1; i < n_mesh; i++)
     {
       // Set up Dim again
       get_dim_helper(problem_pt, mesh_pt[i], external_face_mesh_pt[i]);
  
       if (Dim != old_dim)
       {
         std::ostringstream error_stream;
         error_stream << "Inconsistency: Mesh  " << i << " has Dim=" << Dim
                      << "while mesh 0 has Dim=" << old_dim << std::endl;
         throw OomphLibError(
           error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
       }
     }
 #endif
  
     if (Dim > 2)
     {
       std::ostringstream error_stream;
       error_stream << "The elements within the two interacting meshes have a\n"
                    << "dimension not equal to 1 or 2.\n"
                    << "The multi-domain method will not work in this case.\n"
                    << "The dimension is: " << Dim << "\n";
       throw OomphLibError(
         error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
  
     // Now do the actual work for all meshes simultaneously
     aux_setup_multi_domain_interaction<EXT_ELEMENT, FACE_ELEMENT_GEOM_OBJECT>(
       problem_pt,
       mesh_pt,
       external_mesh_pt,
       interaction_index,
       external_face_mesh_pt);
   }

References Global_Variables::Dim, get_dim_helper(), i, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and Use_bulk_element_as_external.

◆ setup_multi_domain_interaction() [2/3]

template<class EXT_ELEMENT >

void oomph::Multi_domain_functions::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 mesh_pt and external_mesh_pt occupy the same physical space, and the elements in external_mesh_pt act as "external elements" for the ElementWithExternalElements in mesh_pt (but not vice versa):

mesh_pt points to the mesh of ElemenWithExternalElements for which the interaction is set up.
external_mesh_pt points to the mesh that contains the elements of type EXT_ELEMENT that act as "external elements" for the ElementWithExternalElements in \ mesh_pt.
The interaction_index parameter defaults to zero and must be otherwise set by the user if there is more than one mesh that provides sources for the Mesh pointed to by mesh_pt.

Function to set up the one-way multi-domain interaction for problems where the meshes pointed to by mesh_pt and external_mesh_pt occupy the same physical space, and the elements in external_mesh_pt act as "external elements" for the ElementWithExternalElements in mesh_pt (but not vice versa):

mesh_pt points to the mesh of ElemenWithExternalElements for which the interaction is set up.
external_mesh_pt points to the mesh that contains the elements of type EXT_ELEMENT that act as "external elements" for the ElementWithExternalElements in mesh_pt.
The interaction_index parameter defaults to zero and must be otherwise set by the user if there is more than one mesh that provides sources for the Mesh pointed to by mesh_pt.

   {
     // Bulk elements must not be external elements in this case
     Use_bulk_element_as_external = false;
  
     // Call the auxiliary function with GEOM_OBJECT=EXT_ELEMENT
     // and EL_DIM_EUL=EL_DIM_LAG=dimension returned from helper function
     get_dim_helper(problem_pt, mesh_pt, external_mesh_pt);
  
     if (Dim > 3)
     {
       std::ostringstream error_stream;
       error_stream << "The elements within the two interacting meshes have a\n"
                    << "dimension not equal to 1, 2 or 3.\n"
                    << "The multi-domain method will not work in this case.\n"
                    << "The dimension is: " << Dim << "\n";
       throw OomphLibError(
         error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
  
     // Wrapper for each dimension (template parameter)
     aux_setup_multi_domain_interaction<EXT_ELEMENT, EXT_ELEMENT>(
       problem_pt, mesh_pt, external_mesh_pt, interaction_index);
   }

References Global_Variables::Dim, get_dim_helper(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and Use_bulk_element_as_external.

Referenced by MortaringValidationProblem< ELEMENT, NON_MORTAR_ELEMENT >::actions_after_adapt(), oomph::RefineableTetgenMesh< ELEMENT >::adapt(), ContactProblem< ELEMENT >::complete_problem_setup(), MortaringValidationProblem< ELEMENT, NON_MORTAR_ELEMENT >::MortaringValidationProblem(), oomph::ProjectionProblem< PROJECTABLE_ELEMENT >::project(), and setup_bulk_elements_adjacent_to_face_mesh().

◆ setup_multi_domain_interaction() [3/3]

template<class EXT_ELEMENT , class FACE_ELEMENT_GEOM_OBJECT >

void oomph::Multi_domain_functions::setup_multi_domain_interaction	(	Problem *	problem_pt,
		Mesh *const &	mesh_pt,
		Mesh *const &	external_mesh_pt,
		Mesh *const &	external_face_mesh_pt,
		const unsigned &	interaction_index = `0`
	)

Function to set up the one-way multi-domain interaction for FSI-like problems.

mesh_pt points to the mesh of ElemenWithExternalElements for which the interaction is set up. In an FSI example, this mesh would contain the FSIWallElements (either beam/shell elements or the FSISolidTractionElements that apply the traction to a "bulk" solid mesh that is loaded by the fluid.)
external_mesh_pt points to the mesh that contains the elements of type EXT_ELEMENT that provide the "source" for the ElementWithExternalElements. In an FSI example, this mesh would contain the "bulk" fluid elements.
external_face_mesh_pt points to the mesh of FaceElements attached to the external_mesh_pt. The mesh pointed to by external_face_mesh_pt has the same dimension as mesh_pt. The elements contained in external_face_mesh_pt are of type FACE_ELEMENT_GEOM_OBJECT. In an FSI example, these elements are usually the FaceElementAsGeomObjects (templated by the type of the "bulk" fluid elements to which they are attached) that define the FSI boundary of the fluid domain.
The interaction_index parameter defaults to zero and must otherwise be set by the user if there is more than one mesh that provides "external elements" for the Mesh pointed to by mesh_pt (e.g. in the case when a beam or shell structure is loaded by fluid from both sides.)

   {
     // Bulk elements must be external elements in this case
     Use_bulk_element_as_external = true;
  
     // Call the auxiliary routine with GEOM_OBJECT=FACE_ELEMENT_GEOM_OBJECT
     // and EL_DIM_LAG=Dim, EL_DIM_EUL=Dim+1
     get_dim_helper(problem_pt, mesh_pt, external_face_mesh_pt);
  
     if (Dim > 2)
     {
       std::ostringstream error_stream;
       error_stream << "The elements within the two interacting meshes have a\n"
                    << "dimension not equal to 1 or 2.\n"
                    << "The multi-domain method will not work in this case.\n"
                    << "The dimension is: " << Dim << "\n";
       throw OomphLibError(
         error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
  
     // Call the function that actually does the work
     aux_setup_multi_domain_interaction<EXT_ELEMENT, FACE_ELEMENT_GEOM_OBJECT>(
       problem_pt,
       mesh_pt,
       external_mesh_pt,
       interaction_index,
       external_face_mesh_pt);
   }

References Global_Variables::Dim, get_dim_helper(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and Use_bulk_element_as_external.

◆ setup_multi_domain_interactions()

template<class ELEMENT_0 , class ELEMENT_1 >

void oomph::Multi_domain_functions::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.

   {
     // Delete all the current external halo(ed) element and node storage
     first_mesh_pt->delete_all_external_storage();
     second_mesh_pt->delete_all_external_storage();
  
     // Call setup_multi_domain_interaction in both directions
     setup_multi_domain_interaction<ELEMENT_1>(
       problem_pt, first_mesh_pt, second_mesh_pt, first_interaction);
  
     setup_multi_domain_interaction<ELEMENT_0>(
       problem_pt, second_mesh_pt, first_mesh_pt, second_interaction);
   }

References oomph::Mesh::delete_all_external_storage().

Referenced by RefineableConvectionProblem< NST_ELEMENT, AD_ELEMENT >::actions_after_adapt(), RefineableDDConvectionProblem< NST_ELEMENT, AD_ELEMENT >::actions_after_adapt(), ConvectionProblem< NST_ELEMENT, AD_ELEMENT >::actions_after_distribute(), RefineableConvectionProblem< NST_ELEMENT, AD_ELEMENT >::actions_after_distribute(), and DDConvectionProblem< NST_ELEMENT, AD_ELEMENT >::DDConvectionProblem().

Variable Documentation

◆ Accept_failed_locate_zeta_in_setup_multi_domain_interaction

bool oomph::Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_domain_interaction = false

Boolean to indicate that failure in setup multi domain functions is acceptable; defaults to false. If set to true external element pointers are set to null for those elements for which external elements couldn't be located.

Referenced by aux_setup_multi_domain_interaction(), MortaringHelpers::setup_constraint_elements_at_nodes(), and MortaringValidationProblem< ELEMENT, NON_MORTAR_ELEMENT >::setup_multi_domain_interaction().

◆ Allow_use_of_halo_elements_as_external_elements

bool oomph::Multi_domain_functions::Allow_use_of_halo_elements_as_external_elements = true

Boolean to indicate if we're allowed to use halo elements as external elements. Can drastically reduce the number of external halo elements – currently not aware of any problems therefore set to true by default but retention of this flag allows easy return to previous implementation.

Referenced by locate_zeta_for_local_coordinates().

◆ Allow_use_of_halo_elements_as_external_elements_for_projection

bool oomph::Multi_domain_functions::Allow_use_of_halo_elements_as_external_elements_for_projection = true

Indicate whether we are allowed to use halo elements as external elements for projection, possibly only required in parallel unstructured mesh generation during the projection stage. Default set to true

◆ Counter_for_flat_packed_doubles

unsigned oomph::Multi_domain_functions::Counter_for_flat_packed_doubles

Counter used when processing vector of flat-packed doubles – this is really "private" data, declared here to avoid having to pass it (and the associated array) between the various helper functions

◆ Counter_for_flat_packed_unsigneds

unsigned oomph::Multi_domain_functions::Counter_for_flat_packed_unsigneds

Counter used when processing vector of flat-packed unsigneds – this is really "private" data, declared here to avoid having to pass it (and the associated array) between the various helper functions

◆ Dim

unsigned oomph::Multi_domain_functions::Dim

Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in flat-packed form.

Referenced by first_closer_than_second(), get_dim_helper(), and locate_zeta_for_local_coordinates().

◆ Doc_boundary_coordinate_file

std::ofstream oomph::Multi_domain_functions::Doc_boundary_coordinate_file

Output file to document the boundary coordinate along the mesh boundary of the bulk mesh during call to setup_bulk_elements_adjacent_to_face_mesh(...)

Referenced by PseudoElasticCollapsibleChannelProblem< FLUID_ELEMENT, SOLID_ELEMENT >::PseudoElasticCollapsibleChannelProblem(), setup_bulk_elements_adjacent_to_face_mesh(), and UnstructuredFSIProblem< FLUID_ELEMENT, SOLID_ELEMENT >::UnstructuredFSIProblem().

◆ Doc_full_stats

bool oomph::Multi_domain_functions::Doc_full_stats = false

Boolean to indicate whether to output further info during setup_multi_domain_interaction() routines

Boolean to indicate whether to document further info (to screen) during setup_multi_domain_interaction() routines

◆ Doc_stats

bool oomph::Multi_domain_functions::Doc_stats = false

Boolean to indicate whether to output basic info during setup_multi_domain_interaction() routines

Boolean to indicate whether to document basic info (to screen) during setup_multi_domain_interaction() routines

Referenced by RefineableConvectionProblem< NST_ELEMENT, AD_ELEMENT >::RefineableConvectionProblem().

◆ Doc_timings

bool oomph::Multi_domain_functions::Doc_timings = false

Boolean to indicate whether to doc timings or not.

◆ External_element_located

Vector< Vector< unsigned > > oomph::Multi_domain_functions::External_element_located

Lookup scheme for whether a local element's integration point has had an external element assigned to it – essentially boolean. External_element_located[e][ipt] = {0,1} if external element for ipt-th integration in local element e {has not, has} been found. Used locally to ensure that we're not searching for the same elements over and over again when we go around the spirals.

Referenced by aux_setup_multi_domain_interaction(), clean_up(), and locate_zeta_for_local_coordinates().

◆ Flat_packed_doubles

Vector< double > oomph::Multi_domain_functions::Flat_packed_doubles

Vector of flat-packed doubles to be communicated with other processors

Referenced by clean_up().

◆ Flat_packed_located_coordinates

Vector< double > oomph::Multi_domain_functions::Flat_packed_located_coordinates

Vector of flat-packed local coordinates for zeta tuples that have been located

Referenced by clean_up().

◆ Flat_packed_unsigneds

Vector< unsigned > oomph::Multi_domain_functions::Flat_packed_unsigneds

Vector of flat-packed unsigneds to be communicated with other processors – this is really "private" data, declared here to avoid having to pass the array between the various helper functions

Referenced by clean_up().

◆ Flat_packed_zetas_not_found_locally

Vector< double > oomph::Multi_domain_functions::Flat_packed_zetas_not_found_locally

Vector of flat-packed zeta coordinates for which the external element could not be found during current local search. These will be sent to the next processor in the ring-like parallel search. The zeta coordinates come in groups of Dim (scalar) coordinates.

Referenced by aux_setup_multi_domain_interaction(), clean_up(), and locate_zeta_for_local_coordinates().

◆ Located_element_status

Vector< unsigned > oomph::Multi_domain_functions::Located_element_status

Vector to indicate (to another processor) whether a located element (that will have to represented as an external halo element on that processor) should be newly created on that processor (2), already exists on that processor (1), or is not on the current processor either (0).

Referenced by clean_up().

◆ Proc_id_plus_one_of_external_element

Vector< int > oomph::Multi_domain_functions::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 which the i-th zeta coordinate tuple received from elsewhere (in the order in which these are stored in Received_flat_packed_zetas_to_be_found) was located; it's zero if it wasn't found during the current stage of the ring-like parallel search.

Referenced by clean_up().

◆ Received_flat_packed_zetas_to_be_found

Vector< double > oomph::Multi_domain_functions::Received_flat_packed_zetas_to_be_found

Vector of flat-packed zeta coordinates for which the external element could not be found on another processor and for which we're currently searching here. Whatever can't be found here, gets written into Flat_packed_zetas_not_found_locally and then passed on to the next processor during the ring-like parallel search. The zeta coordinates come in groups of Dim (scalar) coordinates.

Referenced by clean_up().

◆ Use_bulk_element_as_external

bool oomph::Multi_domain_functions::Use_bulk_element_as_external = false

Boolean to indicate when to use the bulk element as the external element. Defaults to false, you must have set up FaceElements properly first in order for it to work

Referenced by aux_setup_multi_domain_interaction(), locate_zeta_for_local_coordinates(), and setup_multi_domain_interaction().

◆ Zeta_coords_for_further_away_comparison

Vector<double> oomph::Multi_domain_functions::Zeta_coords_for_further_away_comparison

Vector of zeta coordinates that we're currently trying to locate; used in sorting of bin entries in further_away() comparison function

Referenced by first_closer_than_second().

Enumerations

Functions

Variables

Enumeration Type Documentation

◆ anonymous enum

Function Documentation

◆ aux_setup_multi_domain_interaction() [1/2]

◆ aux_setup_multi_domain_interaction() [2/2]

◆ clean_up()

◆ first_closer_than_second()

◆ get_dim_helper()

◆ locate_zeta_for_local_coordinates()

◆ setup_bulk_elements_adjacent_to_face_mesh() [1/2]

◆ setup_bulk_elements_adjacent_to_face_mesh() [2/2]

◆ setup_multi_domain_interaction() [1/3]

◆ setup_multi_domain_interaction() [2/3]

◆ setup_multi_domain_interaction() [3/3]

◆ setup_multi_domain_interactions()

Variable Documentation

◆ Accept_failed_locate_zeta_in_setup_multi_domain_interaction

◆ Allow_use_of_halo_elements_as_external_elements

◆ Allow_use_of_halo_elements_as_external_elements_for_projection

◆ Counter_for_flat_packed_doubles

◆ Counter_for_flat_packed_unsigneds

◆ Dim

◆ Doc_boundary_coordinate_file

◆ Doc_full_stats

◆ Doc_stats

◆ Doc_timings

◆ External_element_located

◆ Flat_packed_doubles

◆ Flat_packed_located_coordinates

◆ Flat_packed_unsigneds

◆ Flat_packed_zetas_not_found_locally

◆ Located_element_status

◆ Proc_id_plus_one_of_external_element

◆ Received_flat_packed_zetas_to_be_found

◆ Use_bulk_element_as_external

◆ Zeta_coords_for_further_away_comparison