From 96421e00b3484d0bd029e46126b04879ee9139f9 Mon Sep 17 00:00:00 2001 From: vbd Date: Tue, 11 Sep 2007 13:01:48 +0000 Subject: [PATCH] staffan : * doc updates * change namespace IntersectorTetra * added 2nd constructor to IntersectorTetra --- src/INTERP_KERNEL/Interpolation3D.cxx | 4 - src/INTERP_KERNEL/IntersectorTetra.cxx | 61 +++++--- src/INTERP_KERNEL/IntersectorTetra.hxx | 207 ++++++++++++++++++------- 3 files changed, 190 insertions(+), 82 deletions(-) diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx index b3d592477..f243e979d 100644 --- a/src/INTERP_KERNEL/Interpolation3D.cxx +++ b/src/INTERP_KERNEL/Interpolation3D.cxx @@ -113,7 +113,6 @@ namespace MEDMEM // create empty maps for all source elements matrix.resize(numTargetElems); - int totalFiltered = 0; #ifdef USE_RECURSIVE_BBOX_FILTER @@ -192,7 +191,6 @@ namespace MEDMEM } } - totalFiltered += intersector.filtered; } else // recursion @@ -330,7 +328,6 @@ namespace MEDMEM } } - totalFiltered += intersector.filtered; } #endif @@ -346,7 +343,6 @@ namespace MEDMEM delete targetElems[i]; } - LOG(2, "Intersector filtered out " << totalFiltered << " elements" << std::endl); } diff --git a/src/INTERP_KERNEL/IntersectorTetra.cxx b/src/INTERP_KERNEL/IntersectorTetra.cxx index 33baf6eec..eb2b089ef 100644 --- a/src/INTERP_KERNEL/IntersectorTetra.cxx +++ b/src/INTERP_KERNEL/IntersectorTetra.cxx @@ -15,38 +15,60 @@ using namespace MEDMEM; using namespace MED_EN; using namespace INTERP_UTILS; -namespace MEDMEM +/// Smallest volume of the intersecting elements in the transformed space that will be returned as non-zero. +/// Since the scale is always the same in the transformed space (the target tetrahedron is unitary), this number is independent of the scale of the meshes. +#define SPARSE_TRUNCATION_LIMIT 1.0e-14 + +namespace INTERP_UTILS { + /** - * Constructor + * Constructor creating object from target cell global number * * @param srcMesh mesh containing the source elements * @param targetMesh mesh containing the target elements * @param targetCell global number of the target cell + * */ - IntersectorTetra::IntersectorTetra(const MESH& srcMesh, const MESH& targetMesh, int targetCell) - : _srcMesh(srcMesh), filtered(0) + IntersectorTetra::IntersectorTetra(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, int targetCell) + : _srcMesh(srcMesh) { 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 + // get array of corners of target tetraeder const double* tetraCorners[4]; for(int i = 0 ; i < 4 ; ++i) { tetraCorners[i] = getCoordsOfNode(i + 1, targetCell, targetMesh); } - - // create AffineTransform from tetrahedron - _t = new TetraAffineTransform( tetraCorners ); + // create the affine transform + createAffineTransform(tetraCorners); + + } + + /** + * Constructor creating object from the four corners of the tetrahedron. + * + * @param srcMesh mesh containing the source elements + * @param tetraCorners array of four pointers to double[3] arrays containing the coordinates of the + * corners of the tetrahedron + */ + IntersectorTetra::IntersectorTetra(const MEDMEM::MESH& srcMesh, const double** tetraCorners) + : _srcMesh(srcMesh) + { + // create the affine transform + createAffineTransform(tetraCorners); } /** * Destructor * + * Deletes _t and the coordinates (double[3] vectors) in _nodes + * */ IntersectorTetra::~IntersectorTetra() { @@ -69,13 +91,12 @@ namespace MEDMEM * The class will cache the intermediary calculations of transformed nodes of source cells and volumes associated * with triangulated faces to avoid having to recalculate these. * - * @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 element global number of the source element (1 <= srcCell < # source cells) */ double IntersectorTetra::intersectSourceCell(int element) { - //{ could be done on outside + //{ could be done on outside? // check if we have planar tetra element if(_t->determinant() == 0.0) { @@ -103,16 +124,12 @@ namespace MEDMEM { // we could store mapping local -> global numbers too, but not sure it is worth it const int globalNodeNum = getGlobalNumberOfNode(i, element, _srcMesh); - if(_nodes.find(globalNodeNum) == _nodes.end()) + if(_nodes.find(globalNodeNum) == _nodes.end()) { - calculateNode(globalNodeNum); + const double* node = calculateNode(globalNodeNum); } - checkIsOutside(_nodes[globalNodeNum], isOutside); - - // local caching of globalNodeNum - // not sure this is efficient - // globalNodeNumbers[i] = globalNodeNum; + checkIsOutside(_nodes[globalNodeNum], isOutside); } // halfspace filtering check @@ -145,7 +162,7 @@ namespace MEDMEM assert(locNodeNum >= 1); assert(locNodeNum <= cellModel.getNumberOfNodes()); - faceNodes[j-1] = getGlobalNumberOfNode(locNodeNum, element, _srcMesh);//globalNodeNumbers[locNodeNum]; + faceNodes[j-1] = getGlobalNumberOfNode(locNodeNum, element, _srcMesh); } switch(faceType) @@ -219,14 +236,10 @@ namespace MEDMEM } } } - else - { - ++filtered; - } // reset if it is very small to keep the matrix sparse // is this a good idea? - if(epsilonEqual(totalVolume, 0.0, 1.0e-14)) + if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT)) { totalVolume = 0.0; } diff --git a/src/INTERP_KERNEL/IntersectorTetra.hxx b/src/INTERP_KERNEL/IntersectorTetra.hxx index df639fad6..6ebf41588 100644 --- a/src/INTERP_KERNEL/IntersectorTetra.hxx +++ b/src/INTERP_KERNEL/IntersectorTetra.hxx @@ -16,38 +16,127 @@ using __gnu_cxx::hash_map; namespace INTERP_UTILS { + /** + * Class representing a triangular face, used as key in caching hash map in IntersectorTetra. + * + */ class TriangleFaceKey { public: - int _nodes[3]; - int _hashVal; + /** + * Constructor + * Sorts the given nodes (so that the order in which they are passed does not matter) and + * calculates a hash value for the key. + * + * @param node1 global number of the first node of the face + * @param node2 global number of the second node of the face + * @param node3 global number of the third node of the face + */ TriangleFaceKey(int node1, int node2, int node3) { sort3Ints(_nodes, node1, node2, node3); _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29; } - bool operator==(const TriangleFaceKey& rhs) const + /** + * Equality comparison operator. + * Compares this TriangleFaceKey object to another and determines if they represent the same face. + * + * @param key TriangleFaceKey with which to compare + * @return true if key has the same three nodes as this object, false if not + */ + bool operator==(const TriangleFaceKey& key) const { - return _nodes[0] == rhs._nodes[0] && _nodes[1] == rhs._nodes[1] && _nodes[2] == rhs._nodes[2]; + return _nodes[0] == key._nodes[0] && _nodes[1] == key._nodes[1] && _nodes[2] == key._nodes[2]; } + /** + * Returns a hash value for the object, based on its three nodes. + * This value is not unique for each face. + * + * @return a hash value for the object + */ int hashVal() const { return _hashVal; } - + inline void sort3Ints(int* sorted, int node1, int node2, int node3); + + private: + /// global numbers of the three nodes, sorted in ascending order + int _nodes[3]; + + /// hash value for the object, calculated in the constructor + int _hashVal; }; + + /** + * Method to sort three integers in ascending order + * + * @param sorted int[3] array in which to store the result + * @param x1 first integer + * @param x2 second integer + * @param x3 third integer + */ + inline void TriangleFaceKey::sort3Ints(int* sorted, int x1, int x2, int x3) + { + if(x1 < x2) + { + if(x1 < x3) + { + // x1 is min + sorted[0] = x1; + sorted[1] = x2 < x3 ? x2 : x3; + sorted[2] = x2 < x3 ? x3 : x2; + } + else + { + // x3, x1, x2 + sorted[0] = x3; + sorted[1] = x1; + sorted[2] = x2; + } + } + else // x2 < x1 + { + if(x2 < x3) + { + // x2 is min + sorted[0] = x2; + sorted[1] = x1 < x3 ? x1 : x3; + sorted[2] = x1 < x3 ? x3 : x1; + } + else + { + // x3, x2, x1 + sorted[0] = x3; + sorted[1] = x2; + sorted[2] = x1; + } + } + } } namespace __gnu_cxx { + + /** + * Template specialization of __gnu_cxx::hash function object for use with a __gnu_cxx::hash_map + * with TriangleFaceKey as key class. + * + */ template<> class hash { public: + /** + * Operator() that returns the precalculated hashvalue of a TriangleFaceKey object. + * + * @param key a TriangleFaceKey object + * @return an integer hash value for key + */ int operator()(const INTERP_UTILS::TriangleFaceKey& key) const { return key.hashVal(); @@ -58,24 +147,31 @@ namespace __gnu_cxx namespace INTERP_UTILS { - - /** - * Class calculating the volume of intersection of individual 3D elements. + /** + * Class calculating the volume of intersection between a tetrahedral target element and + * source elements with triangular or quadratilateral faces. * */ class IntersectorTetra { - public : - IntersectorTetra(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, int targetCell); + public: + IntersectorTetra(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, int targetCell); + + IntersectorTetra(const MEDMEM::MESH& srcMesh, const double** tetraCorners); + ~IntersectorTetra(); double intersectSourceCell(int srcCell); - - mutable int filtered; private: + + // member functions + inline void createAffineTransform(const double* corners); + inline void checkIsOutside(const double* pt, bool* isOutside) const; + inline void calculateNode(int globalNodeNum); + inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key); /// disallow copying IntersectorTetra(const IntersectorTetra& t); @@ -83,19 +179,42 @@ namespace INTERP_UTILS /// disallow assignment IntersectorTetra& operator=(const IntersectorTetra& t); + // member variables + /// affine transform associated with this target element TetraAffineTransform* _t; + /// hash_map relating node numbers to transformed nodes, used for caching hash_map< int, double*> _nodes; - hash_map< TriangleFaceKey, double > _volumes; - inline void checkIsOutside(const double* pt, bool* isOutside) const; - inline void calculateNode(int globalNodeNum); - inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key); + /// hash_map relating triangular faces to calculated volume contributions, used for caching + hash_map< TriangleFaceKey, double > _volumes; + /// reference to the source mesh const MEDMEM::MESH& _srcMesh; }; + /** + * Creates the affine transform _t from the corners of the tetrahedron. Used by the constructors + * + * @param corners double*[4] array containing pointers to four double[3] arrays with the + * coordinates of the corners of the tetrahedron + */ + inline void IntersectorTetra::createAffineTransform(const double** corners) + { + // create AffineTransform from tetrahedron + _t = new TetraAffineTransform( tetraCorners ); + } + + /** + * Function used to filter out elements by checking if they belong to one of the halfspaces + * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables + * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace, + * the corresponding array element will be true if and only if it was true when the method was called and pt lies in the halfspace. + * + * @param pt double[3] containing the coordiantes of a transformed point + * @param isOutside bool[8] which indicate the results of earlier checks. + */ inline void IntersectorTetra::checkIsOutside(const double* pt, bool* isOutside) const { isOutside[0] = isOutside[0] && (pt[0] <= 0.0); @@ -108,10 +227,20 @@ namespace INTERP_UTILS isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] >= 1.0); } + /** + * Calculates the transformed node with a given global node number. + * Gets the coordinates for the node in _srcMesh with the given global number and applies TetraAffineTransform + * _t to it. Stores the result in the cache _nodes. The non-existance of the node in _nodes should be verified before + * calling. + * + * @param globalNodeNum global node number of the node in the mesh _srcMesh + * + */ inline void IntersectorTetra::calculateNode(int globalNodeNum) { assert(globalNodeNum >= 0); assert(globalNodeNum < _srcMesh.getNumberOfNodes()); + const double* node = &(_srcMesh.getCoordinates(MED_EN::MED_FULL_INTERLACE)[3*globalNodeNum]); double* transformedNode = new double[3]; assert(transformedNode != 0); @@ -120,50 +249,20 @@ namespace INTERP_UTILS _nodes[globalNodeNum] = transformedNode; } + /** + * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in . + * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation. + * + * @param tri triangle for which to calculate the volume contribution + * @param key key associated with the face + */ inline void IntersectorTetra::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key) { const double vol = tri.calculateIntersectionVolume(); _volumes.insert(make_pair(key, vol)); } - inline void TriangleFaceKey::sort3Ints(int* sorted, int node1, int node2, int node3) - { - if(node1 < node2) - { - if(node1 < node3) - { - // node 1 is min - sorted[0] = node1; - sorted[1] = node2 < node3 ? node2 : node3; - sorted[2] = node2 < node3 ? node3 : node2; - } - else - { - // node3 , node1, node2 - sorted[0] = node3; - sorted[1] = node1; - sorted[2] = node2; - } - } - else // x2 < x1 - { - if(node2 < node3) - { - // node 2 is min - sorted[0] = node2; - sorted[1] = node1 < node3 ? node1 : node3; - sorted[2] = node1 < node3 ? node3 : node1; - } - else - { - // node 3, node 2, node 1 - sorted[0] = node3; - sorted[1] = node2; - sorted[2] = node1; - } - } - } - + }; -- 2.39.2