Back to index

salome-med  6.5.0
MEDFileMeshElt.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "MEDFileMeshElt.hxx"
00021 
00022 #include "MEDCouplingUMesh.hxx"
00023 
00024 #include "InterpKernelAutoPtr.hxx"
00025 #include "CellModel.hxx"
00026 
00027 #include <iostream>
00028 
00029 extern med_geometry_type typmai3[32];
00030 
00031 using namespace ParaMEDMEM;
00032 
00033 MEDFileUMeshPerType *MEDFileUMeshPerType::New(med_idt fid, const char *mName, int dt, int it, int mdim, med_geometry_type geoElt, INTERP_KERNEL::NormalizedCellType geoElt2)
00034 {
00035   med_entity_type whichEntity;
00036   if(!isExisting(fid,mName,dt,it,geoElt,whichEntity))
00037     return 0;
00038   return new MEDFileUMeshPerType(fid,mName,dt,it,mdim,geoElt,geoElt2,whichEntity);
00039 }
00040 
00041 bool MEDFileUMeshPerType::isExisting(med_idt fid, const char *mName, int dt, int it, med_geometry_type geoElt, med_entity_type& whichEntity)
00042 {
00043   static const med_entity_type entities[3]={MED_CELL,MED_DESCENDING_FACE,MED_DESCENDING_EDGE};
00044   int nbOfElt=0;
00045   for(int i=0;i<3;i++)
00046     {
00047       med_bool changement,transformation;
00048       int tmp=MEDmeshnEntity(fid,mName,dt,it,entities[i],geoElt,MED_CONNECTIVITY,MED_NODAL,
00049                              &changement,&transformation);
00050       if(tmp>nbOfElt)
00051         {
00052           nbOfElt=tmp;
00053           whichEntity=entities[i];
00054           if(i>0)
00055             std::cerr << "WARNING : MEDFile has been detected to be no compilant with MED 3 : Please change entity in MEDFile for geotype " <<  geoElt << std::endl;
00056         }
00057     }
00058   return nbOfElt>0;
00059 }
00060 
00061 int MEDFileUMeshPerType::getDim() const
00062 {
00063   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
00064   return cm.getDimension();
00065 }
00066 
00067 MEDFileUMeshPerType::MEDFileUMeshPerType(med_idt fid, const char *mName, int dt, int it, int mdim, med_geometry_type geoElt, INTERP_KERNEL::NormalizedCellType type,
00068                                          med_entity_type entity):_type(type),_entity(entity)
00069 {
00070   med_bool changement,transformation;
00071   int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_CONNECTIVITY,MED_NODAL,
00072                                  &changement,&transformation);
00073   if(type!=INTERP_KERNEL::NORM_POLYGON && type!=INTERP_KERNEL::NORM_POLYHED)
00074     {
00075       loadFromStaticType(fid,mName,dt,it,mdim,curNbOfElem,geoElt,type,entity);
00076       return;
00077     }
00078   if(type==INTERP_KERNEL::NORM_POLYGON)
00079     {
00080       loadPolyg(fid,mName,dt,it,mdim,curNbOfElem,geoElt,entity);
00081       return;
00082     }
00083   //if(type==INTERP_KERNEL::NORM_POLYHED)
00084   loadPolyh(fid,mName,dt,it,mdim,curNbOfElem,geoElt,entity);
00085 }
00086 
00087 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,
00088                                              med_entity_type entity)
00089 {
00090   _conn=DataArrayInt::New();
00091   int nbOfNodesPerCell=(geoElt%100);
00092   _conn->alloc((nbOfNodesPerCell+1)*curNbOfElem,1);
00093   _conn_index=DataArrayInt::New();
00094   _conn_index->alloc(curNbOfElem+1,1);
00095   INTERP_KERNEL::AutoPtr<int> connTab=new int[(nbOfNodesPerCell)*curNbOfElem];
00096   _num=DataArrayInt::New();
00097   _num->alloc(curNbOfElem,1);
00098   _fam=DataArrayInt::New();
00099   _fam->alloc(curNbOfElem,1);
00100   med_bool changement,transformation;
00101   INTERP_KERNEL::AutoPtr<char> noms=new char[MED_SNAME_SIZE*curNbOfElem+1];
00102   MEDmeshElementConnectivityRd(fid,mName,dt,it,entity,geoElt,MED_NODAL,MED_FULL_INTERLACE,connTab);
00103   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_FAMILY_NUMBER,MED_NODAL,&changement,&transformation)>0)
00104     {
00105       if(MEDmeshEntityFamilyNumberRd(fid,mName,dt,it,entity,geoElt,_fam->getPointer())!=0)
00106         std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
00107     }
00108   else
00109     std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
00110   if(MEDmeshnEntity(fid,mName,dt,it,entity,geoElt,MED_NUMBER,MED_NODAL,&changement,&transformation)>0)
00111     {
00112       _num=DataArrayInt::New();
00113       _num->alloc(curNbOfElem,1);
00114       if(MEDmeshEntityNumberRd(fid,mName,dt,it,entity,geoElt,_num->getPointer())!=0)
00115         _num=0;
00116     }
00117   else
00118     _num=0;
00119   int *w1=_conn->getPointer();
00120   int *w2=_conn_index->getPointer();
00121   *w2++=0;
00122   const int *wi=connTab;
00123   for(int i=0;i<curNbOfElem;i++,wi+=nbOfNodesPerCell,w2++)
00124     {
00125       *w1++=(int)type;
00126       w1=std::transform(wi,wi+nbOfNodesPerCell,w1,std::bind2nd(std::plus<int>(),-1));
00127       *w2=w2[-1]+nbOfNodesPerCell+1;
00128     }
00129 }
00130 
00131 void MEDFileUMeshPerType::loadPolyg(med_idt fid, const char *mName, int dt, int it, int mdim, int arraySize, med_geometry_type geoElt,
00132                                     med_entity_type entity)
00133 {
00134   med_bool changement,transformation;
00135   med_int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,entity,MED_POLYGON,MED_INDEX_NODE,MED_NODAL,&changement,&transformation)-1;
00136   _conn_index=DataArrayInt::New();
00137   _conn_index->alloc(curNbOfElem+1,1);
00138   _conn=DataArrayInt::New();
00139   _conn->alloc(arraySize+curNbOfElem,1);
00140   _num=DataArrayInt::New();
00141   _num->alloc(curNbOfElem,1);
00142   _fam=DataArrayInt::New();
00143   _fam->alloc(curNbOfElem,1);
00144   INTERP_KERNEL::AutoPtr<int> locConn=new int[arraySize];
00145   MEDmeshPolygonRd(fid,mName,dt,it,MED_CELL,MED_NODAL,_conn_index->getPointer(),locConn);
00146   int *w1=_conn->getPointer();
00147   int *w2=_conn_index->getPointer();
00148   const int *wi=locConn;
00149   for(int i=0;i<curNbOfElem;i++,w2++)
00150     {
00151       *w1++=(int)INTERP_KERNEL::NORM_POLYGON;
00152       const int *wi2=wi+(w2[1]-w2[0]);
00153       w1=std::transform(wi,wi2,w1,std::bind2nd(std::plus<int>(),-1));
00154       wi=wi2;
00155       *w2=*w2-1+i;
00156     }
00157   *w2=*w2-1+curNbOfElem;
00158   if(MEDmeshnEntity(fid,mName,dt,it,entity,MED_POLYGON,MED_FAMILY_NUMBER,MED_NODAL,&changement,&transformation)>0)
00159     {
00160       if(MEDmeshEntityFamilyNumberRd(fid,mName,dt,it,entity,MED_POLYGON,_fam->getPointer())!=0)
00161         std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
00162     }
00163   else
00164     std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
00165   if(MEDmeshnEntity(fid,mName,dt,it,entity,MED_POLYGON,MED_NUMBER,MED_NODAL,&changement,&transformation)>0)
00166     {
00167       if(MEDmeshEntityNumberRd(fid,mName,dt,it,entity,MED_POLYGON,_num->getPointer())!=0)
00168         _num=0;
00169     }
00170   else
00171     _num=0;
00172 }
00173 
00174 void MEDFileUMeshPerType::loadPolyh(med_idt fid, const char *mName, int dt, int it, int mdim, int connFaceLgth, med_geometry_type geoElt,
00175                                     med_entity_type entity)
00176 {
00177   med_bool changement,transformation;
00178   med_int indexFaceLgth=MEDmeshnEntity(fid,mName,dt,it,MED_CELL,MED_POLYHEDRON,MED_INDEX_NODE,MED_NODAL,&changement,&transformation);
00179   int curNbOfElem=MEDmeshnEntity(fid,mName,dt,it,MED_CELL,MED_POLYHEDRON,MED_INDEX_FACE,MED_NODAL,&changement,&transformation)-1;
00180   INTERP_KERNEL::AutoPtr<int> index=new int[curNbOfElem+1];
00181   INTERP_KERNEL::AutoPtr<int> indexFace=new int[indexFaceLgth];
00182   INTERP_KERNEL::AutoPtr<int> locConn=new int[connFaceLgth];
00183   _fam=DataArrayInt::New();
00184   _fam->alloc(curNbOfElem,1);
00185   MEDmeshPolyhedronRd(fid,mName,dt,it,MED_CELL,MED_NODAL,index,indexFace,locConn);
00186   if(MEDmeshnEntity(fid,mName,dt,it,entity,MED_POLYHEDRON,MED_FAMILY_NUMBER,MED_NODAL,&changement,&transformation)>0)
00187     {
00188       if(MEDmeshEntityFamilyNumberRd(fid,mName,dt,it,entity,MED_POLYHEDRON,_fam->getPointer())!=0)
00189         std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
00190     }
00191   else
00192     std::fill(_fam->getPointer(),_fam->getPointer()+curNbOfElem,0);
00193   int arraySize=connFaceLgth;
00194   for(int i=0;i<curNbOfElem;i++)
00195     arraySize+=index[i+1]-index[i]-1;
00196   _conn=DataArrayInt::New();
00197   _conn->alloc(arraySize+curNbOfElem,1);
00198   int *wFinalConn=_conn->getPointer();
00199   _conn_index=DataArrayInt::New();
00200   _conn_index->alloc(curNbOfElem+1,1);
00201   int *finalIndex=_conn_index->getPointer();
00202   finalIndex[0]=0;
00203   for(int i=0;i<curNbOfElem;i++)
00204     {
00205       *wFinalConn++=(int)INTERP_KERNEL::NORM_POLYHED;
00206       finalIndex[i+1]=finalIndex[i]+index[i+1]-index[i]-1+indexFace[index[i+1]-1]-indexFace[index[i]-1]+1;
00207       wFinalConn=std::transform(locConn+indexFace[index[i]-1]-1,locConn+indexFace[index[i]]-1,wFinalConn,std::bind2nd(std::plus<int>(),-1));
00208       for(int j=index[i];j<index[i+1]-1;j++)
00209         {
00210           *wFinalConn++=-1;
00211           wFinalConn=std::transform(locConn+indexFace[j]-1,locConn+indexFace[j+1]-1,wFinalConn,std::bind2nd(std::plus<int>(),-1));
00212         }
00213     }
00214   _num=DataArrayInt::New();
00215   _num->alloc(curNbOfElem,1);
00216   if(MEDmeshnEntity(fid,mName,dt,it,entity,MED_POLYHEDRON,MED_NUMBER,MED_NODAL,&changement,&transformation)>0)
00217     {
00218       if(MEDmeshEntityNumberRd(fid,mName,dt,it,entity,MED_POLYHEDRON,_num->getPointer())!=0)
00219         _num=0;
00220     }
00221   else
00222     _num=0;
00223 }
00224 
00225 void MEDFileUMeshPerType::write(med_idt fid, const char *mname, int mdim, const MEDCouplingUMesh *m, const DataArrayInt *fam, const DataArrayInt *num)
00226 {
00227   int nbOfCells=m->getNumberOfCells();
00228   if(nbOfCells<1)
00229     return ;
00230   int dt,it;
00231   double timm=m->getTime(dt,it);
00232   INTERP_KERNEL::NormalizedCellType ikt=m->getTypeOfCell(0);
00233   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(ikt);
00234   med_geometry_type curMedType=typmai3[(int)ikt];
00235   const int *conn=m->getNodalConnectivity()->getConstPointer();
00236   const int *connI=m->getNodalConnectivityIndex()->getConstPointer();
00237   if(ikt!=INTERP_KERNEL::NORM_POLYGON && ikt!=INTERP_KERNEL::NORM_POLYHED)
00238     {
00239       int nbNodesPerCell=cm.getNumberOfNodes();
00240       INTERP_KERNEL::AutoPtr<int> tab=new int[nbNodesPerCell*nbOfCells];
00241       int *w=tab;
00242       for(int i=0;i<nbOfCells;i++)
00243         w=std::transform(conn+connI[i]+1,conn+connI[i+1],w,std::bind2nd(std::plus<int>(),1));
00244       MEDmeshElementConnectivityWr(fid,mname,dt,it,timm,MED_CELL,curMedType,MED_NODAL,MED_FULL_INTERLACE,nbOfCells,tab);
00245     }
00246   else
00247     {
00248       if(ikt==INTERP_KERNEL::NORM_POLYGON)
00249         {
00250           INTERP_KERNEL::AutoPtr<int> tab1=new int[nbOfCells+1];
00251           INTERP_KERNEL::AutoPtr<int> tab2=new int[m->getMeshLength()];
00252           int *wI=tab1; *wI=1;
00253           int *w=tab2;
00254           for(int i=0;i<nbOfCells;i++,wI++)
00255             {
00256               wI[1]=wI[0]+connI[i+1]-connI[i]-1;
00257               w=std::transform(conn+connI[i]+1,conn+connI[i+1],w,std::bind2nd(std::plus<int>(),1));
00258             }
00259           MEDmeshPolygonWr(fid,mname,dt,it,timm,MED_CELL,MED_NODAL,nbOfCells+1,tab1,tab2);
00260         }
00261       else
00262         {
00263           int meshLgth=m->getMeshLength();
00264           int nbOfFaces=std::count(conn,conn+meshLgth,-1)+nbOfCells;
00265           INTERP_KERNEL::AutoPtr<int> tab1=new int[nbOfCells+1];
00266           int *w1=tab1; *w1=1;
00267           INTERP_KERNEL::AutoPtr<int> tab2=new int[nbOfFaces+1];
00268           int *w2=tab2; *w2=1;
00269           INTERP_KERNEL::AutoPtr<int> bigtab=new int[meshLgth-nbOfCells];
00270           int *bt=bigtab;
00271           for(int i=0;i<nbOfCells;i++,w1++)
00272             {
00273               int nbOfFaces2=0;
00274               for(const int *w=conn+connI[i]+1;w!=conn+connI[i+1];w2++)
00275                 {
00276                   const int *wend=std::find(w,conn+connI[i+1],-1);
00277                   bt=std::transform(w,wend,bt,std::bind2nd(std::plus<int>(),1));
00278                   int nbOfNode=std::distance(w,wend);
00279                   w2[1]=w2[0]+nbOfNode;
00280                   if(wend!=conn+connI[i+1])
00281                     w=wend+1;
00282                   else
00283                     w=wend;
00284                   nbOfFaces2++;
00285                 }
00286               w1[1]=w1[0]+nbOfFaces2;
00287             }
00288           MEDmeshPolyhedronWr(fid,mname,dt,it,timm,MED_CELL,MED_NODAL,nbOfCells+1,tab1,nbOfFaces+1,tab2,bigtab);
00289         }
00290     }
00291   if(fam)
00292     MEDmeshEntityFamilyNumberWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,fam->getConstPointer());
00293   if(num)
00294     MEDmeshEntityNumberWr(fid,mname,dt,it,MED_CELL,curMedType,nbOfCells,num->getConstPointer());
00295 }