From: abn Date: Fri, 17 Feb 2017 09:21:39 +0000 (+0100) Subject: On the road to conformize3D ... faces implemented and new SkyLine array. X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=39e10933b1d75b238807eea0e935b9844d6a27ce;p=tools%2Fmedcoupling.git On the road to conformize3D ... faces implemented and new SkyLine array. --- diff --git a/src/INTERP_KERNEL/CellModel.hxx b/src/INTERP_KERNEL/CellModel.hxx index c581eadb2..cf839aa9e 100644 --- a/src/INTERP_KERNEL/CellModel.hxx +++ b/src/INTERP_KERNEL/CellModel.hxx @@ -32,7 +32,7 @@ namespace INTERP_KERNEL class DiameterCalculator; /*! - * This class descibes all static elements (different from polygons and polyhedron) 3D, 2D and 1D. + * This class describes all static elements (different from polygons and polyhedron) 3D, 2D and 1D. */ class CellModel { diff --git a/src/INTERP_KERNEL/PlanarIntersector.hxx b/src/INTERP_KERNEL/PlanarIntersector.hxx index 104eb187b..89fa1a156 100644 --- a/src/INTERP_KERNEL/PlanarIntersector.hxx +++ b/src/INTERP_KERNEL/PlanarIntersector.hxx @@ -60,7 +60,7 @@ namespace INTERP_KERNEL void getRealSourceCoordinatesPermute(ConnType icellS, int offset, std::vector& coordsS); void getRealCoordinates(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, std::vector& coordsT, std::vector& coordsS, int& orientation); double getValueRegardingOption(double val) const; - static void rotate3DTriangle( double* PP1, double*PP2, double*PP3, + static void Rotate3DTriangle( double* PP1, double*PP2, double*PP3, TranslationRotationMatrix& rotation_matrix); protected: const ConnType *_connectT; diff --git a/src/INTERP_KERNEL/PlanarIntersector.txx b/src/INTERP_KERNEL/PlanarIntersector.txx index a99e7a53d..92099bac1 100644 --- a/src/INTERP_KERNEL/PlanarIntersector.txx +++ b/src/INTERP_KERNEL/PlanarIntersector.txx @@ -385,7 +385,7 @@ namespace INTERP_KERNEL { TranslationRotationMatrix rotation; //rotate3DTriangle(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2], rotation); - rotate3DTriangle(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2, rotation); + Rotate3DTriangle(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2, rotation); for (int i=0; i - void PlanarIntersector::rotate3DTriangle(double* PP1, double*PP2, double*PP3, + void PlanarIntersector::Rotate3DTriangle(double* PP1, double*PP2, double*PP3, TranslationRotationMatrix& rotation_matrix) { //initializes diff --git a/src/INTERP_KERNEL/TranslationRotationMatrix.hxx b/src/INTERP_KERNEL/TranslationRotationMatrix.hxx index 3a227ef79..d5e79ecef 100644 --- a/src/INTERP_KERNEL/TranslationRotationMatrix.hxx +++ b/src/INTERP_KERNEL/TranslationRotationMatrix.hxx @@ -119,7 +119,43 @@ namespace INTERP_KERNEL rotate_vector(P); } - + static void Rotate3DTriangle(const double* PP1,const double*PP2,const double*PP3, TranslationRotationMatrix& rotation_matrix) + { + //initializes + rotation_matrix.translate(PP1); + + double P2w[3]; + double P3w[3]; + P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2]; + P3w[0]=PP3[0]; P3w[1]=PP3[1];P3w[2]=PP3[2]; + + // translating to set P1 at the origin + for (int i=0; i<3; i++) + { + P2w[i]-=PP1[i]; + P3w[i]-=PP1[i]; + } + + // rotating to set P2 on the Oxy plane + TranslationRotationMatrix A; + A.rotate_x(P2w); + A.rotate_vector(P3w); + rotation_matrix.multiply(A); + + //rotating to set P2 on the Ox axis + TranslationRotationMatrix B; + B.rotate_z(P2w); + B.rotate_vector(P3w); + rotation_matrix.multiply(B); + + //rotating to set P3 on the Oxy plane + TranslationRotationMatrix C; + C.rotate_x(P3w); + rotation_matrix.multiply(C); + } + + + private: static const double EPS; static const unsigned ROT_SIZE=9; diff --git a/src/INTERP_KERNEL/VectorUtils.hxx b/src/INTERP_KERNEL/VectorUtils.hxx index 375c36d4b..7ac0582bd 100644 --- a/src/INTERP_KERNEL/VectorUtils.hxx +++ b/src/INTERP_KERNEL/VectorUtils.hxx @@ -145,6 +145,19 @@ namespace INTERP_KERNEL // return std::fabs(x - y) < errTol; } + + /** + * Test whether two 3D vectors are colinear. The two vectors are expected to be of unit norm (not checked) + * Implemented by checking that the norm of the cross product is null. + */ + inline bool isColinear3D(const double *v1, const double *v2, const double eps = DEFAULT_ABS_TOL) + { + double cros[3]; + cross(v1, v2, cros); + return epsilonEqual(dot(cros, cros), 0.0, eps); + } + + /** * Compares doubles using a relative tolerance * This is suitable mainly for comparing larger values to each other. Before performing the relative test, diff --git a/src/MEDCoupling/MEDCouplingSkyLineArray.cxx b/src/MEDCoupling/MEDCouplingSkyLineArray.cxx index c0d1d52e1..dcd45b152 100644 --- a/src/MEDCoupling/MEDCouplingSkyLineArray.cxx +++ b/src/MEDCoupling/MEDCouplingSkyLineArray.cxx @@ -24,7 +24,7 @@ using namespace MEDCoupling; MEDCouplingSkyLineArray::MEDCouplingSkyLineArray(): - _index( DataArrayInt::New() ), _values( DataArrayInt::New() ), _sub_values( DataArrayInt::New() ) + _index( DataArrayInt::New() ), _values( DataArrayInt::New() ), _super_index( DataArrayInt::New() ) { } @@ -58,24 +58,124 @@ MEDCouplingSkyLineArray* MEDCouplingSkyLineArray::New( DataArrayInt* index, Data MEDCouplingSkyLineArray* MEDCouplingSkyLineArray::New( const MEDCouplingSkyLineArray & other ) { MEDCouplingSkyLineArray* ret = new MEDCouplingSkyLineArray(); + ret->_super_index = other._super_index; ret->_index = other._index; ret->_values = other._values; - ret->_sub_values = other._sub_values; return ret; } +/**! Build a three level SkyLine array from the dynamic connectivity of a dynamic mesh (i.e. containing only + * polyhedrons or polygons). + * The input arrays are deep copied, contrary to the other ctors. + */ +MEDCouplingSkyLineArray * MEDCouplingSkyLineArray::BuildFromPolyhedronConn( const DataArrayInt* c, const DataArrayInt* cI ) +{ + using namespace std; + + MEDCouplingSkyLineArray* ret = new MEDCouplingSkyLineArray(); + + const int * cP(c->begin()), * cIP(cI->begin()); + int prev = -1; + if (c->getNbOfElems() != *(cI->end()-1)) + throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::BuildFromDynamicConn: misformatted connectivity (wrong nb of tuples)!"); + for (int i=0; i < cI->getNbOfElems(); i++) + { + int j = cIP[i]; + if (cIP[i] < prev) + throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::BuildFromDynamicConn: misformatted connectivity (indices not monotonic ascending)!"); + prev = cIP[i]; + if (i!=cI->getNbOfElems()-1) + if (cP[j] != INTERP_KERNEL::NORM_POLYHED) + throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::BuildFromDynamicConn: connectivity containing other types than POLYHED!"); + } + + vector superIdx, idx, vals; + int cnt = 0, cnt2 = 0; + superIdx.reserve(cI->getNbOfElems()); + superIdx.push_back(0); + idx.push_back(0); + vals.resize(c->getNbOfElems()); // too much because of the type and the -1, but still better than push_back(). + for (int i=0; i < cI->getNbOfElems()-1; i++) + { + int start = cIP[i]+1, end = cIP[i+1]; + int * work = vals.data() + cnt; + const int * w = cP+start; + const int * w2 = find(w, cP+end, -1); + while (w2 != cP+end) + { + copy(w, w2, work); + int d = distance(w, w2); + cnt += d; work +=d; + idx.push_back(cnt); cnt2++; + w = w2+1; // skip the -1 + w2 = find(w, cP+end, -1); + } + copy(w, cP+end, work); + cnt += distance(w, cP+end); + idx.push_back(cnt); cnt2++; + superIdx.push_back(cnt2); + } + ret->_super_index->alloc(superIdx.size(),1); + copy(superIdx.begin(), superIdx.end(), ret->_super_index->getPointer()); + ret->_index->alloc(idx.size(),1); + copy(idx.begin(), idx.end(), ret->_index->getPointer()); + ret->_values->alloc(cnt,1); + copy(vals.begin(), vals.begin()+cnt, ret->_values->getPointer()); + + return ret; +} + +/** + * Convert a three-level SkyLineArray into a polyhedral connectivity. + * The super-packs are interpreted as cell description, and the packs represent the face connectivity. + */ +void MEDCouplingSkyLineArray::convertToPolyhedronConn( MCAuto& c, MCAuto& cI) const +{ + // TODO: in this case an iterator would be nice + using namespace std; + + checkSuperIndex("convertToPolyhedronConn"); + + const int * siP(_super_index->begin()), * iP(_index->begin()), *vP(_values->begin()); + int cnt = 0; + cI->alloc(_super_index->getNbOfElems(),1); // same number of super packs as number of cells + std::vector cVec, cIVec; + MCAuto dsi = _index->deltaShiftIndex(); + int sz = dsi->accumulate(0) + dsi->getNbOfElems(); // think about it: one slot for the type, -1 at the end of each face of the cell + cVec.reserve(sz); + int * cVecP(cVec.data()); + c->alloc(sz, 1); + + for (int i=0; i < _super_index->getNbOfElems()-1; i++) + { + cIVec.push_back(cnt); + int endId = siP[i+1]; + cVecP[cnt++] = INTERP_KERNEL::NORM_POLYHED; + for (int j=siP[i]; j < endId; j++) + { + int startId2 = iP[j], endId2 = iP[j+1]; + copy(vP+startId2, vP+endId2, cVecP+cnt); + cnt += endId2-startId2; + if(j != endId-1) + cVecP[cnt++] = -1; + } + } + cIVec.push_back(cnt); + copy(cIVec.begin(), cIVec.end(), cI->getPointer()); + copy(cVec.begin(), cVec.begin()+cnt, c->getPointer()); +} std::size_t MEDCouplingSkyLineArray::getHeapMemorySizeWithoutChildren() const { - return _index->getHeapMemorySizeWithoutChildren()+_values->getHeapMemorySizeWithoutChildren()+_sub_values->getHeapMemorySizeWithoutChildren(); + return _index->getHeapMemorySizeWithoutChildren()+_values->getHeapMemorySizeWithoutChildren()+_super_index->getHeapMemorySizeWithoutChildren(); } std::vector MEDCouplingSkyLineArray::getDirectChildrenWithNull() const { std::vector ret; + ret.push_back(_super_index); ret.push_back(_index); ret.push_back(_values); - ret.push_back(_sub_values); return ret; } @@ -90,6 +190,20 @@ void MEDCouplingSkyLineArray::set( DataArrayInt* index, DataArrayInt* value ) else _values = DataArrayInt::New(); } +void MEDCouplingSkyLineArray::set3( DataArrayInt* superIndex, DataArrayInt* index, DataArrayInt* value ) +{ + _super_index=superIndex; + if ( (DataArrayInt*)_super_index ) _super_index->incrRef(); + else _super_index = DataArrayInt::New(); + set(index, value); +} + +DataArrayInt* MEDCouplingSkyLineArray::getSuperIndexArray() const +{ + return ((MEDCouplingSkyLineArray*)this)->_super_index; +} + + DataArrayInt* MEDCouplingSkyLineArray::getIndexArray() const { return ((MEDCouplingSkyLineArray*)this)->_index; @@ -100,19 +214,49 @@ DataArrayInt* MEDCouplingSkyLineArray::getValuesArray() const return ((MEDCouplingSkyLineArray*)this)->_values; } +void MEDCouplingSkyLineArray::checkSuperIndex(const std::string& func) const +{ + if (!_super_index->getNbOfElems()) + { + std::ostringstream oss; + oss << "MEDCouplingSkyLineArray::"<< func << ": not a three level SkyLineArray! Method is not available for two-level SkyLineArray."; + throw INTERP_KERNEL::Exception(oss.str()); + } +} + std::string MEDCouplingSkyLineArray::simpleRepr() const { std::ostringstream oss; - oss << "MEDCouplingSkyLineArray" << std::endl; - oss << " Nb of items: " << getNumberOf() << std::endl; + oss << "MEDCouplingSkyLineArray (" << this << ")" << std::endl; + MCAuto super_index = _super_index->deepCopy(); + if (_super_index->getNbOfElems()) + oss << " Nb of super-packs: " << getSuperNumberOf() << std::endl; + else + { + super_index->alloc(2,1); + super_index->setIJSilent(0,0,0); + super_index->setIJSilent(1,0,_index->getNbOfElems()-1); + } + oss << " Nb of packs: " << getNumberOf() << std::endl; oss << " Nb of values: " << getLength() << std::endl; - oss << " Index:" << std::endl; + + if (_super_index->getNbOfElems()) + { + oss << " Super-indices:" << std::endl; + oss << " "; + const int * i = _super_index->begin(); + for ( ; i != _super_index->end(); ++i ) + oss << *i << " "; + oss << std::endl; + } + + oss << " Indices:" << std::endl; oss << " "; const int * i = _index->begin(); for ( ; i != _index->end(); ++i ) oss << *i << " "; oss << std::endl; - oss << " Value:" << std::endl; + oss << " Values:" << std::endl; oss << " "; const int * v = _values->begin(); int cnt = 0; @@ -129,3 +273,109 @@ std::string MEDCouplingSkyLineArray::simpleRepr() const return oss.str(); } + +void MEDCouplingSkyLineArray::getPackSafe(const int absolutePackId, std::vector & pack) const +{ + if(absolutePackId < 0 || absolutePackId >= _index->getNbOfElems()) + throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::getPackSafe: invalid index!"); + const int * iP(_index->begin()), *vP(_values->begin()); + int sz = iP[absolutePackId+1]-iP[absolutePackId]; + pack.resize(sz); + std::copy(vP+iP[absolutePackId], vP+iP[absolutePackId+1],pack.begin()); +} + + +/**! + * For each given super-pack ID, provide the sub-index of the first matching pack. If no matching pack is found for the + * given super-pack -1 is returned. + * \param[in] superPackIndices the list of super-packs that should be inspected + * \param[in] pack the pack that the function is looking for in each of the provided super-pack + * \param[out] a vector of int, having the same size as superPackIndices and containing for each inspected super-pack + * the index of the first matching pack, or -1 if none found. + */ +void MEDCouplingSkyLineArray::findPackIds(const std::vector & superPackIndices, const std::vector & pack, + std::vector& out) const +{ + using namespace std; + + checkSuperIndex("findPackIds"); + out.resize(superPackIndices.size()); + int i = 0; + const int * siP(_super_index->begin()), * iP(_index->begin()), *vP(_values->begin()); + for(vector::const_iterator it=superPackIndices.begin(); it!=superPackIndices.end(); ++it, i++) + { + out[i] = -1; + const int sPackIdx = *it; + // for each pack + for (int idx=siP[sPackIdx], j=0; idx < siP[sPackIdx+1]; idx++, j++) + if (equal(&vP[iP[idx]], &vP[iP[idx+1]], pack.begin())) + { + out[i] = j; + break; + } + } +} + +/**! + * Delete pack number j in super-pack number i. + * \param[in] i is the super-pack number + * \param[in] j is the pack index inside the super-pack i. + */ +void MEDCouplingSkyLineArray::deletePack(const int i, const int j) +{ + using namespace std; + + checkSuperIndex("deletePack"); + + int * vP = _values->getPointer(); + int * siP(_super_index->getPointer()), *iP(_index->getPointer()); + const int start = iP[siP[i]+j], end = iP[siP[i]+j+1]; + // _values + copy(vP+end, vP+_values->getNbOfElems(), vP+start); + _values->reAlloc(_values->getNbOfElems() - (end-start)); + + // _index + int nt = _index->getNbOfElems(); + copy(iP+siP[i]+j+1, iP+nt, iP+siP[i]+j); + _index->reAlloc(nt-1); iP = _index->getPointer(); // better not forget this ... + for(int ii = siP[i]; ii < nt-1; ii++) + iP[ii] -= (end-start); + + // _super_index + for(int ii = i+1; ii < _super_index->getNbOfElems(); ii++) + (siP[ii])--; +} + +/**! + * Insert a new pack in super-pack at index i. The pack is inserted at the end of the pack list of the chosen super-pack. + */ +void MEDCouplingSkyLineArray::pushBackPack(const int i, const std::vector & pack) +{ + using namespace std; + + checkSuperIndex("pushBackPack"); + + int *siP(_super_index->getPointer()), *iP(_index->getPointer()); + const int sz(pack.size()); + + // _values + _values->reAlloc(_values->getNbOfElems()+sz); + int * vPE(_values->getPointer()+_values->getNbOfElems()); + int *vP(_values->getPointer()); + copy(vP+iP[siP[i+1]], vPE-sz, vP+iP[siP[i+1]]+sz); + // insert pack + copy(pack.begin(), pack.end(), vP+iP[siP[i+1]]); + + // _index + int nt = _index->getNbOfElems(); + _index->reAlloc(nt+1); iP = _index->getPointer(); + copy(iP+siP[i+1]+1, iP+nt, iP+siP[i+1]+2); + iP[siP[i+1]+1] = iP[siP[i+1]] + pack.size(); + for(int ii = siP[i+1]+2; ii < nt+1; ii++) + iP[ii] += sz; + + // _super_index + for(int ii = i+1; ii < _super_index->getNbOfElems(); ii++) + (siP[ii])++; +} + diff --git a/src/MEDCoupling/MEDCouplingSkyLineArray.hxx b/src/MEDCoupling/MEDCouplingSkyLineArray.hxx index dd943a126..7ee80ecd0 100644 --- a/src/MEDCoupling/MEDCouplingSkyLineArray.hxx +++ b/src/MEDCoupling/MEDCouplingSkyLineArray.hxx @@ -30,13 +30,29 @@ namespace MEDCoupling { /**! - * Class allowing the easy manipulation of the indexed array format, where the first array is a set of offsets to - * be used in the second array, to extract packs of values. + * Class allowing the easy manipulation of the indexed array format, where the first array (the "indices") is a set of offsets to + * be used in the second array (the "values"), to extract packs of values. * - * This class allows to pursuie this logic up to 3 levels, i.e. the first array points to packs in the second, which - * itself points to identifiers in the third array. + * Example: + * index = [0,2,7,10] + * values = [1,24,2,33,6,7,10,11,9,28] + * which describes 3 packs of (integer) values : [1,24] and [2,33,6,7,10] and [11,9,28] + * + * Thus the index array is always monotic ascendant. + * + * This class allows to pursue this logic up to 3 levels, i.e. the first array (the "indices") points to packs in the + * second (the "values"), which itself points to identifiers in a third array (the "sub-values"). * * This particularly useful for connectivity of pure polygonal/polyhedral meshes. + * + * Example: + * super-index = [0,1,3] + * index = [0,3,6,10] + * values = [28,1,4,2,35,8,9,10,1,12] + * which represent two 3 packs and two super-packs. The first super-pack is [[28,1,4]] and has only one pack [28,1,4]. + * The second super-pack is [[2,35,8], [9,10,1,12]] and has two packs [2,35,8] and [9,10,1,12]. + * Note that contrary to index, the integers in super-index are interpreted as being inclusive: the first super-pack + * goes from offset 0 to offset 1, inclusive. This is not the same for index, where the upper bound is exclusive. */ class MEDCOUPLING_EXPORT MEDCouplingSkyLineArray : public RefCountObject { @@ -46,30 +62,45 @@ namespace MEDCoupling static MEDCouplingSkyLineArray * New( DataArrayInt* index, DataArrayInt* value ); static MEDCouplingSkyLineArray * New( const MEDCouplingSkyLineArray & other ); + static MEDCouplingSkyLineArray * BuildFromPolyhedronConn( const DataArrayInt* c, const DataArrayInt* cI ); + std::size_t getHeapMemorySizeWithoutChildren() const; std::vector getDirectChildrenWithNull() const; void set( DataArrayInt* index, DataArrayInt* value ); + void set3( DataArrayInt* superIndex, DataArrayInt* index, DataArrayInt* value ); + int getSuperNumberOf() const { return _super_index->getNbOfElems()-1; } int getNumberOf() const { return _index->getNbOfElems()-1; } int getLength() const { return _values->getNbOfElems(); } + + const int* getSuperIndex() const { return _super_index->begin(); } const int* getIndex() const { return _index->begin(); } const int* getValues() const { return _values->begin(); } + DataArrayInt* getSuperIndexArray() const; DataArrayInt* getIndexArray() const; DataArrayInt* getValuesArray() const; std::string simpleRepr() const; -// replaceWithPackFromOther() + void getPackSafe(const int absolutePackId, std::vector & pack) const; + void findPackIds(const std::vector & superPackIndices, const std::vector & pack, std::vector& out) const; + + void deletePack(const int i, const int j); + void pushBackPack(const int i, const std::vector & pack); + + void convertToPolyhedronConn( MCAuto& c, MCAuto& cI) const; private: MEDCouplingSkyLineArray(); ~MEDCouplingSkyLineArray(); + void checkSuperIndex(const std::string& func) const; + + MCAuto _super_index; MCAuto _index; MCAuto _values; - MCAuto _sub_values; }; } # endif diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 1a3e947e0..6192f0565 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -231,6 +231,7 @@ namespace MEDCoupling MEDCOUPLING_EXPORT DataArrayDouble *computePlaneEquationOf3DFaces() const; MEDCOUPLING_EXPORT DataArrayInt *conformize2D(double eps); MEDCOUPLING_EXPORT DataArrayInt *colinearize2D(double eps); + MEDCOUPLING_EXPORT DataArrayInt *conformize3D(double eps); MEDCOUPLING_EXPORT int split2DCells(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *midOpt=0, const DataArrayInt *midOptI=0); MEDCOUPLING_EXPORT static MEDCouplingUMesh *Build0DMeshFromCoords(DataArrayDouble *da); MEDCOUPLING_EXPORT static MCAuto Build1DMeshFromCoords(DataArrayDouble *da); diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index 8fdf01882..49bbc7bb9 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx @@ -20,6 +20,7 @@ #include "MEDCouplingUMesh.hxx" #include "MEDCoupling1GTUMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" #include "CellModel.hxx" #include "VolSurfUser.txx" #include "InterpolationUtils.hxx" @@ -33,6 +34,9 @@ #include "InterpKernelGeo2DEdgeLin.hxx" #include "InterpKernelGeo2DEdgeArcCircle.hxx" #include "InterpKernelGeo2DQuadraticPolygon.hxx" +#include "TranslationRotationMatrix.hxx" +#include "VectorUtils.hxx" +#include "MEDCouplingSkyLineArray.hxx" #include #include @@ -1994,4 +1998,164 @@ DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps) return ret.retn(); } +/*! + * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty). + * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) ! + * This method performs a conformization of \b this. + * + * Only polyhedron cells are supported. You can call convertAllToPoly() + * + * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too. + * This method expects that all nodes in \a this are not closer than \a eps. + * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method. + * + * \param [in] eps the relative error to detect merged edges. + * \return DataArrayInt * - The list of cellIds in \a this that have been subdivided. If empty, nothing changed in \a this (as if it were a const method). The array is a newly allocated array + * that the user is expected to deal with. + * + * \throw If \a this is not coherent. + * \throw If \a this has not spaceDim equal to 3. + * \throw If \a this has not meshDim equal to 3. + * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly() + */ +DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) +{ + using namespace std; + + static const int SPACEDIM=3; + checkConsistencyLight(); + if(getSpaceDimension()!=3 || getMeshDimension()!=3) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!"); + if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first."); + MCAuto desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New()); + MCAuto mDesc(buildDescendingConnectivity(desc,descI,revDesc,revDescI)); + const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer()); + int *c(getNodalConnectivity()->getPointer()),*cI(getNodalConnectivityIndex()->getPointer()); + MCAuto connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex())); + const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin()); + MCAuto connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivity(), + mDesc->getNodalConnectivityIndex())); + const double * cooDesc(mDesc->_coords->begin()); + MCAuto bboxArr(mDesc->getBoundingBoxForBBTree()); + const double *bbox(bboxArr->begin()),*coords(getCoords()->begin()); + int nCell=getNumberOfCells(), nDescCell=mDesc->getNumberOfCells(); + + std::vector addCoo; + BBTree myTree(bbox,0,0,nDescCell,-eps); + // Surfaces - handle biggest first + MCAuto surfF = mDesc->getMeasureField(true); + DataArrayDouble * surfs = surfF->getArray(); + // Normal field + MCAuto normalsF = mDesc->buildOrthogonalField(); + DataArrayDouble * normals = normalsF->getArray(); + const double * normalsP = normals->getConstPointer(); + + // Sort faces by decreasing surface: + vector> S; + for(int i=0;i < surfs->getNumberOfTuples();i++){ + pair p = make_pair(surfs->getIJ(i, 0), i); + S.push_back(p); + } + sort(S.rbegin(),S.rend()); // reverse sort + vector hit(nDescCell); + fill(hit.begin(), hit.end(), false); + vector hitPoly; // the final result: which 3D cells have been modified. + + // STEP 1 -- faces + for( vector>::const_iterator it = S.begin(); it != S.end(); it++) + { + int sIdx = (*it).second; + if (hit[sIdx]) continue; + + vector candidates, cands2; + myTree.getIntersectingElems(bbox+sIdx*2*SPACEDIM,candidates); + // Keep only candidates whose normal matches the normal of current face + for(vector::const_iterator it=candidates.begin();it!=candidates.end();it++) + { + bool col = INTERP_KERNEL::isColinear3D(normalsP + (*it)*sIdx, normalsP + *(it)*SPACEDIM, eps); + if (*it != sIdx && col) + cands2.push_back(*it); + } + if (!cands2.size()) + continue; + + // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later + INTERP_KERNEL::TranslationRotationMatrix rotation; + INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(cooDesc+SPACEDIM*(cDesc[cIDesc[sIdx]+1]), + cooDesc+SPACEDIM*(cDesc[cIDesc[sIdx]+2]), + cooDesc+SPACEDIM*(cDesc[cIDesc[sIdx]+3]), rotation); + + MCAuto mPartRef(buildPartOfMySelf(&sIdx, &sIdx, false)); + MCAuto mPartCand(buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); + double * cooPartRef(mPartCand->_coords->getPointer()); + double * cooPartCand(mPartCand->_coords->getPointer()); + for (int ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) + rotation.transform_vector(cooPartRef+SPACEDIM*ii); + for (int ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++) + rotation.transform_vector(cooPartCand+SPACEDIM*ii); + double ze_z = cooPartRef[2]; + + MCAuto baryPart = mPartCand->computeCellCenterOfMass(); + vector compo; compo.push_back(2); + MCAuto baryPartZ = baryPart->keepSelectedComponents(compo); + MCAuto idsGoodPlane = baryPartZ->findIdsInRange(ze_z-eps, ze_z+eps); + if (!idsGoodPlane->getNumberOfTuples()) + continue; + + MCAuto cc(DataArrayInt::New()), ccI(DataArrayInt::New()); + mPartRef->getCellsContainingPoints(baryPart->begin(), baryPart->getNumberOfTuples(), eps, cc, ccI); + + // Localize faces in 2D thanks to barycenters + if (!cc->getNumberOfTuples()) + continue; + MCAuto dsi = ccI->deltaShiftIndex(); + + { + MCAuto tmp = dsi->findIdsInRange(0, 2); + if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples()) + throw INTERP_KERNEL::Exception("Non expected non-conformity. Only simple (=partition-like) non-conformities are handled."); + } + // TODO check matching total surface! + + MCAuto ids = dsi->findIdsEqual(1); + // Boundary face: + if (!ids->getNumberOfTuples()) + continue; + + for (const int * ii = ids->begin(); ii != ids->end(); ii++) + hit[cands2[idsGoodPlane->getIJ(*ii,0)]] = true; + + // Current face connectivity + const int * sIdxConn = cDesc + cIDesc[sIdx] + 1; + + // For all polyhedrons using this face, replace connectivity: + vector polyIndices, subPacksIds, facePack; + for (int ii=*(revDescIP+sIdx); ii < *(revDescIP+sIdx+1); ii++) + polyIndices.push_back(*(revDescP+ii)); + connSla->findPackIds(polyIndices, facePack, subPacksIds); + // Deletion of old faces + int jj=0; + for (vector::const_iterator it=polyIndices.begin(); it!=polyIndices.end(); ++it, ++jj) + connSla->deletePack(*it, subPacksIds[jj]); + // Insertion of new faces: + for (const int * ii = ids->begin(); ii != ids->end(); ii++) + { + // Build pack from the face to insert: + int faceIdx = cands2[idsGoodPlane->getIJ(*ii,0)]; + vector facePack; + connSlaDesc->getPackSafe(faceIdx, facePack); + // Insert it in all hit polyhedrons: + for (vector::const_iterator it=polyIndices.begin(); it!=polyIndices.end(); ++it) + connSla->pushBackPack(*it, facePack); + } + } + + // STEP 2 -- edges + // TODO + + MCAuto ret(DataArrayInt::New()); + return ret.retn(); +} + diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py index bb4c1e205..28d9f2d2c 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py @@ -4166,6 +4166,51 @@ class MEDCouplingBasicsTest5(unittest.TestCase): pass + def testMEDCouplingSkyLineArrayThreeLevels(self): + # [[28,1,4]] , [[2,35,8], [9,10,1,12]] + superi = DataArrayInt([ 0,1,3 ]) + index = DataArrayInt ([ 0,3,6,10 ]) + value = DataArrayInt ([ 28,1,4,2,35,8,9,10,1,12 ]) + + sla0 = MEDCouplingSkyLineArray() + self.assertEqual( -1, sla0.getSuperNumberOf() ) + self.assertEqual( -1, sla0.getNumberOf() ) + self.assertEqual( 0, sla0.getLength() ) + sla0.set3( superi.deepCopy(), index.deepCopy(), value.deepCopy() ) + self.assertTrue( superi.isEqual( sla0.getSuperIndexArray() )) + + pack = sla0.getPackSafe(2) + self.assertEqual([9,10,1,12], pack) + ids = sla0.findPackIds([0,1], [9,10,1,12]) + self.assertEqual([-1,1], ids) + + sla0.deletePack(1, 0) + si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray() + self.assertEqual([28,1,4,9,10,1,12], val.getValues()) + self.assertEqual([0,3,7], idx.getValues()) + self.assertEqual([0,1,2], si.getValues()) + + sla0.pushBackPack(0, [3,2,1,0]) + si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray() + self.assertEqual([0,2,3], si.getValues()) + self.assertEqual([0,3,7,11], idx.getValues()) + self.assertEqual([28,1,4,3,2,1,0, 9,10,1,12], val.getValues()) + + # Build connectivity from POLYHED connectivity + cI = [0,16,41] + c = [NORM_POLYHED, 1,2,3,-1, 2,3,4,-1, 3,4,5,-1, 4,5,6, + NORM_POLYHED, 7,8,9,10,-1, 9,10,11,12,-1, 3,4,5,6,-1, 5,6,7,8,-1, 9,10,11,12] + sla0 = MEDCouplingSkyLineArray.BuildFromPolyhedronConn(DataArrayInt(c), DataArrayInt(cI)) + si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray() + self.assertEqual([0,4,9], si.getValues()) + self.assertEqual([0,3,6,9,12,16,20,24,28,32], idx.getValues()) + self.assertEqual([1,2,3, 2,3,4, 3,4,5, 4,5,6, + 7,8,9,10, 9,10,11,12, 3,4,5,6, 5,6,7,8, 9,10,11,12], val.getValues()) + c1, cI1 = sla0.convertToPolyhedronConn() + self.assertEqual(c1.getValues(), c) + self.assertEqual(cI1.getValues(), cI) + pass + def testMEDCouplingUMeshgenerateGraph(self): # cartesian mesh 3x3 arr=DataArrayDouble(4) ; arr.iota() diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index bad2f7127..88a57ba8c 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -410,6 +410,7 @@ using namespace INTERP_KERNEL; %newobject MEDCoupling::DenseMatrix::__mul__; %newobject MEDCoupling::MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell; %newobject MEDCoupling::MEDCouplingGaussLocalization::buildRefCell; +%newobject MEDCoupling::MEDCouplingSkyLineArray::BuildFromPolyhedronConn; %feature("unref") MEDCouplingPointSet "$this->decrRef();" %feature("unref") MEDCouplingMesh "$this->decrRef();" @@ -1178,11 +1179,21 @@ namespace MEDCoupling class MEDCouplingSkyLineArray { public: + static MEDCouplingSkyLineArray *BuildFromPolyhedronConn( const DataArrayInt* c, const DataArrayInt* cI ) throw(INTERP_KERNEL::Exception); + void set( DataArrayInt* index, DataArrayInt* value ); + void set3( DataArrayInt* superIndex, DataArrayInt* index, DataArrayInt* value ); + + int getSuperNumberOf() const; int getNumberOf() const; int getLength() const; + + DataArrayInt* getSuperIndexArray() const; DataArrayInt* getIndexArray() const; DataArrayInt* getValuesArray() const; + + void deletePack(const int i, const int j) throw(INTERP_KERNEL::Exception); + %extend { MEDCouplingSkyLineArray() throw(INTERP_KERNEL::Exception) @@ -1209,7 +1220,40 @@ namespace MEDCoupling { return self->simpleRepr(); } - + + PyObject *getPackSafe(int absolutePackId) const throw(INTERP_KERNEL::Exception) + { + std::vector ret; + self->getPackSafe(absolutePackId,ret); + return convertIntArrToPyList2(ret); + } + + PyObject *findPackIds(PyObject *superPackIndices, PyObject *pack) const throw(INTERP_KERNEL::Exception) + { + std::vector vpack, vspIdx, out; + convertPyToNewIntArr3(superPackIndices,vspIdx); + convertPyToNewIntArr3(pack,vpack); + self->findPackIds(vspIdx,vpack, out); + return convertIntArrToPyList2(out); + } + + void pushBackPack(const int i, PyObject *pack) throw(INTERP_KERNEL::Exception) + { + std::vector vpack; + convertPyToNewIntArr3(pack,vpack); + self->pushBackPack(i,vpack); + } + + PyObject *convertToPolyhedronConn() const throw(INTERP_KERNEL::Exception) + { + MCAuto d0=DataArrayInt::New(); + MCAuto d1=DataArrayInt::New(); + self->convertToPolyhedronConn(d0,d1); + PyObject *ret=PyTuple_New(2); + PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(d0.retn()),SWIGTYPE_p_MEDCoupling__DataArrayInt, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(d1.retn()),SWIGTYPE_p_MEDCoupling__DataArrayInt, SWIG_POINTER_OWN | 0 )); + return ret; + } } }; }