Extraire_plan.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 <Extraire_surface.h>
#include <Probleme_base.h>
#include <Extraire_plan.h>
#include <Equation_base.h>
#include <NettoieNoeuds.h>
#include <Domaine_VF.h>
#include <EChaine.h>
#include <SChaine.h>
#include <Domaine.h>
#include <Param.h>
 
Implemente_instanciable(Extraire_plan,"Extraire_plan",Interprete_geometrique_base);
// XD extraire_plan interprete extraire_plan BRACE This keyword extracts a plane mesh named domain_name (this domain
// XD_CONT should have been declared before) from the mesh of the pb_name problem. The plane can be either a triangle
// XD_CONT (defined by the keywords Origine, Point1, Point2 and Triangle), either a regular quadrangle (with keywords
// XD_CONT Origine, Point1 and Point2), or either a generalized quadrangle (with keywords Origine, Point1, Point2,
// XD_CONT Point3). The keyword Epaisseur specifies the thickness of volume around the plane which contains the faces of
// XD_CONT the extracted mesh. The keyword via_extraire_surface will create a plan and use Extraire_surface algorithm.
// XD_CONT Inverse_condition_element keyword then will be used in the case where the plane is a boundary not well
// XD_CONT oriented, and avec_certains_bords_pour_extraire_surface is the option related to the Extraire_surface option
// XD_CONT named avec_certains_bords. Point1, Point2 and Triangle), either a regular quadrangle (with keywords Origine,
// XD_CONT Point1 and Point2), or either a generalized quadrangle (with keywords Origine, Point1, Point2, Point3). The
// XD_CONT keyword Epaisseur specifies the thickness of volume around the plane which contains the faces of the
// XD_CONT extracted mesh. The keyword via_extraire_surface will create a plan and use Extraire_surface algorithm.
// XD_CONT Inverse_condition_element keyword then will be used in the case where the plane is a boundary not well
// XD_CONT oriented, and avec_certains_bords_pour_extraire_surface is the option related to the Extraire_surface option
// XD_CONT named avec_certains_bords.
 
Sortie& Extraire_plan::printOn(Sortie& os) const { return Interprete::printOn(os); }
 
Entree& Extraire_plan::readOn(Entree& is) { return Interprete::readOn(is); }
 
void calcul_normal_norme(const ArrOfDouble& org, const ArrOfDouble& point1, const ArrOfDouble& point2, ArrOfDouble& normal)
{
  calcul_normal(org.addr(), point1.addr(),  point2.addr(), normal.addr());
  normal/=norme_array(normal);
}
 
