Aretes.cpp

/****************************************************************************
* Copyright (c) 2022, 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 <Aretes.h>
#include <Domaine_VDF.h>
 
Implemente_instanciable(Aretes,"Aretes",Objet_U);
 
 
// printOn and readOn
 
Sortie& Aretes::printOn(Sortie& s ) const
{
  s << "Aretes " << finl;
  s << faces_ << finl;
  s << type1_ << finl;
  s << type2_ << finl;
  return s;
}
 
Entree& Aretes::readOn(Entree& s )
{
  return s ;
}
 
/*! @brief Assigns faces f1, f2, f3, f4 to edge numero.
 *
 * @brief The edge type2_ is defined as:
 *  type2_ = -1 for a corner edge
 *  type2_ = 0 for a boundary edge
 *  type2_ = 1 for a mixed edge
 *  type2_ = 2 for an internal edge
 *  The edge type1_ dir is defined as:
 *  type1_ = 0 for an XY edge
 *  type1_ = 1 for an XZ edge
 *  type1_ = 2 for a YZ edge
 *  In 2D there are only XY edges.
 *
 * @param numero_a Edge index (incremented in place).
 * @param dir Edge direction type.
 * @param type Edge category type.
 * @param nb_face Total number of faces.
 * @param f1 First face index.
 * @param f2 Second face index.
 * @param f3 Third face index.
 * @param f4 Fourth face index.
 * @param est_une_plaque Array indicating which faces are plate faces.
 */
void Aretes::affecter(int& numero_a, int dir, int type, int nb_face,
                      int f1, int f2, int f3, int f4, const ArrOfInt& est_une_plaque)
{
  int nb_plaques = 0;
  nb_plaques += f1>=0 ? est_une_plaque[f1] : 0;
  nb_plaques += f2>=0 ? est_une_plaque[f2] : 0;
  nb_plaques += f3>=0 ? est_une_plaque[f3] : 0;
  nb_plaques += f4>=0 ? est_une_plaque[f4] : 0;
  int coin = -1 ;
  int bord = 0 ;
  if (type>0)
    {
      if (nb_plaques!=0) return;
      if( (f1<nb_face) || (f2<nb_face) || (f3<nb_face) || (f4<nb_face) )
        numero_a++;
      else
        return;
    }
  if (type==bord)
    {
      if( (f1<nb_face) || (f2<nb_face) || (f3<nb_face) )
        numero_a++;
      else
        return;
    }
  if (type==coin)
    {
      if( (f1<nb_face) && (f2<nb_face) && (f3<nb_face) && (f4<nb_face) )
        numero_a++;
      else
        return;
    }
 
  faces_(numero_a, 0)=f1;
  faces_(numero_a, 1)=f2;
  faces_(numero_a, 2)=f3;
  faces_(numero_a, 3)=f4;
  type1_(numero_a)=dir;
  type2_(numero_a)=type;
 
  if(type >0)
    {
      assert(faces_(numero_a, 0) !=-1);
      assert(faces_(numero_a, 1) !=-1);
      assert(faces_(numero_a, 2) !=-1);
      assert(faces_(numero_a, 3) !=-1);
    }
}
 
