Partitionneur_Tranche.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 <Reordonner_faces_periodiques.h>
#include <Partitionneur_Tranche.h>
#include <TRUSTArray.h>
#include <Domaine.h>
#include <Param.h>
 
Implemente_instanciable_32_64(Partitionneur_Tranche_32_64, "Partitionneur_Tranche", Partitionneur_base_32_64<_T_>);
// XD partitionneur_tranche partitionneur_deriv tranche INHERITS_BRACE This algorithm will create a geometrical
// XD_CONT partitionning by slicing the mesh in the two or three axis directions, based on the geometric center of each
// XD_CONT mesh element. nz must be given if dimension=3. Each slice contains the same number of elements (slices don\'t
// XD_CONT have the same geometrical width, and for VDF meshes, slice boundaries are generally not flat except if the
// XD_CONT number of mesh elements in each direction is an exact multiple of the number of slices). First, nx slices in
// XD_CONT the X direction are created, then each slice is split in ny slices in the Y direction, and finally, each part
// XD_CONT is split in nz slices in the Z direction. The resulting number of parts is nx*ny*nz. If one particular
// XD_CONT direction has been declared periodic, the default slicing (0, 1, 2, ..., n-1)is replaced by (0, 1, 2, ...
// XD_CONT n-1, 0), each of the two \'0\' slices having twice less elements than the other slices.
// XD attr tranches listentierf tranches OPT Partitioned by nx in the X direction, ny in the Y direction, nz in the Z
// XD_CONT direction. Works only for structured meshes. No warranty for unstructured meshes.
 
 
template <typename _SIZE_>
Sortie& Partitionneur_Tranche_32_64<_SIZE_>::printOn(Sortie& os) const
{
  Cerr << "Partitionneur_Tranche_32_64<_SIZE_>::printOn invalid\n" << finl;
  Process::exit();
  return os;
}
 
/*! @brief The syntax is { Tranches nx ny [ nz ] }
 */
template <typename _SIZE_>
void Partitionneur_Tranche_32_64<_SIZE_>::set_param(Param& param) const
{
  if (!ref_domaine_)
    {
      Cerr << " Error: the domain has not been associated" << finl;
      Process::exit();
    }
  param.ajouter_arr_size_predefinie("tranches",&nb_tranches_,Param::REQUIRED);
 
}
 
/*! @brief The syntax is { Tranches nx ny [ nz ] }
 */
template <typename _SIZE_>
void Partitionneur_Tranche_32_64<_SIZE_>::validate_params() const
{
 
  // used to be done in readOn
  if (min_array(nb_tranches_)<1)
    {
      Cerr << "Error for the cutting domain tool (Tranche) specifications : " <<finl;
      Cerr<< " the number of slice must be greater than 0 for each direction. " << finl;
      Process::exit();
    }
}
 
/*! @brief First initialization step of the partitioner: associates a domain.
 *
 */
template <typename _SIZE_>
void Partitionneur_Tranche_32_64<_SIZE_>::associer_domaine(const Domaine_t& domaine)
{
  ref_domaine_ = domaine;
  nb_tranches_.resize(domaine.dimension);
  nb_tranches_ = -1;
}
 
/*! @brief Second initialization step: defines the number of slabs.
 *
 * (readOn can be used instead).
 *
 */
template <typename _SIZE_>
void Partitionneur_Tranche_32_64<_SIZE_>::initialiser(const ArrOfInt& nb_tranches)
{
 
  assert(ref_domaine_);
  assert(nb_tranches.size_array() == nb_tranches_.size_array());
  assert(min_array(nb_tranches) > 0);
 
  nb_tranches_ = nb_tranches;
}
 
/*! @brief Fills the directions_perio array from the names of the periodic boundaries.
 *
 * For 0 <= i < Objet_U::dimension, directions_perio[i] is 1 if
 *   there exists a periodic boundary whose delta vector is directed in direction i.
 *
 */
