From: vbd Date: Wed, 29 Aug 2007 12:53:46 +0000 (+0000) Subject: staffan : commit for ParaMEDMEM integration X-Git-Tag: trio_trio_coupling~62 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=5c543cd3746ef65394e1b44ed80c0c5d094b2a19;p=tools%2Fmedcoupling.git staffan : commit for ParaMEDMEM integration --- diff --git a/src/INTERP_KERNEL/BoundingBox.cxx b/src/INTERP_KERNEL/BoundingBox.cxx index 7bf4b12b0..2a4e8d921 100644 --- a/src/INTERP_KERNEL/BoundingBox.cxx +++ b/src/INTERP_KERNEL/BoundingBox.cxx @@ -2,15 +2,10 @@ #include #include +#include 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 @@ -22,9 +17,12 @@ namespace INTERP_UTILS * */ 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]; @@ -44,12 +42,17 @@ namespace INTERP_UTILS } /** - * 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]); @@ -64,26 +67,26 @@ namespace INTERP_UTILS */ 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) @@ -145,6 +148,10 @@ namespace INTERP_UTILS } } + /* + * Prints the coordinates of the box to std::cout + * + */ void BoundingBox::dumpCoords() const { std::cout << "[xmin, xmax] = [" << _coords[XMIN] << ", " << _coords[XMAX] << "]" << " | "; @@ -153,15 +160,21 @@ namespace INTERP_UTILS 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; } } diff --git a/src/INTERP_KERNEL/BoundingBox.hxx b/src/INTERP_KERNEL/BoundingBox.hxx index 03bc074b4..4abe00d0a 100644 --- a/src/INTERP_KERNEL/BoundingBox.hxx +++ b/src/INTERP_KERNEL/BoundingBox.hxx @@ -12,61 +12,17 @@ namespace INTERP_UTILS { 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); @@ -77,14 +33,12 @@ namespace INTERP_UTILS 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 diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx index b3ef55ad7..318d58c29 100644 --- a/src/INTERP_KERNEL/Interpolation3D.cxx +++ b/src/INTERP_KERNEL/Interpolation3D.cxx @@ -31,7 +31,6 @@ namespace MEDMEM */ Interpolation3D::Interpolation3D() { - // not implemented } /** @@ -40,7 +39,6 @@ namespace MEDMEM */ Interpolation3D::~Interpolation3D() { - // not implemented } /** @@ -57,69 +55,81 @@ namespace MEDMEM //} 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 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(); @@ -141,7 +151,7 @@ namespace MEDMEM } // 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. @@ -157,9 +167,9 @@ namespace MEDMEM 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 @@ -171,14 +181,17 @@ namespace MEDMEM for(vector::const_iterator iter = currNode->getTargetRegion().getBeginElements() ; iter != currNode->getTargetRegion().getEndElements() ; ++iter) { + const int targetIdx = (*iter)->getIndex(); const double vol = intersector->intersectCells(srcIdx, targetIdx); - if(!epsilonEqual(vol, 0.0)) + //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 { @@ -200,10 +213,10 @@ namespace MEDMEM // 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(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; @@ -229,6 +242,7 @@ namespace MEDMEM 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); @@ -297,6 +311,7 @@ namespace MEDMEM for(vector::const_iterator iter = intersectElems.begin() ; iter != intersectElems.end() ; ++iter) { +#if 0 const int srcIdx = *iter + 1; const double vol = intersector->intersectCells(srcIdx, targetIdx); @@ -304,6 +319,7 @@ namespace MEDMEM { matrix[srcIdx - 1].insert(make_pair(targetIdx, vol)); } +#endif } } @@ -323,92 +339,4 @@ namespace MEDMEM } - /** - * 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 triangles; - srcElement.triangulate(triangles, T); - - double volume = 0.0; - - // std::cout << "num triangles = " << triangles.size() << std::endl; - int i = 0; - for(vector::iterator iter = triangles.begin() ; iter != triangles.end(); ++iter) - { - std::cout << std::endl << "= > Triangle " << ++i << std::endl; - iter->dumpCoords(); - volume += iter->calculateIntersectionVolume(); - } - - std::cout << "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 } diff --git a/src/INTERP_KERNEL/Interpolation3D.hxx b/src/INTERP_KERNEL/Interpolation3D.hxx index 3d2d21705..8bdf51b40 100644 --- a/src/INTERP_KERNEL/Interpolation3D.hxx +++ b/src/INTERP_KERNEL/Interpolation3D.hxx @@ -9,7 +9,6 @@ #include #include -// NOTE : Method comments should be moved to source file during the implementation phase // typedefs typedef std::vector< std::map< int, double > > IntersectionMatrix; @@ -36,73 +35,15 @@ namespace MEDMEM 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); - }; }; diff --git a/src/INTERP_KERNEL/Intersector3D.cxx b/src/INTERP_KERNEL/Intersector3D.cxx index 37fbe0aa0..a41ed2ac7 100644 --- a/src/INTERP_KERNEL/Intersector3D.cxx +++ b/src/INTERP_KERNEL/Intersector3D.cxx @@ -1,25 +1,61 @@ #include "Intersector3D.hxx" - +#include "TetraAffineTransform.hxx" +#include "TransformedTriangle.hxx" +#include "MeshUtils.hxx" +#include "VectorUtils.hxx" #include "Log.hxx" +#include +#include + +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]; @@ -28,7 +64,7 @@ namespace MEDMEM tetraCorners[i] = getCoordsOfNode(i + 1, targetCell, _targetMesh); } - // create AffineTransform + // create AffineTransform from tetrahedron TetraAffineTransform T( tetraCorners ); // check if we have planar tetra element @@ -38,56 +74,156 @@ namespace MEDMEM 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 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::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& triangles, const TetraAffineTransform& T) const +#ifdef OPTIMIZE_ // optimized version + void Intersector3D::triangulate(const medGeometryElement type, const int element, std::vector& 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& 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) @@ -96,10 +232,8 @@ namespace MEDMEM 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) @@ -110,20 +244,17 @@ namespace MEDMEM 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) { @@ -132,17 +263,36 @@ namespace MEDMEM 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 - - + }; diff --git a/src/INTERP_KERNEL/Intersector3D.hxx b/src/INTERP_KERNEL/Intersector3D.hxx index 2d1375bcd..8642d2cf8 100644 --- a/src/INTERP_KERNEL/Intersector3D.hxx +++ b/src/INTERP_KERNEL/Intersector3D.hxx @@ -4,37 +4,65 @@ #include "MEDMEM_define.hxx" #include "MEDMEM_Mesh.hxx" -#include "TetraAffineTransform.hxx" -#include "TransformedTriangle.hxx" -#include "MeshUtils.hxx" #include "Intersector.hxx" #include -using namespace MEDMEM; -using namespace MED_EN; -using namespace INTERP_UTILS; +namespace 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& 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& 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; + } + }; diff --git a/src/INTERP_KERNEL/Log.hxx b/src/INTERP_KERNEL/Log.hxx index c3adaf883..86591b9be 100644 --- a/src/INTERP_KERNEL/Log.hxx +++ b/src/INTERP_KERNEL/Log.hxx @@ -2,23 +2,25 @@ #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 -#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 diff --git a/src/INTERP_KERNEL/Makefile.in b/src/INTERP_KERNEL/Makefile.in index 1ab573047..dab294301 100644 --- a/src/INTERP_KERNEL/Makefile.in +++ b/src/INTERP_KERNEL/Makefile.in @@ -88,20 +88,38 @@ LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome CXXFLAGS+=@CXXTMPDPTHFLAGS@ $(MED2_INCLUDES) CPPFLAGS+=$(BOOST_CPPFLAGS) -CXXFLAGS+=-fprofile-arcs -ftest-coverage -CPPFLAGS+=-fprofile-arcs -ftest-coverage +# 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} @@ -109,6 +127,8 @@ ifeq ($(MED_WITH_KERNEL),yes) LDFLAGSFORBIN+= ${KERNEL_LDFLAGS} -lSALOMELocalTrace -lSALOMEBasics endif +LDFLAGSFORBIN+= -pg + LIBSFORBIN=$(BOOSTLIBS) $(MPI_LIBS) LIBS= diff --git a/src/INTERP_KERNEL/MeshElement.cxx b/src/INTERP_KERNEL/MeshElement.cxx index d13585ad3..3069662c7 100644 --- a/src/INTERP_KERNEL/MeshElement.cxx +++ b/src/INTERP_KERNEL/MeshElement.cxx @@ -7,23 +7,24 @@ #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) @@ -48,61 +49,67 @@ namespace INTERP_UTILS } } - /** - * 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; } diff --git a/src/INTERP_KERNEL/MeshElement.hxx b/src/INTERP_KERNEL/MeshElement.hxx index b08e9d74e..56f3fce53 100644 --- a/src/INTERP_KERNEL/MeshElement.hxx +++ b/src/INTERP_KERNEL/MeshElement.hxx @@ -15,105 +15,50 @@ using namespace MED_EN; 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& 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 }; }; diff --git a/src/INTERP_KERNEL/MeshRegion.cxx b/src/INTERP_KERNEL/MeshRegion.cxx index 1dacdf2b4..16f874092 100644 --- a/src/INTERP_KERNEL/MeshRegion.cxx +++ b/src/INTERP_KERNEL/MeshRegion.cxx @@ -29,7 +29,8 @@ namespace INTERP_UTILS } /** - * 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 * @@ -45,7 +46,7 @@ namespace INTERP_UTILS { 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); @@ -75,6 +76,7 @@ namespace INTERP_UTILS */ 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 @@ -97,29 +99,53 @@ namespace INTERP_UTILS 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::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::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(); diff --git a/src/INTERP_KERNEL/MeshRegion.hxx b/src/INTERP_KERNEL/MeshRegion.hxx index 4311ec60a..b190154ce 100644 --- a/src/INTERP_KERNEL/MeshRegion.hxx +++ b/src/INTERP_KERNEL/MeshRegion.hxx @@ -20,36 +20,12 @@ namespace INTERP_UTILS { 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; @@ -61,9 +37,11 @@ namespace INTERP_UTILS 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 _elements; + + /// BoundingBox containing all the nodes of all the elements in the region. BoundingBox* _box; }; diff --git a/src/INTERP_KERNEL/MeshUtils.hxx b/src/INTERP_KERNEL/MeshUtils.hxx new file mode 100644 index 000000000..498ea0aa5 --- /dev/null +++ b/src/INTERP_KERNEL/MeshUtils.hxx @@ -0,0 +1,54 @@ +#ifndef __MESH_UTILS_HXX__ +#define __MESH_UTILS_HXX__ + +#include "MEDMEM_define.hxx" +#include + +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(type) - 300; + } + +}; + + + + + +#endif diff --git a/src/INTERP_KERNEL/RegionNode.cxx b/src/INTERP_KERNEL/RegionNode.cxx index 6f965e91d..d49b3d9cd 100644 --- a/src/INTERP_KERNEL/RegionNode.cxx +++ b/src/INTERP_KERNEL/RegionNode.cxx @@ -17,7 +17,6 @@ namespace INTERP_UTILS RegionNode::~RegionNode() { } - /** * Accessor to source region diff --git a/src/INTERP_KERNEL/RegionNode.hxx b/src/INTERP_KERNEL/RegionNode.hxx index 24caeadca..b8352c368 100644 --- a/src/INTERP_KERNEL/RegionNode.hxx +++ b/src/INTERP_KERNEL/RegionNode.hxx @@ -10,43 +10,23 @@ namespace INTERP_UTILS * 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 }; diff --git a/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx b/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx index 9276d95e1..cc365c224 100644 --- a/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx +++ b/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx @@ -9,16 +9,16 @@ #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; @@ -56,15 +56,20 @@ bool Interpolation3DTest::areCompatitable(const IntersectionMatrix& m1, const In 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; } @@ -91,15 +96,20 @@ bool Interpolation3DTest::testSymmetric(const IntersectionMatrix& m1, const Inte 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; } @@ -116,15 +126,20 @@ bool Interpolation3DTest::testDiagonal(const IntersectionMatrix& m) const 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; } @@ -173,7 +188,7 @@ void Interpolation3DTest::calcIntersectionMatrix(const char* mesh1path, const ch } -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 << "=============================" ); @@ -186,33 +201,32 @@ void Interpolation3DTest::intersectMeshes(const char* mesh1path, const char* mes 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); diff --git a/src/INTERP_KERNEL/Test/Interpolation3DTest.hxx b/src/INTERP_KERNEL/Test/Interpolation3DTest.hxx index 7de5c6a24..91a76c6b7 100644 --- a/src/INTERP_KERNEL/Test/Interpolation3DTest.hxx +++ b/src/INTERP_KERNEL/Test/Interpolation3DTest.hxx @@ -28,13 +28,11 @@ class Interpolation3DTest : public CppUnit::TestFixture CPPUNIT_TEST( tetraSimpleHalfstripOnly ); CPPUNIT_TEST( generalTetra ); CPPUNIT_TEST( trickyTetra1 ); - CPPUNIT_TEST( inconsistentTetra ); // multi - element - - + CPPUNIT_TEST( tetraComplexIncluded ); CPPUNIT_TEST( dividedUnitTetraSimplerReflexive ); CPPUNIT_TEST( dividedUnitTetraReflexive ); @@ -50,6 +48,7 @@ class Interpolation3DTest : public CppUnit::TestFixture CPPUNIT_TEST( moderateBoxEvenSmallerReflexive ); CPPUNIT_TEST( tinyBoxReflexive ); + CPPUNIT_TEST( simpleHexaBox ); CPPUNIT_TEST_SUITE_END(); @@ -211,7 +210,12 @@ public: { 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; @@ -226,7 +230,8 @@ private: 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; diff --git a/src/INTERP_KERNEL/Test/Makefile.in b/src/INTERP_KERNEL/Test/Makefile.in index 5e08e1d13..f9d22c043 100644 --- a/src/INTERP_KERNEL/Test/Makefile.in +++ b/src/INTERP_KERNEL/Test/Makefile.in @@ -34,10 +34,7 @@ VPATH=.:@srcdir@:@top_srcdir@/idl @COMMENCE@ # header files -EXPORT_HEADERS = CppUnitTest.hxx \ - TransformedTriangleTest.hxx \ - TransformedTriangleIntersectTest.hxx \ - Interpolation3DTest.hxx +EXPORT_HEADERS = TestingUtils.hxx # Libraries targets @@ -49,10 +46,10 @@ LIB_CLIENT_IDL = # Executables targets -BIN = TestInterpKernel +BIN = PerfTest + +BIN_SRC = -BIN_SRC = CppUnitTest.cxx TransformedTriangleTest.cxx TransformedTriangleIntersectTest.cxx \ -Interpolation3DTest.cxx BIN_CLIENT_IDL = @@ -63,18 +60,35 @@ CXXFLAGS+=@CXXTMPDPTHFLAGS@ $(MED2_INCLUDES) $(HDF5_INCLUDES) -I$(top_srcdir)/sr 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} @@ -87,10 +101,13 @@ LIBS = @LIBS@ @CPPUNIT_LIBS@ 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@ diff --git a/src/INTERP_KERNEL/Test/PerfTest.cxx b/src/INTERP_KERNEL/Test/PerfTest.cxx new file mode 100644 index 000000000..958152527 --- /dev/null +++ b/src/INTERP_KERNEL/Test/PerfTest.cxx @@ -0,0 +1,25 @@ +#include "Interpolation3D.hxx" +#include "TestingUtils.hxx" +#include +#include + +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; + +} diff --git a/src/INTERP_KERNEL/Test/TestInterpKernel.cxx b/src/INTERP_KERNEL/Test/TestInterpKernel.cxx index 4fa722002..39793a3e3 100644 --- a/src/INTERP_KERNEL/Test/TestInterpKernel.cxx +++ b/src/INTERP_KERNEL/Test/TestInterpKernel.cxx @@ -25,9 +25,9 @@ // --- 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 diff --git a/src/INTERP_KERNEL/Test/TestingUtils.hxx b/src/INTERP_KERNEL/Test/TestingUtils.hxx new file mode 100644 index 000000000..23ff398c6 --- /dev/null +++ b/src/INTERP_KERNEL/Test/TestingUtils.hxx @@ -0,0 +1,289 @@ +#ifndef _TESTING_UTILS_HXX_ +#define _TESTING_UTILS_HXX_ + +#include "Interpolation3D.hxx" +#include "MEDMEM_Mesh.hxx" + +#include +#include +#include +#include +#include + +#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 volumes; + for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter) + { + for(map::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::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::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 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::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::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 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::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 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 diff --git a/src/INTERP_KERNEL/TetraAffineTransform.cxx b/src/INTERP_KERNEL/TetraAffineTransform.cxx index f2f127313..2f3baddc9 100644 --- a/src/INTERP_KERNEL/TetraAffineTransform.cxx +++ b/src/INTERP_KERNEL/TetraAffineTransform.cxx @@ -1 +1,324 @@ #include "TetraAffineTransform.hxx" +#include "VectorUtils.hxx" + +#include +#include + +#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]]; + } + + +}; diff --git a/src/INTERP_KERNEL/TetraAffineTransform.hxx b/src/INTERP_KERNEL/TetraAffineTransform.hxx index bc0bc5712..03d2159a6 100644 --- a/src/INTERP_KERNEL/TetraAffineTransform.hxx +++ b/src/INTERP_KERNEL/TetraAffineTransform.hxx @@ -1,307 +1,49 @@ #ifndef __TETRA_AFFINE_TRANSFORM_HXX__ #define __TETRA_AFFINE_TRANSFORM_HXX__ -#include -#include - -#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; }; diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index fd135e74d..3cea450fd 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -14,55 +14,6 @@ #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 { @@ -93,7 +44,7 @@ private: 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]))); @@ -135,15 +86,6 @@ namespace INTERP_UTILS _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]; @@ -160,6 +102,11 @@ namespace INTERP_UTILS // initialise rest of data preCalculateDoubleProducts(); + +#ifdef OPTIMIZE + preCalculateTriangleSurroundsEdge(); +#endif + preCalculateTripleProducts(); } @@ -187,14 +134,13 @@ namespace INTERP_UTILS */ 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 @@ -206,7 +152,6 @@ namespace INTERP_UTILS 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"); @@ -217,7 +162,6 @@ namespace INTERP_UTILS // normalize sign = sign > 0.0 ? 1.0 : -1.0; - LOG(2, "-- Calculating intersection polygons ... "); calculateIntersectionPolygons(); @@ -250,16 +194,16 @@ namespace INTERP_UTILS #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 << "***********"); @@ -287,13 +231,19 @@ namespace INTERP_UTILS 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); @@ -320,10 +270,96 @@ namespace INTERP_UTILS 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)) { @@ -340,6 +376,7 @@ namespace INTERP_UTILS _polygonB.push_back(ptB); LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B"); } + } } @@ -379,6 +416,46 @@ namespace INTERP_UTILS } } } +#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)) @@ -403,6 +480,8 @@ namespace INTERP_UTILS LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B"); } } + +#endif } // inclusion tests @@ -599,11 +678,11 @@ namespace INTERP_UTILS #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; @@ -656,14 +735,14 @@ namespace INTERP_UTILS 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 diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 63a95813d..f506bcc1d 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -9,19 +9,23 @@ #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 { @@ -33,6 +37,41 @@ 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 @@ -43,7 +82,6 @@ namespace INTERP_UTILS friend class ::TransformedTriangleTest; friend class ::TransformedTriangleIntersectTest; - friend class ::TransformedTriangleCalcVolumeTest; /** @@ -77,7 +115,6 @@ namespace INTERP_UTILS double calculateIntersectionVolume(); - // temporary debug method void dumpCoords(); @@ -110,11 +147,11 @@ namespace INTERP_UTILS /// 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; @@ -124,7 +161,7 @@ namespace INTERP_UTILS 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); @@ -132,11 +169,11 @@ namespace INTERP_UTILS 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 /////////////// @@ -164,7 +201,7 @@ namespace INTERP_UTILS 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; @@ -172,11 +209,11 @@ namespace INTERP_UTILS 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; @@ -185,19 +222,39 @@ namespace INTERP_UTILS //////////////////////////////////////////////////////////////////////////////////// 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 _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 ///////////////// @@ -245,8 +302,33 @@ namespace INTERP_UTILS // 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 diff --git a/src/INTERP_KERNEL/TransformedTriangle_inline.hxx b/src/INTERP_KERNEL/TransformedTriangle_inline.hxx new file mode 100644 index 000000000..9b805b529 --- /dev/null +++ b/src/INTERP_KERNEL/TransformedTriangle_inline.hxx @@ -0,0 +1,191 @@ +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; + + } diff --git a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx index 49c152e24..6a1aed56e 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx @@ -68,11 +68,63 @@ namespace INTERP_UTILS 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]) * @@ -83,7 +135,7 @@ namespace INTERP_UTILS { 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]) @@ -124,6 +176,7 @@ namespace INTERP_UTILS } } +#ifndef OPTIMIZE /** * Tests if the given segment of the triangle intersects the given facet of the tetrahedron. * (Grandy, eq. [19]) @@ -136,6 +189,7 @@ namespace INTERP_UTILS { return testFacetSurroundsSegment(seg, facet) && testSegmentIntersectsFacet(seg, facet); } +#endif /** * Calculates the point of intersection between the given segment of the triangle @@ -192,16 +246,8 @@ namespace INTERP_UTILS */ 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)) @@ -212,7 +258,8 @@ namespace INTERP_UTILS { 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 @@ -339,14 +386,7 @@ namespace INTERP_UTILS */ 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] = @@ -357,6 +397,7 @@ namespace INTERP_UTILS 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) { @@ -372,6 +413,7 @@ namespace INTERP_UTILS return false; } } +#endif // check segment intersect a facet for(int i = 0 ; i < 3 ; ++i) @@ -386,7 +428,7 @@ namespace INTERP_UTILS 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]) @@ -398,6 +440,7 @@ namespace INTERP_UTILS { return testTriangleSurroundsRay( corner ) && testSurfaceAboveCorner( corner ); } +#endif /** * Tests if the given segment of the triangle intersects the half-strip above the @@ -422,23 +465,9 @@ namespace INTERP_UTILS 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 @@ -457,7 +486,7 @@ namespace INTERP_UTILS // 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] ); @@ -488,15 +517,8 @@ namespace INTERP_UTILS // 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); @@ -539,27 +561,6 @@ namespace INTERP_UTILS // readjust index since O is not used const int cornerIdx = static_cast(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] = @@ -572,7 +573,8 @@ namespace INTERP_UTILS 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 @@ -585,6 +587,8 @@ namespace INTERP_UTILS 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); @@ -615,6 +619,7 @@ namespace INTERP_UTILS } +#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) @@ -642,6 +647,7 @@ namespace INTERP_UTILS 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) @@ -684,7 +690,7 @@ namespace INTERP_UTILS return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0; } - +#endif //////////////////////////////////////////////////////////////////////////////////// diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index dd3081a5e..49de4a361 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -69,7 +69,8 @@ namespace INTERP_UTILS _doubleProducts[8*seg + dp] = calcUnstableC(seg, dp); } } - + + std::map distances; // -- (1) for each segment : check that double products satisfy Grandy, [46] // -- and make corrections if not @@ -77,7 +78,7 @@ namespace INTERP_UTILS { if(!areDoubleProductsConsistent(seg)) { - std::map distances; + LOG(4, "inconsistent! "); for(TetraCorner corner = O ; corner <= Z ; corner = TetraCorner(corner + 1)) @@ -101,6 +102,7 @@ namespace INTERP_UTILS } #endif } + distances.clear(); } @@ -151,8 +153,8 @@ namespace INTERP_UTILS 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); @@ -175,6 +177,7 @@ namespace INTERP_UTILS // 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 @@ -193,7 +196,7 @@ namespace INTERP_UTILS _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) @@ -238,16 +241,19 @@ namespace INTERP_UTILS 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 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 anglesForRows; for(int row = 1 ; row < 4 ; ++row) { @@ -257,7 +263,11 @@ namespace INTERP_UTILS 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); @@ -291,6 +301,7 @@ namespace INTERP_UTILS _validTP[corner] = false; } + anglesForRows.clear(); } @@ -349,6 +360,7 @@ namespace INTERP_UTILS } +#ifndef OPTIMIZE /** * Returns the stable double product c_{xy}^{pq} * @@ -365,7 +377,6 @@ namespace INTERP_UTILS 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 @@ -384,6 +395,7 @@ namespace INTERP_UTILS 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 @@ -410,6 +422,8 @@ namespace INTERP_UTILS 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 diff --git a/src/INTERP_KERNEL/VectorUtils.hxx b/src/INTERP_KERNEL/VectorUtils.hxx index aa434e3b1..c969b9af0 100644 --- a/src/INTERP_KERNEL/VectorUtils.hxx +++ b/src/INTERP_KERNEL/VectorUtils.hxx @@ -9,11 +9,21 @@ #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) @@ -22,6 +32,12 @@ namespace INTERP_UTILS } } + /* + * 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); @@ -29,6 +45,13 @@ namespace INTERP_UTILS 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]; @@ -36,34 +59,55 @@ namespace INTERP_UTILS 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 @@ -73,9 +117,9 @@ namespace INTERP_UTILS 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; } };