void Aretes::calculer_centre_de_gravite(Domaine_VDF& domaine)
{
  const IntTab& so = domaine.face_sommets();
  const DoubleTab& co = domaine.domaine().les_sommets();
  DoubleTab& xa_ = domaine.xa();
  int i,j,k;
  int f0=-1,f1=-1,s00,s01,s10,s11;
  int type;
  int nb_aretes = faces_.dimension(0);
  //const IntVect& orient =  domaine.orientation();
  // Compute the edge coordinates
  if(dimension==2)
    {
      xa_.resize(nb_aretes,2);
      for(i=0; i<nb_aretes; i++)
        {
          type = type2_(i);
          if((type==2)||(type==1)) // internal or mixed edge
            {
              f0 = faces_(i,0);
              s00 = so(f0,0);
              s01 = so(f0,1);
              f1 = faces_(i,1);
              s10 = so(f1,0);
              s11 = so(f1,1);
              if((s00==s10)||(s00==s11))
                {
                  xa_(i,0) = co(s00,0);
                  xa_(i,1) = co(s00,1);
                }
              else if((s01==s10)||(s01==s11))
                {
                  xa_(i,0) = co(s01,0);
                  xa_(i,1) = co(s01,1);
                }
              else
                {
                  Cerr<<"Error: no common vertices found"<<finl;
                  exit();
                }
            }
          else if((type == 0)||(type == -1)) // boundary or corner edge
            {
              for(j=0; j<4; j++)
                {
                  f0=faces_(i,j);
                  if(f0!=-1) break;
                }
              assert(j<4);
              k=j+1;
              s00 = so(f0,0);
              s01 = so(f0,1);
              for(j=k; j<4; j++)
                {
                  f1=faces_(i,j);
                  if(f1!=-1) break;
                }
              assert(j<4);
              s10 = so(f1,0);
              s11 = so(f1,1);
              if((s00==s10)||(s00==s11))
                {
                  xa_(i,0) = co(s00,0);
                  xa_(i,1) = co(s00,1);
                }
              else if((s01==s10)||(s01==s11))
                {
                  xa_(i,0) = co(s01,0);
                  xa_(i,1) = co(s01,1);
                }
              else
                {
                  Cerr<<"Error: no common vertices found"<<finl;
                  exit();
                }
            }
        }
    }
  else if(dimension==3)
    {
      xa_.resize(nb_aretes,3);
      for(i=0; i<nb_aretes; i++)
        {
          type = type2_(i);
          if((type==2)||(type==1)) // internal or mixed edge
            {
              f0 = faces_(i,0);
              f1 = faces_(i,1);
              int s0,s1,s0j;
              int deux;
              for(j=0; j<4; j++)
                {
                  s0j=so(f0,j);
                  for(k=0; k<4; k++)
                    if(s0j==so(f1,k)) break;
                  if(k<4) break;
                }
              assert(j<4);
              s0 = s0j;
              deux = j+1;
              for(j=deux; j<4; j++)
                {
                  s0j=so(f0,j);
                  for(k=0; k<4; k++)
                    if(s0j==so(f1,k)) break;
                  if(k<4) break;
                }
              assert(j<4);
              s1 = s0j;
              xa_(i,0) = (co(s0,0)+co(s1,0))/2.0;
              xa_(i,1) = (co(s0,1)+co(s1,1))/2.0;
              xa_(i,2) = (co(s0,2)+co(s1,2))/2.0;
            }
          else if((type == 0)||(type == -1)) // boundary or corner edge
            {
              int f,fdeux;
              int s0,s1,s0j;
              int deux;
              for(f=0; f<4; f++)
                {
                  f0=faces_(i,f);
                  if(f0!=-1) break;
                }
              assert(f<4);
              fdeux=f+1;
              for(f=fdeux; f<4; f++)
                {
                  f1=faces_(i,f);
                  if(f1!=-1) break;
                }
              for(j=0; j<4; j++)
                {
                  s0j=so(f0,j);
                  for(k=0; k<4; k++)
                    if(s0j==so(f1,k)) break;
                  if(k<4) break;
                }
              assert(j<4);
              s0 = s0j;
              deux = j+1;
              for(j=deux; j<4; j++)
                {
                  s0j=so(f0,j);
                  for(k=0; k<4; k++)
                    if(s0j==so(f1,k)) break;
                  if(k<4) break;
                }
              assert(j<4);
              s1 = s0j;
              xa_(i,0) = (co(s0,0)+co(s1,0))/2.0;
              xa_(i,1) = (co(s0,1)+co(s1,1))/2.0;
              xa_(i,2) = (co(s0,2)+co(s1,2))/2.0;
            }
        }
    }
}
 
/*! @brief Resizes the arrays.
 *
 */
void Aretes::dimensionner(int n)
{
  faces_.resize(n,4);
  type1_.resize(n);
  type2_.resize(n) ;
}
 
/*! @brief Called by trier. Swaps edges a1 and a2.
 *
 */
void Aretes::swap(int a1, int a2)
{
  int tmp;
  tmp = faces_(a1, 0);
  faces_(a1, 0) = faces_(a2, 0);
  faces_(a2, 0)=tmp;
  tmp = faces_(a1, 1);
  faces_(a1, 1) = faces_(a2, 1);
  faces_(a2, 1)=tmp;
  tmp = faces_(a1, 2);
  faces_(a1, 2) = faces_(a2, 2);
  faces_(a2, 2)=tmp;
  tmp = faces_(a1, 3);
  faces_(a1, 3) = faces_(a2, 3);
  faces_(a2, 3)=tmp;
  tmp = type1_(a1);
  type1_(a1)=type1_(a2);
  type1_(a2)=tmp;
  tmp = type2_(a1);
  type2_(a1)=type2_(a2);
  type2_(a2)=tmp;
}
/*! @brief Reorders the edge array: first corner edges (they have only two faces),
 *
 *  then boundary edges (they have three faces, two of which are boundary),
 *  then mixed edges (they have four faces, two of which are boundary),
 *  then internal edges (they have four internal faces).
 *
 */
