Declarer_bord_perio.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 <Declarer_bord_perio.h>
#include <Format_Post_base.h>
#include <Synonyme_info.h>
#include <Octree_Double.h>
#include <Domaine_bord.h>
#include <ArrOfBit.h>
#include <Param.h>
 
Implemente_instanciable_32_64(Declarer_bord_perio_32_64,"Declarer_bord_perio",Interprete_geometrique_base_32_64<_T_>);
Add_synonym(Declarer_bord_perio,"Corriger_frontiere_periodique");
Add_synonym(Declarer_bord_perio_64,"Corriger_frontiere_periodique_64");
 
// XD declarer_bord_perio interprete corriger_frontiere_periodique BRACE The Declarer_bord_perio keyword is mandatory to
// XD_CONT first define the periodic boundaries, to reorder the faces and eventually fix unaligned nodes of these
// XD_CONT boundaries. Faces on one side of the periodic domain are put first, then the faces on the opposite side, in
// XD_CONT the same order. It must be run in sequential before mesh splitting.
// XD attr domaine chaine domaine REQ Name of domain.
// XD attr bord chaine bord REQ the name of the boundary (which must contain two opposite sides of the domain)
// XD attr direction list direction OPT defines the periodicity direction vector (a vector that points from one node on
// XD_CONT one side to the opposite node on the other side). This vector must be given if the automatic algorithm fails,
// XD_CONT that is:NL2 - when the node coordinates are not perfectly periodic NL2 - when the periodic direction is not
// XD_CONT aligned with the normal vector of the boundary faces
// XD attr fichier_post chaine fichier_post OPT .
 
template <typename _SIZE_>
Entree& Declarer_bord_perio_32_64<_SIZE_>::readOn(Entree& is)
{
  return is;
}
 
template <typename _SIZE_>
Sortie& Declarer_bord_perio_32_64<_SIZE_>::printOn(Sortie& os) const
{
  return os;
}
 
/*! @brief Main routine recording the periodic boundaries in the domain and performing the appropriate reordering of the faces
 */
template <typename _SIZE_>
void Declarer_bord_perio_32_64<_SIZE_>::adapt_som_and_faces()
{
  Domaine_t& dom = this->domaine();
  const int_t dim = dom.coord_sommets().dimension(1);
  int direction_perio_set;
  if (direction_perio_.size_array() == 0)
    {
      direction_perio_set=0;
      Cerr << "No direction given, searching periodicity direction automatically:" << finl;
      Reordonner_faces_periodiques_32_64<_SIZE_>::chercher_direction_perio(direction_perio_, dom, nom_bord_);
    }
  else
    {
      direction_perio_set=1;
      if (direction_perio_.size_array() != dim)
        {
          Cerr << "Error in Declarer_bord_perio::interpreter: direction should be of size "
               << dim << finl;
          this->exit();
        }
    }
  Cerr << "Searching and moving periodicity nodes for domain " << dom.le_nom() << " boundary " << nom_bord_ << finl;
  corriger_coordonnees_sommets_perio();
  Bord_t& bord = dom.bord(nom_bord_);
  IntTab_t& faces = bord.faces().les_sommets();
  const double epsilon = Objet_U::precision_geom;
  const int ok = Reordonner_faces_periodiques_32_64<_SIZE_>::reordonner_faces_periodiques(dom, faces, direction_perio_, epsilon);
  if (!ok)
    {
      if (direction_perio_set)
        Cerr << "May be you give a wrong vector for the DIRECTION option" << finl;
      else
        Cerr << "May be you could use the DIRECTION option to specify the vector of periodicity" << finl;
      Cerr << "for the keyword:" << finl;
      Cerr << "Declarer_bord_perio { Domaine " << dom.le_nom() << " Bord " << nom_bord_ << " } " << finl;
      this->exit();
    }
}
 
template <typename _SIZE_>
Entree& Declarer_bord_perio_32_64<_SIZE_>::interpreter_(Entree& is)
{
  bool declare_only = false;
 
  Nom nom_dom;
  Param param(this->que_suis_je());
  param.ajouter("domaine", &nom_dom, Param::REQUIRED);
  param.ajouter("bord", &nom_bord_, Param::REQUIRED);
  param.ajouter("direction", &direction_perio_);
  param.ajouter("fichier_post", &nom_fichier_post_);
  param.ajouter_flag("declare_only", &declare_only);
  param.lire_avec_accolades_depuis(is);
 
  if (this->nproc() > 1)
    {
      Cerr << "Error in Declarer_bord_perio::interpreter():\n"
           << " this function must be run in sequential before mesh splitting." << finl;
      this->barrier();
      this->exit();
    }
 
  this->associer_domaine(nom_dom);
 
  if(!declare_only)
    adapt_som_and_faces();
 
  // Register 'bord' in the list of periodic boundary of the domain:
  Domaine_t& dom = this->domaine();
  dom.bords_perio().add(nom_bord_);
 
  return is;
}
 
 
/*! @brief Pour chaque sommet du bord periodique, on cherche son sommet oppose dans la direction + ou - vecteur_perio et dans un rayon search_radius.
 *
 *   Parmi le couple de sommets forme, on deplace celui qui se trouve en +vecteur_perio
 *   pour le mettre en face de l'autre sommet.
 *   Exit en cas d'erreur (si les sommets sont trop eloignes de leur sommet associe)
 */
