Raccord_distant_homogene.cpp

/****************************************************************************
* Copyright (c) 2026, CEA
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*****************************************************************************/
 
#include <Raccord_distant_homogene.h>
#include <TRUSTTab.h>
#include <Schema_Comm.h>
#include <ArrOfBit.h>
#include <Domaine_dis_base.h>
#include <Domaine_VF.h>
#include <TRUSTTabs.h>
#include <Domaine.h>
#include <communications.h>
 
#include <medcoupling++.h>
#ifdef MEDCOUPLING_
#include <MEDCouplingMemArray.hxx>
 
using namespace MEDCoupling;
#endif
 
Implemente_instanciable_32_64(Raccord_distant_homogene_32_64,"Raccord_distant_homogene",Raccord_distant_32_64<_T_>);
 
/*! @brief Simple call to: Raccord_distant::printOn(Sortie& ) (+ finl)
 *
 * @param (Sortie& s) an output stream
 * @return (Sortie&) the modified output stream
 */
template <typename _SIZE_>
Sortie& Raccord_distant_homogene_32_64<_SIZE_>::printOn(Sortie& s ) const
{
  Raccord_distant_32_64<_SIZE_>::printOn(s) << finl;
  s << (int)1 <<finl;
  s << tab_envoi << finl;
  s << tab_recep << finl;
  s << nom_frontiere_voisine_ << finl;
  s << est_initialise_ << finl;
  s << e_ << finl;
  return s;
}
 
/*! @brief Simple call to: Raccord_distant::readOn(Entree& )
 *
 * @param (Entree& s) an input stream
 * @return (Entree&) the modified input stream
 */
template <typename _SIZE_>
Entree& Raccord_distant_homogene_32_64<_SIZE_>::readOn(Entree& s)
{
  Raccord_distant_32_64<_SIZE_>::readOn(s) ;
  int flag;
  s >> flag;
  if(flag==1)
    {
      s >> tab_envoi;
      s >> tab_recep;
      s >> nom_frontiere_voisine_;
      s >> est_initialise_;
      s >> e_;
    }
  return s;
}
 
 
/*! @brief Returns in DoubleTab x the trace on the distant connector of DoubleTab y located at elements of the distant domain
 */
template<>
void Raccord_distant_homogene_32_64<int>::trace_elem_distant(const DoubleTab& y, DoubleTab& x) const
{
  assert(est_initialise());
  const int nb_compo_= y.line_size();
 
  const IntTab& send_data = Tab_Envoi();
  const ArrOfInt& recv_data = Tab_Recep();
  const Faces& lesfaces = faces();
  const int n1 = send_data.dimension(0);
  const int n2 = recv_data.size_array();
 
  // Resize x if not done yet
  if (x.size_array()==0 && n2!=0)
    x.resize(n2, nb_compo_);
  else if (x.dimension(0) != n2)
    {
      Cerr << "Call to Frontiere::trace_elem with a DoubleTab x not located on boundary faces." << finl;
      Process::exit();
    }
  Schema_Comm schema;
  schema.set_send_recv_pe_list(send_pe_list_, recv_pe_list_);
  schema.begin_comm();
  for (int i = 0; i < n1; i++)
    {
      const int pe_dest = send_data(i, 0);
      const int face  = send_data(i, 1);
      int item = lesfaces.voisin(face,0);
      if (item == -1)
        item = lesfaces.voisin(face,1);
 
      for (int j=0; j<nb_compo_; j++)
        schema.send_buffer(pe_dest) << y(item,j);
    }
  schema.echange_taille_et_messages();
 
  for (int i = 0; i < n2; i++)
    {
      const int pe_source = recv_data[i];
      for(int j=0; j<nb_compo_; j++)
        schema.recv_buffer(pe_source) >> x(i,j);
    }
  schema.end_comm();
}
 
 
/*! @brief Returns in DoubleTab x the trace on the distant connector of DoubleTab y located at faces
 *
 */
