Champ_front_lu.cpp

/****************************************************************************
* Copyright (c) 2024, 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 <Champ_front_lu.h>
#include <Domaine.h>
#include <LecFicDiffuse.h>
#include <Interprete.h>
#include <Domaine_VF.h>
 
Implemente_instanciable(Champ_front_lu,"Champ_front_lu",Ch_front_var_stationnaire);
// XD champ_front_lu front_field_base champ_front_lu NO_BRACE boundary field which is given from data issued from a read
// XD_CONT file. The format of this file has to be the same that the one generated by Ecrire_fichier_xyz_valeur NL2
// XD_CONT Example for K and epsilon quantities to be defined for inlet condition in a boundary named \'entree\':
// XD_CONT NL2entree frontiere_ouverte_K_Eps_impose Champ_Front_lu dom 2pb_K_EPS_PERIO_1006.306198.dat
// XD attr domaine ref_domaine domain REQ Name of domain
// XD attr dim entier dim REQ number of components
// XD attr file chaine file REQ path for the read file
 
 
/*! @brief Print the field to an output stream.
 *
 * Prints the field size and the (constant) value on
 *     the boundary.
 *
 * @param (Sortie& os) an output stream
 * @return (Sortie&) the modified output stream
 */
Sortie& Champ_front_lu::printOn(Sortie& os) const
{
  const DoubleTab& tab=valeurs();
  os << tab.size() << " ";
  for(int i=0; i<tab.size(); i++)
    os << tab(0,i);
  return os;
}
 
 
/*! @brief Read the field from an input stream.
 *
 * Format:
 *       Champ_front_lu nb_compo vrel_1 ... [vrel_i]
 *
 * @param (Entree& is) an input stream
 * @return (Entree& is) the modified input stream
 */
Entree& Champ_front_lu::readOn(Entree& is)
{
  int dim;
  //Motcle motlu;
  Nom nom;
  is >> nom;
 
  mon_domaine = ref_cast(Domaine, Interprete::objet(nom));
 
 
 
  dim=lire_dimension(is,que_suis_je());
  fixer_nb_comp(dim);
  is >> nom_fic;
  // Reading values from file fic
  return is;
}
 
void erreur_chfront(const Nom& nom_fic, const int dim, const Nom& nom_domaine)
{
  Cerr << "****************************************************************************************************************" << finl;
  Cerr << "Error! It seems that several faces into the file " << nom_fic << " does not match any faces of the current mesh." << finl;
  Cerr << finl;
  Cerr << "First, check that the file " << nom_fic << " can be used to define a boundary field on the specified boundary of the current mesh." << finl;
  Cerr << "Second, try to switch to a different algorithm in your data file (find the nearest geometric face of the mesh for each face defined in the " << nom_fic << " file):" << finl;
  Cerr << "Replace:" <<finl;
  Cerr << "Champ_front_lu " << nom_domaine << " " << dim << " " << nom_fic << finl;
  Nom un((Nom)dim),zero((Nom)dim);
  for (int i=0; i<dim; i++)
    {
      un+=" 1";
      zero+=" 0";
    }
  Cerr << "By:"<<finl;
  Nom chaine="moyenne_imposee connexion_approchee fichier";
  Cerr << "Champ_front_recyclage { pb_champ_evaluateur problem_name field_name "<<dim<<finl;
  Cerr << "                        "<<chaine<<" "<<nom_fic<<finl;
  Cerr << "                        ampli_moyenne_imposee "<<un<<finl;
  Cerr << "                        ampli_moyenne_recyclee "<<zero<<finl;
  Cerr << "                        ampli_fluctuation "<<zero<<" }" <<finl;
  Process::exit();
}
 
