--- /dev/null
+#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
#include "MEDMEM_Field.hxx"
-
-
#include <cmath>
#include <map>
#include <sstream>
--- /dev/null
+#ifndef __INTERPOLATION_3D_HXX__
+#define __INTERPOLATION_3D_HXX__
+
+// MEDMEM includes
+#include "Interpolation.hxx"
+
+
+// standard includes
+#include <vector>
+#include <map>
+
+// 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
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}
--- /dev/null
+#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
--- /dev/null
+#ifndef __MESH_REGION_HXX__
+#define __MESH_REGION_HXX__
+
+#include <vector>
+
+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<MeshElement*> _elements;
+ BoundingBox _box;
+
+ };
+
+};
+
+
+#endif
--- /dev/null
+#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
--- /dev/null
+#include "TransformedTriangle.hxx"
+#include <iostream>
+#include <fstream>
+#include <cassert>
+
+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
--- /dev/null
+#ifndef __TRANSFORMED_TRIANGLE_HXX__
+#define __TRANSFORMED_TRIANGLE_HXX__
+
+#include <vector>
+
+#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<double*> _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
