Convert_ICoCoTrioField.cpp

/****************************************************************************
* Copyright (c) 2025, 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 <Convert_ICoCoTrioField.h>
#include <ICoCoTrioField.h>
#include <ICoCoMEDDoubleField.hxx>
#include <Domaine.h>
#include <Champ_Generique_base.h>
#include <Domaine_VF.h>
#include <PE_Groups.h>
#include <Comm_Group.h>
#include <Polyedre.h>
 
void affecte_double_avec_doubletab(double** p, const ArrOfDouble& trio)
{
  *p=new double[trio.size_array()];
  memcpy(*p,trio.addr(),trio.size_array()*sizeof(double));
}
 
void affecte_int_avec_inttab(int** p, const ArrOfInt& trio)
{
  int sz=trio.size_array();
  *p=new int[sz];
  memcpy(*p,trio.addr(),sz*sizeof(int));
}
 
void build_triofield(const Champ_Generique_base& ch, ICoCo::TrioField& afield)
{
  build_triomesh(ch.get_ref_domaine_dis_base(), afield, ch.get_localisation() == Entity::NODE, ch.get_localisation() == Entity::FACE);
 
  afield.setName(ch.le_nom().getString());
  afield._time1 = afield._time2 = ch.get_time(), afield._itnumber = 0;
 
  /* copy of field values */
  afield._has_field_ownership = true;
  OWN_PTR(Champ_base) espace_stockage;
  const Champ_base& champ_ecriture = ch.get_champ(espace_stockage);
  const DoubleTab& vals = champ_ecriture.valeurs();
  afield._nb_field_components = vals.nb_dim() > 1 ? vals.dimension(1) : 1;
  affecte_double_avec_doubletab(&afield._field, vals);
}
 
void build_triofield(const Champ_base& ch, const Domaine_dis_base& dom_dis, ICoCo::TrioField& afield)
{
  build_triomesh(dom_dis, afield, 0, 0); // XXX seulement elem...
 
  afield.setName(ch.le_nom().getString());
  afield._time1 = afield._time2 = 0.0, afield._itnumber = 0;
 
  /* copy of field values */
  afield._has_field_ownership = true;
  const DoubleTab& vals = ch.valeurs();
  afield._nb_field_components = vals.nb_dim() > 1 ? vals.dimension(1) : 1;
  affecte_double_avec_doubletab(&afield._field, vals);
}
 
void build_triomesh(const Domaine_dis_base& dom_dis, ICoCo::TrioField& afield, int type, int loc_faces)
{
  const Domaine_VF& zvf = ref_cast(Domaine_VF, dom_dis);
  const Domaine& dom = dom_dis.domaine();
  afield.clear();
  afield._type = type;
 
  /* node array: copy from the domain */
  const DoubleTab& coord = dom.les_sommets();
  afield._space_dim = dom.dimension;
  afield._nbnodes = coord.dimension(0);
  affecte_double_avec_doubletab(&afield._coords, coord);
 
  /* dimension of the domain elements */
  Motcle type_elem_ = dom_dis.domaine().type_elem()->que_suis_je();
  Motcle type_elem(type_elem_);
  type_elem.prefix("_AXI");
  if (type_elem != Motcle(type_elem_))
    {
      if (type_elem == "QUADRILATERE_2D")
        type_elem = "SEGMENT_2D";
      if (type_elem == "RECTANGLE_2D")
        type_elem = "RECTANGLE";
    }
  if ((type_elem == "RECTANGLE") || (type_elem == "QUADRANGLE") || (type_elem == "TRIANGLE") || (type_elem == "TRIANGLE_3D") || (type_elem == "QUADRANGLE_3D") || (type_elem == "POLYGONE") || (type_elem == "POLYGONE_3D"))
    afield._mesh_dim=2;
  else if ((type_elem == "HEXAEDRE") || (type_elem == "HEXAEDRE_VEF") || (type_elem == "POLYEDRE") || (type_elem == "PRISME") || (type_elem == "TETRAEDRE"))
    afield._mesh_dim=3;
  else if ((type_elem == "SEGMENT_2D") || (type_elem == "SEGMENT"))
    afield._mesh_dim=1;
  else
    {
      Cerr << "build_triofield: " << type_elem<< " not coded" <<finl;
      Process::exit();
    }
 
  /* elements: those of the domain if the field is at nodes/elements, or faces if the field is at faces */
  if (loc_faces) afield._mesh_dim--;
  afield._nb_elems = loc_faces ? zvf.nb_faces() : dom_dis.domaine().nb_elem();
  if (loc_faces || type_elem != "POLYEDRE") //face mesh -> connectivity = face_sommets
    {
      const IntTab& conn = loc_faces ? dom_dis.face_sommets() : dom_dis.domaine().les_elems();
      //the only way to prevent polygons from being mistaken for quadrilaterals is to have a connectivity array of width > 4...
      afield._nodes_per_elem = std::max(conn.dimension(1), type_elem == "POLYGONE" || type_elem == "POLYGONE_3D"  || type_elem == "POLYEDRE"  ? (int) 5 : 0);
      afield._connectivity = new int[afield._nb_elems * afield._nodes_per_elem];
      for (int i = 0; i < afield._nb_elems; i++)
        for (int j = 0; j < afield._nodes_per_elem; j++)
          afield._connectivity[afield._nodes_per_elem * i + j] = j < conn.dimension(1) ? conn(i, j) : -1;
    }
  else //polyhedron mesh -> connectivity in MEDCoupling format, to be built manually
    {
      const Polyedre& poly = ref_cast(Polyedre, dom.type_elem().valeur());
      const ArrOfInt& e_fi = poly.getPolyhedronIndex(), &f_si = poly.getFacesIndex(), &sl = poly.getNodes();
      const IntTab& e_s = dom.les_elems();
      int e, f, s, nef_max = 0, nfs_max = 0; //max nb faces/elem and nodes/face
      for (e = 0; e + 1 < e_fi.size_array(); e++) nef_max = std::max(nef_max, e_fi(e + 1) - e_fi(e));
      for (f = 0; f + 1 < f_si.size_array(); f++) nfs_max = std::max(nfs_max, f_si(f + 1) - f_si(f));
      afield._nodes_per_elem = std::max(nef_max * (nfs_max + 1), (int) 9); //one -1 after each face: at least 9 to avoid butterfly cells
      int *p = afield._connectivity = new int[afield._nb_elems * afield._nodes_per_elem];
      for (e = 0; e < afield._nb_elems; e++)
        {
          /* insertion of the connectivity of each face, followed by a -1 */
          for (f = e_fi(e); f < e_fi(e + 1); f++, *p = -1, p++)
            for (s = f_si(f); s < f_si(f + 1); s++)
              *p = e_s(e, sl(s)), p++;
          /* fill with -1 up to the next row */
          for ( ; p < afield._connectivity + (e + 1) * afield._nodes_per_elem; p++) *p = -1;
        }
    }
}
 
