* @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)
LOG(2, "-- Calculating intersection polygons ... ");
calculateIntersectionAndProjectionPolygons();
- double barycenter[3];
// calculate volume under A
double volA = 0.0;
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
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);
-
}
/**
_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]);
crossprod<3>( triangle.getCorner(P),triangle.getCorner(Q),triangle.getCorner(R), triNormal);
const std::vector<double*> * pPolygonA = &triangle.getPolygonA();
+ const double * baryA = triangle.getBarycenterA();
+ std::copy(baryA, baryA+3, _barycenterA);
+
if ( pPolygonA->size() < 3 )
{
if ( !epsilonEqual( triNormal[_ZZ], 0 ))
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
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<nbFaceInit; i++, face_it++);
+
+ for (int i=nbFaceInit; i < (int)_faces.size(); i++)
+ {
+ const auto& face = *(face_it++);
+ _facesBary.emplace_back(std::vector<double>(3, 0.0));
+ _polyNormals.emplace_back(std::vector<double>(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
}
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<double> coords;
+ std::vector<mcIdType> 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<double*>::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<mcIdType,INTERP_KERNEL::ALL_C_MODE>(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;
}
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;
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<double*>::iterator it = polygon.begin() ; it != polygon.end() ; ++it)
+ for(auto it = polygon.begin() ; it != polygon.end() ; ++it)
{
delete[] *it;
*it = 0;
