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
  23 #include "MEDCouplingUMesh.hxx"
  24
  25 #include "InterpKernelAutoPtr.hxx"
  26 #include "CellModel.hxx"
  27
  28 #include <iostream>
  29
  30 extern med_geometry_type typmai3[32];
  31
  32 using namespace ParaMEDMEM;
  33
  34 MEDFileUMeshPerType *MEDFileUMeshPerType::New(med_idt fid, const char *mName, int dt, int it, int mdim, med_geometry_type geoElt, INTERP_KERNEL::NormalizedCellType geoElt2)
  35 {
  36   med_entity_type whichEntity;
  37   if(!isExisting(fid,mName,dt,it,geoElt,whichEntity))
  38     return 0;
  39   return new MEDFileUMeshPerType(fid,mName,dt,it,mdim,geoElt,geoElt2,whichEntity);
  40 }
  41
  42 bool MEDFileUMeshPerType::isExisting(med_idt fid, const char *mName, int dt, int it, med_geometry_type geoElt, med_entity_type& whichEntity)
  43 {
  44   static const med_entity_type entities[3]={MED_CELL,MED_DESCENDING_FACE,MED_DESCENDING_EDGE};
  45   int nbOfElt=0;
  46   for(int i=0;i<3;i++)
  47     {
  48       med_bool changement,transformation;
  49       int tmp=MEDmeshnEntity(fid,mName,dt,it,entities[i],geoElt,MED_CONNECTIVITY,MED_NODAL,
  50                              &changement,&transformation);
  51       if(tmp>nbOfElt)
  52         {
  53           nbOfElt=tmp;
  54           whichEntity=entities[i];
  55           if(i>0)
  56             std::cerr << "WARNING : MEDFile has been detected to be no compilant with MED 3 : Please change entity in MEDFile for geotype " <<  geoElt << std::endl;
  57         }
  58     }
  59   return nbOfElt>0;
  60 }
  61
  62 int MEDFileUMeshPerType::getDim() const
  63 {
  64   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
  65   return cm.getDimension();
  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):_type(type),_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   if(type!=INTERP_KERNEL::NORM_POLYGON && type!=INTERP_KERNEL::NORM_POLYHED)
  75     {
  76       loadFromStaticType(fid,mName,dt,it,mdim,curNbOfElem,geoElt,type,entity);
  77       return;
  78     }
  79   if(type==INTERP_KERNEL::NORM_POLYGON)
  80     {
  81       loadPolyg(fid,mName,dt,it,mdim,curNbOfElem,geoElt,entity);
  82       return;
  83     }
  84   //if(type==INTERP_KERNEL::NORM_POLYHED)
  85   loadPolyh(fid,mName,dt,it,mdim,curNbOfElem,geoElt,entity);
  86 }
  87
  88 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,
  89                                              med_entity_type entity)
  90 {
  91   _conn=DataArrayInt::New();
  92   int nbOfNodesPerCell=(geoElt%100);
  93   _conn->alloc((nbOfNodesPerCell+1)*curNbOfElem,1);
  94   _conn_index=DataArrayInt::New();
  95   _conn_index->alloc(curNbOfElem+1,1);
  96   INTERP_KERNEL::AutoPtr<int> connTab=new int[(nbOfNodesPerCell)*curNbOfElem];
  97   MEDmeshElementConnectivityRd(fid,mName,dt,it,entity,geoElt,MED_NODAL,MED_FULL_INTERLACE,connTab);
  98   loadCommonPart(fid,mName,dt,it,mdim,curNbOfElem,geoElt,entity);
  99   int *w1=_conn->getPointer();
 100   int *w2=_conn_index->getPointer();
 101   *w2++=0;
 102   const int *wi=connTab;
 103   for(int i=0;i<curNbOfElem;i++,wi+=nbOfNodesPerCell,w2++)
 104     {
 105       *w1++=(int)type;
 106       w1=std::transform(wi,wi+nbOfNodesPerCell,w1,std::bind2nd(std::plus<int>(),-1));
 107       *w2=w2[-1]+nbOfNodesPerCell+1;
 108     }
 109 }
 110
 111 void MEDFileUMeshPerType::loadCommonPart(med_idt fid, const char *mName, int dt, int it, int mdim, int curNbOfElem, med_geometry_type geoElt,
 112                                          med_entity_type entity)
 113 {
 114   med_bool changement,transformation;
 115   _fam=DataArrayInt::New();
 116   _fam->alloc(curNbOfElem,1);
 117   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_FAMILY_NUMBER,MED_NODAL,&changement,&transformation)>0)
 118     {
 119       if(MEDmeshEntityFamilyNumberRd(fid,mName,dt,it,entity,geoElt,_fam->getPointer())!=0)
 120         std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
 121     }
 122   else
 123     std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
 124   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_NUMBER,MED_NODAL,&changement,&transformation)>0)
 125     {
 126       _num=DataArrayInt::New();
 127       _num->alloc(curNbOfElem,1);
 128       if(MEDmeshEntityNumberRd(fid,mName,dt,it,entity,geoElt,_num->getPointer())!=0)
 129         _num=0;
 130     }
 131   else
 132     _num=0;
 133   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_NAME,MED_NODAL,&changement,&transformation)>0)
 134     {
 135       _names=DataArrayAsciiChar::New();
 136       _names->alloc(curNbOfElem+1,MED_SNAME_SIZE);//not a bug to avoid the memory corruption due to last \0 at the end
 137       if(MEDmeshEntityNameRd(fid,mName,dt,it,entity,geoElt,_names->getPointer())!=0)
 138         _names=0;
 139       else
 140         _names->reAlloc(curNbOfElem);//not a bug to avoid the memory corruption due to last \0 at the end
 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)
 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   _conn_index=DataArrayInt::New();
 150   _conn_index->alloc(curNbOfElem+1,1);
 151   _conn=DataArrayInt::New();
 152   _conn->alloc(arraySize+curNbOfElem,1);
 153   INTERP_KERNEL::AutoPtr<int> locConn=new int[arraySize];
 154   MEDmeshPolygonRd(fid,mName,dt,it,MED_CELL,MED_NODAL,_conn_index->getPointer(),locConn);
 155   int *w1=_conn->getPointer();
 156   int *w2=_conn_index->getPointer();
 157   const int *wi=locConn;
 158   for(int i=0;i<curNbOfElem;i++,w2++)
 159     {
 160       *w1++=(int)INTERP_KERNEL::NORM_POLYGON;
 161       const int *wi2=wi+(w2[1]-w2[0]);
 162       w1=std::transform(wi,wi2,w1,std::bind2nd(std::plus<int>(),-1));
 163       wi=wi2;
 164       *w2=*w2-1+i;
 165     }
 166   *w2=*w2-1+curNbOfElem;
 167   loadCommonPart(fid,mName,dt,it,mdim,curNbOfElem,MED_POLYGON,entity);
 168 }
 169
 170 void MEDFileUMeshPerType::loadPolyh(med_idt fid, const char *mName, int dt, int it, int mdim, int connFaceLgth, med_geometry_type geoElt,
 171                                     med_entity_type entity)
 172 {
 173   med_bool changement,transformation;
 174   med_int indexFaceLgth=MEDmeshnEntity(fid,mName,dt,it,MED_CELL,MED_POLYHEDRON,MED_INDEX_NODE,MED_NODAL,&changement,&transformation);
 175   int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,MED_CELL,MED_POLYHEDRON,MED_INDEX_FACE,MED_NODAL,&changement,&transformation)-1;
 176   INTERP_KERNEL::AutoPtr<int> index=new int[curNbOfElem+1];
 177   INTERP_KERNEL::AutoPtr<int> indexFace=new int[indexFaceLgth];
 178   INTERP_KERNEL::AutoPtr<int> locConn=new int[connFaceLgth];
 179   MEDmeshPolyhedronRd(fid,mName,dt,it,MED_CELL,MED_NODAL,index,indexFace,locConn);
 180   int arraySize=connFaceLgth;
 181   for(int i=0;i<curNbOfElem;i++)
 182     arraySize+=index[i+1]-index[i]-1;
 183   _conn=DataArrayInt::New();
 184   _conn->alloc(arraySize+curNbOfElem,1);
 185   int *wFinalConn=_conn->getPointer();
 186   _conn_index=DataArrayInt::New();
 187   _conn_index->alloc(curNbOfElem+1,1);
 188   int *finalIndex=_conn_index->getPointer();
 189   finalIndex[0]=0;
 190   for(int i=0;i<curNbOfElem;i++)
 191     {
 192       *wFinalConn++=(int)INTERP_KERNEL::NORM_POLYHED;
 193       finalIndex[i+1]=finalIndex[i]+index[i+1]-index[i]-1+indexFace[index[i+1]-1]-indexFace[index[i]-1]+1;
 194       wFinalConn=std::transform(locConn+indexFace[index[i]-1]-1,locConn+indexFace[index[i]]-1,wFinalConn,std::bind2nd(std::plus<int>(),-1));
 195       for(int j=index[i];j<index[i+1]-1;j++)
 196         {
 197           *wFinalConn++=-1;
 198           wFinalConn=std::transform(locConn+indexFace[j]-1,locConn+indexFace[j+1]-1,wFinalConn,std::bind2nd(std::plus<int>(),-1));
 199         }
 200     }
 201   loadCommonPart(fid,mName,dt,it,mdim,curNbOfElem,MED_POLYHEDRON,entity);
 202 }
 203
 204 void MEDFileUMeshPerType::write(med_idt fid, const char *mname, int mdim, const MEDCouplingUMesh *m, const DataArrayInt *fam, const DataArrayInt *num, const DataArrayAsciiChar *names)
 205 {
 206   int nbOfCells=m->getNumberOfCells();
 207   if(nbOfCells<1)
 208     return ;
 209   int dt,it;
 210   double timm=m->getTime(dt,it);
 211   INTERP_KERNEL::NormalizedCellType ikt=m->getTypeOfCell(0);
 212   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(ikt);
 213   med_geometry_type curMedType=typmai3[(int)ikt];
 214   const int *conn=m->getNodalConnectivity()->getConstPointer();
 215   const int *connI=m->getNodalConnectivityIndex()->getConstPointer();
 216   if(ikt!=INTERP_KERNEL::NORM_POLYGON && ikt!=INTERP_KERNEL::NORM_POLYHED)
 217     {
 218       int nbNodesPerCell=cm.getNumberOfNodes();
 219       INTERP_KERNEL::AutoPtr<int> tab=new int[nbNodesPerCell*nbOfCells];
 220       int *w=tab;
 221       for(int i=0;i<nbOfCells;i++)
 222         w=std::transform(conn+connI[i]+1,conn+connI[i+1],w,std::bind2nd(std::plus<int>(),1));
 223       MEDmeshElementConnectivityWr(fid,mname,dt,it,timm,MED_CELL,curMedType,MED_NODAL,MED_FULL_INTERLACE,nbOfCells,tab);
 224     }
 225   else
 226     {
 227       if(ikt==INTERP_KERNEL::NORM_POLYGON)
 228         {
 229           INTERP_KERNEL::AutoPtr<int> tab1=new int[nbOfCells+1];
 230           INTERP_KERNEL::AutoPtr<int> tab2=new int[m->getMeshLength()];
 231           int *wI=tab1; *wI=1;
 232           int *w=tab2;
 233           for(int i=0;i<nbOfCells;i++,wI++)
 234             {
 235               wI[1]=wI[0]+connI[i+1]-connI[i]-1;
 236               w=std::transform(conn+connI[i]+1,conn+connI[i+1],w,std::bind2nd(std::plus<int>(),1));
 237             }
 238           MEDmeshPolygonWr(fid,mname,dt,it,timm,MED_CELL,MED_NODAL,nbOfCells+1,tab1,tab2);
 239         }
 240       else
 241         {
 242           int meshLgth=m->getMeshLength();
 243           int nbOfFaces=std::count(conn,conn+meshLgth,-1)+nbOfCells;
 244           INTERP_KERNEL::AutoPtr<int> tab1=new int[nbOfCells+1];
 245           int *w1=tab1; *w1=1;
 246           INTERP_KERNEL::AutoPtr<int> tab2=new int[nbOfFaces+1];
 247           int *w2=tab2; *w2=1;
 248           INTERP_KERNEL::AutoPtr<int> bigtab=new int[meshLgth-nbOfCells];
 249           int *bt=bigtab;
 250           for(int i=0;i<nbOfCells;i++,w1++)
 251             {
 252               int nbOfFaces2=0;
 253               for(const int *w=conn+connI[i]+1;w!=conn+connI[i+1];w2++)
 254                 {
 255                   const int *wend=std::find(w,conn+connI[i+1],-1);
 256                   bt=std::transform(w,wend,bt,std::bind2nd(std::plus<int>(),1));
 257                   int nbOfNode=std::distance(w,wend);
 258                   w2[1]=w2[0]+nbOfNode;
 259                   if(wend!=conn+connI[i+1])
 260                     w=wend+1;
 261                   else
 262                     w=wend;
 263                   nbOfFaces2++;
 264                 }
 265               w1[1]=w1[0]+nbOfFaces2;
 266             }
 267           MEDmeshPolyhedronWr(fid,mname,dt,it,timm,MED_CELL,MED_NODAL,nbOfCells+1,tab1,nbOfFaces+1,tab2,bigtab);
 268         }
 269     }
 270   if(fam)
 271     MEDmeshEntityFamilyNumberWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,fam->getConstPointer());
 272   if(num)
 273     MEDmeshEntityNumberWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,num->getConstPointer());
 274   if(names)
 275     {
 276       if(names->getNumberOfComponents()!=MED_SNAME_SIZE)
 277         {
 278           std::ostringstream oss; oss << "MEDFileUMeshPerType::write : expected a name field on cells with number of components set to " << MED_SNAME_SIZE;
 279           oss << " ! The array has " << names->getNumberOfComponents() << " components !";
 280           throw INTERP_KERNEL::Exception(oss.str().c_str());
 281         }
 282       MEDmeshEntityNameWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,names->getConstPointer());
 283     }
 284 }