From 21ef6e5429f72ddfa9e37484f867ee434b4e2098 Mon Sep 17 00:00:00 2001 From: vbd Date: Mon, 30 Jul 2007 11:59:33 +0000 Subject: [PATCH] staffan : changes before recreating directory tree --- src/INTERP_KERNEL/BoundingBox.cxx | 134 ++++++++++++++++++ src/INTERP_KERNEL/BoundingBox.hxx | 6 +- src/INTERP_KERNEL/InterpolationUtils.hxx | 2 +- src/INTERP_KERNEL/MeshElement.cxx | 114 +++++++++++++++ src/INTERP_KERNEL/MeshElement.hxx | 63 ++++++++ src/INTERP_KERNEL/MeshRegion.cxx | 74 ++++++++++ .../TransformedTriangle_math.cxx | 6 +- 7 files changed, 393 insertions(+), 6 deletions(-) create mode 100644 src/INTERP_KERNEL/BoundingBox.cxx create mode 100644 src/INTERP_KERNEL/MeshElement.cxx create mode 100644 src/INTERP_KERNEL/MeshRegion.cxx diff --git a/src/INTERP_KERNEL/BoundingBox.cxx b/src/INTERP_KERNEL/BoundingBox.cxx new file mode 100644 index 000000000..9ac02709c --- /dev/null +++ b/src/INTERP_KERNEL/BoundingBox.cxx @@ -0,0 +1,134 @@ +#include "BoundingBox.hxx" + +#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 + * 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::BoundingBox(const double** pts, const int numPts); + { + using namespace std; + assert(numPts > 1); + + // initialize with first two points + const double* pt1 = pts[0]; + const double* pt2 = pts[1]; + + for(BoxCoord c = XMIN ; c <= ZMIN ; c = BoxCoord(c + 1)) + { + _coords[c] = min(pt1[c], pt2[c]); + _coords[c + 3] = max(pt1[c + 3], pt2[c + 3]); + } + + for(int i = 2 ; i < numPts ; ++i) + { + updateWithPoint(pts[i]); + } + } + /** + * Constructor creating box from union of two boxes + * + */ + BoundingBox::BoundingBox(const BoundingBox& box1, const BoundingBox& box1) + { + 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]); + } + } + + /** + * Destructor + * + */ + BoundingBox::~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 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 + if(_coords[c] > otherMaxCoord + || _coords[c + 3] < otherMinCoord) + { + return true; + } + } + return false; + } + + /** + * Sets a coordinate of the box to a given value. + * + * @param coord coordinate to set + * @param value new value for coordinate + * + */ + void BoundingBox::setCoordinate(const BoxCoord coord, double value) + { + _coords[coord] = value; + } + + /** + * Gets a coordinate of the box + * + * @param coord coordinate to set + * @returns value of coordinate + * + */ + double BoundingBox::getCoordinate(const BoxCoord coord) const + { + return _coords[coord]; + } + + /** + * Updates the bounding box to include a given point + * + * @param pt point to be included + * + */ + void BoundingBox::updateWithPoint(const double* pt) + { + using namespace std; + + for(BoxCoord c = XMIN ; c <= ZMIN ; c = BoxCoord(c + 1)) + { + const double ptVal = pt[c]; + + // update min and max coordinates + _coords[c] = min(_coords[c], ptVal); + _coords[c + 3] = max(_coords[c + 3], ptVal); + + } + } + +}; diff --git a/src/INTERP_KERNEL/BoundingBox.hxx b/src/INTERP_KERNEL/BoundingBox.hxx index 8a6662d01..daee9bf2d 100644 --- a/src/INTERP_KERNEL/BoundingBox.hxx +++ b/src/INTERP_KERNEL/BoundingBox.hxx @@ -11,14 +11,14 @@ namespace INTERP_UTILS class BoundingBox { public: - enum BoxCoord { XMIN = 0, XMAX = 1, YMIN = 2, YMAX = 3, ZMIN = 4, ZMAX = 5 }; + 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}) {} + //BoundingBox() : coords({0.0, 0.0, 0.0, 0.0, 0.0, 0.0}) {} /** * Constructor creating box from an array of the points corresponding @@ -69,6 +69,8 @@ namespace INTERP_UTILS */ double getCoordinate(const BoxCoord coord) const; + void updateWithPoint(const double* pt); + private: double[6] _coords; diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index 59917a12d..317dd7eca 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -457,7 +457,7 @@ namespace INTERP_UTILS /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ -* /* fonction pour reconstituer un polygone convexe à partir */ + /* fonction pour reconstituer un polygone convexe à partir */ /* d'un nuage de point. */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ diff --git a/src/INTERP_KERNEL/MeshElement.cxx b/src/INTERP_KERNEL/MeshElement.cxx new file mode 100644 index 000000000..52b2c2245 --- /dev/null +++ b/src/INTERP_KERNEL/MeshElement.cxx @@ -0,0 +1,114 @@ +#include "MeshElement.hxx" + +namespace INTERP_UTILS +{ + + /** + * Constructor + * + * @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 MEDEN::medGeometryElement type, const MEDMEM::Mesh& mesh) + : _index(index - 1), _mesh(&mesh), _type(type), _box(0) + { + // get coordinates of vertices + const int numNodes = getNumberNodes(); + const double* vertices[numNodes]; + + for(int i = 0 ; i < numNodes ; ++i) + { + vertices[i] = getCoordsOfNode(i + 1); + } + + // create bounding box + _box = new BoundingBox(vertices, numNodes); + + } + + /** + * Destructor + * + */ + MeshElement::~MeshElement() + { + if(_box) + { + delete _box; + } + } + + /** + * 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 MeshElement::isElementIncludedIn(const MeshElement& otherElement) const + { + // not implemented + return false; + } + + /** + * 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 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 + */ + bool MeshElement::isElementTriviallyDisjointWith(const MeshElement& otherElement) const + { + // not implemented + return false; + } + + /** + * Returns the number of nodes of this element + * + * @returns the number of nodes of this element + */ + int MeshElement::getNumberNodes() const + { + assert(_type > 300); + assert(_type < 400); + + // int(type) = yxx where y is dimension of entity (here 3) + // and xx is the number of nodes of the element + return static_cast(_type) - 300; + } + + /** + * Returns the coordinates of a node of this element + * (1 <= node <= #nodes) + * + * @param node the node for which the coordinates are sought + * @returns pointer to an array of 3 doubles containing the coordinates + */ + const double* MeshElement::getCoordsOfNode(int node) const + { + const int nodeOffset = node - 1; + const int elemIdx = mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index]; + return &(_mesh->getCoordinates(MED_FULL_INTERLACE)[elemIdx + nodeOffset]); + } + + /** + * Triangulate the faces of this element and apply an affine Transform to the triangles + * + * @param triangles vector in which triangles are stored + * @param T affine transform that is applied to the nodes of the triangles + */ + void MeshElement::triangulate(std::vector& triangles, const TetraAffineTransform& T) const + { + // not implemented + } + +}; diff --git a/src/INTERP_KERNEL/MeshElement.hxx b/src/INTERP_KERNEL/MeshElement.hxx index 040977b17..cc060a529 100644 --- a/src/INTERP_KERNEL/MeshElement.hxx +++ b/src/INTERP_KERNEL/MeshElement.hxx @@ -68,6 +68,69 @@ namespace INTERP_UTILS */ const double* getCoordsOfNode(int node); + /** + * 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); + + /** + * 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); + private: const int _index; const MEDMEM::Mesh* _mesh; diff --git a/src/INTERP_KERNEL/MeshRegion.cxx b/src/INTERP_KERNEL/MeshRegion.cxx new file mode 100644 index 000000000..93145bf34 --- /dev/null +++ b/src/INTERP_KERNEL/MeshRegion.cxx @@ -0,0 +1,74 @@ +#include "MeshRegion.hxx" + +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. + * + */ + class MeshRegion + { + 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) + { + _elements.push_back(element); + + if(_box == 0) + { + // get coordinates of elements + _box = new BoundingBox(); + } + } + + /** + * 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; + + }; + +}; + + +#endif diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index b3c9f7c4b..2916da10b 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -29,9 +29,9 @@ namespace INTERP_UTILS C_YH, C_XH, C_XY // Z }; - const double TransformedTriangle::MACH_EPS = 1.0e-15; - const double TransformedTriangle::MULT_PREC_F = 4.0*TransformedTriangle::MACH_EPS; - const double TransformedTriangle::THRESHOLD_F = 20.0; + const long double TransformedTriangle::MACH_EPS = 1.0e-15; + const long double TransformedTriangle::MULT_PREC_F = 4.0*TransformedTriangle::MACH_EPS; + const long double TransformedTriangle::THRESHOLD_F = 20.0; const double TransformedTriangle::TRIPLE_PRODUCT_ANGLE_THRESHOLD = 0.1; -- 2.39.2