Entree& Extraire_plan::interpreter_(Entree& is)
{
  Nom nom_pb;
  Nom nom_dom;
  ArrOfDouble origine,point1,point2,point3;
  bool triangle = false;
  double epaisseur = -123.;
  Param param(que_suis_je());
  param.ajouter("domaine",&nom_dom,Param::REQUIRED); // XD_ADD_P ref_domaine
  // XD_CONT domain name
  param.ajouter("probleme",&nom_pb,Param::REQUIRED); // XD_ADD_P ref_Pb_base
  // XD_CONT pb_name
  param.ajouter("origine",&origine,Param::REQUIRED); // XD_ADD_P list
  // XD_CONT not_set
  param.ajouter("point1",&point1,Param::REQUIRED); // XD_ADD_P list
  // XD_CONT not_set
  param.ajouter("point2",&point2,Param::REQUIRED); // XD_ADD_P list
  // XD_CONT not_set
  param.ajouter("point3",&point3); // XD_ADD_P list
  // XD_CONT not_set
  param.ajouter_flag("triangle",&triangle); // XD_ADD_P rien
  // XD_CONT not_set
  param.ajouter("epaisseur",&epaisseur,Param::REQUIRED); // XD_ADD_P floattant
  // XD_CONT thickness
 
  bool via_extraire_surface = false;
  bool inverse_condition_element = false;
  param.ajouter_flag("via_extraire_surface",&via_extraire_surface); // XD_ADD_P rien
  // XD_CONT not_set
  param.ajouter_flag("inverse_condition_element",&inverse_condition_element); // XD_ADD_P rien
  // XD_CONT not_set
 
  Noms bords;
 
  param.ajouter("avec_certains_bords_pour_extraire_surface",&bords); // XD_ADD_P listchaine
  // XD_CONT name of boundaries to include when extracting plan
 
  param.lire_avec_accolades_depuis(is);
 
  associer_domaine(nom_dom);
 
  if (point3.size_array()==0)
    {
      point3.resize_array(3);
      for (int dir=0; dir<3; dir++)
        if (triangle)
          point3[dir]=(point1[dir]+point2[dir])/2.;
        else
          point3[dir]=point1[dir]+point2[dir]-origine[dir];
    }
 
 
  if (via_extraire_surface)
    {
      // determination de la normale
      ArrOfDouble normal(3);
      calcul_normal_norme(origine,point1,point2,normal);
 
      double dx=0;
      for (int dir=0; dir<3; dir++)
        dx+=normal[dir]*origine[dir];
      Noms axes(3);
      axes[0]="x";
      axes[1]="y";
      axes[2]="z";
 
      SChaine out;
      out << " { domaine "<<nom_dom <<" probleme "<<nom_pb <<finl;
      out<<" condition_elements (" <<-dx ;
      for (int dir=0; dir<3; dir++)
        if (normal[dir])
          out<<"+"<<axes[dir]<<"*("<<normal[dir]<<")";
      if (inverse_condition_element)
        out<<")_ge_0"<<finl;
      else
        out<<")_le_0"<<finl;
      out<<" condition_faces (";
      for (int ori=0; ori<2; ori++)
 
        {
          ArrOfDouble& A=(ori==0?origine:point3);
          ArrOfDouble& D=(ori!=0?origine:point3);
          ArrOfDouble normal2(3),prov(3);
          for (int te=0; te<2; te++)
            {
              ArrOfDouble& B=(te==0?point1:point2);
              //ArrOfDouble& C=(ori!=0?point1:point2);
              prov=normal;
              prov+=A;
              calcul_normal_norme(A,B,prov,normal2);
              //prov=REF;
              //prov-=A;
              double dxbis=0;
              double ref=0;
              out<<"((";
              for (int dir=0; dir<3; dir++)
                {
                  dxbis+=normal2[dir]*D[dir];
                  ref+=normal2[dir]*A[dir];
                  if (normal2[dir])
                    out << axes[dir]<<"*("<< normal2[dir]<<")+";
                }
              out<<"("<<-ref<<"))";
              //Cout<<"INFO " << dx-ref<<finl;
              if ((dxbis-ref)>0)
                out <<"_ge_0)*";
              else
                out <<"_le_0)*";
            }
        }
      out<<"1)"<<finl;
      //out<<" avec_les_bords "<<finl;
      if (bords.size())
        out<<"avec_certains_bords " <<bords<<finl;
      out<<"}" << finl;
      Cout<<" Creation of the plan through Extraire_surface "<< out.get_str()<<finl;;
      EChaine in(out.get_str());
      Extraire_surface extraction;
      extraction.interpreter(in);
      // abort();
      return is;
    }
 
  Domaine& dom=domaine();
  // on recupere le pb
  if(! sub_type(Probleme_base, objet(nom_pb)))
    {
      Cerr << nom_pb << " is of type " << objet(nom_pb).que_suis_je() << finl;
      Cerr << "and not of type Probleme_base" << finl;
      exit();
    }
  Probleme_base& pb=ref_cast(Probleme_base, objet(nom_pb));
  const Domaine_VF& domaine_vf=ref_cast(Domaine_VF,pb.domaine_dis());
  dom.les_sommets()=domaine_vf.domaine().les_sommets();
  const DoubleTab& coord=dom.les_sommets();
  const Nom& type_elem=domaine_vf.domaine().type_elem()->que_suis_je();
  if (type_elem==Motcle("Tetraedre"))
    dom.typer("Triangle");
  else
    {
      Cerr<<que_suis_je()<<"  not coded for this type of elements "<<finl;
      exit();
    }
 
  // creation d'un domaine pipo pour pouvoir chercher les faces
  Domaine domaine_test;
  ArrOfDouble normal(3);
  {
    if (triangle)
      domaine_test.typer("Prisme");
    else
      domaine_test.typer("Hexaedre_vef");
 
    DoubleTab& somm_hexa=domaine_test.les_sommets();
    int nb_som=4;
    if (triangle) nb_som=3;
    somm_hexa.resize(nb_som*2,3);
 
    for (int dir=0; dir<3; dir++)
      {
        somm_hexa(0,dir)=origine[dir];
        somm_hexa(1,dir)=point1[dir];
        somm_hexa(2,dir)=point2[dir];
        somm_hexa(3,dir)=point3[dir];
        for (int t=0; t<nb_som; t++)
          somm_hexa(t+nb_som,dir)=somm_hexa(t,dir);
      }
    // il faut maintenant epaissir l'hexa dans la direction normale
 
    calcul_normal_norme(origine,point1,point2,normal);
    double d_epaisseur=epaisseur/2.;
    for (int t=0; t<nb_som; t++)
      for (int dir=0; dir<3; dir++)
        {
          somm_hexa(t,dir)=somm_hexa(t,dir)-d_epaisseur*normal[dir];
          somm_hexa(nb_som+t,dir)=somm_hexa(nb_som+t,dir)+d_epaisseur*normal[dir];
        }
    IntTab& elem_test=domaine_test.les_elems();
    elem_test.resize(1,nb_som*2);
    for (int s=0; s<nb_som*2; s++) elem_test(0,s)=s;
    /*
      SFichier es("hex.geom");
      es<<test;
    */
  }
  const DoubleTab& xv =domaine_vf.xv();
 
  int nbfaces=domaine_vf.nb_faces();
 
  ArrOfInt marq(nbfaces);
  // on marque les joints
  //const Joints& joints=domaine_vf.face_joints();
  int nbjoints=domaine_vf.nb_joints();
 
  for(int njoint=0; njoint<nbjoints; njoint++)
    {
      const Joint& joint_temp = domaine_vf.joint(njoint);
      int pe_voisin=joint_temp.PEvoisin();
      if (pe_voisin<me())
        {
          const IntTab& indices_faces_joint = joint_temp.joint_item(JOINT_ITEM::FACE).renum_items_communs();
          const int nb_faces_j = indices_faces_joint.dimension(0);
          for (int j = 0; j < nb_faces_j; j++)
            {
              int face_de_joint = indices_faces_joint(j, 1);
              marq[face_de_joint] = -1;
            }
        }
    }
  int nb_t=0;
 
  for (int fac=0; fac<nbfaces; fac++)
    {
      if (domaine_test.chercher_elements(xv(fac,0),xv(fac,1),xv(fac,2))==0)
        if (marq[fac]!=-1)
          {
            // tester si item_commun....
            marq[fac]=1;
            nb_t++;
          }
    }
  ArrOfDouble point0b(3),point1b(3),point2b(3);
  Cerr<<"Number of elements of the new domain "<<nb_t<<finl;
  IntTab& les_elems=dom.les_elems();
  les_elems.resize(nb_t,3);
  const IntTab& face_sommets=domaine_vf.face_sommets();
  int nb=0;
  for (int fac=0; fac<nbfaces; fac++)
    if (marq[fac]==1)
      {
        //Cerr<<fac <<" ";
        for (int s=0; s<3; s++)
          les_elems(nb,s)=face_sommets(fac,s);
        // on calcule la normale
        ArrOfDouble normal_b(3);
        for (int i=0; i<3; i++)
          {
            point0b[i]=coord(les_elems(nb,0),i);
            point1b[i]=coord(les_elems(nb,1),i);
            point2b[i]=coord(les_elems(nb,2),i);
          }
        calcul_normal(point0b.addr(),point1b.addr(),point2b.addr(),normal_b.addr());
        if (dotproduct_array(normal,normal_b)<0)
          {
            // si normal a l'envers on inverse les deux sommets
            les_elems(nb,1)=face_sommets(fac,2);
            les_elems(nb,2)=face_sommets(fac,1);
          }
        nb++;
      }
  Cerr<<finl;;
  assert(nb==nb_t);
  NettoieNoeuds::nettoie(dom);
 
  return is;
}