src/MEDLoader/MEDFileMeshElt.cxx

   1 // Copyright (C) 2007-2014  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, or (at your option) any later version.
   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[34];
  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 std::size_t MEDFileUMeshPerType::getHeapMemorySizeWithoutChildren() const
  44 {
  45   return 0;
  46 }
  47
  48 std::vector<const BigMemoryObject *> MEDFileUMeshPerType::getDirectChildren() const
  49 {
  50   std::vector<const BigMemoryObject *> ret;
  51   if((const MEDCoupling1GTUMesh *)_m)
  52     ret.push_back((const MEDCoupling1GTUMesh *)_m);
  53   if((const DataArrayInt *)_num)
  54     ret.push_back((const DataArrayInt *)_num);
  55   if((const DataArrayInt *)_fam)
  56     ret.push_back((const DataArrayInt *)_fam);
  57   if((const DataArrayAsciiChar *)_names)
  58     ret.push_back((const DataArrayAsciiChar *)_names);
  59   return ret;
  60 }
  61
  62 bool MEDFileUMeshPerType::isExisting(med_idt fid, const char *mName, int dt, int it, med_geometry_type geoElt, med_entity_type& whichEntity)
  63 {
  64   static const med_entity_type entities[3]={MED_CELL,MED_DESCENDING_FACE,MED_DESCENDING_EDGE};
  65   int nbOfElt=0;
  66   for(int i=0;i<3;i++)
  67     {
  68       med_bool changement,transformation;
  69       int tmp=MEDmeshnEntity(fid,mName,dt,it,entities[i],geoElt,MED_CONNECTIVITY,MED_NODAL,
  70           &changement,&transformation);
  71       if(tmp>nbOfElt)
  72         {
  73           nbOfElt=tmp;
  74           whichEntity=entities[i];
  75           if(i>0)
  76             std::cerr << "WARNING : MEDFile has been detected to be no compilant with MED 3 : Please change entity in MEDFile for geotype " <<  geoElt << std::endl;
  77         }
  78     }
  79   return nbOfElt>0;
  80 }
  81
  82 int MEDFileUMeshPerType::getDim() const
  83 {
  84   return _m->getMeshDimension();
  85 }
  86
  87 MEDFileUMeshPerType::MEDFileUMeshPerType(med_idt fid, const char *mName, int dt, int it, int mdim, med_geometry_type geoElt, INTERP_KERNEL::NormalizedCellType type,
  88                                          med_entity_type entity, MEDFileMeshReadSelector *mrs):_entity(entity)
  89 {
  90   med_bool changement,transformation;
  91   int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_CONNECTIVITY,MED_NODAL,
  92       &changement,&transformation);
  93   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
  94   if(!cm.isDynamic())
  95     {
  96       loadFromStaticType(fid,mName,dt,it,mdim,curNbOfElem,geoElt,type,entity,mrs);
  97       return;
  98     }
  99   if(type==INTERP_KERNEL::NORM_POLYGON || type==INTERP_KERNEL::NORM_QPOLYG)
 100     {
 101       loadPolyg(fid,mName,dt,it,mdim,curNbOfElem,geoElt,entity,mrs);
 102       return;
 103     }
 104   //if(type==INTERP_KERNEL::NORM_POLYHED)
 105   loadPolyh(fid,mName,dt,it,mdim,curNbOfElem,geoElt,entity,mrs);
 106 }
 107
 108 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,
 109                                              med_entity_type entity, MEDFileMeshReadSelector *mrs)
 110 {
 111   _m=MEDCoupling1SGTUMesh::New(mName,type);
 112   MEDCoupling1SGTUMesh *mc(dynamic_cast<MEDCoupling1SGTUMesh *>((MEDCoupling1GTUMesh *)_m));
 113   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
 114   int nbOfNodesPerCell=mc->getNumberOfNodesPerCell();
 115   conn->alloc(nbOfNodesPerCell*curNbOfElem,1);
 116   MEDmeshElementConnectivityRd(fid,mName,dt,it,entity,geoElt,MED_NODAL,MED_FULL_INTERLACE,conn->getPointer());
 117   std::transform(conn->begin(),conn->end(),conn->getPointer(),std::bind2nd(std::plus<int>(),-1));
 118   mc->setNodalConnectivity(conn);
 119   loadCommonPart(fid,mName,dt,it,mdim,curNbOfElem,geoElt,entity,mrs);
 120 }
 121
 122 void MEDFileUMeshPerType::loadCommonPart(med_idt fid, const char *mName, int dt, int it, int mdim, int curNbOfElem, med_geometry_type geoElt,
 123                                          med_entity_type entity, MEDFileMeshReadSelector *mrs)
 124 {
 125   med_bool changement,transformation;
 126   _fam=0;
 127   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_FAMILY_NUMBER,MED_NODAL,&changement,&transformation)>0)
 128     {
 129       if(!mrs || mrs->isCellFamilyFieldReading())
 130         {
 131           _fam=DataArrayInt::New();
 132           _fam->alloc(curNbOfElem,1);
 133           if(MEDmeshEntityFamilyNumberRd(fid,mName,dt,it,entity,geoElt,_fam->getPointer())!=0)
 134             std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
 135         }
 136     }
 137   _num=0;
 138   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_NUMBER,MED_NODAL,&changement,&transformation)>0)
 139     {
 140       if(!mrs || mrs->isCellNumFieldReading())
 141         {
 142           _num=DataArrayInt::New();
 143           _num->alloc(curNbOfElem,1);
 144           if(MEDmeshEntityNumberRd(fid,mName,dt,it,entity,geoElt,_num->getPointer())!=0)
 145             _num=0;
 146         }
 147     }
 148   _names=0;
 149   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_NAME,MED_NODAL,&changement,&transformation)>0)
 150     {
 151       if(!mrs || mrs->isCellNameFieldReading())
 152         {
 153           _names=DataArrayAsciiChar::New();
 154           _names->alloc(curNbOfElem+1,MED_SNAME_SIZE);//not a bug to avoid the memory corruption due to last \0 at the end
 155           if(MEDmeshEntityNameRd(fid,mName,dt,it,entity,geoElt,_names->getPointer())!=0)
 156             _names=0;
 157           else
 158             _names->reAlloc(curNbOfElem);//not a bug to avoid the memory corruption due to last \0 at the end
 159         }
 160     }
 161 }
 162
 163 void MEDFileUMeshPerType::loadPolyg(med_idt fid, const char *mName, int dt, int it, int mdim, int arraySize, med_geometry_type geoElt,
 164                                     med_entity_type entity, MEDFileMeshReadSelector *mrs)
 165 {
 166   med_bool changement,transformation;
 167   med_int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_INDEX_NODE,MED_NODAL,&changement,&transformation)-1;
 168   _m=MEDCoupling1DGTUMesh::New(mName,geoElt==MED_POLYGON?INTERP_KERNEL::NORM_POLYGON:INTERP_KERNEL::NORM_QPOLYG);
 169   MEDCoupling1DGTUMesh *mc(dynamic_cast<MEDCoupling1DGTUMesh *>((MEDCoupling1GTUMesh *)_m));
 170   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New()),connI(DataArrayInt::New());
 171   conn->alloc(arraySize,1); connI->alloc(curNbOfElem+1,1);
 172   MEDmeshPolygon2Rd(fid,mName,dt,it,MED_CELL,geoElt,MED_NODAL,connI->getPointer(),conn->getPointer());
 173   std::transform(conn->begin(),conn->end(),conn->getPointer(),std::bind2nd(std::plus<int>(),-1));
 174   std::transform(connI->begin(),connI->end(),connI->getPointer(),std::bind2nd(std::plus<int>(),-1));
 175   mc->setNodalConnectivity(conn,connI);
 176   loadCommonPart(fid,mName,dt,it,mdim,curNbOfElem,geoElt,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   _m=MEDCoupling1DGTUMesh::New(mName,INTERP_KERNEL::NORM_POLYHED);
 186   MEDCoupling1DGTUMesh *mc(dynamic_cast<MEDCoupling1DGTUMesh *>((MEDCoupling1GTUMesh *)_m));
 187   INTERP_KERNEL::AutoPtr<int> index=new int[curNbOfElem+1];
 188   INTERP_KERNEL::AutoPtr<int> indexFace=new int[indexFaceLgth];
 189   INTERP_KERNEL::AutoPtr<int> locConn=new int[connFaceLgth];
 190   MEDmeshPolyhedronRd(fid,mName,dt,it,MED_CELL,MED_NODAL,index,indexFace,locConn);
 191   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New()),connI(DataArrayInt::New());
 192   int arraySize=connFaceLgth;
 193   for(int i=0;i<curNbOfElem;i++)
 194     arraySize+=index[i+1]-index[i]-1;
 195   conn=DataArrayInt::New();
 196   conn->alloc(arraySize,1);
 197   int *wFinalConn=conn->getPointer();
 198   connI->alloc(curNbOfElem+1,1);
 199   int *finalIndex(connI->getPointer());
 200   finalIndex[0]=0;
 201   for(int i=0;i<curNbOfElem;i++)
 202     {
 203       finalIndex[i+1]=finalIndex[i]+index[i+1]-index[i]-1+indexFace[index[i+1]-1]-indexFace[index[i]-1];
 204       wFinalConn=std::transform(locConn+indexFace[index[i]-1]-1,locConn+indexFace[index[i]]-1,wFinalConn,std::bind2nd(std::plus<int>(),-1));
 205       for(int j=index[i];j<index[i+1]-1;j++)
 206         {
 207           *wFinalConn++=-1;
 208           wFinalConn=std::transform(locConn+indexFace[j]-1,locConn+indexFace[j+1]-1,wFinalConn,std::bind2nd(std::plus<int>(),-1));
 209         }
 210     }
 211   mc->setNodalConnectivity(conn,connI);
 212   loadCommonPart(fid,mName,dt,it,mdim,curNbOfElem,MED_POLYHEDRON,entity,mrs);
 213 }
 214
 215 void MEDFileUMeshPerType::Write(med_idt fid, const std::string& mname, int mdim, const MEDCoupling1GTUMesh *m, const DataArrayInt *fam, const DataArrayInt *num, const DataArrayAsciiChar *names)
 216 {
 217   int nbOfCells=m->getNumberOfCells();
 218   if(nbOfCells<1)
 219     return ;
 220   int dt,it;
 221   double timm=m->getTime(dt,it);
 222   INTERP_KERNEL::NormalizedCellType ikt=m->getTypeOfCell(0);
 223   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(ikt);
 224   med_geometry_type curMedType=typmai3[(int)ikt];
 225   if(!cm.isDynamic())
 226     {
 227       const MEDCoupling1SGTUMesh *m0(dynamic_cast<const MEDCoupling1SGTUMesh *>(m));
 228       if(!m0)
 229         throw INTERP_KERNEL::Exception("MEDFileUMeshPerType::Write : internal error #1 !");
 230       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr(m0->getNodalConnectivity()->deepCpy());
 231       std::transform(arr->begin(),arr->end(),arr->getPointer(),std::bind2nd(std::plus<int>(),1));
 232       MEDmeshElementConnectivityWr(fid,mname.c_str(),dt,it,timm,MED_CELL,curMedType,MED_NODAL,MED_FULL_INTERLACE,nbOfCells,arr->begin());
 233     }
 234   else
 235     {
 236       const MEDCoupling1DGTUMesh *m0(dynamic_cast<const MEDCoupling1DGTUMesh *>(m));
 237       if(!m0)
 238         throw INTERP_KERNEL::Exception("MEDFileUMeshPerType::Write : internal error #2 !");
 239       if(ikt==INTERP_KERNEL::NORM_POLYGON || ikt==INTERP_KERNEL::NORM_QPOLYG)
 240         {
 241           MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr(m0->getNodalConnectivity()->deepCpy()),arrI(m0->getNodalConnectivityIndex()->deepCpy());
 242           std::transform(arr->begin(),arr->end(),arr->getPointer(),std::bind2nd(std::plus<int>(),1));
 243           std::transform(arrI->begin(),arrI->end(),arrI->getPointer(),std::bind2nd(std::plus<int>(),1));
 244           MEDmeshPolygon2Wr(fid,mname.c_str(),dt,it,timm,MED_CELL,ikt==INTERP_KERNEL::NORM_POLYGON?MED_POLYGON:MED_POLYGON2,MED_NODAL,nbOfCells+1,arrI->begin(),arr->begin());
 245         }
 246       else
 247         {
 248           const int *conn(m0->getNodalConnectivity()->begin()),*connI(m0->getNodalConnectivityIndex()->begin());
 249           int meshLgth=m0->getNodalConnectivityLength();
 250           int nbOfFaces=std::count(conn,conn+meshLgth,-1)+nbOfCells;
 251           INTERP_KERNEL::AutoPtr<int> tab1=new int[nbOfCells+1];
 252           int *w1=tab1; *w1=1;
 253           INTERP_KERNEL::AutoPtr<int> tab2=new int[nbOfFaces+1];
 254           int *w2=tab2; *w2=1;
 255           INTERP_KERNEL::AutoPtr<int> bigtab=new int[meshLgth];
 256           int *bt=bigtab;
 257           for(int i=0;i<nbOfCells;i++,w1++)
 258             {
 259               int nbOfFaces2=0;
 260               for(const int *w=conn+connI[i];w!=conn+connI[i+1];w2++)
 261                 {
 262                   const int *wend=std::find(w,conn+connI[i+1],-1);
 263                   bt=std::transform(w,wend,bt,std::bind2nd(std::plus<int>(),1));
 264                   int nbOfNode=std::distance(w,wend);
 265                   w2[1]=w2[0]+nbOfNode;
 266                   if(wend!=conn+connI[i+1])
 267                     w=wend+1;
 268                   else
 269                     w=wend;
 270                   nbOfFaces2++;
 271                 }
 272               w1[1]=w1[0]+nbOfFaces2;
 273             }
 274           MEDmeshPolyhedronWr(fid,mname.c_str(),dt,it,timm,MED_CELL,MED_NODAL,nbOfCells+1,tab1,nbOfFaces+1,tab2,bigtab);
 275         }
 276     }
 277   if(fam)
 278     MEDmeshEntityFamilyNumberWr(fid,mname.c_str(),dt,it,MED_CELL,curMedType,nbOfCells,fam->getConstPointer());
 279   if(num)
 280     MEDmeshEntityNumberWr(fid,mname.c_str(),dt,it,MED_CELL,curMedType,nbOfCells,num->getConstPointer());
 281   if(names)
 282     {
 283       if(names->getNumberOfComponents()!=MED_SNAME_SIZE)
 284         {
 285           std::ostringstream oss; oss << "MEDFileUMeshPerType::write : expected a name field on cells with number of components set to " << MED_SNAME_SIZE;
 286           oss << " ! The array has " << names->getNumberOfComponents() << " components !";
 287           throw INTERP_KERNEL::Exception(oss.str().c_str());
 288         }
 289       MEDmeshEntityNameWr(fid,mname.c_str(),dt,it,MED_CELL,curMedType,nbOfCells,names->getConstPointer());
 290     }
 291 }