From: Anthony Geay Date: Fri, 27 Jan 2017 14:04:48 +0000 (+0100) Subject: A first fast implementation of 2D Voronoi X-Git-Tag: V8_3_0a2~30 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=331685311f395c0b0654223c1cffddd254ed49cb;p=tools%2Fmedcoupling.git A first fast implementation of 2D Voronoi --- diff --git a/src/MEDCoupling/CMakeLists.txt b/src/MEDCoupling/CMakeLists.txt index 3a70181ea..4c66e6dd6 100644 --- a/src/MEDCoupling/CMakeLists.txt +++ b/src/MEDCoupling/CMakeLists.txt @@ -63,6 +63,7 @@ SET(medcoupling_SOURCES MEDCouplingMatrix.cxx MEDCouplingPartDefinition.cxx MEDCouplingSkyLineArray.cxx + MEDCouplingVoronoi.cxx ../ICoCo/ICoCoField.cxx ../ICoCo/ICoCoMEDField.cxx ) diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx index 4c52fb881..92d8e863e 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx @@ -25,7 +25,8 @@ #include "MEDCouplingUMesh.hxx" #include "MEDCouplingTimeDiscretization.hxx" #include "MEDCouplingFieldDiscretization.hxx" -#include "MCAuto.hxx" +#include "MCAuto.txx" +#include "MEDCouplingVoronoi.hxx" #include "MEDCouplingNatureOfField.hxx" #include "InterpKernelAutoPtr.hxx" @@ -2202,6 +2203,19 @@ bool MEDCouplingFieldDouble::simplexize(int policy) return true; } +/*! + * This method makes the hypothesis that \a this is a Gauss field. This method returns a newly created field on cells with same number of tuples than \a this. + * Each Gauss points in \a this is replaced by a polygon or polyhedron cell with associated region following Voronoi algorithm. + */ +MCAuto MEDCouplingFieldDouble::voronoize(double eps) const +{ + checkConsistencyLight(); + const MEDCouplingMesh *mesh(getMesh()); + if(mesh->getSpaceDimension()==2 && mesh->getSpaceDimension()==2) + return voronoize2D(eps); + throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::voronoize : only 2D is supported for the moment !"); +} + /*! * This is expected to be a 3 components vector field on nodes (if not an exception will be thrown). \a this is also expected to lie on a MEDCouplingPointSet mesh. * Finaly \a this is also expected to be consistent. @@ -3068,6 +3082,65 @@ std::string MEDCouplingFieldDouble::WriteVTK(const std::string& fileName, const return ret; } +MCAuto MEDCouplingFieldDouble::voronoize2D(double eps) const +{ + checkConsistencyLight(); + const MEDCouplingMesh *inpMesh(getMesh()); + int nbCells(inpMesh->getNumberOfCells()); + const MEDCouplingFieldDiscretization *disc(getDiscretization()); + const MEDCouplingFieldDiscretizationGauss *disc2(dynamic_cast(disc)); + if(!disc2) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::voronoize2D : Not a ON_GAUSS_PT field"); + int nbLocs(disc2->getNbOfGaussLocalization()); + std::vector< MCAuto > cells(nbCells); + for(int i=0;igetGaussLocalization(i)); + if(gl.getDimension()!=2) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::voronoize2D : not a 2D one !"); + MCAuto mesh(gl.buildRefCell()); + const std::vector& coo(gl.getGaussCoords()); + MCAuto coo2(DataArrayDouble::NewFromStdVector(coo)); + coo2->rearrange(2); + // + MCAuto coo3(MEDCouplingUMesh::Build0DMeshFromCoords(coo2)); + // + MCAuto vorCellsForCurDisc(Voronoize2D(mesh,coo2,eps)); + std::vector ids; + MCAuto ptsInReal; + disc2->getCellIdsHavingGaussLocalization(i,ids); + { + MCAuto tmp4(inpMesh->buildUnstructured()); + MCAuto subMesh(tmp4->buildPartOfMySelf(&ids[0],&ids[0]+ids.size())); + ptsInReal=gl.localizePtsInRefCooForEachCell(vorCellsForCurDisc->getCoords(),subMesh); + } + int nbPtsPerCell(vorCellsForCurDisc->getNumberOfNodes()); + for(std::size_t i=0;i elt(vorCellsForCurDisc->clone(false)); + MCAuto coo(ptsInReal->selectByTupleIdSafeSlice(i*nbPtsPerCell,(i+1)*nbPtsPerCell,1)); + elt->setCoords(coo); + cells[ids[i]]=elt; + } + } + std::vector< const MEDCouplingUMesh * > cellsPtr(VecAutoToVecOfCstPt(cells)); + MCAuto outMesh(MEDCouplingUMesh::MergeUMeshes(cellsPtr)); + MCAuto onCells(MEDCouplingFieldDouble::New(ON_CELLS)); + onCells->setMesh(outMesh); + { + MCAuto arr(getArray()->deepCopy()); + onCells->setArray(arr); + } + onCells->setTimeUnit(getTimeUnit()); + { + int b,c; + double a(getTime(b,c)); + onCells->setTime(a,b,c); + } + onCells->setName(getName()); + return onCells; +} + MEDCouplingTimeDiscretization *MEDCouplingFieldDouble::timeDiscr() { MEDCouplingTimeDiscretizationTemplate *ret(_time_discr); @@ -3089,3 +3162,4 @@ const MEDCouplingTimeDiscretization *MEDCouplingFieldDouble::timeDiscr() const throw INTERP_KERNEL::Exception("Field Double Null invalid type of time discr !"); return retc; } + diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.hxx b/src/MEDCoupling/MEDCouplingFieldDouble.hxx index 20243324a..26d5e4802 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.hxx @@ -114,6 +114,7 @@ namespace MEDCoupling MEDCOUPLING_EXPORT bool zipConnectivity(int compType, double epsOnVals=1e-15); MEDCOUPLING_EXPORT MEDCouplingFieldDouble *extractSlice3D(const double *origin, const double *vec, double eps) const; MEDCOUPLING_EXPORT bool simplexize(int policy); + MEDCOUPLING_EXPORT MCAuto voronoize(double eps) const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *computeVectorFieldCyl(const double center[3], const double vect[3]) const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *doublyContractedProduct() const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *determinant() const; @@ -166,6 +167,7 @@ namespace MEDCoupling MEDCouplingFieldDouble(const MEDCouplingFieldTemplate& ft, TypeOfTimeDiscretization td); MEDCouplingFieldDouble(const MEDCouplingFieldDouble& other, bool deepCopy); MEDCouplingFieldDouble(NatureOfField n, MEDCouplingTimeDiscretization *td, MEDCouplingFieldDiscretization *type); + MCAuto voronoize2D(double eps) const; MEDCouplingTimeDiscretization *timeDiscr(); const MEDCouplingTimeDiscretization *timeDiscr() const; }; diff --git a/src/MEDCoupling/MEDCouplingVoronoi.cxx b/src/MEDCoupling/MEDCouplingVoronoi.cxx new file mode 100644 index 000000000..7ec022a57 --- /dev/null +++ b/src/MEDCoupling/MEDCouplingVoronoi.cxx @@ -0,0 +1,218 @@ +// Copyright (C) 2007-2017 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay (EDF R&D) + +#include "MEDCouplingVoronoi.hxx" +#include "MEDCoupling1GTUMesh.hxx" +#include "MEDCouplingCMesh.hxx" +#include "MCAuto.txx" + +using namespace MEDCoupling; + +MCAuto ComputeBigCellFrom(const double pt1[2], const double pt2[2], const std::vector& bbox, double eps) +{ + static const double FACT=1.2; + MCAuto m(MEDCouplingCMesh::New()); + MCAuto arr1(DataArrayDouble::New()); arr1->alloc(2,1) ; arr1->setIJ(0,0,bbox[0]); arr1->setIJ(1,0,bbox[1]); + MCAuto arr2(DataArrayDouble::New()); arr2->alloc(2,1) ; arr2->setIJ(0,0,bbox[2]); arr2->setIJ(1,0,bbox[3]); + m->setCoords(arr1,arr2); + static const double PT[2]={0.,0.}; + m->scale(PT,FACT); + MCAuto mu(m->buildUnstructured()); + double l(std::max(bbox[1]-bbox[0],bbox[3]-bbox[2])); + double middle[2]={(pt1[0]+pt2[0])/2.,(pt1[1]+pt2[1])/2.}; + double v[2]={pt1[0],pt1[1]}; + DataArrayDouble::Rotate2DAlg(middle,M_PI/2,1,v,v); + v[0]=middle[0]-v[0]; v[1]=middle[1]-v[1]; + { + double nor(sqrt(v[0]*v[0]+v[1]*v[1])); + v[0]/=nor; v[1]/=nor; + } + MCAuto line(MEDCouplingUMesh::New("line",1)); + { + MCAuto coo(DataArrayDouble::New()); coo->alloc(2,2); + 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]); + line->setCoords(coo); + } + line->allocateCells(); + static const int CONN[2]={0,1}; + line->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,CONN); + MCAuto sp2,sp1; + { + DataArrayInt *cellNb1(0),*cellNb2(0); + MEDCouplingUMesh *sp2Pt(0),*sp1Pt(0); + MEDCouplingUMesh::Intersect2DMeshWith1DLine(mu,line,eps,sp2Pt,sp1Pt,cellNb1,cellNb2); + sp1=sp1Pt; sp2=sp2Pt; + MCAuto cellNb10(cellNb1),cellNb20(cellNb2); + } + std::vector ccp; + sp2->getCellsContainingPoint(pt1,eps,ccp); + if(ccp.size()!=1) + throw INTERP_KERNEL::Exception("ComputeBigCellFrom : expected single element !"); + MCAuto ret(sp2->buildPartOfMySelfSlice(ccp[0],ccp[0]+1,1,true)); + ret->zipCoords(); + return ret; +} + +MCAuto MergeVorCells(const std::vector< MCAuto >& vcs, double eps) +{ + std::size_t sz(vcs.size()); + if(sz<1) + throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !"); + if(sz==1) + return vcs[0]; + MCAuto p; + { + std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs)); + p=MEDCouplingUMesh::MergeUMeshes(vcsBis); + } + p->zipCoords(); + { + bool dummy; int dummy2; + MCAuto dummy3(p->mergeNodes(eps,dummy,dummy2)); + } + MCAuto edgeToKeep; + MCAuto p0; + { + MCAuto d(DataArrayInt::New()),di(DataArrayInt::New()),rd(DataArrayInt::New()),rdi(DataArrayInt::New()); + p0=p->buildDescendingConnectivity(d,di,rd,rdi); + MCAuto dsi(rdi->deltaShiftIndex()); + edgeToKeep=dsi->findIdsEqual(1); + } + MCAuto skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end())); + skinOfRes->zipCoords(); + if(skinOfRes->getNumberOfCells()!=skinOfRes->getNumberOfNodes()) + throw INTERP_KERNEL::Exception("MergeVorCells : result of merge looks bad !"); + MCAuto d(skinOfRes->orderConsecutiveCells1D()); + MCAuto skinOfRes2; + { + MCAuto part(skinOfRes->buildPartOfMySelf(d->begin(),d->end())); + skinOfRes2=MEDCoupling1SGTUMesh::New(part); + } + MCAuto c(skinOfRes2->getNodalConnectivity()->deepCopy()); + c->circularPermutation(1); + c->rearrange(2); + std::vector< MCAuto > vdi(c->explodeComponents()); + if(!vdi[0]->isEqual(*vdi[1])) + throw INTERP_KERNEL::Exception("MergeVorCells : internal error !"); + MCAuto m(MEDCouplingUMesh::New("",2)); + m->setCoords(skinOfRes2->getCoords()); + m->allocateCells(); + m->insertNextCell(INTERP_KERNEL::NORM_POLYGON,vdi[0]->getNumberOfTuples(),vdi[0]->begin()); + return m; +} + +MCAuto MEDCoupling::Voronoize2D(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) +{ + if(!m || !points) + throw INTERP_KERNEL::Exception("Voronoize2D : null pointer !"); + m->checkConsistencyLight(); + points->checkAllocated(); + if(m->getMeshDimension()!=2 || m->getSpaceDimension()!=2 || points->getNumberOfComponents()!=2) + throw INTERP_KERNEL::Exception("Voronoize2D : spacedim must be equal to 2 and meshdim also equal to 2 !"); + if(m->getNumberOfCells()!=1) + throw INTERP_KERNEL::Exception("Voronoize2D : mesh is expected to have only one cell !"); + int nbPts(points->getNumberOfTuples()); + if(nbPts<1) + throw INTERP_KERNEL::Exception("Voronoize2D : at least one point expected !"); + std::vector bbox(4); + m->getBoundingBox(&bbox[0]); + std::vector< MCAuto > l0(1,MCAuto(m->deepCopy())); + const double *pts(points->begin()); + for(int i=1;i vorTess; + { + std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0)); + vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis); + } + { + bool dummy; + int newNbNodes; + MCAuto dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes)); + } + std::vector polygsToIterOn; + const double *pt(pts+i*2); + vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn); + if(polygsToIterOn.size()<1) + throw INTERP_KERNEL::Exception("Voronoize2D : presence of a point outside the given cell !"); + std::set elemsToDo,elemsDone; elemsToDo.insert(polygsToIterOn[0]); + std::vector< MCAuto > newVorCells; + while(!elemsToDo.empty()) + { + int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly); + const double *seed(pts+2*poly); + MCAuto cell(ComputeBigCellFrom(pt,seed,bbox,eps)); + MCAuto tile(l0[poly]); + tile->zipCoords(); + MCAuto a; + MCAuto b,c; + { + DataArrayInt *bPtr(0),*cPtr(0); + a=MEDCouplingUMesh::Intersect2DMeshes(tile,cell,eps,bPtr,cPtr); + b=bPtr; c=cPtr; + } + MCAuto part(c->findIdsEqual(-1)); + if(part->getNumberOfTuples()!=1) + throw INTERP_KERNEL::Exception("Voronoize2D : internal error"); + MCAuto newVorCell; + { + MCAuto tmp(part->buildComplement(a->getNumberOfCells())); + newVorCell=a->buildPartOfMySelf(tmp->begin(),tmp->end()); + } + newVorCell->zipCoords(); + MCAuto modifiedCell(a->buildPartOfMySelf(part->begin(),part->end())); + modifiedCell->zipCoords(); + l0[poly]=modifiedCell; + // + MCAuto ids; + { + DataArrayInt *tmp(0); + bool sta(a->getCoords()->areIncludedInMe(cell->getCoords(),eps,tmp)); + ids=tmp; + if(!sta) + throw INTERP_KERNEL::Exception("Voronoize2D : internal error 2 !"); + } + MCAuto newCoords; + { + MCAuto tmp(ids->buildComplement(a->getNumberOfNodes())); + newCoords=a->getCoords()->selectByTupleId(tmp->begin(),tmp->end()); + } + const double *cPtr(newCoords->begin()); + for(int i=0;igetNumberOfTuples();i++,cPtr+=2) + { + std::set zeCandidates; + { + std::vector zeCandidatesTmp; + vorTess->getCellsContainingPoint(cPtr,eps,zeCandidatesTmp); + zeCandidates.insert(zeCandidatesTmp.begin(),zeCandidatesTmp.end()); + } + std::set tmp,newElementsToDo; + std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp,tmp.begin())); + std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp.begin(),tmp.end(),std::inserter(newElementsToDo,newElementsToDo.begin())); + elemsToDo=newElementsToDo; + } + newVorCells.push_back(newVorCell); + } + l0.push_back(MergeVorCells(newVorCells,eps)); + } + std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0)); + MCAuto ret(MEDCouplingUMesh::MergeUMeshes(l0Bis)); + return ret; +} diff --git a/src/MEDCoupling/MEDCouplingVoronoi.hxx b/src/MEDCoupling/MEDCouplingVoronoi.hxx new file mode 100644 index 000000000..06dc32180 --- /dev/null +++ b/src/MEDCoupling/MEDCouplingVoronoi.hxx @@ -0,0 +1,33 @@ +// Copyright (C) 2007-2017 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay (EDF R&D) + +#ifndef __MEDCOUPLINGVORONOI_HXX__ +#define __MEDCOUPLINGVORONOI_HXX__ + +#include "MEDCoupling.hxx" + +#include "MEDCouplingUMesh.hxx" + +namespace MEDCoupling +{ + MCAuto Voronoize2D(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps); +} + +#endif diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index d4ad97677..5f7c1943e 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -212,6 +212,7 @@ using namespace INTERP_KERNEL; %newobject MEDCoupling::MEDCouplingFieldDouble::findIdsInRange; %newobject MEDCoupling::MEDCouplingFieldDouble::buildSubPart; %newobject MEDCoupling::MEDCouplingFieldDouble::buildSubPartRange; +%newobject MEDCoupling::MEDCouplingFieldDouble::voronoize; %newobject MEDCoupling::MEDCouplingFieldDouble::__getitem__; %newobject MEDCoupling::MEDCouplingFieldDouble::__neg__; %newobject MEDCoupling::MEDCouplingFieldDouble::__add__; @@ -3937,7 +3938,13 @@ namespace MEDCoupling self->reprQuickOverview(oss); return oss.str(); } - + + MEDCouplingFieldDouble *voronoize(double eps) const throw(INTERP_KERNEL::Exception) + { + MCAuto ret(self->voronoize(eps)); + return ret.retn(); + } + MEDCouplingFieldDouble *computeVectorFieldCyl(PyObject *center, PyObject *vector) const { const char msg[]="Python wrap of MEDCouplingFieldDouble::computeVectorFieldCyl : ";