TroisDto2D.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 <TroisDto2D.h>
#include <Interprete_bloc.h>
#include <Scatter.h>
 
Implemente_instanciable_32_64(TroisDto2D_32_64, "Extract_2D_from_3D", Interprete_geometrique_base_32_64<_T_>);
// XD extract_2d_from_3d interprete extract_2d_from_3d INHERITS_BRACE Keyword to extract a 2D mesh by selecting a
// XD_CONT boundary of the 3D mesh. To generate a 2D axisymmetric mesh prefer Extract_2Daxi_from_3D keyword.
// XD attr dom3D ref_domaine dom3D REQ Domain name of the 3D mesh
// XD attr bord chaine bord REQ Boundary name. This boundary becomes the new 2D mesh and all the boundaries, in 3D,
// XD_CONT attached to the selected boundary, give their name to the new boundaries, in 2D.
// XD attr dom2D chaine dom2D REQ Domain name of the new 2D mesh
 
 
template <typename _SIZE_>
Sortie& TroisDto2D_32_64<_SIZE_>::printOn(Sortie& os) const
{
  return Interprete::printOn(os);
}
 
 
template <typename _SIZE_>
Entree& TroisDto2D_32_64<_SIZE_>::readOn(Entree& is)
{
  return Interprete::readOn(is);
}
 
/*! @brief Fonction principale de l'interprete Extract_2D_from_3D Structure du jeu de donnee (en dimension 3) :
 *
 *     Extract_2D_from_3D dom3D bord3D dom2D
 *
 * @param (Entree& is) un flot d'entree
 * @return (Entree&) le flot d'entree
 * @throws l'objet a mailler n'est pas du type Domaine
 */
template <typename _SIZE_>
Entree& TroisDto2D_32_64<_SIZE_>::interpreter_(Entree& is)
{
  if(Objet_U::dimension==2)
    {
      Cerr << finl;
      Cerr << "You can not use " << this->que_suis_je() << " on a case already 2D" << finl;
      Process::exit();
    }
  Objet_U::dimension=3; // waooo!! violent :-)
  int test_axi=0;
  if (this->que_suis_je()=="Extract_2Daxi_from_3D")
    {
      test_axi=1;
      Cerr<<"We generate a 2D axi domain"<<finl;
    }
  Nom nom_bord, nom_dom2D;
  this->associer_domaine(is);
  is >> nom_bord >> nom_dom2D;
  Cerr << "Extraction of boundary " << nom_bord
       << " of domain " << this->domaine().le_nom()
       << " --> " << nom_dom2D << finl;
  {
    // Si l'objet nom_dom2D n'existe pas encore, on
    // ajoute un nouvel objet a l'interprete bloc courant
    Interprete_bloc& interp = Interprete_bloc::interprete_courant();
    if (interp.objet_global_existant(nom_dom2D))
      {
        if (sub_type(Domaine_t, interp.objet_global(nom_dom2D)))
          {
            Cerr << "Domain " << nom_dom2D
                 << " already exists, writing to this object." << finl;
          }
        else
          {
            Cerr << "Error : object " << nom_dom2D
                 << " already exists and is not a domain." << finl;
            Process::exit();
          }
      }
    else
      {
        Cerr << "Declaring a new domain with name " << nom_dom2D << finl;
        DerObjU ob;
        ob.typer("Domaine");
        interp.ajouter(nom_dom2D, ob);
      }
  }
  Domaine_t& dom2D = ref_cast(Domaine_t, this->objet(nom_dom2D));
  const Domaine_t& dom3D = this->domaine();
  const Bord_t& bord3D=dom3D.bord(nom_bord);
  Scatter::uninit_sequential_domain(dom2D);
  extraire_2D(dom3D, dom2D, bord3D,nom_bord, test_axi);
  Scatter::init_sequential_domain(dom2D);
  return is;
 
}
 
