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 Main function of the Extract_2D_from_3D interpreter. Data set structure (in dimension 3):
 *
 *     Extract_2D_from_3D dom3D bord3D dom2D
 *
 * @param is Input stream.
 * @return The input stream.
 * @throws The object to mesh is not of 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;
  {
    // If the object nom_dom2D does not yet exist,
    // add a new object to the current block interpreter
    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;
 
  // Build 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;
  // Build 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)
    {
      // General algorithm
      // The boundary must be planar (not verified), parallel to an axis or not
      // Warning: the coordinates of the 2D domain are in the new frame
      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.;
 
      // determine an orthonormal frame (A;IJK) on the boundary to extract
      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);
 
      // compute I = AB/norm(AB)
      Ix=xb-xa;
      Iy=yb-ya;
      Iz=zb-za;
      norm = sqrt(Ix*Ix+Iy*Iy+Iz*Iz);
      Ix/=norm;
      Iy/=norm;
      Iz/=norm;
 
      // compute K = AB^AC/norm(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;
 
      //  compute J = K^I
      Jx = Ky*Iz-Iy*Kz;
      Jy = Kz*Ix-Iz*Kx;
      Jz = Kx*Iy-Ky*Ix;
 
      // dot products OA.I and OA.J
      scalI = xa*Ix+ya*Iy+za*Iz;
      scalJ = xa*Jx+ya*Jy+za*Jz;
 
      // compute the coordinates of the boundary points in the new 2D frame (A;I;J)
      // by performing a change of frame from the 3D frame to (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);
// do not reduce to 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
    {
      // Cutting algorithm used by Xprepro in VDF
      // Warning: limited to a planar boundary parallel to an axis
      // The coordinates of the 2D domain are preserved
      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;
  // Retrieve the boundaries:
  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");
          // create the 2D boundary faces here!
          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);
          // Loop over the faces of the 2D domain
          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)
                        {
                          // If 3 vertices of the face belong to the boundary
                          // there is a problem: add all 3 edges and look for
                          // the duplicates that will need to be removed...
                          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)
            {
              // Search for duplicates (O(n^2) algorithm)
              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;
                        }
                    }
                }
              // Remove duplicates
              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;
            }
        }
    }
 
  // Retrieve the connectors:
  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");
          // create the 2D boundary faces here!
          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