src/MEDCoupling/MEDCouplingVoronoi.cxx

   1 // Copyright (C) 2007-2017  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, 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 R&D)
  20
  21 #include "MEDCouplingVoronoi.hxx"
  22 #include "MEDCoupling1GTUMesh.hxx"
  23 #include "MEDCouplingCMesh.hxx"
  24 #include "MCAuto.txx"
  25
  26 using namespace MEDCoupling;
  27
  28 MCAuto<MEDCouplingUMesh> ComputeBigCellFrom(const double pt1[2], const double pt2[2], const std::vector<double>& bbox, double eps)
  29 {
  30   static const double FACT=1.2;
  31   MCAuto<MEDCouplingCMesh> m(MEDCouplingCMesh::New());
  32   MCAuto<DataArrayDouble> arr1(DataArrayDouble::New()); arr1->alloc(2,1) ; arr1->setIJ(0,0,bbox[0]); arr1->setIJ(1,0,bbox[1]);
  33   MCAuto<DataArrayDouble> arr2(DataArrayDouble::New()); arr2->alloc(2,1) ; arr2->setIJ(0,0,bbox[2]); arr2->setIJ(1,0,bbox[3]);
  34   m->setCoords(arr1,arr2);
  35   static const double PT[2]={0.,0.};
  36   m->scale(PT,FACT);
  37   MCAuto<MEDCouplingUMesh> mu(m->buildUnstructured());
  38   double l(std::max(bbox[1]-bbox[0],bbox[3]-bbox[2]));
  39   double middle[2]={(pt1[0]+pt2[0])/2.,(pt1[1]+pt2[1])/2.};
  40   double v[2]={pt1[0],pt1[1]};
  41   DataArrayDouble::Rotate2DAlg(middle,M_PI/2,1,v,v);
  42   v[0]=middle[0]-v[0]; v[1]=middle[1]-v[1];
  43   {
  44     double nor(sqrt(v[0]*v[0]+v[1]*v[1]));
  45     v[0]/=nor; v[1]/=nor;
  46   }
  47   MCAuto<MEDCouplingUMesh> line(MEDCouplingUMesh::New("line",1));
  48   {
  49     MCAuto<DataArrayDouble> coo(DataArrayDouble::New()); coo->alloc(2,2);
  50     coo->setIJ(0,0,middle[0]-2.*l*v[0]); coo->setIJ(0,1,middle[1]-2.*l*v[1]); coo->setIJ(1,0,middle[0]+2.*l*v[0]); coo->setIJ(1,1,middle[1]+2.*l*v[1]);
  51     line->setCoords(coo);
  52   }
  53   line->allocateCells();
  54   static const int CONN[2]={0,1};
  55   line->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,CONN);
  56   MCAuto<MEDCouplingUMesh> sp2,sp1;
  57   {
  58     DataArrayInt *cellNb1(0),*cellNb2(0);
  59     MEDCouplingUMesh *sp2Pt(0),*sp1Pt(0);
  60     MEDCouplingUMesh::Intersect2DMeshWith1DLine(mu,line,eps,sp2Pt,sp1Pt,cellNb1,cellNb2);
  61     sp1=sp1Pt; sp2=sp2Pt;
  62     MCAuto<DataArrayInt> cellNb10(cellNb1),cellNb20(cellNb2);
  63   }
  64   std::vector<int> ccp;
  65   sp2->getCellsContainingPoint(pt1,eps,ccp);
  66   if(ccp.size()!=1)
  67     throw INTERP_KERNEL::Exception("ComputeBigCellFrom : expected single element !");
  68   MCAuto<MEDCouplingUMesh> ret(sp2->buildPartOfMySelfSlice(ccp[0],ccp[0]+1,1,true));
  69   ret->zipCoords();
  70   return ret;
  71 }
  72
  73 MCAuto<MEDCouplingUMesh> MergeVorCells(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
  74 {
  75   std::size_t sz(vcs.size());
  76   if(sz<1)
  77     throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
  78   if(sz==1)
  79     return vcs[0];
  80   MCAuto<MEDCouplingUMesh> p;
  81   {
  82     std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
  83     p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
  84   }
  85   p->zipCoords();
  86   {
  87     bool dummy; int dummy2;
  88     MCAuto<DataArrayInt> dummy3(p->mergeNodes(eps,dummy,dummy2));
  89   }
  90   MCAuto<DataArrayInt> edgeToKeep;
  91   MCAuto<MEDCouplingUMesh> p0;
  92   {
  93     MCAuto<DataArrayInt> d(DataArrayInt::New()),di(DataArrayInt::New()),rd(DataArrayInt::New()),rdi(DataArrayInt::New());
  94     p0=p->buildDescendingConnectivity(d,di,rd,rdi);
  95     MCAuto<DataArrayInt> dsi(rdi->deltaShiftIndex());
  96     edgeToKeep=dsi->findIdsEqual(1);
  97   }
  98   MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
  99   skinOfRes->zipCoords();
 100   if(skinOfRes->getNumberOfCells()!=skinOfRes->getNumberOfNodes())
 101     throw INTERP_KERNEL::Exception("MergeVorCells : result of merge looks bad !");
 102   MCAuto<DataArrayInt> d(skinOfRes->orderConsecutiveCells1D());
 103   MCAuto<MEDCoupling1SGTUMesh> skinOfRes2;
 104   {
 105     MCAuto<MEDCouplingUMesh> part(skinOfRes->buildPartOfMySelf(d->begin(),d->end()));
 106     skinOfRes2=MEDCoupling1SGTUMesh::New(part);
 107   }
 108   MCAuto<DataArrayInt> c(skinOfRes2->getNodalConnectivity()->deepCopy());
 109   c->circularPermutation(1);
 110   c->rearrange(2);
 111   std::vector< MCAuto<DataArrayInt> > vdi(c->explodeComponents());
 112   if(!vdi[0]->isEqual(*vdi[1]))
 113     throw INTERP_KERNEL::Exception("MergeVorCells : internal error !");
 114   MCAuto<MEDCouplingUMesh> m(MEDCouplingUMesh::New("",2));
 115   m->setCoords(skinOfRes2->getCoords());
 116   m->allocateCells();
 117   m->insertNextCell(INTERP_KERNEL::NORM_POLYGON,vdi[0]->getNumberOfTuples(),vdi[0]->begin());
 118   return m;
 119 }
 120
 121 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronoize2D(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps)
 122 {
 123   if(!m || !points)
 124     throw INTERP_KERNEL::Exception("Voronoize2D : null pointer !");
 125   m->checkConsistencyLight();
 126   points->checkAllocated();
 127   if(m->getMeshDimension()!=2 || m->getSpaceDimension()!=2 || points->getNumberOfComponents()!=2)
 128     throw INTERP_KERNEL::Exception("Voronoize2D : spacedim must be equal to 2 and meshdim also equal to 2 !");
 129   if(m->getNumberOfCells()!=1)
 130     throw INTERP_KERNEL::Exception("Voronoize2D : mesh is expected to have only one cell !");
 131   int nbPts(points->getNumberOfTuples());
 132   if(nbPts<1)
 133     throw INTERP_KERNEL::Exception("Voronoize2D : at least one point expected !");
 134   std::vector<double> bbox(4);
 135   m->getBoundingBox(&bbox[0]);
 136   std::vector< MCAuto<MEDCouplingUMesh> > l0(1,MCAuto<MEDCouplingUMesh>(m->deepCopy()));
 137   const double *pts(points->begin());
 138   for(int i=1;i<nbPts;i++)
 139     {
 140       MCAuto<MEDCouplingUMesh> vorTess;
 141       {
 142         std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 143         vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
 144       }
 145       {
 146         bool dummy;
 147         int newNbNodes;
 148         MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
 149       }
 150       std::vector<int> polygsToIterOn;
 151       const double *pt(pts+i*2);
 152       vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
 153       if(polygsToIterOn.size()<1)
 154         throw INTERP_KERNEL::Exception("Voronoize2D : presence of a point outside the given cell !");
 155       std::set<int> elemsToDo,elemsDone; elemsToDo.insert(polygsToIterOn[0]);
 156       std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
 157       while(!elemsToDo.empty())
 158         {
 159           int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly);
 160           const double *seed(pts+2*poly);
 161           MCAuto<MEDCouplingUMesh> cell(ComputeBigCellFrom(pt,seed,bbox,eps));
 162           MCAuto<MEDCouplingUMesh> tile(l0[poly]);
 163           tile->zipCoords();
 164           MCAuto<MEDCouplingUMesh> a;
 165           MCAuto<DataArrayInt> b,c;
 166           {
 167             DataArrayInt *bPtr(0),*cPtr(0);
 168             a=MEDCouplingUMesh::Intersect2DMeshes(tile,cell,eps,bPtr,cPtr);
 169             b=bPtr; c=cPtr;
 170           }
 171           MCAuto<DataArrayInt> part(c->findIdsEqual(-1));
 172           if(part->getNumberOfTuples()!=1)
 173             throw INTERP_KERNEL::Exception("Voronoize2D : internal error");
 174           MCAuto<MEDCouplingUMesh> newVorCell;
 175           {
 176             MCAuto<DataArrayInt> tmp(part->buildComplement(a->getNumberOfCells()));
 177             newVorCell=a->buildPartOfMySelf(tmp->begin(),tmp->end());
 178           }
 179           newVorCell->zipCoords();
 180           MCAuto<MEDCouplingUMesh> modifiedCell(a->buildPartOfMySelf(part->begin(),part->end()));
 181           modifiedCell->zipCoords();
 182           l0[poly]=modifiedCell;
 183           //
 184           MCAuto<DataArrayInt> ids;
 185           {
 186             DataArrayInt *tmp(0);
 187             bool sta(a->getCoords()->areIncludedInMe(cell->getCoords(),eps,tmp));
 188             ids=tmp;
 189             if(!sta)
 190               throw INTERP_KERNEL::Exception("Voronoize2D : internal error 2 !");
 191           }
 192           MCAuto<DataArrayDouble> newCoords;
 193           {
 194             MCAuto<DataArrayInt> tmp(ids->buildComplement(a->getNumberOfNodes()));
 195             newCoords=a->getCoords()->selectByTupleId(tmp->begin(),tmp->end());
 196           }
 197           const double *cPtr(newCoords->begin());
 198           for(int i=0;i<newCoords->getNumberOfTuples();i++,cPtr+=2)
 199             {
 200               std::set<int> zeCandidates;
 201               {
 202                 std::vector<int> zeCandidatesTmp;
 203                 vorTess->getCellsContainingPoint(cPtr,eps,zeCandidatesTmp);
 204                 zeCandidates.insert(zeCandidatesTmp.begin(),zeCandidatesTmp.end());
 205               }
 206               std::set<int> tmp,newElementsToDo;
 207               std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp,tmp.begin()));
 208               std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp.begin(),tmp.end(),std::inserter(newElementsToDo,newElementsToDo.begin()));
 209               elemsToDo=newElementsToDo;
 210             }
 211           newVorCells.push_back(newVorCell);
 212         }
 213       l0.push_back(MergeVorCells(newVorCells,eps));
 214     }
 215   std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 216   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
 217   return ret;
 218 }