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   return _m->getMeshDimension();
  66 }
  67
  68 MEDFileUMeshPerType::MEDFileUMeshPerType(med_idt fid, const char *mName, int dt, int it, int mdim, med_geometry_type geoElt, INTERP_KERNEL::NormalizedCellType type,
  69                                          med_entity_type entity, MEDFileMeshReadSelector *mrs):_entity(entity)
  70 {
  71   med_bool changement,transformation;
  72   int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_CONNECTIVITY,MED_NODAL,
  73                                  &changement,&transformation);
  74   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
  75   if(!cm.isDynamic())
  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   _m=MEDCoupling1SGTUMesh::New(mName,type);
  93   MEDCoupling1SGTUMesh *mc(dynamic_cast<MEDCoupling1SGTUMesh *>((MEDCoupling1GTUMesh *)_m));
  94   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
  95   int nbOfNodesPerCell=mc->getNumberOfNodesPerCell();
  96   conn->alloc(nbOfNodesPerCell*curNbOfElem,1);
  97   MEDmeshElementConnectivityRd(fid,mName,dt,it,entity,geoElt,MED_NODAL,MED_FULL_INTERLACE,conn->getPointer());
  98   std::transform(conn->begin(),conn->end(),conn->getPointer(),std::bind2nd(std::plus<int>(),-1));
  99   mc->setNodalConnectivity(conn);
 100   loadCommonPart(fid,mName,dt,it,mdim,curNbOfElem,geoElt,entity,mrs);
 101 }
 102
 103 void MEDFileUMeshPerType::loadCommonPart(med_idt fid, const char *mName, int dt, int it, int mdim, int curNbOfElem, med_geometry_type geoElt,
 104                                          med_entity_type entity, MEDFileMeshReadSelector *mrs)
 105 {
 106   med_bool changement,transformation;
 107   _fam=0;
 108   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_FAMILY_NUMBER,MED_NODAL,&changement,&transformation)>0)
 109     {
 110       if(!mrs || mrs->isCellFamilyFieldReading())
 111         {
 112           _fam=DataArrayInt::New();
 113           _fam->alloc(curNbOfElem,1);
 114           if(MEDmeshEntityFamilyNumberRd(fid,mName,dt,it,entity,geoElt,_fam->getPointer())!=0)
 115             std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
 116         }
 117     }
 118   _num=0;
 119   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_NUMBER,MED_NODAL,&changement,&transformation)>0)
 120     {
 121       if(!mrs || mrs->isCellNumFieldReading())
 122         {
 123           _num=DataArrayInt::New();
 124           _num->alloc(curNbOfElem,1);
 125           if(MEDmeshEntityNumberRd(fid,mName,dt,it,entity,geoElt,_num->getPointer())!=0)
 126             _num=0;
 127         }
 128     }
 129   _names=0;
 130   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_NAME,MED_NODAL,&changement,&transformation)>0)
 131     {
 132       if(!mrs || mrs->isCellNameFieldReading())
 133         {
 134           _names=DataArrayAsciiChar::New();
 135           _names->alloc(curNbOfElem+1,MED_SNAME_SIZE);//not a bug to avoid the memory corruption due to last \0 at the end
 136           if(MEDmeshEntityNameRd(fid,mName,dt,it,entity,geoElt,_names->getPointer())!=0)
 137             _names=0;
 138           else
 139             _names->reAlloc(curNbOfElem);//not a bug to avoid the memory corruption due to last \0 at the end
 140         }
 141     }
 142 }
 143
 144 void MEDFileUMeshPerType::loadPolyg(med_idt fid, const char *mName, int dt, int it, int mdim, int arraySize, med_geometry_type geoElt,
 145                                     med_entity_type entity, MEDFileMeshReadSelector *mrs)
 146 {
 147   med_bool changement,transformation;
 148   med_int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,entity,MED_POLYGON,MED_INDEX_NODE,MED_NODAL,&changement,&transformation)-1;
 149   _m=MEDCoupling1DGTUMesh::New(mName,INTERP_KERNEL::NORM_POLYGON);
 150   MEDCoupling1DGTUMesh *mc(dynamic_cast<MEDCoupling1DGTUMesh *>((MEDCoupling1GTUMesh *)_m));
 151   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New()),connI(DataArrayInt::New());
 152   conn->alloc(arraySize,1); connI->alloc(curNbOfElem+1,1);
 153   MEDmeshPolygonRd(fid,mName,dt,it,MED_CELL,MED_NODAL,connI->getPointer(),conn->getPointer());
 154   std::transform(conn->begin(),conn->end(),conn->getPointer(),std::bind2nd(std::plus<int>(),-1));
 155   std::transform(connI->begin(),connI->end(),connI->getPointer(),std::bind2nd(std::plus<int>(),-1));
 156   mc->setNodalConnectivity(conn,connI);
 157   loadCommonPart(fid,mName,dt,it,mdim,curNbOfElem,MED_POLYGON,entity,mrs);
 158 }
 159
 160 void MEDFileUMeshPerType::loadPolyh(med_idt fid, const char *mName, int dt, int it, int mdim, int connFaceLgth, med_geometry_type geoElt,
 161                                     med_entity_type entity, MEDFileMeshReadSelector *mrs)
 162 {
 163   med_bool changement,transformation;
 164   med_int indexFaceLgth=MEDmeshnEntity(fid,mName,dt,it,MED_CELL,MED_POLYHEDRON,MED_INDEX_NODE,MED_NODAL,&changement,&transformation);
 165   int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,MED_CELL,MED_POLYHEDRON,MED_INDEX_FACE,MED_NODAL,&changement,&transformation)-1;
 166   _m=MEDCoupling1DGTUMesh::New(mName,INTERP_KERNEL::NORM_POLYHED);
 167   MEDCoupling1DGTUMesh *mc(dynamic_cast<MEDCoupling1DGTUMesh *>((MEDCoupling1GTUMesh *)_m));
 168   INTERP_KERNEL::AutoPtr<int> index=new int[curNbOfElem+1];
 169   INTERP_KERNEL::AutoPtr<int> indexFace=new int[indexFaceLgth];
 170   INTERP_KERNEL::AutoPtr<int> locConn=new int[connFaceLgth];
 171   MEDmeshPolyhedronRd(fid,mName,dt,it,MED_CELL,MED_NODAL,index,indexFace,locConn);
 172   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New()),connI(DataArrayInt::New());
 173   int arraySize=connFaceLgth;
 174   for(int i=0;i<curNbOfElem;i++)
 175     arraySize+=index[i+1]-index[i]-1;
 176   conn=DataArrayInt::New();
 177   conn->alloc(arraySize,1);
 178   int *wFinalConn=conn->getPointer();
 179   connI->alloc(curNbOfElem+1,1);
 180   int *finalIndex(connI->getPointer());
 181   finalIndex[0]=0;
 182   for(int i=0;i<curNbOfElem;i++)
 183     {
 184       finalIndex[i+1]=finalIndex[i]+index[i+1]-index[i]-1+indexFace[index[i+1]-1]-indexFace[index[i]-1];
 185       wFinalConn=std::transform(locConn+indexFace[index[i]-1]-1,locConn+indexFace[index[i]]-1,wFinalConn,std::bind2nd(std::plus<int>(),-1));
 186       for(int j=index[i];j<index[i+1]-1;j++)
 187         {
 188           *wFinalConn++=-1;
 189           wFinalConn=std::transform(locConn+indexFace[j]-1,locConn+indexFace[j+1]-1,wFinalConn,std::bind2nd(std::plus<int>(),-1));
 190         }
 191     }
 192   mc->setNodalConnectivity(conn,connI);
 193   loadCommonPart(fid,mName,dt,it,mdim,curNbOfElem,MED_POLYHEDRON,entity,mrs);
 194 }
 195
 196 void MEDFileUMeshPerType::Write(med_idt fid, const char *mname, int mdim, const MEDCoupling1GTUMesh *m, const DataArrayInt *fam, const DataArrayInt *num, const DataArrayAsciiChar *names)
 197 {
 198   int nbOfCells=m->getNumberOfCells();
 199   if(nbOfCells<1)
 200     return ;
 201   int dt,it;
 202   double timm=m->getTime(dt,it);
 203   INTERP_KERNEL::NormalizedCellType ikt=m->getTypeOfCell(0);
 204   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(ikt);
 205   med_geometry_type curMedType=typmai3[(int)ikt];
 206   if(!cm.isDynamic())
 207     {
 208       const MEDCoupling1SGTUMesh *m0(dynamic_cast<const MEDCoupling1SGTUMesh *>(m));
 209       if(!m0)
 210         throw INTERP_KERNEL::Exception("MEDFileUMeshPerType::Write : internal error #1 !");
 211       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr(m0->getNodalConnectivity()->deepCpy());
 212       std::transform(arr->begin(),arr->end(),arr->getPointer(),std::bind2nd(std::plus<int>(),1));
 213       MEDmeshElementConnectivityWr(fid,mname,dt,it,timm,MED_CELL,curMedType,MED_NODAL,MED_FULL_INTERLACE,nbOfCells,arr->begin());
 214     }
 215   else
 216     {
 217       const MEDCoupling1DGTUMesh *m0(dynamic_cast<const MEDCoupling1DGTUMesh *>(m));
 218       if(!m0)
 219         throw INTERP_KERNEL::Exception("MEDFileUMeshPerType::Write : internal error #2 !");
 220       if(ikt==INTERP_KERNEL::NORM_POLYGON)
 221         {
 222           MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr(m0->getNodalConnectivity()->deepCpy()),arrI(m0->getNodalConnectivityIndex()->deepCpy());
 223           std::transform(arr->begin(),arr->end(),arr->getPointer(),std::bind2nd(std::plus<int>(),1));
 224           std::transform(arrI->begin(),arrI->end(),arrI->getPointer(),std::bind2nd(std::plus<int>(),1));
 225           MEDmeshPolygonWr(fid,mname,dt,it,timm,MED_CELL,MED_NODAL,nbOfCells+1,arrI->begin(),arr->begin());
 226         }
 227       else
 228         {
 229           const int *conn(m0->getNodalConnectivity()->begin()),*connI(m0->getNodalConnectivityIndex()->begin());
 230           int meshLgth=m0->getNodalConnectivityLength();
 231           int nbOfFaces=std::count(conn,conn+meshLgth,-1)+nbOfCells;
 232           INTERP_KERNEL::AutoPtr<int> tab1=new int[nbOfCells+1];
 233           int *w1=tab1; *w1=1;
 234           INTERP_KERNEL::AutoPtr<int> tab2=new int[nbOfFaces+1];
 235           int *w2=tab2; *w2=1;
 236           INTERP_KERNEL::AutoPtr<int> bigtab=new int[meshLgth];
 237           int *bt=bigtab;
 238           for(int i=0;i<nbOfCells;i++,w1++)
 239             {
 240               int nbOfFaces2=0;
 241               for(const int *w=conn+connI[i];w!=conn+connI[i+1];w2++)
 242                 {
 243                   const int *wend=std::find(w,conn+connI[i+1],-1);
 244                   bt=std::transform(w,wend,bt,std::bind2nd(std::plus<int>(),1));
 245                   int nbOfNode=std::distance(w,wend);
 246                   w2[1]=w2[0]+nbOfNode;
 247                   if(wend!=conn+connI[i+1])
 248                     w=wend+1;
 249                   else
 250                     w=wend;
 251                   nbOfFaces2++;
 252                 }
 253               w1[1]=w1[0]+nbOfFaces2;
 254             }
 255           MEDmeshPolyhedronWr(fid,mname,dt,it,timm,MED_CELL,MED_NODAL,nbOfCells+1,tab1,nbOfFaces+1,tab2,bigtab);
 256         }
 257     }
 258   if(fam)
 259     MEDmeshEntityFamilyNumberWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,fam->getConstPointer());
 260   if(num)
 261     MEDmeshEntityNumberWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,num->getConstPointer());
 262   if(names)
 263     {
 264       if(names->getNumberOfComponents()!=MED_SNAME_SIZE)
 265         {
 266           std::ostringstream oss; oss << "MEDFileUMeshPerType::write : expected a name field on cells with number of components set to " << MED_SNAME_SIZE;
 267           oss << " ! The array has " << names->getNumberOfComponents() << " components !";
 268           throw INTERP_KERNEL::Exception(oss.str().c_str());
 269         }
 270       MEDmeshEntityNameWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,names->getConstPointer());
 271     }
 272 }