#include <iostream>
#include <algorithm>
+#include <cassert>
namespace INTERP_UTILS
{
-
- /**
- * Default constructor
- *
- */
- //BoundingBox() : coords({0.0, 0.0, 0.0, 0.0, 0.0, 0.0}) {}
/**
* Constructor creating box from an array of the points corresponding
*
*/
BoundingBox::BoundingBox(const double** pts, const int numPts)
+ :_coords(new double[6])
{
using namespace std;
+ assert(_coords != 0);
assert(numPts > 1);
+
// initialize with first two points
const double* pt1 = pts[0];
}
/**
- * Constructor creating box from union of two boxes
+ * Constructor creating box from union of two boxes, resulting in a box that enclose both of them
*
+ * @param box1 the first box
+ * @param box2 the second box
*/
BoundingBox::BoundingBox(const BoundingBox& box1, const BoundingBox& box2)
+ : _coords(new double[6])
{
using namespace std;
+ assert(_coords != 0);
+
for(BoxCoord c = XMIN ; c <= ZMIN ; c = BoxCoord(c + 1))
{
_coords[c] = min(box1._coords[c], box2._coords[c]);
*/
BoundingBox::~BoundingBox()
{
+ delete[] _coords;
}
/**
* Determines if the intersection with a given box is empty
*
- * @param box box with which intersection is tested
+ * @param box BoundingBox with which intersection is tested
* @returns true if intersection between boxes is empty, false if not
*/
bool BoundingBox::isDisjointWith(const BoundingBox& box) const
{
for(BoxCoord c = XMIN ; c <= ZMIN ; c = BoxCoord(c + 1))
{
-
const double otherMinCoord = box.getCoordinate(c);
const double otherMaxCoord = box.getCoordinate(BoxCoord(c + 3));
// boxes are disjoint if there exists a direction in which the
// minimum coordinate of one is greater than the maximum coordinate of the other
- //? stable version
+ //? more stable version
/*const double tol = 1.0e-2*_coords[c];
if(_coords[c] > otherMaxCoord + tol
|| _coords[c + 3] < otherMinCoord - tol)
}
}
+ /*
+ * Prints the coordinates of the box to std::cout
+ *
+ */
void BoundingBox::dumpCoords() const
{
std::cout << "[xmin, xmax] = [" << _coords[XMIN] << ", " << _coords[XMAX] << "]" << " | ";
std::cout << std::endl;
}
+ /*
+ * Checks if the box is valid, which it is if its minimum coordinates are
+ * smaller than its maximum coordinates in all directions.
+ *
+ * @returns true if the box is valid, false if not
+ */
bool BoundingBox::isValid() const
{
bool valid = true;
for(BoxCoord c = XMIN ; c < ZMIN ; c = BoxCoord(c + 1))
{
- if(_coords[c] >= _coords[c + 3])
+ if(_coords[c] > _coords[c + 3])
{
std::cout << "+++ Error in BoundingBox |: coordinate " << c << " is invalid : "
- <<_coords[c] << " >= " << _coords[c+3] << std::endl;
+ <<_coords[c] << " > " << _coords[c+3] << std::endl;
valid = false;
}
}
{
public:
enum BoxCoord { XMIN = 0, YMIN = 1, ZMIN = 2, XMAX = 3, YMAX = 4, ZMAX = 5 };
-
-
- /**
- * Default constructor
- *
- */
- //BoundingBox() : coords({0.0, 0.0, 0.0, 0.0, 0.0, 0.0}) {}
-
- /**
- * Constructor creating box from an array of the points corresponding
- * to the vertices of the element.
- * Each point is represented by an array of three doubles.
- *
- * @param pts array of points
- * @param numPts number of vertices
- *
- */
+
BoundingBox(const double** pts, const int numPts);
- /**
- * Constructor creating box from union of two boxes
- *
- */
BoundingBox(const BoundingBox& box1, const BoundingBox& box2);
- /**
- * Destructor
- *
- */
~BoundingBox();
- /**
- * Determines if the intersection with a given box is empty
- *
- * @param box box with which intersection is tested
- * @returns true if intersection between boxes is empty, false if not
- */
bool isDisjointWith(const BoundingBox& box) const;
- /**
- * Sets a coordinate of the box to a given value.
- *
- * @param coord coordinate to set
- * @param value new value for coordinate
- *
- */
void setCoordinate(const BoxCoord coord, double value);
- /**
- * Gets a coordinate of the box
- *
- * @param coord coordinate to set
- * @returns value of coordinate
- *
- */
double getCoordinate(const BoxCoord coord) const;
void updateWithPoint(const double* pt);
bool isValid() const;
- double _coords[6];
-
-
+ /// vector containing the coordinates of the box
+ /// interlaced in the order XMIN, YMIN, ZMIN, XMAX, YMAX, ZMAX
+ double* _coords;
};
- // typedef BoundingBox::BoxCoord BoxCoord;
-
};
#endif
*/
Interpolation3D::Interpolation3D()
{
- // not implemented
}
/**
*/
Interpolation3D::~Interpolation3D()
{
- // not implemented
}
/**
//} 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 do more sanity checking here - eliminating meshes that
+ // cannot be treated
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.
- *
+ * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh.
+ * The calculation is done in two steps. First a filtering process reduces the number of pairs of elements for which the
+ * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
+ * volume of intersection is calculated by an object of type Intersector3D for the remaining pairs, and entered into the
+ * intersection matrix.
+ *
+ * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
+ * The length of the vector is equal to the number of source elements - for each source element there is a map, regardless
+ * of whether the element intersects any target elements or not. But in the maps there are only entries for those target elements
+ * which have a non-zero intersection volume with the source element. The vector has indices running from
+ * 0 to (#source elements - 1), meaning that the map for source element i is stored at index i - 1. In the maps, however,
+ * the indexing is more natural : the intersection volume of the source element with target element j is found at matrix[i-1][j].
+ *
+
* @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)
+ void Interpolation3D::calculateIntersectionVolumes(const MESH& srcMesh, const 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 intersector element
+ Intersector3D* intersector = new Intersector3D(srcMesh, targetMesh);
// 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);
- Intersector3D* intersector = new Intersector3D(srcMesh, targetMesh);
-
LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
- // create empty maps for all source elements
- matrix.resize(numSrcElems);
-
-
MeshElement* srcElems[numSrcElems];
MeshElement* targetElems[numTargetElems];
-
+
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 < 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);
}
+
+ // create empty maps for all source elements
+ matrix.resize(numSrcElems);
#ifdef USE_RECURSIVE_BBOX_FILTER
-
+
+ /*
+ * 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 MeshElement contains a bounding box and the global number of the corresponding element in the mesh.
+ */
+
// 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
+ // its target region with all the target mesh elements whose bounding box
+ // intersects that of the source region
RegionNode* firstNode = new RegionNode();
}
// 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
+ // right part of the source region of the current node (created 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.
if(currNode->getSrcRegion().getNumberOfElements() == 1)
{
+ // calculate volumes
LOG(4, " - One element");
- // volume calculation
MeshElement* srcElement = *(currNode->getSrcRegion().getBeginElements());
// NB : srcElement indices are from 0 .. numSrcElements - 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))
+ //if(!epsilonEqual(vol, 0.0))
+ if(vol != 0.0)
{
volumes->insert(make_pair(targetIdx, vol));
LOG(2, "Result : V (" << srcIdx << ", " << targetIdx << ") = " << matrix[srcIdx - 1][targetIdx]);
}
}
+
}
else // recursion
{
// 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
+ // Anyway, it basically chooses the next axis, cyclically
axis = (axis != BoundingBox::ZMAX) ? static_cast<BoundingBox::BoxCoord>(axis + 1) : BoundingBox::XMAX;
- // add target elements of curr node that overlap the two new nodes
+ // add target elements of current node that overlap the source regions of the new nodes
LOG(5, " -- Adding target elements");
int numLeftElements = 0;
int numRightElements = 0;
LOG(5, "Left target region has " << numLeftElements << " elements and right target region has "
<< numRightElements << " elements");
+ // push new nodes on stack
if(numLeftElements != 0)
{
nodes.push(leftNode);
for(vector<int>::const_iterator iter = intersectElems.begin() ; iter != intersectElems.end() ; ++iter)
{
+#if 0
const int srcIdx = *iter + 1;
const double vol = intersector->intersectCells(srcIdx, targetIdx);
{
matrix[srcIdx - 1].insert(make_pair(targetIdx, vol));
}
+#endif
}
}
}
- /**
- * 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;
- }
-
- // 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);
-
- 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 << "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 )
-
- //? 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) ;
-
- }
-#endif
}
#include <vector>
#include <map>
-// NOTE : Method comments should be moved to source file during the implementation phase
// typedefs
typedef std::vector< std::map< int, double > > IntersectionMatrix;
public :
- /**
- * Default constructor
- *
- */
Interpolation3D();
-
-
- /**
- * Destructor
- *
- */
virtual ~Interpolation3D();
- /**
- * 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
- */
virtual IntersectionMatrix interpol_maillages(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh);
private:
- /**
- * 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 SearchNodes, 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 calculateIntersectionVolumes(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, IntersectionMatrix& matrix);
- /**
- * 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 calculateIntersectionVolume(const INTERP_UTILS::MeshElement& srcElement, const INTERP_UTILS::MeshElement& targetElement);
-
};
};
#include "Intersector3D.hxx"
-
+#include "TetraAffineTransform.hxx"
+#include "TransformedTriangle.hxx"
+#include "MeshUtils.hxx"
+#include "VectorUtils.hxx"
#include "Log.hxx"
+#include <cmath>
+#include <cassert>
+
+using namespace MEDMEM;
+using namespace MED_EN;
+using namespace INTERP_UTILS;
+
namespace MEDMEM
{
-
- Intersector3D::Intersector3D(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh)
+ /*
+ * Constructor
+ *
+ * @param srcMesh mesh containing the source elements
+ * @param targetMesh mesh containing the target elements
+ *
+ */
+ Intersector3D::Intersector3D(const MESH& srcMesh, const MESH& targetMesh)
: _srcMesh(srcMesh), _targetMesh(targetMesh)
{
}
+ /*
+ * Destructor
+ *
+ */
Intersector3D::~Intersector3D()
{
}
+ /*
+ * Calculates the volume of intersection of an element in the source mesh and an element
+ * in the target mesh. The method is based on the algorithm of Grandy. It first calculates the transformation
+ * that takes the target tetrahedron into the unit tetrahedron. After that, the
+ * faces of the source element are triangulated and the calculated transformation is applied
+ * to each triangle. The algorithm of Grandy, implemented in INTERP_UTILS::TransformedTriangle is used
+ * to calculate the contribution to the volume from each triangle. The volume returned is the sum of these contributions
+ * divided by the determinant of the transformation.
+ *
+ * @pre The element in _targetMesh referenced by targetCell is of type MED_TETRA4.
+ * @param srcCell global number of the source element (1 <= srcCell < # source cells)
+ * @param targetCell global number of the target element (1 <= targetCell < # target cells) - this element must be a tetrahedron
+ */
double Intersector3D::intersectCells(int srcCell, int targetCell)
{
-
+ // get type of cell
const medGeometryElement srcType = _srcMesh.getElementType(MED_CELL, srcCell);
- // const medGeometryElement targetType = _targetMesh.getElementType(MED_CELL, targetCell);
+ const medGeometryElement targetType = _targetMesh.getElementType(MED_CELL, targetCell);
+
+ // maybe we should do something more civilized here
+ assert(targetType == MED_TETRA4);
// get array of points of target tetraeder
const double* tetraCorners[4];
tetraCorners[i] = getCoordsOfNode(i + 1, targetCell, _targetMesh);
}
- // create AffineTransform
+ // create AffineTransform from tetrahedron
TetraAffineTransform T( tetraCorners );
// check if we have planar tetra element
LOG(2, "Planar tetra -- volume 0");
return 0.0;
}
- // triangulate source element faces (assumed tetraeder for the time being)
- // do nothing
+
+ // triangulate source element faces
vector<TransformedTriangle> triangles;
- triangulate(srcType, srcCell, _srcMesh, triangles, T);
+ triangulate(srcType, srcCell, triangles, T);
double volume = 0.0;
LOG(4, "num triangles = " << triangles.size());
int i = 0;
+
+ // add up volumes
+ // ? potential stability issue - maybe volumes should be added in increasing order of magnitude?
for(vector<TransformedTriangle>::iterator iter = triangles.begin() ; iter != triangles.end(); ++iter)
{
LOG(2, "= > Triangle " << ++i);
-#ifdef LOG_ACTIVE
- if(LOG_LEVEL >= 2)
- {
- iter->dumpCoords();
- }
+
+#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
+ // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation
+ // that should be used (which is equivalent to dividing by the determinant)
return std::abs(1.0 / T.determinant() * volume) ;
-
}
/**
- * Triangulate the faces of an element and apply an affine Transform to the triangles
+ * Triangulate the faces of a source 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
+ * @param type type of the element
+ * @param element global number of the element
+ * @param triangles vector in which the transformed 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
+#ifdef OPTIMIZE_ // optimized version
+ void Intersector3D::triangulate(const medGeometryElement type, const int element, std::vector<TransformedTriangle>& triangles, const TetraAffineTransform& T) const
{
+
// get cell model for the element
CELLMODEL cellModel(type);
assert(cellModel.getDimension() == 3);
+ assert(element >= 1);
+ assert(element <= _srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS));
+
+ // try to avoid reallocations with push_back()
+ // 2 * the number of faces should be plenty
+ triangles.reserve(2 * cellModel.getNumberOfConstituents(1));
+
+ // loop over faces
+ const int numFaces = cellModel.getNumberOfConstituents(1);
+ double* transformedNodeList[numFaces];
+ medGeometryElement faceTypes[numFaces];
+ bool isOutsideTetra = true;
+
+ for(int i = 1 ; i <= numFaces ; ++i)
+ {
+ faceTypes[i - 1] = cellModel.getConstituentType(1, i);
+ CELLMODEL faceModel(faceTypes[i-1]);
+
+ assert(faceModel.getDimension() == 2);
+ // assert(faceModel.getNumberOfNodes() == 3);
+
+ transformedNodeList[i - 1] = new double[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);
- // start index in connectivity array for cell
- // const int cellIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index] - 1;
+ assert(localNodeNumber >= 1);
+ assert(localNodeNumber <= cellModel.getNumberOfNodes());
+
+ const double* node = getCoordsOfNode(localNodeNumber, element, _srcMesh);
+
+ // transform
+ //{ not totally efficient since we transform each node once per face
+ T.apply(&transformedNodeList[i - 1][3*(j-1)], node);
+ //T.apply(&transformedNodes[3*(j-1)], node);
+
+ LOG(4, "Node " << localNodeNumber << " = " << vToStr(node) << " transformed to "
+ << vToStr(&transformedNodeList[i - 1][3*(j-1)]));
+
+ }
- // assert(cellIdx >= 0);
- //assert(cellIdx < _mesh->getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS));
+ isOutsideTetra = isOutsideTetra && isTriangleOutsideTetra(transformedNodeList[i - 1]);
+ }
+
+ if(!isOutsideTetra)
+ {
+ for(int i = 1 ; i <= numFaces ; ++i)
+ {
+
+ // create transformed triangles from face
+ switch(faceTypes[i - 1])
+ {
+ case MED_TRIA3:
+ // simple take the triangle as it is
+ triangles.push_back(TransformedTriangle(&transformedNodeList[i - 1][0], &transformedNodeList[i - 1][3], &transformedNodeList[i - 1][6]));
+ break;
+
+ case MED_QUAD4:
+ // simple triangulation of faces along a diagonal :
+ //
+ // 2 ------ 3
+ // | / |
+ // | / |
+ // | / |
+ // | / |
+ // | / |
+ // | / |
+ // 1 ------ 4
+ //
+ //? not sure if this always works
+
+ // local nodes 1, 2, 3
+ triangles.push_back(TransformedTriangle(&transformedNodeList[i - 1][0], &transformedNodeList[i - 1][3], &transformedNodeList[i - 1][6]));
+
+ // local nodes 1, 3, 4
+ triangles.push_back(TransformedTriangle(&transformedNodeList[i - 1][0], &transformedNodeList[i - 1][6], &transformedNodeList[i - 1][9]));
+
+ break;
+
+ default:
+ std::cout << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." << std::endl;
+ assert(false);
+ }
+ }
+ }
+ for(int i = 1 ; i <= numFaces ; ++i)
+ {
+ delete[] transformedNodeList[i-1];
+ }
+ }
+
+#else // un-optimized version
+
+ void Intersector3D::triangulate(const medGeometryElement type, const int element, std::vector<TransformedTriangle>& triangles, const TetraAffineTransform& T) const
+ {
+ // get cell model for the element
+ CELLMODEL cellModel(type);
+
+ assert(element >= 1);
+ assert(element <= _srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS));
+ assert(cellModel.getDimension() == 3);
// loop over faces
for(int i = 1 ; i <= cellModel.getNumberOfConstituents(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)
assert(localNodeNumber >= 1);
assert(localNodeNumber <= cellModel.getNumberOfNodes());
- const double* node = getCoordsOfNode(localNodeNumber, element, mesh);
+ const double* node = getCoordsOfNode(localNodeNumber, element, _srcMesh);
// 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)]));
+ << vToStr(&transformedNodes[3*(j-1)]));
}
- // to be removed
- assert(faceType == MED_TRIA3);
-
// create transformed triangles from face
switch(faceType)
{
break;
// add other cases here to treat hexahedra, pyramides, etc
-
+ case MED_QUAD4:
+ // simple triangulation of faces along a diagonal :
+ //
+ // 2 ------ 3
+ // | / |
+ // | / |
+ // | / |
+ // | / |
+ // | / |
+ // | / |
+ // 1 ------ 4
+ //
+ //? not sure if this always works
+
+ // local nodes 1, 2, 3
+ triangles.push_back(TransformedTriangle(&transformedNodes[0], &transformedNodes[3], &transformedNodes[6]));
+
+ // local nodes 1, 3, 4
+ triangles.push_back(TransformedTriangle(&transformedNodes[0], &transformedNodes[6], &transformedNodes[9]));
+
+ break;
default:
std::cout << "+++ Error : Only elements with triangular faces are supported at the moment." << std::endl;
;
}
}
-
-
}
+
+#endif
-
-
+
};
#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 INTERP_UTILS
+{
+ class TransformedTriangle;
+ class TetraAffineTransform;
+};
namespace MEDMEM
{
-
+ /*
+ * Class calculating the volume of intersection of individual 3D elements.
+ *
+ */
class Intersector3D : public Intersector
{
public :
- Intersector3D(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh);
+ Intersector3D(const MESH& srcMesh, const MESH& targetMesh);
- ~Intersector3D();
+ virtual ~Intersector3D();
- double intersectCells(int srcCell, int targetCell);
+ virtual 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;
+
+ inline bool isTriangleOutsideTetra(double* points) const;
+
+ void triangulate(const MED_EN::medGeometryElement type, const int element, std::vector<INTERP_UTILS::TransformedTriangle>& triangles, const INTERP_UTILS::TetraAffineTransform& T) const;
+
+ const MESH& _srcMesh;
+ const MESH& _targetMesh;
};
+
+ inline bool Intersector3D::isTriangleOutsideTetra(double* points) const
+ {
+ const bool leftX = points[0] <= 0.0 && points[3] <= 0.0 && points[6] <= 0.0;
+ const bool rightX = points[0] >= 1.0 && points[3] >= 1.0 && points[6] >= 1.0;
+ const bool leftY = points[1] <= 0.0 && points[4] <= 0.0 && points[7] <= 0.0;
+ const bool rightY = points[1] >= 1.0 && points[4] >= 1.0 && points[7] >= 1.0;
+ const bool leftZ = points[2] <= 0.0 && points[5] <= 0.0 && points[8] <= 0.0;
+ const bool rightZ = points[2] >= 1.0 && points[5] >= 1.0 && points[8] >= 1.0;
+ const double h[3] =
+ {
+ 1 - points[0] - points[1] - points[2],
+ 1 - points[3] - points[4] - points[5],
+ 1 - points[6] - points[7] - points[8]
+ };
+ const bool leftH = h[0] <= 0.0 && h[1] <= 0.0 && h[2] <= 0.0;
+ const bool rightH = h[0] >= 1.0 && h[1] >= 1.0 && h[2] >= 1.0;
+
+ return leftX || rightX || leftY || rightY || leftZ || rightZ || leftH || rightH;
+ }
+
};
#define _LOG_H_
/* Simple pre-processor logging utility.
- * Replaces LOG( x ) with "if(x <= LOG_LEVEL) std::cout" when logging is active (LOG_ACTIVE defined)
- * x is the level at which the message should be logged.
+ * Replaces LOG( lvl, x ) with "if(lvl <= LOG_LEVEL) std::cout << x << std::endl" when logging is active
+ * (LOG_LEVEL > 0 is defined). x is the level at which the message should be logged - if it is smaller or equal to
+ * LOG_LEVEL (which can be defined at compile-time for each file by passing option -DLOG_LEVEL=x to gcc)
+ * than the message is logged.
*
*
*
*/
+/// define LOG_LEVEL here if it is not already defined
+#ifndef LOG_LEVEL
+#define LOG_LEVEL 0
+#endif
-
-#define LOG_ACTIVE
-
-#ifdef LOG_ACTIVE
+#if LOG_LEVEL > 0
#include <iostream>
-#define LOG_LEVEL 1
-
+/// write message msg to std::cout if x <= LOG_LEVEL
#define LOG(x, msg) if(x <= LOG_LEVEL) std::cout << msg << std::endl;
#else
CXXFLAGS+=@CXXTMPDPTHFLAGS@ $(MED2_INCLUDES)
CPPFLAGS+=$(BOOST_CPPFLAGS)
-CXXFLAGS+=-fprofile-arcs -ftest-coverage
-CPPFLAGS+=-fprofile-arcs -ftest-coverage
+# optimization
+#CXXFLAGS+= -O1 -DOPTIMIZE
+#CPPFLAGS+= -O1 -DOPTIMIZE
+CXXFLAGS+= -DOPTIMIZE -O2
+CPPFLAGS+= -DOPTIMIZE -O2
+
+# for log
+CXXFLAGS+= -DLOG_LEVEL=0
+CPPFLAGS+= -DLOG_LEVEL=0
+
+# for gcov
+#CXXFLAGS+=-fprofile-arcs -ftest-coverage
+#CPPFLAGS+=-fprofile-arcs -ftest-coverage
+
+#for gprof
+CXXFLAGS+=-pg
+CPPFLAGS+=-pg
#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
+# for gcov
+#LDFLAGS+= -lgcov
+LDFLAGS+= -pg
#LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS)
# change motivated by the bug KERNEL4778.
LDFLAGSFORBIN+= -lm $(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 -lmedmem $(BOOST_LIBS) -lutil
#LDFLAGSFORBIN+= -lm $(HDF5_LIBS) -lmedmem $(BOOST_LIBS) -lutil
+
ifeq ($(MED_WITH_KERNEL),yes)
CPPFLAGS+= ${KERNEL_CXXFLAGS}
CXXFLAGS+= ${KERNEL_CXXFLAGS}
LDFLAGSFORBIN+= ${KERNEL_LDFLAGS} -lSALOMELocalTrace -lSALOMEBasics
endif
+LDFLAGSFORBIN+= -pg
+
LIBSFORBIN=$(BOOSTLIBS) $(MPI_LIBS)
LIBS=
#include "BoundingBox.hxx"
-
namespace INTERP_UTILS
{
/**
* Constructor
*
+ * @param index global number of element in the mesh
* @param mesh mesh that the element belongs to
* @param type geometric type of the element
- * @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), _box(0), _type(type)
{
// get coordinates of vertices
const int numNodes = getNumberOfNodesForType(type);
-
+
+ assert(numNodes >= 3);
+
const double* vertices[numNodes];
for(int i = 0 ; i < numNodes ; ++i)
}
}
- /**
- * Determines if this element is in the interior of another element
- * by calculating the triple products for each point of this element with respect
- * to all the faces of the other object (faces must be triangulated ... ) If all triple
- * products have the same sign, then the element is in the interior of the other element
+ /*
+ * Accessor to global number
*
- * @param otherElement the supposedly enclosing element
- * @returns true if this element is enclosed in the other element, false if not
+ * @return global number of the element
*/
- bool MeshElement::isElementIncludedIn(const MeshElement& otherElement) const
+ int MeshElement::getIndex() const
{
- // not implemented
- return false;
- }
+ return _index;
+ }
+
+ /*
+ * Accessor to bounding box
+ *
+ * @return pointer to bounding box of the element
+ */
+ const BoundingBox* MeshElement::getBoundingBox() const
+ {
+ return _box;
+ }
- /**
- * Determines whether the intersection of this element is trivially empty. This is done by checking for each
- * face of one element if it is such that all the vertices of the other element is on the same side of this face.
- * If there is such a face, then the intersection is trivially empty. If there is no such face, we do not know if
- * the intersection is empty.
+ /*
+ * Accessor to the type of the element
*
- * @pre The elements are convex. If this is not true, we return false.
- * @param otherElement the element believed to be disjoint with this one
- * @returns true if the two elements are convex and there exists a face of this element such as described
- * above, false otherwise
+ * @return type of the element
*/
- bool MeshElement::isElementTriviallyDisjointWith(const MeshElement& otherElement) const
- {
- // not implemented
- return false;
- }
+ MED_EN::medGeometryElement MeshElement::getType() const
+ {
+ return _type;
+ }
-
- int MeshElement::getIndex() const
- {
- return _index;
- }
-
- void MeshElement::dumpCoords() const
+ /////////////////////////////////////////////////////////////////////
+ /// ElementBBoxOrder /////////////
+ /////////////////////////////////////////////////////////////////////
+ /*
+ * Constructor
+ *
+ * @param coord BoundingBox coordinate (XMIN, XMAX, etc) on which to base the ordering
+ */
+ ElementBBoxOrder::ElementBBoxOrder(BoundingBox::BoxCoord coord)
+ : _coord(coord)
{
- std::cout << "Bounding box of element " << _index << " is " << std::endl;
- _box->dumpCoords();
}
-
- /// ElementBBoxOrder
+ /*
+ * Comparison operator based on the bounding boxes of the elements
+ *
+ * @returns true if the coordinate _coord of the bounding box of elem1 is
+ * strictly smaller than that of the bounding box of elem2
+ */
bool ElementBBoxOrder::operator()( MeshElement* elem1, MeshElement* elem2)
{
+ const BoundingBox* box1 = elem1->getBoundingBox();
+ const BoundingBox* box2 = elem2->getBoundingBox();
+
assert(elem1 != 0);
assert(elem2 != 0);
- assert(elem1->_box != 0);
- assert(elem2->_box != 0);
+ assert(box1 != 0);
+ assert(box2 != 0);
- const double coord1 = elem1->_box->getCoordinate(_coord);
- const double coord2 = elem2->_box->getCoordinate(_coord);
+ const double coord1 = box1->getCoordinate(_coord);
+ const double coord2 = box2->getCoordinate(_coord);
return coord1 < coord2;
}
namespace INTERP_UTILS
{
- // class TransformedTriangle;
- // class TetraAffineTransform;
-
/**
* Class representing a single element of a mesh together with its bounding box.
- * It permits access to the nodes of the element, to test for disunion and inclusion with other elements
- * and to triangulate its faces.
+ * It permits access to the element's global number, type and bounding box and allows
+ * easy bounding box intersection tests between MeshElements and collections of MeshElement (MeshRegions)
*/
class MeshElement
{
-
-
- friend class ElementBBoxOrder;
- friend class MeshRegion;
public:
- /**
- * Constructor
- *
- * @param mesh mesh that the element belongs to
- * @param type geometric type of the element
- * @param index connectivity index of element in the mesh
- */
MeshElement(const int index, const MED_EN::medGeometryElement type, const MEDMEM::MESH& mesh);
- /**
- * Destructor
- *
- */
~MeshElement();
-
- /**
- * Determines if this element is in the interior of another element
- * by calculating the triple products for each point of this element with respect
- * to all the faces of the other object (faces must be triangulated ... ) If all triple
- * products have the same sign, then the element is in the interior of the other element
- *
- * @param otherElement the supposedly enclosing element
- * @returns true if this element is enclosed in the other element, false if not
- */
- bool isElementIncludedIn(const MeshElement& otherElement) const;
-
- /**
- * Determines whether the intersection of this element is trivially empty. This is done by checking for each
- * face of one element if it is such that all the vertices of the other element is on the same side of this face.
- * If there is such a face, then the intersection is trivially empty. If there is no such face, we do not know if
- * the intersection is empty.
- *
- * @pre The elements are convex. If this is no true, we return false.
- * @param otherElement the element believed to be disjoint with this one
- * @returns true if the two elements are convex and there exists a face of this element such as described
- * above, false otherwise
- */
- bool isElementTriviallyDisjointWith(const MeshElement& otherElement) 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;
int getIndex() const;
-
- void dumpCoords() const;
- const BoundingBox* getBoundingBox() const
- {
- return _box;
- }
+ const BoundingBox* getBoundingBox() const;
- MED_EN::medGeometryElement getType() const
- {
- return _type;
- }
+ MED_EN::medGeometryElement getType() const;
private:
- const int _index;
- BoundingBox* _box; // should not change after initialisation
- MED_EN::medGeometryElement _type;
+ const int _index; /// global number of the element
+ BoundingBox* _box; /// bounding box of the element - does not change after having been initialised
+ MED_EN::medGeometryElement _type;/// type of the element
};
-
+ /*
+ * Class defining an order for MeshElements based on their bounding boxes.
+ * The order defined between two elements is that between a given coordinate of
+ * their bounding boxes. For instance, if the order is based on YMIN, an element whose boxes
+ * has a smaller YMIN is sorted before one with a larger YMIN.
+ *
+ */
class ElementBBoxOrder
{
public :
- ElementBBoxOrder(BoundingBox::BoxCoord coord)
- : _coord(coord)
- {
- }
+ ElementBBoxOrder(BoundingBox::BoxCoord coord);
bool operator()(MeshElement* elem1, MeshElement* elem2);
private :
- BoundingBox::BoxCoord _coord;
+ BoundingBox::BoxCoord _coord; /// BoundingBox coordinate (XMIN, XMAX, etc) on which to base the ordering
};
};
}
/**
- * Adds an element to the region, updating the bounding box.
+ * Adds an element to the region, updating the bounding box. If the bounding box does not yet
+ * exist, it is created here. This creation is delayed to make it possible to have empty MeshRegions
*
* @param element pointer to element to add to region
*
{
const double* pts[numNodes];
- // get coordinates of elements
+ // get coordinates of the nodes of the element
for(int i = 0 ; i < numNodes ; ++i)
{
pts[i] = getCoordsOfNode(i + 1, elemIdx, mesh);
*/
void MeshRegion::split(MeshRegion& region1, MeshRegion& region2, BoundingBox::BoxCoord coord, const MEDMEM::MESH& mesh)
{
+ // create ordering
ElementBBoxOrder cmp(coord);
// sort elements by their bounding boxes
region2.addElement(*iter, mesh);
++iter;
}
-
- // assert(std::abs(region1._elements.size() - region2._elements.size()) < 2);
-
}
+ /*
+ * Determines if a given element can intersect the elements of this region by
+ * testing whether the bounding box of the region intersects the bounding box of the element.
+ * Note that the test is only true in one direction : if the bounding boxes are disjoint, the
+ * element cannot intersect any of the elements in the region, but if they are not disjoint, the
+ * element may or may not do so.
+ *
+ * @param elem Element with which to test for disjoint-ness
+ * @return true if the bounding box of the element is disjoint with the bounding box of the region, false otherwise
+ */
bool MeshRegion::isDisjointWithElementBoundingBox(const MeshElement& elem) const
{
+ const BoundingBox* elemBox = elem.getBoundingBox();
+
assert(_box != 0);
- assert(elem._box != 0);
+ assert(elemBox != 0);
- return _box->isDisjointWith(*(elem._box));
+ return _box->isDisjointWith(*elemBox);
}
+ /*
+ * Accessor to beginning of elements vector
+ *
+ * @return constant iterator pointing at the beginning of the vector or elements
+ */
std::vector<MeshElement*>::const_iterator MeshRegion::getBeginElements() const
{
return _elements.begin();
}
+ /*
+ * Accessor to end of elements vector
+ *
+ * @return constant iterator pointing at the end of the vector or elements
+ */
std::vector<MeshElement*>::const_iterator MeshRegion::getEndElements() const
{
return _elements.end();
}
-
+
+ /*
+ * Gives information on how many elements are contained in the region.
+ *
+ * @return the number of elements contained in the region
+ */
int MeshRegion::getNumberOfElements() const
{
return _elements.size();
{
public:
- /**
- * Default constructor
- *
- */
MeshRegion();
-
- /**
- * Destructor
- *
- */
+
~MeshRegion();
- /**
- * Adds an element to the region, updating the bounding box.
- *
- * @param element pointer to element to add to region
- *
- */
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
- * coordinate along the axis are smaller than the middle of the bounding box of this region in region1. The
- * rest of the elements are copied to region2.
- *
- * @param region1 region in which to store one half of this region
- * @param region2 region in which to store the other of this region
- * @param axis axis along which to split the region
- *
- */
void split(MeshRegion& region1, MeshRegion& region2, BoundingBox::BoxCoord coord, const MEDMEM::MESH& mesh);
bool isDisjointWithElementBoundingBox(const MeshElement& elem) const;
int getNumberOfElements() const;
private:
- // Vector of pointers to elements. NB : these pointers are not owned by the region object, and are thus
- // neither allocated or liberated in this class. The elements must therefore be allocated and liberated outside this class
+ /// Vector of pointers to elements. NB : these pointers are not owned by the region object, and are thus
+ /// neither allocated or liberated in this class. The elements must therefore be allocated and liberated outside this class
std::vector<MeshElement*> _elements;
+
+ /// BoundingBox containing all the nodes of all the elements in the region.
BoundingBox* _box;
};
--- /dev/null
+#ifndef __MESH_UTILS_HXX__
+#define __MESH_UTILS_HXX__
+
+#include "MEDMEM_define.hxx"
+#include <cassert>
+
+using namespace MEDMEM;
+using namespace MED_EN;
+
+namespace INTERP_UTILS
+{
+
+ /**
+ * Returns the coordinates of a node of an element
+ * (1 <= node <= #nodes)
+ *
+ * @param node the node for which the coordinates are sought
+ * @param element an element of the mesh
+ * @param mesh a mesh
+ * @return pointer to an array of 3 doubles containing the coordinates of the node
+ */
+ inline const double* getCoordsOfNode(int node, int element, const MESH& mesh)
+ {
+ assert(node >= 1);
+ assert(node <= mesh.getNumberOfNodes());
+ const int nodeOffset = node - 1;
+ const int elemIdx = mesh.getConnectivityIndex(MED_NODAL, MED_CELL)[element - 1] - 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]);
+ }
+
+ /**
+ * Returns the number of nodes of a given type of element
+ *
+ * @param type a type of element
+ * @return the number of nodes in elements of the type
+ */
+ inline int getNumberOfNodesForType(const medGeometryElement type)
+ {
+ 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;
+ }
+
+};
+
+
+
+
+
+#endif
RegionNode::~RegionNode()
{
}
-
/**
* Accessor to source region
* Class containing a tuplet of a source region and a target region.
* This is used as the object to put on the stack in the depth-first search
* in the bounding-box filtering process.
- *
*/
class RegionNode
{
public:
- /**
- * Default constructor
- *
- */
RegionNode();
- /**
- * Destructor
- *
- */
~RegionNode();
- /**
- * Accessor to source region
- *
- * @return reference to source region
- */
MeshRegion& getSrcRegion();
- /**
- * Accessor to target region
- *
- * @return reference to target region
- */
MeshRegion& getTargetRegion();
-
private:
- MeshRegion _srcRegion;
- MeshRegion _targetRegion;
+ MeshRegion _srcRegion; /// source region
+ MeshRegion _targetRegion; /// target region
};
#include "VectorUtils.hxx"
-// Levels :
+// levels :
// 1 - titles and volume results
-// 2 - empty
-// 3 - symmetry / diagonal results
-// 4 - intersection matrix output
+// 2 - symmetry / diagonal results and intersection matrix output
+// 3 - empty
+// 4 - empty
// 5 - misc
#include "Log.hxx"
-#define VOL_PREC 1.0e-6
+//#define VOL_PREC 1.0e-6
using namespace MEDMEM;
using namespace std;
int j = iter2->first;
if(m2.at(j-1).count(i+1) == 0)
{
- if(!epsilonEqualRelative(iter2->second, 0.0, VOL_PREC))
+ if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
{
- LOG(3, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
+ LOG(2, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
+ LOG(2, "(" << i << ", " << j << ") fails");
compatitable = false;
}
}
}
++i;
}
+ if(!compatitable)
+ {
+ LOG(1, "*** matrices are not compatitable");
+ }
return compatitable;
}
const double v2 = theMap[i + 1];
if(v1 != v2)
{
- LOG(3, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 );
+ LOG(2, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 );
if(!epsilonEqualRelative(v1, v2, VOL_PREC))
{
+ LOG(2, "(" << i << ", " << j << ") fails");
isSymmetric = false;
}
}
}
++i;
}
+ if(!isSymmetric)
+ {
+ LOG(1, "*** matrices are not symmetric");
+ }
return isSymmetric;
}
const double vol = iter2->second;
if(vol != 0.0 && (i != j))
{
- LOG(3, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" );
- if(!epsilonEqualRelative(vol, 0.0, VOL_PREC))
+ LOG(2, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" );
+ if(!epsilonEqual(vol, 0.0, VOL_PREC))
{
+ LOG(2, "(" << i << ", " << j << ") fails");
isDiagonal = false;
}
}
}
++i;
}
+ if(!isDiagonal)
+ {
+ LOG(1, "*** matrix is not diagonal");
+ }
return isDiagonal;
}
}
-void Interpolation3DTest::intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec) const
+void Interpolation3DTest::intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec, bool doubleTest) const
{
LOG(1, std::endl << std::endl << "=============================" );
IntersectionMatrix matrix1;
calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
-#ifdef LOG_ACTIVE
- if(LOG_LEVEL >= 4 )
- {
- dumpIntersectionMatrix(matrix1);
- }
+#if LOG_LEVEL >= 2
+ dumpIntersectionMatrix(matrix1);
#endif
std::cout.precision(16);
const double vol1 = sumVolume(matrix1);
- if(isTestReflexive)
+ if(!doubleTest)
{
LOG(1, "vol = " << vol1 <<" correctVol = " << correctVol );
CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1));
- CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
+
+ if(isTestReflexive)
+ {
+ CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
+ }
}
else
{
+
IntersectionMatrix matrix2;
calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);
-#ifdef LOG_ACTIVE
- if(LOG_LEVEL >= 4 )
- {
- dumpIntersectionMatrix(matrix2);
- }
+#if LOG_LEVEL >= 2
+ dumpIntersectionMatrix(matrix2);
#endif
const double vol2 = sumVolume(matrix2);
CPPUNIT_TEST( tetraSimpleHalfstripOnly );
CPPUNIT_TEST( generalTetra );
CPPUNIT_TEST( trickyTetra1 );
-
CPPUNIT_TEST( inconsistentTetra );
// multi - element
-
-
+
CPPUNIT_TEST( tetraComplexIncluded );
CPPUNIT_TEST( dividedUnitTetraSimplerReflexive );
CPPUNIT_TEST( dividedUnitTetraReflexive );
CPPUNIT_TEST( moderateBoxEvenSmallerReflexive );
CPPUNIT_TEST( tinyBoxReflexive );
+ CPPUNIT_TEST( simpleHexaBox );
CPPUNIT_TEST_SUITE_END();
{
intersectMeshes("meshes/TinyBox.med", "TinyBox", "meshes/TinyBox.med", "TinyBox", 979200);
}
-
+
+ void simpleHexaBox()
+ {
+ intersectMeshes("meshes/BoxHexa1.med", "BoxHexa1", "meshes/BoxTetra2.med", "BoxTetra2", 65250, 1.0e-5, false);
+ }
+
private:
Interpolation3D* interpolator;
void dumpIntersectionMatrix(const IntersectionMatrix& m) const;
- void intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec = 1.0e-6) const;
+ // 1.0e-5 here is due to limited precision of "correct" volumes calculated in Salome
+ void intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec = 1.0e-5, bool doubleTest = true) const;
void calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const;
@COMMENCE@
# header files
-EXPORT_HEADERS = CppUnitTest.hxx \
- TransformedTriangleTest.hxx \
- TransformedTriangleIntersectTest.hxx \
- Interpolation3DTest.hxx
+EXPORT_HEADERS = TestingUtils.hxx
# Libraries targets
# Executables targets
-BIN = TestInterpKernel
+BIN = PerfTest
+
+BIN_SRC =
-BIN_SRC = CppUnitTest.cxx TransformedTriangleTest.cxx TransformedTriangleIntersectTest.cxx \
-Interpolation3DTest.cxx
BIN_CLIENT_IDL =
CPPFLAGS+=$(BOOST_CPPFLAGS) $(MED2_INCLUDES) $(HDF5_INCLUDES) -I$(top_srcdir)/src/INTERP_KERNEL
#include CppUnit
-CXXFLAGS+= -I/usr/include/cppunit
-CPPFLAGS+= -I/usr/include/cppunit
+#CXXFLAGS+= -I/usr/include/cppunit
+#CPPFLAGS+= -I/usr/include/cppunit
+
+# for log
+CXXFLAGS+= -DLOG_LEVEL=1
+CPPFLAGS+= -DLOG_LEVEL=1
+
+# for gcov
+#CXXFLAGS+=-fprofile-arcs -ftest-coverage
+#CPPFLAGS+=-fprofile-arcs -ftest-coverage
+
+#for gprof
+CXXFLAGS+=-pg
+CPPFLAGS+=-pg
#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) -lmedmem -lutil -lmed
+# for gcov
+#LDFLAGS+= -lgcov
+LDFLAGS+= -pg
+
#LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS)
# change motivated by the bug KERNEL4778.
#LDFLAGSFORBIN+= -lm $(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 -lmedmem $(BOOST_LIBS) -lutil -linterpkernel
+
ifeq ($(MED_WITH_KERNEL),yes)
CPPFLAGS+= ${KERNEL_CXXFLAGS}
CXXFLAGS+= ${KERNEL_CXXFLAGS}
LDFLAGSFORBIN += $(LDFLAGS) -lm $(MED2_LIBS) $(HDF5_LIBS) \
-L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome -lmed_V2_1 -lmedmem \
- -lcppunit -linterpkernel
+ -linterpkernel
#LDFLAGSFORBIN += $(LDFLAGS) -lm $(HDF5_LIBS) \
# -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome \
# -lcppunit -linterpkernel
-UNIT_TEST_PROG = TestInterpKernel
+
+LDFLAGSFORBIN += -pg
+
+UNIT_TEST_PROG = #PerfTest
@CONCLUDE@
--- /dev/null
+#include "Interpolation3D.hxx"
+#include "TestingUtils.hxx"
+#include <cassert>
+#include <string>
+
+int main(int argc, char** argv)
+{
+ assert(argc == 3);
+
+ // load meshes
+ const string mesh1 = argv[1];
+ const string mesh2 = argv[2];
+
+ const string mesh1path = "meshes/" + mesh1 + ".med";
+ const string mesh2path = "meshes/" + mesh2 + ".med";
+
+ IntersectionMatrix m;
+
+ calcIntersectionMatrix(mesh1path.c_str(), mesh1.c_str(), mesh2path.c_str(), mesh2.c_str(), m);
+
+ //dumpIntersectionMatrix(m);
+
+ return 0;
+
+}
// --- Registers the fixture into the 'registry'
//CPPUNIT_TEST_SUITE_REGISTRATION( Interpolation3DTestMultiElement );
-CPPUNIT_TEST_SUITE_REGISTRATION( Interpolation3DTest );
+//CPPUNIT_TEST_SUITE_REGISTRATION( Interpolation3DTest );
CPPUNIT_TEST_SUITE_REGISTRATION( TransformedTriangleIntersectTest );
-CPPUNIT_TEST_SUITE_REGISTRATION( TransformedTriangleTest );
+//CPPUNIT_TEST_SUITE_REGISTRATION( TransformedTriangleTest );
//CPPUNIT_TEST_SUITE_REGISTRATION( TestBogusClass );
// --- generic Main program from KERNEL_SRC/src/Basics/Test
--- /dev/null
+#ifndef _TESTING_UTILS_HXX_
+#define _TESTING_UTILS_HXX_
+
+#include "Interpolation3D.hxx"
+#include "MEDMEM_Mesh.hxx"
+
+#include <iostream>
+#include <map>
+#include <vector>
+#include <cmath>
+#include <algorithm>
+
+#include "VectorUtils.hxx"
+
+// levels :
+// 1 - titles and volume results
+// 2 - symmetry / diagonal results and intersection matrix output
+// 3 - empty
+// 4 - empty
+// 5 - misc
+#include "Log.hxx"
+
+using namespace MEDMEM;
+using namespace std;
+using namespace INTERP_UTILS;
+using namespace MED_EN;
+
+
+double sumVolume(const IntersectionMatrix& m)
+{
+
+ vector<double> volumes;
+ for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
+ {
+ for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
+ {
+ volumes.push_back(iter2->second);
+ // vol += std::abs(iter2->second);
+ }
+ }
+
+ // sum in ascending order to avoid rounding errors
+
+ sort(volumes.begin(), volumes.end());
+ const double vol = accumulate(volumes.begin(), volumes.end(), 0.0);
+
+ return vol;
+}
+
+#if 0
+
+bool areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2)
+{
+ bool compatitable = true;
+ int i = 0;
+ for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
+ {
+ for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
+ {
+ int j = iter2->first;
+ if(m2.at(j-1).count(i+1) == 0)
+ {
+ if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
+ {
+ LOG(2, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
+ LOG(2, "(" << i << ", " << j << ") fails");
+ compatitable = false;
+ }
+ }
+ }
+ ++i;
+ }
+ if(!compatitable)
+ {
+ LOG(1, "*** matrices are not compatitable");
+ }
+ return compatitable;
+}
+
+bool testSymmetric(const IntersectionMatrix& m1, const IntersectionMatrix& m2)
+{
+
+ int i = 0;
+ bool isSymmetric = true;
+
+ LOG(1, "Checking symmetry src - target" );
+ isSymmetric = isSymmetric & areCompatitable(m1, m2) ;
+ LOG(1, "Checking symmetry target - src" );
+ isSymmetric = isSymmetric & areCompatitable(m2, m1);
+
+ for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
+ {
+ for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
+ {
+ int j = iter2->first;
+ const double v1 = iter2->second;
+ //if(m2[j - 1].count(i+1) > 0)
+ // {
+ map<int, double> theMap = m2.at(j-1);
+ const double v2 = theMap[i + 1];
+ if(v1 != v2)
+ {
+ LOG(2, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 );
+ if(!epsilonEqualRelative(v1, v2, VOL_PREC))
+ {
+ LOG(2, "(" << i << ", " << j << ") fails");
+ isSymmetric = false;
+ }
+ }
+ }
+ ++i;
+ }
+ if(!isSymmetric)
+ {
+ LOG(1, "*** matrices are not symmetric");
+ }
+ return isSymmetric;
+}
+
+bool testDiagonal(const IntersectionMatrix& m)
+{
+ LOG(1, "Checking if matrix is diagonal" );
+ int i = 1;
+ bool isDiagonal = true;
+ for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
+ {
+ for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
+ {
+ int j = iter2->first;
+ const double vol = iter2->second;
+ if(vol != 0.0 && (i != j))
+ {
+ LOG(2, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" );
+ if(!epsilonEqual(vol, 0.0, VOL_PREC))
+ {
+ LOG(2, "(" << i << ", " << j << ") fails");
+ isDiagonal = false;
+ }
+ }
+ }
+ ++i;
+ }
+ if(!isDiagonal)
+ {
+ LOG(1, "*** matrix is not diagonal");
+ }
+ return isDiagonal;
+}
+
+#endif
+
+void dumpIntersectionMatrix(const IntersectionMatrix& m)
+{
+ int i = 0;
+ std::cout << "Intersection matrix is " << endl;
+ for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
+ {
+ for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
+ {
+
+ std::cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << endl;
+
+ }
+ ++i;
+ }
+ std::cout << "Sum of volumes = " << sumVolume(m) << std::endl;
+}
+
+std::pair<int,int> countNumberOfMatrixEntries(const IntersectionMatrix& m)
+{
+
+ int numElems = 0;
+ int numNonZero = 0;
+ for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
+ {
+ numElems += iter->size();
+ for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
+ {
+ if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
+ {
+ ++numNonZero;
+ }
+ }
+ }
+ return std::make_pair(numElems, numNonZero);
+}
+
+
+void calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m)
+{
+ const string dataDir = getenv("DATA_DIR");
+
+ LOG(1, std::endl << "=== -> intersecting src = " << mesh1 << ", target = " << mesh2 );
+
+ LOG(5, "Loading " << mesh1 << " from " << mesh1path);
+ const MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1);
+ const int numSrcElems = sMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
+ LOG(1, "Source mesh has " << numSrcElems << " elements");
+
+
+ LOG(5, "Loading " << mesh2 << " from " << mesh2path);
+ const MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
+ const int numTargetElems = tMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
+
+ LOG(1, "Target mesh has " << numTargetElems << " elements");
+
+ Interpolation3D* interpolator = new Interpolation3D();
+
+ m = interpolator->interpol_maillages(sMesh, tMesh);
+
+ std::pair<int, int> eff = countNumberOfMatrixEntries(m);
+ LOG(1, eff.first << " of " << numTargetElems * numSrcElems << " intersections calculated : ratio = "
+ << double(eff.first) / double(numTargetElems * numSrcElems));
+ LOG(1, eff.second << " non-zero elements of " << eff.first << " total : filter efficiency = "
+ << double(eff.second) / double(eff.first));
+
+ delete interpolator;
+
+ LOG(1, "Intersection calculation done. " << std::endl );
+
+}
+
+
+
+
+
+
+
+
+#if 0
+void intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec, bool doubleTest)
+{
+ LOG(1, std::endl << std::endl << "=============================" );
+
+ using std::string;
+ const string path1 = string(mesh1path) + string(mesh1);
+ const string path2 = string(mesh2path) + string(mesh2);
+
+ const bool isTestReflexive = (path1.compare(path2) == 0);
+
+ IntersectionMatrix matrix1;
+ calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
+
+#if LOG_LEVEL >= 2
+ dumpIntersectionMatrix(matrix1);
+#endif
+
+ std::cout.precision(16);
+
+ const double vol1 = sumVolume(matrix1);
+
+ if(!doubleTest)
+ {
+ LOG(1, "vol = " << vol1 <<" correctVol = " << correctVol );
+ // CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1));
+
+ if(isTestReflexive)
+ {
+ // CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
+ }
+ }
+ else
+ {
+
+ IntersectionMatrix matrix2;
+ calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);
+
+#if LOG_LEVEL >= 2
+ dumpIntersectionMatrix(matrix2);
+#endif
+
+ const double vol2 = sumVolume(matrix2);
+
+ LOG(1, "vol1 = " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol );
+
+ // CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(vol1, correctVol));
+ // CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec * std::max(vol2, correctVol));
+ // CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec * std::max(vol1, vol2));
+ // CPPUNIT_ASSERT_EQUAL_MESSAGE("Symmetry test failed", true, testSymmetric(matrix1, matrix2));
+ }
+
+}
+
+
+
+#endif
+
+
+#endif
#include "TetraAffineTransform.hxx"
+#include "VectorUtils.hxx"
+
+#include <cmath>
+#include <iostream>
+
+#include "Log.hxx"
+
+namespace INTERP_UTILS
+{
+ /////////////////////////////////////////////////////////////////////////////////////////
+ /// PUBLIC INTERFACE METHODS //////////////
+ /////////////////////////////////////////////////////////////////////////////////////////
+
+ /*
+ * Constructor
+ * Create the TetraAffineTransform object from the tetrahedron
+ * with corners specified in pts. If the tetrahedron is degenerate or almost degenerate,
+ * construction succeeds, but the determinant of the transform is set to 0.
+ *
+ * @param pts a 4x3 matrix containing 4 points (pts[0], ..., pts[3]) of 3 coordinates each
+ */
+ TetraAffineTransform::TetraAffineTransform(const double** pts)
+ {
+
+ LOG(2,"Creating transform from tetraeder : ");
+ LOG(2, vToStr(pts[0]) << ", " << vToStr(pts[1]) << ", " << vToStr(pts[2]) << ", " << vToStr(pts[3]));
+
+ // three last points -> linear transform
+ for(int i = 0; i < 3 ; ++i)
+ {
+ for(int j = 0 ; j < 3 ; ++j)
+ {
+ // NB we insert columns, not rows
+ _linearTransform[3*j + i] = (pts[i+1])[j] - (pts[0])[j];
+ }
+ }
+
+ calculateDeterminant();
+
+ 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
+ if(epsilonEqual(_determinant, 0.0))
+ {
+ _determinant = 0.0;
+ return;
+ }
+
+ // we need the inverse transform
+ invertLinearTransform();
+
+ // first point -> translation
+ // calculate here because translation takes place in "transformed space",
+ // or in other words b = -A*O where A is the linear transform
+ // and O is the position vector of the point that is mapped onto the origin
+ for(int i = 0 ; i < 3 ; ++i)
+ {
+ _translation[i] = -(_linearTransform[3*i]*(pts[0])[0] + _linearTransform[3*i+1]*(pts[0])[1] + _linearTransform[3*i+2]*(pts[0])[2]) ;
+ }
+
+ // precalculate determinant (again after inversion of transform)
+ calculateDeterminant();
+
+#ifdef INVERSION_SELF_CHECK
+ // debugging : check that applying the inversed transform to the original points
+ // gives us the unit tetrahedron
+ 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]);
+ 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));
+ }
+ }
+
+ LOG(4, " ok");
+#endif
+ }
+
+ /*
+ * Calculates the transform of point srcPt and stores the result in destPt.
+ * If destPt == srcPt, then srcPt is overwritten safely.
+ *
+ *
+ * @param destPt double[3] in which to store the transformed point
+ * @param srcPt double[3] containing coordinates of points to be transformed
+ *
+ */
+ void TetraAffineTransform::apply(double* destPt, const double* srcPt) const
+ {
+ double* dest = destPt;
+
+ // are we self-allocating ?
+ const bool selfAllocation = (destPt == srcPt);
+
+ if(selfAllocation)
+ {
+ // alloc temporary memory
+ dest = new double[3];
+
+ LOG(6, "Info : Self-affectation in TetraAffineTransform::apply");
+ }
+
+ for(int i = 0 ; i < 3 ; ++i)
+ {
+ // matrix - vector multiplication
+ dest[i] = _linearTransform[3*i] * srcPt[0] + _linearTransform[3*i + 1] * srcPt[1] + _linearTransform[3*i + 2] * srcPt[2];
+
+ // translation
+ dest[i] += _translation[i];
+ }
+
+ if(selfAllocation)
+ {
+ // copy result back to destPt
+ for(int i = 0 ; i < 3 ; ++i)
+ {
+ destPt[i] = dest[i];
+ }
+ delete[] dest;
+ }
+ }
+
+ /*
+ * Gives the determinant of the linear part of the transform
+ *
+ * @return determinant of the transform
+ *
+ */
+ double TetraAffineTransform::determinant() const
+ {
+ return _determinant;
+ }
+
+ /*
+ * Outputs to std::cout the matrix A and the vector b
+ * of the transform Ax + b
+ *
+ */
+ void TetraAffineTransform::dump() const
+ {
+ using namespace std;
+
+ std::cout << "A = " << std::endl << "[";
+ for(int i = 0; i < 3; ++i)
+ {
+ std::cout << _linearTransform[3*i] << ", " << _linearTransform[3*i + 1] << ", " << _linearTransform[3*i + 2];
+ if(i != 2 ) std::cout << endl;
+ }
+ std::cout << "]" << endl;
+
+ std::cout << "b = " << "[" << _translation[0] << ", " << _translation[1] << ", " << _translation[2] << "]" << endl;
+ }
+
+ /////////////////////////////////////////////////////////////////////////////////////////
+ /// PRIVATE METHODS //////////////
+ /////////////////////////////////////////////////////////////////////////////////////////
+
+ /*
+ * Calculates the inverse of the matrix A, stored in _linearTransform
+ * by LU-factorization and substitution
+ *
+ */
+ void TetraAffineTransform::invertLinearTransform()
+ {
+ //{ we copy the matrix for the lu-factorization
+ // maybe inefficient
+ double lu[9];
+ for(int i = 0 ; i < 9; ++i)
+ {
+ lu[i] = _linearTransform[i];
+ }
+
+ // calculate LU factorization
+ int idx[3];
+ factorizeLU(lu, idx);
+
+ // calculate inverse by forward and backward substitution
+ // store in _linearTransform
+ // NB _linearTransform cannot be overwritten with lu in the loop
+ for(int i = 0 ; i < 3 ; ++i)
+ {
+ // form standard base vector i
+ const double b[3] =
+ {
+ int(i == 0),
+ int(i == 1),
+ int(i == 2)
+ };
+
+ LOG(6, "b = [" << b[0] << ", " << b[1] << ", " << b[2] << "]");
+
+ double y[3];
+ forwardSubstitution(y, lu, b, idx);
+
+ double x[3];
+ backwardSubstitution(x, lu, y, idx);
+
+ // copy to _linearTransform matrix
+ // NB : this is a column operation, so we cannot
+ // do this directly when we calculate x
+ for(int j = 0 ; j < 3 ; j++)
+ {
+ _linearTransform[3*j + i] = x[idx[j]];
+ }
+ }
+ }
+
+ /*
+ * Updates the member _determinant of the matrix A of the transformation.
+ *
+ */
+ void TetraAffineTransform::calculateDeterminant()
+ {
+ const double subDet[3] =
+ {
+ _linearTransform[4] * _linearTransform[8] - _linearTransform[5] * _linearTransform[7],
+ _linearTransform[3] * _linearTransform[8] - _linearTransform[5] * _linearTransform[6],
+ _linearTransform[3] * _linearTransform[7] - _linearTransform[4] * _linearTransform[6]
+ };
+
+ _determinant = _linearTransform[0] * subDet[0] - _linearTransform[1] * subDet[1] + _linearTransform[2] * subDet[2];
+ }
+
+ /*
+ * Calculates the LU-factorization of the matrix A (_linearTransform)
+ * and stores it in lu. Since partial pivoting is used, there are
+ * row swaps. This is represented by the index permutation vector idx : to access element
+ * (i,j) of lu, use lu[3*idx[i] + j]
+ *
+ * @param lu double[9] in which to store LU-factorization
+ * @param idx int[3] in which to store row permutation vector
+ */
+ void TetraAffineTransform::factorizeLU(double* lu, int* idx) const
+ {
+ // 3 x 3 LU factorization
+ // initialise idx
+ for(int i = 0 ; i < 3 ; ++i)
+ {
+ idx[i] = i;
+ }
+
+ for(int k = 0; k < 2 ; ++k)
+ {
+
+ // find pivot
+ int i = k;
+ double max = std::abs(lu[3*idx[k] + k]);
+ int row = i;
+ while(i < 3)
+ {
+ if(std::abs(lu[3*idx[i] + k]) > max)
+ {
+ max = fabs(lu[3*idx[i] + k]);
+ row = i;
+ }
+ ++i;
+ }
+
+ // swap rows in index vector
+ int tmp = idx[k];
+ idx[k] = idx[row];
+ idx[row] = tmp;
+
+ // calculate row
+ for(int j = k + 1 ; j < 3 ; ++j)
+ {
+ // l_jk = u_jk / u_kk
+ lu[3*idx[j] + k] /= lu[3*idx[k] + k];
+ for(int s = k + 1; s < 3 ; ++s)
+ {
+ // case s = k will always become zero, and then be replaced by
+ // the l-value
+ // there's no need to calculate this explicitly
+
+ // u_js = u_js - l_jk * u_ks
+ lu[3*idx[j] + s] -= lu[3*idx[j] + k] * lu[3*idx[k] + s] ;
+ }
+ }
+ }
+ }
+
+ /*
+ * Solves the system Lx = b, where L is lower unit-triangular (ones on the diagonal)
+ *
+ * @param x double[3] in which the solution is stored
+ * @param lu double[9] containing the LU-factorization
+ * @param b double[3] containing the right-hand side
+ * @param idx int[3] containing the permutation vector associated with lu
+ *
+ */
+ void TetraAffineTransform::forwardSubstitution(double* x, const double* lu, const double* b, const int* idx) const
+ {
+ // divisions are not carried out explicitly because lu does not store
+ // the diagonal values of L, which are all 1
+ x[idx[0]] = ( b[idx[0]] );// / lu[3*idx[0]]
+ x[idx[1]] = ( b[idx[1]] - lu[3*idx[1]] * x[idx[0]] ); // / lu[3*idx[1] + 1];
+ x[idx[2]] = ( b[idx[2]] - lu[3*idx[2]] * x[idx[0]] - lu[3*idx[2] + 1] * x[idx[1]] ); // / lu[3*idx[2] + 2];
+ }
+
+ /*
+ * Solves the system Ux = b, where U is upper-triangular
+ *
+ * @param x double[3] in which the solution is stored
+ * @param lu double[9] containing the LU-factorization
+ * @param b double[3] containing the right-hand side
+ * @param idx int[3] containing the permutation vector associated with lu
+ *
+ */
+ void TetraAffineTransform::backwardSubstitution(double* x, const double* lu, const double* b, const int* idx) const
+ {
+ x[idx[2]] = ( b[idx[2]] ) / lu[3*idx[2] + 2];
+ x[idx[1]] = ( b[idx[1]] - lu[3*idx[1] + 2] * x[idx[2]] ) / lu[3*idx[1] + 1];
+ x[idx[0]] = ( b[idx[0]] - lu[3*idx[0] + 1] * x[idx[1]] - lu[3*idx[0] + 2] * x[idx[2]] ) / lu[3*idx[0]];
+ }
+
+
+};
#ifndef __TETRA_AFFINE_TRANSFORM_HXX__
#define __TETRA_AFFINE_TRANSFORM_HXX__
-#include <math.h>
-#include <iostream>
-
-#include "VectorUtils.hxx"
-
-#include "Log.hxx"
-
-#undef NULL_COORD_CORRECTION
+#undef INVERSION_SELF_CHECK // debugging : check that calculated inverse is correct
namespace INTERP_UTILS
{
+ /*
+ * Class representing an affine transformation x -> Ax + b that transforms a given tetrahedron
+ * into the unit tetrahedron.
+ *
+ */
class TetraAffineTransform
{
-
public:
- TetraAffineTransform(const double** pts)
- {
-
- LOG(2,"Creating transform from tetraeder : ");
- LOG(2, vToStr(pts[0]) << ", " << vToStr(pts[1]) << ", " << vToStr(pts[2]) << ", " << vToStr(pts[3]));
-
-#if 0
- do {
-#endif
- // three last points -> linear transform
- for(int i = 0; i < 3 ; ++i)
- {
- for(int j = 0 ; j < 3 ; ++j)
- {
- // NB we insert columns, not rows
- _linearTransform[3*j + i] = (pts[i+1])[j] - (pts[0])[j];
- }
- }
-
- calculateDeterminant();
-
- 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
- if(epsilonEqual(_determinant, 0.0))
- {
- _determinant = 0.0;
- return;
- }
-
-#if 0
-
- //assert(_determinant != 0.0);
-
- if(_determinant < 0.0)
- {
- // swap two columns
- const double* tmp = pts[1];
- pts[1] = pts[2];
- pts[2] = tmp;
- }
-
- } while(_determinant < 0.0); // should do at most one more loop
-#endif
-
- // we need the inverse transform
- invertLinearTransform();
-
- // first point -> translation
- // calculate here because translation takes place in "transformed space",
- // or in other words b = -A*O where A is the linear transform
- // and O is the position vector of the point that is mapped onto the origin
- for(int i = 0 ; i < 3 ; ++i)
- {
- _translation[i] = -(_linearTransform[3*i]*(pts[0])[0] + _linearTransform[3*i+1]*(pts[0])[1] + _linearTransform[3*i+2]*(pts[0])[2]) ;
- }
-
- // precalculate determinant (again after inversion of transform)
- calculateDeterminant();
-
- // self-check
- 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]);
- 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));
- }
- }
-
- LOG(4, " ok");
- }
-
- void apply(double* destPt, const double* srcPt) const
- {
- double* dest = destPt;
-
- const bool selfAllocation = (destPt == srcPt);
-
-
-
- if(selfAllocation)
- {
- // alloc temporary memory
- dest = new double[3];
-
- LOG(6, "Oops! self-affectation");
- }
-
- for(int i = 0 ; i < 3 ; ++i)
- {
- // matrix - vector multiplication
- dest[i] = _linearTransform[3*i] * srcPt[0] + _linearTransform[3*i + 1] * srcPt[1] + _linearTransform[3*i + 2] * srcPt[2];
+ TetraAffineTransform(const double** pts);
- // translation
- dest[i] += _translation[i];
+ void apply(double* destPt, const double* srcPt) const;
-#ifdef NULL_COORD_CORRECTION
- if(epsilonEqual(dest[i], 0.0, 1.0e-15))
- {
- dest[i] = 0.0;
- }
-#endif
-
- }
-
- if(selfAllocation)
- {
- // copy result back to destPt
- for(int i = 0 ; i < 3 ; ++i)
- {
- destPt[i] = dest[i];
- }
- delete[] dest;
- }
- }
+ double determinant() const;
- double determinant() const
- {
- return _determinant;
- }
-
- void dump()
- {
- using namespace std;
-
- std::cout << "A = " << std::endl << "[";
- for(int i = 0; i < 3; ++i)
- {
- cout << _linearTransform[3*i] << ", " << _linearTransform[3*i + 1] << ", " << _linearTransform[3*i + 2];
- if(i != 2 ) cout << endl;
- }
- cout << "]" << endl;
-
- cout << "b = " << "[" << _translation[0] << ", " << _translation[1] << ", " << _translation[2] << "]" << endl;
- }
+ void dump() const;
private:
- void invertLinearTransform()
- {
- //{ we copy the matrix for the lu-factorization
- // maybe inefficient
-
- double lu[9];
- for(int i = 0 ; i < 9; ++i)
- {
- lu[i] = _linearTransform[i];
- }
-
- // calculate LU factorization
- int idx[3];
- // double parity = 1.0;
- factorizeLU(lu, idx/*, &parity*/);
-
- //_determinant = parity / (lu[3*idx[0]] * lu[3*idx[1] + 1] * lu[3*idx[2] + 2]);
- //std::cout << "lu-determinant = " << _determinant << std::endl;
-
- // calculate inverse by forward and backward substitution
- // store in _linearTransform
- // NB _linearTransform cannot be overwritten with lu in the loop
- for(int i = 0 ; i < 3 ; ++i)
- {
- // form standard base vector i
- const double b[3] =
- {
- int(i == 0),
- int(i == 1),
- int(i == 2)
- };
-
- LOG(6, "b = [" << b[0] << ", " << b[1] << ", " << b[2] << "]");
-
- double y[3];
- forwardSubstitution(y, lu, b, idx);
-
- double x[3];
- backwardSubstitution(x, lu, y, idx);
+ void invertLinearTransform();
- // copy to _linearTransform matrix
- // NB : this is a column operation, so we cannot
- // do this directly when we calculate x
- for(int j = 0 ; j < 3 ; j++)
- {
- _linearTransform[3*j + i] = x[idx[j]];
- }
-
- }
-
- }
-
- void calculateDeterminant()
- {
- const double subDet[3] =
- {
- _linearTransform[4] * _linearTransform[8] - _linearTransform[5] * _linearTransform[7],
- _linearTransform[3] * _linearTransform[8] - _linearTransform[5] * _linearTransform[6],
- _linearTransform[3] * _linearTransform[7] - _linearTransform[4] * _linearTransform[6]
- };
-
- _determinant = _linearTransform[0] * subDet[0] - _linearTransform[1] * subDet[1] + _linearTransform[2] * subDet[2];
- }
+ void calculateDeterminant();
/////////////////////////////////////////////////
/// Auxiliary methods for inverse calculation ///
/////////////////////////////////////////////////
- void factorizeLU(double* lu, int* idx/*, double* parity*/) const
- {
- // 3 x 3 LU factorization
- // initialise idx
- for(int i = 0 ; i < 3 ; ++i)
- {
- idx[i] = i;
- }
-
- for(int k = 0; k < 2 ; ++k)
- {
-
- // find pivot
- int i = k;
- double max = fabs(lu[3*idx[k] + k]);
- int row = i;
- while(i < 3)
- {
- if(fabs(lu[3*idx[i] + k]) > max)
- {
- max = fabs(lu[3*idx[i] + k]);
- row = i;
- }
- ++i;
- }
-
- // swap rows
- int tmp = idx[k];
- idx[k] = idx[row];
- idx[row] = tmp;
- /* if(k != row)
- *parity = -(*parity);
- */
-
- // calculate row
- for(int j = k + 1 ; j < 3 ; ++j)
- {
- // l_jk = u_jk / u_kk
- lu[3*idx[j] + k] /= lu[3*idx[k] + k];
- for(int s = k + 1; s < 3 ; ++s)
- {
- // case s = k will always become zero, and then be replaced by
- // the l-value
- // there's no need to calculate this explicitly
-
- // u_js = u_js - l_jk * u_ks
- lu[3*idx[j] + s] -= lu[3*idx[j] + k] * lu[3*idx[k] + s] ;
- }
- }
- }
-
-
- }
+ void factorizeLU(double* lu, int* idx) const;
- void forwardSubstitution(double* x, const double* lu, const double* b, const int* idx) const
- {
- x[idx[0]] = ( b[idx[0]] );// / lu[3*idx[0]];
- x[idx[1]] = ( b[idx[1]] - lu[3*idx[1]] * x[idx[0]] ); // / lu[3*idx[1] + 1];
- x[idx[2]] = ( b[idx[2]] - lu[3*idx[2]] * x[idx[0]] - lu[3*idx[2] + 1] * x[idx[1]] ); // / lu[3*idx[2] + 2];
- }
-
- void backwardSubstitution(double* x, const double* lu, const double* b, const int* idx) const
- {
- x[idx[2]] = ( b[idx[2]] ) / lu[3*idx[2] + 2];
- x[idx[1]] = ( b[idx[1]] - lu[3*idx[1] + 2] * x[idx[2]] ) / lu[3*idx[1] + 1];
- x[idx[0]] = ( b[idx[0]] - lu[3*idx[0] + 1] * x[idx[1]] - lu[3*idx[0] + 2] * x[idx[2]] ) / lu[3*idx[0]];
- }
+ void forwardSubstitution(double* x, const double* lu, const double* b, const int* idx) const;
+ void backwardSubstitution(double* x, const double* lu, const double* b, const int* idx) const;
-
- double _translation[3];
+ /// The affine transformation Ax + b is represented with _linearTransformation containing the elements of
+ /// A in row-first ordering and _translation containing the elements of b
double _linearTransform[9];
+ double _translation[3];
+ /// The determinant of the matrix A is calculated at the construction of the object and cached.
double _determinant;
};
#undef COORDINATE_CORRECTION 1.0e-15
-class CircularSortOrder
-{
-public:
-
- enum CoordType { XY, XZ, YZ };
-
- CircularSortOrder(const double* barycenter, const double* pt0, const CoordType type)
- : _pt0( pt0 )
- {
- // get the indices to use in determinant
- _aIdx = (type == YZ) ? 2 : 0;
- _bIdx = (type == XY) ? 1 : 2;
-
- _a = barycenter[_aIdx] - pt0[_aIdx];
- _b = barycenter[_bIdx] - pt0[_bIdx];
- //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)
- {
- //{ calculations could be optimised here, avoiding intermediary affectations
- const double diff[4] =
- {
- pt1[_aIdx] - _pt0[_aIdx], pt1[_bIdx] - _pt0[_bIdx], // pt1 - pt0
- pt2[_aIdx] - _pt0[_aIdx], pt2[_bIdx] - _pt0[_bIdx], // pt2 - pt0
- };
-
- // We need to check if one of the points is equal to pt0,
- // since pt0 should always end up at the beginning
- // is pt1 == pt0 ? -> if yes, pt1 < pt2
- if(diff[0] == 0.0 && diff[1] == 0.0)
- return true;
- // is pt2 == pt0 ? -> if yes pt2 < pt1
- if(diff[2] == 0.0 && diff[3] == 0.0)
- return false;
-
- // normal case
- const double det1 = _a*diff[1] - _b*diff[0];
- const double det2 = _a*diff[3] - _b*diff[2];
-
- return det1 < det2;
- }
-
-private:
- int _aIdx, _bIdx;
- double _a, _b;
- const double* _pt0;
-};
class ProjectedCentralCircularSortOrder
{
class Vector3Cmp
{
- public:
+public:
bool operator()(double* const& pt1, double* const& pt2)
{
LOG(6, "points are equal ? : " << int((pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2])));
_coords[5*Q + 3] = 1 - q[0] - q[1] - q[2];
_coords[5*R + 3] = 1 - r[0] - r[1] - r[2];
- // order of substractions give different results ?
- /*
- _coords[5*P + 3] = 1 - (p[0] - (p[1] - p[2]));
- _coords[5*Q + 3] = 1 - (q[0] - (q[1] - q[2]));
- 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]));
- */
-
-
// H coordinate
_coords[5*P + 4] = 1 - p[0] - p[1];
_coords[5*Q + 4] = 1 - q[0] - q[1];
// initialise rest of data
preCalculateDoubleProducts();
+
+#ifdef OPTIMIZE
+ preCalculateTriangleSurroundsEdge();
+#endif
+
preCalculateTripleProducts();
}
*/
double TransformedTriangle::calculateIntersectionVolume()
{
- // check first that we are not below z - plane
-
+ // check first that we are not below z - plane
if(isTriangleBelowTetraeder())
{
LOG(2, " --- Triangle is below tetraeder - V = 0.0");
return 0.0;
}
-
+
// get the sign of the volume - equal to the sign of the z-component of the normal
// of the triangle, u_x * v_y - u_y * v_x, where u = q - p and v = r - p
// if it is zero, the triangle is perpendicular to the z - plane and so the volume is zero
double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
-
if(sign == 0.0)
{
LOG(2, " --- Triangle is perpendicular to z-plane - V = 0.0");
// normalize
sign = sign > 0.0 ? 1.0 : -1.0;
-
LOG(2, "-- Calculating intersection polygons ... ");
calculateIntersectionPolygons();
#ifdef MERGE_CALC
if((!mergeCalculation) && _polygonB.size() > 2)
#else
- if(!isTriangleInPlaneOfFacet(XYZ) && _polygonB.size() > 2)
+ if(!isTriangleInPlaneOfFacet(XYZ) && _polygonB.size() > 2)
#endif
- {
- LOG(2, "---- Treating polygon B ... ");
+ {
+ LOG(2, "---- Treating polygon B ... ");
- calculatePolygonBarycenter(B, barycenter);
- sortIntersectionPolygon(B, barycenter);
- volB = calculateVolumeUnderPolygon(B, barycenter);
- LOG(2, "Volume is " << sign * volB);
- }
+ calculatePolygonBarycenter(B, barycenter);
+ sortIntersectionPolygon(B, barycenter);
+ volB = calculateVolumeUnderPolygon(B, barycenter);
+ LOG(2, "Volume is " << sign * volB);
+ }
LOG(2, "volA + volB = " << sign * (volA + volB) << std::endl << "***********");
assert(_polygonA.size() == 0);
assert(_polygonB.size() == 0);
+#ifdef OPTIMIZE
+ // avoid reallocations in push_back() by pre-allocating enough memory
+ // we should never have more than 20 points
+ _polygonA.reserve(20);
+ _polygonB.reserve(20);
+#endif
+
// -- surface intersections
// surface - edge
for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
{
if(testSurfaceEdgeIntersection(edge))
{
- //{ we only really need to calculate the point once
double* ptA = new double[3];
calcIntersectionPtSurfaceEdge(edge, ptA);
_polygonA.push_back(ptA);
LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B");
}
}
+#ifdef OPTIMIZE
+
+ // -- segment intersections
+ for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
+ {
+
+ bool isZero[NO_DP];
+
+ // check beforehand which double-products are zero
+ for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
+ {
+ isZero[dp] = (calcStableC(seg, dp) == 0.0);
+ }
+
+ // segment - facet
+ for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
+ {
+ // is this test worth it?
+ const bool doTest =
+ !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] &&
+ !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
+ !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
+
+ if(doTest && testSegmentFacetIntersection(seg, facet))
+ {
+ double* ptA = new double[3];
+ calcIntersectionPtSegmentFacet(seg, facet, ptA);
+ _polygonA.push_back(ptA);
+ LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
+ if(facet == XYZ)
+ {
+ double* ptB = new double[3];
+ copyVector3(ptA, ptB);
+ _polygonB.push_back(ptB);
+ LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B");
+ }
+
+ }
+ }
+
+ // segment - edge
+ for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
+ {
+ const DoubleProduct edge_dp = DoubleProduct(edge);
+
+ if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge))
+ {
+ double* ptA = new double[3];
+ calcIntersectionPtSegmentEdge(seg, edge, ptA);
+ _polygonA.push_back(ptA);
+ LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
+ if(edge >= XY)
+ {
+ double* ptB = new double[3];
+ copyVector3(ptA, ptB);
+ _polygonB.push_back(ptB);
+ }
+ }
+ }
+
+ // segment - corner
+ for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
+ {
+ const bool doTest =
+ isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
+ isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
+ isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
+
+ if(doTest && testSegmentCornerIntersection(seg, corner))
+ {
+ double* ptA = new double[3];
+ copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
+ _polygonA.push_back(ptA);
+ LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
+ if(corner != O)
+ {
+ double* ptB = new double[3];
+ _polygonB.push_back(ptB);
+ copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
+ LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B");
+ }
+ }
+ }
+#else
+
// -- segment intersections
for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
{
+
// segment - facet
for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
{
_polygonB.push_back(ptB);
LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B");
}
+
}
}
}
}
}
+#endif
+
+#ifdef OPTIMIZE
+
+ // segment - ray
+ for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
+ {
+ if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*corner]] && testSegmentRayIntersection(seg, corner))
+ {
+ double* ptB = new double[3];
+ copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
+ _polygonB.push_back(ptB);
+ LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B");
+ }
+ }
+
+ // segment - halfstrip
+ for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
+ {
+
+#if 0
+ const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX
+ const bool doTest =
+ !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] &&
+ !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]];
+
+
+ if(doTest && testSegmentHalfstripIntersection(seg, edge))
+#endif
+ if(testSegmentHalfstripIntersection(seg, edge))
+ {
+ double* ptB = new double[3];
+ calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
+ _polygonB.push_back(ptB);
+ LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
+ }
+ }
+
+
+#else
// segment - ray
for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
}
}
+
+#endif
}
// inclusion tests
#ifdef EPS_TESTING
if(!epsilonEqualRelative(_coords[5*c + coord], 0.0, TEST_EPS * _coords[5*c + coord]))
#else
- if(_coords[5*c + coord] != 0.0)
+ if(_coords[5*c + coord] != 0.0)
#endif
- {
- return false;
- }
+ {
+ return false;
+ }
}
return true;
return true;
}
-void TransformedTriangle::dumpCoords()
-{
- std::cout << "Coords : ";
- for(int i = 0 ; i < 3; ++i)
- {
- std::cout << vToStr(&_coords[5*i]) << ",";
- }
- std::cout << std::endl;
-}
+ void TransformedTriangle::dumpCoords()
+ {
+ std::cout << "Coords : ";
+ for(int i = 0 ; i < 3; ++i)
+ {
+ std::cout << vToStr(&_coords[5*i]) << ",";
+ }
+ std::cout << std::endl;
+ }
}; // NAMESPACE
#define TEST_EPS 1.0e-14
#endif
+#ifdef OPTIMIZE
+#define OPT_INLINE inline
+#else
+#define OPT_INLINE
+#endif
+
// 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;
-class TransformedTriangleCalcVolumeTest;
-
namespace INTERP_UTILS
{
*
* Reference : J. Grandy, "Conservative Remapping and Region Overlays by Intersecting Arbitrary Polyhedra",
* Journal of Computational Physics (1999)
+ */
+
+ /* ! READ ME FIRST !
+ * OVERVIEW of how the class works : (details can be found in the commentaries above each method)
+ *
+ * Constructor :
+ * The constructor takes as arguments three pointers to double[3] vectors holding the transformed
+ * coordinates of the corners of the triangle. It copies there coordinates and then proceeds to pre-calculating certain
+ * entities used in the intersection calculation : the double products, triple products and the values of the function E
+ * (Grandy, [53]).
+ *
+ * calculateIntersectionVolume() :
+ * This is the only method in the public interface. It calculates the volume under the intersection polygons
+ * between the triangle and the unit tetrahedron, as described in Grandy, pp. 435-447. It does this by first calculating the
+ * intersection polygons A and B, with the method calculateIntersectionPolygons(). It then calculates the barycenter of each
+ * polygon in calculatePolygonBarycenter(), and sorts their points in a circular order around the barycenter in
+ * sortIntersecionPolygon(). The sorting is done with STL sort, using the order defined in the class
+ * ProjectedCentralCircularSortOrder. The volume under each polygon is then calculated with calculateVolumeUnderPolygon(), which
+ * implements formula [34] in Grandy.
+ *
+ * calculateIntersectionPolygons() :
+ * This method goes through all the possible ways in which the triangle can intersect the tetrahedron and tests for these
+ * types of intersections in accordance with the formulas described in Grandy. These tests are implemented in the test* - methods.
+ * The formulas in the article are stated for one case each only, while the calculation must take into account all cases.
+ * To this end, a number of tables, implemented as static const arrays of different types, are used. The tables
+ * mainly contain values of the different enumeration types described at the beginning of the class interface. For example,
+ * the formula Grandy gives for the segment-halfstrip intersection tests ([30]) is for use with the halfstrip above the zx edge.
+ * For the other two halfstrips (above the xy and yz edges), other double products are used.
+ * These double products are stored in the table DP_FOR_HALFSTRIP_INTERSECTION, which permits to treat
+ * all the edges equally, avoiding switch() - statements. It is the careful choice of order of the enumeration types that makes this
+ * possible. Notably, there is a correspondance between the TetraEdge type and the DoubleProduct type (see Grandy, tatble III) that
+ * is used throughout the code, permitting statements such as DoubleProduct(my_edge) to work.
+ * When an intersection point has been detected it is calculated with a corresponding calc* - method in the cases where it
+ * is not known directly. It is then added to the polygon A and/or B as necessary.
+ *
*
*/
class TransformedTriangle
friend class ::TransformedTriangleTest;
friend class ::TransformedTriangleIntersectTest;
- friend class ::TransformedTriangleCalcVolumeTest;
/**
double calculateIntersectionVolume();
- // temporary debug method
void dumpCoords();
/// Intersection test methods and intersection point calculations ////////
////////////////////////////////////////////////////////////////////////////////////
- bool testSurfaceEdgeIntersection(const TetraEdge edge) const;
+ OPT_INLINE bool testSurfaceEdgeIntersection(const TetraEdge edge) const;
void calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const;
- bool testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const;
+ OPT_INLINE bool testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const;
void calcIntersectionPtSegmentFacet(const TriSegment seg, const TetraFacet facet, double* pt) const;
bool testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const ;
- bool testSurfaceRayIntersection(const TetraCorner corner) const;
+ OPT_INLINE bool testSurfaceRayIntersection(const TetraCorner corner) const;
bool testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edg);
bool testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const;
- bool testCornerInTetrahedron(const TriCorner corner) const;
+ OPT_INLINE bool testCornerInTetrahedron(const TriCorner corner) const;
- bool testCornerOnXYZFacet(const TriCorner corner) const;
+ OPT_INLINE bool testCornerOnXYZFacet(const TriCorner corner) const;
- bool testCornerAboveXYZFacet(const TriCorner corner) const;
+ OPT_INLINE bool testCornerAboveXYZFacet(const TriCorner corner) const;
////////////////////////////////////////////////////////////////////////////////////
/// Utility methods used in intersection tests ///////////////
bool areDoubleProductsConsistent(const TriSegment seg) const;
- void resetDoubleProducts(const TriSegment seg, const TetraCorner corner);
+ OPT_INLINE void resetDoubleProducts(const TriSegment seg, const TetraCorner corner);
double calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const;
double calculateAngleEdgeTriangle(const TetraEdge edge) const;
- double calcStableC(const TriSegment seg, const DoubleProduct dp) const;
+ OPT_INLINE double calcStableC(const TriSegment seg, const DoubleProduct dp) const;
- double calcStableT(const TetraCorner corner) const;
+ OPT_INLINE double calcStableT(const TetraCorner corner) const;
- double calcUnstableC(const TriSegment seg, const DoubleProduct dp) const;
+ OPT_INLINE double calcUnstableC(const TriSegment seg, const DoubleProduct dp) const;
double calcTByDevelopingRow(const TetraCorner corner, const int row = 1, const bool project = false) const;
////////////////////////////////////////////////////////////////////////////////////
private:
- // storage :
+ // order :
// [ p_x, p_y, p_z, p_h, p_H, q_x, q_y, q_z, q_h, q_H, r_x, r_y, r_z, r_h, r_H ]
double _coords[15];
+
+ /// flags showing whether the double and triple products have been precalculated for this class
+ bool _isDoubleProductsCalculated, _isTripleProductsCalculated;
- bool _isDoubleProductsCalculated, _isTripleProductsCalculated;
+ /// array containing the 24 double products
+ /// order : c^PQ_YZ, ... ,cPQ_10, ... c^QR_YZ, ... c^RP_YZ
+ /// following order in enumeration DoubleProduct
double _doubleProducts[24];
+
+ /// array containing the 4 triple products
+ /// order : t_O, t_X, t_Y, t_Z
double _tripleProducts[4];
+
+ /// arrays holding the points in the two intersection polygons A and B
+ /// these points are allocated in calculateIntersectionPolygons() and liberated in the destructor
std::vector<double*> _polygonA, _polygonB;
+
+ /// vectors holding the coordinates of the barycenters of the polygons A and B
+ /// these points are calculated in calculatePolygonBarycenter
double _barycenterA[3], _barycenterB[3];
// used for debugging
bool _validTP[4];
+#ifdef OPTIMIZE
+ void preCalculateTriangleSurroundsEdge();
+ bool _triangleSurroundsEdgeCache[NO_TET_EDGE];
+ bool _isOutsideTetra;
+#endif
+
////////////////////////////////////////////////////////////////////////////////////
/// Constants /////////////////
// correspondance edge - corners
static const TetraCorner CORNERS_FOR_EDGE[12];
+ // correspondance edge - facets
+ // facets shared by each edge
+ static const TetraFacet FACET_FOR_EDGE[12];
+
+ // correspondance edge - corners
+ static const TetraEdge EDGES_FOR_CORNER[12];
+
+ // double products used in segment-halfstrip test
+ static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12];
+
+ // double products used in segment - ray test
+ static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21];
+
+
+ inline bool isTriangleOutsideTetra(void) const;
+
};
+
+
+ // include definitions of inline methods
+#ifdef OPTIMIZE
+
+#include "TransformedTriangle_inline.hxx"
+
+#endif
+
};
#endif
--- /dev/null
+inline void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const TetraCorner corner)
+{
+ // set the three corresponding double products to 0.0
+ static const DoubleProduct DOUBLE_PRODUCTS[12] =
+ {
+ C_YZ, C_ZX, C_XY, // O
+ C_YZ, C_ZH, C_YH, // X
+ C_ZX, C_ZH, C_XH, // Y
+ C_XY, C_YH, C_XH // Z
+ };
+
+ for(int i = 0 ; i < 3 ; ++i) {
+ const DoubleProduct dp = DOUBLE_PRODUCTS[3*corner + i];
+
+ LOG(6, std::endl << "resetting inconsistent dp :" << dp << " for corner " << corner);
+ _doubleProducts[8*seg + dp] = 0.0;
+ };
+}
+
+/**
+ * Returns the stable double product c_{xy}^{pq}
+ *
+ * @pre The stable double products have been calculated with preCalculateDoubleProducts.
+ * @param seg segment of triangle
+ * @param dp double product sought
+ *
+ * @returns stabilised double product c_{xy}^{pq}
+ *
+ */
+inline double TransformedTriangle::calcStableC(const TriSegment seg, const DoubleProduct dp) const
+{
+ return _doubleProducts[8*seg + dp];
+}
+
+/**
+ * Returns the stable triple product t_X for a given corner
+ * The triple product gives the signed volume of the tetrahedron between
+ * this corner and the triangle PQR. These triple products have been calculated
+ * in a way to avoid problems with cancellation.
+ *
+ * @pre double products have already been calculated
+ * @pre triple products have already been calculated
+ * @param corner corner for which the triple product is calculated
+ * @returns triple product associated with corner (see Grandy, eqs. [50]-[52])
+ */
+inline double TransformedTriangle::calcStableT(const TetraCorner corner) const
+{
+ assert(_isTripleProductsCalculated);
+ assert(_validTP[corner]);
+ return _tripleProducts[corner];
+}
+
+/**
+ * Calculates the given double product c_{xy}^{pq} = x_p*y_q - y_p*x_q for a
+ * a segment PQ of the triangle. This method does not compensate for
+ * precision errors.
+ *
+ * @param seg segment of triangle
+ * @param dp double product sought
+ *
+ * @returns double product c_{xy}^{pq}
+ *
+ */
+inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp) const
+{
+
+ // find the points of the triangle
+ // 0 -> P, 1 -> Q, 2 -> R
+ const int pt1 = seg;
+ const int pt2 = (seg + 1) % 3;
+
+ // find offsets
+ const int off1 = DP_OFFSET_1[dp];
+ const int off2 = DP_OFFSET_2[dp];
+
+ return _coords[5*pt1 + off1] * _coords[5*pt2 + off2] - _coords[5*pt1 + off2] * _coords[5*pt2 + off1];
+}
+
+inline void TransformedTriangle::preCalculateTriangleSurroundsEdge()
+{
+ for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
+ {
+ _triangleSurroundsEdgeCache[edge] = testTriangleSurroundsEdge(edge);
+ }
+}
+
+/**
+ * Tests if the given edge of the tetrahedron intersects the triangle PQR. (Grandy, eq [17])
+ *
+ * @param edge edge of tetrahedron
+ * @returns true if PQR intersects the edge, and the edge is not in the plane of the triangle.
+ */
+inline bool TransformedTriangle::testSurfaceEdgeIntersection(const TetraEdge edge) const
+{
+ return _triangleSurroundsEdgeCache[edge] && testEdgeIntersectsTriangle(edge);
+}
+
+/**
+ * Tests if the given segment of the triangle intersects the given facet of the tetrahedron.
+ * (Grandy, eq. [19])
+ *
+ * @param seg segment of the triangle
+ * @param facet facet of the tetrahedron
+ * @returns true if the segment intersects the facet
+ */
+inline bool TransformedTriangle::testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const
+{
+ return testFacetSurroundsSegment(seg, facet) && testSegmentIntersectsFacet(seg, facet);
+}
+
+/**
+ * Tests if the triangle PQR intersects the ray pointing in the upwards z - direction
+ * from the given corner of the tetrahedron. (Grandy eq. [29])
+ *
+ * @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z)
+ * @returns true if the upwards ray from the corner intersects the triangle.
+ */
+inline bool TransformedTriangle::testSurfaceRayIntersection(const TetraCorner corner) const
+{
+ return testTriangleSurroundsRay( corner ) && testSurfaceAboveCorner( corner );
+}
+
+/**
+ * Tests if the given corner of triangle PQR lies in the interior of the unit tetrahedron
+ * (0 <= x_p, y_p, z_p, h_p <= 1)
+ *
+ * @param corner corner of the triangle PQR
+ * @returns true if the corner lies in the interior of the unit tetrahedron.
+ */
+inline bool TransformedTriangle::testCornerInTetrahedron(const TriCorner corner) const
+{
+ const double pt[4] =
+ {
+ _coords[5*corner], // x
+ _coords[5*corner + 1], // y
+ _coords[5*corner + 2], // z
+ _coords[5*corner + 3] // z
+ };
+
+ for(int i = 0 ; i < 4 ; ++i)
+ {
+ if(pt[i] < 0.0 || pt[i] > 1.0)
+ {
+ return false;
+ }
+ }
+ return true;
+}
+
+/**
+ * Tests if the given corner of triangle PQR lies on the facet h = 0 (the XYZ facet)
+ * (0 <= x_p, y_p, z_p <= 1 && h_p = 0)
+ *
+ * @param corner corner of the triangle PQR
+ * @returns true if the corner lies on the facet h = 0
+ */
+inline bool TransformedTriangle::testCornerOnXYZFacet(const TriCorner corner) const
+ {
+ const double pt[4] =
+ {
+ _coords[5*corner], // x
+ _coords[5*corner + 1], // y
+ _coords[5*corner + 2], // z
+ _coords[5*corner + 3] // h
+ };
+
+ if(pt[3] != 0.0)
+ {
+ return false;
+ }
+
+ for(int i = 0 ; i < 3 ; ++i)
+ {
+ if(pt[i] < 0.0 || pt[i] > 1.0)
+ {
+ return false;
+ }
+ }
+ return true;
+ }
+
+inline bool TransformedTriangle::testCornerAboveXYZFacet(const TriCorner corner) const
+ {
+ const double x = _coords[5*corner];
+ const double y = _coords[5*corner + 1];
+ const double h = _coords[5*corner + 3];
+ const double H = _coords[5*corner + 4];
+
+ return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0;
+
+ }
Y, Z, // YZ
Z, X // ZX
};
+
+ // correspondance edge - facets
+ // facets shared by each edge
+ const TransformedTriangle::TetraFacet TransformedTriangle::FACET_FOR_EDGE[12] =
+ {
+ OXY, OZX, // OX
+ OXY, OYZ, // OY
+ OZX, OYZ, // OZ
+ OXY, XYZ, // XY
+ OYZ, XYZ, // YZ
+ OZX, XYZ // ZX
+ };
+
+ // edges meeting at a given corner
+ const TransformedTriangle::TetraEdge TransformedTriangle::EDGES_FOR_CORNER[12] =
+ {
+ OX, OY, OZ, // O
+ OX, XY, ZX, // X
+ OY, XY, YZ, // Y
+ OZ, ZX, YZ // Z
+ };
+
+ // NB : some uncertainty whether these last are correct
+ const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_HALFSTRIP_INTERSECTION[12] =
+ {
+ C_10, C_01, C_ZH, C_10, // XY
+ C_01, C_XY, C_XH, C_01, // YZ
+ C_XY, C_10, C_YH, C_XY // ZX
+ };
+
+ // double products to use in segment-ray test
+ // dp 1 -> cond 1
+ // dp 2-7 -> cond 3
+#if SEG_RAY_TABLE==1
+ const TransformedTriangle::DoubleProduct TransformedTriangle::DP_SEGMENT_RAY_INTERSECTION[21] =
+ {
+ C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
+ C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_10, // Y
+ C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_01 // Z
+ };
+#else
+ const TransformedTriangle::DoubleProduct TransformedTriangle::DP_SEGMENT_RAY_INTERSECTION[21] =
+ {
+ C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
+ C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_YZ, // Y
+ C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_ZX // Z
+ };
+#endif
+
+
+
////////////////////////////////////////////////////////////////////////////////////
/// Intersection test methods and intersection point calculations ////////
////////////////////////////////////////////////////////////////////////////////////
+#ifndef OPTIMIZE
/**
* Tests if the given edge of the tetrahedron intersects the triangle PQR. (Grandy, eq [17])
*
{
return testTriangleSurroundsEdge(edge) && testEdgeIntersectsTriangle(edge);
}
-
+#endif
/**
* Calculates the point of intersection between the given edge of the tetrahedron and the
* triangle PQR. (Grandy, eq [22])
}
}
+#ifndef OPTIMIZE
/**
* Tests if the given segment of the triangle intersects the given facet of the tetrahedron.
* (Grandy, eq. [19])
{
return testFacetSurroundsSegment(seg, facet) && testSegmentIntersectsFacet(seg, facet);
}
+#endif
/**
* Calculates the point of intersection between the given segment of the triangle
*/
bool TransformedTriangle::testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const
{
- // facets shared by each edge
- static const TetraFacet FACET_FOR_EDGE[12] =
- {
- OXY, OZX, // OX
- OXY, OYZ, // OY
- OZX, OYZ, // OZ
- OXY, XYZ, // XY
- OYZ, XYZ, // YZ
- OZX, XYZ // ZX
- };
+
+#ifndef OPTIMIZE // in this case, we have already checked if the double product is zero
#ifdef EPS_TESTING
if(!epsilonEqual(calcStableC(seg,DoubleProduct( edge )), 0.0, TEST_EPS))
{
return false;
}
- else
+ else
+#endif // OPTIMIZE
{
// check condition that the double products for one of the two
// facets adjacent to the edge has a positive product
*/
bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const
{
- // edges meeting at a given corner
- static const TetraEdge EDGES_FOR_CORNER[12] =
- {
- OX, OY, OZ, // O
- OX, XY, ZX, // X
- OY, XY, YZ, // Y
- OZ, ZX, YZ // Z
- };
+
// facets meeting at a given corner
static const TetraFacet FACETS_FOR_CORNER[12] =
OZX, XYZ, OYZ // Z
};
+#ifndef OPTIMIZE // if optimized, we have already checked that the double products are zero
// check double products are zero
for(int i = 0 ; i < 3 ; ++i)
{
return false;
}
}
+#endif
// check segment intersect a facet
for(int i = 0 ; i < 3 ; ++i)
return false;
}
-
+#ifndef OPTIMIZE
/**
* Tests if the triangle PQR intersects the ray pointing in the upwards z - direction
* from the given corner of the tetrahedron. (Grandy eq. [29])
{
return testTriangleSurroundsRay( corner ) && testSurfaceAboveCorner( corner );
}
+#endif
/**
* Tests if the given segment of the triangle intersects the half-strip above the
C_01, C_XY, C_XH, C_01, // YZ
C_XY, C_10, C_YH, C_XY // ZX
};
- /*static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12] =
- {
- C_10, C_01, C_ZH, C_10, // XY
- C_01, C_XY, C_XH, C_ZX, // YZ
- C_XY, C_10, C_YH, C_XY // ZX
- };
- */
- /*static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12] =
- {
- C_10, C_01, C_ZH, C_XY, // XY
- C_01, C_XY, C_XH, C_XY, // YZ
- C_XY, C_10, C_YH, C_XY // ZX
- };
- */
-
+
// facets to use in second condition (S_m)
- static TetraFacet FACET_FOR_HALFSTRIP_INTERSECTION[3] =
+ static const TetraFacet FACET_FOR_HALFSTRIP_INTERSECTION[3] =
{
NO_TET_FACET, // XY -> special case : test with plane H = 0
OYZ, // YZ
// special case : facet H = 0
- bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet);
+ const bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet);
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] );
// for edge AB : (x,y,z)* = (1-alpha) * A + alpha * B
// where alpha = cB / (cB - cA)
- const DoubleProduct DP_FOR_EDGE[6] =
- {
- C_10, C_01, // XY
- C_01, C_XY, // YZ
- C_XY, C_10 // ZX
- };
-
- const double cA = calcStableC(seg, DP_FOR_EDGE[2*edgeIndex]);
- const double cB = calcStableC(seg, DP_FOR_EDGE[2*edgeIndex + 1]);
+ const double cA = calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex]);
+ const double cB = calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 1]);
assert(cA != cB);
const double alpha = cA / (cA - cB);
// readjust index since O is not used
const int cornerIdx = static_cast<int>(corner) - 1;
- // double products to use in test
- // dp 1 -> cond 1
- // dp 2-7 -> cond 3
- //? NB : last two rows are not completely understood and may contain errors
-#if SEG_RAY_TABLE==1
- static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21] =
- {
- C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
- C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_10, // Y
- C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_01 // Z
- };
-#else
-
- static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21] =
- {
- C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
- C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_YZ, // Y
- C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_ZX // Z
- };
-#endif
-
// facets to use
//? not sure this is correct
static const TetraFacet FIRST_FACET_SEGMENT_RAY_INTERSECTION[3] =
const DoubleProduct dp0 = DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx];
const double cVal0 = calcStableC(seg, dp0);
-
+
+#ifndef OPTIMIZE // in this case we have already checked that the double product is zero
//? epsilon-equality here?
// cond. 1
#ifdef EPS_TESTING
return false;
}
+#endif
+
// cond 2
const bool cond21 = testSegmentIntersectsFacet(seg, FIRST_FACET_SEGMENT_RAY_INTERSECTION[cornerIdx]);
const bool cond22 = (corner == Z) ? testSegmentIntersectsFacet(seg, OYZ) : testSegmentIntersectsHPlane(seg);
}
+#ifndef OPTIMIZE
/**
* Tests if the given corner of triangle PQR lies in the interior of the unit tetrahedron
* (0 <= x_p, y_p, z_p, h_p <= 1)
return true;
}
+
/**
* Tests if the given corner of triangle PQR lies on the facet h = 0 (the XYZ facet)
* (0 <= x_p, y_p, z_p <= 1 && h_p = 0)
return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0;
}
-
+#endif
////////////////////////////////////////////////////////////////////////////////////
_doubleProducts[8*seg + dp] = calcUnstableC(seg, dp);
}
}
-
+
+ std::map<double, TetraCorner> distances;
// -- (1) for each segment : check that double products satisfy Grandy, [46]
// -- and make corrections if not
{
if(!areDoubleProductsConsistent(seg))
{
- std::map<double, TetraCorner> distances;
+
LOG(4, "inconsistent! ");
for(TetraCorner corner = O ; corner <= Z ; corner = TetraCorner(corner + 1))
}
#endif
}
+ distances.clear();
}
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 );
+ LOG(2, "for seg " << seg << " consistency " << term1 + term2 + term3 );
+ LOG(2, "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);
// return epsilonEqual(term1 + term2 + term3, 0.0);
}
+#ifndef OPTIMIZE
void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const TetraCorner corner)
{
// set the three corresponding double products to 0.0
_doubleProducts[8*seg + dp] = 0.0;
};
}
-
+#endif OPTIMIZE
double TransformedTriangle::calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const
{
// NB uses fact that TriSegment <=> TriCorner that is first point of segment (PQ <=> P)
void TransformedTriangle::preCalculateTripleProducts(void)
{
if(_isTripleProductsCalculated)
- return;
+ {
+ return;
+ }
+
+ // find edge / row to use -> that whose edge makes the smallest angle to the triangle
+ // use a map to find the minimum
+ std::map<double, int> anglesForRows;
LOG(4, "Precalculating triple products" );
for(TetraCorner corner = O ; corner <= Z ; corner = TetraCorner(corner + 1))
{
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
- std::map<double, int> anglesForRows;
for(int row = 1 ; row < 4 ; ++row)
{
TetraEdge edge = TetraEdge(dp);
// use edge only if it is surrounded by the surface
+#ifdef OPTIMIZE
+ if( _triangleSurroundsEdgeCache[edge] )
+#else
if( testTriangleSurroundsEdge(edge) )
+#endif
{
// -- calculate angle between edge and PQR
const double angle = calculateAngleEdgeTriangle(edge);
_validTP[corner] = false;
}
+ anglesForRows.clear();
}
}
+#ifndef OPTIMIZE
/**
* Returns the stable double product c_{xy}^{pq}
*
return _doubleProducts[8*seg + dp];
}
-
/**
* Returns the stable triple product t_X for a given corner
* The triple product gives the signed volume of the tetrahedron between
return _tripleProducts[corner];
}
+
/**
* Calculates the given double product c_{xy}^{pq} = x_p*y_q - y_p*x_q for a
* a segment PQ of the triangle. This method does not compensate for
return _coords[5*pt1 + off1] * _coords[5*pt2 + off2] - _coords[5*pt1 + off2] * _coords[5*pt2 + off1];
}
+#endif
+
/**
* Calculates triple product associated with the given corner of tetrahedron, developing
* the determinant by the given row. The triple product gives the signed volume of
#define VOL_PREC 1.0e-6
#define DEFAULT_REL_TOL 1.0e-6
-#define DEFAULT_ABS_TOL 1.0e-12
+#define DEFAULT_ABS_TOL 5.0e-12
namespace INTERP_UTILS
{
-
+ ///////////////////////////////////////////////////////////////////////
+ // Math operations for vectors represented by double[3] - arrays ////
+ ///////////////////////////////////////////////////////////////////////
+
+ /*
+ * Copies a double[3] vector from src to dest
+ *
+ * @param src source vector
+ * @param dest destination vector
+ *
+ */
inline void copyVector3(const double* src, double* dest)
{
for(int i = 0 ; i < 3 ; ++i)
}
}
+ /*
+ * Creates a string representation of a double[3] vector
+ *
+ * @param pt a 3-vector
+ * @return a string of the form [x, y, z]
+ */
inline const std::string vToStr(const double* pt)
{
std::stringstream ss(std::ios::out);
return ss.str();
}
+ /*
+ * Calculates the cross product of two double[3] - vectors.
+ *
+ * @param v1 vector v1
+ * @param v2 vector v2
+ * @param res vector in which to store the result v1 x v2. It should not be one of v1 and v2.
+ */
inline void cross(const double* v1, const double* v2,double* res)
{
res[0] = v1[1]*v2[2] - v1[2]*v2[1];
res[2] = v1[0]*v2[1] - v1[1]*v2[0];
}
+ /*
+ * Calculates the dot product of two double[3] - vectors
+ *
+ * @param v1 vector v1
+ * @param v2 vector v2
+ * @return dot (scalar) product v1.v2
+ */
inline double dot(const double* v1, const double* v2)
{
return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
}
+ /*
+ * Calculates norm of a double[3] vector
+ *
+ * @param v a vector v
+ * @return euclidean norm of v
+ */
inline double norm(const double* v)
{
return sqrt(dot(v,v));
}
- inline double angleBetweenVectors(const double* v1, const double* v2, const double* n)
- {
- const double denominator = dot(v1, v2);
- double v3[3];
-
- cross(v1, v2, v3);
- const double numerator = dot(n, v3);
- return atan2(numerator, denominator);
-
- }
-
- /// Should be used for comparisons to zero
+ /*
+ * Compares doubles using an absolute tolerance
+ * This is suitable mainly for comparisons with 0.0
+ *
+ * @param x first value
+ * @param y second value
+ * @param errTol maximum allowed absolute difference that is to be treated as equality
+ * @returns true if |x - y| < errTol, false otherwise
+ */
inline bool epsilonEqual(const double x, const double y, const double errTol = DEFAULT_ABS_TOL)
{
- return std::abs(x - y) <= errTol;
+ return y < x ? x - y < errTol : y - x < errTol;
+ // return std::abs(x - y) < errTol;
}
- // Should be used for comparisons between numbers that could be large
+ /*
+ * Compares doubles using a relative tolerance
+ * This is suitable mainly for comparing larger values to each other. Before performing the relative test,
+ * an absolute test is performed to guard from problems when comparing to 0.0
+ *
+ * @param x first value
+ * @param y second value
+ * @param relTol maximum allowed relative difference that is to be treated as equality
+ * @param absTol maximum allowed absolute difference that is to be treated as equality
+ * @returns true if |x - y| <= absTol or |x - y|/max(|x|,|y|) <= relTol, false otherwise
+ */
inline bool epsilonEqualRelative(const double x, const double y, const double relTol = DEFAULT_REL_TOL, const double absTol = DEFAULT_ABS_TOL)
{
// necessary for comparing values close to zero
return true;
}
- const double relError = std::abs((x - y) / std::max(x, y));
+ const double relError = std::abs((x - y) / std::max(std::abs(x), std::abs(y)));
- return relError <= relTol;
+ return relError < relTol;
}
};