int Champ_front_lu::initialiser(double temps, const Champ_Inc_base& inco)
{
 
  if (!Ch_front_var_stationnaire::initialiser(temps,inco))
    return 0;
 
  int dim=nb_comp();
  Domaine& domaine=mon_domaine.valeur();
  const Frontiere_dis_base& fr_dis=frontiere_dis();
  const Frontiere& frontiere=fr_dis.frontiere();
  const int nb_faces=frontiere.nb_faces();
  const int ndeb=frontiere.num_premiere_face();
  DoubleTab& les_val=valeurs();
 
  LecFicDiffuse fic;
  int bin=is_bin(nom_fic);
  fic.set_bin(bin);
  fic.ouvrir(nom_fic);
  Nom format;
  if (bin)
    fic >> format;
  else
    format="ascii";
 
  Cerr << "Reading of the " << format << " file : " << nom_fic << finl;
 
  DoubleVect point(dimension);
  DoubleVect val_lu(dim);
  int nb_val_lues;
  fic >> nb_val_lues;
 
  // Verification (the file may contain MORE faces than the boundary, e.g. in periodic cases)
  trustIdType nb_faces_tot = mp_sum(nb_faces);
  if ((trustIdType)nb_val_lues<nb_faces_tot)
    {
      Cerr << "This file contains " << nb_val_lues << " faces which is not enough compare to" << finl;
      Cerr << "the faces number (" << nb_faces_tot << ") of the boundary " << frontiere.le_nom() << finl;
      exit();
    }
 
  // PQ 20/06/06 : to enable a "fast" search of faces
  //               we rely on mon_dom.chercher_elements(x,y,z)
  //               (method inspired by EcritureLectureSpecial::lecture_special)
  //
  //               Restriction --> face fields only
 
  int nb_elems = domaine.nb_elem();
  int nbfacelem=domaine.nb_faces_elem();
  const Domaine_VF& zvf = ref_cast(Domaine_VF, domaine_dis());
  const IntTab& elem_faces=zvf.elem_faces();
  const DoubleTab& xv=zvf.xv();
  IntTab compteur(nb_faces);
  compteur=0;
 
  int k;
  double x,y,z=0;
  int pourcent=0,tmp;
 
  for (int p=0; p<nb_val_lues; p++)
    {
      tmp=(p*10)/nb_val_lues;
      if (tmp>pourcent)
        {
          pourcent=tmp;
          Cerr<<"\rProcess "<<Process::me()<<" has read "<<pourcent*10<<"% of values"<<flush;
        }
      for (k=0; k<dimension; k++)
        fic >> point[k];
      for (k=0; k<dim; k++)
        fic >>val_lu[k];
 
      x=point[0];
      y=point[1];
      if (Objet_U::dimension==3) z=point[2];
      int Case=-1;
      // First algorithm (the fastest) searching for the element
      // containing the face centre of gravity
      int num_elem=domaine.chercher_elements(x,y,z);
      if ((num_elem!=-1) && (num_elem<nb_elems))
        {
          // Find the face of the element
          for (int fac=0; fac<nbfacelem; fac++)
            {
              int face=elem_faces(num_elem,fac);
              int ok=1;
              for (int dir=0; dir<dimension; dir++)
                {
                  ok=ok&&(est_egal(point[dir],xv(face,dir)));
                }
              if (ok)
                {
                  Case=face;
                  break;
                }
            }
        }
      else
        {
          // Second algorithm (slower) used when the point is outside the domain
          // In that case, we scan the faces of the boundary
          for (int face=0; face<nb_faces; face++)
            if (compteur[face] == 0)
              {
                int ok=1;
                for (int dir=0; dir<dimension; dir++)
                  ok=ok&&(est_egal(point[dir],xv(ndeb+face,dir)));
                if (ok)
                  {
                    Case=ndeb+face;
                    break;
                  }
              }
        }
      // Face found
      if (Case!=-1)
        {
          int face_loc=Case-ndeb;
          if ((face_loc<0) || (face_loc>=nb_faces))
            {
              Cerr << " The face " << ndeb << " found doesn't belong to the boundary sought " << finl;
              exit();
            }
          compteur[face_loc] += 1;
          for (k=0; k<dim; k++)
            les_val(face_loc,k) = val_lu[k];
        }
    }
  Cerr<<"\rProcess "<<Process::me()<<" has read 100% of values"<<finl;
 
  // We verify that all faces of the boundary have been assigned
  bool err=false;
  for (int face=0; face<nb_faces; face++)
    if (compteur[face] != 1)
      {
        Cerr.precision(12);
        if (compteur[face] == 0)
          Cerr << "No value found for the face " << face << "/" << nb_faces-1 << " of the boundary " << frontiere.le_nom() << " on the processor " << Process::me() << ".";
        else
          Cerr << "Several values found for the face " << face << "/" << nb_faces-1 << " of the boundary " << frontiere.le_nom() << " on the processor " << Process::me() << ".";
        Cerr << " Coordinates: " << xv(ndeb+face,0) << " " << xv(ndeb+face,1);
        if (dimension==3) Cerr << " " << xv(ndeb+face,2);
        Cerr << finl;
        err=true;
      }
  if (Process::mp_or(err))
    erreur_chfront(nom_fic, dim, mon_domaine->le_nom());
  les_val.echange_espace_virtuel();
  return 1;
}
 
/*! @brief Returns the format of the read file
 *
 * @param (Nom name) an input stream
 * @return (int)
 */
int Champ_front_lu::is_bin(Nom name)
{
  LecFicDiffuse fichier(name);
  Nom format;
  fichier >> format;
 
#ifdef INT_is_64_
  // In case of INT_64, first Nom can be INT64 ...
  if(format=="INT64")
    format="binary";
#endif
 
  if(format=="binary")
    return 1;
  else
    return 0;
}