From 1b2ef0d972bf731593b7d4b966795e0de79e3a7b Mon Sep 17 00:00:00 2001 From: bph Date: Tue, 4 Oct 2011 08:56:36 +0000 Subject: [PATCH] =?utf8?q?Fichiers=20modifi=C3=83=C2=A9s=20pour=20l'optimi?= =?utf8?q?sation?= MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit --- src/INTERP_KERNEL/TransformedTriangle.cxx | 423 +++++++----------- src/INTERP_KERNEL/TransformedTriangle.hxx | 46 +- .../UnitTetraIntersectionBary.cxx | 308 ++++++------- .../UnitTetraIntersectionBary.hxx | 2 + 4 files changed, 313 insertions(+), 466 deletions(-) diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index 68fd7a08b..ef9cd3a34 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -135,10 +135,10 @@ namespace INTERP_KERNEL preCalculateTriangleSurroundsEdge(); preCalculateTripleProducts(); -#if 1 //CS_ALM _indiceA = 0; _indiceB = 0; -#endif + kA=0; + kB=0; } @@ -152,28 +152,10 @@ namespace INTERP_KERNEL { // delete elements of polygons -#if 0 //CS_ALM - for(std::vector::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it) - { - delete[] *it; - } -#else - //_polygonA.clear(); -#endif - - - -#if 0 //CS_ALM - for(std::vector::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it) - { - delete[] *it; - } -#else - //_polygonB.clear(); - _indiceA = 0; _indiceB = 0; -#endif + kA=0; + kB=0; } /** @@ -222,11 +204,7 @@ namespace INTERP_KERNEL // calculate volume under A double volA = 0.0; -#if 0 //CS_ALM - if(_polygonA.size() > 2) -#else if(_indiceA > 2) -#endif { LOG(2, "---- Treating polygon A ... "); calculatePolygonBarycenter(A, barycenter); @@ -237,11 +215,7 @@ namespace INTERP_KERNEL double volB = 0.0; // if triangle is not in h = 0 plane, calculate volume under B -#if 0 //CS_ALM - if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ)) -#else if(_indiceB > 2 && !isTriangleInPlaneOfFacet(XYZ)) -#endif { LOG(2, "---- Treating polygon B ... "); @@ -271,62 +245,37 @@ namespace INTERP_KERNEL */ void TransformedTriangle::calculateIntersectionPolygons() { -#if 1 //CS_ALM assert(_indiceA == 0); assert(_indiceB == 0); -#endif // avoid reallocations in push_back() by pre-allocating enough memory // we should never have more than 20 points -#if 0 //CS_ALM - _polygonA.reserve(20); - _polygonB.reserve(20); -#else - //_polygonA.reserve(60); - //_polygonB.reserve(60); - -#endif // -- surface intersections // surface - edge - - for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1)) { if(testSurfaceEdgeIntersection(edge)) { -#if 0 //CS_ALM - double* ptA = new double[3]; - calcIntersectionPtSurfaceEdge(edge, ptA); - _polygonA.push_back(ptA); - LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A "); -#else double ptA[3] = {0.0,0.0,0.0}; calcIntersectionPtSurfaceEdge(edge, ptA); + for (int i = 0; i < 3; i++) { - _polygonA[_indiceA][i] = ptA[i]; + _polygonA[kA++] = ptA[i]; } + _indiceA += 1; - LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A "); -#endif + LOG(1,"Surface-edge : " << vToStr(ptA) << " added to A "); if(edge >= XY) { -#if 0 //CS_ALM - double* ptB = new double[3]; - copyVector3(ptA, ptB); - _polygonB.push_back(ptB); - LOG(3,"Surface-edge : " << vToStr(ptB) << " added to B "); -#else - double ptB[3]= {0.0,0.0,0.0}; - copyVector3(ptA, ptB); for (int i = 0; i < 3; i++) { - _polygonB[_indiceB][i] = ptB[i]; + _polygonB[kB++] = ptA[i]; } + _indiceB += 1; - LOG(3,"Surface-edge : " << vToStr(ptB) << " added to B "); -#endif + LOG(1,"Surface-edge : " << vToStr(ptA) << " added to B "); } @@ -338,21 +287,13 @@ namespace INTERP_KERNEL { if(testSurfaceRayIntersection(corner)) { -#if 0 //CS_ALM - double* ptB = new double[3]; - copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); - _polygonB.push_back(ptB); - LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B"); -#else - double ptB[3]= {0.0,0.0,0.0} ; - copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); for (int i = 0; i < 3; i++) { + _polygonB[kB++] = COORDS_TET_CORNER[3 * corner+i]; - _polygonB[_indiceB][i] = ptB[i]; } + _indiceB += 1; - LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B"); -#endif + //LOG(1,"Surface-ray : " << vToStr(pt) << " added to B"); } @@ -381,40 +322,26 @@ namespace INTERP_KERNEL if(doTest && testSegmentFacetIntersection(seg, facet)) { -#if 0 //CS_ALM - double* ptA = new double[3]; - calcIntersectionPtSegmentFacet(seg, facet, ptA); - _polygonA.push_back(ptA); - LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A"); -#else double ptA[3] = {0.0,0.0,0.0}; calcIntersectionPtSegmentFacet(seg, facet, ptA); for (int i = 0; i < 3; i++) { - _polygonA[_indiceA][i] = ptA[i]; + _polygonA[kA++] = ptA[i]; } + _indiceA += 1; - LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A"); -#endif + LOG(1,"Segment-facet : " << vToStr(ptA) << " added to A"); if(facet == XYZ) { -#if 0 //CS_ALM - double* ptB = new double[3]; - copyVector3(ptA, ptB); - _polygonB.push_back(ptB); - LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B"); -#else - double ptB[3]= {0.0,0.0,0.0}; - copyVector3(ptA, ptB); for (int i = 0; i < 3; i++) { - _polygonB[_indiceB][i] = ptB[i]; + _polygonB[kB++] = ptA[i]; } + _indiceB += 1; - LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B"); -#endif + LOG(1,"Segment-facet : " << vToStr(ptA) << " added to B"); } @@ -428,38 +355,24 @@ namespace INTERP_KERNEL if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge)) { -#if 0 //CS_ALM - double* ptA = new double[3]; - calcIntersectionPtSegmentEdge(seg, edge, ptA); - _polygonA.push_back(ptA); - LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A"); - if(edge >= XY) - { - double* ptB = new double[3]; - copyVector3(ptA, ptB); - _polygonB.push_back(ptB); - } -#else double ptA[3]= {0.0,0.0,0.0} ; calcIntersectionPtSegmentEdge(seg, edge, ptA); for (int i = 0; i < 3; i++) { - _polygonA[_indiceA][i] = ptA[i]; + _polygonA[kA++] = ptA[i]; } + _indiceA += 1; - LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A"); + LOG(1,"Segment-edge : " << vToStr(ptA) << " added to A"); if(edge >= XY) { - double ptB[3]= {0.0,0.0,0.0}; - copyVector3(ptA, ptB); for (int i = 0; i < 3; i++) { - - _polygonB[_indiceB][i] = ptB[i]; + _polygonB[kB++] = ptA[i]; } + _indiceB += 1; } -#endif } } @@ -473,39 +386,23 @@ namespace INTERP_KERNEL if(doTest && testSegmentCornerIntersection(seg, corner)) { -#if 0 //CS_ALM - double* ptA = new double[3]; - copyVector3(&COORDS_TET_CORNER[3 * corner], ptA); - _polygonA.push_back(ptA); - LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A"); - if(corner != O) - { - double* ptB = new double[3]; - _polygonB.push_back(ptB); - copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); - LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B"); - } -#else - double ptA[3]= {0.0,0.0,0.0}; - copyVector3(&COORDS_TET_CORNER[3 * corner], ptA); for (int i = 0; i < 3; i++) { + _polygonA[kA++] = COORDS_TET_CORNER[3 * corner+i]; - _polygonA[_indiceA][i] = ptA[i]; } _indiceA += 1; - LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A"); + + LOG(1,"Segment-corner : " << vToStr(ptA) << " added to A"); if(corner != O) { - double ptB[3]= {0.0,0.0,0.0}; - - copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); for (int i = 0; i < 3; i++) { - _polygonB[_indiceB][i] = ptB[i]; + _polygonB[kB++] = COORDS_TET_CORNER[3 * corner+i]; + } + _indiceB += 1; - LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B"); + LOG(1,"Segment-corner : " << vToStr(ptB) << " added to B"); } -#endif } } @@ -514,21 +411,13 @@ namespace INTERP_KERNEL { if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner)) { -#if 0 //CS_ALM - double* ptB = new double[3]; - copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); - _polygonB.push_back(ptB); - LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B"); -#else - double ptB[3]= {0.0,0.0,0.0}; - copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); for (int i = 0; i < 3; i++) { - - _polygonB[_indiceB][i] = ptB[i]; + _polygonB[kB++] = COORDS_TET_CORNER[3 * corner+i]; + } + _indiceB += 1; - LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B"); -#endif + LOG(1,"Segment-ray : " << vToStr(ptB) << " added to B"); } } @@ -547,20 +436,14 @@ namespace INTERP_KERNEL #endif if(testSegmentHalfstripIntersection(seg, edge)) { -#if 0 //CS_ALM - double* ptB = new double[3]; - calcIntersectionPtSegmentHalfstrip(seg, edge, ptB); - _polygonB.push_back(ptB); - LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B"); -#else double ptB[3]= {0.0,0.0,0.0}; calcIntersectionPtSegmentHalfstrip(seg, edge, ptB); for (int i = 0; i < 3; i++) { - _polygonB[_indiceB][i] = ptB[i]; + _polygonB[kB++] = ptB[i]; } + _indiceB += 1; - LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B"); -#endif + LOG(1,"Segment-halfstrip : " << vToStr(ptB) << " added to B"); } } } @@ -572,68 +455,53 @@ namespace INTERP_KERNEL // tetrahedron if(testCornerInTetrahedron(corner)) { -#if 0 //CS_ALM - double* ptA = new double[3]; - copyVector3(&_coords[5*corner], ptA); - _polygonA.push_back(ptA); - LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A"); -#else - double ptA[3]= {0.0,0.0,0.0}; - copyVector3(&_coords[5*corner], ptA); - for (int i = 0; i < 3; i++) { - _polygonA[_indiceA][i] = ptA[i]; + for (int i = 0; i < 3; i++) { + _polygonA[kA++] = _coords[5*corner+i]; + } + _indiceA += 1; - LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A"); -#endif - } + LOG(1,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A"); + } // XYZ - plane if(testCornerOnXYZFacet(corner)) { -#if 0 //CS_ALM - double* ptB = new double[3]; - copyVector3(&_coords[5*corner], ptB); - _polygonB.push_back(ptB); - LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B"); -#else - double ptB[3]= {0.0,0.0,0.0}; - copyVector3(&_coords[5*corner], ptB); - for (int i = 0; i < 3; i++) { - - _polygonB[_indiceB][i] = ptB[i]; + for (int i = 0; i < 3; i++) { + _polygonB[kB++] = _coords[5*corner+i]; + } + _indiceB += 1; - LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B"); -#endif + LOG(1,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B"); } // projection on XYZ - facet if(testCornerAboveXYZFacet(corner)) { -#if 0 //CS_ALM - double* ptB = new double[3]; - copyVector3(&_coords[5*corner], ptB); - ptB[2] = 1 - ptB[0] - ptB[1]; - assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0)); - _polygonB.push_back(ptB); - LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B"); -#else double ptB[3]= {0.0,0.0,0.0}; copyVector3(&_coords[5*corner], ptB); ptB[2] = 1 - ptB[0] - ptB[1]; assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0)); for (int i = 0; i < 3; i++) { - _polygonB[_indiceB][i] = ptB[i]; + _polygonB[kB++] = ptB[i]; } + _indiceB += 1; - LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B"); -#endif + LOG(1,"Projection XYZ-plane : " << vToStr(ptB) << " added to B"); } } + for (int i = 0; i < _indiceA; i++) { + _polygonA2[i] = &_polygonA[i*3]; + } + + for (int i = 0; i < _indiceB; i++) { + _polygonB2[i] = &_polygonB[i*3]; + } + } /** @@ -645,44 +513,37 @@ namespace INTERP_KERNEL * @param barycenter array of three doubles where barycenter is stored * */ -#if 0 //CS_ALM void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter) -#else - void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter) -#endif { LOG(3,"--- Calculating polygon barycenter"); // get the polygon points -#if 0 //CS_ALM - std::vector& polygon = (poly == A) ? _polygonA : _polygonB; -#else - VEC3 polygon[20]; + double polygon[20*3]; + double *polygon2[20]; if (poly == A) { - for (int i=0; i < _indiceA; i++) { - for (int j = 0; j < 3; j++) { - polygon[i][j] = _polygonA[i][j]; - } + for (int i=0; i < _indiceA*3; i++) { + polygon[i] = _polygonA[i]; + } + for (int i=0; i < _indiceA; i++) { + polygon2[i] = &polygon[i*3]; + } } else { - for (int i=0; i < _indiceB; i++) { - for (int j = 0; j < 3; j++) { - polygon[i][j] = _polygonB[i][j]; - } + for (int i=0; i < _indiceB*3; i++) { + polygon[i] = _polygonB[i]; + } - } - -#endif + for (int i=0; i < _indiceB; i++) { + polygon2[i] = &polygon[i*3]; + } + } + + // calculate barycenter -#if 0 //CS_ALM - const int m = polygon.size(); -#else const int m = (poly == A) ? _indiceA : _indiceB ; -#endif - for(int j = 0 ; j < 3 ; ++j) { barycenter[j] = 0.0; @@ -692,12 +553,13 @@ namespace INTERP_KERNEL { for(int i = 0 ; i < m ; ++i) { - - const double* pt = polygon[i]; - for(int j = 0 ; j < 3 ; ++j) - { - barycenter[j] += pt[j] / double(m); - } + + const double* pt = polygon2[i]; + for(int j = 0 ; j < 3 ; ++j) + { + barycenter[j] += pt[j] / double(m); + } + } } LOG(3,"Barycenter is " << vToStr(barycenter)); @@ -722,32 +584,38 @@ namespace INTERP_KERNEL typedef SortOrder::CoordType CoordType; // get the polygon points -#if 0 //CS_ALM - std::vector& polygon = (poly == A) ? _polygonA : _polygonB; - if(polygon.size() == 0) - return; -#else - VEC3 polygon[20]; + double* polygon2[20]; + if (poly == A) { + /*for (int i=0; i < _indiceA*3; i++) { + polygon[i] = _polygonA[i]; + + }*/ for (int i=0; i < _indiceA; i++) { - for (int j = 0; j < 3; j++) { - polygon[i][j] = _polygonA[i][j]; - } - } + polygon2[i] = new double[3]; + polygon2[i][0] = _polygonA[i*3]; + polygon2[i][1] = _polygonA[i*3 + 1]; + polygon2[i][2] = _polygonA[i*3 + 2]; + // polygon2[i] = &_polygonA[i*3]; + } } else { + /*for (int i=0; i < _indiceB*3; i++) { + + polygon[i] = _polygonB[i]; + }*/ for (int i=0; i < _indiceB; i++) { - for (int j = 0; j < 3; j++) { - polygon[i][j] = _polygonB[i][j]; - } + // polygon2[i] = &_polygonB[i*3]; + polygon2[i] = new double[3]; + polygon2[i][0] = _polygonB[i*3]; + polygon2[i][1] = _polygonB[i*3 + 1]; + polygon2[i][2] = _polygonB[i*3 + 2]; } } int taillePoly = (poly == A) ? _indiceA : _indiceB; if (taillePoly == 0) return; -#endif - // determine type of sorting CoordType type = SortOrder::XY; @@ -766,7 +634,7 @@ namespace INTERP_KERNEL { type = SortOrder::YZ; } - } + } // create order object SortOrder order(barycenter, type); @@ -774,24 +642,43 @@ namespace INTERP_KERNEL // sort vector with this object // NB : do not change place of first object, with respect to which the order // is defined -#if 0 //CS_ALM - sort((polygon.begin()), polygon.end(), order); -#else - //double* pol; - //std::sort( pol, pol+60, order); - std::sort( polygon, polygon+(taillePoly*3), order); -#endif + std::sort( polygon2, polygon2+taillePoly, order); + if (poly == A) { + for (int i=0; i < _indiceA; i++) { + // _polygonA[i*3] = *polygon2[i]; + _polygonA[i*3] = polygon2[i][0]; + _polygonA[i*3 + 1] = polygon2[i][1]; + _polygonA[i*3 + 2] = polygon2[i][2]; + + } + } + else { + for (int i=0; i < _indiceB; i++) { + // _polygonB[i*3] = *polygon2[i]; + _polygonB[i*3] = polygon2[i][0]; + _polygonB[i*3 + 1] = polygon2[i][1]; + _polygonB[i*3 + 2] = polygon2[i][2]; + + } + } + LOG(3,"Sorted polygon is "); -#if 0 //CS_ALM - for(size_t i = 0 ; i < polygon.size() ; ++i) -#else for(size_t i = 0 ; i < (size_t)taillePoly ; ++i) -#endif + { - LOG(3,vToStr(polygon[i])); + if (poly == A) { + LOG(1,vToStr(&_polygonA[i*3])); + } + else { + LOG(1,vToStr(&_polygonB[i*3])); + } } +// for (int i = 0; i < taillePoly; i++) { +// delete [] polygon2[i]; +// } +// delete polygon2; } /** @@ -810,42 +697,33 @@ namespace INTERP_KERNEL LOG(2,"--- Calculating volume under polygon"); // get the polygon points -#if 0 //CS_ALM - std::vector& polygon = (poly == A) ? _polygonA : _polygonB; - int taillePoly = polygon.size() -#else - VEC3 polygon[20]; - if (poly == A) { + + double *polygon2[20]; + if (poly == A) { + for (int i=0; i < _indiceA; i++) { - for (int j = 0; j < 3; j++) { - polygon[i][j] = _polygonA[i][j]; - } - } + polygon2[i] = &_polygonA[i*3]; + } } else { + for (int i=0; i < _indiceB; i++) { - for (int j = 0; j < 3; j++) { - polygon[i][j] = _polygonB[i][j]; - } - } - } + polygon2[i] = &_polygonB[i*3]; + } + } int taillePoly = (poly == A) ? _indiceA : _indiceB; - -#endif double vol = 0.0; const int m = taillePoly; for(int i = 0 ; i < m ; ++i) { - - const double* ptCurr = polygon[i]; // pt "i" - const double* ptNext = polygon[(i + 1) % m]; // pt "i+1" (pt m == pt 0) - + const double* ptCurr = polygon2[i]; // pt "i" + const double* ptNext = polygon2[(i + 1) % m]; // pt "i+1" (pt m == pt 0) const double factor1 = ptCurr[2] + ptNext[2] + barycenter[2]; const double factor2 = ptCurr[0]*(ptNext[1] - barycenter[1]) @@ -853,6 +731,7 @@ namespace INTERP_KERNEL + barycenter[0]*(ptCurr[1] - ptNext[1]); vol += (factor1 * factor2) / 6.0; } + LOG(2,"Abs. Volume is " << vol); return vol; diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 2f9c99650..1dd749eb4 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -99,7 +99,7 @@ namespace INTERP_KERNEL */ class INTERPKERNEL_EXPORT TransformedTriangle { - + public: @@ -132,11 +132,7 @@ namespace INTERP_KERNEL /// Double products /// NB : order corresponds to TetraEdges (Grandy, table III) enum DoubleProduct { C_YZ = 0, C_ZX, C_XY, C_ZH, C_XH, C_YH, C_01, C_10, NO_DP }; -#if 1 //CS_ALM - /// - typedef double VEC3[3]; - typedef VEC3 POLYGON_TYPE[20]; -#endif + TransformedTriangle(double* p, double* q, double* r); ~TransformedTriangle(); @@ -147,12 +143,13 @@ namespace INTERP_KERNEL // Queries of member values used by UnitTetraIntersectionBary const double* getCorner(TriCorner corner) const { return _coords + 5*corner; } -#if 0 //CS_ALM - const std::vector& getPolygonA() const { return _polygonA; } -#else - //const std::vector& getPolygonA() const { return _polygonA; } - const POLYGON_TYPE& getPolygonA() const { return _polygonA; } -#endif + + + const double* getPolygonA() const { return _polygonA; } + int getIndiceA() const {return _indiceA;} + + int getIndiceB() const {return _indiceB;} + double getVolume() const { return _volume; } @@ -287,27 +284,24 @@ namespace INTERP_KERNEL /// Vector holding the points of the intersection polygon A. /// these points are allocated in calculateIntersectionPolygons() and liberated in the destructor -#if 0 //CS_ALM - std::vector _polygonA; -#else - //double _polygonA[60]; - //std::vector _polygonA; - VEC3 _polygonA[20]; + + double _polygonA[20*3]; + double *_polygonA2[20]; int _indiceA; int _indiceB; + int kA; + int kB; + -#endif /// Vector holding the points of the intersection polygon B. /// These points are allocated in calculateIntersectionPolygons() and liberated in the destructor -#if 0 //CS_ALM - std::vector _polygonB; -#else - //double _polygonB[60]; - //std::vector _polygonB; - VEC3 _polygonB[20]; -#endif + + double _polygonB[20*3]; + double *_polygonB2[20]; + + /// Array holding the coordinates of the barycenter of the polygon A /// This point is calculated in calculatePolygonBarycenter diff --git a/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx b/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx index 9f2df0e4f..96ceb753d 100644 --- a/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx +++ b/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx @@ -83,12 +83,11 @@ namespace INTERP_KERNEL double triNormal[3], polyNormal[3]; crossprod<3>( triangle.getCorner(P),triangle.getCorner(Q),triangle.getCorner(R), triNormal); -#if 0 //CS_ALM - const std::vector * pPolygonA = &triangle.getPolygonA(); -#else - const std::vector * pPolygonA = &triangle.getPolygonA(); -#endif - if ( pPolygonA->size() < 3 ) + + const double* pPolygonA = triangle.getPolygonA(); + + + if ( triangle.getIndiceA() < 3 ) { if ( !epsilonEqual( triNormal[_Z], 0 )) return; // not vertical triangle does not intersect the unit tetra @@ -96,32 +95,35 @@ namespace INTERP_KERNEL // Vertical triangle. Use inherited methods of TransformedTriangle to // calculate intesection polygon *((TransformedTriangle*)this) = triangle; // copy triangle fields - _polygonA.clear(); - _polygonB.clear(); + this->_indiceA = 0; + this->_indiceB= 0; + this->kA = 0 ; + this->kB = 0; + calculateIntersectionPolygons(); - if (this->_polygonA.size() < 3) + if (this->_indiceA < 3) return; calculatePolygonBarycenter(A, _barycenterA); sortIntersectionPolygon(A, _barycenterA); - pPolygonA = & _polygonA; + + pPolygonA = _polygonA; + } + // check if polygon orientation is same as the one of triangle -#if 0 //CS_ALM - std::vector::const_iterator p = pPolygonA->begin(), pEnd = pPolygonA->end(); -#else - std::vector::const_iterator p = pPolygonA->begin(), pEnd = pPolygonA->end(); - /*int taillePoly = pPolygonA->size(); - double p[3]; - double pEnd[3]; - - for (int j = 0; j < 3; j++){ - p[j] = pPolygonA->operator[](j); - }*/ - -#endif - //#ifdef DMP_UNITTETRAINTERSECTIONBARY + for (int i = 0; i < _indiceA; i ++) + _polygonA2[i] = &_polygonA[i*3]; + + for (int i = 0; i < _indiceB; i ++) + _polygonB2[i] = &_polygonB[i*3]; + + const double* p = _polygonA2[0]; + const double* pEnd = _polygonA2[_indiceA]; + + +#ifdef DMP_UNITTETRAINTERSECTIONBARY std::cout.precision(18); std::cout << "**** int polygon() " << std::endl; while ( p != pEnd ) @@ -131,17 +133,13 @@ namespace INTERP_KERNEL p += 3; } p = pPolygonA->begin(); - //#endif +#endif -#if 0 //CS_ALM - double* p1 = *p++; - double* p2 = *p; - while ( samePoint( p1, p2 ) && ++p != pEnd ) - p2 = *p; -#else - -#endif + const double* p1 = p++; + const double* p2 = p; + while ( samePoint( p1, p2 ) && ++p != pEnd ) + p2 = p; if ( p == pEnd ) { @@ -151,9 +149,11 @@ namespace INTERP_KERNEL clearPolygons(); return; } - double* p3 = *p; + + const double* p3 = p; while (( samePoint( p2, p3 ) || samePoint( p1, p3 )) && ++p != pEnd ) - p3 = *p; + p3 = p; + if ( p == pEnd ) { #ifdef DMP_UNITTETRAINTERSECTIONBARY @@ -171,22 +171,25 @@ namespace INTERP_KERNEL _faces.push_back( std::vector< double* > () ); std::vector< double* >& faceCorner = _faces.back(); - faceCorner.resize( pPolygonA->size()/* + 1*/ ); + faceCorner.resize( _indiceA/* + 1*/ ); int i = 0; if ( reverse ) { - std::vector::const_reverse_iterator polyF = pPolygonA->rbegin(), polyEnd; - for ( polyEnd = pPolygonA->rend(); polyF != polyEnd; ++i, ++polyF ) - if ( i==0 || !samePoint( *polyF, faceCorner[i-1] )) { -#if 0 //CS_ALM - copyVector3( *polyF, faceCorner[i] = new double[3] ); -#else + + const double* polyF = _polygonA2[_indiceA]; + const double * polyEnd; + //for ( polyEnd = _polygonA2[0]; polyF != polyEnd; ++i, ++polyF ) + for ( polyEnd = _polygonA2[0]; polyF != polyEnd; ++i, --polyF ) + if ( i==0 || !samePoint( polyF, faceCorner[i-1] )) { + + + double var[3]= {0.0,0.0,0.0}; - faceCorner[i] = var; - copyVector3( *polyF, faceCorner[i] ); - //faceCorner[i] = var; -#endif + copyVector3(faceCorner[i] ,var); + copyVector3( polyF, faceCorner[i] ); + + } else { --i; } @@ -196,16 +199,17 @@ namespace INTERP_KERNEL } else { - std::vector::const_iterator polyF = pPolygonA->begin(), polyEnd; - for ( polyEnd = pPolygonA->end(); polyF != polyEnd; ++i, ++polyF ) - if ( i==0 || !samePoint( *polyF, faceCorner[i-1] )) { -#if 0 //CS_ALM - copyVector3( *polyF, faceCorner[i] = new double[3] ); -#else + + //const double* polyF = pPolygonA; + const double* polyF = _polygonA2[0]; + const double * polyEnd; + for ( polyEnd = _polygonA2[_indiceA]; polyF != polyEnd; ++i, ++polyF ) + if ( i==0 || !samePoint( polyF, faceCorner[i-1] )) { + double var[3]= {0.0,0.0,0.0}; - faceCorner[i] = var; - copyVector3( *polyF, faceCorner[i] ); -#endif + copyVector3(faceCorner[i] ,var); + copyVector3( polyF, faceCorner[i] ); + } else { --i; } @@ -213,12 +217,22 @@ namespace INTERP_KERNEL if ( i < 3 ) { clearPolygons(); // free memory of _polygonA - _polygonA = faceCorner; + + int facesize = faceCorner.size(); + int k = 0; + for (int i = 0; i < facesize*3; i++,k++) { + for (int j = 0; j < 3; j++) { + _polygonA[k++] = faceCorner[i][j]; + } + + } + _faces.pop_back(); } else { - if ( i < (int)pPolygonA->size() ) + + if ( i < (int)_indiceA ) faceCorner.resize( i ); if ( _polyNormals.empty() ) @@ -291,9 +305,9 @@ namespace INTERP_KERNEL double bary[] = { 0, 0, 0 }; // base polygon bary -#if 1 //CS_ALM + int polySize = polygon.size(); -#endif + for ( v = polygon.begin(); v != vEnd ; ++v ) { double* p = *v; @@ -301,30 +315,22 @@ namespace INTERP_KERNEL bary[1] += p[1]; bary[2] += p[2]; } -#if 0 //CS_ALM - bary[0] /= polygon.size(); - bary[1] /= polygon.size(); - bary[2] /= polygon.size(); -#else + bary[0] /= polySize; bary[1] /= polySize; bary[2] /= polySize; -#endif + // pyramid volume double vol = 0; -#if 0 //CS_ALM - for ( int i = 0; i < (int)polygon.size(); ++i ) -#else + for ( int i = 0; i < polySize; ++i ) -#endif + { double* p1 = polygon[i]; -#if 0 //CS_ALM - double* p2 = polygon[(i+1)%polygon.size()]; -#else + double* p2 = polygon[(i+1)%polySize]; -#endif + vol += std::fabs( calculateVolumeForTetra( p1, p2, bary, P )); } @@ -374,15 +380,10 @@ namespace INTERP_KERNEL { std::vector< double* >& polygon = *f; double coordSum[3] = {0,0,0}; -#if 1 //CS_ALM + int polySize = polygon.size(); -#endif -#if 0 //CS_ALM - for ( int i = 0; i < (int)polygon.size(); ++i ) -#else for ( int i = 0; i < polySize; ++i ) -#endif { double* p = polygon[i]; coordSum[0] += p[0]; @@ -394,13 +395,9 @@ namespace INTERP_KERNEL if ( epsilonEqual( coordSum[j], 0.0 )) sideAdded[j] = ++nbAddedSides != 0 ; } -#if 0 //CS_ALM - if ( !sideAdded[3] && - ( epsilonEqual( (coordSum[0]+coordSum[1]+coordSum[2]) / polygon.size(), 1. ))) { -#else + if ( !sideAdded[3] && ( epsilonEqual( (coordSum[0]+coordSum[1]+coordSum[2]) / polySize, 1. ))) { -#endif sideAdded[3] = ++nbAddedSides != 0 ; } } @@ -412,11 +409,9 @@ namespace INTERP_KERNEL // --------------------------------------------------------------------------------- int nbIntersectPolygs = _faces.size(); -#if 1 //CS_ALM + std::vector< double* > * sideFaces[ 4 ]; // future polygons on sides of tetra -#else - std::vector< double[3] > * sideFaces[ 4 ]; // future polygons on sides of tetra -#endif + for ( int i = 0; i < NB_TETRA_SIDES; ++i ) { sideFaces[ i ]=0; @@ -430,18 +425,16 @@ namespace INTERP_KERNEL for ( int iF = 0; iF < nbIntersectPolygs; ++f, ++iF ) // loop on added intersection polygons { std::vector< double* >& polygon = *f; -#if 1 //CS_ALM + int polySize = (int)polygon.size(); -#endif + for ( int i = 0; i < polySize; ++i ) { // segment ends double* p1 = polygon[i]; -#if 0 //CS_ALM - double* p2 = polygon[(i+1)%polygon.size()]; -#else + double* p2 = polygon[(i+1)%polySize]; -#endif + bool p1OnSide, p2OnSide;//, onZeroSide = false; for ( int j = 0; j < 3; ++j ) { @@ -452,19 +445,14 @@ namespace INTERP_KERNEL if ( p1OnSide && p2OnSide ) { // segment p1-p2 is on j-th orthogonal side of tetra -#if 0 //CS_ALM - sideFaces[j]->push_back( new double[3] ); - copyVector3( p1, sideFaces[j]->back() ); - sideFaces[j]->push_back( new double[3] ); - copyVector3( p2, sideFaces[j]->back() ); -#else + double var[3]= {0.0,0.0,0.0}; sideFaces[j]->push_back( var ); copyVector3( p1, sideFaces[j]->back() ); double var2[3]= {0.0,0.0,0.0}; sideFaces[j]->push_back( var2 ); copyVector3( p2, sideFaces[j]->back() ); -#endif + //break; } @@ -474,19 +462,14 @@ namespace INTERP_KERNEL epsilonEqual( p1[_X] + p1[_Y] + p1[_Z], 1.0 ) && epsilonEqual( p2[_X] + p2[_Y] + p2[_Z], 1.0 )) { -#if 0 //CS_ALM - sideFaces[3]->push_back( new double[3] ); - copyVector3( p1, sideFaces[3]->back() ); - sideFaces[3]->push_back( new double[3] ); - copyVector3( p2, sideFaces[3]->back() ); -#else + double var[3]= {0.0,0.0,0.0}; sideFaces[3]->push_back( var ); copyVector3( p1, sideFaces[3]->back() ); double var2[3]= {0.0,0.0,0.0}; sideFaces[3]->push_back( var2 ); copyVector3( p2, sideFaces[3]->back() ); -#endif + } } } @@ -502,14 +485,10 @@ namespace INTERP_KERNEL else { std::vector< double* >& sideFace = *sideFaces[i]; -#if 0 //CS_ALM - for ( int i = 0; i < sideFace.size(); ++i ) - { -#else + int faceSize = sideFace.size(); for ( int i = 0; i < faceSize; ++i ) { -#endif double* p = sideFace[i]; std::cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << std::endl; @@ -696,12 +675,10 @@ namespace INTERP_KERNEL { if ( !cutOffCorners[ ic ] && ic != excludeCorner ) { -#if 0 //CS_ALM - sideFace.push_back( new double[3] ); -#else + double var[3]= {0.0,0.0,0.0}; sideFace.push_back( var ); -#endif + copyVector3( origin, sideFace.back() ); if ( ic ) sideFace.back()[ ic-1 ] = 1.0; @@ -747,9 +724,18 @@ namespace INTERP_KERNEL if ( face.size() >= 3 ) { clearPolygons(); // free memory of _polygonA - _polygonA = face; + + int k = 0; + for (unsigned int i = 0; i < face.size(); i++,k++) { + for (int j = 0; j < 3; j++) { + _polygonA[k++] = face[i][j]; + } + _indiceA++; + } face.clear(); - face.reserve( _polygonA.size() ); + face.reserve( _indiceA ); + + if ( iF >= nbIntersectPolygs ) { // sort points of side faces calculatePolygonBarycenter( A, _barycenterA ); @@ -757,25 +743,42 @@ namespace INTERP_KERNEL sortIntersectionPolygon( A, _barycenterA ); } // exclude equal points - std::vector< double* >::iterator v = _polygonA.begin(), vEnd = _polygonA.end(); - face.push_back( *v ); - *v = 0; - for ( ++v; v != vEnd; ++v ) + + for (int i = 0; i < _indiceA; i++) { + _polygonA2[i] = &_polygonA[i*3]; + } + double* v = _polygonA2[0]; + double *vEnd = _polygonA2[_indiceA]; + face.push_back( v ); + + + + + *v = 0; + for ( ++v; v != vEnd; ++v ) { - double* pPrev = face.back(); - double* p = *v; + double* pPrev = face.back(); + double* p = v; if ( !samePoint( p, pPrev )) - { - face.push_back( p ); - *v = 0; - } - } + { + face.push_back( p ); + *v = 0; + } + } + } if ( face.size() < 3 ) { // size could decrease clearPolygons(); // free memory of _polygonA - _polygonA = face; - face.clear(); + + int k = 0; + for (unsigned int i = 0; i < face.size(); i++,k++) { + for (int j = 0; j < 3; j++) { + _polygonA[k++] = face[i][j]; + } + } + face.clear(); + } else { @@ -826,44 +829,13 @@ namespace INTERP_KERNEL void UnitTetraIntersectionBary::clearPolygons(bool andFaces) { -#if 0 //CS_ALM - for(std::vector::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it) - { - delete[] *it; - *it = 0; - } - - for(std::vector::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it) - { - - delete[] *it; - - *it = 0; - } -#endif - - _polygonA.clear(); - _polygonB.clear(); + _indiceA=0; + _indiceB=0; + kA=0; + kB=0; - if ( andFaces ) - { - std::list< std::vector< double* > >::iterator f = this->_faces.begin(), fEnd = this->_faces.end(); -#if 0 //CS_ALM - for ( ; f != fEnd; ++f ) - { - std::vector< double* >& polygon = *f; - for(std::vector::iterator it = polygon.begin() ; it != polygon.end() ; ++it) - { - - delete[] *it; - *it = 0; - } - } -#endif - this->_faces.clear(); - } } //================================================================================ @@ -874,7 +846,7 @@ namespace INTERP_KERNEL UnitTetraIntersectionBary::~UnitTetraIntersectionBary() { - clearPolygons(/*andFaces=*/true ); + //clearPolygons(/*andFaces=*/true ); } } diff --git a/src/INTERP_KERNEL/UnitTetraIntersectionBary.hxx b/src/INTERP_KERNEL/UnitTetraIntersectionBary.hxx index 726dfdef7..064095d24 100644 --- a/src/INTERP_KERNEL/UnitTetraIntersectionBary.hxx +++ b/src/INTERP_KERNEL/UnitTetraIntersectionBary.hxx @@ -32,7 +32,9 @@ namespace INTERP_KERNEL { + class UnitTetraIntersectionBary : protected TransformedTriangle + { public: INTERPKERNEL_EXPORT UnitTetraIntersectionBary(bool isTetraInversed=false); -- 2.39.2