src/MEDLoader/MEDFileMeshElt.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "MEDFileMeshElt.hxx"
  22 #include "MEDFileMeshReadSelector.hxx"
  23
  24 #include "MEDCouplingUMesh.hxx"
  25
  26 #include "InterpKernelAutoPtr.hxx"
  27 #include "CellModel.hxx"
  28
  29 #include <iostream>
  30
  31 extern med_geometry_type typmai3[32];
  32
  33 using namespace ParaMEDMEM;
  34
  35 MEDFileUMeshPerType *MEDFileUMeshPerType::New(med_idt fid, const char *mName, int dt, int it, int mdim, med_geometry_type geoElt, INTERP_KERNEL::NormalizedCellType geoElt2, MEDFileMeshReadSelector *mrs)
  36 {
  37   med_entity_type whichEntity;
  38   if(!isExisting(fid,mName,dt,it,geoElt,whichEntity))
  39     return 0;
  40   return new MEDFileUMeshPerType(fid,mName,dt,it,mdim,geoElt,geoElt2,whichEntity,mrs);
  41 }
  42
  43 bool MEDFileUMeshPerType::isExisting(med_idt fid, const char *mName, int dt, int it, med_geometry_type geoElt, med_entity_type& whichEntity)
  44 {
  45   static const med_entity_type entities[3]={MED_CELL,MED_DESCENDING_FACE,MED_DESCENDING_EDGE};
  46   int nbOfElt=0;
  47   for(int i=0;i<3;i++)
  48     {
  49       med_bool changement,transformation;
  50       int tmp=MEDmeshnEntity(fid,mName,dt,it,entities[i],geoElt,MED_CONNECTIVITY,MED_NODAL,
  51                              &changement,&transformation);
  52       if(tmp>nbOfElt)
  53         {
  54           nbOfElt=tmp;
  55           whichEntity=entities[i];
  56           if(i>0)
  57             std::cerr << "WARNING : MEDFile has been detected to be no compilant with MED 3 : Please change entity in MEDFile for geotype " <<  geoElt << std::endl;
  58         }
  59     }
  60   return nbOfElt>0;
  61 }
  62
  63 int MEDFileUMeshPerType::getDim() const
  64 {
  65   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
  66   return cm.getDimension();
  67 }
  68
  69 MEDFileUMeshPerType::MEDFileUMeshPerType(med_idt fid, const char *mName, int dt, int it, int mdim, med_geometry_type geoElt, INTERP_KERNEL::NormalizedCellType type,
  70                                          med_entity_type entity, MEDFileMeshReadSelector *mrs):_type(type),_entity(entity)
  71 {
  72   med_bool changement,transformation;
  73   int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_CONNECTIVITY,MED_NODAL,
  74                                  &changement,&transformation);
  75   if(type!=INTERP_KERNEL::NORM_POLYGON && type!=INTERP_KERNEL::NORM_POLYHED)
  76     {
  77       loadFromStaticType(fid,mName,dt,it,mdim,curNbOfElem,geoElt,type,entity,mrs);
  78       return;
  79     }
  80   if(type==INTERP_KERNEL::NORM_POLYGON)
  81     {
  82       loadPolyg(fid,mName,dt,it,mdim,curNbOfElem,geoElt,entity,mrs);
  83       return;
  84     }
  85   //if(type==INTERP_KERNEL::NORM_POLYHED)
  86   loadPolyh(fid,mName,dt,it,mdim,curNbOfElem,geoElt,entity,mrs);
  87 }
  88
  89 void MEDFileUMeshPerType::loadFromStaticType(med_idt fid, const char *mName, int dt, int it, int mdim, int curNbOfElem, med_geometry_type geoElt, INTERP_KERNEL::NormalizedCellType type,
  90                                              med_entity_type entity, MEDFileMeshReadSelector *mrs)
  91 {
  92   _conn=DataArrayInt::New();
  93   int nbOfNodesPerCell=(geoElt%100);
  94   _conn->alloc((nbOfNodesPerCell+1)*curNbOfElem,1);
  95   _conn_index=DataArrayInt::New();
  96   _conn_index->alloc(curNbOfElem+1,1);
  97   INTERP_KERNEL::AutoPtr<int> connTab=new int[(nbOfNodesPerCell)*curNbOfElem];
  98   MEDmeshElementConnectivityRd(fid,mName,dt,it,entity,geoElt,MED_NODAL,MED_FULL_INTERLACE,connTab);
  99   loadCommonPart(fid,mName,dt,it,mdim,curNbOfElem,geoElt,entity,mrs);
 100   int *w1=_conn->getPointer();
 101   int *w2=_conn_index->getPointer();
 102   *w2++=0;
 103   const int *wi=connTab;
 104   for(int i=0;i<curNbOfElem;i++,wi+=nbOfNodesPerCell,w2++)
 105     {
 106       *w1++=(int)type;
 107       w1=std::transform(wi,wi+nbOfNodesPerCell,w1,std::bind2nd(std::plus<int>(),-1));
 108       *w2=w2[-1]+nbOfNodesPerCell+1;
 109     }
 110 }
 111
 112 void MEDFileUMeshPerType::loadCommonPart(med_idt fid, const char *mName, int dt, int it, int mdim, int curNbOfElem, med_geometry_type geoElt,
 113                                          med_entity_type entity, MEDFileMeshReadSelector *mrs)
 114 {
 115   med_bool changement,transformation;
 116   _fam=0;
 117   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_FAMILY_NUMBER,MED_NODAL,&changement,&transformation)>0)
 118     {
 119       if(!mrs || mrs->isCellFamilyFieldReading())
 120         {
 121           _fam=DataArrayInt::New();
 122           _fam->alloc(curNbOfElem,1);
 123           if(MEDmeshEntityFamilyNumberRd(fid,mName,dt,it,entity,geoElt,_fam->getPointer())!=0)
 124             std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
 125         }
 126     }
 127   _num=0;
 128   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_NUMBER,MED_NODAL,&changement,&transformation)>0)
 129     {
 130       if(!mrs || mrs->isCellNumFieldReading())
 131         {
 132           _num=DataArrayInt::New();
 133           _num->alloc(curNbOfElem,1);
 134           if(MEDmeshEntityNumberRd(fid,mName,dt,it,entity,geoElt,_num->getPointer())!=0)
 135             _num=0;
 136         }
 137     }
 138   _names=0;
 139   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_NAME,MED_NODAL,&changement,&transformation)>0)
 140     {
 141       if(!mrs || mrs->isCellNameFieldReading())
 142         {
 143           _names=DataArrayAsciiChar::New();
 144           _names->alloc(curNbOfElem+1,MED_SNAME_SIZE);//not a bug to avoid the memory corruption due to last \0 at the end
 145           if(MEDmeshEntityNameRd(fid,mName,dt,it,entity,geoElt,_names->getPointer())!=0)
 146             _names=0;
 147           else
 148             _names->reAlloc(curNbOfElem);//not a bug to avoid the memory corruption due to last \0 at the end
 149         }
 150     }
 151 }
 152
 153 void MEDFileUMeshPerType::loadPolyg(med_idt fid, const char *mName, int dt, int it, int mdim, int arraySize, med_geometry_type geoElt,
 154                                     med_entity_type entity, MEDFileMeshReadSelector *mrs)
 155 {
 156   med_bool changement,transformation;
 157   med_int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,entity,MED_POLYGON,MED_INDEX_NODE,MED_NODAL,&changement,&transformation)-1;
 158   _conn_index=DataArrayInt::New();
 159   _conn_index->alloc(curNbOfElem+1,1);
 160   _conn=DataArrayInt::New();
 161   _conn->alloc(arraySize+curNbOfElem,1);
 162   INTERP_KERNEL::AutoPtr<int> locConn=new int[arraySize];
 163   MEDmeshPolygonRd(fid,mName,dt,it,MED_CELL,MED_NODAL,_conn_index->getPointer(),locConn);
 164   int *w1=_conn->getPointer();
 165   int *w2=_conn_index->getPointer();
 166   const int *wi=locConn;
 167   for(int i=0;i<curNbOfElem;i++,w2++)
 168     {
 169       *w1++=(int)INTERP_KERNEL::NORM_POLYGON;
 170       const int *wi2=wi+(w2[1]-w2[0]);
 171       w1=std::transform(wi,wi2,w1,std::bind2nd(std::plus<int>(),-1));
 172       wi=wi2;
 173       *w2=*w2-1+i;
 174     }
 175   *w2=*w2-1+curNbOfElem;
 176   loadCommonPart(fid,mName,dt,it,mdim,curNbOfElem,MED_POLYGON,entity,mrs);
 177 }
 178
 179 void MEDFileUMeshPerType::loadPolyh(med_idt fid, const char *mName, int dt, int it, int mdim, int connFaceLgth, med_geometry_type geoElt,
 180                                     med_entity_type entity, MEDFileMeshReadSelector *mrs)
 181 {
 182   med_bool changement,transformation;
 183   med_int indexFaceLgth=MEDmeshnEntity(fid,mName,dt,it,MED_CELL,MED_POLYHEDRON,MED_INDEX_NODE,MED_NODAL,&changement,&transformation);
 184   int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,MED_CELL,MED_POLYHEDRON,MED_INDEX_FACE,MED_NODAL,&changement,&transformation)-1;
 185   INTERP_KERNEL::AutoPtr<int> index=new int[curNbOfElem+1];
 186   INTERP_KERNEL::AutoPtr<int> indexFace=new int[indexFaceLgth];
 187   INTERP_KERNEL::AutoPtr<int> locConn=new int[connFaceLgth];
 188   MEDmeshPolyhedronRd(fid,mName,dt,it,MED_CELL,MED_NODAL,index,indexFace,locConn);
 189   int arraySize=connFaceLgth;
 190   for(int i=0;i<curNbOfElem;i++)
 191     arraySize+=index[i+1]-index[i]-1;
 192   _conn=DataArrayInt::New();
 193   _conn->alloc(arraySize+curNbOfElem,1);
 194   int *wFinalConn=_conn->getPointer();
 195   _conn_index=DataArrayInt::New();
 196   _conn_index->alloc(curNbOfElem+1,1);
 197   int *finalIndex=_conn_index->getPointer();
 198   finalIndex[0]=0;
 199   for(int i=0;i<curNbOfElem;i++)
 200     {
 201       *wFinalConn++=(int)INTERP_KERNEL::NORM_POLYHED;
 202       finalIndex[i+1]=finalIndex[i]+index[i+1]-index[i]-1+indexFace[index[i+1]-1]-indexFace[index[i]-1]+1;
 203       wFinalConn=std::transform(locConn+indexFace[index[i]-1]-1,locConn+indexFace[index[i]]-1,wFinalConn,std::bind2nd(std::plus<int>(),-1));
 204       for(int j=index[i];j<index[i+1]-1;j++)
 205         {
 206           *wFinalConn++=-1;
 207           wFinalConn=std::transform(locConn+indexFace[j]-1,locConn+indexFace[j+1]-1,wFinalConn,std::bind2nd(std::plus<int>(),-1));
 208         }
 209     }
 210   loadCommonPart(fid,mName,dt,it,mdim,curNbOfElem,MED_POLYHEDRON,entity,mrs);
 211 }
 212
 213 void MEDFileUMeshPerType::write(med_idt fid, const char *mname, int mdim, const MEDCouplingUMesh *m, const DataArrayInt *fam, const DataArrayInt *num, const DataArrayAsciiChar *names)
 214 {
 215   int nbOfCells=m->getNumberOfCells();
 216   if(nbOfCells<1)
 217     return ;
 218   int dt,it;
 219   double timm=m->getTime(dt,it);
 220   INTERP_KERNEL::NormalizedCellType ikt=m->getTypeOfCell(0);
 221   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(ikt);
 222   med_geometry_type curMedType=typmai3[(int)ikt];
 223   const int *conn=m->getNodalConnectivity()->getConstPointer();
 224   const int *connI=m->getNodalConnectivityIndex()->getConstPointer();
 225   if(ikt!=INTERP_KERNEL::NORM_POLYGON && ikt!=INTERP_KERNEL::NORM_POLYHED)
 226     {
 227       int nbNodesPerCell=cm.getNumberOfNodes();
 228       INTERP_KERNEL::AutoPtr<int> tab=new int[nbNodesPerCell*nbOfCells];
 229       int *w=tab;
 230       for(int i=0;i<nbOfCells;i++)
 231         w=std::transform(conn+connI[i]+1,conn+connI[i+1],w,std::bind2nd(std::plus<int>(),1));
 232       MEDmeshElementConnectivityWr(fid,mname,dt,it,timm,MED_CELL,curMedType,MED_NODAL,MED_FULL_INTERLACE,nbOfCells,tab);
 233     }
 234   else
 235     {
 236       if(ikt==INTERP_KERNEL::NORM_POLYGON)
 237         {
 238           INTERP_KERNEL::AutoPtr<int> tab1=new int[nbOfCells+1];
 239           INTERP_KERNEL::AutoPtr<int> tab2=new int[m->getMeshLength()];
 240           int *wI=tab1; *wI=1;
 241           int *w=tab2;
 242           for(int i=0;i<nbOfCells;i++,wI++)
 243             {
 244               wI[1]=wI[0]+connI[i+1]-connI[i]-1;
 245               w=std::transform(conn+connI[i]+1,conn+connI[i+1],w,std::bind2nd(std::plus<int>(),1));
 246             }
 247           MEDmeshPolygonWr(fid,mname,dt,it,timm,MED_CELL,MED_NODAL,nbOfCells+1,tab1,tab2);
 248         }
 249       else
 250         {
 251           int meshLgth=m->getMeshLength();
 252           int nbOfFaces=std::count(conn,conn+meshLgth,-1)+nbOfCells;
 253           INTERP_KERNEL::AutoPtr<int> tab1=new int[nbOfCells+1];
 254           int *w1=tab1; *w1=1;
 255           INTERP_KERNEL::AutoPtr<int> tab2=new int[nbOfFaces+1];
 256           int *w2=tab2; *w2=1;
 257           INTERP_KERNEL::AutoPtr<int> bigtab=new int[meshLgth-nbOfCells];
 258           int *bt=bigtab;
 259           for(int i=0;i<nbOfCells;i++,w1++)
 260             {
 261               int nbOfFaces2=0;
 262               for(const int *w=conn+connI[i]+1;w!=conn+connI[i+1];w2++)
 263                 {
 264                   const int *wend=std::find(w,conn+connI[i+1],-1);
 265                   bt=std::transform(w,wend,bt,std::bind2nd(std::plus<int>(),1));
 266                   int nbOfNode=std::distance(w,wend);
 267                   w2[1]=w2[0]+nbOfNode;
 268                   if(wend!=conn+connI[i+1])
 269                     w=wend+1;
 270                   else
 271                     w=wend;
 272                   nbOfFaces2++;
 273                 }
 274               w1[1]=w1[0]+nbOfFaces2;
 275             }
 276           MEDmeshPolyhedronWr(fid,mname,dt,it,timm,MED_CELL,MED_NODAL,nbOfCells+1,tab1,nbOfFaces+1,tab2,bigtab);
 277         }
 278     }
 279   if(fam)
 280     MEDmeshEntityFamilyNumberWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,fam->getConstPointer());
 281   if(num)
 282     MEDmeshEntityNumberWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,num->getConstPointer());
 283   if(names)
 284     {
 285       if(names->getNumberOfComponents()!=MED_SNAME_SIZE)
 286         {
 287           std::ostringstream oss; oss << "MEDFileUMeshPerType::write : expected a name field on cells with number of components set to " << MED_SNAME_SIZE;
 288           oss << " ! The array has " << names->getNumberOfComponents() << " components !";
 289           throw INTERP_KERNEL::Exception(oss.str().c_str());
 290         }
 291       MEDmeshEntityNameWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,names->getConstPointer());
 292     }
 293 }