From f9cdf861759639bfd122e8246c43ef9514e47daf Mon Sep 17 00:00:00 2001 From: bph Date: Fri, 5 Aug 2011 14:50:09 +0000 Subject: [PATCH] *** empty log message *** --- src/INTERP_KERNEL/CurveIntersectorP0P1.txx | 11 +- .../PolyhedronIntersectorP0P0.txx | 6 +- .../PolyhedronIntersectorP0P1.txx | 23 ++ .../PolyhedronIntersectorP1P0.txx | 24 +- .../PolyhedronIntersectorP1P0Bary.txx | 24 +- src/INTERP_KERNEL/SplitterTetra.hxx | 39 ++- src/INTERP_KERNEL/SplitterTetra.txx | 138 ++++++-- src/INTERP_KERNEL/TransformedTriangle.cxx | 300 +++++++++++++++++- src/INTERP_KERNEL/TransformedTriangle.hxx | 29 +- .../TransformedTriangleIntersect.cxx | 1 + .../UnitTetraIntersectionBary.cxx | 151 ++++++++- 11 files changed, 699 insertions(+), 47 deletions(-) diff --git a/src/INTERP_KERNEL/CurveIntersectorP0P1.txx b/src/INTERP_KERNEL/CurveIntersectorP0P1.txx index d6bed1d79..803fe02d2 100644 --- a/src/INTERP_KERNEL/CurveIntersectorP0P1.txx +++ b/src/INTERP_KERNEL/CurveIntersectorP0P1.txx @@ -60,9 +60,18 @@ namespace INTERP_KERNEL ::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) { std::vector segmentsT; - BASE_INTERSECTOR::getDualSegments( icellT, BASE_INTERSECTOR::_meshT, segmentsT); + BASE_INTERSECTOR::getDualSegments( icellT, + BASE_INTERSECTOR::_meshT, segmentsT); +#if 1 //CS_ALM + int taille = (int)segmentsT.size(); +#endif +#if 0 //CS_ALM for ( int t = 0; t < (int)segmentsT.size(); ++t ) { +#else + for ( int t = 0; t < taille; ++t ) + { +#endif typename MyMatrix::value_type& resRow = res[ OTT::ind2C( segmentsT[t]._nodeId )]; for(typename std::vector::const_iterator iter=icellsS.begin(); iter!=icellsS.end(); iter++) diff --git a/src/INTERP_KERNEL/PolyhedronIntersectorP0P0.txx b/src/INTERP_KERNEL/PolyhedronIntersectorP0P0.txx index 260a4c110..40ca69aba 100644 --- a/src/INTERP_KERNEL/PolyhedronIntersectorP0P0.txx +++ b/src/INTERP_KERNEL/PolyhedronIntersectorP0P0.txx @@ -82,11 +82,13 @@ namespace INTERP_KERNEL for(typename std::vector::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++) { double volume = 0.; - for(typename std::vector*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) + for(typename std::vector*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) + { volume += (*iter)->intersectSourceCell(*iterCellS); + } if(volume!=0.) res[targetCell].insert(std::make_pair(OTT::indFC(*iterCellS), volume)); - } + } _split.releaseArrays(); } } diff --git a/src/INTERP_KERNEL/PolyhedronIntersectorP0P1.txx b/src/INTERP_KERNEL/PolyhedronIntersectorP0P1.txx index 4dbf72dc9..c7dff8a71 100644 --- a/src/INTERP_KERNEL/PolyhedronIntersectorP0P1.txx +++ b/src/INTERP_KERNEL/PolyhedronIntersectorP0P1.txx @@ -80,15 +80,38 @@ namespace INTERP_KERNEL int nbOfNodesT=Intersector3D::_target_mesh.getNumberOfNodesOfElement(OTT::indFC(targetCell)); releaseArrays(); _split.splitTargetCell(targetCell,nbOfNodesT,_tetra); + + for(typename std::vector::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++) { +#if 0 //CS_ALM + typename MyMeshType::MyConnType cellSrc = *iterCellS; + const MyMeshType& _src_mesh = Intersector3DP0P1::_src_mesh; + //int cellSrcIdx = OTT::indFC(cellSrc); + NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT::indFC(cellSrc)); + const CellModel& cellModelCell=CellModel::GetCellModel(normCellType); + unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT::indFC(cellSrc)) : cellModelCell.getNumberOfNodes(); + // halfspace filtering + bool isOutside[8] = {true, true, true, true, true, true, true, true}; + bool isTargetOutside = false; + + // calculate the coordinates of the nodes + int *cellNodes=new int[nbOfNodes4Type]; + +#endif for(typename std::vector*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) { (*iter)->splitIntoDualCells(subTetras); + for(int i=0;i<24;i++) { SplitterTetra *tmp=subTetras[i]; +#if 1 //CS_ALM double volume = tmp->intersectSourceCell(*iterCellS); +#else + double volume = tmp->intersectSourceCell(*iterCellS,normCellType, cellModelCell, nbOfNodes4Type,cellNodes, isOutside, + isTargetOutside); +#endif if(volume!=0.) { typename MyMatrix::value_type& resRow=res[tmp->getId(0)]; diff --git a/src/INTERP_KERNEL/PolyhedronIntersectorP1P0.txx b/src/INTERP_KERNEL/PolyhedronIntersectorP1P0.txx index 43ca2840a..3bef18613 100644 --- a/src/INTERP_KERNEL/PolyhedronIntersectorP1P0.txx +++ b/src/INTERP_KERNEL/PolyhedronIntersectorP1P0.txx @@ -83,16 +83,38 @@ namespace INTERP_KERNEL typename MyMatrix::value_type& resRow=res[targetCell]; for(typename std::vector::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++) { + releaseArrays(); - int nbOfNodesS=Intersector3D::_src_mesh.getNumberOfNodesOfElement(OTT::indFC(*iterCellS)); + int nbOfNodesS=Intersector3DP1P0::_src_mesh.getNumberOfNodesOfElement(OTT::indFC(*iterCellS)); _split.splitTargetCell(*iterCellS,nbOfNodesS,_tetra); +#if 0 //CS_ALM + typename MyMeshType::MyConnType cellSrc = *iterCellS; + const MyMeshType& _src_mesh = Intersector3DP1P0::_src_mesh; + //int cellSrcIdx = OTT::indFC(cellSrc); + NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT::indFC(cellSrc)); + const CellModel& cellModelCell=CellModel::GetCellModel(normCellType); + unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT::indFC(cellSrc)) : cellModelCell.getNumberOfNodes(); + // halfspace filtering + bool isOutside[8] = {true, true, true, true, true, true, true, true}; + bool isTargetOutside = false; + + // calculate the coordinates of the nodes + int *cellNodes=new int[nbOfNodes4Type]; + +#endif + for(typename std::vector*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) { (*iter)->splitIntoDualCells(subTetras); for(int i=0;i<24;i++) { SplitterTetra *tmp=subTetras[i]; +#if 1 //CS_ALM double volume = tmp->intersectSourceCell(targetCell); +#else + double volume = tmp->intersectSourceCell(targetCell, normCellType, cellModelCell, nbOfNodes4Type,cellNodes, isOutside, + isTargetOutside); +#endif ConnType sourceNode=tmp->getId(0); if(volume!=0.) { diff --git a/src/INTERP_KERNEL/PolyhedronIntersectorP1P0Bary.txx b/src/INTERP_KERNEL/PolyhedronIntersectorP1P0Bary.txx index d3ae2560f..a40a9ae0b 100644 --- a/src/INTERP_KERNEL/PolyhedronIntersectorP1P0Bary.txx +++ b/src/INTERP_KERNEL/PolyhedronIntersectorP1P0Bary.txx @@ -105,11 +105,33 @@ namespace INTERP_KERNEL // intersect a source tetrahedron with each target tetrahedron: get intersection volume and barycenter double baryCentre[SPACEDIM], total_baryCentre[3] = { 0., 0., 0.}; double interVolume = 0; +#if 0 //CS_ALM + typename MyMeshType::MyConnType cellSrc = *iterCellS; + //int cellSrcIdx = OTT::indFC(cellSrc); + const MyMeshType& _src_mesh = Intersector3DP1P0Bary::_src_mesh; + NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT::indFC(cellSrc)); + const CellModel& cellModelCell=CellModel::GetCellModel(normCellType); + unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT::indFC(cellSrc)) : cellModelCell.getNumberOfNodes(); + // halfspace filtering + bool isOutside[8] = {true, true, true, true, true, true, true, true}; + bool isTargetOutside = false; + + // calculate the coordinates of the nodes + int *cellNodes=new int[nbOfNodes4Type]; + +#endif for(typename std::vector*>::iterator iterTetraT = _tetra.begin(); iterTetraT != _tetra.end(); ++iterTetraT) { SplitterTetra *tmp=*iterTetraT; + tmp->clearVolumesCache(); - double volume = tmp->intersectSourceCell(*iterCellS, baryCentre); +#if 1 //CS_ALM + double volume = tmp->intersectSourceCell(*iterCellS, + baryCentre); +#else + double volume = tmp->intersectSourceCell(*iterCellS,normCellType, cellModelCell, nbOfNodes4Type,cellNodes, isOutside, + isTargetOutside,baryCentre); +#endif if ( volume > 0 ) { interVolume += volume; diff --git a/src/INTERP_KERNEL/SplitterTetra.hxx b/src/INTERP_KERNEL/SplitterTetra.hxx index 5127c6d84..b21dceb9b 100644 --- a/src/INTERP_KERNEL/SplitterTetra.hxx +++ b/src/INTERP_KERNEL/SplitterTetra.hxx @@ -24,6 +24,8 @@ #include "TetraAffineTransform.hxx" #include "InterpolationOptions.hxx" #include "InterpKernelHashMap.hxx" +#include "CellModel.hxx" + #include #include @@ -36,6 +38,7 @@ namespace INTERP_KERNEL * \brief Class representing a triangular face, used as key in caching hash map in SplitterTetra. * */ + #if 1 //CS_ALM class TriangleFaceKey { public: @@ -87,7 +90,7 @@ namespace INTERP_KERNEL /// hash value for the object, calculated in the constructor int _hashVal; }; - + /** * Method to sort three integers in ascending order * @@ -153,8 +156,9 @@ namespace INTERP_KERNEL return key.hashVal(); } }; +#endif } - +//#endif namespace INTERP_KERNEL { @@ -172,7 +176,12 @@ namespace INTERP_KERNEL ~SplitterTetra(); +#if 1 //CS_ALM double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0); +#else + double intersectSourceCell(typename MyMeshType::MyConnType srcCell,NormalizedCellType normCellType, const CellModel& cellModelCell, unsigned nbOfNodes4Type, + int* cellNodes, bool* isOutside, bool isTargetOutside, double* baryCentre=0); +#endif double intersectTetra(const double** tetraCorners); @@ -189,7 +198,11 @@ namespace INTERP_KERNEL inline void createAffineTransform(const double** corners); inline void checkIsOutside(const double* pt, bool* isOutside) const; inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum); +#if 1//CS_ALM inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key); +#else + inline double calculateVolume(TransformedTriangle& tri/*CS_ALM , const TriangleFaceKey& key*/); +#endif /// disallow copying @@ -205,12 +218,14 @@ namespace INTERP_KERNEL /// HashMap relating node numbers to transformed nodes, used for caching HashMap< int , double* > _nodes; +#if 1//CS_ALM /// HashMap relating triangular faces to calculated volume contributions, used for caching HashMap< TriangleFaceKey, double // #ifdef WIN32 // , hash_compare // #endif > _volumes; +#endif /// reference to the source mesh const MyMeshType& _src_mesh; @@ -269,7 +284,11 @@ namespace INTERP_KERNEL inline void SplitterTetra::calculateNode(typename MyMeshType::MyConnType globalNodeNum) { const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum; +#if 1 //CS_ALM double* transformedNode = new double[MyMeshType::MY_SPACEDIM]; +#else + double transformedNode[MyMeshType::MY_SPACEDIM]; +#endif assert(transformedNode != 0); _t->apply(transformedNode, node); _nodes[globalNodeNum] = transformedNode; @@ -282,12 +301,20 @@ namespace INTERP_KERNEL * @param tri triangle for which to calculate the volume contribution * @param key key associated with the face */ +#if 1 //CS_ALM template inline void SplitterTetra::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key) { const double vol = tri.calculateIntersectionVolume(); _volumes.insert(std::make_pair(key, vol)); } +#else + template + inline double SplitterTetra::calculateVolume(TransformedTriangle& tri /*CS_ALM , const TriangleFaceKey& key*/) + { + return tri.calculateIntersectionVolume(); + } +#endif template class SplitterTetra2 @@ -355,7 +382,11 @@ namespace INTERP_KERNEL inline const double* SplitterTetra2::getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node) { // replace "at()" with [] for unsafe but faster access +#if 1 //CS_ALM return _nodes.at(node); +#else + return _nodes[node]; +#endif } /** @@ -368,7 +399,11 @@ namespace INTERP_KERNEL template const double* SplitterTetra2::getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId) { +#if 1 //CS_ALM const double *ret=_nodes.at(node); +#else + const double *ret=_nodes[node]; +#endif if(node<8) nodeId=_node_ids[node]; else diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx index 0372bc4a6..6c5bfa5c0 100644 --- a/src/INTERP_KERNEL/SplitterTetra.txx +++ b/src/INTERP_KERNEL/SplitterTetra.txx @@ -133,7 +133,9 @@ namespace INTERP_KERNEL template void SplitterTetra::clearVolumesCache() { +#if 1 //CS_ALM _volumes.clear(); +#endif } /*! @@ -175,9 +177,26 @@ namespace INTERP_KERNEL * * @param element global number of the source element in C mode. */ + + /*const NormalizedCellType polyType, + const int polyNodesNbr, + const int *const polyNodes, + const double *const *const polyCoords, + std::set& listOfTetraFacesTreated*/ +#if 1 //CS_ALM template double SplitterTetra::intersectSourceCell(typename MyMeshType::MyConnType element, double* baryCentre) +#else + template + double SplitterTetra::intersectSourceCell(typename MyMeshType::MyConnType element,NormalizedCellType normCellType, + const CellModel& cellModelCell, + unsigned nbOfNodes4Type, + int* cellNodes, + bool* isOutside, + bool isTargetOutside, + double* baryCentre) +#endif { typedef typename MyMeshType::MyConnType ConnType; const NumberingPolicy numPol=MyMeshType::My_numPol; @@ -190,6 +209,7 @@ namespace INTERP_KERNEL return 0.0; } + //#if 1 //CS_ALM // get type of cell NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT::indFC(element)); const CellModel& cellModelCell=CellModel::GetCellModel(normCellType); @@ -200,6 +220,7 @@ namespace INTERP_KERNEL // calculate the coordinates of the nodes int *cellNodes=new int[nbOfNodes4Type]; + //#else for(int i = 0;i<(int)nbOfNodes4Type;++i) { // we could store mapping local -> global numbers too, but not sure it is worth it @@ -215,6 +236,7 @@ namespace INTERP_KERNEL checkIsOutside(_nodes[globalNodeNum], isOutside); } + //#endif // halfspace filtering check // NB : might not be beneficial for caching of triangles @@ -223,9 +245,13 @@ namespace INTERP_KERNEL if(isOutside[i]) { isTargetOutside = true; +#if 0 //CS_ALM + break; +#endif } } + double totalVolume = 0.0; if(!isTargetOutside) @@ -264,20 +290,42 @@ namespace INTERP_KERNEL case NORM_TRI3: { // create the face key - TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]); +#if 1 //CS_ALM + TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]); // calculate the triangle if needed + //CS_ALM : test optimizing if(_volumes.find(key) == _volumes.end()) { TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]); - calculateVolume(tri, key); + calculateVolume(tri, key); totalVolume += _volumes[key]; + if ( baryCentre ) baryCalculator.addSide( tri ); - } else { + } else { // count negative as face has reversed orientation totalVolume -= _volumes[key]; - } + } + +#else + //CS_ALMTriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]); + + // calculate the triangle if needed + //CS_ALM : test optimizing + //CS_ALM if(_volumes.find(key) == _volumes.end()) + //CS_ALM { + TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]); + //CS_ALM calculateVolume(tri, key); + //CS_ALM totalVolume += _volumes[key]; + totalVolume += calculateVolume(tri); + if ( baryCentre ) + baryCalculator.addSide( tri ); + //CS_ALM } else { + //CS_ALM // count negative as face has reversed orientation + //CS_ALM totalVolume -= _volumes[key]; + //CS_ALM } +#endif } break; @@ -299,30 +347,62 @@ namespace INTERP_KERNEL // calculate the triangles if needed // local nodes 1, 2, 3 - TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]); - if(_volumes.find(key1) == _volumes.end()) - { +#if 1 // CS_ALM + TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]); + if(_volumes.find(key1) == _volumes.end()) + { TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]); - calculateVolume(tri, key1); - totalVolume += _volumes[key1]; - } else { - // count negative as face has reversed orientation - totalVolume -= _volumes[key1]; + calculateVolume(tri, key1); + totalVolume += _volumes[key1]; + + } else { + // count negative as face has reversed orientation + totalVolume -= _volumes[key1]; } // local nodes 1, 3, 4 - TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]); + TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]); if(_volumes.find(key2) == _volumes.end()) - { + { TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[2]], _nodes[faceNodes[3]]); calculateVolume(tri, key2); - totalVolume += _volumes[key2]; - } + totalVolume += _volumes[key2]; + + } else { - // count negative as face has reversed orientation + // count negative as face has reversed orientation totalVolume -= _volumes[key2]; } +#else + //CS_ALM TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]); + //CS_ALM if(_volumes.find(key1) == _volumes.end()) + { + TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]); + //CS_ALM calculateVolume(tri, key1); + //CS_ALM totalVolume += _volumes[key1]; + totalVolume += calculateVolume(tri); + //CS_ALM } else { + //CS_ALM // count negative as face has reversed orientation + //CS_ALM totalVolume -= _volumes[key1]; + } + + // local nodes 1, 3, 4 + //CS_ALM TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]); + //CS_ALM if(_volumes.find(key2) == _volumes.end()) + { + TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[2]], _nodes[faceNodes[3]]); + //CS_ALM calculateVolume(tri, key2); + //CS_ALM totalVolume += _volumes[key2]; + + totalVolume += calculateVolume(tri); + } + //CS_ALM else + //CS_ALM { + //CS_ALM // count negative as face has reversed orientation + //CS_ALM totalVolume -= _volumes[key2]; + //CS_ALM } +#endif } break; @@ -331,17 +411,35 @@ namespace INTERP_KERNEL int nbTria = nbFaceNodes - 2; // split polygon into nbTria triangles for ( int iTri = 0; iTri < nbTria; ++iTri ) { - TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1+iTri], faceNodes[2+iTri]); +#if 1 //CS_ALM + TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1+iTri], faceNodes[2+iTri]); if(_volumes.find(key) == _volumes.end()) { TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1+iTri]], _nodes[faceNodes[2+iTri]]); calculateVolume(tri, key); totalVolume += _volumes[key]; + + } else { totalVolume -= _volumes[key]; } +#else + //CS_ALMTriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1+iTri], faceNodes[2+iTri]); + //CS_ALM if(_volumes.find(key) == _volumes.end()) + //CS_ALM { + TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1+iTri]], _nodes[faceNodes[2+iTri]]); + //CS_ALMcalculateVolume(tri, key); + //CS_ALMtotalVolume += _volumes[key]; + totalVolume += calculateVolume(tri); + + //CS_ALM } + //CS_ALM else + //CS_ALM { + //CS_ALM totalVolume -= _volumes[key]; + //CS_ALM } +#endif } } break; @@ -358,7 +456,11 @@ namespace INTERP_KERNEL _t->reverseApply( baryCentre, baryCentre ); } } +#if 1 //CS_ALM delete [] cellNodes; +#else + +#endif // reset if it is very small to keep the matrix sparse // is this a good idea? if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT)) diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index 661d71e65..68fd7a08b 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -135,6 +135,10 @@ namespace INTERP_KERNEL preCalculateTriangleSurroundsEdge(); preCalculateTripleProducts(); +#if 1 //CS_ALM + _indiceA = 0; + _indiceB = 0; +#endif } @@ -147,14 +151,29 @@ namespace INTERP_KERNEL TransformedTriangle::~TransformedTriangle() { // 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 } /** @@ -203,7 +222,11 @@ 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); @@ -214,7 +237,11 @@ 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 ... "); @@ -244,28 +271,63 @@ namespace INTERP_KERNEL */ void TransformedTriangle::calculateIntersectionPolygons() { - assert(_polygonA.size() == 0); - assert(_polygonB.size() == 0); +#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]; + + } + _indiceA += 1; + LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A "); +#endif + 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]; + + } + _indiceB += 1; + LOG(3,"Surface-edge : " << vToStr(ptB) << " added to B "); +#endif + } } @@ -276,10 +338,23 @@ 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[_indiceB][i] = ptB[i]; + } + _indiceB += 1; + LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B"); +#endif + + } } @@ -306,16 +381,41 @@ 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); + _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]; + + } + _indiceA += 1; + LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A"); +#endif + + 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]; + + } + _indiceB += 1; + LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B"); +#endif + } } @@ -328,6 +428,7 @@ 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); @@ -338,6 +439,27 @@ namespace INTERP_KERNEL 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]; + + } + _indiceA += 1; + LOG(3,"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]; + } + _indiceB += 1; + } +#endif } } @@ -351,6 +473,7 @@ 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); @@ -362,6 +485,27 @@ namespace INTERP_KERNEL 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[_indiceA][i] = ptA[i]; + } + _indiceA += 1; + LOG(3,"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]; + } + _indiceB += 1; + LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B"); + } +#endif } } @@ -370,10 +514,21 @@ 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]; + } + _indiceB += 1; + LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B"); +#endif } } @@ -392,10 +547,20 @@ 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]; + } + _indiceB += 1; + LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B"); +#endif } } } @@ -407,30 +572,64 @@ 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]; + } + _indiceA += 1; + LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A"); +#endif } // 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]; + } + _indiceB += 1; + LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B"); +#endif } // 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]; + } + _indiceB += 1; + LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B"); +#endif } } @@ -446,15 +645,43 @@ 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]; + if (poly == A) { + for (int i=0; i < _indiceA; i++) { + for (int j = 0; j < 3; j++) { + polygon[i][j] = _polygonA[i][j]; + } + } + } + else { + for (int i=0; i < _indiceB; i++) { + for (int j = 0; j < 3; j++) { + polygon[i][j] = _polygonB[i][j]; + } + } + } + +#endif + // 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) { @@ -465,6 +692,7 @@ namespace INTERP_KERNEL { for(int i = 0 ; i < m ; ++i) { + const double* pt = polygon[i]; for(int j = 0 ; j < 3 ; ++j) { @@ -494,10 +722,32 @@ 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]; + if (poly == A) { + for (int i=0; i < _indiceA; i++) { + for (int j = 0; j < 3; j++) { + polygon[i][j] = _polygonA[i][j]; + } + } + } + else { + for (int i=0; i < _indiceB; i++) { + for (int j = 0; j < 3; j++) { + polygon[i][j] = _polygonB[i][j]; + } + } + } + + int taillePoly = (poly == A) ? _indiceA : _indiceB; + if (taillePoly == 0) + return; +#endif + // determine type of sorting CoordType type = SortOrder::XY; @@ -516,7 +766,7 @@ namespace INTERP_KERNEL { type = SortOrder::YZ; } - } + } // create order object SortOrder order(barycenter, type); @@ -524,10 +774,20 @@ 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 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])); } @@ -550,13 +810,39 @@ 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) { + for (int i=0; i < _indiceA; i++) { + for (int j = 0; j < 3; j++) { + polygon[i][j] = _polygonA[i][j]; + } + } + } + else { + for (int i=0; i < _indiceB; i++) { + for (int j = 0; j < 3; j++) { + polygon[i][j] = _polygonB[i][j]; + } + } + } + + + + int taillePoly = (poly == A) ? _indiceA : _indiceB; + +#endif double vol = 0.0; - const int m = polygon.size(); + 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) diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 84877b9b9..2f9c99650 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -132,7 +132,11 @@ 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(); @@ -143,8 +147,12 @@ 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 double getVolume() const { return _volume; } @@ -181,6 +189,7 @@ namespace INTERP_KERNEL inline bool testSurfaceEdgeIntersection(const TetraEdge edge) const; + void calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const; inline bool testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const; @@ -278,11 +287,27 @@ 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]; + int _indiceA; + int _indiceB; + +#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 /// Array holding the coordinates of the barycenter of the polygon A /// This point is calculated in calculatePolygonBarycenter diff --git a/src/INTERP_KERNEL/TransformedTriangleIntersect.cxx b/src/INTERP_KERNEL/TransformedTriangleIntersect.cxx index faf7c9b61..94dc1e762 100644 --- a/src/INTERP_KERNEL/TransformedTriangleIntersect.cxx +++ b/src/INTERP_KERNEL/TransformedTriangleIntersect.cxx @@ -148,6 +148,7 @@ namespace INTERP_KERNEL * @param edge edge of tetrahedron * @param pt array of three doubles in which to store the coordinates of the intersection point */ + void TransformedTriangle::calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const { assert(edge < H01); diff --git a/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx b/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx index 6aecef649..9f2df0e4f 100644 --- a/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx +++ b/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx @@ -83,7 +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 ) { if ( !epsilonEqual( triNormal[_Z], 0 )) @@ -103,21 +107,42 @@ namespace INTERP_KERNEL } // 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(); -#ifdef DMP_UNITTETRAINTERSECTIONBARY +#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 std::cout.precision(18); std::cout << "**** int polygon() " << std::endl; while ( p != pEnd ) { - double* pp = *p++; + double* pp = p; std::cout << pEnd - p << ": ( " << pp[0] << ", " << pp[1] << ", " << pp[2] << " )" << std::endl; + 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 + if ( p == pEnd ) { #ifdef DMP_UNITTETRAINTERSECTIONBARY @@ -142,8 +167,10 @@ namespace INTERP_KERNEL if (_isTetraInversed) reverse = !reverse; // store polygon + _faces.push_back( std::vector< double* > () ); std::vector< double* >& faceCorner = _faces.back(); + faceCorner.resize( pPolygonA->size()/* + 1*/ ); int i = 0; @@ -151,10 +178,18 @@ namespace INTERP_KERNEL { 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 ( i==0 || !samePoint( *polyF, faceCorner[i-1] )) { +#if 0 //CS_ALM copyVector3( *polyF, faceCorner[i] = new double[3] ); - else +#else + double var[3]= {0.0,0.0,0.0}; + faceCorner[i] = var; + copyVector3( *polyF, faceCorner[i] ); + //faceCorner[i] = var; +#endif + } else { --i; + } polyNormal[0] *= -1.; polyNormal[1] *= -1.; polyNormal[2] *= -1.; @@ -163,10 +198,17 @@ namespace INTERP_KERNEL { 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 ( i==0 || !samePoint( *polyF, faceCorner[i-1] )) { +#if 0 //CS_ALM copyVector3( *polyF, faceCorner[i] = new double[3] ); - else +#else + double var[3]= {0.0,0.0,0.0}; + faceCorner[i] = var; + copyVector3( *polyF, faceCorner[i] ); +#endif + } else { --i; + } } if ( i < 3 ) { @@ -249,6 +291,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; @@ -256,16 +301,30 @@ 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 )); } @@ -315,7 +374,15 @@ 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]; @@ -327,9 +394,15 @@ 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]) / polygon.size(), 1. ))) + ( epsilonEqual( (coordSum[0]+coordSum[1]+coordSum[2]) / polySize, 1. ))) { +#endif sideAdded[3] = ++nbAddedSides != 0 ; + } } if ( nbAddedSides == NB_TETRA_SIDES ) return nbAddedSides; @@ -339,8 +412,11 @@ 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; @@ -354,11 +430,18 @@ namespace INTERP_KERNEL for ( int iF = 0; iF < nbIntersectPolygs; ++f, ++iF ) // loop on added intersection polygons { std::vector< double* >& polygon = *f; - for ( int i = 0; i < (int)polygon.size(); ++i ) +#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 ) { @@ -369,10 +452,20 @@ 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; } } @@ -381,10 +474,19 @@ 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 } } } @@ -400,8 +502,15 @@ 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; } @@ -587,7 +696,12 @@ 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; @@ -712,31 +826,42 @@ 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; + { + + delete[] *it; + *it = 0; } + for(std::vector::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it) { + delete[] *it; + *it = 0; } +#endif _polygonA.clear(); _polygonB.clear(); + 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(); } } -- 2.39.2