template <typename _SIZE_>
void Partitionneur_Tranche_32_64<_SIZE_>::chercher_direction_perio(const Domaine_t& domaine,
                                                                   ArrOfInt& directions_perio)
{
  Cerr << "Search of periodic directions of domain " << domaine.le_nom() << finl;
  const int dim = Objet_U::dimension;
  const Noms& liste_bords_perio = domaine.bords_perio();
 
  directions_perio.resize_array(dim);
  directions_perio = 0;
  for (auto& itr : liste_bords_perio)
    {
      int num_bord = 0;
      const Nom& nom_bord = itr;
      const int nb_bords = domaine.nb_bords();
      while (num_bord < nb_bords)
        {
          if (domaine.bord(num_bord).le_nom() == nom_bord)
            break;
          num_bord++;
        }
      if (num_bord == nb_bords)
        {
          Cerr << "Error in Partitionneur_Tranche_32_64<_SIZE_>::chercher_direction_perio\n"
               << " boundary not found : " << nom_bord << finl;
          Process::exit();
        }
      ArrOfDouble delta;
      ArrOfDouble erreur;
      const double epsilon = Objet_U::precision_geom;
      Reordonner_faces_periodiques_32_64<_SIZE_>::check_faces_periodiques(domaine.bord(num_bord), delta, erreur);
      int count = 0;
      for (int j = 0; j < dim; j++)
        {
          if (std::abs(delta[j]) > epsilon)
            {
              directions_perio[j]++;
              Cerr << " Boundary : " << nom_bord << " periodic direction : " << j << finl;
              count++;
            }
        }
      if (count != 1)
        {
          Cerr << "Error in Partitionneur_Tranche_32_64<_SIZE_>::chercher_direction_perio" << finl;
          Cerr << " periodic direction not found for the boundary " << nom_bord << finl;
          Cerr << " Vector delta found between faces twin : " << delta << finl;
          Cerr << " with a maximum error : " << erreur << finl;
          Cerr << "TIP: Try to add 'declare_only' flag into declare_bord_perio block or switch to metis tool" << finl;
          Process::exit();
        }
    }
  if (max_array(directions_perio) > 1)
    {
      Cerr << "Error in Partitionneur_Tranche_32_64<_SIZE_>::chercher_direction_perio" << finl;
      Cerr << " several boundaries have the same periodic direction" << finl;
      Process::exit();
    }
}
 
