const int numSrcElems = srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
const int numTargetElems = targetMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
+ std::cout << "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements " << std::endl;
+
// create empty maps for all source elements
matrix.resize(numSrcElems);
{
RegionNode* currNode = nodes.top();
nodes.pop();
- // std::cout << "Popping node " << std::endl;
+ std::cout << "Popping node " << std::endl;
if(currNode->getSrcRegion().getNumberOfElements() == 1)
{
- // std::cout << " - One element" << std::endl;
+ std::cout << " - One element" << std::endl;
// volume calculation
MeshElement* srcElement = *(currNode->getSrcRegion().getBeginElements());
iter != currNode->getTargetRegion().getEndElements() ; ++iter)
{
const double vol = calculateIntersectionVolume(*srcElement, **iter);
- if(vol != 0.0)
+ // if(vol != 0.0)
{
const int targetIdx = (*iter)->getIndex();
else // recursion
{
- // std::cout << " - Recursion" << std::endl;
+ std::cout << " - Recursion" << std::endl;
RegionNode* leftNode = new RegionNode();
RegionNode* rightNode = new RegionNode();
currNode->getSrcRegion().split(leftNode->getSrcRegion(), rightNode->getSrcRegion(), axis);
+ std::cout << "After split, left src region has " << leftNode->getSrcRegion().getNumberOfElements() <<
+ " elements and right src region has " << rightNode->getSrcRegion().getNumberOfElements() << " elements" << std::endl;
+
// 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 != BoundingBox::ZMAX) ? static_cast<BoundingBox::BoxCoord>(axis + 1) : BoundingBox::XMAX;
// add target elements of curr node that overlap the two new nodes
- // std::cout << " -- Adding target elements" << std::endl;
+ std::cout << " -- Adding target elements" << std::endl;
int numLeftElements = 0;
int numRightElements = 0;
for(vector<MeshElement*>::const_iterator iter = currNode->getTargetRegion().getBeginElements() ;
iter != currNode->getTargetRegion().getEndElements() ; ++iter)
{
- // std::cout << " --- New target node" << std::endl;
+ //std::cout << " --- New target node" << std::endl;
if(!leftNode->getSrcRegion().isDisjointWithElementBoundingBox(**iter))
{
++numRightElements;
}
+
}
+ std::cout << "Left target region has " << numLeftElements << " elements and right target region has " << numRightElements << " elements" << std::endl;
+
if(numLeftElements != 0)
{
nodes.push(leftNode);
}
delete currNode;
- // std::cout << "Next iteration. Nodes left : " << nodes.size() << std::endl;
+ std::cout << "Next iteration. Nodes left : " << nodes.size() << std::endl;
}
// create AffineTransform
TetraAffineTransform T( tetraCorners );
- // T.dump();
+ std::cout << "Transform : " << std::endl;
+ T.dump();
+ std::cout << std::endl;
// triangulate source element faces (assumed tetraeder for the time being)
// do nothing
double volume = 0.0;
// std::cout << "num triangles = " << triangles.size() << std::endl;
-
+ int i = 0;
for(vector<TransformedTriangle>::iterator iter = triangles.begin() ; iter != triangles.end(); ++iter)
{
- // std::cout << std::endl << "= > Triangle " << ++i << std::endl;
+ std::cout << std::endl << "= > Triangle " << ++i << std::endl;
iter->dumpCoords();
volume += iter->calculateIntersectionVolume();
}
// typedefs
typedef std::vector< std::map< int, double > > IntersectionMatrix;
-
-#ifdef TESTING_INTERP_KERNEL
-class Interpolation3DTest;
-#endif
-
namespace INTERP_UTILS
{
class MeshElement;
class MeshRegion;
};
+class Interpolation3DTest;
+
namespace MEDMEM
{
/**
class Interpolation3D : public Interpolation
{
-#ifdef TESTING_INTERP_KERNEL
friend class ::Interpolation3DTest;
-#endif
+
public :
EXPORT_HEADERS = \
+Interpolation.hxx\
+Interpolation3D.hxx\
+RegionNode.hxx\
+MeshRegion.hxx\
+MeshElement.hxx\
+BoundingBox.hxx\
+TetraAffineTransform.hxx\
+TranslationRotationMatrix.hxx\
Interpolation2D.hxx\
Interpolation3DSurf.hxx\
-Interpolation.hxx\
InterpolationUtils.hxx\
-TranslationRotationMatrix.hxx\
BBTree.H
# Libraries targets
TransformedTriangle.cxx\
TransformedTriangle_intersect.cxx\
TransformedTriangle_math.cxx\
-Interpolation3DSurf.cxx\
-Interpolation2D.cxx\
+MeshElement.cxx\
+MeshRegion.cxx\
+RegionNode.cxx\
+BoundingBox.cxx\
+TetraAffineTransform.cxx\
+Interpolation3D.cxx
+#Interpolation3DSurf.cxx\
+#Interpolation2D.cxx\
# Executables targets
BIN =
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)
-# 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
+#LDFLAGS+= $(HDF5_LIBS) $(STDLIB) -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 $(HDF5_LIBS) -lmedmem $(BOOST_LIBS) -lutil
ifeq ($(MED_WITH_KERNEL),yes)
CPPFLAGS+= ${KERNEL_CXXFLAGS}
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]);
+ const int elemIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index] - 1;
+ return &(_mesh->getCoordinates(MED_FULL_INTERLACE)[elemIdx + 3*nodeOffset]);
}
/**
*/
void MeshElement::triangulate(std::vector<TransformedTriangle>& triangles, const TetraAffineTransform& T) const
{
+ std::cout << "Triangulating element " << getIndex() << std::endl;
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];
-
+
+ std::cout << "elements has faces at indices " << beginFaceIdx << " to " << endFaceIdx << std::endl;
+ assert(endFaceIdx - beginFaceIdx == 4);
+
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];
-
+ std::cout << "Face " << faceIdx << " with nodes in [" << beginNodeIdx << "," << endNodeIdx << "[" << std::endl;
assert(endNodeIdx - beginNodeIdx == 3);
-
+
+ double transformedPts[9];
+
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]);
+ _mesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL, MED_FACE, MED_ALL_ELEMENTS)[beginNodeIdx + j - 1] - 1;
+ const double* pt = &(_mesh->getCoordinates(MED_FULL_INTERLACE)[3*connNodeIdx]);
+
+ //const double* pt = getCoordsOfNode(j + 1);
// transform
T.apply(&transformedPts[3*j], pt);
+ std::cout << "Transforming : " << vToStr(pt) << " -> " << vToStr(&transformedPts[3*j]) << std::endl;
}
triangles.push_back(TransformedTriangle(&transformedPts[0], &transformedPts[3], &transformedPts[6]));
_box = new BoundingBox(pts, numNodes);
+ } else {
+ const int numNodes = element->getNumberNodes();
+
+ for(int i = 1 ; i <= numNodes ; ++i)
+ {
+ const double* pt = element->getCoordsOfNode(i);
+ _box->updateWithPoint(pt);
+ }
}
}
# header files
EXPORT_HEADERS = CppUnitTest.hxx \
TransformedTriangleTest.hxx \
- TransformedTriangleIntersectTest.hxx
+ TransformedTriangleIntersectTest.hxx \
+ Interpolation3DTest.hxx
# Libraries targets
BIN = TestInterpKernel
-BIN_SRC = CppUnitTest.cxx TransformedTriangleTest.cxx TransformedTriangleIntersectTest.cxx
+BIN_SRC = CppUnitTest.cxx TransformedTriangleTest.cxx TransformedTriangleIntersectTest.cxx \
+Interpolation3DTest.cxx
BIN_CLIENT_IDL =
CXXFLAGS+= -I/usr/include/cppunit
CPPFLAGS+= -I/usr/include/cppunit
-# 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
+#LDFLAGS+= $(HDF5_LIBS) $(STDLIB) -lmedmem -lutil -lmed
#LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS)
# change motivated by the bug KERNEL4778.
LDFLAGSFORBIN += $(LDFLAGS) -lm $(MED2_LIBS) $(HDF5_LIBS) \
-L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome -lmed_V2_1 -lmedmem \
-lcppunit -linterpkernel
-
+#LDFLAGSFORBIN += $(LDFLAGS) -lm $(HDF5_LIBS) \
+# -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome \
+# -lcppunit -linterpkernel
UNIT_TEST_PROG = TestInterpKernel
@CONCLUDE@
#include <math.h>
#include <iostream>
-#ifdef TESTING_INTERP_KERNEL
-class TetraAffineTransformTest;
-#endif
+#include "VectorUtils.hxx"
namespace INTERP_UTILS
{
class TetraAffineTransform
{
-#ifdef TESTING_INTERP_KERNEL
- friend class ::TetraAffineTransformTest;
-#endif
+
public:
TetraAffineTransform(const double** pts)
{
+
+ std::cout << "Creating transform from tetraeder : " << std::endl;
+ std::cout << vToStr(pts[0]) << ", " << vToStr(pts[1]) << ", " << vToStr(pts[2]) << ", " << vToStr(pts[3]) << ", " << std::endl;
+
#if 0
do {
#endif
// precalculate determinant (again after inversion of transform)
calculateDeterminant();
+
+ std::cout << "*Self-check : Applying transformation to original points ... ";
+ for(int i = 0; i < 4 ; ++i)
+ {
+ double v[3];
+ apply(v, pts[i]);
+ // std::cout << vToStr(v) << std::endl;
+ for(int j = 0; j < 3; ++j)
+ {
+ assert(epsilonEqual(v[j], (3*i+j == 3 || 3*i+j == 7 || 3*i+j == 11 ) ? 1.0 : 0.0));
+ }
+ }
+
+ std::cout << " ok" << std::endl;
}
void apply(double* destPt, const double* srcPt) const
std::cout << "A = " << std::endl << "[";
for(int i = 0; i < 3; ++i)
{
- cout << _linearTransform[3*i] << ", " << _linearTransform[3*i + 1] << ", " << _linearTransform[3*i + 1];
- if(i !=2 ) cout << endl;
+ cout << _linearTransform[3*i] << ", " << _linearTransform[3*i + 1] << ", " << _linearTransform[3*i + 2];
+ if(i != 2 ) cout << endl;
}
cout << "]" << endl;
cout << "b = " << "[" << _translation[0] << ", " << _translation[1] << ", " << _translation[2] << "]" << endl;
}
+
private:
void invertLinearTransform()
if(isTriangleBelowTetraeder())
{
+ std::cout << std::endl << "Triangle is below tetraeder - V = 0.0" << std::endl << std::endl ;
return 0.0;
}
if(sign == 0.0)
{
+ std::cout << std::endl << "Triangle is perpendicular to z-plane - V = 0.0" << std::endl << std::endl;
return 0.0;
}
{
double* ptB = new double[3];
copyVector3(&_coords[5*corner], ptB);
- _polygonA.push_back(ptB);
+ _polygonB.push_back(ptB);
std::cout << "Inclusion XYZ-plane : " << vToStr(ptB) << " added to B" << std::endl;
}
+
+ // projection on XYZ - facet
+ if(testCornerAboveXYZFacet(corner))
+ {
+ double* ptB = new double[3];
+ copyVector3(&_coords[5*corner], ptB);
+ ptB[2] = 1 - ptB[0] - ptB[1];
+ assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
+ _polygonB.push_back(ptB);
+ std::cout << "Projection XYZ-plane : " << vToStr(ptB) << " added to B" << std::endl;
+ }
+
}
+
}
/**
#include <vector>
-#ifdef TESTING_INTERP_KERNEL
+
class TransformedTriangleTest;
class TransformedTriangleIntersectTest;
class TransformedTriangleCalcVolumeTest;
-#endif
+
namespace INTERP_UTILS
{
public:
-#ifdef TESTING_INTERP_KERNEL
+
friend class ::TransformedTriangleTest;
friend class ::TransformedTriangleIntersectTest;
friend class ::TransformedTriangleCalcVolumeTest;
-#endif
+
/**
* Enumerations representing the different geometric elements of the unit tetrahedron
bool testCornerInTetrahedron(const TriCorner corner) const;
bool testCornerOnXYZFacet(const TriCorner corner) const;
-
+
+ bool testCornerAboveXYZFacet(const TriCorner corner) const;
////////////////////////////////////////////////////////////////////////////////////
/// Utility methods used in intersection tests ///////////////
bool testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const;
- bool testSegmentIntersectsH(const TriSegment seg);
+ bool testSegmentIntersectsHPlane(const TriSegment seg) const;
bool testSurfaceAboveCorner(const TetraCorner corner) const;
// coordinates of corners of tetrahedron
static const double COORDS_TET_CORNER[12];
+
+ // indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION
+ // for the calculation of the coordinates (x,y,z) of the intersection points
+ // for Segment-Facet and Segment-Edge intersections
+ static const int TransformedTriangle::DP_INDEX[12];
+
+ // correspondance edge - corners
+ static const TetraCorner CORNERS_FOR_EDGE[12];
};
#include <fstream>
#include <cassert>
#include <cmath>
+#include <limits>
namespace INTERP_UTILS
{
C_YH, C_XH, C_XY // Z
};
- const long double TransformedTriangle::MACH_EPS = 1.0e-15;
- const long double TransformedTriangle::MULT_PREC_F = 4.0*TransformedTriangle::MACH_EPS;
+ //const double TransformedTriangle::MACH_EPS = 1.0e-15;
+ const long double TransformedTriangle::MACH_EPS = std::numeric_limits<double>::epsilon();
+ 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;
+ // coordinates of corner ptTetCorner
+ /* const double TransformedTriangle::COORDS_TET_CORNER[12] =
+ {
+ 0.0, 0.0, 0.0, // O
+ 1.0, 0.0, 0.0, // X
+ 0.0, 1.0, 0.0, // Y
+ 0.0, 0.0, 1.0 // Z
+ };
+ */
////////////////////////////////////////////////////////////////////////////////////
/// Double and triple product calculations ///////////////
////////////////////////////////////////////////////////////////////////////////////
// coordinates of corner
const double ptTetCorner[3] =
{
- corner == TT::X ? 1.0 : 0.0,
- corner == TT::Y ? 1.0 : 0.0,
- corner == TT::Z ? 1.0 : 0.0
+ COORDS_TET_CORNER[3*corner ],
+ COORDS_TET_CORNER[3*corner + 1],
+ COORDS_TET_CORNER[3*corner + 2]
};
// dist^2 = ( PQ x CP )^2 / |PQ|^2 where C is the corner point
// difference vectors
const double diffPQ[3] = { ptQ[0] - ptP[0], ptQ[1] - ptP[1], ptQ[2] - ptP[2] };
const double diffCornerP[3] = { ptP[0] - ptTetCorner[0], ptP[1] - ptTetCorner[1], ptP[2] - ptTetCorner[2] };
-
+
+ //{ use functions in VectorUtils for this
+
// cross product of difference vectors
const double cross[3] =
{
for(int i = 0 ; i < 3 ; ++i) {
DoubleProduct dp = DOUBLE_PRODUCTS[3*min_corner + i];
- // std::cout << std::endl << "in code inconsistent dp :" << dp << std::endl;
-
+ // std::cout << std::endl << "inconsistent dp :" << dp << std::endl;
_doubleProducts[8*seg + dp] = 0.0;
}
const int off1 = DP_OFFSET_1[dp];
const int off2 = DP_OFFSET_2[dp];
- const double term1 = _coords[5*pt1 + off1] * _coords[5*pt2 + off2];
- const double term2 = _coords[5*pt1 + off2] * _coords[5*pt2 + off1];
- const double absTerm1 = (term1 > 0.0 ? term1 : -term1);
- const double absTerm2 = (term2 > 0.0 ? term2 : -term2);
+ const long double term1 = static_cast<long double>(_coords[5*pt1 + off1]) * static_cast<long double>(_coords[5*pt2 + off2]);
+ const long double term2 = static_cast<long double>(_coords[5*pt1 + off2]) * static_cast<long double>(_coords[5*pt2 + off1]);
- const double delta = MULT_PREC_F*(absTerm1 + absTerm2);
- const double absDoubleProduct = _doubleProducts[8*seg + dp] > 0.0 ? _doubleProducts[8*seg + dp] : -_doubleProducts[8*seg + dp];
+ const long double delta = MULT_PREC_F * ( std::abs(term1) + std::abs(term2) );
- if(absDoubleProduct < THRESHOLD_F*delta)
+ if( std::abs(static_cast<long double>(_doubleProducts[8*seg + dp])) < THRESHOLD_F*delta )
{
- //std::cout << "Double product " << 8*seg+dp << " = " << absDoubleProduct << " is imprecise, reset to 0.0" << std::endl;
+ if(_doubleProducts[8*seg + dp] != 0.0)
+ {
+ std::cout << "Double product for (seg,dp) = (" << seg << ", " << dp << ") = " << std::abs(_doubleProducts[8*seg + dp]) << " is imprecise, reset to 0.0" << std::endl;
+ }
_doubleProducts[8*seg + dp] = 0.0;
+
}
}
}
if(_isTripleProductsCalculated)
return;
+ std::cout << "Precalculating triple products" << std::endl;
for(TetraCorner corner = O ; corner <= Z ; corner = TetraCorner(corner + 1))
{
+ std::cout << "- Triple product for corner " << corner << std::endl;
bool isGoodRowFound = false;
// find edge / row to use
// find normal to PQR - cross PQ and PR
const double pq[3] =
{
- _coords[5*Q] - _coords[5*P],
+ _coords[5*Q] - _coords[5*P],
_coords[5*Q + 1] - _coords[5*P + 1],
_coords[5*Q + 2] - _coords[5*P + 2]
};
{
_tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, false);
}
+ _validTP[corner] = true;
}
else
{
// this value will not be used
// we set it to whatever
+ // std::cout << "Triple product not calculated for corner " << corner << std::endl;
_tripleProducts[corner] = -3.14159265;
+ _validTP[corner] = false;
+
}
}
double TransformedTriangle::calcStableT(const TetraCorner corner) const
{
assert(_isTripleProductsCalculated);
+ // assert(_validTP[corner]);
return _tripleProducts[corner];
}
#include <string>
#include <sstream>
-#include <math.h>
+#include <cmath>
+#include <numeric>
namespace INTERP_UTILS
{
- void copyVector3(const double* src, double* dest)
+ inline void copyVector3(const double* src, double* dest)
{
for(int i = 0 ; i < 3 ; ++i)
{
}
}
- const std::string vToStr(const double* pt)
+ inline const std::string vToStr(const double* pt)
{
std::stringstream ss(std::ios::out);
ss << "[" << pt[0] << ", " << pt[1] << ", " << pt[2] << "]";
return ss.str();
}
- void cross(const double* v1, const double* v2,double* res)
+ inline void cross(const double* v1, const double* v2,double* res)
{
res[0] = v1[1]*v2[2] - v1[2]*v2[1];
res[1] = v1[2]*v2[0] - v1[0]*v2[2];
res[2] = v1[0]*v2[1] - v1[1]*v2[0];
}
- double dot(const double* v1, const double* v2)
+ inline double dot(const double* v1, const double* v2)
{
return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
}
- double norm(const double* v)
+ inline double norm(const double* v)
{
return sqrt(dot(v,v));
}
- double angleBetweenVectors(const double* v1, const double* v2, const double* n)
+ inline double angleBetweenVectors(const double* v1, const double* v2, const double* n)
{
const double denominator = dot(v1, v2);
double v3[3];
}
+ inline bool epsilonEqual(const double x, const double y, const double errTol = 1.0e-12)
+ {
+ return std::abs(x - y) < errTol;
+ }
};
