From 731cf01250b97fd57e8a635ddad6655a4a88715e Mon Sep 17 00:00:00 2001 From: vbd Date: Thu, 5 Jul 2007 14:34:30 +0000 Subject: [PATCH] ronnas : Added files pertaining to project "calculation of intersection volumes in 3D" Most .hxx files contain only the class interfaces for the moment. Implementation has begun for the TransformedTriangle class, whose code has been split into three .cxx - files --- src/INTERP_KERNEL/BoundingBox.hxx | 81 +++++++++ src/INTERP_KERNEL/Interpolation2D.hxx | 2 - src/INTERP_KERNEL/Interpolation3D.hxx | 99 +++++++++++ src/INTERP_KERNEL/Makefile.in | 20 ++- src/INTERP_KERNEL/MeshElement.hxx | 84 ++++++++++ src/INTERP_KERNEL/MeshRegion.hxx | 60 +++++++ src/INTERP_KERNEL/RegionNode.hxx | 49 ++++++ src/INTERP_KERNEL/TransformedTriangle.cxx | 136 +++++++++++++++ src/INTERP_KERNEL/TransformedTriangle.hxx | 194 ++++++++++++++++++++++ 9 files changed, 716 insertions(+), 9 deletions(-) create mode 100644 src/INTERP_KERNEL/BoundingBox.hxx create mode 100644 src/INTERP_KERNEL/Interpolation3D.hxx create mode 100644 src/INTERP_KERNEL/MeshElement.hxx create mode 100644 src/INTERP_KERNEL/MeshRegion.hxx create mode 100644 src/INTERP_KERNEL/RegionNode.hxx create mode 100644 src/INTERP_KERNEL/TransformedTriangle.cxx create mode 100644 src/INTERP_KERNEL/TransformedTriangle.hxx diff --git a/src/INTERP_KERNEL/BoundingBox.hxx b/src/INTERP_KERNEL/BoundingBox.hxx new file mode 100644 index 000000000..8a6662d01 --- /dev/null +++ b/src/INTERP_KERNEL/BoundingBox.hxx @@ -0,0 +1,81 @@ +#ifndef __BOUNDING_BOX_HXX__ +#define __BOUNDING_BOX_HXX__ + +namespace INTERP_UTILS +{ + + /** + * Class representing the bounding box of a number of points. + * + */ + class BoundingBox + { + public: + enum BoxCoord { XMIN = 0, XMAX = 1, YMIN = 2, YMAX = 3, ZMIN = 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& box1) {} + + /** + * 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; + + private: + double[6] _coords; + + }; + + // typedef BoundingBox::BoxCoord BoxCoord; + +}; + +#endif diff --git a/src/INTERP_KERNEL/Interpolation2D.hxx b/src/INTERP_KERNEL/Interpolation2D.hxx index 2a46aab77..a66814fef 100755 --- a/src/INTERP_KERNEL/Interpolation2D.hxx +++ b/src/INTERP_KERNEL/Interpolation2D.hxx @@ -5,8 +5,6 @@ #include "MEDMEM_Field.hxx" - - #include #include #include diff --git a/src/INTERP_KERNEL/Interpolation3D.hxx b/src/INTERP_KERNEL/Interpolation3D.hxx new file mode 100644 index 000000000..fba170134 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation3D.hxx @@ -0,0 +1,99 @@ +#ifndef __INTERPOLATION_3D_HXX__ +#define __INTERPOLATION_3D_HXX__ + +// MEDMEM includes +#include "Interpolation.hxx" + + +// standard includes +#include +#include + +// NOTE : Method comments should be moved to source file during the implementation phase + +// typedefs +typedef std::vector< std::map< int, double > > IntersectionMatrix; + +namespace MEDMEM +{ + /** + * Class used to calculate the volumes of intersection between the elements of two 3D meshes. + * + */ + class Interpolation3D : public Interpolation + { + + 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 Element srcElement&, const Element targetElement); + + }; + +}; + +#endif diff --git a/src/INTERP_KERNEL/Makefile.in b/src/INTERP_KERNEL/Makefile.in index 42f8b876b..17d330ba9 100644 --- a/src/INTERP_KERNEL/Makefile.in +++ b/src/INTERP_KERNEL/Makefile.in @@ -51,31 +51,37 @@ BBTree.H LIB=libinterpkernel.la LIB_SRC = \ +TransformedTriangle.cxx\ +TransformedTriangle_intersect.cxx\ +TransformedTriangle_math.cxx\ +Interpolation3DSurf.cxx\ Interpolation2D.cxx\ -Interpolation3DSurf.cxx - - # Executables targets -BIN = +BIN = BIN_SRC = BIN_SERVER_IDL = BIN_CLIENT_IDL = -TEST_PROGS = make_plane test_INTERPOL_2D +TEST_PROGS = make_plane test_INTERPOL_2D LDFLAGS+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome CXXFLAGS+=@CXXTMPDPTHFLAGS@ $(MED2_INCLUDES) -CPPFLAGS+=$(BOOST_CPPFLAGS) +CPPFLAGS+=$(BOOST_CPPFLAGS) + +# enable testing +CXXFLAGS+= -DTESTING_INTERP_KERNEL +CPPFLAGS+= -DTESTING_INTERP_KERNEL + #LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 $(STDLIB) -lmedmem -lutil #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 $(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 -lmedmem $(BOOST_LIBS) -lutil ifeq ($(MED_WITH_KERNEL),yes) CPPFLAGS+= ${KERNEL_CXXFLAGS} diff --git a/src/INTERP_KERNEL/MeshElement.hxx b/src/INTERP_KERNEL/MeshElement.hxx new file mode 100644 index 000000000..040977b17 --- /dev/null +++ b/src/INTERP_KERNEL/MeshElement.hxx @@ -0,0 +1,84 @@ +#ifndef __MESH_ELEMENT_HXX__ +#define __MESH_ELEMENT_HXX__ + +namespace INTERP_UTILS +{ + + /** + * 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. + */ + 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); + + private: + const int _index; + const MEDMEM::Mesh* _mesh; + const MEDEN::medGeometryElement _type; + BoundingBox _box; // should not change after initialisation + + }; + + // to be added : triangulation methods + +}; + + +#endif diff --git a/src/INTERP_KERNEL/MeshRegion.hxx b/src/INTERP_KERNEL/MeshRegion.hxx new file mode 100644 index 000000000..fdbb0df15 --- /dev/null +++ b/src/INTERP_KERNEL/MeshRegion.hxx @@ -0,0 +1,60 @@ +#ifndef __MESH_REGION_HXX__ +#define __MESH_REGION_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. + * + */ + class MeshRegion + { + 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); + + /** + * 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(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/RegionNode.hxx b/src/INTERP_KERNEL/RegionNode.hxx new file mode 100644 index 000000000..e21ea864f --- /dev/null +++ b/src/INTERP_KERNEL/RegionNode.hxx @@ -0,0 +1,49 @@ +#ifndef __REGION_NODE_HXX__ +#define __REGION_NODE_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 + * in the bounding-box filtering process. + * + */ + class RegionNode + { + public: + + /** + * Default constructor + * + */ + RegionNode(); + + /** + * Destructor + * + */ + ~RegionNode(); + + /** + * + */ + MeshRegion& getSrcRegion() const; + + /** + * + */ + MeshRegion& getTargetRegion() const; + + + private: + + MeshRegion _srcRegion; + MeshRegion _targetRegion; + + }; + +}; + + +#endif diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx new file mode 100644 index 000000000..91e751540 --- /dev/null +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -0,0 +1,136 @@ +#include "TransformedTriangle.hxx" +#include +#include +#include + +namespace INTERP_UTILS +{ + //////////////////////////////////////////////////////////////////////////// + /// PUBLIC //////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////// + + /** + * Constructor + * + * The coordinates are copied to the internal member variables + * (slower but safer - investigate if this is necessary) + * + * @param p array of three doubles containing coordinates of P + * @param q array of three doubles containing coordinates of Q + * @param r array of three doubles containing coordinates of R + */ + TransformedTriangle::TransformedTriangle(double* p, double* q, double* r) + : _isDoubleProductsCalculated(false), _isTripleProductsCalculated(false) + { + + for(int i = 0 ; i < 3 ; ++i) + { + // xyz coordinates + _coords[5*P + i] = p[i]; + _coords[5*Q + i] = q[i]; + _coords[5*R + i] = r[i]; + } + + // h coordinate + _coords[5*P + 3] = 1 - p[0] - p[1] - p[2]; + _coords[5*Q + 3] = 1 - q[0] - q[1] - q[2]; + _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]; + _coords[5*R + 4] = 1 - r[0] - r[1]; + + // initialise rest of data + preCalculateDoubleProducts(); + preCalculateTripleProducts(); + + } + + TransformedTriangle::~TransformedTriangle() + { + } + + /** + * Calculates the volume of intersection between the triangle and the + * unit tetrahedron. + * + * @returns volume of intersection of this triangle with unit tetrahedron, + * as described in Grandy + * + */ + double TransformedTriangle::calculateIntersectionVolume() + { + // not implemented + return 0.0; + } + + //////////////////////////////////////////////////////////////////////////// + /// PRIVATE //////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////// + + //////////////////////////////////////////////////////////////////////////////////// + /// High-level methods called directly by calculateIntersectionVolume() //////// + //////////////////////////////////////////////////////////////////////////////////// + /** + * Calculates the intersection polygons A and B, performing the intersection tests + * and storing the corresponding points in the vectors _polygonA and _polygonB + * + * @post _polygonA contains the intersection polygon A and _polygonB contains the + * intersection polygon B. + * + */ + void TransformedTriangle::calculateIntersectionPolygons() + { + // not implemented + } + + /** + * Calculates the barycenters of the given intersection polygon. + * + * @pre the intersection polygons have been calculated with calculateIntersectionPolygons() + * + * @param poly one of the two intersection polygons + * @param barycenter array of three doubles where barycenter is stored + * + */ + void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter) + { + // not implemented + } + + /** + * Sorts the given intersection polygon in circular order around its barycenter. + + * @pre the intersection polygons have been calculated with calculateIntersectionPolygons() + * @post the vertices in _polygonA and _polygonB are sorted in circular order around their + * respective barycenters + * + * @param poly one of the two intersection polygons + * @param barycenter array of three doubles with the coordinates of the barycenter + * + */ + void TransformedTriangle::sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter) + { + // not implemented + } + + /** + * Calculates the volume between the given polygon and the z = 0 plane. + * + * @pre the intersection polygones have been calculated with calculateIntersectionPolygons(), + * and they have been sorted in circular order with sortIntersectionPolygons(void) + * + * @param poly one of the two intersection polygons + * @param barycenter array of three doubles with the coordinates of the barycenter + * @returns the volume between the polygon and the z = 0 plane + * + */ + double TransformedTriangle::calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter) + { + // not implemented + return -3.14; + } + + +}; // NAMESPACE diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx new file mode 100644 index 000000000..9281486ed --- /dev/null +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -0,0 +1,194 @@ +#ifndef __TRANSFORMED_TRIANGLE_HXX__ +#define __TRANSFORMED_TRIANGLE_HXX__ + +#include + +#ifdef TESTING_INTERP_KERNEL +class TransformedTriangleTest; +#endif + +namespace INTERP_UTILS +{ + + /** + * Class representing one of the faces of the triangulated source polyhedron after having been transformed + * with the affine transform that takes the target tetrahedron to the unit tetrahedron. It contains the + * logic for calculating the volume of intersection between the triangle and the unit tetrahedron. + * + * Reference : J. Grandy, "Conservative Remapping and Region Overlays by Intersecting Arbitrary Polyhedra", + * Journal of Computational Physics (1999) + * + */ + class TransformedTriangle + { + //#ifdef TESTING_INTERP_KERNEL + + //#endif + + public: +#ifdef TESTING_INTERP_KERNEL + friend class ::TransformedTriangleTest; +#endif + /** + * Enumerations representing the different geometric elements of the unit tetrahedron + * and the triangle. + */ + /// Corners of tetrahedron + enum TetraCorner { O = 0, X, Y, Z }; + + /// Edges of tetrahedron + enum TetraEdge { OX = 0, OY, OZ, XY, YZ, ZX, H01, H10 }; + + /// Facets (faces) of tetrahedron + enum TetraFacet { OYZ = 0, OZX, OXY, XYZ }; + + /// Corners of triangle + enum TriCorner { P = 0, Q, R }; + + /// Segments (edges) of triangle + enum TriSegment { PQ = 0, QR, RP }; + + /// Polygons + enum IntersectionPolygon{ A = 0, B }; + + /// Double products + /// NB : order corresponds to TetraEdges (Grandy, table III) + enum DoubleProduct { C_YZ = 0, C_ZX, C_XY, C_ZH, C_XH, C_YH, C_01, C_10 }; + + TransformedTriangle(double* p, double* q, double* r); + ~TransformedTriangle(); + + double calculateIntersectionVolume(); + + private: + + //////////////////////////////////////////////////////////////////////////////////// + /// High-level methods called directly by calculateIntersectionVolume() //////// + //////////////////////////////////////////////////////////////////////////////////// + + void calculateIntersectionPolygons(); + + void calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter); + + void sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter); + + double calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter); + + + //////////////////////////////////////////////////////////////////////////////////// + /// Intersection test methods and intersection point calculations //////// + //////////////////////////////////////////////////////////////////////////////////// + + bool testSurfaceEdgeIntersection(const TetraEdge edge) const; + + void calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const; + + bool testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const; + + void calcIntersectionPtSegmentFacet(const TriSegment seg, const TetraFacet facet, double* pt) const; + + bool testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const; + + void calcIntersectionPtSegmentEdge(const TriSegment seg, const TetraEdge edge, double* pt) const ; + + bool testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const ; + + bool testSurfaceRayIntersection(const TetraCorner corner) const; + + bool testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edg); + + void calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const; + + bool testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const; + + bool testCornerInTetrahedron(const TriCorner corner) const; + + bool testCornerOnXYZFacet(const TriCorner corner) const; + + + //////////////////////////////////////////////////////////////////////////////////// + /// Utility methods used in intersection tests /////////////// + //////////////////////////////////////////////////////////////////////////////////// + + bool testTriangleSurroundsEdge(const TetraEdge edge) const; + + bool testEdgeIntersectsTriangle(const TetraEdge edge) const; + + bool testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const; + + bool testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const; + + bool testSurfaceAboveCorner(const TetraCorner corner); + + //////////////////////////////////////////////////////////////////////////////////// + /// Double and triple product calculations /////////////// + //////////////////////////////////////////////////////////////////////////////////// + + void preCalculateDoubleProducts(void); + + void preCalculateTripleProducts(void); + + double calcStableC(const TriSegment seg, const DoubleProduct dp) const; + + double calcStableT(const TetraCorner corner) const; + + double calcUnstableC(const TriSegment seg, const DoubleProduct dp) const; + + double calcTByDevelopingRow(const TetraCorner corner, const int row = 1, const bool project = false) const; + + //////////////////////////////////////////////////////////////////////////////////// + /// Member variables /////////////// + //////////////////////////////////////////////////////////////////////////////////// + private: + + // storage : + // [ 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]; + + bool _isDoubleProductsCalculated, _isTripleProductsCalculated; + double _doubleProducts[24]; + double _tripleProducts[4]; + std::vector _polygonA, _polygonB; + double _barycenterA[3], _barycenterB[3]; + + //////////////////////////////////////////////////////////////////////////////////// + /// Constants ///////////////// + //////////////////////////////////////////////////////////////////////////////////// + + // offsets : 0 -> x, 1 -> y, 2 -> z, 3 -> h, 4 -> H + // corresponds to order of double products in DoubleProduct + // so that offset[C_*] gives the right coordinate + static const int DP_OFFSET_1[8]; + static const int DP_OFFSET_2[8]; + + // the coordinates used in the expansion of triple products by a given row + // in constellation (corner, row-1) + // (0,1,2,3) <=> (x,y,z,h) + static const int COORDINATE_FOR_DETERMINANT_EXPANSION[12]; + + // contains the edge of the double product used when + // expanding the triple product determinant associated with each corner + // by a given row + static const DoubleProduct DP_FOR_DETERMINANT_EXPANSION[12]; + + // values used to decide how imprecise the double products + // should be to set them to 0.0 + static const double MACH_EPS; // machine epsilon + static const double MULT_PREC_F; // precision of multiplications (Grandy : f) + static const double THRESHOLD_F; // threshold for zeroing (Grandy : F/f) + + static const double TRIPLE_PRODUCT_ANGLE_THRESHOLD; + + // correspondance facet - double product + // Grandy, table IV + static const DoubleProduct DP_FOR_SEG_FACET_INTERSECTION[12]; + + // signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION + static const double SIGN_FOR_SEG_FACET_INTERSECTION[12]; + + + }; + +}; + +#endif -- 2.39.2