template<>
void Raccord_distant_homogene_32_64<int>::trace_face_distant(const DoubleTab& y, DoubleTab& x) const
{
  assert(est_initialise());
  int n, N = y.line_size();
 
  const IntTab& send_data = Tab_Envoi();
  const ArrOfInt& recv_data = Tab_Recep();
  const int n1 = send_data.dimension(0);
  const int n2 = recv_data.size_array();
 
  // Resize x if not done yet
  if (x.size_array()==0 && n2!=0)
    {
      if (y.nb_dim() == 2) x.resize(n2, y.dimension(1));
      else if (y.nb_dim() == 3) x.resize(n2, y.dimension(1), y.dimension(2));
      else if (y.nb_dim() == 4) x.resize(n2, y.dimension(1), y.dimension(2), y.dimension(3));
    }
 
  Schema_Comm schema;
  schema.set_send_recv_pe_list(send_pe_list_, recv_pe_list_);
  schema.begin_comm();
  for (int i = 0; i < n1; i++)
    {
      const int pe_dest = send_data(i, 0);
      const int face  = num_premiere_face() + send_data(i, 1);
      for (n = 0; n < N; n++) schema.send_buffer(pe_dest) << y.addr()[N * face + n];
    }
  schema.echange_taille_et_messages();
  for (int i = 0; i < n2; i++)
    {
      const int pe_source = recv_data[i];
      for (n = 0; n < N; n++) schema.recv_buffer(pe_source) >> x.addr()[N * i + n];
    }
  schema.end_comm();
}
 
 
/*! @brief Returns in DoubleVect x the trace on the distant connector of DoubleVect y located at faces of the distant connector
 *
 */
template<>
void Raccord_distant_homogene_32_64<int>::trace_face_distant(const DoubleVect& y, DoubleVect& x) const
{
  // Verify that the array is indeed indexed on the connector faces
  assert(y.size()==nb_faces());
  assert(est_initialise());
  const IntTab& send_data = Tab_Envoi();
  const ArrOfInt& recv_data = Tab_Recep();
  const int n1 = send_data.dimension(0);
  const int n2 = recv_data.size_array();
 
  // Resize x if not done yet
  if (x.size_array()==0 && n2!=0)
    x.resize(n2);
 
  Schema_Comm schema;
  schema.set_send_recv_pe_list(send_pe_list_, recv_pe_list_);
  schema.begin_comm();
  for (int i = 0; i < n1; i++)
    {
      const int pe_dest = send_data(i, 0);
      const int face  = send_data(i, 1);
      schema.send_buffer(pe_dest) << y(face);
    }
  schema.echange_taille_et_messages();
  for (int i = 0; i < n2; i++)
    {
      const int pe_source = recv_data[i];
      schema.recv_buffer(pe_source) >> x(i);
    }
  schema.end_comm();
}
 
// 64 bits version of these should never be called:
template<typename _SIZE_>
void Raccord_distant_homogene_32_64<_SIZE_>::trace_elem_distant(const DoubleTab& y, DoubleTab& x) const
{
  assert(false);
  throw;
}
 
template<typename _SIZE_>
void Raccord_distant_homogene_32_64<_SIZE_>::trace_face_distant(const DoubleTab& y, DoubleTab& x) const
{
  assert(false);
  throw;
}
 
template<typename _SIZE_>
void Raccord_distant_homogene_32_64<_SIZE_>::trace_face_distant(const DoubleVect& y, DoubleVect& x) const
{
  assert(false);
  throw;
}
 
 
template <typename _SIZE_>
void Raccord_distant_homogene_32_64<_SIZE_>::completer()
{
  const IntTab& send_data = Tab_Envoi();
  const ArrOfInt& recv_data = Tab_Recep();
  const int n1 = send_data.dimension(0);
  const int n2 = recv_data.size_array();
 
  const int nbproc = Process::nproc();
  ArrOfBit flags(nbproc);
  flags = 0;
  for (int i = 0; i < n1; i++)
    {
      const int pe = send_data(i, 0);
      if (!flags.testsetbit(pe))
        send_pe_list_.append_array(pe);
    }
  flags = 0;
  for (int i = 0; i < n2; i++)
    {
      const int pe = recv_data[i];
      if (!flags.testsetbit(pe))
        recv_pe_list_.append_array(pe);
    }
  send_pe_list_.ordonne_array();
  recv_pe_list_.ordonne_array();
}
 
