X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FMEDCoupling%2FMEDCouplingUMesh.cxx;h=697b876ab41c2ccf18c9ebb0565217f4e7f06206;hb=fc22b4cd63404700f73e09be7507a11b838e55bc;hp=2757c52aec249dd9cab5e70cea4c85b534480186;hpb=c1a99376f5e8a6804d29ea0b42057ac87fd9890c;p=tools%2Fmedcoupling.git diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 2757c52ae..697b876ab 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -16,7 +16,7 @@ // // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -// Author : Anthony Geay (CEA/DEN) +// Author : Anthony Geay (EDF R&D) #include "MEDCouplingUMesh.hxx" #include "MEDCouplingCMesh.hxx" @@ -41,6 +41,7 @@ #include "InterpKernelGeo2DEdgeLin.hxx" #include "InterpKernelGeo2DEdgeArcCircle.hxx" #include "InterpKernelGeo2DQuadraticPolygon.hxx" +#include "OrientationInverter.hxx" #include "MEDCouplingUMesh_internal.hxx" #include @@ -55,8 +56,8 @@ using namespace MEDCoupling; double MEDCouplingUMesh::EPS_FOR_POLYH_ORIENTATION=1.e-14; /// @cond INTERNAL -const INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::MEDMEM_ORDER[N_MEDMEM_ORDER] = { INTERP_KERNEL::NORM_POINT1, INTERP_KERNEL::NORM_SEG2, INTERP_KERNEL::NORM_SEG3, INTERP_KERNEL::NORM_SEG4, INTERP_KERNEL::NORM_POLYL, INTERP_KERNEL::NORM_TRI3, INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_TRI7, INTERP_KERNEL::NORM_QUAD8, INTERP_KERNEL::NORM_QUAD9, INTERP_KERNEL::NORM_POLYGON, INTERP_KERNEL::NORM_QPOLYG, INTERP_KERNEL::NORM_TETRA4, INTERP_KERNEL::NORM_PYRA5, INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXA8, INTERP_KERNEL::NORM_HEXGP12, INTERP_KERNEL::NORM_TETRA10, INTERP_KERNEL::NORM_PYRA13, INTERP_KERNEL::NORM_PENTA15, INTERP_KERNEL::NORM_HEXA20, INTERP_KERNEL::NORM_HEXA27, INTERP_KERNEL::NORM_POLYHED }; -const int MEDCouplingUMesh::MEDCOUPLING2VTKTYPETRADUCER[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}; +const INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::MEDMEM_ORDER[N_MEDMEM_ORDER] = { INTERP_KERNEL::NORM_POINT1, INTERP_KERNEL::NORM_SEG2, INTERP_KERNEL::NORM_SEG3, INTERP_KERNEL::NORM_SEG4, INTERP_KERNEL::NORM_POLYL, INTERP_KERNEL::NORM_TRI3, INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_TRI7, INTERP_KERNEL::NORM_QUAD8, INTERP_KERNEL::NORM_QUAD9, INTERP_KERNEL::NORM_POLYGON, INTERP_KERNEL::NORM_QPOLYG, INTERP_KERNEL::NORM_TETRA4, INTERP_KERNEL::NORM_PYRA5, INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXA8, INTERP_KERNEL::NORM_HEXGP12, INTERP_KERNEL::NORM_TETRA10, INTERP_KERNEL::NORM_PYRA13, INTERP_KERNEL::NORM_PENTA15, INTERP_KERNEL::NORM_PENTA18, INTERP_KERNEL::NORM_HEXA20, INTERP_KERNEL::NORM_HEXA27, INTERP_KERNEL::NORM_POLYHED }; +const int MEDCouplingUMesh::MEDCOUPLING2VTKTYPETRADUCER[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,32,-1,25,42,36,4}; /// @endcond MEDCouplingUMesh *MEDCouplingUMesh::New() @@ -589,18 +590,15 @@ void MEDCouplingUMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const { checkFullyDefined(); - int nbOfNodes=getNumberOfNodes(); + int nbOfNodes(getNumberOfNodes()); int *revNodalIndxPtr=(int *)malloc((nbOfNodes+1)*sizeof(int)); revNodalIndx->useArray(revNodalIndxPtr,true,C_DEALLOC,nbOfNodes+1,1); std::fill(revNodalIndxPtr,revNodalIndxPtr+nbOfNodes+1,0); - const int *conn=_nodal_connec->getConstPointer(); - const int *connIndex=_nodal_connec_index->getConstPointer(); - int nbOfCells=getNumberOfCells(); - int nbOfEltsInRevNodal=0; + const int *conn(_nodal_connec->begin()),*connIndex(_nodal_connec_index->begin()); + int nbOfCells(getNumberOfCells()),nbOfEltsInRevNodal(0); for(int eltId=0;eltId=0)//for polyhedrons { @@ -838,10 +836,10 @@ void MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(const DataArrayInt *desc, cons { if(!desc || !descIndx || !revDesc || !revDescIndx) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ComputeNeighborsOfCellsAdv some input array is empty !"); - const int *descPtr=desc->getConstPointer(); - const int *descIPtr=descIndx->getConstPointer(); - const int *revDescPtr=revDesc->getConstPointer(); - const int *revDescIPtr=revDescIndx->getConstPointer(); + const int *descPtr=desc->begin(); + const int *descIPtr=descIndx->begin(); + const int *revDescPtr=revDesc->begin(); + const int *revDescIPtr=revDescIndx->begin(); // int nbCells=descIndx->getNumberOfTuples()-1; MCAuto out0=DataArrayInt::New(); @@ -863,6 +861,35 @@ void MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(const DataArrayInt *desc, cons neighborsIndx=out1.retn(); } +/*! + * Explodes \a this into edges whatever its dimension. + */ +MCAuto MEDCouplingUMesh::explodeIntoEdges(MCAuto& desc, MCAuto& descIndex, MCAuto& revDesc, MCAuto& revDescIndx) const +{ + checkFullyDefined(); + int mdim(getMeshDimension()); + desc=DataArrayInt::New(); descIndex=DataArrayInt::New(); revDesc=DataArrayInt::New(); revDescIndx=DataArrayInt::New(); + MCAuto mesh1D; + switch(mdim) + { + case 3: + { + mesh1D=explode3DMeshTo1D(desc,descIndex,revDesc,revDescIndx); + break; + } + case 2: + { + mesh1D=buildDescendingConnectivity(desc,descIndex,revDesc,revDescIndx); + break; + } + default: + { + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::computeNeighborsOfNodes : Mesh dimension supported are [3,2] !"); + } + } + return mesh1D; +} + /*! * \b WARNING this method do the assumption that connectivity lies on the coordinates set. * For speed reasons no check of this will be done. This method calls @@ -877,13 +904,15 @@ void MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(const DataArrayInt *desc, cons * The number of tuples is equal to the last values in \b neighborsIndx. * \param [out] neighborsIdx is an array of size this->getNumberOfCells()+1 newly allocated and should * be dealt by the caller. This arrays allow to use the first output parameter \b neighbors. + * + * \sa MEDCouplingUMesh::computeEnlargedNeighborsOfNodes */ void MEDCouplingUMesh::computeNeighborsOfNodes(DataArrayInt *&neighbors, DataArrayInt *&neighborsIdx) const { checkFullyDefined(); int mdim(getMeshDimension()),nbNodes(getNumberOfNodes()); MCAuto desc(DataArrayInt::New()),descIndx(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescIndx(DataArrayInt::New()); - MCAuto mesh1D; + MCConstAuto mesh1D; switch(mdim) { case 3: @@ -898,8 +927,7 @@ void MEDCouplingUMesh::computeNeighborsOfNodes(DataArrayInt *&neighbors, DataArr } case 1: { - mesh1D=const_cast(this); - mesh1D->incrRef(); + mesh1D.takeRef(this); break; } default: @@ -922,6 +950,45 @@ void MEDCouplingUMesh::computeNeighborsOfNodes(DataArrayInt *&neighbors, DataArr neighborsIdx=descIndx.retn(); } +/*! + * Computes enlarged neighbors for each nodes in \a this. The behavior of this method is close to MEDCouplingUMesh::computeNeighborsOfNodes except that the neighborhood of each node is wider here. + * A node j is considered to be in the neighborhood of i if and only if there is a cell in \a this containing in its nodal connectivity both i and j. + * This method is useful to find ghost cells of a part of a mesh with a code based on fields on nodes. + * + * \sa MEDCouplingUMesh::computeNeighborsOfNodes + */ +void MEDCouplingUMesh::computeEnlargedNeighborsOfNodes(MCAuto &neighbors, MCAuto& neighborsIdx) const +{ + checkFullyDefined(); + int nbOfNodes(getNumberOfNodes()); + const int *conn(_nodal_connec->begin()),*connIndex(_nodal_connec_index->begin()); + int nbOfCells(getNumberOfCells()); + std::vector< std::set > st0(nbOfNodes); + for(int eltId=0;eltId s(strtNdlConnOfCurCell,endNdlConnOfCurCell); s.erase(-1); //for polyhedrons + for(std::set::const_iterator iter2=s.begin();iter2!=s.end();iter2++) + st0[*iter2].insert(s.begin(),s.end()); + } + neighborsIdx=DataArrayInt::New(); neighborsIdx->alloc(nbOfNodes+1,1); neighborsIdx->setIJ(0,0,0); + { + int *neighIdx(neighborsIdx->getPointer()); + for(std::vector< std::set >::const_iterator it=st0.begin();it!=st0.end();it++,neighIdx++) + neighIdx[1]=neighIdx[0]+(*it).size()-1; + } + neighbors=DataArrayInt::New(); neighbors->alloc(neighborsIdx->back(),1); + { + const int *neighIdx(neighborsIdx->begin()); + int *neigh(neighbors->getPointer()),nodeId(0); + for(std::vector< std::set >::const_iterator it=st0.begin();it!=st0.end();it++,neighIdx++,nodeId++) + { + std::set s(*it); s.erase(nodeId); + std::copy(s.begin(),s.end(),neigh+*neighIdx); + } + } +} + /*! * Converts specified cells to either polygons (if \a this is a 2D mesh) or * polyhedrons (if \a this is a 3D mesh). The cells to convert are specified by an @@ -954,7 +1021,7 @@ void MEDCouplingUMesh::convertToPolyTypes(const int *cellIdsToConvertBg, const i int nbOfCells(getNumberOfCells()); if(dim==2) { - const int *connIndex=_nodal_connec_index->getConstPointer(); + const int *connIndex=_nodal_connec_index->begin(); int *conn=_nodal_connec->getPointer(); for(const int *iter=cellIdsToConvertBg;iter!=cellIdsToConvertEnd;iter++) { @@ -1917,19 +1984,19 @@ void MEDCouplingUMesh::setPartOfMySelf(const int *cellIdsBg, const int *cellIdsE oss << ", whereas other mesh dimension is set equal to " << otherOnSameCoordsThanThis.getMeshDimension() << " !"; throw INTERP_KERNEL::Exception(oss.str()); } - int nbOfCellsToModify=(int)std::distance(cellIdsBg,cellIdsEnd); + std::size_t nbOfCellsToModify(std::distance(cellIdsBg,cellIdsEnd)); if(nbOfCellsToModify!=otherOnSameCoordsThanThis.getNumberOfCells()) { std::ostringstream oss; oss << "MEDCouplingUMesh::setPartOfMySelf : cells ids length (" << nbOfCellsToModify << ") do not match the number of cells of other mesh (" << otherOnSameCoordsThanThis.getNumberOfCells() << ") !"; throw INTERP_KERNEL::Exception(oss.str()); } - int nbOfCells=getNumberOfCells(); - bool easyAssign=true; - const int *connI=_nodal_connec_index->getConstPointer(); - const int *connIOther=otherOnSameCoordsThanThis._nodal_connec_index->getConstPointer(); + std::size_t nbOfCells(getNumberOfCells()); + bool easyAssign(true); + const int *connI(_nodal_connec_index->begin()); + const int *connIOther=otherOnSameCoordsThanThis._nodal_connec_index->begin(); for(const int *it=cellIdsBg;it!=cellIdsEnd && easyAssign;it++,connIOther++) { - if(*it>=0 && *it=0 && *it<(int)nbOfCells) { easyAssign=(connIOther[1]-connIOther[0])==(connI[*it+1]-connI[*it]); } @@ -1967,7 +2034,7 @@ void MEDCouplingUMesh::setPartOfMySelfSlice(int start, int end, int step, const throw INTERP_KERNEL::Exception(oss.str()); } int nbOfCellsToModify=DataArray::GetNumberOfItemGivenBESRelative(start,end,step,"MEDCouplingUMesh::setPartOfMySelfSlice : "); - if(nbOfCellsToModify!=otherOnSameCoordsThanThis.getNumberOfCells()) + if(nbOfCellsToModify!=(int)otherOnSameCoordsThanThis.getNumberOfCells()) { std::ostringstream oss; oss << "MEDCouplingUMesh::setPartOfMySelfSlice : cells ids length (" << nbOfCellsToModify << ") do not match the number of cells of other mesh (" << otherOnSameCoordsThanThis.getNumberOfCells() << ") !"; throw INTERP_KERNEL::Exception(oss.str()); @@ -2296,7 +2363,7 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On DAInt neighIInit00(tmp11); // Neighbor information of the mesh WITH the crack (some neighbors are removed): DataArrayInt *idsTmp=0; - bool b=m01->areCellsIncludedIn(&otherDimM1OnSameCoords,2,idsTmp); + m01->areCellsIncludedIn(&otherDimM1OnSameCoords,2,idsTmp); DAInt ids(idsTmp); // In the neighbor information remove the connection between high dimension cells and its low level constituents which are part // of the frontier given in parameter (i.e. the cells of low dimension from the group delimiting the crack): @@ -2323,7 +2390,7 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On // Connex zone without the crack (to compute the next seed really) int dnu; DAInt connexCheck = MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(&seed, &seed+1, neighInit00,neighIInit00, -1, dnu); - int cnt = 0; + std::size_t cnt(0); for (int * ptr = connexCheck->getPointer(); cnt < connexCheck->getNumberOfTuples(); ptr++, cnt++) hitCells->setIJ(*ptr,0,1); // Connex zone WITH the crack (to identify cells lying on either part of the crack) @@ -2708,11 +2775,10 @@ DataArrayInt *MEDCouplingUMesh::getCellsInBoundingBox(const INTERP_KERNEL::Direc * \return INTERP_KERNEL::NormalizedCellType - enumeration item describing the cell type. * \throw If \a cellId is invalid. Valid range is [0, \a this->getNumberOfCells() ). */ -INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::getTypeOfCell(int cellId) const +INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::getTypeOfCell(std::size_t cellId) const { - const int *ptI=_nodal_connec_index->getConstPointer(); - const int *pt=_nodal_connec->getConstPointer(); - if(cellId>=0 && cellId<(int)_nodal_connec_index->getNbOfElems()-1) + const int *ptI(_nodal_connec_index->begin()),*pt(_nodal_connec->begin()); + if(cellId<_nodal_connec_index->getNbOfElems()-1) return (INTERP_KERNEL::NormalizedCellType) pt[ptI[cellId]]; else { @@ -2754,13 +2820,11 @@ DataArrayInt *MEDCouplingUMesh::giveCellsWithType(INTERP_KERNEL::NormalizedCellT /*! * Returns nb of cells having the geometric type \a type. No throw if no cells in \a this has the geometric type \a type. */ -int MEDCouplingUMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const +std::size_t MEDCouplingUMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const { - const int *ptI=_nodal_connec_index->getConstPointer(); - const int *pt=_nodal_connec->getConstPointer(); - int nbOfCells=getNumberOfCells(); - int ret=0; - for(int i=0;ibegin()),*pt(_nodal_connec->begin()); + std::size_t nbOfCells(getNumberOfCells()),ret(0); + for(std::size_t i=0;igetNumberOfCells() ). */ -void MEDCouplingUMesh::getNodeIdsOfCell(int cellId, std::vector& conn) const +void MEDCouplingUMesh::getNodeIdsOfCell(std::size_t cellId, std::vector& conn) const { - const int *ptI=_nodal_connec_index->getConstPointer(); - const int *pt=_nodal_connec->getConstPointer(); + const int *ptI(_nodal_connec_index->begin()),*pt(_nodal_connec->begin()); for(const int *w=pt+ptI[cellId]+1;w!=pt+ptI[cellId+1];w++) if(*w>=0) conn.push_back(*w); @@ -2977,14 +3040,14 @@ void MEDCouplingUMesh::setConnectivity(DataArrayInt *conn, DataArrayInt *connInd * Copy constructor. If 'deepCopy' is false \a this is a shallow copy of other. * If 'deeCpy' is true all arrays (coordinates and connectivities) are deeply copied. */ -MEDCouplingUMesh::MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCopy):MEDCouplingPointSet(other,deepCopy),_mesh_dim(other._mesh_dim), +MEDCouplingUMesh::MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCpy):MEDCouplingPointSet(other,deepCpy),_mesh_dim(other._mesh_dim), _nodal_connec(0),_nodal_connec_index(0), _types(other._types) { if(other._nodal_connec) - _nodal_connec=other._nodal_connec->performCopyOrIncrRef(deepCopy); + _nodal_connec=other._nodal_connec->performCopyOrIncrRef(deepCpy); if(other._nodal_connec_index) - _nodal_connec_index=other._nodal_connec_index->performCopyOrIncrRef(deepCopy); + _nodal_connec_index=other._nodal_connec_index->performCopyOrIncrRef(deepCpy); } MEDCouplingUMesh::~MEDCouplingUMesh() @@ -3010,7 +3073,7 @@ void MEDCouplingUMesh::computeTypes() * \return int - the number of cells in \a this mesh. * \throw If the nodal connectivity of cells is not defined. */ -int MEDCouplingUMesh::getNumberOfCells() const +std::size_t MEDCouplingUMesh::getNumberOfCells() const { if(_nodal_connec_index) return _nodal_connec_index->getNumberOfTuples()-1; @@ -3815,7 +3878,7 @@ void MEDCouplingUMesh::project1D(const double *pt, const double *v, double eps, MCAuto f=buildDirectionVectorField(); const double *fPtr=f->getArray()->getConstPointer(); double tmp[3]; - for(int i=0;igetNumberOfComponents()!=spaceDim) + if((int)pts->getNumberOfComponents()!=spaceDim) { std::ostringstream oss; oss << "MEDCouplingUMesh::distanceToPoints : input pts DataArrayDouble has " << pts->getNumberOfComponents() << " components whereas it should be equal to " << spaceDim << " (mesh spaceDimension) !"; throw INTERP_KERNEL::Exception(oss.str()); @@ -4788,6 +4851,27 @@ void MEDCouplingUMesh::orientCorrectlyPolyhedrons() updateTime(); } +/*! + * This method invert orientation of all cells in \a this. + * After calling this method the absolute value of measure of cells in \a this are the same than before calling. + * This method only operates on the connectivity so coordinates are not touched at all. + */ +void MEDCouplingUMesh::invertOrientationOfAllCells() +{ + checkConnectivityFullyDefined(); + std::set gts(getAllGeoTypes()); + int *conn(_nodal_connec->getPointer()); + const int *conni(_nodal_connec_index->begin()); + for(std::set::const_iterator gt=gts.begin();gt!=gts.end();gt++) + { + INTERP_KERNEL::AutoCppPtr oi(INTERP_KERNEL::OrientationInverter::BuildInstanceFrom(*gt)); + MCAuto cwt(giveCellsWithType(*gt)); + for(const int *it=cwt->begin();it!=cwt->end();it++) + oi->operate(conn+conni[*it]+1,conn+conni[*it+1]); + } + updateTime(); +} + /*! * Finds and fixes incorrectly oriented linear extruded volumes (INTERP_KERNEL::NORM_HEXA8, * INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXGP12 etc) to respect the MED convention @@ -5333,6 +5417,8 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTreeFast() const DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arcDetEps) const { checkFullyDefined(); + INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(arcDetEps); + int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells()); if(spaceDim!=2 || mDim!=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!"); @@ -5344,7 +5430,6 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arc { const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[*connI])); int sz(connI[1]-connI[0]-1); - INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=arcDetEps; std::vector nodes(sz); INTERP_KERNEL::QuadraticPolygon *pol(0); for(int j=0;j ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim); double *bbox(ret->getPointer()); const double *coords(_coords->begin()); @@ -5389,7 +5475,6 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic(double arc { const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[*connI])); int sz(connI[1]-connI[0]-1); - INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=arcDetEps; std::vector nodes(sz); INTERP_KERNEL::Edge *edge(0); for(int j=0;j=3) + if(nodalI[1]-nodalI[0]>=4) { + double aa[3]={coor[nodal[nodalI[0]+1+1]*3+0]-coor[nodal[nodalI[0]+1+0]*3+0], + coor[nodal[nodalI[0]+1+1]*3+1]-coor[nodal[nodalI[0]+1+0]*3+1], + coor[nodal[nodalI[0]+1+1]*3+2]-coor[nodal[nodalI[0]+1+0]*3+2]} + ,bb[3]={coor[nodal[nodalI[0]+1+2]*3+0]-coor[nodal[nodalI[0]+1+0]*3+0], + coor[nodal[nodalI[0]+1+2]*3+1]-coor[nodal[nodalI[0]+1+0]*3+1], + coor[nodal[nodalI[0]+1+2]*3+2]-coor[nodal[nodalI[0]+1+0]*3+2]}; + double cc[3]={aa[1]*bb[2]-aa[2]*bb[1],aa[2]*bb[0]-aa[0]*bb[2],aa[0]*bb[1]-aa[1]*bb[0]}; for(int j=0;j<3;j++) { int nodeId(nodal[nodalI[0]+1+j]); @@ -6404,14 +6496,34 @@ DataArrayDouble *MEDCouplingUMesh::computePlaneEquationOf3DFaces() const throw INTERP_KERNEL::Exception(oss.str()); } } + if(sqrt(cc[0]*cc[0]+cc[1]*cc[1]+cc[2]*cc[2])>1e-7) + { + INTERP_KERNEL::inverseMatrix(matrix,4,matrix2); + retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15]; + } + else + { + if(nodalI[1]-nodalI[0]==4) + { + std::ostringstream oss; oss << "MEDCouplingUMesh::computePlaneEquationOf3DFaces : cell" << i << " : Presence of The 3 colinear points !"; + throw INTERP_KERNEL::Exception(oss.str()); + } + // + double dd[3]={0.,0.,0.}; + for(int offset=nodalI[0]+1;offset()); + int nbOfNodesInCell(nodalI[1]-nodalI[0]-1); + std::transform(dd,dd+3,dd,std::bind2nd(std::multiplies(),1./(double)nbOfNodesInCell)); + std::copy(dd,dd+3,matrix+4*2); + INTERP_KERNEL::inverseMatrix(matrix,4,matrix2); + retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15]; + } } else { std::ostringstream oss; oss << "MEDCouplingUMesh::computePlaneEquationOf3DFaces : invalid 2D cell #" << i << " ! Must be constitued by more than 3 nodes !"; throw INTERP_KERNEL::Exception(oss.str()); } - INTERP_KERNEL::inverseMatrix(matrix,4,matrix2); - retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15]; } return ret.retn(); } @@ -7615,7 +7727,7 @@ bool MEDCouplingUMesh::RemoveIdsFromIndexedArrays(const int *idsToRemoveBg, cons previousArrI=*arrIPtr; *arrIPtr=(int)arrOut.size(); } - if(arr->getNumberOfTuples()==(int)arrOut.size()) + if(arr->getNumberOfTuples()==arrOut.size()) return false; arr->alloc((int)arrOut.size(),1); std::copy(arrOut.begin(),arrOut.end(),arr->getPointer());