From 6373ab7d5d25de2da322e17e8bf06afc07104826 Mon Sep 17 00:00:00 2001 From: vbd Date: Tue, 31 Jul 2007 07:17:59 +0000 Subject: [PATCH] staffan : re-implementation of BoundingBox, MeshElement, MeshRegion, RegionNode * not yet tested --- src/INTERP_KERNEL/BoundingBox.cxx | 9 +- src/INTERP_KERNEL/BoundingBox.hxx | 4 +- src/INTERP_KERNEL/Interpolation3D.cxx | 6 +- src/INTERP_KERNEL/MeshElement.cxx | 66 ++++++++++- src/INTERP_KERNEL/MeshElement.hxx | 103 +++++++---------- src/INTERP_KERNEL/MeshRegion.cxx | 159 +++++++++++++++++--------- src/INTERP_KERNEL/MeshRegion.hxx | 16 ++- src/INTERP_KERNEL/RegionNode.cxx | 41 +++++++ src/INTERP_KERNEL/RegionNode.hxx | 11 +- 9 files changed, 281 insertions(+), 134 deletions(-) create mode 100644 src/INTERP_KERNEL/RegionNode.cxx diff --git a/src/INTERP_KERNEL/BoundingBox.cxx b/src/INTERP_KERNEL/BoundingBox.cxx index 9ac02709c..bb4f17fe5 100644 --- a/src/INTERP_KERNEL/BoundingBox.cxx +++ b/src/INTERP_KERNEL/BoundingBox.cxx @@ -20,7 +20,7 @@ namespace INTERP_UTILS * @param numPts number of vertices * */ - BoundingBox::BoundingBox(const double** pts, const int numPts); + BoundingBox::BoundingBox(const double** pts, const int numPts) { using namespace std; assert(numPts > 1); @@ -44,12 +44,13 @@ namespace INTERP_UTILS * Constructor creating box from union of two boxes * */ - BoundingBox::BoundingBox(const BoundingBox& box1, const BoundingBox& box1) + BoundingBox::BoundingBox(const BoundingBox& box1, const BoundingBox& box2) { + using namespace std; for(BoxCoord c = XMIN ; c <= ZMIN ; c = BoxCoord(c + 1)) { - _coords[c] = min(box1[c], box2[c]); - _coords[c + 3] = max(box1[c + 3], box2[c + 3]); + _coords[c] = min(box1._coords[c], box2._coords[c]); + _coords[c + 3] = max(box1._coords[c + 3], box2._coords[c + 3]); } } diff --git a/src/INTERP_KERNEL/BoundingBox.hxx b/src/INTERP_KERNEL/BoundingBox.hxx index daee9bf2d..64297d906 100644 --- a/src/INTERP_KERNEL/BoundingBox.hxx +++ b/src/INTERP_KERNEL/BoundingBox.hxx @@ -35,7 +35,7 @@ namespace INTERP_UTILS * Constructor creating box from union of two boxes * */ - BoundingBox(const BoundingBox& box1, const BoundingBox& box1) {} + BoundingBox(const BoundingBox& box1, const BoundingBox& box2); /** * Destructor @@ -72,7 +72,7 @@ namespace INTERP_UTILS void updateWithPoint(const double* pt); private: - double[6] _coords; + double _coords[6]; }; diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx index 6caec3d0f..69f55d794 100644 --- a/src/INTERP_KERNEL/Interpolation3D.cxx +++ b/src/INTERP_KERNEL/Interpolation3D.cxx @@ -166,7 +166,7 @@ namespace MEDMEM } else // recursion { - typedef MeshElement::ZoneBBoxComparison Cmp; + // std::cout << " - Recursion" << std::endl; RegionNode* leftNode = new RegionNode(); @@ -174,14 +174,14 @@ namespace MEDMEM // split current source region //} decide on axis - static Cmp::BoxAxis axis = Cmp::X; + static BoundingBox::BoxCoord axis = BoundingBox::XMAX; currNode->getSrcRegion().split(leftNode->getSrcRegion(), rightNode->getSrcRegion(), axis); // ugly hack to avoid problem with enum which does not start at 0 // I guess I ought to implement ++ for it instead ... // Anyway, it basically chooses the next axis, circually - axis = (axis != Cmp::Z) ? static_cast(axis + 1) : Cmp::X; + axis = (axis != BoundingBox::ZMAX) ? static_cast(axis + 1) : BoundingBox::XMAX; // add target elements of curr node that overlap the two new nodes // std::cout << " -- Adding target elements" << std::endl; diff --git a/src/INTERP_KERNEL/MeshElement.cxx b/src/INTERP_KERNEL/MeshElement.cxx index 52b2c2245..1c4d22c54 100644 --- a/src/INTERP_KERNEL/MeshElement.cxx +++ b/src/INTERP_KERNEL/MeshElement.cxx @@ -1,5 +1,8 @@ #include "MeshElement.hxx" +#include "TetraAffineTransform.hxx" +#include "TransformedTriangle.hxx" + namespace INTERP_UTILS { @@ -10,7 +13,7 @@ namespace INTERP_UTILS * @param type geometric type of the element * @param index global number of element in the mesh */ - MeshElement::MeshElement(const int index, const MEDEN::medGeometryElement type, const MEDMEM::Mesh& mesh) + MeshElement::MeshElement(const int index, const MED_EN::medGeometryElement type, const MEDMEM::MESH& mesh) : _index(index - 1), _mesh(&mesh), _type(type), _box(0) { // get coordinates of vertices @@ -96,7 +99,7 @@ namespace INTERP_UTILS const double* MeshElement::getCoordsOfNode(int node) const { const int nodeOffset = node - 1; - const int elemIdx = mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index]; + const int elemIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index]; return &(_mesh->getCoordinates(MED_FULL_INTERLACE)[elemIdx + nodeOffset]); } @@ -108,7 +111,64 @@ namespace INTERP_UTILS */ void MeshElement::triangulate(std::vector& triangles, const TetraAffineTransform& T) const { - // not implemented + switch(_type) + { + case MED_TETRA4 : + // triangles == faces + const int beginFaceIdx = _mesh->getConnectivityIndex(MED_DESCENDING, MED_CELL)[_index]; + const int endFaceIdx = _mesh->getConnectivityIndex(MED_DESCENDING, MED_CELL)[_index + 1]; + + for(int i = beginFaceIdx ; i < endFaceIdx ; ++i) // loop over faces of element + { + const int faceIdx = _mesh->getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL, MED_ALL_ELEMENTS)[i - 1]; + const int beginNodeIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_FACE)[faceIdx - 1]; + const int endNodeIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_FACE)[faceIdx]; + + double transformedPts[9]; + + assert(endNodeIdx - beginNodeIdx == 3); + + for(int j = 0 ; j < 3 ; ++j) // loop over nodes of face + { + // { optimise here using getCoordinatesForNode ? + // could maybe use the connNodeIdx directly and only transform each node once + // instead of once per face + const int connNodeIdx = + _mesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL, MED_FACE, MED_ALL_ELEMENTS)[beginNodeIdx + j - 1]; + const double* pt = &(_mesh->getCoordinates(MED_FULL_INTERLACE)[connNodeIdx - 1]); + + // transform + T.apply(&transformedPts[3*j], pt); + } + + triangles.push_back(TransformedTriangle(&transformedPts[0], &transformedPts[3], &transformedPts[6])); + } + + break; + default: + break; + } + } + + int MeshElement::getIndex() const + { + return _index + 1; + } + + + /// ElementBBoxOrder + bool ElementBBoxOrder::operator()( MeshElement* elem1, MeshElement* elem2) + { + assert(elem1 != 0); + assert(elem2 != 0); + assert(elem1->_box != 0); + assert(elem2->_box != 0); + + const double coord1 = elem1->_box->getCoordinate(_coord); + const double coord2 = elem2->_box->getCoordinate(_coord); + + return coord1 < coord2; } + }; diff --git a/src/INTERP_KERNEL/MeshElement.hxx b/src/INTERP_KERNEL/MeshElement.hxx index cc060a529..ab339da23 100644 --- a/src/INTERP_KERNEL/MeshElement.hxx +++ b/src/INTERP_KERNEL/MeshElement.hxx @@ -1,9 +1,23 @@ #ifndef __MESH_ELEMENT_HXX__ #define __MESH_ELEMENT_HXX__ +#include "MEDMEM_define.hxx" +#include "MEDMEM_Mesh.hxx" + +#include "BoundingBox.hxx" + +using namespace MEDMEM; +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 @@ -11,63 +25,13 @@ namespace INTERP_UTILS */ class MeshElement { - 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 MEDEN::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; - /** - * Returns the number of nodes of this element - * - * @returns the number of nodes of this element - */ - int getNumberNodes() const; - - /** - * Returns the coordinates of a node of this element - * (1 <= node <= #nodes) - * - * @param node the node for which the coordinates are sought - * @returns pointer to an array of 3 doubles containing the coordinates - */ - const double* getCoordsOfNode(int node); + friend class ElementBBoxOrder; + friend class MeshRegion; + public: + /** * Constructor * @@ -75,7 +39,7 @@ namespace INTERP_UTILS * @param type geometric type of the element * @param index connectivity index of element in the mesh */ - MeshElement(const int index, const MEDEN::medGeometryElement type, const MEDMEM::Mesh& mesh); + MeshElement(const int index, const MED_EN::medGeometryElement type, const MEDMEM::MESH& mesh); /** * Destructor @@ -121,7 +85,7 @@ namespace INTERP_UTILS * @param node the node for which the coordinates are sought * @returns pointer to an array of 3 doubles containing the coordinates */ - const double* getCoordsOfNode(int node); + const double* getCoordsOfNode(int node) const; /** * Triangulate the faces of this element and apply an affine Transform to the triangles @@ -129,19 +93,34 @@ namespace INTERP_UTILS * @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); + void triangulate(std::vector& triangles, const TetraAffineTransform& T) const; + + int getIndex() const; private: const int _index; - const MEDMEM::Mesh* _mesh; - const MEDEN::medGeometryElement _type; - BoundingBox _box; // should not change after initialisation + const MEDMEM::MESH* _mesh; + const MED_EN::medGeometryElement _type; + BoundingBox* _box; // should not change after initialisation }; - // to be added : triangulation methods -}; + class ElementBBoxOrder + { + public : + + ElementBBoxOrder(BoundingBox::BoxCoord coord) + : _coord(coord) + { + } + + bool operator()(MeshElement* elem1, MeshElement* elem2); + + private : + BoundingBox::BoxCoord _coord; + }; +}; #endif diff --git a/src/INTERP_KERNEL/MeshRegion.cxx b/src/INTERP_KERNEL/MeshRegion.cxx index 93145bf34..6a58263fd 100644 --- a/src/INTERP_KERNEL/MeshRegion.cxx +++ b/src/INTERP_KERNEL/MeshRegion.cxx @@ -1,74 +1,121 @@ #include "MeshRegion.hxx" +#include "MeshElement.hxx" + +//#include + namespace INTERP_UTILS { + /** - * Class representing a set of elements in a mesh together with their bounding box. - * It permits to split itself in two, which is used in the depth-first search filtering process. - * + * Default constructor + * */ - class MeshRegion + MeshRegion::MeshRegion() + : _box(0) { - public: - - /** - * Default constructor - * - */ - MeshRegion::MeshRegion() - : _box(0) - { - } + } - /** - * Destructor - * - */ - MeshRegion::~MeshRegion() - { - if(_box != 0) - { - delete _box; - } - - /** - * Adds an element to the region, updating the bounding box. - * - * @param element pointer to element to add to region - * - */ - void MeshRegion::addElement(MeshElement* const element) + /** + * Destructor + * + */ + MeshRegion::~MeshRegion() + { + if(_box != 0) { - _elements.push_back(element); + delete _box; + } + } + + /** + * Adds an element to the region, updating the bounding box. + * + * @param element pointer to element to add to region + * + */ + void MeshRegion::addElement(MeshElement* const element) + { + _elements.push_back(element); - if(_box == 0) + if(_box == 0) + { + const int numNodes = element->getNumberNodes(); + const double* pts[numNodes]; + + // get coordinates of elements + for(int i = 0 ; i < numNodes ; ++i) { - // get coordinates of elements - _box = new BoundingBox(); + pts[i] = element->getCoordsOfNode(i + 1); } + + _box = new BoundingBox(pts, numNodes); + + } + } + + /** + * 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 coord coordinate of BoundingBox to use when splitting the region + * + */ + void MeshRegion::split(MeshRegion& region1, MeshRegion& region2, BoundingBox::BoxCoord coord) + { + ElementBBoxOrder cmp(coord); + + // sort elements by their bounding boxes + std::sort(_elements.begin(), _elements.end(), cmp); + + // put the first half of the elements in region1 and the + // rest in region2 + std::vector::const_iterator iter = _elements.begin(); + int elemCount = 0; + + while(elemCount < static_cast(_elements.size() / 2)) + { + region1.addElement(*iter); + ++iter; + ++elemCount; + } + + while(iter != _elements.end()) + { + region2.addElement(*iter); + ++iter; } - /** - * 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 MeshRegion::split(Region& region1, Region& region2, int axis) 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 - std::vector _elements; - BoundingBox* _box; - - }; + // assert(std::abs(region1._elements.size() - region2._elements.size()) < 2); + + } + + bool MeshRegion::isDisjointWithElementBoundingBox(const MeshElement& elem) const + { + assert(_box != 0); + assert(elem._box != 0); + + return _box->isDisjointWith(*(elem._box)); + } + std::vector::const_iterator MeshRegion::getBeginElements() const + { + return _elements.begin(); + } + + std::vector::const_iterator MeshRegion::getEndElements() const + { + return _elements.end(); + } + + int MeshRegion::getNumberOfElements() const + { + return _elements.size(); + } }; -#endif + diff --git a/src/INTERP_KERNEL/MeshRegion.hxx b/src/INTERP_KERNEL/MeshRegion.hxx index fdbb0df15..62081aaad 100644 --- a/src/INTERP_KERNEL/MeshRegion.hxx +++ b/src/INTERP_KERNEL/MeshRegion.hxx @@ -3,8 +3,12 @@ #include +#include "BoundingBox.hxx" + namespace INTERP_UTILS { + class MeshElement; + /** * Class representing a set of elements in a mesh together with their bounding box. * It permits to split itself in two, which is used in the depth-first search filtering process. @@ -44,13 +48,21 @@ namespace INTERP_UTILS * @param axis axis along which to split the region * */ - void split(Region& region1, Region& region2, int axis) const; + void split(MeshRegion& region1, MeshRegion& region2, BoundingBox::BoxCoord coord); + + bool isDisjointWithElementBoundingBox(const MeshElement& elem) const; + + std::vector::const_iterator getBeginElements() const; + + std::vector::const_iterator getEndElements() const; + + int getNumberOfElements() const; private: // Vector of pointers to elements. NB : these pointers are not owned by the region object, and are thus // neither allocated or liberated in this class. The elements must therefore be allocated and liberated outside this class std::vector _elements; - BoundingBox _box; + BoundingBox* _box; }; diff --git a/src/INTERP_KERNEL/RegionNode.cxx b/src/INTERP_KERNEL/RegionNode.cxx new file mode 100644 index 000000000..6f965e91d --- /dev/null +++ b/src/INTERP_KERNEL/RegionNode.cxx @@ -0,0 +1,41 @@ +#include "RegionNode.hxx" + +namespace INTERP_UTILS +{ + /** + * Default constructor + * + */ + RegionNode::RegionNode() + { + } + + /** + * Destructor + * + */ + RegionNode::~RegionNode() + { + } + + + /** + * Accessor to source region + * + * @return reference to source region + */ + MeshRegion& RegionNode::getSrcRegion() + { + return _srcRegion; + } + /** + * Accessor to target region + * + * @return reference to target region + */ + MeshRegion& RegionNode::getTargetRegion() + { + return _targetRegion; + } + +}; diff --git a/src/INTERP_KERNEL/RegionNode.hxx b/src/INTERP_KERNEL/RegionNode.hxx index e21ea864f..24caeadca 100644 --- a/src/INTERP_KERNEL/RegionNode.hxx +++ b/src/INTERP_KERNEL/RegionNode.hxx @@ -1,8 +1,11 @@ #ifndef __REGION_NODE_HXX__ #define __REGION_NODE_HXX__ +#include "MeshRegion.hxx" + 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 @@ -26,14 +29,18 @@ namespace INTERP_UTILS ~RegionNode(); /** + * Accessor to source region * + * @return reference to source region */ - MeshRegion& getSrcRegion() const; + MeshRegion& getSrcRegion(); /** + * Accessor to target region * + * @return reference to target region */ - MeshRegion& getTargetRegion() const; + MeshRegion& getTargetRegion(); private: -- 2.39.2