#ifndef NO_MEDFIELD
#include <MEDCouplingUMesh.hxx>
#include <MEDCouplingFieldDouble.hxx>
#include <MCAuto.hxx>
#include <string.h>
#include <vector>
 
 
using ICoCo::TrioField;
using ICoCo::MEDDoubleField;
using std::vector;
 
 
/*!
 * @brief Build a MEDDoubleField from a TrioField.
 *
 * This method is non const only due to this->_field that can be modified (to point to the same domaine than returned object).
 * So \b warning, to access to \a this->_field only when the returned object is alive.
 * @param triofield the TrioField to convert
 * @return the corresponding MEDDoubleField
 */
MEDDoubleField build_medfield(TrioField& triofield)
{
  MEDCoupling::MCAuto<MEDCoupling::MEDCouplingUMesh> mesh(MEDCoupling::MEDCouplingUMesh::New("",triofield._mesh_dim));
  MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> coo(MEDCoupling::DataArrayDouble::New());
  coo->alloc(triofield._nbnodes,triofield._space_dim);
  mesh->setCoords(coo);
  double *ptr(coo->getPointer());
  std::copy(triofield._coords,triofield._coords+triofield._space_dim*triofield._nbnodes,ptr);
  mesh->allocateCells(triofield._nb_elems);
  INTERP_KERNEL::NormalizedCellType elemtype;
  switch(triofield._mesh_dim)
    {
    case 0 :
      {
        switch (triofield._nodes_per_elem)
          {
          case 0: // cas field vide
            elemtype=INTERP_KERNEL::NORM_SEG2; // to avoid warning
            break;
          default:
            throw INTERP_KERNEL::Exception("incompatible Trio field - wrong nb of nodes per elem");
          }
        break;
      }
    case 1:
      {
        switch (triofield._nodes_per_elem)
          {
          case 2:
            elemtype=INTERP_KERNEL::NORM_SEG2;
            break;
          default:
            throw INTERP_KERNEL::Exception("incompatible Trio field - wrong nb of nodes per elem");
          }
        break;
      }
    case 2:
      {
        switch (triofield._nodes_per_elem)
          {
          case 3:
            elemtype=INTERP_KERNEL::NORM_TRI3;
            break;
          case 4 :
            elemtype=INTERP_KERNEL::NORM_QUAD4;
            break;
          default:
            elemtype=INTERP_KERNEL::NORM_POLYGON;
          }
        break;
      }
    case 3:
      {
        switch (triofield._nodes_per_elem)
          {
          case 4:
            elemtype=INTERP_KERNEL::NORM_TETRA4;
            break;
          case 8 :
            elemtype=INTERP_KERNEL::NORM_HEXA8;
            break;
          default:
            elemtype=INTERP_KERNEL::NORM_POLYHED;
          }
        break;
      default:
        throw INTERP_KERNEL::Exception("incompatible Trio field - wrong mesh dimension");
      }
    }
  //creating a connectivity table that complies to MED (1 indexing) <- actually no
  //and passing it to _mesh
  MEDCoupling::MCAuto<MEDCoupling::MEDCouplingFieldDouble> field;
  int *conn(new int[triofield._nodes_per_elem]);
  for (int i=0; i<triofield._nb_elems; i++)
    {
 
      for(int j=0; j<triofield._nodes_per_elem; j++)
        {
          conn[j]=triofield._connectivity[i*triofield._nodes_per_elem+j];
        }
      if (elemtype==INTERP_KERNEL::NORM_QUAD4)
        {
          // different node ordering in TRUST
          int tmp=conn[3];
          conn[3]=conn[2];
          conn[2]=tmp;
        }
      if (elemtype==INTERP_KERNEL::NORM_HEXA8)
        {
          // different node ordering in TRUST
          int tmp=conn[3];
          conn[3]=conn[2];
          conn[2]=tmp;
          tmp=conn[7];
          conn[7]=conn[6];
          conn[6]=tmp;
        }
      int size = triofield._nodes_per_elem;
      while (conn[size - 1] == -1) size--; //remove trailing -1 values from connectivity
#if INT_is_64_ == 2
      // Convert int into mcIdType
      std::vector<mcIdType> conn_mc(conn, conn+size);
      mesh->insertNextCell(elemtype,size,conn_mc.data());
#else
      mesh->insertNextCell(elemtype,size,conn);
#endif
    }
  delete [] conn;
  mesh->finishInsertingCells();
  //
 
 
  std::vector<int> cells;
  // if ((mesh->getSpaceDimension() == 2 || mesh->getSpaceDimension() == 3) && mesh->getMeshDimension() == 2)
  //   mesh->checkButterflyCells(cells);
  if (!cells.empty())
    {
      Cerr<<" cells are butterflyed "<<cells[0]<<finl;
      PE_Groups::groupe_TRUST().abort();
      Process::abort();
      Process::exit();
    }
  //field on the sending end
  int nb_case=triofield.nb_values();
  if (triofield._type==0)
    {
      field =  MEDCoupling::MEDCouplingFieldDouble::New(MEDCoupling::ON_CELLS,MEDCoupling::ONE_TIME);
    }
  else
    {
      field =  MEDCoupling::MEDCouplingFieldDouble::New(MEDCoupling::ON_NODES,MEDCoupling::ONE_TIME );
    }
  field->setMesh(mesh);
  field->setNature(MEDCoupling::IntensiveMaximum);
  MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> fieldArr(MEDCoupling::DataArrayDouble::New());
  fieldArr->alloc(field->getNumberOfTuplesExpected(),triofield._nb_field_components);
  field->setName(triofield.getName());
  std::string meshName("SupportOf_");
  meshName+=triofield.getName();
  mesh->setName(meshName);
  field->setTime(triofield._time1,0,triofield._itnumber);
  if (triofield._field!=0)
    {
      for (int i =0; i<nb_case; i++)
        for (int j=0; j<triofield._nb_field_components; j++)
          {
            fieldArr->setIJ(i,j,triofield._field[i*triofield._nb_field_components+j]);
          }
    }
  //field on the receiving end
  else
    {
      // the trio field points to the pointer inside the MED field
      triofield._field=fieldArr->getPointer();
      for (int i=0; i<triofield._nb_field_components*nb_case; i++)
        triofield._field[i]=0.0;
    }
  field->setArray(fieldArr);
  return MEDDoubleField(field);
}
MEDDoubleField build_medfield(const Champ_Generique_base& ch)
{
  TrioField fl;
  build_triofield(ch, fl);
  return build_medfield(fl);
}
 
 
#else
namespace ICoCo
{
class MEDDoubleField
{
};
}
ICoCo::MEDDoubleField build_medfield(ICoCo::TrioField& ch)
{
  Cerr<<"Version compiled without MEDCoupling"<<finl;
  Process::exit();
  throw;
}
ICoCo::MEDDoubleField build_medfield(const Champ_Generique_base& ch)
{
  Cerr<<"Version compiled without MEDCoupling"<<finl;
  Process::exit();
  throw;
}
#endif