From 07229eb54326cae50efc660a2a11be2691996cbe Mon Sep 17 00:00:00 2001 From: abn Date: Mon, 5 Feb 2024 15:34:39 +0100 Subject: [PATCH] [TetraIntersec] Fixing P1P0 barycentric computation + using proper barycenter computation (not a mere point average ...) for the center of mass of the intersection polygon. --- src/INTERP_KERNEL/TransformedTriangle.cxx | 26 +-- src/INTERP_KERNEL/TransformedTriangle.hxx | 3 +- .../UnitTetraIntersectionBary.cxx | 160 +++++++++++------- .../UnitTetraIntersectionBary.hxx | 9 +- 4 files changed, 123 insertions(+), 75 deletions(-) diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index b3478ae9c..909e9869b 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -105,7 +105,9 @@ namespace INTERP_KERNEL * @param r array of three doubles containing coordinates of R */ TransformedTriangle::TransformedTriangle(double* p, double* q, double* r) - : _is_double_products_calculated(false), _is_triple_products_calculated(false), _volume(0) + : _is_double_products_calculated(false), + _is_triple_products_calculated(false), + _volume(0) { for(int i = 0 ; i < 3 ; ++i) @@ -196,7 +198,6 @@ namespace INTERP_KERNEL LOG(2, "-- Calculating intersection polygons ... "); calculateIntersectionAndProjectionPolygons(); - double barycenter[3]; // calculate volume under A double volA = 0.0; @@ -208,12 +209,13 @@ namespace INTERP_KERNEL for(const auto& pt: _polygonA) LOG(3,vToStr(pt)); #endif - calculatePolygonBarycenter(A, barycenter); - sortIntersectionPolygon(A, barycenter); - volA = calculateVolumeUnderPolygon(A, barycenter); + calculatePolygonBarycenter(A, _barycenterA); // _barycenterA is re-used later in UnitTetraIntersector + sortIntersectionPolygon(A, _barycenterA); + volA = calculateVolumeUnderPolygon(A, _barycenterA); LOG(2, "Volume is " << sign * volA); } + double barycenterB[3]; double volB = 0.0; // if triangle is not in h = 0 plane, calculate volume under B if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ)) // second condition avoids double counting in case triangle fully included in h=0 facet @@ -224,22 +226,21 @@ namespace INTERP_KERNEL for(const auto& pt: _polygonB) LOG(3,vToStr(pt)); #endif - calculatePolygonBarycenter(B, barycenter); - sortIntersectionPolygon(B, barycenter); - volB = calculateVolumeUnderPolygon(B, barycenter); + calculatePolygonBarycenter(B, barycenterB); + sortIntersectionPolygon(B, barycenterB); + volB = calculateVolumeUnderPolygon(B, barycenterB); LOG(2, "Volume is " << sign * volB); } #if LOG_LEVEL >= 2 LOG(2, "############ Triangle :") dumpCoords(); - LOG(2, "vol A = " << volA); + LOG(2, "vol A = " << std::setprecision(15) << volA); LOG(2, "vol B = " << volB); LOG(2, "TOTAL = " << sign*(volA+volB)); #endif return _volume = sign * (volA + volB); - } /** @@ -264,9 +265,8 @@ namespace INTERP_KERNEL _volume = 0.; if(_polygonA.size() > 2) { - double barycenter[3]; - calculatePolygonBarycenter(A, barycenter); - sortIntersectionPolygon(A, barycenter); + calculatePolygonBarycenter(A, _barycenterA); + sortIntersectionPolygon(A, _barycenterA); const std::size_t nbPoints = _polygonA.size(); for(std::size_t i = 0 ; i < nbPoints ; ++i) tat->reverseApply(_polygonA[i], _polygonA[i]); diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 1c7aa1fd2..4e5980bb8 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -144,6 +144,7 @@ namespace INTERP_KERNEL const double* getCorner(TriCorner corner) const { return _coords + 5*corner; } const std::vector& getPolygonA() const { return _polygonA; } double getVolume() const { return _volume; } + const double* getBarycenterA() const { return _barycenterA; } protected: TransformedTriangle() { } @@ -272,7 +273,7 @@ namespace INTERP_KERNEL /// calculated volume for use of UnitTetraIntersectionBary double _volume; - + /** * Calls TransformedTriangle::testTriangleSurroundsEdge for edges OX to ZX and stores the result in * member variable array_triangleSurroundsEdgeCache. diff --git a/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx b/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx index 3494350a7..f4b007e4a 100644 --- a/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx +++ b/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx @@ -84,6 +84,9 @@ namespace INTERP_KERNEL crossprod<3>( triangle.getCorner(P),triangle.getCorner(Q),triangle.getCorner(R), triNormal); const std::vector * pPolygonA = &triangle.getPolygonA(); + const double * baryA = triangle.getBarycenterA(); + std::copy(baryA, baryA+3, _barycenterA); + if ( pPolygonA->size() < 3 ) { if ( !epsilonEqual( triNormal[_ZZ], 0 )) @@ -173,15 +176,19 @@ namespace INTERP_KERNEL clearPolygons(); // free memory of _polygonA _polygonA = faceCorner; _faces.pop_back(); + _facesBary.pop_back(); } else { if ( i < (int)pPolygonA->size() ) faceCorner.resize( i ); - + // Store normal and barycenter of the face: if ( _polyNormals.empty() ) _polyNormals.reserve(4); _polyNormals.push_back( std::vector< double >( polyNormal, polyNormal+3 )); + if (_facesBary.empty()) + _facesBary.reserve(4); + _facesBary.push_back( std::vector< double >( _barycenterA, _barycenterA+3 )); } #ifdef DMP_UNITTETRAINTERSECTIONBARY @@ -196,15 +203,69 @@ namespace INTERP_KERNEL clearPolygons(); } + /** Compute the face (rough) barycenter and normal vector for the last faces added by + * method addSideFaces() + * For the other first faces, information is already there. + */ + void UnitTetraIntersectionBary::completeFaceBarycentersAndNormals() + { + int nbFaceInit = (int)_facesBary.size(); + + int idx = 0, cnt = 0; + // Position ourselves at the right place in the list of faces: + auto face_it = _faces.begin(); + for(int i=0; i(3, 0.0)); + _polyNormals.emplace_back(std::vector(3, 0.0)); + const int npts = (int)face.size(); + if (face.empty() || npts < 3) continue; + + // Barycenter computation + auto & v = _facesBary.back(); + for (const auto& p: face) + for(int d=0; d<3; d++) + v[d] += p[d]; + for(int d=0; d<3; d++) + v[d] /= (double)npts; + + // Face normal computation + auto & n = _polyNormals.back(); + crossprod<3>(face[0], face[1], face[2], n.data()); + } + } + + /** Rough element barycenter computation. This is really just used to check whether a face of the + * (always convex) polyhedron is well oriented. + */ + void UnitTetraIntersectionBary::roughBarycenter(double * bary) + { + bary[0] = bary[1] = bary[2] = 0.0; + int idx = -1, cnt = 0; + for (const auto& face : _faces) + { + idx++; + if (face.empty()) continue; + for(int d=0; d<3; d++) + bary[d] += _facesBary[idx][d]; + cnt++; + } + for(int d=0; d<3; d++) + bary[d] /= (double)cnt; + } + //================================================================================ /*! * \brief Computes and returns coordinates of barycentre */ //================================================================================ - bool UnitTetraIntersectionBary::getBary(double* baryCenter) { baryCenter[0] = baryCenter[1] = baryCenter[2] = -1.0; + if ( addSideFaces() < NB_TETRA_SIDES ) { // tetra is not intersected @@ -217,75 +278,54 @@ namespace INTERP_KERNEL } return false; } - // Algo: - // - pick up one point P among the summits of the polyhedron - // - for each face of the polyhedron which does not contain the point : - // - compute the barycenter of the volume obtained by forming the "pyramid" with - // the face as a base and point P as a summit - // - compute the volume of the "pyramid" - // - Add up all barycenter positions weighting them with the volumes. - baryCenter[0] = baryCenter[1] = baryCenter[2] = 0.; + // After call to addSideFaces(), new faces were created - some face barycenters and normals are missing: + completeFaceBarycentersAndNormals(); - std::list< std::vector< double* > >::iterator f = _faces.begin(), fEnd = _faces.end(); - double * PP = f->at(0); + // Store coords and connectivity of intersection polyhedron to be able to call 'barycenterOfPolyhedron' + std::vector coords; + std::vector connec; + mcIdType cnt = 0; - for ( ++f; f != fEnd; ++f ) - { - std::vector< double* >& polygon = *f; - if ( polygon.empty() ) - continue; + // just an average of the bary of the faces,which are themselves not cleanly computed + double roughBary[3]; + roughBarycenter(roughBary); - bool pBelongsToPoly = false; - std::vector::iterator v = polygon.begin(), vEnd = polygon.end(); - for ( ; !pBelongsToPoly && v != vEnd; ++v ) - pBelongsToPoly = samePoint( PP, *v ); - if ( pBelongsToPoly ) - continue; + int i = 0; + for (auto it = _faces.begin(); it != _faces.end(); it++, i++) + { + const auto& face = *it; + if ((*it).size() < 3) continue; // Degenerate faces ... - // Compute the barycenter of the volume. Barycenter of pyramid is on line - // ( barycenter of polygon -> PP ) with 1/4 of pyramid height from polygon. + const auto& face_bary = _facesBary[i]; + const auto& norm = _polyNormals[i]; - double bary[] = { 0, 0, 0 }; + // Is the normal oriented in the same way as the vector joining the polyhedron (rough) barycenter and the face bary? + double barysV[3] = {0.,0.,0.}; + for (int d = 0; d<3; d++) barysV[d] = face_bary[d]-roughBary[d]; + const bool reverse = dotprod<3>(norm.data(), barysV) < 0.0; - // base polygon bary - for ( v = polygon.begin(); v != vEnd ; ++v ) + mcIdType cnt_prev = (mcIdType)connec.size(); + for(const auto& p: *it) { - double* p = *v; - bary[0] += p[0]; - bary[1] += p[1]; - bary[2] += p[2]; + for (int i=0; i < 3; i++) + coords.push_back(p[i]); + connec.push_back(cnt++); } - bary[0] /= (int)polygon.size(); - bary[1] /= (int)polygon.size(); - bary[2] /= (int)polygon.size(); + if(reverse) // TODO: optim: should store directly backward: + std::reverse(connec.data()+cnt_prev, connec.data()+connec.size()); + connec.push_back(-1); + } + connec.pop_back(); // Remove last -1: - // pyramid volume - double vol = 0; - for ( int i = 0; i < (int)polygon.size(); ++i ) - { - double* p1 = polygon[i]; - double* p2 = polygon[(i+1)%polygon.size()]; - vol += std::fabs( calculateVolumeForTetra( p1, p2, bary, PP )); - } + // + // Proper center of mass computation, even if several points aligned on an edge for example: + // + barycenterOfPolyhedron(connec.data(),connec.size(),coords.data(),baryCenter); - // put bary on the line ( barycenter of polygon -> PP ) and multiply by volume - baryCenter[0] += ( bary[0] * 0.75 + PP[0] * 0.25 ) * vol; - baryCenter[1] += ( bary[1] * 0.75 + PP[1] * 0.25 ) * vol; - baryCenter[2] += ( bary[2] * 0.75 + PP[2] * 0.25 ) * vol; - } if ( _int_volume < 0. ) _int_volume = -_int_volume; - baryCenter[0] /= _int_volume; - baryCenter[1] /= _int_volume; - baryCenter[2] /= _int_volume; -#ifdef DMP_UNITTETRAINTERSECTIONBARY - std::cout.precision(5); - std::cout << "**** Barycenter " << baryCenter[0] <<", "<< baryCenter[1] <<", "<< baryCenter[2] - << "\t **** Volume " << _int_volume << std::endl; - std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; -#endif return true; } @@ -310,7 +350,7 @@ namespace INTERP_KERNEL bool sideAdded[NB_TETRA_SIDES] = { false, false, false, false }; int nbAddedSides = 0; - std::list< std::vector< double* > >::iterator f = _faces.begin(), fEnd = _faces.end(); + auto f = _faces.begin(), fEnd = _faces.end(); for ( ; f != fEnd; ++f ) { std::vector< double* >& polygon = *f; @@ -727,11 +767,11 @@ namespace INTERP_KERNEL if ( andFaces ) { - std::list< std::vector< double* > >::iterator f = this->_faces.begin(), fEnd = this->_faces.end(); + auto f = this->_faces.begin(), fEnd = this->_faces.end(); for ( ; f != fEnd; ++f ) { std::vector< double* >& polygon = *f; - for(std::vector::iterator it = polygon.begin() ; it != polygon.end() ; ++it) + for(auto it = polygon.begin() ; it != polygon.end() ; ++it) { delete[] *it; *it = 0; diff --git a/src/INTERP_KERNEL/UnitTetraIntersectionBary.hxx b/src/INTERP_KERNEL/UnitTetraIntersectionBary.hxx index 596f46b2e..497cbacc2 100644 --- a/src/INTERP_KERNEL/UnitTetraIntersectionBary.hxx +++ b/src/INTERP_KERNEL/UnitTetraIntersectionBary.hxx @@ -66,12 +66,19 @@ namespace INTERP_KERNEL void clearPolygons(bool andFaces=false); + void completeFaceBarycentersAndNormals(); + + void roughBarycenter(double * bary); + /// volume of intersection double _int_volume; - /// faces of intersection polyhedron + /// Faces of intersection polyhedron std::list< std::vector< double* > > _faces; + /// Normals of faces: std::vector< std::vector< double > > _polyNormals; + /// Barycenter of faces: + std::vector< std::vector< double > > _facesBary; bool _isTetraInversed; }; -- 2.39.2