void Aretes::trier(int& nb_aretes_coin, int& nb_aretes_bord,
                   int& nb_aretes_mixte, int& nb_aretes_interne)
{
  //
  nb_aretes_coin=nb_aretes_bord=nb_aretes_mixte=nb_aretes_interne=0;
  int coin = -1 ;
  int bord = 0 ;
  int mixte = 1 ;
  int interne = 2 ;
  int nb_aretes = type1_.size();
  int courante=0;
  int arete;
  while( (courante<nb_aretes)&&(type2_(courante)==coin) )
    {
      courante++;
      nb_aretes_coin++;
    }
  for(arete=courante; arete<nb_aretes; arete++)
    {
      if(type2_(arete)==coin)
        {
          swap(arete, courante);
          while( (courante<nb_aretes)&&(type2_(courante)==coin) )
            {
              courante++;
              nb_aretes_coin++;
            }
          //arete=courante;
        }
    }
  while( (courante<nb_aretes)&&(type2_(courante)==bord) )
    {
      courante++;
      nb_aretes_bord++;
    }
  for(arete=courante; arete<nb_aretes; arete++)
    {
      if(type2_(arete)==bord)
        {
          swap(arete, courante);
          assert(type2_(courante) == bord);
          while( (courante<nb_aretes)&&(type2_(courante)==bord) )
            {
              courante++;
              nb_aretes_bord++;
            }
          //arete=courante;
        }
    }
  while( (courante<nb_aretes)&&(type2_(courante)==mixte) )
    {
      courante++;
      nb_aretes_mixte++;
    }
  for(arete=courante; arete<nb_aretes; arete++)
    {
      if(type2_(arete)==mixte)
        {
          swap(arete, courante);
          assert(faces_(courante, 0) !=-1);
          assert(faces_(courante, 1) !=-1);
          assert(faces_(courante, 2) !=-1);
          assert(faces_(courante, 3) !=-1);
          while( (courante<nb_aretes)&&(type2_(courante)==mixte) )
            {
              courante++;
              nb_aretes_mixte++;
            }
          //arete=courante;
        }
    }
  while( (courante<nb_aretes)&&(type2_(courante)==interne) )
    {
      courante++;
      nb_aretes_interne++;
    }
  for(arete=courante; arete<nb_aretes; arete++)
    {
      if(type2_(arete)==interne)
        {
          swap(arete, courante);
          assert(faces_(courante, 0) !=-1);
          assert(faces_(courante, 1) !=-1);
          assert(faces_(courante, 2) !=-1);
          assert(faces_(courante, 3) !=-1);
          while( (courante<nb_aretes)&&(type2_(courante)==interne) )
            {
              courante++;
              nb_aretes_interne++;
            }
          //arete=courante;
        }
    }
}
 
void Aretes::trier_pour_debog(int& nb_aretes_coin, int& nb_aretes_bord,
                              int& nb_aretes_mixte, int& nb_aretes_interne,const DoubleTab&
                              xv)
{
 
  ArrOfDouble XVref(dimension),XVref2(dimension);
  int nb_aretes = type1_.size();
  int arete;
  for (int boucle=0; boucle<3; boucle++)
    {
      int deb=nb_aretes - nb_aretes_interne;
      int fin=nb_aretes;
      if (boucle==1)
        {
          fin=deb;
          deb-=nb_aretes_mixte;
        }
      else if (boucle==2)
        {
          fin=deb;
          deb-=  nb_aretes_bord;
        }
      for ( arete=deb; arete<fin; arete++)
        {
          if ((boucle==2)&&(type2_(arete)!=0))
            {
              Cerr<<"gros pb "<<arete<<finl;
              exit();
            }
          int ref=faces_(arete,0);
          for (int i=0; i<dimension; i++) XVref[i]=xv(ref,i);
          int marq=-1;
          int ref2;
          for (int arete2=arete; arete2<fin; arete2++)
            {
              ref2=faces_(arete2,0);
              for (int i=0; i<dimension; i++) XVref2[i]=xv(ref2,i);
              int test=-1;
              //              int testsa=-1;
              if (dimension==2)
                {
                  if (sup_strict(XVref2[1],XVref[1])) test=1;
                  else if ((XVref2[1]==XVref[1])&&(sup_strict(XVref2[0],XVref[0]))) test=0;
                }
              else if (dimension==3)
                {
                  if (XVref2[2]>XVref[2]) test=2;
                  else if  (XVref2[2]==XVref[2])
                    {
                      if (XVref2[1]>XVref[1]) test=1;
                      else if ((XVref2[1]==XVref[1])&&(XVref2[0]>XVref[0])) test=0;
                    }
                }
              if (test!=-1)
                {
                  marq=arete2;
                  XVref=XVref2;
                }
            }
          if (marq!=-1) swap(arete,marq);
        }
    }
  for (arete=nb_aretes - nb_aretes_interne; arete<nb_aretes*0; arete++)
    {
      Cerr <<me()<<" arete "<<arete<<" "<< faces_(arete,0)<<" "<< faces_(arete,1)<<" "<< faces_(arete,2)<<" "<< faces_(arete,3)<<finl;
    }
}