From 2094626def1d9fcae0697d1b410b4e75fb47c4c5 Mon Sep 17 00:00:00 2001 From: vbd Date: Thu, 23 Aug 2007 14:20:58 +0000 Subject: [PATCH] staffan : * factored out intersection code into class Intersector3D * added alternative filtering with BBTree to Interpolation3D * replaced std::cout with LOG - macros * added equalEpsilonRelative for better checking of equality of doubles --- src/INTERP_KERNEL/BoundingBox.cxx | 10 +- src/INTERP_KERNEL/BoundingBox.hxx | 2 + src/INTERP_KERNEL/Interpolation3D.cxx | 634 ++++++++++-------- src/INTERP_KERNEL/Interpolation3D.hxx | 2 +- src/INTERP_KERNEL/Intersector3D.cxx | 148 ++++ src/INTERP_KERNEL/Intersector3D.hxx | 41 ++ src/INTERP_KERNEL/Makefile.in | 13 +- src/INTERP_KERNEL/MeshElement.cxx | 129 +--- src/INTERP_KERNEL/MeshElement.hxx | 49 +- src/INTERP_KERNEL/MeshRegion.cxx | 23 +- src/INTERP_KERNEL/MeshRegion.hxx | 8 +- src/INTERP_KERNEL/TetraAffineTransform.hxx | 22 +- src/INTERP_KERNEL/TransformedTriangle.cxx | 75 ++- src/INTERP_KERNEL/TransformedTriangle.hxx | 10 + .../TransformedTriangle_intersect.cxx | 44 +- .../TransformedTriangle_math.cxx | 69 +- 16 files changed, 739 insertions(+), 540 deletions(-) create mode 100644 src/INTERP_KERNEL/Intersector3D.cxx create mode 100644 src/INTERP_KERNEL/Intersector3D.hxx diff --git a/src/INTERP_KERNEL/BoundingBox.cxx b/src/INTERP_KERNEL/BoundingBox.cxx index 3817afb3a..7bf4b12b0 100644 --- a/src/INTERP_KERNEL/BoundingBox.cxx +++ b/src/INTERP_KERNEL/BoundingBox.cxx @@ -145,6 +145,13 @@ namespace INTERP_UTILS } } + void BoundingBox::dumpCoords() const + { + std::cout << "[xmin, xmax] = [" << _coords[XMIN] << ", " << _coords[XMAX] << "]" << " | "; + std::cout << "[ymin, ymax] = [" << _coords[YMIN] << ", " << _coords[YMAX] << "]" << " | "; + std::cout << "[zmin, zmax] = [" << _coords[ZMIN] << ", " << _coords[ZMAX] << "]"; + std::cout << std::endl; + } bool BoundingBox::isValid() const { @@ -153,7 +160,8 @@ namespace INTERP_UTILS { if(_coords[c] >= _coords[c + 3]) { - std::cout << "BoundingBox |: coordinate " << c << " is invalid : " << _coords[c] << " >= " << _coords[c+3] << std::endl; + std::cout << "+++ Error in BoundingBox |: coordinate " << c << " is invalid : " + <<_coords[c] << " >= " << _coords[c+3] << std::endl; valid = false; } } diff --git a/src/INTERP_KERNEL/BoundingBox.hxx b/src/INTERP_KERNEL/BoundingBox.hxx index 49640555e..03bc074b4 100644 --- a/src/INTERP_KERNEL/BoundingBox.hxx +++ b/src/INTERP_KERNEL/BoundingBox.hxx @@ -71,6 +71,8 @@ namespace INTERP_UTILS void updateWithPoint(const double* pt); + void dumpCoords() const; + private: bool isValid() const; diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx index 8d96f5d45..faf19c113 100644 --- a/src/INTERP_KERNEL/Interpolation3D.cxx +++ b/src/INTERP_KERNEL/Interpolation3D.cxx @@ -8,334 +8,408 @@ #include "TetraAffineTransform.hxx" #include "TransformedTriangle.hxx" #include "VectorUtils.hxx" +#include "Intersector3D.hxx" +#include "Log.hxx" + using namespace INTERP_UTILS; using namespace MEDMEM; using namespace MED_EN; +#define USE_RECURSIVE_BBOX_FILTER + +#ifndef USE_RECURSIVE_BBOX_FILTER +#include "BBTree.H" +#endif + namespace MEDMEM { - /** - * Default constructor - * - */ - Interpolation3D::Interpolation3D() - { - // not implemented - } + /** + * Default constructor + * + */ + Interpolation3D::Interpolation3D() + { + // not implemented + } - /** - * Destructor - * - */ - Interpolation3D::~Interpolation3D() - { - // not implemented - } - - /** - * Main interface method of the class, which verifies that the meshes are valid for the calculation, - * and then returns the matrix of intersection volumes. - * - * @param srcMesh 3-dimensional source mesh - * @param targetMesh 3-dimesional target mesh, containing only tetraedra - * @returns vector containing for each element i of the source mesh, a map giving for each element j - * of the target mesh which i intersects, the volume of the intersection - */ - IntersectionMatrix Interpolation3D::interpol_maillages(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh) - { - //} it seems wasteful to make a copy here - IntersectionMatrix matrix; + /** + * Destructor + * + */ + Interpolation3D::~Interpolation3D() + { + // not implemented + } + + /** + * Main interface method of the class, which verifies that the meshes are valid for the calculation, + * and then returns the matrix of intersection volumes. + * + * @param srcMesh 3-dimensional source mesh + * @param targetMesh 3-dimesional target mesh, containing only tetraedra + * @returns vector containing for each element i of the source mesh, a map giving for each element j + * of the target mesh which i intersects, the volume of the intersection + */ + IntersectionMatrix Interpolation3D::interpol_maillages(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh) + { + //} it seems wasteful to make a copy here + IntersectionMatrix matrix; - // we should maybe do more sanity checking here - eliminating meshes that - // are too complicated + // we should maybe do more sanity checking here - eliminating meshes that + // are too complicated - calculateIntersectionVolumes(srcMesh, targetMesh, matrix); - return matrix; - } + calculateIntersectionVolumes(srcMesh, targetMesh, matrix); + return matrix; + } - /** - * Performs a depth-first search over srcMesh, using bounding boxes to recursively eliminate the elements of targetMesh - * which cannot intersect smaller and smaller regions of srcMesh. At each level, each region is divided in two, forming - * a binary search tree with leaves consisting of only one element of the source mesh together with the elements of the - * target mesh that can intersect it. The recursion is implemented with a stack of RegionNodes, each one containing a - * source region and a target region. Each region has an associated bounding box and a vector of pointers to the elements - * that belong to it. Each element contains a bounding box, an element type and an index in the MEDMEM ConnectivityIndex array. - * - * @param srcMesh 3-dimensional source mesh - * @param targetMesh 3-dimesional target mesh, containing only tetraedra - * @param matrix vector of maps in which the result is stored - * - */ - void Interpolation3D::calculateIntersectionVolumes(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, IntersectionMatrix& matrix) - { - // calculate descending connectivities - srcMesh.calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL); - targetMesh.calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL); - - // create MeshElement objects corresponding to each element of the two meshes + /** + * Performs a depth-first search over srcMesh, using bounding boxes to recursively eliminate the elements of targetMesh + * which cannot intersect smaller and smaller regions of srcMesh. At each level, each region is divided in two, forming + * a binary search tree with leaves consisting of only one element of the source mesh together with the elements of the + * target mesh that can intersect it. The recursion is implemented with a stack of RegionNodes, each one containing a + * source region and a target region. Each region has an associated bounding box and a vector of pointers to the elements + * that belong to it. Each element contains a bounding box, an element type and an index in the MEDMEM ConnectivityIndex array. + * + * @param srcMesh 3-dimensional source mesh + * @param targetMesh 3-dimesional target mesh, containing only tetraedra + * @param matrix vector of maps in which the result is stored + * + */ + void Interpolation3D::calculateIntersectionVolumes(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, IntersectionMatrix& matrix) + { + // calculate descending connectivities + // srcMesh.calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL); + // targetMesh.calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL); + + // create MeshElement objects corresponding to each element of the two meshes - const int numSrcElems = srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS); - const int numTargetElems = targetMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS); + const int numSrcElems = srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS); + const int numTargetElems = targetMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS); - // std::cout << "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements " << std::endl; + Intersector3D* intersector = new Intersector3D(srcMesh, targetMesh); - // create empty maps for all source elements - matrix.resize(numSrcElems); + LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements "); - MeshElement* srcElems[numSrcElems]; - MeshElement* targetElems[numTargetElems]; + // create empty maps for all source elements + matrix.resize(numSrcElems); + + + MeshElement* srcElems[numSrcElems]; + MeshElement* targetElems[numTargetElems]; - std::map indices; + std::map indices; - for(int i = 0 ; i < numSrcElems ; ++i) - { - //const int index = srcMesh.getConnectivityIndex(MED_NODAL, MED_CELL)[i]; - const medGeometryElement type = srcMesh.getElementType(MED_CELL, i + 1); - srcElems[i] = new MeshElement(i + 1, type, srcMesh); + for(int i = 0 ; i < numSrcElems ; ++i) + { + //const int index = srcMesh.getConnectivityIndex(MED_NODAL, MED_CELL)[i]; + const medGeometryElement type = srcMesh.getElementType(MED_CELL, i + 1); + srcElems[i] = new MeshElement(i + 1, type, srcMesh); - } - - for(int i = 0 ; i < numTargetElems ; ++i) - { - // const int index = targetMesh.getConnectivityIndex(MED_NODAL, MED_CELL)[i]; - const medGeometryElement type = targetMesh.getElementType(MED_CELL, i + 1); - targetElems[i] = new MeshElement(i + 1, type, targetMesh); - } + } + + for(int i = 0 ; i < numTargetElems ; ++i) + { + // const int index = targetMesh.getConnectivityIndex(MED_NODAL, MED_CELL)[i]; + const medGeometryElement type = targetMesh.getElementType(MED_CELL, i + 1); + targetElems[i] = new MeshElement(i + 1, type, targetMesh); + } + +#ifdef USE_RECURSIVE_BBOX_FILTER - // create initial RegionNode and fill up its source region with all the source mesh elements and - // its target region with all the target mesh elements whose bbox intersects that of the source region + // create initial RegionNode and fill up its source region with all the source mesh elements and + // its target region with all the target mesh elements whose bbox intersects that of the source region - RegionNode* firstNode = new RegionNode(); + RegionNode* firstNode = new RegionNode(); - MeshRegion& srcRegion = firstNode->getSrcRegion(); - - for(int i = 0 ; i < numSrcElems ; ++i) - { - srcRegion.addElement(srcElems[i]); - } - - MeshRegion& targetRegion = firstNode->getTargetRegion(); - - for(int i = 0 ; i < numTargetElems ; ++i) - { - if(!srcRegion.isDisjointWithElementBoundingBox( *(targetElems[i]) )) - { - targetRegion.addElement(targetElems[i]); - } - } - - // using a stack, descend recursively, creating at each step two new RegionNodes having as source region the left and - // right part of the source region of the current node (determined using MeshRegion::split()) and as target region all the - // elements of the target mesh whose bounding box intersects the corresponding part - // Continue until the source region contains only one element, at which point the intersection volumes are - // calculated with all the remaining target mesh elements and stored in the matrix if they are non-zero. - - stack nodes; - nodes.push(firstNode); - - while(!nodes.empty()) - { - RegionNode* currNode = nodes.top(); - nodes.pop(); - // std::cout << "Popping node " << std::endl; - - if(currNode->getSrcRegion().getNumberOfElements() == 1) - { - // std::cout << " - One element" << std::endl; - - // volume calculation - MeshElement* srcElement = *(currNode->getSrcRegion().getBeginElements()); + MeshRegion& srcRegion = firstNode->getSrcRegion(); + + for(int i = 0 ; i < numSrcElems ; ++i) + { + srcRegion.addElement(srcElems[i], srcMesh); + } + + MeshRegion& targetRegion = firstNode->getTargetRegion(); + + for(int i = 0 ; i < numTargetElems ; ++i) + { + if(!srcRegion.isDisjointWithElementBoundingBox( *(targetElems[i]) )) + { + targetRegion.addElement(targetElems[i], targetMesh); + } + } + + // using a stack, descend recursively, creating at each step two new RegionNodes having as source region the left and + // right part of the source region of the current node (determined using MeshRegion::split()) and as target region all the + // elements of the target mesh whose bounding box intersects the corresponding part + // Continue until the source region contains only one element, at which point the intersection volumes are + // calculated with all the remaining target mesh elements and stored in the matrix if they are non-zero. + + stack nodes; + nodes.push(firstNode); + + while(!nodes.empty()) + { + RegionNode* currNode = nodes.top(); + nodes.pop(); + LOG(4, "Popping node "); + + if(currNode->getSrcRegion().getNumberOfElements() == 1) + { + LOG(4, " - One element"); + + // volume calculation + MeshElement* srcElement = *(currNode->getSrcRegion().getBeginElements()); - // NB : srcElement indices are from 0 .. numSrcElements - 1 - // targetElement indicies from 1 .. numTargetElements - // maybe this is not ideal ... - const int srcIdx = srcElement->getIndex() - 1; - std::map< int, double >* volumes = &(matrix[srcIdx]); - - for(vector::const_iterator iter = currNode->getTargetRegion().getBeginElements() ; - iter != currNode->getTargetRegion().getEndElements() ; ++iter) - { - const double vol = calculateIntersectionVolume(*srcElement, **iter); - if(!epsilonEqual(vol, 0.0, 1.0e-10)) + // NB : srcElement indices are from 0 .. numSrcElements - 1 + // targetElement indicies from 1 .. numTargetElements + // maybe this is not ideal ... + const int srcIdx = srcElement->getIndex(); + std::map< int, double >* volumes = &(matrix[srcIdx - 1]); + + for(vector::const_iterator iter = currNode->getTargetRegion().getBeginElements() ; + iter != currNode->getTargetRegion().getEndElements() ; ++iter) + { + const int targetIdx = (*iter)->getIndex(); + const double vol = intersector->intersectCells(srcIdx, targetIdx); + if(!epsilonEqual(vol, 0.0, 1.0e-10)) // if(vol != 0.0) - { - const int targetIdx = (*iter)->getIndex(); - - volumes->insert(make_pair(targetIdx, vol)); - // std::cout << "Result : V (" << srcIdx << ", " << targetIdx << ") = " << matrix[srcIdx][targetIdx] << std::endl; - } - } - } - else // recursion - { - - // std::cout << " - Recursion" << std::endl; - - RegionNode* leftNode = new RegionNode(); - RegionNode* rightNode = new RegionNode(); + { + volumes->insert(make_pair(targetIdx, vol)); + LOG(2, "Result : V (" << srcIdx << ", " << targetIdx << ") = " << matrix[srcIdx - 1][targetIdx]); + } + } + } + else // recursion + { + + LOG(4, " - Recursion"); + + RegionNode* leftNode = new RegionNode(); + RegionNode* rightNode = new RegionNode(); - // split current source region - //} decide on axis - static BoundingBox::BoxCoord axis = BoundingBox::XMAX; + // split current source region + //} decide on axis + static BoundingBox::BoxCoord axis = BoundingBox::XMAX; - currNode->getSrcRegion().split(leftNode->getSrcRegion(), rightNode->getSrcRegion(), axis); - - // std::cout << "After split, left src region has " << leftNode->getSrcRegion().getNumberOfElements() << - // " elements and right src region has " << rightNode->getSrcRegion().getNumberOfElements() << " elements" << std::endl; - - // ugly hack to avoid problem with enum which does not start at 0 - // I guess I ought to implement ++ for it instead ... - // Anyway, it basically chooses the next axis, circually - axis = (axis != BoundingBox::ZMAX) ? static_cast(axis + 1) : BoundingBox::XMAX; - - // add target elements of curr node that overlap the two new nodes - // std::cout << " -- Adding target elements" << std::endl; - int numLeftElements = 0; - int numRightElements = 0; - for(vector::const_iterator iter = currNode->getTargetRegion().getBeginElements() ; - iter != currNode->getTargetRegion().getEndElements() ; ++iter) - { - //// std::cout << " --- New target node" << std::endl; + currNode->getSrcRegion().split(leftNode->getSrcRegion(), rightNode->getSrcRegion(), axis, srcMesh); + + LOG(5, "After split, left src region has " << leftNode->getSrcRegion().getNumberOfElements() + << " elements and right src region has " << rightNode->getSrcRegion().getNumberOfElements() + << " elements"); + + // ugly hack to avoid problem with enum which does not start at 0 + // I guess I ought to implement ++ for it instead ... + // Anyway, it basically chooses the next axis, circually + axis = (axis != BoundingBox::ZMAX) ? static_cast(axis + 1) : BoundingBox::XMAX; + + // add target elements of curr node that overlap the two new nodes + LOG(5, " -- Adding target elements"); + int numLeftElements = 0; + int numRightElements = 0; + for(vector::const_iterator iter = currNode->getTargetRegion().getBeginElements() ; + iter != currNode->getTargetRegion().getEndElements() ; ++iter) + { + LOG(6, " --- New target node"); - if(!leftNode->getSrcRegion().isDisjointWithElementBoundingBox(**iter)) - { - leftNode->getTargetRegion().addElement(*iter); - ++numLeftElements; - } + if(!leftNode->getSrcRegion().isDisjointWithElementBoundingBox(**iter)) + { + leftNode->getTargetRegion().addElement(*iter, targetMesh); + ++numLeftElements; + } - if(!rightNode->getSrcRegion().isDisjointWithElementBoundingBox(**iter)) - { - rightNode->getTargetRegion().addElement(*iter); - ++numRightElements; - } + if(!rightNode->getSrcRegion().isDisjointWithElementBoundingBox(**iter)) + { + rightNode->getTargetRegion().addElement(*iter, targetMesh); + ++numRightElements; + } - } - - // std::cout << "Left target region has " << numLeftElements << " elements and right target region has " << numRightElements << " elements" << std::endl; - - if(numLeftElements != 0) - { - nodes.push(leftNode); - } - else - { - delete leftNode; - } - - if(numRightElements != 0) - { - nodes.push(rightNode); - } - else - { - delete rightNode; - } - } + } + + LOG(5, "Left target region has " << numLeftElements << " elements and right target region has " + << numRightElements << " elements"); + + if(numLeftElements != 0) + { + nodes.push(leftNode); + } + else + { + delete leftNode; + } + + if(numRightElements != 0) + { + nodes.push(rightNode); + } + else + { + delete rightNode; + } + } - delete currNode; - // std::cout << "Next iteration. Nodes left : " << nodes.size() << std::endl; - } + delete currNode; + LOG(4, "Next iteration. Nodes left : " << nodes.size()); + } + + +#else // use BBTree + // create BBTree structure + // - get bounding boxes + double bboxes[6 * numSrcElems]; + int srcElemIdx[numSrcElems]; + for(int i = 0; i < numSrcElems ; ++i) + { + // get source bboxes in right order + const BoundingBox* box = srcElems[i]->getBoundingBox(); + bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN); + bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX); + bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN); + bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX); + bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN); + bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX); + srcElemIdx[i] = srcElems[i]->getIndex() - 1; + } + + BBTree<3> tree(bboxes, srcElemIdx, 0, numSrcElems); - // free allocated memory - for(int i = 0 ; i < numSrcElems ; ++i) - { - delete srcElems[i]; - } - for(int i = 0 ; i < numTargetElems ; ++i) - { - delete targetElems[i]; - } - - } - - /** - * Calculates volume of intersection between an element of the source mesh and an element of the target mesh. - * The calculation passes by the following steps : - * a) test if srcElement is in the interior of targetElement -> if yes, return volume of srcElement - * --> test by call to Element::isElementIncludedIn - * b) test if targetElement is in the interior of srcElement -> if yes return volume of targetElement - * --> test by call to Element::isElementIncludedIn - * c) test if srcElement and targetElement are disjoint -> if yes return 0.0 [this test is only possible if srcElement is convex] - * --> test by call to Element::isElementTriviallyDisjointWith - * d) (1) find transformation M that takes the target element (a tetraeder) to the unit tetraeder - * --> call calculateTransform() - * (2) divide srcElement in triangles - * --> call triangulateElement() - * (3) for each triangle in element, - * (i) find polygones --> calculateIntersectionPolygones() - * (ii) calculate volume under polygones --> calculateVolumeUnderPolygone() - * (iii) add volume to sum - * (4) return det(M)*sumVolume (det(M) = 6*vol(targetElement)) - * - * @param srcElement - * @param targetElement - * @returns volume of intersection between srcElement and targetElement - * - */ - double Interpolation3D::calculateIntersectionVolume(const MeshElement& srcElement, const MeshElement& targetElement) + // for each target element, get source elements with which to calculate intersection + // - calculate intersection by calling intersectCells + for(int i = 0; i < numTargetElems; ++i) { + const BoundingBox* box = targetElems[i]->getBoundingBox(); + const int targetIdx = targetElems[i]->getIndex(); + + // get target bbox in right order + double targetBox[6]; + targetBox[0] = box->getCoordinate(BoundingBox::XMIN); + targetBox[1] = box->getCoordinate(BoundingBox::XMAX); + targetBox[2] = box->getCoordinate(BoundingBox::YMIN); + targetBox[3] = box->getCoordinate(BoundingBox::YMAX); + targetBox[4] = box->getCoordinate(BoundingBox::ZMIN); + targetBox[5] = box->getCoordinate(BoundingBox::ZMAX); - //std::cout << "Intersecting elems " << srcElement.getIndex() << " and " << targetElement.getIndex() << std::endl; - // (a), (b) and (c) not yet implemented + vector intersectElems; - // (d) : without fine-level filtering (a) - (c) for the time being - // std::cout << "------------------" << std::endl; - // std::cout << "Source : "; - srcElement.dumpCoords(); - // std::cout << "Target : "; - targetElement.dumpCoords(); + tree.getIntersectingElems(targetBox, intersectElems); - // get array of points of target tetraeder - const double* tetraCorners[4]; - for(int i = 0 ; i < 4 ; ++i) + for(vector::const_iterator iter = intersectElems.begin() ; iter != intersectElems.end() ; ++iter) { - tetraCorners[i] = targetElement.getCoordsOfNode(i + 1); + const int srcIdx = *iter + 1; + const double vol = intersector->intersectCells(srcIdx, targetIdx); + + if(!epsilonEqual(vol, 0.0, 1.0e-10)) + { + matrix[srcIdx - 1].insert(make_pair(targetIdx, vol)); + } } + } + +#endif + + // free allocated memory + for(int i = 0 ; i < numSrcElems ; ++i) + { + delete srcElems[i]; + } + for(int i = 0 ; i < numTargetElems ; ++i) + { + delete targetElems[i]; + } - // create AffineTransform - TetraAffineTransform T( tetraCorners ); + delete intersector; + + } + + /** + * Calculates volume of intersection between an element of the source mesh and an element of the target mesh. + * The calculation passes by the following steps : + * a) test if srcElement is in the interior of targetElement -> if yes, return volume of srcElement + * --> test by call to Element::isElementIncludedIn + * b) test if targetElement is in the interior of srcElement -> if yes return volume of targetElement + * --> test by call to Element::isElementIncludedIn + * c) test if srcElement and targetElement are disjoint -> if yes return 0.0 [this test is only possible if srcElement is convex] + * --> test by call to Element::isElementTriviallyDisjointWith + * d) (1) find transformation M that takes the target element (a tetraeder) to the unit tetraeder + * --> call calculateTransform() + * (2) divide srcElement in triangles + * --> call triangulateElement() + * (3) for each triangle in element, + * (i) find polygones --> calculateIntersectionPolygones() + * (ii) calculate volume under polygones --> calculateVolumeUnderPolygone() + * (iii) add volume to sum + * (4) return det(M)*sumVolume (det(M) = 6*vol(targetElement)) + * + * @param srcElement + * @param targetElement + * @returns volume of intersection between srcElement and targetElement + * + */ +#if 0 + double Interpolation3D::calculateIntersectionVolume(const MeshElement& srcElement, const MeshElement& targetElement) + { + + //std::cout << "Intersecting elems " << srcElement.getIndex() << " and " << targetElement.getIndex() << std::endl; + // (a), (b) and (c) not yet implemented + + // (d) : without fine-level filtering (a) - (c) for the time being + // std::cout << "------------------" << std::endl; + std::cout << "Source : "; + srcElement.dumpCoords(); + std::cout << "Target : "; + targetElement.dumpCoords(); + + // get array of points of target tetraeder + const double* tetraCorners[4]; + for(int i = 0 ; i < 4 ; ++i) + { + tetraCorners[i] = targetElement.getCoordsOfNode(i + 1); + } + + // create AffineTransform + TetraAffineTransform T( tetraCorners ); - // check if we have planar tetra element - if(epsilonEqual(T.determinant(), 0.0, 1.0e-16)) - { - // tetra is planar - // std::cout << "Planar tetra -- volume 0" << std::endl; - return 0.0; - } + // check if we have planar tetra element + if(epsilonEqual(T.determinant(), 0.0, 1.0e-16)) + { + // tetra is planar + // std::cout << "Planar tetra -- volume 0" << std::endl; + return 0.0; + } - // std::cout << "Transform : " << std::endl; - // T.dump(); - // std::cout << std::endl; + // std::cout << "Transform : " << std::endl; + // T.dump(); + // std::cout << std::endl; - // triangulate source element faces (assumed tetraeder for the time being) - // do nothing - vector triangles; - srcElement.triangulate(triangles, T); + // triangulate source element faces (assumed tetraeder for the time being) + // do nothing + vector triangles; + srcElement.triangulate(triangles, T); - double volume = 0.0; + double volume = 0.0; - // std::cout << "num triangles = " << triangles.size() << std::endl; - int i = 0; - for(vector::iterator iter = triangles.begin() ; iter != triangles.end(); ++iter) - { - // std::cout << std::endl << "= > Triangle " << ++i << std::endl; - iter->dumpCoords(); - volume += iter->calculateIntersectionVolume(); - } + // std::cout << "num triangles = " << triangles.size() << std::endl; + int i = 0; + for(vector::iterator iter = triangles.begin() ; iter != triangles.end(); ++iter) + { + std::cout << std::endl << "= > Triangle " << ++i << std::endl; + iter->dumpCoords(); + volume += iter->calculateIntersectionVolume(); + } - // std::cout << "Volume = " << volume << ", det= " << T.determinant() << std::endl; + std::cout << "Volume = " << volume << ", det= " << T.determinant() << std::endl; - //? trying without abs to see if the sign of the determinant will always cancel that of the volume - //? but maybe we should take abs( det ( T ) ) or abs ( 1 / det * vol ) + //? trying without abs to see if the sign of the determinant will always cancel that of the volume + //? but maybe we should take abs( det ( T ) ) or abs ( 1 / det * vol ) - //? fault in article, Grandy, [8] : it is the determinant of the inverse transformation that - // should be used + //? fault in article, Grandy, [8] : it is the determinant of the inverse transformation that + // should be used - return std::abs(1.0 / T.determinant() * volume) ; + return std::abs(1.0 / T.determinant() * volume) ; - } + } +#endif } diff --git a/src/INTERP_KERNEL/Interpolation3D.hxx b/src/INTERP_KERNEL/Interpolation3D.hxx index eaf5ece8b..3d2d21705 100644 --- a/src/INTERP_KERNEL/Interpolation3D.hxx +++ b/src/INTERP_KERNEL/Interpolation3D.hxx @@ -101,7 +101,7 @@ namespace MEDMEM * @returns volume of intersection between srcElement and targetElement * */ - double calculateIntersectionVolume(const INTERP_UTILS::MeshElement& srcElement, const INTERP_UTILS::MeshElement& targetElement); + // double calculateIntersectionVolume(const INTERP_UTILS::MeshElement& srcElement, const INTERP_UTILS::MeshElement& targetElement); }; diff --git a/src/INTERP_KERNEL/Intersector3D.cxx b/src/INTERP_KERNEL/Intersector3D.cxx new file mode 100644 index 000000000..37fbe0aa0 --- /dev/null +++ b/src/INTERP_KERNEL/Intersector3D.cxx @@ -0,0 +1,148 @@ +#include "Intersector3D.hxx" + + +#include "Log.hxx" + +namespace MEDMEM +{ + + Intersector3D::Intersector3D(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh) + : _srcMesh(srcMesh), _targetMesh(targetMesh) + { + } + + Intersector3D::~Intersector3D() + { + } + + double Intersector3D::intersectCells(int srcCell, int targetCell) + { + + const medGeometryElement srcType = _srcMesh.getElementType(MED_CELL, srcCell); + // const medGeometryElement targetType = _targetMesh.getElementType(MED_CELL, targetCell); + + // get array of points of target tetraeder + const double* tetraCorners[4]; + for(int i = 0 ; i < 4 ; ++i) + { + tetraCorners[i] = getCoordsOfNode(i + 1, targetCell, _targetMesh); + } + + // create AffineTransform + TetraAffineTransform T( tetraCorners ); + + // check if we have planar tetra element + if(T.determinant() == 0.0) + { + // tetra is planar + LOG(2, "Planar tetra -- volume 0"); + return 0.0; + } + // triangulate source element faces (assumed tetraeder for the time being) + // do nothing + vector triangles; + triangulate(srcType, srcCell, _srcMesh, triangles, T); + + double volume = 0.0; + + LOG(4, "num triangles = " << triangles.size()); + int i = 0; + for(vector::iterator iter = triangles.begin() ; iter != triangles.end(); ++iter) + { + LOG(2, "= > Triangle " << ++i); +#ifdef LOG_ACTIVE + if(LOG_LEVEL >= 2) + { + iter->dumpCoords(); + } +#endif + volume += iter->calculateIntersectionVolume(); + } + + LOG(2, "Volume = " << volume << ", det= " << T.determinant()); + + //? trying without abs to see if the sign of the determinant will always cancel that of the volume + //? but maybe we should take abs( det ( T ) ) or abs ( 1 / det * vol ) + + //? fault in article, Grandy, [8] : it is the determinant of the inverse transformation that + // should be used + return std::abs(1.0 / T.determinant() * volume) ; + + } + + /** + * Triangulate the faces of an element and apply an affine Transform to the triangles + * + * @param triangles vector in which triangles are stored + * @param T affine transform that is applied to the nodes of the triangles + */ + void Intersector3D::triangulate(const medGeometryElement type, const int element, const MEDMEM::MESH& mesh, std::vector& triangles, const TetraAffineTransform& T) const + { + // get cell model for the element + CELLMODEL cellModel(type); + + assert(cellModel.getDimension() == 3); + + // start index in connectivity array for cell + // const int cellIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index] - 1; + + // assert(cellIdx >= 0); + //assert(cellIdx < _mesh->getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS)); + + // loop over faces + for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i) + { + medGeometryElement faceType = cellModel.getConstituentType(1, i); + CELLMODEL faceModel(faceType); + + assert(faceModel.getDimension() == 2); + assert(faceModel.getNumberOfNodes() == 3); + + double transformedNodes[3 * faceModel.getNumberOfNodes()]; + + + // loop over nodes of face + for(int j = 1; j <= faceModel.getNumberOfNodes(); ++j) + { + // offset of node from cellIdx + int localNodeNumber = cellModel.getNodeConstituent(1, i, j); + + assert(localNodeNumber >= 1); + assert(localNodeNumber <= cellModel.getNumberOfNodes()); + + const double* node = getCoordsOfNode(localNodeNumber, element, mesh); + + // transform + //{ not totally efficient since we transform each node once per face + T.apply(&transformedNodes[3*(j-1)], node); + + LOG(4, "Node " << localNodeNumber << " = " << vToStr(node) << " transformed to " + << vToStr(&transformedNodes[3*(j-1)])); + + } + + // to be removed + assert(faceType == MED_TRIA3); + + // create transformed triangles from face + switch(faceType) + { + case MED_TRIA3: + triangles.push_back(TransformedTriangle(&transformedNodes[0], &transformedNodes[3], &transformedNodes[6])); + break; + + // add other cases here to treat hexahedra, pyramides, etc + + default: + std::cout << "+++ Error : Only elements with triangular faces are supported at the moment." << std::endl; + ; + } + } + + + } + + + + +}; diff --git a/src/INTERP_KERNEL/Intersector3D.hxx b/src/INTERP_KERNEL/Intersector3D.hxx new file mode 100644 index 000000000..2d1375bcd --- /dev/null +++ b/src/INTERP_KERNEL/Intersector3D.hxx @@ -0,0 +1,41 @@ +#ifndef __INTERSECTOR_3D_HXX__ +#define __INTERSECTOR_3D_HXX__ + +#include "MEDMEM_define.hxx" +#include "MEDMEM_Mesh.hxx" + +#include "TetraAffineTransform.hxx" +#include "TransformedTriangle.hxx" +#include "MeshUtils.hxx" +#include "Intersector.hxx" + +#include + +using namespace MEDMEM; +using namespace MED_EN; +using namespace INTERP_UTILS; + +namespace MEDMEM +{ + + class Intersector3D : public Intersector + { + + public : + Intersector3D(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh); + + ~Intersector3D(); + + double intersectCells(int srcCell, int targetCell); + + void triangulate(const medGeometryElement type, const int element, const MEDMEM::MESH& mesh, std::vector& triangles, const TetraAffineTransform& T) const; + + private: + const MEDMEM::MESH& _srcMesh; + const MEDMEM::MESH& _targetMesh; + + }; +}; + + +#endif diff --git a/src/INTERP_KERNEL/Makefile.in b/src/INTERP_KERNEL/Makefile.in index 394aa7cc2..1ab573047 100644 --- a/src/INTERP_KERNEL/Makefile.in +++ b/src/INTERP_KERNEL/Makefile.in @@ -50,7 +50,10 @@ TranslationRotationMatrix.hxx\ Interpolation2D.hxx\ Interpolation3DSurf.hxx\ InterpolationUtils.hxx\ -BBTree.H +BBTree.H\ +MeshUtils.hxx\ +Intersector3D.hxx\ +Log.hxx # Libraries targets @@ -65,7 +68,9 @@ MeshRegion.cxx\ RegionNode.cxx\ BoundingBox.cxx\ TetraAffineTransform.cxx\ -Interpolation3D.cxx +Interpolation3D.cxx\ +Intersector3D.cxx\ +#Log.cxx #Interpolation3DSurf.cxx\ #Interpolation2D.cxx\ @@ -83,12 +88,14 @@ LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome CXXFLAGS+=@CXXTMPDPTHFLAGS@ $(MED2_INCLUDES) CPPFLAGS+=$(BOOST_CPPFLAGS) - +CXXFLAGS+=-fprofile-arcs -ftest-coverage +CPPFLAGS+=-fprofile-arcs -ftest-coverage #LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 $(STDLIB) -lmedmem -lutil #LDFLAGS+= $(HDF5_LIBS) $(STDLIB) -lutil +LDFLAGS+= -lgcov #LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. diff --git a/src/INTERP_KERNEL/MeshElement.cxx b/src/INTERP_KERNEL/MeshElement.cxx index 32077a5c7..d13585ad3 100644 --- a/src/INTERP_KERNEL/MeshElement.cxx +++ b/src/INTERP_KERNEL/MeshElement.cxx @@ -3,6 +3,10 @@ #include "TetraAffineTransform.hxx" #include "TransformedTriangle.hxx" #include "MEDMEM_CellModel.hxx" +#include "MeshUtils.hxx" +#include "BoundingBox.hxx" + + namespace INTERP_UTILS { @@ -15,16 +19,16 @@ namespace INTERP_UTILS * @param index global number of element in the mesh */ MeshElement::MeshElement(const int index, const MED_EN::medGeometryElement type, const MEDMEM::MESH& mesh) - : _index(index - 1), _mesh(&mesh), _type(type), _box(0) + : _index(index), _box(0), _type(type) { // get coordinates of vertices - const int numNodes = getNumberNodes(); + const int numNodes = getNumberOfNodesForType(type); const double* vertices[numNodes]; for(int i = 0 ; i < numNodes ; ++i) { - vertices[i] = getCoordsOfNode(i + 1); + vertices[i] = getCoordsOfNode(i + 1, index, mesh); } // create bounding box @@ -76,126 +80,18 @@ namespace INTERP_UTILS return false; } - /** - * Returns the number of nodes of this element - * - * @returns the number of nodes of this element - */ - int MeshElement::getNumberNodes() const - { - assert(_type > 300); - assert(_type < 400); - - // int(type) = yxx where y is dimension of entity (here 3) - // and xx is the number of nodes of the element - return static_cast(_type) - 300; - } - - /** - * Returns the coordinates of a node of this element - * (1 <= node <= #nodes) - * - * @param node the node for which the coordinates are sought - * @returns pointer to an array of 3 doubles containing the coordinates - */ - const double* MeshElement::getCoordsOfNode(int node) const - { - assert(node >= 1); - assert(node <= getNumberNodes()); - const int nodeOffset = node - 1; - const int elemIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index] - 1; - const int connIdx = _mesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL, MED_CELL, MED_ALL_ELEMENTS)[elemIdx + nodeOffset] - 1; - return &(_mesh->getCoordinates(MED_FULL_INTERLACE)[3*connIdx]); - } - /** - * Triangulate the faces of this element and apply an affine Transform to the triangles - * - * @param triangles vector in which triangles are stored - * @param T affine transform that is applied to the nodes of the triangles - */ - void MeshElement::triangulate(std::vector& triangles, const TetraAffineTransform& T) const + int MeshElement::getIndex() const { - // get cell model for the element - CELLMODEL cellModel(_type); - - assert(cellModel.getDimension() == 3); - - // start index in connectivity array for cell - // const int cellIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index] - 1; - - // assert(cellIdx >= 0); - //assert(cellIdx < _mesh->getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS)); - - // loop over faces - for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i) - { - medGeometryElement faceType = cellModel.getConstituentType(1, i); - CELLMODEL faceModel(faceType); - - assert(faceModel.getDimension() == 2); - assert(faceModel.getNumberOfNodes() == 3); - - double transformedNodes[3 * faceModel.getNumberOfNodes()]; - - - // loop over nodes of face - for(int j = 1; j <= faceModel.getNumberOfNodes(); ++j) - { - // offset of node from cellIdx - int localNodeNumber = cellModel.getNodeConstituent(1, i, j); - - assert(localNodeNumber >= 1); - assert(localNodeNumber <= cellModel.getNumberOfNodes()); - - const double* node = getCoordsOfNode(localNodeNumber); - - // transform - //{ not totally efficient since we transform each node once per face - T.apply(&transformedNodes[3*(j-1)], node); - - // std::cout << "Node " << localNodeNumber << " = " << vToStr(node) << " transformed to " << vToStr(&transformedNodes[3*(j-1)]) << std::endl; - - } - - // to be removed - assert(faceType == MED_TRIA3); - - // create transformed triangles from face - - switch(faceType) - { - case MED_TRIA3: - // std::cout << std::endl<< "** Adding triangle " << i << std::endl; - triangles.push_back(TransformedTriangle(&transformedNodes[0], &transformedNodes[3], &transformedNodes[6])); - break; - - // add other cases here to treat hexahedra, pyramides, etc - - default: - // std::cout << "Only elements with triangular faces are supported at the moment." << std::endl; - ; - } - } - - + return _index; } - int MeshElement::getIndex() const + void MeshElement::dumpCoords() const { - return _index + 1; + std::cout << "Bounding box of element " << _index << " is " << std::endl; + _box->dumpCoords(); } - void MeshElement::dumpCoords() const - { - // std::cout << "Element " << _index + 1 << " has nodes " << std::endl; - for(int i = 1 ; i <= getNumberNodes() ; ++i) - { - // std::cout << vToStr(getCoordsOfNode(i)) << ", "; - } - // std::cout << std::endl; - } - /// ElementBBoxOrder bool ElementBBoxOrder::operator()( MeshElement* elem1, MeshElement* elem2) @@ -211,5 +107,4 @@ namespace INTERP_UTILS return coord1 < coord2; } - }; diff --git a/src/INTERP_KERNEL/MeshElement.hxx b/src/INTERP_KERNEL/MeshElement.hxx index f58c36d0f..b08e9d74e 100644 --- a/src/INTERP_KERNEL/MeshElement.hxx +++ b/src/INTERP_KERNEL/MeshElement.hxx @@ -6,17 +6,17 @@ #include "BoundingBox.hxx" -using namespace MEDMEM; -using namespace MED_EN; +using namespace MEDMEM; +using namespace MED_EN; namespace INTERP_UTILS { - class TransformedTriangle; - class TetraAffineTransform; + // class TransformedTriangle; + // class TetraAffineTransform; /** * Class representing a single element of a mesh together with its bounding box. @@ -71,40 +71,32 @@ namespace INTERP_UTILS */ bool isElementTriviallyDisjointWith(const MeshElement& otherElement) const; - /** - * Returns the number of nodes of this element - * - * @returns the number of nodes of this element - */ - int getNumberNodes() const; - - /** - * Returns the coordinates of a node of this element - * (1 <= node <= #nodes) - * - * @param node the node for which the coordinates are sought - * @returns pointer to an array of 3 doubles containing the coordinates - */ - const double* getCoordsOfNode(int node) const; - /** * Triangulate the faces of this element and apply an affine Transform to the triangles * * @param triangles vector in which triangles are stored * @param T affine transform that is applied to the nodes of the triangles */ - void triangulate(std::vector& triangles, const TetraAffineTransform& T) const; + // void triangulate(std::vector& triangles, const TetraAffineTransform& T) const; int getIndex() const; void dumpCoords() const; + + const BoundingBox* getBoundingBox() const + { + return _box; + } + + MED_EN::medGeometryElement getType() const + { + return _type; + } private: const int _index; - const MEDMEM::MESH* _mesh; - const MED_EN::medGeometryElement _type; BoundingBox* _box; // should not change after initialisation - + MED_EN::medGeometryElement _type; }; @@ -112,17 +104,20 @@ namespace INTERP_UTILS class ElementBBoxOrder { public : - + ElementBBoxOrder(BoundingBox::BoxCoord coord) : _coord(coord) { } - + bool operator()(MeshElement* elem1, MeshElement* elem2); - + private : BoundingBox::BoxCoord _coord; }; }; + + + #endif diff --git a/src/INTERP_KERNEL/MeshRegion.cxx b/src/INTERP_KERNEL/MeshRegion.cxx index 128b046c4..1dacdf2b4 100644 --- a/src/INTERP_KERNEL/MeshRegion.cxx +++ b/src/INTERP_KERNEL/MeshRegion.cxx @@ -2,7 +2,7 @@ #include "MeshElement.hxx" -//#include +#include "MeshUtils.hxx" namespace INTERP_UTILS { @@ -34,29 +34,30 @@ namespace INTERP_UTILS * @param element pointer to element to add to region * */ - void MeshRegion::addElement(MeshElement* const element) + void MeshRegion::addElement(MeshElement* const element, const MEDMEM::MESH& mesh) { _elements.push_back(element); + + const int numNodes = getNumberOfNodesForType(element->getType()); + const int elemIdx = element->getIndex(); if(_box == 0) { - const int numNodes = element->getNumberNodes(); const double* pts[numNodes]; - + // get coordinates of elements for(int i = 0 ; i < numNodes ; ++i) { - pts[i] = element->getCoordsOfNode(i + 1); + pts[i] = getCoordsOfNode(i + 1, elemIdx, mesh); } _box = new BoundingBox(pts, numNodes); } else { - const int numNodes = element->getNumberNodes(); - for(int i = 1 ; i <= numNodes ; ++i) + for(int i = 0 ; i < numNodes ; ++i) { - const double* pt = element->getCoordsOfNode(i); + const double* pt = getCoordsOfNode(i + 1, elemIdx, mesh); _box->updateWithPoint(pt); } } @@ -72,7 +73,7 @@ namespace INTERP_UTILS * @param coord coordinate of BoundingBox to use when splitting the region * */ - void MeshRegion::split(MeshRegion& region1, MeshRegion& region2, BoundingBox::BoxCoord coord) + void MeshRegion::split(MeshRegion& region1, MeshRegion& region2, BoundingBox::BoxCoord coord, const MEDMEM::MESH& mesh) { ElementBBoxOrder cmp(coord); @@ -86,14 +87,14 @@ namespace INTERP_UTILS while(elemCount < static_cast(_elements.size() / 2)) { - region1.addElement(*iter); + region1.addElement(*iter, mesh); ++iter; ++elemCount; } while(iter != _elements.end()) { - region2.addElement(*iter); + region2.addElement(*iter, mesh); ++iter; } diff --git a/src/INTERP_KERNEL/MeshRegion.hxx b/src/INTERP_KERNEL/MeshRegion.hxx index 62081aaad..4311ec60a 100644 --- a/src/INTERP_KERNEL/MeshRegion.hxx +++ b/src/INTERP_KERNEL/MeshRegion.hxx @@ -5,6 +5,8 @@ #include "BoundingBox.hxx" +#include "MEDMEM_Mesh.hxx" + namespace INTERP_UTILS { class MeshElement; @@ -36,7 +38,7 @@ namespace INTERP_UTILS * @param element pointer to element to add to region * */ - void addElement(MeshElement* const element); + void addElement(MeshElement* const element, const MEDMEM::MESH& mesh); /** * Splits the region in two along the given axis, copying the elements with bounding boxes whose maximum @@ -48,7 +50,7 @@ namespace INTERP_UTILS * @param axis axis along which to split the region * */ - void split(MeshRegion& region1, MeshRegion& region2, BoundingBox::BoxCoord coord); + void split(MeshRegion& region1, MeshRegion& region2, BoundingBox::BoxCoord coord, const MEDMEM::MESH& mesh); bool isDisjointWithElementBoundingBox(const MeshElement& elem) const; @@ -63,7 +65,7 @@ namespace INTERP_UTILS // neither allocated or liberated in this class. The elements must therefore be allocated and liberated outside this class std::vector _elements; BoundingBox* _box; - + }; }; diff --git a/src/INTERP_KERNEL/TetraAffineTransform.hxx b/src/INTERP_KERNEL/TetraAffineTransform.hxx index 6b781bf63..bc0bc5712 100644 --- a/src/INTERP_KERNEL/TetraAffineTransform.hxx +++ b/src/INTERP_KERNEL/TetraAffineTransform.hxx @@ -6,6 +6,8 @@ #include "VectorUtils.hxx" +#include "Log.hxx" + #undef NULL_COORD_CORRECTION namespace INTERP_UTILS @@ -19,8 +21,8 @@ namespace INTERP_UTILS TetraAffineTransform(const double** pts) { - // std::cout << "Creating transform from tetraeder : " << std::endl; - // std::cout << vToStr(pts[0]) << ", " << vToStr(pts[1]) << ", " << vToStr(pts[2]) << ", " << vToStr(pts[3]) << ", " << std::endl; + LOG(2,"Creating transform from tetraeder : "); + LOG(2, vToStr(pts[0]) << ", " << vToStr(pts[1]) << ", " << vToStr(pts[2]) << ", " << vToStr(pts[3])); #if 0 do { @@ -37,7 +39,7 @@ namespace INTERP_UTILS calculateDeterminant(); - // std::cout << "determinant before inverse = " << _determinant << std::endl; + LOG(3, "determinant before inverse = " << _determinant); // check that tetra is non-planar -> determinant is not zero // otherwise set _determinant to zero to signal caller that transformation did not work @@ -78,20 +80,20 @@ namespace INTERP_UTILS calculateDeterminant(); // self-check - // std::cout << "transform determinant is " << _determinant << std::endl; - // std::cout << "*Self-check : Applying transformation to original points ... "; + LOG(4, "transform determinant is " << _determinant); + LOG(4, "*Self-check : Applying transformation to original points ... "); for(int i = 0; i < 4 ; ++i) { double v[3]; apply(v, pts[i]); - // std::cout << vToStr(v) << std::endl; + LOG(4, vToStr(v)) for(int j = 0; j < 3; ++j) { assert(epsilonEqual(v[j], (3*i+j == 3 || 3*i+j == 7 || 3*i+j == 11 ) ? 1.0 : 0.0)); } } - // std::cout << " ok" << std::endl; + LOG(4, " ok"); } void apply(double* destPt, const double* srcPt) const @@ -107,7 +109,7 @@ namespace INTERP_UTILS // alloc temporary memory dest = new double[3]; - //// std::cout << "Oops! self-affectation" << std::endl; + LOG(6, "Oops! self-affectation"); } for(int i = 0 ; i < 3 ; ++i) @@ -147,7 +149,7 @@ namespace INTERP_UTILS { using namespace std; - // std::cout << "A = " << std::endl << "["; + std::cout << "A = " << std::endl << "["; for(int i = 0; i < 3; ++i) { cout << _linearTransform[3*i] << ", " << _linearTransform[3*i + 1] << ", " << _linearTransform[3*i + 2]; @@ -192,7 +194,7 @@ namespace INTERP_UTILS int(i == 2) }; - //std::cout << "b = [" << b[0] << ", " << b[1] << ", " << b[2] << "]" << std::endl; + LOG(6, "b = [" << b[0] << ", " << b[1] << ", " << b[2] << "]"); double y[3]; forwardSubstitution(y, lu, b, idx); diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index 6f2a9650e..49e7e6013 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -5,13 +5,15 @@ #include #include #include -#include "VectorUtils.hxx" #include #include +#include "VectorUtils.hxx" + #undef MERGE_CALC #undef COORDINATE_CORRECTION 1.0e-15 + class CircularSortOrder { public: @@ -27,8 +29,8 @@ public: _a = barycenter[_aIdx] - pt0[_aIdx]; _b = barycenter[_bIdx] - pt0[_bIdx]; - // std::cout << "Creating order of type " << type << " with pt0= " << vToStr(pt0) << std::endl; - // std::cout << "a = " << _a << ", b = " << _b << std::endl; + //LOG(4, "Creating order of type " << type << " with pt0= " << vToStr(pt0)); + //LOG(4, "a = " << _a << ", b = " << _b) } bool operator()(const double* pt1, const double* pt2) @@ -94,7 +96,7 @@ class Vector3Cmp public: bool operator()(double* const& pt1, double* const& pt2) { - // std::cout << "points are equal ? : " << int((pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2])) << std::endl; + LOG(6, "points are equal ? : " << int((pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2]))); return (pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2]); } }; @@ -137,7 +139,7 @@ namespace INTERP_UTILS /* _coords[5*P + 3] = 1 - (p[0] - (p[1] - p[2])); _coords[5*Q + 3] = 1 - (q[0] - (q[1] - q[2])); - std::cout << "old = " << 1 - q[0] - q[1] - q[2] << " calculated = " << 1 - (q[0] - (q[1] - q[2])) << " stored : " << _coords[5*Q + 3] << std::endl; + LOG(6, "old = " << 1 - q[0] - q[1] - q[2] << " calculated = " << 1 - (q[0] - (q[1] - q[2])) << " stored : " << _coords[5*Q + 3]); _coords[5*R + 3] = 1 - (r[0] -(r[1] - r[2])); */ @@ -189,7 +191,7 @@ namespace INTERP_UTILS if(isTriangleBelowTetraeder()) { - // std::cout << std::endl << "Triangle is below tetraeder - V = 0.0" << std::endl << std::endl ; + LOG(2, " --- Triangle is below tetraeder - V = 0.0"); return 0.0; } @@ -207,7 +209,7 @@ namespace INTERP_UTILS if(sign == 0.0) { - // std::cout << std::endl << "Triangle is perpendicular to z-plane - V = 0.0" << std::endl << std::endl; + LOG(2, " --- Triangle is perpendicular to z-plane - V = 0.0"); return 0.0; } @@ -216,7 +218,7 @@ namespace INTERP_UTILS sign = sign > 0.0 ? 1.0 : -1.0; - // std::cout << std::endl << "-- Calculating intersection polygons ... " << std::endl; + LOG(2, "-- Calculating intersection polygons ... "); calculateIntersectionPolygons(); #ifdef MERGE_CALC @@ -236,11 +238,11 @@ namespace INTERP_UTILS double volA = 0.0; if(_polygonA.size() > 2) { - // std::cout << std::endl << "-- Treating polygon A ... " << std::endl; + LOG(2, "---- Treating polygon A ... "); calculatePolygonBarycenter(A, barycenter); sortIntersectionPolygon(A, barycenter); volA = calculateVolumeUnderPolygon(A, barycenter); - // std::cout << "Volume is " << sign * volA << std::endl; + LOG(2, "Volume is " << sign * volA); } double volB = 0.0; @@ -251,16 +253,15 @@ namespace INTERP_UTILS if(!isTriangleInPlaneOfFacet(XYZ) && _polygonB.size() > 2) #endif { - // std::cout << std::endl << "-- Treating polygon B ... " << std::endl; - // std::cout << _coords[5*P + 3] << ", " << _coords[5*Q + 3] << ", " << _coords[5*R+ 3] << std::endl; + LOG(2, "---- Treating polygon B ... "); calculatePolygonBarycenter(B, barycenter); sortIntersectionPolygon(B, barycenter); volB = calculateVolumeUnderPolygon(B, barycenter); - // std::cout << "Volume is " << sign * volB << std::endl; + LOG(2, "Volume is " << sign * volB); } - // std::cout << std::endl << "volA + volB = " << sign * (volA + volB) << std::endl << std::endl; + LOG(2, "volA + volB = " << sign * (volA + volB) << std::endl << "***********"); return sign * (volA + volB); @@ -296,13 +297,13 @@ namespace INTERP_UTILS double* ptA = new double[3]; calcIntersectionPtSurfaceEdge(edge, ptA); _polygonA.push_back(ptA); - // std::cout << "Surface-edge : " << vToStr(ptA) << " added to A " << std::endl; + LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A "); if(edge >= XY) { double* ptB = new double[3]; copyVector3(ptA, ptB); _polygonB.push_back(ptB); - // std::cout << "Surface-edge : " << vToStr(ptB) << " added to B " << std::endl; + LOG(3,"Surface-edge : " << vToStr(ptB) << " added to B "); } } @@ -316,7 +317,7 @@ namespace INTERP_UTILS double* ptB = new double[3]; copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); _polygonB.push_back(ptB); - // std::cout << "Surface-ray : " << vToStr(ptB) << " added to B" << std::endl; + LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B"); } } @@ -331,13 +332,13 @@ namespace INTERP_UTILS double* ptA = new double[3]; calcIntersectionPtSegmentFacet(seg, facet, ptA); _polygonA.push_back(ptA); - // std::cout << "Segment-facet : " << vToStr(ptA) << " added to A" << std::endl; + LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A"); if(facet == XYZ) { double* ptB = new double[3]; copyVector3(ptA, ptB); _polygonB.push_back(ptB); - // std::cout << "Segment-facet : " << vToStr(ptB) << " added to B" << std::endl; + LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B"); } } } @@ -350,7 +351,7 @@ namespace INTERP_UTILS double* ptA = new double[3]; calcIntersectionPtSegmentEdge(seg, edge, ptA); _polygonA.push_back(ptA); - // std::cout << "Segment-edge : " << vToStr(ptA) << " added to A" << std::endl; + LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A"); if(edge >= XY) { double* ptB = new double[3]; @@ -368,13 +369,13 @@ namespace INTERP_UTILS double* ptA = new double[3]; copyVector3(&COORDS_TET_CORNER[3 * corner], ptA); _polygonA.push_back(ptA); - // std::cout << "Segment-corner : " << vToStr(ptA) << " added to A" << std::endl; + 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); - // std::cout << "Segment-corner : " << vToStr(ptB) << " added to B" << std::endl; + LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B"); } } } @@ -387,7 +388,7 @@ namespace INTERP_UTILS double* ptB = new double[3]; copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); _polygonB.push_back(ptB); - // std::cout << "Segment-ray : " << vToStr(ptB) << " added to B" << std::endl; + LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B"); } } @@ -399,7 +400,7 @@ namespace INTERP_UTILS double* ptB = new double[3]; calcIntersectionPtSegmentHalfstrip(seg, edge, ptB); _polygonB.push_back(ptB); - // std::cout << "Segment-halfstrip : " << vToStr(ptB) << " added to B" << std::endl; + LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B"); } } } @@ -414,7 +415,7 @@ namespace INTERP_UTILS double* ptA = new double[3]; copyVector3(&_coords[5*corner], ptA); _polygonA.push_back(ptA); - // std::cout << "Inclusion tetrahedron : " << vToStr(ptA) << " added to A" << std::endl; + LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A"); } // XYZ - plane @@ -423,7 +424,7 @@ namespace INTERP_UTILS double* ptB = new double[3]; copyVector3(&_coords[5*corner], ptB); _polygonB.push_back(ptB); - // std::cout << "Inclusion XYZ-plane : " << vToStr(ptB) << " added to B" << std::endl; + LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B"); } // projection on XYZ - facet @@ -434,7 +435,7 @@ namespace INTERP_UTILS ptB[2] = 1 - ptB[0] - ptB[1]; assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0)); _polygonB.push_back(ptB); - // std::cout << "Projection XYZ-plane : " << vToStr(ptB) << " added to B" << std::endl; + LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B"); } } @@ -452,7 +453,7 @@ namespace INTERP_UTILS */ void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter) { - // std::cout << "--- Calculating polygon barycenter" << std::endl; + LOG(3,"--- Calculating polygon barycenter"); // get the polygon points std::vector& polygon = (poly == A) ? _polygonA : _polygonB; @@ -476,7 +477,7 @@ namespace INTERP_UTILS } } } - // std::cout << "Barycenter is " << vToStr(barycenter) << std::endl; + LOG(3,"Barycenter is " << vToStr(barycenter)); } /** @@ -492,7 +493,7 @@ namespace INTERP_UTILS */ void TransformedTriangle::sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter) { - // std::cout << "--- Sorting polygon ..."<< std::endl; + LOG(3,"--- Sorting polygon ..."); using ::ProjectedCentralCircularSortOrder; typedef ProjectedCentralCircularSortOrder SortOrder; // change is only necessary here and in constructor @@ -529,10 +530,10 @@ namespace INTERP_UTILS //stable_sort((++polygon.begin()), polygon.end(), order); - // std::cout << "Sorted polygon is " << std::endl; + LOG(3,"Sorted polygon is "); for(int i = 0 ; i < polygon.size() ; ++i) { - // std::cout << vToStr(polygon[i]) << std::endl; + LOG(3,vToStr(polygon[i])); } } @@ -550,7 +551,7 @@ namespace INTERP_UTILS */ double TransformedTriangle::calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter) { - // std::cout << "--- Calculating volume under polygon" << std::endl; + LOG(2,"--- Calculating volume under polygon"); // get the polygon points std::vector& polygon = (poly == A) ? _polygonA : _polygonB; @@ -571,7 +572,7 @@ namespace INTERP_UTILS vol += (factor1 * factor2) / 6.0; } - // std::cout << "Abs. Volume is " << vol << std::endl; + LOG(2,"Abs. Volume is " << vol); return vol; } @@ -654,12 +655,12 @@ namespace INTERP_UTILS void TransformedTriangle::dumpCoords() { - // std::cout << "Coords : "; + std::cout << "Coords : "; for(int i = 0 ; i < 3; ++i) { - // std::cout << vToStr(&_coords[5*i]) << ","; + std::cout << vToStr(&_coords[5*i]) << ","; } - // std::cout << std::endl; + std::cout << std::endl; } }; // NAMESPACE diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 6797310eb..55e8c5aec 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -3,6 +3,14 @@ #include +// Levels : +// 1 - overview of algorithm + volume result +// 2 - algorithm detail +// 3 - intersection polygon results detail +// 4 - intersection polygon search detail +// higher -> misc. gory details of calculation +//#define LOG_LEVEL 3 +#include "Log.hxx" class TransformedTriangleTest; class TransformedTriangleIntersectTest; @@ -148,6 +156,8 @@ namespace INTERP_UTILS void preCalculateDoubleProducts(void); + bool areDoubleProductsConsistent(const TriSegment seg) const; + void resetDoubleProducts(const TriSegment seg, const TetraCorner corner); double calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const; diff --git a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx index 300ee72c0..35afa8721 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx @@ -109,15 +109,15 @@ namespace INTERP_UTILS const double alpha = tA / (tA - tB); // calculate point - // std::cout << "corner A = " << corners[0] << " corner B = " << corners[1] << std::endl; - // std::cout << "tA = " << tA << " tB = " << tB << " alpha= " << alpha << std::endl; + LOG(4, "corner A = " << corners[0] << " corner B = " << corners[1] ); + LOG(4, "tA = " << tA << " tB = " << tB << " alpha= " << alpha ); for(int i = 0; i < 3; ++i) { pt[i] = (1 - alpha) * COORDS_TET_CORNER[3*corners[0] + i] + alpha * COORDS_TET_CORNER[3*corners[1] + i]; - // // std::cout << pt[i] << std::endl; - assert(pt[i] >= 0.0); - assert(pt[i] <= 1.0); + LOG(6, pt[i] ); + assert(pt[i] >= 0.0); + assert(pt[i] <= 1.0); } } @@ -171,8 +171,8 @@ namespace INTERP_UTILS const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx]; const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx]; pt[i] = -( sign * calcStableC(seg, dp) ) / s; - // std::cout << "SegmentFacetIntPtCalc : pt[" << i << "] = " << pt[i] << std::endl; - // std::cout << "c(" << seg << ", " << dp << ") = " << sign * calcStableC(seg, dp) << std::endl; + LOG(4, "SegmentFacetIntPtCalc : pt[" << i << "] = " << pt[i] ); + LOG(4, "c(" << seg << ", " << dp << ") = " << sign * calcStableC(seg, dp) ); assert(pt[i] >= 0.0); assert(pt[i] <= 1.0); } @@ -448,8 +448,8 @@ namespace INTERP_UTILS // special case : facet H = 0 bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet); - // std::cout << "Halfstrip tests (" << seg << ", " << edge << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) << std::endl; - // std::cout << "c2 = " << cVals[2] << ", c3 = " << cVals[3] << std::endl; + LOG(4, "Halfstrip tests (" << seg << ", " << edge << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) ); + LOG(4, "c2 = " << cVals[2] << ", c3 = " << cVals[3] ); return (cVals[0]*cVals[1] < 0.0) && cond2 && (cVals[2]*cVals[3] > 0.0); } @@ -506,7 +506,7 @@ namespace INTERP_UTILS }; pt[i] = (1 - alpha) * cornerCoords[0] + alpha * cornerCoords[1]; - // std::cout << pt[i] << std::endl; + LOG(6, pt[i] ); assert(pt[i] >= 0.0); assert(pt[i] <= 1.0); } @@ -524,7 +524,7 @@ namespace INTERP_UTILS bool TransformedTriangle::testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const { assert(corner == X || corner == Y || corner == Z); - // std::cout << "Testing seg - ray intersection for seg = " << seg << ", corner = " << corner << std::endl; + LOG(4, "Testing seg - ray intersection for seg = " << seg << ", corner = " << corner ); // readjust index since O is not used const int cornerIdx = static_cast(corner) - 1; @@ -567,7 +567,7 @@ namespace INTERP_UTILS // cond. 1 if(cVal0 != 0.0) { - // std::cout << "SR fails at cond 1 cVal0 = " << cVal0 << std::endl; + LOG(4, "SR fails at cond 1 cVal0 = " << cVal0 ); return false; } @@ -577,7 +577,7 @@ namespace INTERP_UTILS if(!(cond21 || cond22)) { - // std::cout << "SR fails at cond 2 : cond21 = " << cond21 << ", cond22 = " << cond22 << std::endl; + LOG(4, "SR fails at cond 2 : cond21 = " << cond21 << ", cond22 = " << cond22 ); return false; } @@ -595,7 +595,7 @@ namespace INTERP_UTILS // cond. 3 if(( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) >= 0.0) { - // std::cout << "SR fails at cond 3 : " << (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] << std::endl; + LOG(4, "SR fails at cond 3 : " << (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ); } return ( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) < 0.0; @@ -693,7 +693,7 @@ namespace INTERP_UTILS const double cQR = calcStableC(QR, DoubleProduct(edge)); const double cRP = calcStableC(RP, DoubleProduct(edge)); - // std::cout << "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" << std::endl; + LOG(5, "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" ); // if two or more c-values are zero we disallow x-edge intersection // Grandy, p.446 @@ -702,7 +702,7 @@ namespace INTERP_UTILS if(numZeros >= 2 ) { - // std::cout << "TriangleSurroundsEdge test fails due to too many 0 dp" << std::endl; + LOG(5, "TriangleSurroundsEdge test fails due to too many 0 dp" ); } return (cPQ*cQR >= 0.0) && (cQR*cRP >= 0.0) && (cRP*cPQ >= 0.0) && numZeros < 2; @@ -737,7 +737,7 @@ namespace INTERP_UTILS const double t2 = calcStableT(TRIPLE_PRODUCTS[2*edge + 1]); //? should equality with zero use epsilon? - // std::cout << "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 << std::endl; + LOG(5, "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 ); return (t1*t2 <= 0.0) && (t1 - t2 != 0.0); } @@ -791,7 +791,7 @@ namespace INTERP_UTILS const double coord2 = _coords[5*( (seg + 1) % 3) + facet]; //? should we use epsilon-equality here in second test? - // std::cout << "coord1 : " << coord1 << " coord2 : " << coord2 << std::endl; + LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 ); //return (coord1*coord2 <= 0.0) && epsilonEqual(coord1,coord2); return (coord1*coord2 <= 0.0) && (coord1 != coord2); } @@ -802,7 +802,7 @@ namespace INTERP_UTILS const double coord1 = _coords[5*seg + 4]; const double coord2 = _coords[5*( (seg + 1) % 3) + 4]; //? should we use epsilon-equality here in second test? - // std::cout << "coord1 : " << coord1 << " coord2 : " << coord2 << std::endl; + LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 ); return (coord1*coord2 <= 0.0) && (coord1 != coord2); } @@ -819,8 +819,8 @@ namespace INTERP_UTILS //? is it always YZ here ? //? changed to XY ! const double normal = calcStableC(PQ, C_XY) + calcStableC(QR, C_XY) + calcStableC(RP, C_XY); - //std::cout << "surface above corner " << corner << " : " << "n = " << normal << ", t = [" << calcTByDevelopingRow(corner, 1, false) << ", " << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) << std::endl; - //std::cout << "] - stable : " << calcStableT(corner) << std::endl; + LOG(6, "surface above corner " << corner << " : " << "n = " << normal << ", t = [" << calcTByDevelopingRow(corner, 1, false) << ", " << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) ); + LOG(6, "] - stable : " << calcStableT(corner) ); //? we don't care here if the triple product is "invalid", that is, the triangle does not surround one of the // edges going out from the corner (Grandy [53]) @@ -863,7 +863,7 @@ namespace INTERP_UTILS //? NB here we have no correction for precision - is this good? // Our authority Grandy says nothing - // std::cout << "dp in triSurrRay for corner " << corner << " = [" << cPQ << ", " << cQR << ", " << cRP << "]" << std::endl; + LOG(5, "dp in triSurrRay for corner " << corner << " = [" << cPQ << ", " << cQR << ", " << cRP << "]" ); return ( cPQ*cQR > 0.0 ) && ( cPQ*cRP > 0.0 ); } diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index e003671f6..dd3081a5e 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -12,7 +12,6 @@ #undef SECOND_CORNER_RESET #undef FIXED_DELTA 1.0e-15 - namespace INTERP_UTILS { @@ -76,27 +75,10 @@ namespace INTERP_UTILS // -- and make corrections if not for(TriSegment seg = PQ ; seg <= RP ; seg = TriSegment(seg + 1)) { - const double term1 = _doubleProducts[8*seg + C_YZ] * _doubleProducts[8*seg + C_XH]; - const double term2 = _doubleProducts[8*seg + C_ZX] * _doubleProducts[8*seg + C_YH]; - const double term3 = _doubleProducts[8*seg + C_XY] * _doubleProducts[8*seg + C_ZH]; - - // std::cout << std::endl; - // std::cout << "for seg " << seg << " consistency " << term1 + term2 + term3 << std::endl; - // std::cout << "term1 :" << term1 << " term2 :" << term2 << " term3: " << term3 << std::endl; - - // if(term1 + term2 + term3 != 0.0) - const int num_zero = (term1 == 0.0 ? 1 : 0) + (term2 == 0.0 ? 1 : 0) + (term3 == 0.0 ? 1 : 0); - const int num_neg = (term1 < 0.0 ? 1 : 0) + (term2 < 0.0 ? 1 : 0) + (term3 < 0.0 ? 1 : 0); - - // calculated geometry is inconsistent if we have one of the following cases - // * one term zero and the other two of the same sign - // * two terms zero - // * all terms positive - // * all terms negative - if((num_zero == 1 && num_neg != 1) || num_zero == 2 || (num_neg == 0 && num_zero != 3) || num_neg == 3 ) + if(!areDoubleProductsConsistent(seg)) { std::map distances; - //std::cout << "inconsistent! "<< std::endl; + LOG(4, "inconsistent! "); for(TetraCorner corner = O ; corner <= Z ; corner = TetraCorner(corner + 1)) { @@ -150,8 +132,8 @@ namespace INTERP_UTILS { if(_doubleProducts[8*seg + dp] != 0.0) { - // std::cout << "Double product for (seg,dp) = (" << seg << ", " << dp << ") = " - // << std::abs(_doubleProducts[8*seg + dp]) << " is imprecise, reset to 0.0" << std::endl; + LOG(5, "Double product for (seg,dp) = (" << seg << ", " << dp << ") = " ); + LOG(5, std::abs(_doubleProducts[8*seg + dp]) << " is imprecise, reset to 0.0" ); } _doubleProducts[8*seg + dp] = 0.0; @@ -162,6 +144,37 @@ namespace INTERP_UTILS _isDoubleProductsCalculated = true; } + bool TransformedTriangle::areDoubleProductsConsistent(const TriSegment seg) const + { + const double term1 = _doubleProducts[8*seg + C_YZ] * _doubleProducts[8*seg + C_XH]; + const double term2 = _doubleProducts[8*seg + C_ZX] * _doubleProducts[8*seg + C_YH]; + const double term3 = _doubleProducts[8*seg + C_XY] * _doubleProducts[8*seg + C_ZH]; + + + LOG(6, "for seg " << seg << " consistency " << term1 + term2 + term3 ); + LOG(6, "term1 :" << term1 << " term2 :" << term2 << " term3: " << term3 ); + + const int num_zero = (term1 == 0.0 ? 1 : 0) + (term2 == 0.0 ? 1 : 0) + (term3 == 0.0 ? 1 : 0); + const int num_neg = (term1 < 0.0 ? 1 : 0) + (term2 < 0.0 ? 1 : 0) + (term3 < 0.0 ? 1 : 0); + const int num_pos = (term1 > 0.0 ? 1 : 0) + (term2 > 0.0 ? 1 : 0) + (term3 > 0.0 ? 1 : 0); + + assert( num_zero + num_neg + num_pos == 3 ); + + // calculated geometry is inconsistent if we have one of the following cases + // * one term zero and the other two of the same sign + // * two terms zero + // * all terms positive + // * all terms negative + if(((num_zero == 1 && num_neg != 1) || num_zero == 2 || (num_neg == 0 && num_zero != 3) || num_neg == 3 )) + { + LOG(4, "inconsistent dp found" ); + } + return !((num_zero == 1 && num_neg != 1) || num_zero == 2 || (num_neg == 0 && num_zero != 3) || num_neg == 3 ); + + // return (num_zero == 3) || (num_neg != 0 && num_neg != 3 && num_pos != 0 && num_pos != 3); + // return epsilonEqual(term1 + term2 + term3, 0.0); + } + void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const TetraCorner corner) { // set the three corresponding double products to 0.0 @@ -176,7 +189,7 @@ namespace INTERP_UTILS for(int i = 0 ; i < 3 ; ++i) { const DoubleProduct dp = DOUBLE_PRODUCTS[3*corner + i]; - // std::cout << std::endl << "resetting inconsistent dp :" << dp << " for corner " << corner << std::endl; + LOG(6, std::endl << "resetting inconsistent dp :" << dp << " for corner " << corner); _doubleProducts[8*seg + dp] = 0.0; }; } @@ -227,10 +240,10 @@ namespace INTERP_UTILS if(_isTripleProductsCalculated) return; - // std::cout << "Precalculating triple products" << std::endl; + LOG(4, "Precalculating triple products" ); for(TetraCorner corner = O ; corner <= Z ; corner = TetraCorner(corner + 1)) { - // std::cout << "- Triple product for corner " << corner << std::endl; + LOG(6, "- Triple product for corner " << corner ); // find edge / row to use -> that whose edge makes the smallest angle to the triangle // use a map to find the minimum @@ -273,7 +286,7 @@ namespace INTERP_UTILS { // this value will not be used // we set it to whatever - // std::cout << "Triple product not calculated for corner " << corner << std::endl; + LOG(6, "Triple product not calculated for corner " << corner ); _tripleProducts[corner] = -3.14159265; _validTP[corner] = false; @@ -514,14 +527,14 @@ namespace INTERP_UTILS if( epsilonEqual( p_term + q_term + r_term, 0.0, THRESHOLD_F * delta) ) { - // std::cout << "Reset imprecise triple product for corner " << corner << " to zero" << std::endl; + LOG(4, "Reset imprecise triple product for corner " << corner << " to zero" ); return 0.0; } else { // NB : using plus also for the middle term compensates for a double product // which is inversely ordered - // std::cout << "Triple product for corner " << corner << ", row " << row << " = " << sign*( p_term + q_term + r_term ) << std::endl; + LOG(6, "Triple product for corner " << corner << ", row " << row << " = " << sign*( p_term + q_term + r_term ) ); return sign*( p_term + q_term + r_term ); } -- 2.39.2