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