src/MEDLoader/MEDFileMeshElt.cxx

   1 // Copyright (C) 2007-2012  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   _num=DataArrayInt::New();
  98   _num->alloc(curNbOfElem,1);
  99   _fam=DataArrayInt::New();
 100   _fam->alloc(curNbOfElem,1);
 101   med_bool changement,transformation;
 102   INTERP_KERNEL::AutoPtr<char> noms=new char[MED_SNAME_SIZE*(std::size_t)curNbOfElem+1];
 103   MEDmeshElementConnectivityRd(fid,mName,dt,it,entity,geoElt,MED_NODAL,MED_FULL_INTERLACE,connTab);
 104   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_FAMILY_NUMBER,MED_NODAL,&changement,&transformation)>0)
 105     {
 106       if(MEDmeshEntityFamilyNumberRd(fid,mName,dt,it,entity,geoElt,_fam->getPointer())!=0)
 107         std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
 108     }
 109   else
 110     std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
 111   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_NUMBER,MED_NODAL,&changement,&transformation)>0)
 112     {
 113       _num=DataArrayInt::New();
 114       _num->alloc(curNbOfElem,1);
 115       if(MEDmeshEntityNumberRd(fid,mName,dt,it,entity,geoElt,_num->getPointer())!=0)
 116         _num=0;
 117     }
 118   else
 119     _num=0;
 120   int *w1=_conn->getPointer();
 121   int *w2=_conn_index->getPointer();
 122   *w2++=0;
 123   const int *wi=connTab;
 124   for(int i=0;i<curNbOfElem;i++,wi+=nbOfNodesPerCell,w2++)
 125     {
 126       *w1++=(int)type;
 127       w1=std::transform(wi,wi+nbOfNodesPerCell,w1,std::bind2nd(std::plus<int>(),-1));
 128       *w2=w2[-1]+nbOfNodesPerCell+1;
 129     }
 130 }
 131
 132 void MEDFileUMeshPerType::loadPolyg(med_idt fid, const char *mName, int dt, int it, int mdim, int arraySize, med_geometry_type geoElt,
 133                                     med_entity_type entity)
 134 {
 135   med_bool changement,transformation;
 136   med_int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,entity,MED_POLYGON,MED_INDEX_NODE,MED_NODAL,&changement,&transformation)-1;
 137   _conn_index=DataArrayInt::New();
 138   _conn_index->alloc(curNbOfElem+1,1);
 139   _conn=DataArrayInt::New();
 140   _conn->alloc(arraySize+curNbOfElem,1);
 141   _num=DataArrayInt::New();
 142   _num->alloc(curNbOfElem,1);
 143   _fam=DataArrayInt::New();
 144   _fam->alloc(curNbOfElem,1);
 145   INTERP_KERNEL::AutoPtr<int> locConn=new int[arraySize];
 146   MEDmeshPolygonRd(fid,mName,dt,it,MED_CELL,MED_NODAL,_conn_index->getPointer(),locConn);
 147   int *w1=_conn->getPointer();
 148   int *w2=_conn_index->getPointer();
 149   const int *wi=locConn;
 150   for(int i=0;i<curNbOfElem;i++,w2++)
 151     {
 152       *w1++=(int)INTERP_KERNEL::NORM_POLYGON;
 153       const int *wi2=wi+(w2[1]-w2[0]);
 154       w1=std::transform(wi,wi2,w1,std::bind2nd(std::plus<int>(),-1));
 155       wi=wi2;
 156       *w2=*w2-1+i;
 157     }
 158   *w2=*w2-1+curNbOfElem;
 159   if(MEDmeshnEntity(fid,mName,dt,it,entity,MED_POLYGON,MED_FAMILY_NUMBER,MED_NODAL,&changement,&transformation)>0)
 160     {
 161       if(MEDmeshEntityFamilyNumberRd(fid,mName,dt,it,entity,MED_POLYGON,_fam->getPointer())!=0)
 162         std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
 163     }
 164   else
 165     std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
 166   if(MEDmeshnEntity(fid,mName,dt,it,entity,MED_POLYGON,MED_NUMBER,MED_NODAL,&changement,&transformation)>0)
 167     {
 168       if(MEDmeshEntityNumberRd(fid,mName,dt,it,entity,MED_POLYGON,_num->getPointer())!=0)
 169         _num=0;
 170     }
 171   else
 172     _num=0;
 173 }
 174
 175 void MEDFileUMeshPerType::loadPolyh(med_idt fid, const char *mName, int dt, int it, int mdim, int connFaceLgth, med_geometry_type geoElt,
 176                                     med_entity_type entity)
 177 {
 178   med_bool changement,transformation;
 179   med_int indexFaceLgth=MEDmeshnEntity(fid,mName,dt,it,MED_CELL,MED_POLYHEDRON,MED_INDEX_NODE,MED_NODAL,&changement,&transformation);
 180   int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,MED_CELL,MED_POLYHEDRON,MED_INDEX_FACE,MED_NODAL,&changement,&transformation)-1;
 181   INTERP_KERNEL::AutoPtr<int> index=new int[curNbOfElem+1];
 182   INTERP_KERNEL::AutoPtr<int> indexFace=new int[indexFaceLgth];
 183   INTERP_KERNEL::AutoPtr<int> locConn=new int[connFaceLgth];
 184   _fam=DataArrayInt::New();
 185   _fam->alloc(curNbOfElem,1);
 186   MEDmeshPolyhedronRd(fid,mName,dt,it,MED_CELL,MED_NODAL,index,indexFace,locConn);
 187   if(MEDmeshnEntity(fid,mName,dt,it,entity,MED_POLYHEDRON,MED_FAMILY_NUMBER,MED_NODAL,&changement,&transformation)>0)
 188     {
 189       if(MEDmeshEntityFamilyNumberRd(fid,mName,dt,it,entity,MED_POLYHEDRON,_fam->getPointer())!=0)
 190         std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
 191     }
 192   else
 193     std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
 194   int arraySize=connFaceLgth;
 195   for(int i=0;i<curNbOfElem;i++)
 196     arraySize+=index[i+1]-index[i]-1;
 197   _conn=DataArrayInt::New();
 198   _conn->alloc(arraySize+curNbOfElem,1);
 199   int *wFinalConn=_conn->getPointer();
 200   _conn_index=DataArrayInt::New();
 201   _conn_index->alloc(curNbOfElem+1,1);
 202   int *finalIndex=_conn_index->getPointer();
 203   finalIndex[0]=0;
 204   for(int i=0;i<curNbOfElem;i++)
 205     {
 206       *wFinalConn++=(int)INTERP_KERNEL::NORM_POLYHED;
 207       finalIndex[i+1]=finalIndex[i]+index[i+1]-index[i]-1+indexFace[index[i+1]-1]-indexFace[index[i]-1]+1;
 208       wFinalConn=std::transform(locConn+indexFace[index[i]-1]-1,locConn+indexFace[index[i]]-1,wFinalConn,std::bind2nd(std::plus<int>(),-1));
 209       for(int j=index[i];j<index[i+1]-1;j++)
 210         {
 211           *wFinalConn++=-1;
 212           wFinalConn=std::transform(locConn+indexFace[j]-1,locConn+indexFace[j+1]-1,wFinalConn,std::bind2nd(std::plus<int>(),-1));
 213         }
 214     }
 215   _num=DataArrayInt::New();
 216   _num->alloc(curNbOfElem,1);
 217   if(MEDmeshnEntity(fid,mName,dt,it,entity,MED_POLYHEDRON,MED_NUMBER,MED_NODAL,&changement,&transformation)>0)
 218     {
 219       if(MEDmeshEntityNumberRd(fid,mName,dt,it,entity,MED_POLYHEDRON,_num->getPointer())!=0)
 220         _num=0;
 221     }
 222   else
 223     _num=0;
 224 }
 225
 226 void MEDFileUMeshPerType::write(med_idt fid, const char *mname, int mdim, const MEDCouplingUMesh *m, const DataArrayInt *fam, const DataArrayInt *num)
 227 {
 228   int nbOfCells=m->getNumberOfCells();
 229   if(nbOfCells<1)
 230     return ;
 231   int dt,it;
 232   double timm=m->getTime(dt,it);
 233   INTERP_KERNEL::NormalizedCellType ikt=m->getTypeOfCell(0);
 234   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(ikt);
 235   med_geometry_type curMedType=typmai3[(int)ikt];
 236   const int *conn=m->getNodalConnectivity()->getConstPointer();
 237   const int *connI=m->getNodalConnectivityIndex()->getConstPointer();
 238   if(ikt!=INTERP_KERNEL::NORM_POLYGON && ikt!=INTERP_KERNEL::NORM_POLYHED)
 239     {
 240       int nbNodesPerCell=cm.getNumberOfNodes();
 241       INTERP_KERNEL::AutoPtr<int> tab=new int[nbNodesPerCell*nbOfCells];
 242       int *w=tab;
 243       for(int i=0;i<nbOfCells;i++)
 244         w=std::transform(conn+connI[i]+1,conn+connI[i+1],w,std::bind2nd(std::plus<int>(),1));
 245       MEDmeshElementConnectivityWr(fid,mname,dt,it,timm,MED_CELL,curMedType,MED_NODAL,MED_FULL_INTERLACE,nbOfCells,tab);
 246     }
 247   else
 248     {
 249       if(ikt==INTERP_KERNEL::NORM_POLYGON)
 250         {
 251           INTERP_KERNEL::AutoPtr<int> tab1=new int[nbOfCells+1];
 252           INTERP_KERNEL::AutoPtr<int> tab2=new int[m->getMeshLength()];
 253           int *wI=tab1; *wI=1;
 254           int *w=tab2;
 255           for(int i=0;i<nbOfCells;i++,wI++)
 256             {
 257               wI[1]=wI[0]+connI[i+1]-connI[i]-1;
 258               w=std::transform(conn+connI[i]+1,conn+connI[i+1],w,std::bind2nd(std::plus<int>(),1));
 259             }
 260           MEDmeshPolygonWr(fid,mname,dt,it,timm,MED_CELL,MED_NODAL,nbOfCells+1,tab1,tab2);
 261         }
 262       else
 263         {
 264           int meshLgth=m->getMeshLength();
 265           int nbOfFaces=std::count(conn,conn+meshLgth,-1)+nbOfCells;
 266           INTERP_KERNEL::AutoPtr<int> tab1=new int[nbOfCells+1];
 267           int *w1=tab1; *w1=1;
 268           INTERP_KERNEL::AutoPtr<int> tab2=new int[nbOfFaces+1];
 269           int *w2=tab2; *w2=1;
 270           INTERP_KERNEL::AutoPtr<int> bigtab=new int[meshLgth-nbOfCells];
 271           int *bt=bigtab;
 272           for(int i=0;i<nbOfCells;i++,w1++)
 273             {
 274               int nbOfFaces2=0;
 275               for(const int *w=conn+connI[i]+1;w!=conn+connI[i+1];w2++)
 276                 {
 277                   const int *wend=std::find(w,conn+connI[i+1],-1);
 278                   bt=std::transform(w,wend,bt,std::bind2nd(std::plus<int>(),1));
 279                   int nbOfNode=std::distance(w,wend);
 280                   w2[1]=w2[0]+nbOfNode;
 281                   if(wend!=conn+connI[i+1])
 282                     w=wend+1;
 283                   else
 284                     w=wend;
 285                   nbOfFaces2++;
 286                 }
 287               w1[1]=w1[0]+nbOfFaces2;
 288             }
 289           MEDmeshPolyhedronWr(fid,mname,dt,it,timm,MED_CELL,MED_NODAL,nbOfCells+1,tab1,nbOfFaces+1,tab2,bigtab);
 290         }
 291     }
 292   if(fam)
 293     MEDmeshEntityFamilyNumberWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,fam->getConstPointer());
 294   if(num)
 295     MEDmeshEntityNumberWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,num->getConstPointer());
 296 }