src/MEDCoupling/MEDCouplingStructuredMesh.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 "MEDCouplingStructuredMesh.hxx"
  22 #include "MEDCouplingFieldDouble.hxx"
  23 #include "MEDCouplingMemArray.hxx"
  24 #include "MEDCouplingUMesh.hxx"
  25
  26 #include <numeric>
  27
  28 using namespace ParaMEDMEM;
  29
  30 MEDCouplingStructuredMesh::MEDCouplingStructuredMesh()
  31 {
  32 }
  33
  34 MEDCouplingStructuredMesh::MEDCouplingStructuredMesh(const MEDCouplingStructuredMesh& other, bool deepCopy):MEDCouplingMesh(other)
  35 {
  36 }
  37
  38 MEDCouplingStructuredMesh::~MEDCouplingStructuredMesh()
  39 {
  40 }
  41
  42 std::size_t MEDCouplingStructuredMesh::getHeapMemorySize() const
  43 {
  44   return MEDCouplingMesh::getHeapMemorySize();
  45 }
  46
  47 void MEDCouplingStructuredMesh::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
  48 {
  49   MEDCouplingMesh::copyTinyStringsFrom(other);
  50 }
  51
  52 bool MEDCouplingStructuredMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
  53 {
  54   return MEDCouplingMesh::isEqualIfNotWhy(other,prec,reason);
  55 }
  56
  57 INTERP_KERNEL::NormalizedCellType MEDCouplingStructuredMesh::getTypeOfCell(int cellId) const
  58 {
  59   switch(getMeshDimension())
  60     {
  61     case 3:
  62       return INTERP_KERNEL::NORM_HEXA8;
  63     case 2:
  64       return INTERP_KERNEL::NORM_QUAD4;
  65     case 1:
  66       return INTERP_KERNEL::NORM_SEG2;
  67     default:
  68       throw INTERP_KERNEL::Exception("Unexpected dimension for MEDCouplingCurveLinearMesh::getTypeOfCell !");
  69     }
  70 }
  71
  72 std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingStructuredMesh::getAllGeoTypes() const
  73 {
  74   INTERP_KERNEL::NormalizedCellType ret;
  75   switch(getMeshDimension())
  76     {
  77     case 3:
  78       ret=INTERP_KERNEL::NORM_HEXA8;
  79       break;
  80     case 2:
  81       ret=INTERP_KERNEL::NORM_QUAD4;
  82       break;
  83     case 1:
  84       ret=INTERP_KERNEL::NORM_SEG2;
  85       break;
  86     default:
  87       throw INTERP_KERNEL::Exception("Unexpected dimension for MEDCouplingStructuredMesh::getAllGeoTypes !");
  88     }
  89   std::set<INTERP_KERNEL::NormalizedCellType> ret2;
  90   ret2.insert(ret);
  91   return ret2;
  92 }
  93
  94 int MEDCouplingStructuredMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
  95 {
  96   int ret=getNumberOfCells();
  97   int dim=getMeshDimension();
  98   switch(type)
  99     {
 100     case INTERP_KERNEL::NORM_HEXA8:
 101       if(dim==3)
 102         return ret;
 103     case INTERP_KERNEL::NORM_QUAD4:
 104       if(dim==2)
 105         return ret;
 106     case INTERP_KERNEL::NORM_SEG2:
 107       if(dim==1)
 108         return ret;
 109     default:
 110       throw INTERP_KERNEL::Exception("Unexpected dimension for MEDCouplingStructuredMesh::getTypeOfCell !");
 111     }
 112   return 0;
 113 }
 114
 115 void MEDCouplingStructuredMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const
 116 {
 117   int meshDim=getMeshDimension();
 118   int tmpCell[3],tmpNode[3];
 119   getSplitCellValues(tmpCell);
 120   getSplitNodeValues(tmpNode);
 121   int tmp2[3];
 122   GetPosFromId(cellId,meshDim,tmpCell,tmp2);
 123   switch(meshDim)
 124     {
 125     case 1:
 126       conn.push_back(tmp2[0]); conn.push_back(tmp2[0]+1);
 127       break;
 128     case 2:
 129       conn.push_back(tmp2[1]*tmpCell[1]+tmp2[0]); conn.push_back(tmp2[1]*tmpCell[1]+tmp2[0]+1);
 130       conn.push_back((tmp2[1]+1)*(tmpCell[1]+1)+tmp2[0]+1); conn.push_back((tmp2[1]+1)*(tmpCell[1]+1)+tmp2[0]);
 131       break;
 132     case 3:
 133       conn.push_back(tmp2[1]*tmpCell[1]+tmp2[0]+tmp2[2]*tmpNode[2]); conn.push_back(tmp2[1]*tmpCell[1]+tmp2[0]+1+tmp2[2]*tmpNode[2]);
 134       conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+1+tmp2[2]*tmpNode[2]); conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+tmp2[2]*tmpNode[2]);
 135       conn.push_back(tmp2[1]*tmpCell[1]+tmp2[0]+(tmp2[2]+1)*tmpNode[2]); conn.push_back(tmp2[1]*tmpCell[1]+tmp2[0]+1+(tmp2[2]+1)*tmpNode[2]);
 136       conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+1+(tmp2[2]+1)*tmpNode[2]); conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+(tmp2[2]+1)*tmpNode[2]);
 137       break;
 138     default:
 139       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::getNodeIdsOfCell : big problem spacedim must be in 1,2 or 3 !");
 140     };
 141 }
 142
 143 /*!
 144  * See MEDCouplingUMesh::getDistributionOfTypes for more information
 145  */
 146 std::vector<int> MEDCouplingStructuredMesh::getDistributionOfTypes() const throw(INTERP_KERNEL::Exception)
 147 {
 148   //only one type of cell
 149   std::vector<int> ret(3);
 150   ret[0]=getTypeOfCell(0);
 151   ret[1]=getNumberOfCells();
 152   ret[2]=0; //ret[3*k+2]==0 because it has no sense here
 153   return ret;
 154 }
 155
 156 /*!
 157  * See MEDCouplingUMesh::checkTypeConsistencyAndContig for more information
 158  */
 159 DataArrayInt *MEDCouplingStructuredMesh::checkTypeConsistencyAndContig(const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
 160 {
 161   if(code.empty())
 162     throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkTypeConsistencyAndContig : code is empty, should not !");
 163   std::size_t sz=code.size();
 164   if(sz!=3)
 165     throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkTypeConsistencyAndContig : code should be of size 3 exactly !");
 166
 167   int nbCells=getNumberOfCellsWithType((INTERP_KERNEL::NormalizedCellType)code[0]);
 168   if(code[2]==-1)
 169     {
 170       if(code[1]==nbCells)
 171         return 0;
 172       else
 173         throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkTypeConsistencyAndContig : number of cells mismatch !");
 174     }
 175   else
 176     {
 177       if(code[2]<-1)
 178         throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkTypeConsistencyAndContig : code[2]<-1 mismatch !");
 179       if(code[2]>=(int)idsPerType.size())
 180         throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkTypeConsistencyAndContig : code[2]>size idsPerType !");
 181       return idsPerType[code[2]]->deepCpy();
 182     }
 183 }
 184
 185 /*!
 186  * See MEDCouplingUMesh::splitProfilePerType for more information
 187  */
 188 void MEDCouplingStructuredMesh::splitProfilePerType(const DataArrayInt *profile, std::vector<int>& code, std::vector<DataArrayInt *>& idsInPflPerType, std::vector<DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
 189 {
 190   int nbCells=getNumberOfCells();
 191   code.resize(3);
 192   code[0]=(int)getTypeOfCell(0);
 193   code[1]=nbCells;
 194   code[2]=0;
 195   idsInPflPerType.push_back(profile->deepCpy());
 196   idsPerType.push_back(profile->deepCpy());
 197 }
 198
 199 MEDCouplingUMesh *MEDCouplingStructuredMesh::buildUnstructured() const throw(INTERP_KERNEL::Exception)
 200 {
 201   int meshDim=getMeshDimension();
 202   MEDCouplingUMesh *ret=MEDCouplingUMesh::New(getName(),meshDim);
 203   DataArrayDouble *coords=getCoordinatesAndOwner();
 204   ret->setCoords(coords);
 205   coords->decrRef();
 206   switch(meshDim)
 207     {
 208     case 1:
 209       fill1DUnstructuredMesh(ret);
 210       break;
 211     case 2:
 212       fill2DUnstructuredMesh(ret);
 213       break;
 214     case 3:
 215       fill3DUnstructuredMesh(ret);
 216       break;
 217     default:
 218       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::buildUnstructured : big problem spacedim must be in 1,2 or 3 !");
 219     };
 220   return ret;
 221 }
 222
 223 MEDCouplingMesh *MEDCouplingStructuredMesh::buildPart(const int *start, const int *end) const
 224 {
 225   MEDCouplingUMesh *um=buildUnstructured();
 226   MEDCouplingMesh *ret=um->buildPart(start,end);
 227   um->decrRef();
 228   return ret;
 229 }
 230
 231 MEDCouplingMesh *MEDCouplingStructuredMesh::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
 232 {
 233   MEDCouplingUMesh *um=buildUnstructured();
 234   MEDCouplingMesh *ret=um->buildPartAndReduceNodes(start,end,arr);
 235   um->decrRef();
 236   return ret;
 237 }
 238
 239 DataArrayInt *MEDCouplingStructuredMesh::simplexize(int policy) throw(INTERP_KERNEL::Exception)
 240 {
 241   throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::simplexize : not available for Cartesian mesh !");
 242 }
 243
 244 MEDCouplingFieldDouble *MEDCouplingStructuredMesh::buildOrthogonalField() const
 245 {
 246   if(getMeshDimension()!=2)
 247     throw INTERP_KERNEL::Exception("Expected a MEDCouplingStructuredMesh with meshDim == 2 !");
 248   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
 249   DataArrayDouble *array=DataArrayDouble::New();
 250   int nbOfCells=getNumberOfCells();
 251   array->alloc(nbOfCells,3);
 252   double *vals=array->getPointer();
 253   for(int i=0;i<nbOfCells;i++)
 254     { vals[3*i]=0.; vals[3*i+1]=0.; vals[3*i+2]=1.; }
 255   ret->setArray(array);
 256   array->decrRef();
 257   ret->setMesh(this);
 258   return ret;
 259 }
 260
 261 void MEDCouplingStructuredMesh::fill1DUnstructuredMesh(MEDCouplingUMesh *m) const
 262 {
 263   int nbOfCells=-1;
 264   getNodeGridStructure(&nbOfCells);
 265   nbOfCells--;
 266   DataArrayInt *connI=DataArrayInt::New();
 267   connI->alloc(nbOfCells+1,1);
 268   int *ci=connI->getPointer();
 269   DataArrayInt *conn=DataArrayInt::New();
 270   conn->alloc(3*nbOfCells,1);
 271   ci[0]=0;
 272   int *cp=conn->getPointer();
 273   for(int i=0;i<nbOfCells;i++)
 274     {
 275       cp[3*i]=(int)INTERP_KERNEL::NORM_SEG2;
 276       cp[3*i+1]=i;
 277       cp[3*i+2]=i+1;
 278       ci[i+1]=3*(i+1);
 279     }
 280   m->setConnectivity(conn,connI,true);
 281   conn->decrRef();
 282   connI->decrRef();
 283 }
 284
 285 void MEDCouplingStructuredMesh::fill2DUnstructuredMesh(MEDCouplingUMesh *m) const
 286 {
 287   int ns[2];
 288   getNodeGridStructure(ns);
 289   int n1=ns[0]-1;
 290   int n2=ns[1]-1;
 291   DataArrayInt *connI=DataArrayInt::New();
 292   connI->alloc(n1*n2+1,1);
 293   int *ci=connI->getPointer();
 294   DataArrayInt *conn=DataArrayInt::New();
 295   conn->alloc(5*n1*n2,1);
 296   ci[0]=0;
 297   int *cp=conn->getPointer();
 298   int pos=0;
 299   for(int j=0;j<n2;j++)
 300     for(int i=0;i<n1;i++,pos++)
 301       {
 302         cp[5*pos]=(int)INTERP_KERNEL::NORM_QUAD4;
 303         cp[5*pos+1]=i+1+j*(n1+1);
 304         cp[5*pos+2]=i+j*(n1+1);
 305         cp[5*pos+3]=i+(j+1)*(n1+1);
 306         cp[5*pos+4]=i+1+(j+1)*(n1+1);
 307         ci[pos+1]=5*(pos+1);
 308     }
 309   m->setConnectivity(conn,connI,true);
 310   conn->decrRef();
 311   connI->decrRef();
 312 }
 313
 314 void MEDCouplingStructuredMesh::fill3DUnstructuredMesh(MEDCouplingUMesh *m) const
 315 {
 316   int ns[3];
 317   getNodeGridStructure(ns);
 318   int n1=ns[0]-1;
 319   int n2=ns[1]-1;
 320   int n3=ns[2]-1;
 321   DataArrayInt *connI=DataArrayInt::New();
 322   connI->alloc(n1*n2*n3+1,1);
 323   int *ci=connI->getPointer();
 324   DataArrayInt *conn=DataArrayInt::New();
 325   conn->alloc(9*n1*n2*n3,1);
 326   ci[0]=0;
 327   int *cp=conn->getPointer();
 328   int pos=0;
 329   for(int k=0;k<n3;k++)
 330     for(int j=0;j<n2;j++)
 331       for(int i=0;i<n1;i++,pos++)
 332         {
 333           cp[9*pos]=(int)INTERP_KERNEL::NORM_HEXA8;
 334           int tmp=(n1+1)*(n2+1);
 335           cp[9*pos+1]=i+1+j*(n1+1)+k*tmp;
 336           cp[9*pos+2]=i+j*(n1+1)+k*tmp;
 337           cp[9*pos+3]=i+(j+1)*(n1+1)+k*tmp;
 338           cp[9*pos+4]=i+1+(j+1)*(n1+1)+k*tmp;
 339           cp[9*pos+5]=i+1+j*(n1+1)+(k+1)*tmp;
 340           cp[9*pos+6]=i+j*(n1+1)+(k+1)*tmp;
 341           cp[9*pos+7]=i+(j+1)*(n1+1)+(k+1)*tmp;
 342           cp[9*pos+8]=i+1+(j+1)*(n1+1)+(k+1)*tmp;
 343           ci[pos+1]=9*(pos+1);
 344         }
 345   m->setConnectivity(conn,connI,true);
 346   conn->decrRef();
 347   connI->decrRef();
 348 }
 349
 350 int MEDCouplingStructuredMesh::getCellIdFromPos(int i, int j, int k) const
 351 {
 352   int tmp[3]={i,j,k};
 353   int tmp2[3];
 354   int meshDim=getMeshDimension();
 355   getSplitCellValues(tmp2);
 356   std::transform(tmp,tmp+meshDim,tmp2,tmp,std::multiplies<int>());
 357   return std::accumulate(tmp,tmp+meshDim,0);
 358 }
 359
 360 int MEDCouplingStructuredMesh::getNodeIdFromPos(int i, int j, int k) const
 361 {
 362   int tmp[3]={i,j,k};
 363   int tmp2[3];
 364   int meshDim=getMeshDimension();
 365   getSplitNodeValues(tmp2);
 366   std::transform(tmp,tmp+meshDim,tmp2,tmp,std::multiplies<int>());
 367   return std::accumulate(tmp,tmp+meshDim,0);
 368 }
 369
 370 void MEDCouplingStructuredMesh::GetPosFromId(int nodeId, int meshDim, const int *split, int *res)
 371 {
 372   int work=nodeId;
 373   for(int i=meshDim-1;i>=0;i--)
 374     {
 375       int pos=work/split[i];
 376       work=work%split[i];
 377       res[i]=pos;
 378     }
 379 }