/*! @brief Initialises the distant connector with the boundary and the discretised domain opposite to the distant connector, and the discretised domain of the distant connector
 * Only called from Champ_front* instances, so can remain 32 bits only.
 */
template <>
void Raccord_distant_homogene_32_64<int>::initialise(const Frontiere& opposed_boundary,
                                                     const Domaine_dis_base& opposed_domaine_dis, const Domaine_dis_base& domaine_dis)
{
  Raccord_distant_homogene& raccord_distant = *this;
  if(raccord_distant.est_initialise())
    {
      Cerr << "Remote connection " << raccord_distant.le_nom() << " already initialized." << finl;
      Cerr << "Error, contact TRUST support." << finl;
      Process::exit();
    }
  int dim=Objet_U::dimension;
  raccord_distant.nom_frontiere_voisine()=opposed_boundary.le_nom();
  const Domaine_VF& domaine_dis_vf = ref_cast(Domaine_VF, domaine_dis);
  // Identify faces by their centres of gravity
  int parts = Process::nproc();
  DoubleTabs remote_xv(parts);
  int prem_face2 = raccord_distant.num_premiere_face();
  int nb_face2   = raccord_distant.nb_faces();
  int moi = Process::me();
  remote_xv[moi].resize(nb_face2,dim);
  for (int ind_face=0 ; ind_face < nb_face2 ; ind_face++)
    {
      int face = prem_face2 + ind_face;
      for (int j=0 ; j<dim ; j++)
        remote_xv[moi](ind_face,j) = domaine_dis_vf.xv(face,j);
    }
 
  // Then exchange the arrays of centres of gravity
  // send the arrays
  for (int p = 0; p < parts; p++)
    envoyer_broadcast(remote_xv[p], p);
 
  ArrsOfInt racc_vois(parts);
 
#ifdef MEDCOUPLING_
  // Process the information: each proc knows all XV
  // If the proc holds a portion of the raccord_distant
  int prem_face1 = opposed_boundary.num_premiere_face();
  int nb_face1   = opposed_boundary.nb_faces();
  if (nb_face1>0)
    {
      const Domaine_VF& opposed_domaine_dis_vf = ref_cast(Domaine_VF,opposed_domaine_dis);
      const DoubleTab& local_xv = opposed_domaine_dis_vf.xv();
 
      ArrOfInt& Recep=raccord_distant.Tab_Recep();
      Recep.resize_array(nb_face1);
      //double tolerance = 1e-8; //not very tolerant, are we?
      double tolerance = Objet_U::precision_geom * 100 ; //default value 1e-8 not very tolerant, are we?
 
      //DataArrayDoubles of local xv and all remote_xv (concatenated)
      std::vector<MCAuto<DataArrayDouble> > vxv(parts);
      std::vector<const DataArrayDouble*> cvxv(parts);
      for (int p = 0; p < parts; p++)
        {
          vxv[p] = DataArrayDouble::New();
          vxv[p]->useExternalArrayWithRWAccess(remote_xv[p].addr(), remote_xv[p].dimension(0), remote_xv[p].dimension(1));
          cvxv[p] = vxv[p];
        }
      MCAuto<DataArrayDouble> remote_xvs(DataArrayDouble::Aggregate(cvxv)), local_xvs(DataArrayDouble::New());
      local_xvs->alloc(nb_face1, dim);
      for (int ind_face = 0; ind_face < nb_face1; ind_face++)
        for (int j = 0; j < dim; j++)
          local_xvs->setIJ(ind_face, j, local_xv(prem_face1 + ind_face, j));
 
      //indices of the closest points in remote_xvs for each point of local_xv
      MCAuto<DataArrayIdType> glob_idx(DataArrayIdType::New());
 
      glob_idx = remote_xvs->findClosestTupleId(local_xvs);
 
      // for each face of local_xv: tolerance check, fill array
      for (int ind_face = 0, face1 = prem_face1; ind_face<nb_face1; ind_face++, face1++)
        {
          // convert global index (glob_idx(ind_face)) to the pair (proc, ind_face2)
          int proc = 0;
          mcIdType ind_face2_big = glob_idx->getIJ(ind_face, 0);
          while (ind_face2_big >= remote_xv[proc].dimension(0))
            {
              ind_face2_big -= remote_xv[proc].dimension(0);
              proc++;
            }
          assert(ind_face2_big < std::numeric_limits<int>::max());
          int ind_face2 = (int)ind_face2_big;
          assert(ind_face2 < remote_xv[proc].dimension(0));
 
          // tolerance check
          double distance2 = 0;
          for (int j=0; j<dim; j++)
            {
              double x1=local_xv(face1,j);
              double x2=remote_xv[proc](ind_face2,j);
              distance2 += (x1-x2)*(x1-x2);
            }
          if (distance2 > tolerance * tolerance)
            {
              Cerr << "Warning, there is no remote face found on the local boundary " << opposed_boundary.le_nom() << " for the face local number " << ind_face << " (";
              for (int j=0; j<dim; j++) Cerr << " " << local_xv(face1,j);
              Cerr << " )" << finl << "the nearest face on the remote boundary " << raccord_distant.le_nom() << " is the face " << ind_face2 << " (";
              for (int j=0; j<dim; j++) Cerr << " " << remote_xv[proc](ind_face2,j);
              Cerr << " )" << finl << "which is located at a distance of " << sqrt(distance2) << " and it is above the geometric tolerance " << tolerance << finl;
              Cerr << "Check your mesh or contact TRUST support." << finl;
              Process::exit();
            }
 
          // fill arrays
          Recep[ind_face]=proc;
          racc_vois[proc].append_array(ind_face2);
        }
    }
#else
  Cerr<<"Raccord_distant_homogene_32_64 needs TRUST compiled with MEDCoupling."<<finl;
  exit();
#endif
  ArrsOfInt facteurs(Process::nproc());
  envoyer_all_to_all(racc_vois, facteurs);
 
  IntTab& Envoi = raccord_distant.Tab_Envoi();
  if (nb_face2 > 0)
    Envoi.resize(nb_face2, 2);
  int ind = 0;
  for (int p = 0; p < parts; p++)
    {
      const ArrOfInt& facteur = facteurs[p];
      const int sz = facteur.size_array();
      for (int d=0 ; d<sz ; d++)
        {
          Envoi(ind,0) = p;
          Envoi(ind,1) = facteur[d];
          ind++;
        }
    }
  raccord_distant.e() = 0;
  raccord_distant.est_initialise() = 1;
  raccord_distant.completer();
  Cerr <<"Initialize the remote connection " << domaine_dis.domaine().le_nom() << "/" << raccord_distant.le_nom() << "<-" << opposed_domaine_dis.domaine().le_nom() << "/" << opposed_boundary.le_nom() << finl;
}
 
// 64 bit version should fail
template <typename _SIZE_>
void Raccord_distant_homogene_32_64<_SIZE_>::initialise(const Frontiere_t& opposed_boundary,
                                                        const Domaine_dis_base& opposed_domaine_dis, const Domaine_dis_base& domaine_dis)
{
  assert(false);
  throw;
}
 
 
template class Raccord_distant_homogene_32_64<int>;
#if INT_is_64_ == 2
template class Raccord_distant_homogene_32_64<trustIdType>;
#endif