}
}
+ 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
{
{
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;
}
}
void updateWithPoint(const double* pt);
+ void dumpCoords() const;
+
private:
bool isValid() const;
#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<MeshElement*, int> indices;
+ std::map<MeshElement*, int> 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<RegionNode*> 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<RegionNode*> 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<MeshElement*>::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<MeshElement*>::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<BoundingBox::BoxCoord>(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<MeshElement*>::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<BoundingBox::BoxCoord>(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<MeshElement*>::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<int> 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<int>::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<TransformedTriangle> triangles;
- srcElement.triangulate(triangles, T);
+ // triangulate source element faces (assumed tetraeder for the time being)
+ // do nothing
+ vector<TransformedTriangle> 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<TransformedTriangle>::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<TransformedTriangle>::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
}
* @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);
};
--- /dev/null
+#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<TransformedTriangle> triangles;
+ triangulate(srcType, srcCell, _srcMesh, triangles, T);
+
+ double volume = 0.0;
+
+ LOG(4, "num triangles = " << triangles.size());
+ int i = 0;
+ for(vector<TransformedTriangle>::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<TransformedTriangle>& 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;
+ ;
+ }
+ }
+
+
+ }
+
+
+
+
+};
--- /dev/null
+#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 <vector>
+
+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<TransformedTriangle>& triangles, const TetraAffineTransform& T) const;
+
+ private:
+ const MEDMEM::MESH& _srcMesh;
+ const MEDMEM::MESH& _targetMesh;
+
+ };
+};
+
+
+#endif
Interpolation2D.hxx\
Interpolation3DSurf.hxx\
InterpolationUtils.hxx\
-BBTree.H
+BBTree.H\
+MeshUtils.hxx\
+Intersector3D.hxx\
+Log.hxx
# Libraries targets
RegionNode.cxx\
BoundingBox.cxx\
TetraAffineTransform.cxx\
-Interpolation3D.cxx
+Interpolation3D.cxx\
+Intersector3D.cxx\
+#Log.cxx
#Interpolation3DSurf.cxx\
#Interpolation2D.cxx\
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.
#include "TetraAffineTransform.hxx"
#include "TransformedTriangle.hxx"
#include "MEDMEM_CellModel.hxx"
+#include "MeshUtils.hxx"
+#include "BoundingBox.hxx"
+
+
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
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<int>(_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<TransformedTriangle>& 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)
return coord1 < coord2;
}
-
};
#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.
*/
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<TransformedTriangle>& triangles, const TetraAffineTransform& T) const;
+ // void triangulate(std::vector<TransformedTriangle>& 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;
};
class ElementBBoxOrder
{
public :
-
+
ElementBBoxOrder(BoundingBox::BoxCoord coord)
: _coord(coord)
{
}
-
+
bool operator()(MeshElement* elem1, MeshElement* elem2);
-
+
private :
BoundingBox::BoxCoord _coord;
};
};
+
+
+
#endif
#include "MeshElement.hxx"
-//#include <math.h>
+#include "MeshUtils.hxx"
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);
}
}
* @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);
while(elemCount < static_cast<int>(_elements.size() / 2))
{
- region1.addElement(*iter);
+ region1.addElement(*iter, mesh);
++iter;
++elemCount;
}
while(iter != _elements.end())
{
- region2.addElement(*iter);
+ region2.addElement(*iter, mesh);
++iter;
}
#include "BoundingBox.hxx"
+#include "MEDMEM_Mesh.hxx"
+
namespace INTERP_UTILS
{
class MeshElement;
* @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
* @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;
// neither allocated or liberated in this class. The elements must therefore be allocated and liberated outside this class
std::vector<MeshElement*> _elements;
BoundingBox* _box;
-
+
};
};
#include "VectorUtils.hxx"
+#include "Log.hxx"
+
#undef NULL_COORD_CORRECTION
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 {
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
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
// 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)
{
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];
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);
#include <algorithm>
#include <functional>
#include <iterator>
-#include "VectorUtils.hxx"
#include <math.h>
#include <vector>
+#include "VectorUtils.hxx"
+
#undef MERGE_CALC
#undef COORDINATE_CORRECTION 1.0e-15
+
class CircularSortOrder
{
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)
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]);
}
};
/*
_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]));
*/
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;
}
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;
}
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
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;
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);
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 ");
}
}
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");
}
}
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");
}
}
}
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];
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");
}
}
}
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");
}
}
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");
}
}
}
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
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
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");
}
}
*/
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<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
}
}
}
- // std::cout << "Barycenter is " << vToStr(barycenter) << std::endl;
+ LOG(3,"Barycenter is " << vToStr(barycenter));
}
/**
*/
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
//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]));
}
}
*/
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<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
vol += (factor1 * factor2) / 6.0;
}
- // std::cout << "Abs. Volume is " << vol << std::endl;
+ LOG(2,"Abs. Volume is " << vol);
return vol;
}
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
#include <vector>
+// 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;
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;
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);
}
}
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);
}
// 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);
}
};
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);
}
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<int>(corner) - 1;
// 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;
}
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;
}
// 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;
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
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;
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);
}
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);
}
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);
}
//? 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])
//? 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 );
}
#undef SECOND_CORNER_RESET
#undef FIXED_DELTA 1.0e-15
-
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<double, TetraCorner> distances;
- //std::cout << "inconsistent! "<< std::endl;
+ LOG(4, "inconsistent! ");
for(TetraCorner corner = O ; corner <= Z ; corner = TetraCorner(corner + 1))
{
{
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;
_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
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;
};
}
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
{
// 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;
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 );
}