template <typename _SIZE_>
void Partitionneur_Tranche_32_64<_SIZE_>::construire_partition(BigIntVect_& elem_part, int& nb_parts_tot) const
{
  using DoubleTab_t = DoubleTab_T<_SIZE_>;
 
  assert(ref_domaine_);
  assert(nb_tranches_[0] > 0);
 
  const Domaine_t& dom = ref_domaine_.valeur();
  const int_t nb_elem = dom.nb_elem();
 
  const int dim = Objet_U::dimension;
 
  // For each space dimension, is this direction a periodic direction?
  ArrOfInt directions_perio;
  this->chercher_direction_perio(dom, directions_perio);
 
  // Center of gravity of the elements:
  Cerr << "Calculation of centers of gravity of the elements" << finl;
  DoubleTab_t coord_g;
  dom.calculer_centres_gravite(coord_g);
  assert(coord_g.dimension(0) == nb_elem);
  assert(coord_g.dimension(1) == dim);
 
  Cerr << "Moving of centers of gravity" << finl;
  // Shift coordinates so there is a strict ordering
  // in each direction (in VDF, avoid identical coordinates)
  {
    for (int_t i = 0; i < nb_elem; i++)
      {
        double x = coord_g(i,0);
        double y = coord_g(i,1);
        double z = (dim == 3) ? coord_g(i,2) : 0.;
        coord_g(i,0) = x + y*1e-6 + z*1e-12;
        coord_g(i,1) = y + x*1e-6 + z*1e-12;
        if (dim == 3)
          coord_g(i,2) = z + x*1e-6 + y*1e-12;
      }
  }
 
  // Number of elements in each part
  //  (initialization: one part containing nb_elem elements)
  SmallArrOfTID_T<_SIZE_> nb_elem_part(1);
  nb_elem_part= nb_elem;
 
  // In ascending order of parts, list of elements
  // contained in the part
  // (example:
  //  from 0..nb_elem_part[0]-1                (elements of part 0)
  //  from nb_elem_part[0]..nb_elem_part[0]+nb_elem_part[1]-1 (part 1)
  // Initialization: one large part, all elements
  ArrOfInt_t listes_elem(nb_elem);
  for (int_t i = 0; i < nb_elem; i++)
    listes_elem[i] = i;
 
  SmallArrOfTID_T<_SIZE_> new_nb_elem_part;
 
  // Index used for sorting elements
  ArrOfInt_t index;
 
  // Algorithm: create nb_tranches[0] slabs in the x direction,
  //  each slab having the same number of domain elements.
  //  Then take each slab and split it into nb_tranches[1]
  //  parts in the y direction, also balanced.
  //  Finally, take each resulting part and cut it into nb_tranches[2]
  //  parts in the z direction.
  // Note: if the domain is not a parallelepiped, the partitioning
  //  is not necessarily Cartesian:
  //    ------------------
  //    |       |        |
  //    |-------|        |
  //    |       |--------| <- slab corners do not necessarily coincide
  //    ------------------
 
  for (int direction = 0; direction < dim; direction++)
    {
      // Number of slabs to create in this direction
      const int nb_tranches = nb_tranches_[direction];
 
      int_t index_debut_partie = 0;
      new_nb_elem_part.resize_array(0);
      const int nb_parts = nb_elem_part.size_array();
 
      Cerr << "The mesh contains " << nb_parts << " parts. ";
      Cerr << "Splitting in the direction " << direction << finl;
 
      // Loop over already created parts. Divide each part
      // into nb_tranches sub-parts:
      for (int part = 0; part < nb_parts; part++)
        {
          Cerr << " Splitting of subpart " << part << finl;
          // Extract the element list for this part:
          //  point the index array to a sub-portion of listes_elem
          const int_t nb_elem_partie = nb_elem_part[part];
          assert(index_debut_partie >= 0 && index_debut_partie + nb_elem_partie <= nb_elem);
          index.ref_data(listes_elem.addr()+index_debut_partie, nb_elem_partie);
          // Sort elements of this part in ascending order
          // of the coordinate (index points to a portion of liste_elem, so
          // actually a portion of the liste_elem array is sorted)
          std::sort(index.begin(), index.end(), [&](int_t a, int_t b)
          {
            return ( coord_g(a,direction)<coord_g(b,direction) );
          });
 
          // If this direction is periodic, perform a circular shift by
          //  nb_elem_partie/(nb_tranches*2) elements so that elements at one end
          //  end up in the same part as those at the opposite end
          if (directions_perio[direction])
            {
              ArrOfInt_t copie(nb_elem_partie);
              copie.inject_array(index);
              int_t j = nb_elem_partie/(nb_tranches*2);
              for (int i = 0; i < nb_elem_partie; i++)
                {
                  index[j++] = copie[i];
                  if (j >= nb_elem_partie)
                    j = 0;
                }
            }
          // Build a partition of nb_tranches parts in
          // ascending coordinate order
          {
            int_t n = 0;
            for (int i = 0; i < nb_tranches; i++)
              {
                // All the cast because nb_elem_partie * (i+1) > 2^31 for big meshes
                const long long new_n0 = (long long)nb_elem_partie * (long long)(i+1) / (long long)nb_tranches;
                assert(new_n0 < std::numeric_limits<_SIZE_>::max());
                const int_t new_n = (int_t)new_n0;
                new_nb_elem_part.append_array(new_n - n);
                n = new_n;
              }
          }
          index_debut_partie += nb_elem_partie;
        }
 
      nb_elem_part = new_nb_elem_part;
    }
 
  // Initial partition of the domain: everything in partition 0 at the start
  elem_part.resize(nb_elem);
  elem_part = 0;
  // Build elem_part
  const int nb_parts = nb_elem_part.size_array();
  int_t index_fin = 0, i = 0;
  for (int part = 0; part < nb_parts; part++)
    {
      index_fin += nb_elem_part[part];
      for (; i < index_fin; i++)
        {
          const int_t elem = listes_elem[i];
          elem_part[elem] = part;
        }
    }
 
  if (ref_domaine_->bords_perio().size() > 0)
    this->corriger_bords_avec_liste(dom, 0, elem_part);
 
  Cerr << "Correction elem0 on processor 0" << finl;
  this->corriger_elem0_sur_proc0(elem_part);
}
 
 
template class Partitionneur_Tranche_32_64<int>;
#if INT_is_64_ == 2
template class Partitionneur_Tranche_32_64<trustIdType>;
#endif