template <typename _SIZE_>
void TroisDto2D_32_64<_SIZE_>::extraire_2D(const Domaine_t& dom3D, Domaine_t& dom2D, const Bord_t& bord3D, const Nom& nom_bord, int test_axi)
{
  const DoubleTab_t& coord_sommets3D=dom3D.coord_sommets();
  DoubleTab_t& coord_sommets2D=dom2D.les_sommets();
  const Faces_t& faces3D=bord3D.faces();
  const IntTab_t& les_faces3D=faces3D.les_sommets();
  int_t nb_som3D=coord_sommets3D.dimension(0);
  Cerr << "number of 3D nodes :" << nb_som3D << finl;
  int_t nb_faces3D=les_faces3D.dimension(0);
  Cerr << "number of 3D faces of " << nom_bord << " :" << nb_faces3D << finl;
  int nb_som_fac3D=les_faces3D.dimension_int(1);
  Cerr << "number of nodes per faces in 3D :" << nb_som_fac3D << finl;
 
  // Construction de renum_som3D2D
  ArrOfInt_t renum_som3D2D(nb_som3D);
  for (int_t i=0; i<nb_som3D ; i++)
    renum_som3D2D[i]=-1;
  int_t compteur=0;
  for (int i=0; i<nb_faces3D; i++)
    for (int j=0; j<nb_som_fac3D; j++)
      {
        int_t& tmp=renum_som3D2D[les_faces3D(i,j)];
        if(tmp==-1)
          tmp=compteur++;
      }
 
  int_t nb_som2D=compteur;
  // Construction de renum_som2D3D
  ArrOfInt_t renum_som2D3D(nb_som2D);
 
  for (int_t i=0; i<nb_som3D ; i++)
    {
      int_t j;
      if((j=renum_som3D2D[i])!=-1)
        renum_som2D3D[j]=i;
    }
  coord_sommets2D.resize(nb_som2D,2);
  if (coupe_==2) coord_sommets2D.resize(nb_som2D,3);
  if(dom3D.type_elem()->que_suis_je()=="Hexaedre")
    if (test_axi)
      dom2D.typer("Rectangle_2D_axi");
    else
      dom2D.typer("Rectangle");
  else if(dom3D.type_elem()->que_suis_je()=="Tetraedre")
    dom2D.typer("Triangle");
  else if(dom3D.type_elem()->que_suis_je()=="Hexaedre_VEF")
    dom2D.typer("Quadrangle");
 
  if (coupe_!=1)
    {
      // Algorithme general
      // La frontiere doit etre plane (non verifie), parallele a un axe ou pas
      // Attention: les coordonnees du domaine 2D sont dans le nouveau repere
      double xb,yb,zb;
      double xc,yc,zc;
      double ux,uy,uz;
      double norm;
      double scalI,scalJ;
 
      xa=0.;
      ya=0.;
      za=0.;
      xb=0.;
      yb=0.;
      zb=0.;
      xc=0.;
      yc=0.;
      zc=0;
      ux=0.;
      uy=0.;
      uz=0.;
      norm=0.;
      scalI=0.;
      scalJ=0.;
 
      // on determine un repere orthonorme (A;IJK) sur le bord a extraire
      xa = coord_sommets3D(les_faces3D(0,0),0);
      ya = coord_sommets3D(les_faces3D(0,0),1);
      za = coord_sommets3D(les_faces3D(0,0),2);
 
      xb = coord_sommets3D(les_faces3D(0,1),0);
      yb = coord_sommets3D(les_faces3D(0,1),1);
      zb = coord_sommets3D(les_faces3D(0,1),2);
 
      xc = coord_sommets3D(les_faces3D(0,2),0);
      yc = coord_sommets3D(les_faces3D(0,2),1);
      zc = coord_sommets3D(les_faces3D(0,2),2);
 
      // calcul de I = AB/norme(AB)
      Ix=xb-xa;
      Iy=yb-ya;
      Iz=zb-za;
      norm = sqrt(Ix*Ix+Iy*Iy+Iz*Iz);
      Ix/=norm;
      Iy/=norm;
      Iz/=norm;
 
      // calcul de K = AB^AC/norme(AB^AC)
      ux=xc-xa;
      uy=yc-ya;
      uz=zc-za;
      norm = sqrt(ux*ux+uy*uy+uz*uz);
      ux/=norm;
      uy/=norm;
      uz/=norm;
      Kx = Iy*uz-uy*Iz;
      Ky = Iz*ux-Ix*uz;
      Kz = Ix*uy-Iy*ux;
      norm = sqrt(Kx*Kx+Ky*Ky+Kz*Kz);
      Kx/=norm;
      Ky/=norm;
      Kz/=norm;
 
      //  calcul de J = K^I
      Jx = Ky*Iz-Iy*Kz;
      Jy = Kz*Ix-Iz*Kx;
      Jz = Kx*Iy-Ky*Ix;
 
      // produit scalaire OA.I et OA.J
      scalI = xa*Ix+ya*Iy+za*Iz;
      scalJ = xa*Jx+ya*Jy+za*Jz;
 
      // on calcule les coordonnees des points du bord dans le nouveau repere 2D (A;I;J)
      // pour cela changement de repere du repere 3D vers le repere (A;I;J)
      for(int_t i=0; i<nb_som2D; i++)
        {
          double x = coord_sommets3D(renum_som2D3D[i],0);
          double y = coord_sommets3D(renum_som2D3D[i],1);
          double z = coord_sommets3D(renum_som2D3D[i],2);
          if (coupe_==0)
            {
              coord_sommets2D(i,0) = -scalI + x*Ix + y*Iy + z*Iz;
              coord_sommets2D(i,1) = -scalJ + x*Jx + y*Jy + z*Jz;
            }
          else
            {
              assert(coupe_==2);
// on ne ramene pas en 2D
              coord_sommets2D(i,0)=x;
              coord_sommets2D(i,1)=y;
              coord_sommets2D(i,2)=z;
            }
          //        Cout << "coord 1 = " << coord_sommets2D(i,0) << finl;
          //        Cout << "coord 2 = " << coord_sommets2D(i,1) << finl;
          //        Cout << "coord 3 = " << -xa*Kx-ya*Ky-za*Kz+x*Kx + y*Ky + z*Kz << finl;
        }
    }
  else
    {
      // Algorithme de coupe utilise par Xprepro en VDF
      // Attention: limite a une frontiere plane parallele a un axe
      // Les coordoonnees du domaine 2D sont conservees
      int test=0;
      double precision=DMAXFLOAT;
      for (test=2; test>-1; test--)
        {
          double x1,x2,x3;
          x1=coord_sommets3D(les_faces3D(0,0),test);
          x2=coord_sommets3D(les_faces3D(0,1),test);
          x3=coord_sommets3D(les_faces3D(0,2),test);
          double tmp=std::fabs(x1-x2);
          precision=(tmp<precision && tmp>0?tmp:precision);
          tmp=std::fabs(x3-x2);
          precision=(tmp<precision && tmp>0?tmp:precision);
          if (est_egal(x1,x2)&&est_egal(x2,x3)) break;
        }
      if (test==-1)
        {
          Cerr << "Error into TroisDto2D_32_64.cpp..." << finl;
          Cerr << "May be the " << dom2D.le_nom() << " boundary domain is not plane enough" << finl;
          Cerr << "Try to use: PrecisionGeom " << 100*precision << " in your data file" << finl;
          Cerr << "Or contact TRUST support." << finl;
          Process::exit();
        }
      IntVect dir(2);
      int kk=0;
      for (int jj=0; jj<3; jj++)
        if (test!=jj)
          {
            dir[kk]=jj;
            kk++;
          }
      for(int_t i=0; i<nb_som2D; i++)
        for(int j=0; j<2; j++)
          coord_sommets2D(i,j)=coord_sommets3D(renum_som2D3D[i],dir[j]);
    }
 
  IntTab_t& les_elems2D=dom2D.les_elems();
  les_elems2D=les_faces3D;
  for (int_t i=0; i< nb_faces3D; i++)
    for (int j=0; j< nb_som_fac3D; j++)
      les_elems2D(i,j)=renum_som3D2D[les_elems2D(i,j)];
 
  Objet_U::dimension=2;
  // On recupere les bords :
  for (const auto& itr: dom3D.faces_bord())
    {
      const Frontiere_t& front = itr;
      const Faces_t& faces3D_front=front.faces();
      if ( (front.le_nom()!=nom_bord))
        {
          Bord_t tmp;
          Bord_t& bord2D=dom2D.faces_bord().add(tmp);
          bord2D.nommer(front.le_nom());
          bord2D.associer_domaine(dom2D);
          if (test_axi)
            bord2D.typer_faces("QUADRILATERE_2D_AXI");
          else
            bord2D.typer_faces("segment_2D");
          // creer les faces de bord 2D ici!
          int_t compteur2=0;
          const IntTab_t& faces_sommets=faces3D_front.les_sommets();
          int_t nb_faces=faces_sommets.dimension(0);
          IntTab_t& aretes=bord2D.les_sommets_des_faces();
          IntTab_t& faces_voisins=bord2D.faces().voisins();
 
          aretes.resize(nb_faces, 2);
          // Boucle sur les faces du domaine 2D
          bool doublons=false;
          for(int_t face=0; face < nb_faces; face++)
            {
              int ok=0;
              int_t tmpbis;
              for(int i=0; i<nb_som_fac3D; i++)
                {
                  tmpbis=renum_som3D2D[faces_sommets(face,i)];
                  if (tmpbis!=-1)
                    {
                      ok++;
                      if(ok==1)
                        aretes(compteur2,0)=tmpbis;
                      if(ok==2)
                        {
                          aretes(compteur2,1)=tmpbis;
                          compteur2++;
                        }
                      if(ok==3)
                        {
                          // Si 3 sommets de la face appartienne a la frontiere
                          // il y'a un probleme, on ajoute les 3 aretes et on cherchera
                          // les doublons qu'il faudra eliminer...
                          doublons=true;
                          aretes(compteur2,0)=aretes(compteur2-1,1);
                          aretes(compteur2,1)=tmpbis;
                          compteur2++;
                          aretes(compteur2,0)=aretes(compteur2-1,1);
                          aretes(compteur2,1)=aretes(compteur2-2,0);
                          compteur2++;
                        }
                    }
                }
            }
          if (doublons)
            {
              // Recherche des doublons (algo en n^2)
              for (int_t i=0; i<compteur2; i++)
                {
                  int_t& S10=aretes(i,0);
                  int_t& S11=aretes(i,1);
                  for (int_t j=i+1; j<compteur2; j++)
                    {
                      int_t& S20=aretes(j,0);
                      int_t& S21=aretes(j,1);
                      if ( (S10==S20 && S11==S21) || (S10==S21 && S11==S20) )
                        {
                          S10=-1;
                          S11=-1;
                          S20=-1;
                          S21=-1;
                        }
                    }
                }
              // On supprime les doublons
              for (int_t i=0; i<compteur2; i++)
                if (aretes(i,0)==-1)
                  {
                    for (int_t j=i; j<compteur2-1; j++)
                      {
                        aretes(j,0)=aretes(j+1,0);
                        aretes(j,1)=aretes(j+1,1);
                      }
                    compteur2--;
                  }
            }
          aretes.resize(compteur2, 2);
          faces_voisins.resize(compteur2,2);
          faces_voisins=-1;
          if (compteur2==0)
            {
              Cerr<<"boundary to eliminate " << front.le_nom()<<finl;
              dom2D.faces_bord().suppr(bord2D);
            }
          else
            {
              Cerr<<"boundary maintained " << front.le_nom()<<finl;
            }
        }
    }
 
  // On recupere les raccords :
  for (const auto& itr: dom3D.faces_raccord())
    {
      const Frontiere_t& front = itr;
      const Faces_t& faces3Dfront=front.faces();
      if ( (front.le_nom()!=nom_bord))
        {
          Raccord_t tmpbis;
          Raccord_t& bord2D=dom2D.faces_raccord().add(tmpbis);
          bord2D.typer(itr->le_type());
          bord2D->nommer(front.le_nom());
          bord2D->associer_domaine(dom2D);
          if (test_axi)
            bord2D->typer_faces("QUADRILATERE_2D_AXI");
          else
            bord2D->typer_faces("segment_2D");
          // creer les faces de bord 2D ici!
          int_t compteur2=0;
          const IntTab_t& faces_sommets=faces3Dfront.les_sommets();
          int_t nb_faces=faces_sommets.dimension(0);
          IntTab_t& aretes=bord2D->les_sommets_des_faces();
          IntTab_t& faces_voisins=bord2D->faces().voisins();
 
          aretes.resize(nb_faces, 2);
          for(int_t face=0; face < nb_faces; face++)
            {
              int ok=0;
              int_t tmp;
              for(int i=0; i<nb_som_fac3D; i++)
                {
                  tmp=renum_som3D2D[faces_sommets(face,i)];
                  if (tmp!=-1)
                    {
                      ok++;
                      if(ok==1)
                        aretes(compteur2,0)=tmp;
                      else if(ok==2)
                        {
                          aretes(compteur2,1)=tmp;
                          compteur2++;
                          i=nb_som_fac3D;
                        }
                    }
                }
            }
          aretes.resize(compteur2, 2);
          faces_voisins.resize(compteur2,2);
          faces_voisins=-1;
          if (compteur2==0)
            {
              Cerr<<"connector to eliminate " << front.le_nom()<<finl;
              dom2D.faces_raccord().suppr(bord2D);
            }
          else
            Cerr<<"connector maintained " << front.le_nom()<<finl;
        }
    }
 
  if(dom3D.type_elem()->que_suis_je()!="Hexaedre_VEF")
    if (coupe_!=2)
      dom2D.reordonner();
  Objet_U::dimension=3;
}
 
 
template class TroisDto2D_32_64<int>;
#if INT_is_64_ == 2
template class TroisDto2D_32_64<trustIdType>;
#endif