X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FMEDCoupling%2FMEDCouplingUMesh.cxx;h=7a4aae408c4ec5742cdfd77c61269c1eab2ead3c;hb=659f8c67d0348350e12fde38fe8c4de1ff95dffe;hp=2ba48b592b804217f45e6ad985a10a9790d98f96;hpb=18e56d429e218d5503e4142fe5cde99c80106386;p=tools%2Fmedcoupling.git diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 2ba48b592..7a4aae408 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -1,9 +1,9 @@ -// Copyright (C) 2007-2013 CEA/DEN, EDF R&D +// Copyright (C) 2007-2014 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. +// 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 @@ -58,7 +58,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::New() return new MEDCouplingUMesh; } -MEDCouplingUMesh *MEDCouplingUMesh::New(const char *meshName, int meshDim) +MEDCouplingUMesh *MEDCouplingUMesh::New(const std::string& meshName, int meshDim) { MEDCouplingUMesh *ret=new MEDCouplingUMesh; ret->setName(meshName); @@ -419,12 +419,39 @@ MEDCouplingUMeshCellByTypeEntry *MEDCouplingUMesh::cellsByType() * Returns a set of all cell types available in \a this mesh. * \return std::set - the set of cell types. * \warning this method does not throw any exception even if \a this is not defined. + * \sa MEDCouplingUMesh::getAllGeoTypesSorted */ std::set MEDCouplingUMesh::getAllGeoTypes() const { return _types; } +/*! + * This method returns the sorted list of geometric types in \a this. + * Sorted means in the same order than the cells in \a this. A single entry in return vector means the maximal chunk of consecutive cells in \a this + * having the same geometric type. So a same geometric type can appear more than once if the cells are not sorted per geometric type. + * + * \throw if connectivity in \a this is not correctly defined. + * + * \sa MEDCouplingMesh::getAllGeoTypes + */ +std::vector MEDCouplingUMesh::getAllGeoTypesSorted() const +{ + std::vector ret; + checkConnectivityFullyDefined(); + int nbOfCells(getNumberOfCells()); + if(nbOfCells==0) + return ret; + if(getMeshLength()<1) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getAllGeoTypesSorted : the connectivity in this seems invalid !"); + const int *c(_nodal_connec->begin()),*ci(_nodal_connec_index->begin()); + ret.push_back((INTERP_KERNEL::NormalizedCellType)c[*ci++]); + for(int i=1;igetConstPointer(); const int *connIndex=_nodal_connec_index->getConstPointer(); std::string name="Mesh constituent of "; name+=getName(); - MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingUMesh::New(name.c_str(),getMeshDimension()-SonsGenerator::DELTA); + MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingUMesh::New(name,getMeshDimension()-SonsGenerator::DELTA); ret->setCoords(getCoords()); ret->allocateCells(2*nbOfCells); descIndx->alloc(nbOfCells+1,1); @@ -1822,7 +1849,7 @@ bool MEDCouplingUMesh::areCellsIncludedIn(const MEDCouplingUMesh *other, int com } MEDCouplingAutoRefCountObjectPtr o2n=mesh->zipConnectivityTraducer(compType,nbOfCells); arr=o2n->substr(nbOfCells); - arr->setName(other->getName().c_str()); + arr->setName(other->getName()); int tmp; if(other->getNumberOfCells()==0) return true; @@ -1866,7 +1893,7 @@ bool MEDCouplingUMesh::areCellsIncludedIn2(const MEDCouplingUMesh *other, DataAr } } } - arr2->setName(other->getName().c_str()); + arr2->setName(other->getName()); if(arr2->presenceOfValue(0)) return false; arr=arr2.retn(); @@ -2836,7 +2863,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildSetInstanceFromThis(int spaceDim) const int mdim=getMeshDimension(); if(mdim<0) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSetInstanceFromThis : invalid mesh dimension ! Should be >= 0 !"); - MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingUMesh::New(getName().c_str(),mdim); + MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingUMesh::New(getName(),mdim); MEDCouplingAutoRefCountObjectPtr tmp1,tmp2; bool needToCpyCT=true; if(!_nodal_connec) @@ -2910,7 +2937,7 @@ int MEDCouplingUMesh::getNumberOfNodesInCell(int cellId) const * describing the cell types. * \throw If the coordinates array is not set. * \throw If the nodal connectivity of cells is not defined. - * \sa getAllTypes() + * \sa getAllGeoTypes() */ std::set MEDCouplingUMesh::getTypesOfPart(const int *begin, const int *end) const { @@ -3224,7 +3251,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField(bool isAbs) const name+=getName(); int nbelem=getNumberOfCells(); MEDCouplingAutoRefCountObjectPtr field=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME); - field->setName(name.c_str()); + field->setName(name); MEDCouplingAutoRefCountObjectPtr array=DataArrayDouble::New(); array->alloc(nbelem,1); double *area_vol=array->getPointer(); @@ -3279,7 +3306,7 @@ DataArrayDouble *MEDCouplingUMesh::getPartMeasureField(bool isAbs, const int *be name+=getName(); int nbelem=(int)std::distance(begin,end); MEDCouplingAutoRefCountObjectPtr array=DataArrayDouble::New(); - array->setName(name.c_str()); + array->setName(name); array->alloc(nbelem,1); double *area_vol=array->getPointer(); if(getMeshDimension()!=-1) @@ -3994,6 +4021,9 @@ void MEDCouplingUMesh::DistanceToPoint2DCurveAlg(const double *pt, const int *ce /*! * Finds cells in contact with a ball (i.e. a point with precision). + * For speed reasons, the INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6 and INTERP_KERNEL::NORM_QUAD8 cells are considered as convex cells to detect if a point is IN or OUT. + * If it is not the case, please change their types to INTERP_KERNEL::NORM_POLYGON or INTERP_KERNEL::NORM_QPOLYG before invoking this method. + * * \warning This method is suitable if the caller intends to evaluate only one * point, for more points getCellsContainingPoints() is recommended as it is * faster. @@ -4015,6 +4045,8 @@ int MEDCouplingUMesh::getCellContainingPoint(const double *pos, double eps) cons /*! * Finds cells in contact with a ball (i.e. a point with precision). + * For speed reasons, the INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6 and INTERP_KERNEL::NORM_QUAD8 cells are considered as convex cells to detect if a point is IN or OUT. + * If it is not the case, please change their types to INTERP_KERNEL::NORM_POLYGON or INTERP_KERNEL::NORM_QPOLYG before invoking this method. * \warning This method is suitable if the caller intends to evaluate only one * point, for more points getCellsContainingPoints() is recommended as it is * faster. @@ -4071,6 +4103,7 @@ namespace ParaMEDMEM INTERP_KERNEL::EdgeLin *e1=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[2]].first); INTERP_KERNEL::EdgeLin *e2=new INTERP_KERNEL::EdgeLin(mapp2[bg[2]].first,mapp2[bg[1]].first); INTERP_KERNEL::SegSegIntersector inters(*e1,*e2); + // is the SEG3 degenerated, and thus can be reduced to a SEG2? bool colinearity=inters.areColinears(); delete e1; delete e2; if(colinearity) @@ -4087,11 +4120,14 @@ namespace ParaMEDMEM } /*! - * This method creates a sub mesh in Geometric2D DS. The sub mesh is composed be the sub set of cells in 'candidates' and the global mesh 'mDesc'. - * The input meth 'mDesc' must be so that mDim==1 et spaceDim==3. - * 'mapp' contains a mapping between local numbering in submesh and the global node numbering in 'mDesc'. + * This method creates a sub mesh in Geometric2D DS. The sub mesh is composed by the sub set of cells in 'candidates' taken from + * the global mesh 'mDesc'. + * The input mesh 'mDesc' must be so that mDim==1 and spaceDim==2. + * 'mapp' returns a mapping between local numbering in submesh (represented by a Node*) and the global node numbering in 'mDesc'. */ - INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMesh(const MEDCouplingUMesh *mDesc, const std::vector& candidates, std::map& mapp) throw(INTERP_KERNEL::Exception) + INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMesh(const MEDCouplingUMesh *mDesc, const std::vector& candidates, + std::map& mapp) + throw(INTERP_KERNEL::Exception) { mapp.clear(); std::map > mapp2;//bool is for a flag specifying if node is boundary (true) or only a middle for SEG3. @@ -4136,6 +4172,9 @@ namespace ParaMEDMEM return new INTERP_KERNEL::Node(coo1[2*nodeId],coo1[2*nodeId+1]); } + /** + * Construct a mapping between set of Nodes and the standart MEDCoupling connectivity format (c, cI). + */ void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector& addCoo, const int *desc1Bg, const int *desc1End, const std::vector >& intesctEdges1, /*output*/std::map& mapp, std::map& mappRev) @@ -4165,7 +4204,7 @@ void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const d { elts=DataArrayInt::New(); eltsIndex=DataArrayInt::New(); eltsIndex->alloc(nbOfPoints+1,1); eltsIndex->setIJ(0,0,0); elts->alloc(0,1); int *eltsIndexPtr(eltsIndex->getPointer()); - MEDCouplingAutoRefCountObjectPtr bboxArr(getBoundingBoxForBBTree()); + MEDCouplingAutoRefCountObjectPtr bboxArr(getBoundingBoxForBBTree(eps)); const double *bbox(bboxArr->begin()); int nbOfCells=getNumberOfCells(); const int *conn=_nodal_connec->getConstPointer(); @@ -4184,10 +4223,35 @@ void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const d myTree.getIntersectingElems(bb,candidates); for(std::vector::const_iterator iter=candidates.begin();iter!=candidates.end();iter++) { - int sz=connI[(*iter)+1]-connI[*iter]-1; - if(INTERP_KERNEL::PointLocatorAlgos >::isElementContainsPoint(pos+i*SPACEDIM, - (INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]], - coords,conn+connI[*iter]+1,sz,eps)) + int sz(connI[(*iter)+1]-connI[*iter]-1); + INTERP_KERNEL::NormalizedCellType ct((INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]]); + bool status(false); + if(ct!=INTERP_KERNEL::NORM_POLYGON && ct!=INTERP_KERNEL::NORM_QPOLYG) + status=INTERP_KERNEL::PointLocatorAlgos >::isElementContainsPoint(pos+i*SPACEDIM,ct,coords,conn+connI[*iter]+1,sz,eps); + else + { + if(SPACEDIM!=2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellsContainingPointsAlg : not implemented yet for POLYGON and QPOLYGON in spaceDim 3 !"); + INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps; + INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps; + std::vector nodes(sz); + INTERP_KERNEL::QuadraticPolygon *pol(0); + for(int j=0;jnormalizeMe(b,c); n->applySimilarity(b,c,a); + status=pol->isInOrOut2(n); + delete pol; n->decrRef(); + } + if(status) { eltsIndexPtr[i+1]++; elts->pushBackSilent(*iter); @@ -4199,6 +4263,8 @@ void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const d * Finds cells in contact with several balls (i.e. points with precision). * This method is an extension of getCellContainingPoint() and * getCellsContainingPoint() for the case of multiple points. + * For speed reasons, the INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6 and INTERP_KERNEL::NORM_QUAD8 cells are considered as convex cells to detect if a point is IN or OUT. + * If it is not the case, please change their types to INTERP_KERNEL::NORM_POLYGON or INTERP_KERNEL::NORM_QPOLYG before invoking this method. * \param [in] pos - an array of coordinates of points in full interlace mode : * X0,Y0,Z0,X1,Y1,Z1,... Size of the array must be \a * this->getSpaceDimension() * \a nbOfPoints @@ -5245,7 +5311,7 @@ void MEDCouplingUMesh::tessellate2DCurve(double eps) * and \a this->getMeshDimension() != 3. * \throw If \a policy is not one of the four discussed above. * \throw If the nodal connectivity of cells is not defined. - * \sa MEDCouplingUMesh::tetrahedrize + * \sa MEDCouplingUMesh::tetrahedrize, MEDCoupling1SGTUMesh::sortHexa8EachOther */ DataArrayInt *MEDCouplingUMesh::simplexize(int policy) { @@ -6195,12 +6261,45 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getSkewField() const /*! * This method aggregate the bbox of each cell and put it into bbox parameter. * + * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12) + * For all other cases this input parameter is ignored. * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components. * * \throw If \a this is not fully set (coordinates and connectivity). * \throw If a cell in \a this has no valid nodeId. + * \sa MEDCouplingUMesh::getBoundingBoxForBBTreeFast, MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic */ -DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree() const +DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree(double arcDetEps) const +{ + int mDim(getMeshDimension()),sDim(getSpaceDimension()); + if(mDim!=2 || sDim!=2) + return getBoundingBoxForBBTreeFast(); + else + { + bool presenceOfQuadratic(false); + for(std::set::const_iterator it=_types.begin();it!=_types.end();it++) + { + const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(*it)); + if(cm.isQuadratic()) + presenceOfQuadratic=true; + } + if(presenceOfQuadratic) + return getBoundingBoxForBBTree2DQuadratic(arcDetEps); + else + return getBoundingBoxForBBTreeFast(); + } +} + +/*! + * This method aggregate the bbox of each cell only considering the nodes constituting each cell and put it into bbox parameter. + * So meshes having quadratic cells the computed bounding boxes can be invalid ! + * + * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components. + * + * \throw If \a this is not fully set (coordinates and connectivity). + * \throw If a cell in \a this has no valid nodeId. + */ +DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTreeFast() const { checkFullyDefined(); int spaceDim(getSpaceDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes()); @@ -6239,6 +6338,48 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree() const return ret.retn(); } +/*! + * This method aggregate the bbox regarding foreach 2D cell in \a this the whole shape. So this method is particulary useful for 2D meshes having quadratic cells + * because for this type of cells getBoundingBoxForBBTreeFast method may return invalid bounding boxes. + * + * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12) + * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components. + * \throw If \a this is not fully defined. + * \throw If \a this is not a mesh with meshDimension equal to 2. + * \throw If \a this is not a mesh with spaceDimension equal to 2. + */ +DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arcDetEps) const +{ + checkFullyDefined(); + int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes()); + if(mDim!=2 || spaceDim!=2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic : This method should be applied on mesh with mesh dimension equal to 2 and space dimension also equal to 2!"); + MEDCouplingAutoRefCountObjectPtr ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim); + double *bbox(ret->getPointer()); + const double *coords(_coords->getConstPointer()); + const int *conn(_nodal_connec->getConstPointer()),*connI(_nodal_connec_index->getConstPointer()); + for(int i=0;i nodes(sz); + INTERP_KERNEL::QuadraticPolygon *pol(0); + for(int j=0;jfillBounds(b); delete pol; + bbox[0]=b.getXMin(); bbox[1]=b.getXMax(); bbox[2]=b.getYMin(); bbox[3]=b.getYMax(); + } + return ret.retn(); +} + /// @cond INTERNAL namespace ParaMEDMEMImpl @@ -6281,7 +6422,7 @@ std::vector MEDCouplingUMesh::getDistributionOfTypes() const const int *connI=_nodal_connec_index->getConstPointer(); const int *work=connI; int nbOfCells=getNumberOfCells(); - std::size_t n=getAllTypes().size(); + std::size_t n=getAllGeoTypes().size(); std::vector ret(3*n,-1); //ret[3*k+2]==-1 because it has no sense here std::set types; for(std::size_t i=0;work!=connI+nbOfCells;i++) @@ -6772,7 +6913,7 @@ MEDCoupling1GTUMesh *MEDCouplingUMesh::convertIntoSingleGeoTypeMesh() const if(_types.size()!=1) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convertIntoSingleGeoTypeMesh : current mesh does not contain exactly one geometric type !"); INTERP_KERNEL::NormalizedCellType typ=*_types.begin(); - MEDCouplingAutoRefCountObjectPtr ret=MEDCoupling1GTUMesh::New(getName().c_str(),typ); + MEDCouplingAutoRefCountObjectPtr ret=MEDCoupling1GTUMesh::New(getName(),typ); ret->setCoords(getCoords()); MEDCoupling1SGTUMesh *retC=dynamic_cast((MEDCoupling1GTUMesh*)ret); if(retC) @@ -6952,7 +7093,7 @@ DataArrayInt *MEDCouplingUMesh::convertCellArrayPerGeoType(const DataArrayInt *d const int *conn=_nodal_connec->getConstPointer(); const int *connI=_nodal_connec_index->getConstPointer(); int nbOfCells=getNumberOfCells(); - std::set types=getAllTypes(); + std::set types(getAllGeoTypes()); int *tmp=new int[nbOfCells]; for(std::set::const_iterator iter=types.begin();iter!=types.end();iter++) { @@ -7248,7 +7389,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::Build0DMeshFromCoords(DataArrayDouble *da) if(!da) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Build0DMeshFromCoords : instance of DataArrayDouble must be not null !"); da->checkAllocated(); - MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingUMesh::New(da->getName().c_str(),0); + MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingUMesh::New(da->getName(),0); ret->setCoords(da); int nbOfTuples=da->getNumberOfTuples(); MEDCouplingAutoRefCountObjectPtr c=DataArrayInt::New(); @@ -7521,7 +7662,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::FuseUMeshesOnSameCoords(const std::vectoralloc(curNbOfCells,1); std::copy(o2nPtr+offset,o2nPtr+offset+curNbOfCells,tmp->getPointer()); offset+=curNbOfCells; - tmp->setName(meshes[i]->getName().c_str()); + tmp->setName(meshes[i]->getName()); corr[i]=tmp; } return ret.retn(); @@ -8143,7 +8284,7 @@ void MEDCouplingUMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData int nbOfCells=getNumberOfCells(); if(nbOfCells<=0) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::writeVTK : the unstructured mesh has no cells !"); - static const int PARAMEDMEM2VTKTYPETRADUCER[INTERP_KERNEL::NORM_MAXTYPE+1]={1,3,21,5,9,7,22,34,23,28,-1,-1,-1,-1,10,14,13,-1,12,-1,24,-1,16,27,-1,26,-1,29,-1,-1,25,42,-1,4}; + static const int PARAMEDMEM2VTKTYPETRADUCER[INTERP_KERNEL::NORM_MAXTYPE+1]={1,3,21,5,9,7,22,34,23,28,-1,-1,-1,-1,10,14,13,-1,12,-1,24,-1,16,27,-1,26,-1,29,-1,-1,25,42,36,4}; ofs << " <" << getVTKDataSetType() << ">\n"; ofs << " \n"; ofs << " \n" << pointData << std::endl; @@ -8249,7 +8390,10 @@ std::string MEDCouplingUMesh::getVTKDataSetType() const /*! * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and - * returns a result mesh constituted by polygons. The meshes should be in 2D space. In + * returns a result mesh constituted by polygons. + * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains + * all nodes from m2. + * The meshes should be in 2D space. In * addition, returns two arrays mapping cells of the result mesh to cells of the input * meshes. * \param [in] m1 - the first input mesh which is a partitioned object. @@ -8270,31 +8414,40 @@ std::string MEDCouplingUMesh::getVTKDataSetType() const * \throw If the nodal connectivity of cells is not defined in any of the meshes. * \throw If any of the meshes is not a 2D mesh in 2D space. */ -MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2) +MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, + double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2) { m1->checkFullyDefined(); m2->checkFullyDefined(); if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!"); + + // Step 1: compute all edge intersections (new nodes) std::vector< std::vector > intersectEdge1, colinear2, subDiv2; - MEDCouplingUMesh *m1Desc=0,*m2Desc=0; + MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0; - std::vector addCoo,addCoordsQuadratic; + std::vector addCoo,addCoordsQuadratic; // coordinates of newly created nodes INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps; INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps; - IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,m1Desc,desc1,descIndx1,revDesc1,revDescIndx1, - m2Desc,desc2,descIndx2,revDesc2,revDescIndx2,addCoo); + IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2, + m1Desc,desc1,descIndx1,revDesc1,revDescIndx1, + addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2); revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef(); MEDCouplingAutoRefCountObjectPtr dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2); MEDCouplingAutoRefCountObjectPtr dd5(m1Desc),dd6(m2Desc); + + // Step 2: re-order newly created nodes according to the ordering found in m2 std::vector< std::vector > intersectEdge2; BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2); subDiv2.clear(); dd5=0; dd6=0; + + // Step 3: std::vector cr,crI; //no DataArrayInt because interface with Geometric2D std::vector cNb1,cNb2; //no DataArrayInt because interface with Geometric2D BuildIntersecting2DCellsFromEdges(eps,m1,desc1->getConstPointer(),descIndx1->getConstPointer(),intersectEdge1,colinear2,m2,desc2->getConstPointer(),descIndx2->getConstPointer(),intersectEdge2,addCoo, /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2); - // + + // Step 4: Prepare final result: MEDCouplingAutoRefCountObjectPtr addCooDa=DataArrayDouble::New(); addCooDa->alloc((int)(addCoo.size())/2,2); std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer()); @@ -8315,6 +8468,15 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1 return ret.retn(); } + +/** + * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the + * (newly created) nodes corresponding to the edge intersections. + * Output params: + * @param[out] cr, crI connectivity of the resulting mesh + * @param[out] cNb1, cNb2 correspondance arrays giving for the merged mesh the initial cells IDs in m1 / m2 + * TODO: describe input parameters + */ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1, const std::vector >& intesctEdges1, const std::vector< std::vector >& colinear2, const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector >& intesctEdges2, @@ -8333,6 +8495,7 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo int offset3=offset2+((int)addCoords.size())/2; MEDCouplingAutoRefCountObjectPtr bbox1Arr(m1->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree()); const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin()); + // Here a BBTree on 2D-cells, not on segments: BBTree myTree(bbox2,0,0,m2->getNumberOfCells(),eps); int ncell1=m1->getNumberOfCells(); crI.push_back(0); @@ -8345,7 +8508,9 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo INTERP_KERNEL::QuadraticPolygon pol1; INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]]; const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ); + // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects: MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev); + // pol1 is the full cell from mesh2, in QP format, with all the additional intersecting nodes. pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1, desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1); // @@ -8355,16 +8520,18 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo for(it1.first();!it1.finished();it1.next()) edges1.insert(it1.current()->getPtr()); // - std::map > edgesIn2ForShare; + std::map > edgesIn2ForShare; // common edges std::vector pol2s(candidates2.size()); int ii=0; for(std::vector::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++) { INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]]; const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2); + // Complete mapping with elements coming from the current cell it2 in mesh2: MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev); + // pol2 is the new QP in the final merged result. pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2, - pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2,edgesIn2ForShare); + pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare); } ii=0; for(std::vector::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++) @@ -8374,6 +8541,8 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2); pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2); } + // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched + // by m2 but that we still want to keep in the final result. if(!edges1.empty()) { try @@ -8393,15 +8562,19 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo /*! * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2). - * + * It builds the descending connectivity of the two meshes, and then using a binary tree + * it computes the edge intersections. This results in new points being created : they're stored in addCoo. + * Documentation about parameters colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs(). */ void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, std::vector< std::vector >& intersectEdge1, std::vector< std::vector >& colinear2, std::vector< std::vector >& subDiv2, MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1, - MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2, - std::vector& addCoo) throw(INTERP_KERNEL::Exception) + std::vector& addCoo, + MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2) + throw(INTERP_KERNEL::Exception) { static const int SPACEDIM=2; + // Build desc connectivity desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New(); desc2=DataArrayInt::New(); descIndx2=DataArrayInt::New(); @@ -8414,29 +8587,33 @@ void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, c MEDCouplingAutoRefCountObjectPtr dd9(m1Desc),dd10(m2Desc); const int *c1=m1Desc->getNodalConnectivity()->getConstPointer(); const int *ci1=m1Desc->getNodalConnectivityIndex()->getConstPointer(); + + // Build BB tree of all edges in the tool mesh (second mesh) MEDCouplingAutoRefCountObjectPtr bbox1Arr(m1Desc->getBoundingBoxForBBTree()),bbox2Arr(m2Desc->getBoundingBoxForBBTree()); const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin()); - int ncell1=m1Desc->getNumberOfCells(); - int ncell2=m2Desc->getNumberOfCells(); - intersectEdge1.resize(ncell1); - colinear2.resize(ncell2); - subDiv2.resize(ncell2); + int nDescCell1=m1Desc->getNumberOfCells(); + int nDescCell2=m2Desc->getNumberOfCells(); + intersectEdge1.resize(nDescCell1); + colinear2.resize(nDescCell2); + subDiv2.resize(nDescCell2); BBTree myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps); + std::vector candidates1(1); int offset1=m1->getNumberOfNodes(); int offset2=offset1+m2->getNumberOfNodes(); - for(int i=0;i candidates2; + std::vector candidates2; // edges of mesh2 candidate for intersection myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2); - if(!candidates2.empty()) + if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1 { std::map map1,map2; + // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2); candidates1[0]=i; INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1); - // this following part is to avoid that a some remove nodes (for example due to a merge between pol1 and pol2) can be replaced by a newlt created one - // This trick garanties that Node * are discriminant + // This following part is to avoid that some removed nodes (for example due to a merge between pol1 and pol2) are replaced by a newly created one + // This trick guarantees that Node * are discriminant (i.e. form a unique identifier) std::set nodes; pol1->getAllNodes(nodes); pol2->getAllNodes(nodes); std::size_t szz(nodes.size()); @@ -8445,6 +8622,7 @@ void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, c for(std::size_t iii=0;iiiincrRef(); nodesSafe[iii]=*itt; } // end of protection + // Performs egde cutting: pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo); delete pol2; delete pol1; @@ -8461,13 +8639,18 @@ void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, c * This method has 4 inputs : * - a mesh 'm1' with meshDim==1 and a SpaceDim==2 * - a mesh 'm2' with meshDim==1 and a SpaceDim==2 - * - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm' the splitting node ids in randomly sorted. - * The aim of this method is to sort the splitting nodes, if any, and to put in 'intersectEdge' output paramter based on edges of mesh 'm2' - * \param m1 is expected to be a mesh of meshDimension equal to 1 and spaceDim equal to 2. No check of that is performed by this method. Only present for its coords in case of 'subDiv' shares some nodes of 'm1' + * - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm' the splitting node ids randomly sorted. + * The aim of this method is to sort the splitting nodes, if any, and to put them in 'intersectEdge' output parameter based on edges of mesh 'm2' + * Nodes end up lying consecutively on a cutted edge. + * \param m1 is expected to be a mesh of meshDimension equal to 1 and spaceDim equal to 2. No check of that is performed by this method. + * (Only present for its coords in case of 'subDiv' shares some nodes of 'm1') * \param m2 is expected to be a mesh of meshDimension equal to 1 and spaceDim equal to 2. No check of that is performed by this method. - * \param addCoo input parameter with additionnal nodes linked to intersection of the 2 meshes. + * \param addCoo input parameter with additional nodes linked to intersection of the 2 meshes. + * \param[out] intersectEdge the same content as subDiv, but correclty oriented. */ -void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, const std::vector& addCoo, const std::vector< std::vector >& subDiv, std::vector< std::vector >& intersectEdge) +void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, + const std::vector& addCoo, + const std::vector< std::vector >& subDiv, std::vector< std::vector >& intersectEdge) { int offset1=m1->getNumberOfNodes(); int ncell=m2->getNumberOfCells(); @@ -9296,7 +9479,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildSpreadZonesWithPoly() const std::vector partition=partitionBySpreadZone(); std::vector< MEDCouplingAutoRefCountObjectPtr > partitionAuto; partitionAuto.reserve(partition.size()); std::copy(partition.begin(),partition.end(),std::back_insert_iterator > >(partitionAuto)); - MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingUMesh::New(getName().c_str(),mdim); + MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingUMesh::New(getName(),mdim); ret->setCoords(getCoords()); ret->allocateCells((int)partition.size()); // @@ -9391,10 +9574,13 @@ DataArrayInt *MEDCouplingUMesh::ComputeRangesFromTypeDistribution(const std::vec * an id of old cell producing it. The caller is to delete this array using * decrRef() as it is no more needed. * \return MEDCoupling1SGTUMesh * - the mesh containing only INTERP_KERNEL::NORM_TETRA4 cells. + * + * \warning This method operates on each cells in this independantly ! So it can leads to non conform mesh in returned value ! If you expect to have a conform mesh in output + * the policy PLANAR_FACE_6 should be used on a mesh sorted with MEDCoupling1SGTUMesh::sortHexa8EachOther. * * \throw If \a this is not a 3D mesh (spaceDim==3 and meshDim==3). * \throw If \a this is not fully constituted with linear 3D cells. - * \sa MEDCouplingUMesh::simplexize + * \sa MEDCouplingUMesh::simplexize, MEDCoupling1SGTUMesh::sortHexa8EachOther */ MEDCoupling1SGTUMesh *MEDCouplingUMesh::tetrahedrize(int policy, DataArrayInt *& n2oCells, int& nbOfAdditionalPoints) const { @@ -9403,7 +9589,7 @@ MEDCoupling1SGTUMesh *MEDCouplingUMesh::tetrahedrize(int policy, DataArrayInt *& if(getMeshDimension()!=3 || getSpaceDimension()!=3) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tetrahedrize : only available for mesh with meshdim == 3 and spacedim == 3 !"); int nbOfCells(getNumberOfCells()),nbNodes(getNumberOfNodes()); - MEDCouplingAutoRefCountObjectPtr ret0(MEDCoupling1SGTUMesh::New(getName().c_str(),INTERP_KERNEL::NORM_TETRA4)); + MEDCouplingAutoRefCountObjectPtr ret0(MEDCoupling1SGTUMesh::New(getName(),INTERP_KERNEL::NORM_TETRA4)); MEDCouplingAutoRefCountObjectPtr ret(DataArrayInt::New()); ret->alloc(nbOfCells,1); int *retPt(ret->getPointer()); MEDCouplingAutoRefCountObjectPtr newConn(DataArrayInt::New()); newConn->alloc(0,1);