template <typename _SIZE_>
void Declarer_bord_perio_32_64<_SIZE_>::corriger_coordonnees_sommets_perio()
{
  if (Process::nproc() > 1)
    {
      Cerr << "Error in Declarer_bord_perio_32_64::corriger_coordonnees_sommets_perio\n"
           << " this algorithm is sequential (use it before Decouper)" << finl;
      Process::barrier();
      Process::exit();
    }
 
  Domaine_t& dom = this->domaine();
  Domaine_bord_t domaine_bord;
  domaine_bord.construire_domaine_bord(dom, nom_bord_);
  const ArrOfInt_t& renum_som = domaine_bord.get_renum_som();
  const DoubleTab_t& som_bord = domaine_bord.les_sommets();
  Octree_Double_t octree;
  octree.build_nodes(som_bord, 0 /* do not include virtual nodes */);
 
  const double epsilon_initial = dom.epsilon();
  if (epsilon_initial == 0.)
    {
      Cerr << "Error in Declarer_bord_perio_32_64::corriger_coordonnees_sommets_perio\n"
           << " dom.epsilon = 0." << finl;
      Process::exit();
    }
 
  int error_flag = 0;
 
  const int_t nb_som_bord = som_bord.dimension(0);
  const int dim = static_cast<int>(som_bord.dimension(1));
  // Un tableau contenant le deplacement applique aux sommets (pour postraitement)
  DoubleTab_t delta(nb_som_bord, dim);
  // Un tableau pour detecter les erreurs de periodicite: sommets associes plusieurs fois
  ArrOfBit_t marker(nb_som_bord);
  marker = 0;
  // Tableau temporaire pour l'octree:
  ArrOfInt_t nodes_list;
 
  DoubleTab_t& sommets_src = dom.les_sommets();
  ArrOfDouble coord(dim);
  for (int_t som = 0; som < nb_som_bord; som++)
    {
      if (marker.testsetbit(som))
        continue; // Sommet deja traite
 
      // Cherche le sommet oppose dans les deux directions (-1. et +1.)
      double facteur;
      int_t som2 = -1;
      // Recherche du sommet dans un rayon de plus en plus grand en commencant par epsilon:
      double epsilon = epsilon_initial;
      do
        {
          for (facteur = -1.; facteur < 1.5; facteur += 2.)
            {
              for (int i = 0; i < dim; i++)
                coord[i] = som_bord(som, i) + facteur * direction_perio_[i];
              octree.search_elements_box(coord, epsilon, nodes_list);
              som2 = octree.search_nodes_close_to(coord, som_bord, nodes_list, epsilon);
              if (som2 >= 0)
                break;
            }
          if (som2 >= 0)
            break;
          // Sommet non trouve, on recommence avec un rayon de recherche plus grand
          epsilon *= 4.;
        }
      while (1);
      if (marker.testsetbit(som2))
        {
          Cerr << "Error in Declarer_bord_perio_32_64::corriger_coordonnees_sommets_perio\n"
               << " Coordinate    " << coord
               << " Closest point [ " << som_bord(som2, 0) << " " << som_bord(som2, 0) << " "
               << ((dim==3)?som_bord(som2, 0):0.) << " ] already used for another point" << finl;
          error_flag = 1;
        }
      else
        {
          // On deplace le sommet vers lequel pointe vecteur_perio:
          const int_t som_ref = som;
          const int_t som_deplace = som2;
          // Indice du sommet a deplacer dans le domaine source:
          const int_t s = renum_som[som_deplace];
          // Deplacement
          for (int i = 0; i < dim; i++)
            {
              double old_x = som_bord(som_deplace, i);
              double new_x = som_bord(som_ref, i) + facteur * direction_perio_[i];
              sommets_src(s, i) = new_x;
              delta(som_deplace, i) = new_x - old_x;
            }
        }
    }
 
  if (nom_fichier_post_ != "??")
    {
      if (!std::is_same<_SIZE_, int>::value)
        {
          Cerr << "Declarer_bord_perio_64 - option 'fichier_post' is not implemented yet for 64b domains!" << finl;
          Cerr << "Remove the attribute or contact TRUST support." << finl;
          Process::exit(-1);
        }
#if INT_is_64_ != 2
      //Domaine dom2;
      Cerr << "Writing node displacement into file: " << nom_fichier_post_ << finl;
      OWN_PTR(Format_Post_base) fichier_post;
      fichier_post.typer("FORMAT_POST_LATA");
      fichier_post->initialize(nom_fichier_post_, 1 /* binaire */, "SIMPLE");
      Format_Post_base& post = fichier_post.valeur();
      post.ecrire_entete(0., 0 /*reprise*/, 1 /* premier post */);
      post.ecrire_domaine(domaine_bord, 1 /* premier_post */);
      post.ecrire_temps(0.);
      Noms unites;
      unites.add("m");
      unites.add("m");
      if (dim==3) unites.add("m");
      Noms compos;
      compos.add("dx");
      compos.add("dy");
      if (dim==3) compos.add("dz");
      post.ecrire_champ(domaine_bord,
                        unites, compos, -1 /* ecrire toutes les composantes */,
                        0., /* temps */
                        "vitesse", domaine_bord.le_nom(), "SOM","vector", delta);
      int fin=1;
      post.finir(fin);
#endif
    }
 
  Cerr << "Max of node displacement (m): " << max_abs_array(delta) << finl;
  if (error_flag)
    {
      Cerr << "Could not correct node coordinates. Aborting." << finl;
      Process::exit();
    }
}
 
template class Declarer_bord_perio_32_64<int>;
#if INT_is_64_ == 2
template class Declarer_bord_perio_32_64<trustIdType>;
#endif