--- /dev/null
+#include "BoundingBox.hxx"
+
+#include <algorithm>
+
+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);
+
+ }
+ }
+
+};
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
*/
double getCoordinate(const BoxCoord coord) const;
+ void updateWithPoint(const double* pt);
+
private:
double[6] _coords;
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-* /* fonction pour reconstituer un polygone convexe à partir */
+ /* fonction pour reconstituer un polygone convexe à partir */
/* d'un nuage de point. */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
--- /dev/null
+#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<int>(_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<TransformedTriangle>& triangles, const TetraAffineTransform& T) const
+ {
+ // not implemented
+ }
+
+};
*/
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<TransformedTriangle>& triangles, const TetraAffineTransform& T);
+
private:
const int _index;
const MEDMEM::Mesh* _mesh;
--- /dev/null
+#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<MeshElement*> _elements;
+ BoundingBox* _box;
+
+ };
+
+};
+
+
+#endif
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;
