src/MEDCoupling/MEDCouplingUMesh_internal.hxx

   1 // Copyright (C) 2007-2024  CEA, EDF
   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, or (at your option) any later version.
   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 (EdF)
  20
  21 using namespace MEDCoupling;
  22
  23 class MinusOneSonsGenerator
  24 {
  25 public:
  26   MinusOneSonsGenerator(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
  27   unsigned getNumberOfSons2(const mcIdType *conn, mcIdType lgth) const { return _cm.getNumberOfSons2(conn,lgth); }
  28   unsigned fillSonCellNodalConnectivity2(int sonId, const mcIdType *nodalConn, mcIdType lgth, mcIdType *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillSonCellNodalConnectivity2(sonId,nodalConn,lgth,sonNodalConn,typeOfSon); }
  29   static const int DELTA=1;
  30 private:
  31   const INTERP_KERNEL::CellModel& _cm;
  32 };
  33
  34 class MinusOneSonsGeneratorBiQuadratic
  35 {
  36 public:
  37   MinusOneSonsGeneratorBiQuadratic(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
  38   unsigned getNumberOfSons2(const mcIdType *conn, mcIdType lgth) const { return _cm.getNumberOfSons2(conn,lgth); }
  39   unsigned fillSonCellNodalConnectivity2(int sonId, const mcIdType *nodalConn, mcIdType lgth, mcIdType *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillSonCellNodalConnectivity4(sonId,nodalConn,lgth,sonNodalConn,typeOfSon); }
  40   static const int DELTA=1;
  41 private:
  42   const INTERP_KERNEL::CellModel& _cm;
  43 };
  44
  45 class MinusTwoSonsGenerator
  46 {
  47 public:
  48   MinusTwoSonsGenerator(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
  49   unsigned getNumberOfSons2(const mcIdType *conn, mcIdType lgth) const { return _cm.getNumberOfEdgesIn3D(conn,lgth); }
  50   unsigned fillSonCellNodalConnectivity2(int sonId, const mcIdType *nodalConn, mcIdType lgth, mcIdType *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillSonEdgesNodalConnectivity3D(sonId,nodalConn,lgth,sonNodalConn,typeOfSon); }
  51   static const int DELTA=2;
  52 private:
  53   const INTERP_KERNEL::CellModel& _cm;
  54 };
  55
  56 class MicroEdgesGenerator2D
  57 {
  58 public:
  59   MicroEdgesGenerator2D(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
  60   unsigned getNumberOfSons2(const mcIdType *conn, mcIdType lgth) const { return _cm.getNumberOfMicroEdges(); }
  61   unsigned fillSonCellNodalConnectivity2(int sonId, const mcIdType *nodalConn, mcIdType lgth, mcIdType *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillMicroEdgeNodalConnectivity(sonId,nodalConn,sonNodalConn,typeOfSon); }
  62   static const int DELTA=1;
  63 private:
  64   const INTERP_KERNEL::CellModel& _cm;
  65 };
  66
  67 class MicroEdgesGenerator3D
  68 {
  69 public:
  70   MicroEdgesGenerator3D(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
  71   unsigned getNumberOfSons2(const mcIdType *conn, mcIdType lgth) const { return _cm.getNumberOfMicroEdges(); }
  72   unsigned fillSonCellNodalConnectivity2(int sonId, const mcIdType *nodalConn, mcIdType lgth, mcIdType *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillMicroEdgeNodalConnectivity(sonId,nodalConn,sonNodalConn,typeOfSon); }
  73   static const int DELTA=2;
  74 private:
  75   const INTERP_KERNEL::CellModel& _cm;
  76 };
  77
  78 mcIdType MEDCouplingFastNbrer(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *conn2);
  79 mcIdType MEDCouplingOrientationSensitiveNbrer(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *conn2);
  80
  81 namespace MEDCoupling
  82 {
  83   template<const int SPACEDIMM>
  84   class DummyClsMCUG
  85   {
  86   public:
  87     static const int MY_SPACEDIM=SPACEDIMM;
  88     static const int MY_MESHDIM=8;
  89     typedef mcIdType MyConnType;
  90     static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
  91     // begin
  92     // useless, but for windows compilation ...
  93     const double* getCoordinatesPtr() const { return 0; }
  94     const MyConnType* getConnectivityPtr() const { return 0; }
  95     const MyConnType* getConnectivityIndexPtr() const { return 0; }
  96     INTERP_KERNEL::NormalizedCellType getTypeOfElement(MyConnType) const { return (INTERP_KERNEL::NormalizedCellType)0; }
  97     // end
  98   };
  99 }
 100
 101 template<int SPACEDIM>
 102 void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const double *pos, mcIdType nbOfPoints,
 103                                                    double eps, MCAuto<DataArrayIdType>& elts, MCAuto<DataArrayIdType>& eltsIndex, std::function<bool(INTERP_KERNEL::NormalizedCellType,int)> sensibilityTo2DQuadraticLinearCellsFunc) const
 104 {
 105   // Override precision for this method only:
 106   INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
 107
 108   elts=DataArrayIdType::New(); eltsIndex=DataArrayIdType::New(); eltsIndex->alloc(nbOfPoints+1,1); eltsIndex->setIJ(0,0,0); elts->alloc(0,1);
 109   mcIdType *eltsIndexPtr(eltsIndex->getPointer());
 110   MCAuto<DataArrayDouble> bboxArr(getBoundingBoxForBBTree(eps));
 111   const double *bbox(bboxArr->begin());
 112   mcIdType nbOfCells=getNumberOfCells();
 113   const mcIdType *conn=_nodal_connec->getConstPointer();
 114   const mcIdType *connI=_nodal_connec_index->getConstPointer();
 115   double bb[2*SPACEDIM];
 116   BBTree<SPACEDIM,mcIdType> myTree(&bbox[0],0,0,nbOfCells,-eps);
 117   for(mcIdType i=0;i<nbOfPoints;i++)
 118     {
 119       eltsIndexPtr[i+1]=eltsIndexPtr[i];
 120       for(int j=0;j<SPACEDIM;j++)
 121         {
 122           bb[2*j]=pos[SPACEDIM*i+j];
 123           bb[2*j+1]=pos[SPACEDIM*i+j];
 124         }
 125       std::vector<mcIdType> candidates;
 126       myTree.getIntersectingElems(bb,candidates);
 127       for(std::vector<mcIdType>::const_iterator iter=candidates.begin();iter!=candidates.end();iter++)
 128         {
 129           mcIdType sz(connI[(*iter)+1]-connI[*iter]-1);
 130           INTERP_KERNEL::NormalizedCellType ct((INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]]);
 131           bool status(false);
 132           // [ABN] : point locator algorithms are not impl. for POLY or QPOLY in spaceDim3
 133           if(  SPACEDIM!=2 &&
 134               (ct == INTERP_KERNEL::NORM_POLYGON || sensibilityTo2DQuadraticLinearCellsFunc(ct,_mesh_dim)))
 135             throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellsContainingPointsAlg : not implemented yet for POLYGON and QPOLYGON in spaceDim 3 !");
 136           // Keep calling simple algorithm when this is desired and simple for speed reasons:
 137           if (SPACEDIM == 2 && ct != INTERP_KERNEL::NORM_POLYGON && !sensibilityTo2DQuadraticLinearCellsFunc(ct,_mesh_dim))
 138             status=INTERP_KERNEL::PointLocatorAlgos<DummyClsMCUG<2> >::isElementContainsPointAlgo2DSimple2(pos+i*SPACEDIM,ct,coords,conn+connI[*iter]+1,sz,eps);
 139           else
 140             status=INTERP_KERNEL::PointLocatorAlgos<DummyClsMCUG<SPACEDIM> >::isElementContainsPoint(pos+i*SPACEDIM,ct,coords,conn+connI[*iter]+1,sz,eps);
 141           if(status)
 142             {
 143               eltsIndexPtr[i+1]++;
 144               elts->pushBackSilent(*iter);
 145             }
 146         }
 147     }
 148 }
 149
 150 /*!
 151  * \b WARNING this method do the assumption that connectivity lies on the coordinates set.
 152  * For speed reasons no check of this will be done.
 153  */
 154 template<class SonsGenerator>
 155 MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivityGen(DataArrayIdType *desc, DataArrayIdType *descIndx, DataArrayIdType *revDesc, DataArrayIdType *revDescIndx, DimM1DescNbrer nbrer) const
 156 {
 157   if(!desc || !descIndx || !revDesc || !revDescIndx)
 158     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildDescendingConnectivityGen : present of a null pointer in input !");
 159   checkConnectivityFullyDefined();
 160   mcIdType nbOfCells=getNumberOfCells();
 161   mcIdType nbOfNodes=getNumberOfNodes();
 162   MCAuto<DataArrayIdType> revNodalIndx=DataArrayIdType::New(); revNodalIndx->alloc(nbOfNodes+1,1); revNodalIndx->fillWithZero();
 163   mcIdType *revNodalIndxPtr=revNodalIndx->getPointer();
 164   const mcIdType *conn=_nodal_connec->getConstPointer();
 165   const mcIdType *connIndex=_nodal_connec_index->getConstPointer();
 166   std::string name="Mesh constituent of "; name+=getName();
 167   MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New(name,getMeshDimension()-SonsGenerator::DELTA);
 168   ret->setCoords(getCoords());
 169   ret->allocateCells(2*nbOfCells);
 170   descIndx->alloc(nbOfCells+1,1);
 171   MCAuto<DataArrayIdType> revDesc2(DataArrayIdType::New()); revDesc2->reserve(2*nbOfCells);
 172   mcIdType *descIndxPtr=descIndx->getPointer(); *descIndxPtr++=0;
 173   for(mcIdType eltId=0;eltId<nbOfCells;eltId++,descIndxPtr++)
 174     {
 175       mcIdType pos=connIndex[eltId];
 176       mcIdType posP1=connIndex[eltId+1];
 177       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[pos]);
 178       SonsGenerator sg(cm);
 179       unsigned nbOfSons=sg.getNumberOfSons2(conn+pos+1,posP1-pos-1);
 180       INTERP_KERNEL::AutoPtr<mcIdType> tmp=new mcIdType[posP1-pos];
 181       for(unsigned i=0;i<nbOfSons;i++)
 182         {
 183           INTERP_KERNEL::NormalizedCellType cmsId;
 184           unsigned nbOfNodesSon=sg.fillSonCellNodalConnectivity2(i,conn+pos+1,posP1-pos-1,tmp,cmsId);
 185           for(unsigned k=0;k<nbOfNodesSon;k++)
 186             if(tmp[k]>=0)
 187               revNodalIndxPtr[tmp[k]+1]++;
 188           ret->insertNextCell(cmsId,nbOfNodesSon,tmp);
 189           revDesc2->pushBackSilent(eltId);
 190         }
 191       descIndxPtr[0]=descIndxPtr[-1]+ToIdType(nbOfSons);
 192     }
 193   mcIdType nbOfCellsM1=ret->getNumberOfCells();
 194   std::transform(revNodalIndxPtr+1,revNodalIndxPtr+nbOfNodes+1,revNodalIndxPtr,revNodalIndxPtr+1,std::plus<mcIdType>());
 195   MCAuto<DataArrayIdType> revNodal=DataArrayIdType::New(); revNodal->alloc(revNodalIndx->back(),1);
 196   std::fill(revNodal->getPointer(),revNodal->getPointer()+revNodalIndx->back(),-1);
 197   mcIdType *revNodalPtr=revNodal->getPointer();
 198   const mcIdType *connM1=ret->getNodalConnectivity()->getConstPointer();
 199   const mcIdType *connIndexM1=ret->getNodalConnectivityIndex()->getConstPointer();
 200   for(mcIdType eltId=0;eltId<nbOfCellsM1;eltId++)
 201     {
 202       const mcIdType *strtNdlConnOfCurCell=connM1+connIndexM1[eltId]+1;
 203       const mcIdType *endNdlConnOfCurCell=connM1+connIndexM1[eltId+1];
 204       for(const mcIdType *iter=strtNdlConnOfCurCell;iter!=endNdlConnOfCurCell;iter++)
 205         if(*iter>=0)//for polyhedrons
 206           *std::find_if(revNodalPtr+revNodalIndxPtr[*iter],revNodalPtr+revNodalIndxPtr[*iter+1],std::bind(std::equal_to<mcIdType>(),std::placeholders::_1,-1))=eltId;
 207     }
 208   //
 209   DataArrayIdType *commonCells=0,*commonCellsI=0;
 210   FindCommonCellsAlg(3,0,ret->getNodalConnectivity(),ret->getNodalConnectivityIndex(),revNodal,revNodalIndx,commonCells,commonCellsI);
 211   MCAuto<DataArrayIdType> commonCellsTmp(commonCells),commonCellsITmp(commonCellsI);
 212   const mcIdType *commonCellsPtr(commonCells->getConstPointer()),*commonCellsIPtr(commonCellsI->getConstPointer());
 213   mcIdType newNbOfCellsM1=-1;
 214   MCAuto<DataArrayIdType> o2nM1=DataArrayIdType::ConvertIndexArrayToO2N(nbOfCellsM1,commonCells->begin(),
 215                                                                                                             commonCellsI->begin(),commonCellsI->end(),newNbOfCellsM1);
 216   std::vector<bool> isImpacted(nbOfCellsM1,false);
 217   for(const mcIdType *work=commonCellsI->begin();work!=commonCellsI->end()-1;work++)
 218     for(mcIdType work2=work[0];work2!=work[1];work2++)
 219       isImpacted[commonCellsPtr[work2]]=true;
 220   const mcIdType *o2nM1Ptr=o2nM1->getConstPointer();
 221   MCAuto<DataArrayIdType> n2oM1=o2nM1->invertArrayO2N2N2OBis(newNbOfCellsM1);
 222   const mcIdType *n2oM1Ptr=n2oM1->getConstPointer();
 223   MCAuto<MEDCouplingUMesh> ret2=static_cast<MEDCouplingUMesh *>(ret->buildPartOfMySelf(n2oM1->begin(),n2oM1->end(),true));
 224   ret2->copyTinyInfoFrom(this);
 225   desc->alloc(descIndx->back(),1);
 226   mcIdType *descPtr=desc->getPointer();
 227   const INTERP_KERNEL::CellModel& cmsDft=INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_POINT1);
 228   for(mcIdType i=0;i<nbOfCellsM1;i++,descPtr++)
 229     {
 230       if(!isImpacted[i])
 231         *descPtr=nbrer(o2nM1Ptr[i],0,cmsDft,false,0,0);
 232       else
 233         {
 234           if(i!=n2oM1Ptr[o2nM1Ptr[i]])
 235             {
 236               const INTERP_KERNEL::CellModel& cms=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connM1[connIndexM1[i]]);
 237               *descPtr=nbrer(o2nM1Ptr[i],connIndexM1[i+1]-connIndexM1[i]-1,cms,true,connM1+connIndexM1[n2oM1Ptr[o2nM1Ptr[i]]]+1,connM1+connIndexM1[i]+1);
 238             }
 239           else
 240             *descPtr=nbrer(o2nM1Ptr[i],0,cmsDft,false,0,0);
 241         }
 242     }
 243   revDesc->reserve(newNbOfCellsM1);
 244   revDescIndx->alloc(newNbOfCellsM1+1,1);
 245   mcIdType *revDescIndxPtr=revDescIndx->getPointer(); *revDescIndxPtr++=0;
 246   const mcIdType *revDesc2Ptr=revDesc2->getConstPointer();
 247   for(mcIdType i=0;i<newNbOfCellsM1;i++,revDescIndxPtr++)
 248     {
 249       mcIdType oldCellIdM1=n2oM1Ptr[i];
 250       if(!isImpacted[oldCellIdM1])
 251         {
 252           revDesc->pushBackSilent(revDesc2Ptr[oldCellIdM1]);
 253           revDescIndxPtr[0]=revDescIndxPtr[-1]+1;
 254         }
 255       else
 256         {
 257           for(mcIdType j=commonCellsIPtr[0];j<commonCellsIPtr[1];j++)
 258             revDesc->pushBackSilent(revDesc2Ptr[commonCellsPtr[j]]);
 259           revDescIndxPtr[0]=revDescIndxPtr[-1]+commonCellsIPtr[1]-commonCellsIPtr[0];
 260           commonCellsIPtr++;
 261         }
 262     }
 263   //
 264   return ret2.retn();
 265 }
 266
 267