#include <map>
#include <vector>
#include <cmath>
+#include <algorithm>
#include "VectorUtils.hxx"
-#define VOL_PREC 1.0e-11
+// Levels :
+// 1 - titles and volume results
+// 2 - empty
+// 3 - symmetry / diagonal results
+// 4 - intersection matrix output
+// 5 - misc
+#include "Log.hxx"
+
+
+#define VOL_PREC 1.0e-6
using namespace MEDMEM;
using namespace std;
double Interpolation3DTest::sumVolume(const IntersectionMatrix& m) const
{
- double vol = 0.0;
+
+ vector<double> volumes;
for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
{
for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
{
- vol += fabs(iter2->second);
+ volumes.push_back(iter2->second);
+ // vol += std::abs(iter2->second);
}
}
+
+ // sum in ascending order to avoid rounding errors
+
+ sort(volumes.begin(), volumes.end());
+ const double vol = accumulate(volumes.begin(), volumes.end(), 0.0);
+
return vol;
}
int j = iter2->first;
if(m2.at(j-1).count(i+1) == 0)
{
- if(!epsilonEqual(iter2->second, 0.0))
+ if(!epsilonEqualRelative(iter2->second, 0.0, VOL_PREC))
{
- std::cout << "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " << std::endl;
+ LOG(3, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
compatitable = false;
}
}
int i = 0;
bool isSymmetric = true;
- std::cout << "Checking symmetry src - target" << std::endl;
+ LOG(1, "Checking symmetry src - target" );
isSymmetric = isSymmetric & areCompatitable(m1, m2) ;
- std::cout << "Checking symmetry target - src" << std::endl;
+ LOG(1, "Checking symmetry target - src" );
isSymmetric = isSymmetric & areCompatitable(m2, m1);
for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
const double v2 = theMap[i + 1];
if(v1 != v2)
{
- std::cout << "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 << std::endl;
- if(!epsilonEqual(v1, v2, VOL_PREC))
+ LOG(3, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 );
+ if(!epsilonEqualRelative(v1, v2, VOL_PREC))
{
isSymmetric = false;
}
bool Interpolation3DTest::testDiagonal(const IntersectionMatrix& m) const
{
- std::cout << "Checking if matrix is diagonal" << std::endl;
+ LOG(1, "Checking if matrix is diagonal" );
int i = 1;
bool isDiagonal = true;
for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
const double vol = iter2->second;
if(vol != 0.0 && (i != j))
{
- std::cout << "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" << std::endl;
- if(!epsilonEqual(vol, 0.0, VOL_PREC))
+ LOG(3, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" );
+ if(!epsilonEqualRelative(vol, 0.0, VOL_PREC))
{
isDiagonal = false;
}
void Interpolation3DTest::dumpIntersectionMatrix(const IntersectionMatrix& m) const
{
int i = 0;
- cout << "Intersection matrix is " << endl;
+ std::cout << "Intersection matrix is " << endl;
for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
{
for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
{
- cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << endl;
+ std::cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << endl;
}
++i;
}
- std::cout << "Sum of volumes = " << sumVolume(m) << std::endl << std::endl;
+ std::cout << "Sum of volumes = " << sumVolume(m) << std::endl;
}
void Interpolation3DTest::setUp()
{
const string dataDir = getenv("DATA_DIR");
- std::cout << std::endl << "=== -> intersecting src = " << mesh1 << ", target = " << mesh2 << std::endl;
+ LOG(1, std::endl << "=== -> intersecting src = " << mesh1 << ", target = " << mesh2 );
- std::cout << "Loading " << mesh1 << " from " << mesh1path << endl;
+ LOG(5, "Loading " << mesh1 << " from " << mesh1path);
const MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1);
- std::cout << "Loading " << mesh2 << " from " << mesh2path << endl;
+ LOG(5, "Loading " << mesh2 << " from " << mesh2path);
const MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
m = interpolator->interpol_maillages(sMesh, tMesh);
- std::cout << "Intersection calculation done. " << std::endl << std::endl;
+ LOG(1, "Intersection calculation done. " << std::endl );
+
}
void Interpolation3DTest::intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec) const
{
- std::cout << std::endl << std::endl << "=============================" << std::endl;
+ LOG(1, std::endl << std::endl << "=============================" );
using std::string;
const string path1 = string(mesh1path) + string(mesh1);
IntersectionMatrix matrix1;
calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
- dumpIntersectionMatrix(matrix1);
+
+#ifdef LOG_ACTIVE
+ if(LOG_LEVEL >= 4 )
+ {
+ dumpIntersectionMatrix(matrix1);
+ }
+#endif
+
std::cout.precision(16);
const double vol1 = sumVolume(matrix1);
if(isTestReflexive)
{
- std::cout << "vol = " << vol1 <<" correctVol = " << correctVol << std::endl;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec);
+ LOG(1, "vol = " << vol1 <<" correctVol = " << correctVol );
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1));
CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
}
else
{
IntersectionMatrix matrix2;
calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);
- dumpIntersectionMatrix(matrix2);
+
+#ifdef LOG_ACTIVE
+ if(LOG_LEVEL >= 4 )
+ {
+ dumpIntersectionMatrix(matrix2);
+ }
+#endif
const double vol2 = sumVolume(matrix2);
- std::cout << "vol1 = " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol << std::endl;
+ LOG(1, "vol1 = " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol );
- CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(vol1, correctVol));
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec * std::max(vol2, correctVol));
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec * std::max(vol1, vol2));
CPPUNIT_ASSERT_EQUAL_MESSAGE("Symmetry test failed", true, testSymmetric(matrix1, matrix2));
}
class Interpolation3DTest : public CppUnit::TestFixture
{
- CPPUNIT_TEST_SUITE( Interpolation3DTest );
-
// single - element
- //#if 0
+ CPPUNIT_TEST_SUITE( Interpolation3DTest );
CPPUNIT_TEST( tetraReflexiveUnit );
-
CPPUNIT_TEST( tetraReflexiveGeneral );
-
CPPUNIT_TEST( tetraNudgedSimpler );
CPPUNIT_TEST( tetraNudged );
CPPUNIT_TEST( tetraCorner );
CPPUNIT_TEST( tetraHalfstripOnly2 );
CPPUNIT_TEST( tetraSimpleHalfstripOnly );
CPPUNIT_TEST( generalTetra );
-
CPPUNIT_TEST( trickyTetra1 );
- // multi - element
+ CPPUNIT_TEST( inconsistentTetra );
+
+
+ // multi - element
+
+
CPPUNIT_TEST( tetraComplexIncluded );
CPPUNIT_TEST( dividedUnitTetraSimplerReflexive );
CPPUNIT_TEST( dividedUnitTetraReflexive );
CPPUNIT_TEST( nudgedDividedUnitTetra );
CPPUNIT_TEST( nudgedDividedUnitTetraSimpler );
CPPUNIT_TEST( dividedGenTetra );
- //#endif
CPPUNIT_TEST( boxReflexive );
-
-
CPPUNIT_TEST( boxReflexiveModerate );
- //#if 0
-
CPPUNIT_TEST( tetraBoxes );
-
CPPUNIT_TEST( moderateBoxes );
-
CPPUNIT_TEST( moderateBoxesSmaller );
-
-
CPPUNIT_TEST( moderateBoxSmallReflexive );
-
-
CPPUNIT_TEST( moderateBoxEvenSmallerReflexive );
-
CPPUNIT_TEST( tinyBoxReflexive );
- //#endif
-
CPPUNIT_TEST_SUITE_END();
void tetraSimpleIncluded()
{
- intersectMeshes("meshes/SimpleIncludedTetra.med", "SimpleIncludedTetra", "meshes/SimpleIncludingTetra.med", "SimpleIncludingTetra", 17.0156, 1.0e-4);
+ intersectMeshes("meshes/SimpleIncludedTetra.med", "SimpleIncludedTetra", "meshes/SimpleIncludingTetra.med", "SimpleIncludingTetra", 17.0156);
}
void tetraComplexIncluded()
{
- intersectMeshes("meshes/ComplexIncludedTetra.med", "ComplexIncludedTetra", "meshes/ComplexIncludingTetra.med", "ComplexIncludingTetra", 17.0156, 1.0e-4);
+ intersectMeshes("meshes/ComplexIncludedTetra.med", "ComplexIncludedTetra", "meshes/ComplexIncludingTetra.med", "ComplexIncludingTetra", 17.0156);
}
void tetraHalfstripOnly()
void generalTetra()
{
- intersectMeshes("meshes/GenTetra1.med", "GenTetra1", "meshes/GenTetra2.med", "GenTetra2", 4.91393, 1.0e-5);
+ intersectMeshes("meshes/GenTetra1.med", "GenTetra1", "meshes/GenTetra2.med", "GenTetra2", 4.91393);
}
void trickyTetra1()
{
- intersectMeshes("meshes/UnitTetra.med", "UnitTetra", "meshes/TrickyTetra1.med", "TrickyTetra1", 0.0, 1.0e-8);
+ intersectMeshes("meshes/UnitTetra.med", "UnitTetra", "meshes/TrickyTetra1.med", "TrickyTetra1", 0.0);
+ }
+
+ void inconsistentTetra()
+ {
+ intersectMeshes("meshes/LargeUnitTetra.med", "LargeUnitTetra", "meshes/LargeInconsistentTetra.med", "LargeInconsistent", 7.86231e7);
}
+
void tetraDegenEdge()
{
intersectMeshes("meshes/UnitTetraDegenT.med", "UnitTetraDegenT", "meshes/DegenEdgeXY.med", "DegenEdgeXY", 0.0);
void boxReflexive()
{
- intersectMeshes("meshes/Box3.med", "Box3", "meshes/Box3.med", "Box3", 13.9954, 1.0e-4);
+ intersectMeshes("meshes/Box3.med", "Box3", "meshes/Box3.med", "Box3", 13.9954);
}
void boxReflexiveModerate()
{
- intersectMeshes("meshes/Box1Moderate.med", "Box1Moderate", "meshes/Box1Moderate.med", "Box1Moderate", 1.0e6, 1.0);
+ intersectMeshes("meshes/Box1Moderate.med", "Box1Moderate", "meshes/Box1Moderate.med", "Box1Moderate", 1.0e6);
}
void tetraBoxes()
{
- intersectMeshes("meshes/Box1.med", "Box1", "meshes/Box2.med", "Box2", 124.197, 1.0e-3);
+ intersectMeshes("meshes/Box1.med", "Box1", "meshes/Box2.med", "Box2", 124.197);
}
void moderateBoxes()
{
- intersectMeshes("meshes/Box1Moderate.med", "Box1Moderate", "meshes/Box2Moderate.med", "Box2Moderate", 376856, 1.0);
+ intersectMeshes("meshes/Box1Moderate.med", "Box1Moderate", "meshes/Box2Moderate.med", "Box2Moderate", 376856);
}
void moderateBoxesSmaller()
{
- intersectMeshes("meshes/BoxModSmall1.med", "BoxModSmall1", "meshes/BoxModSmall2.med", "BoxModSmall2", 321853, 1.0);
+ intersectMeshes("meshes/BoxModSmall1.med", "BoxModSmall1", "meshes/BoxModSmall2.med", "BoxModSmall2", 321853);
}
void moderateBoxSmallReflexive()
{
- intersectMeshes("meshes/BoxModSmall1.med", "BoxModSmall1", "meshes/BoxModSmall1.med", "BoxModSmall1", 1.44018e6, 1.0);
+ intersectMeshes("meshes/BoxModSmall1.med", "BoxModSmall1", "meshes/BoxModSmall1.med", "BoxModSmall1", 1.44018e6);
}
void moderateBoxEvenSmallerReflexive()
{
- intersectMeshes("meshes/BoxEvenSmaller1.med", "BoxEvenSmaller1", "meshes/BoxEvenSmaller1.med", "BoxEvenSmaller1", 1.44018e6, 1.0);
+ intersectMeshes("meshes/BoxEvenSmaller1.med", "BoxEvenSmaller1", "meshes/BoxEvenSmaller1.med", "BoxEvenSmaller1", 1.44018e6);
}
void tinyBoxReflexive()
{
- intersectMeshes("meshes/TinyBox.med", "TinyBox", "meshes/TinyBox.med", "TinyBox", 979200, 1.0);
+ intersectMeshes("meshes/TinyBox.med", "TinyBox", "meshes/TinyBox.med", "TinyBox", 979200);
}
private:
#include "Interpolation3DTest.hxx"
// --- Registers the fixture into the 'registry'
-
+//CPPUNIT_TEST_SUITE_REGISTRATION( Interpolation3DTestMultiElement );
CPPUNIT_TEST_SUITE_REGISTRATION( Interpolation3DTest );
CPPUNIT_TEST_SUITE_REGISTRATION( TransformedTriangleIntersectTest );
CPPUNIT_TEST_SUITE_REGISTRATION( TransformedTriangleTest );
-CPPUNIT_TEST_SUITE_REGISTRATION( TestBogusClass );
+//CPPUNIT_TEST_SUITE_REGISTRATION( TestBogusClass );
// --- generic Main program from KERNEL_SRC/src/Basics/Test
#include "TransformedTriangleIntersectTest.hxx"
#include <iostream>
+
+#include "Log.hxx"
+
////////////////////////////////////////////////////////////////////////////////////////////////////////
// Intersection tests
// Each method in this file runs all the intersection tests with some triangle. The goal is to cover all
void TransformedTriangleIntersectTest::testTriangle1()
{
+ LOG(1, "+++++++ Testing triangle 1" );
+
typedef TransformedTriangle TT;
double coords[9] =
void TransformedTriangleIntersectTest::testTriangle2()
{
-
+ LOG(1, "+++++++ Testing triangle 2" );
typedef TransformedTriangle TT;
double coords[9] =
void TransformedTriangleIntersectTest::testTriangle3()
{
-
+ LOG(1, "+++++++ Testing triangle 3" );
typedef TransformedTriangle TT;
double coords[9] =
void TransformedTriangleIntersectTest::testTriangle4()
{
-
+ LOG(1, "+++++++ Testing triangle 4" );
typedef TransformedTriangle TT;
double coords[9] =
void TransformedTriangleIntersectTest::testTriangle5()
{
- // std::cout << std::endl << "+++++++ Testing triangle 5" << std::endl;
+ LOG(1, "+++++++ Testing triangle 5" );
+
typedef TransformedTriangle TT;
double coords[9] =
void TransformedTriangleIntersectTest::testTriangle6()
{
- // std::cout << std::endl << "+++++++ Testing triangle 6" << std::endl;
-
- typedef TransformedTriangle TT;
+ LOG(1, "+++++++ Testing triangle 6" );
- double coords[9] =
+ typedef TransformedTriangle TT;
+
+ double coords[9] =
{
1.5, 0.5, 1.35, // P
0.5, -0.5, 2.1, // Q
-3.0, 3.0, -0.5 // R
};
-
+
TransformedTriangle* tri = new TransformedTriangle(&coords[0], &coords[3], &coords[6]);
// run all intersection tests and ensure that the ones
void TransformedTriangleIntersectTest::testTriangle7()
{
+ LOG(1, "+++++++ Testing triangle 7" );
typedef TransformedTriangle TT;
double coords[9] =
void TransformedTriangleIntersectTest::testTriangle8()
{
-
+ LOG(1, "+++++++ Testing triangle 8" );
typedef TransformedTriangle TT;
double coords[9] =
void TransformedTriangleIntersectTest::testTriangle9()
{
-
+ LOG(1, "+++++++ Testing triangle 9" );
typedef TransformedTriangle TT;
double coords[9] =
void TransformedTriangleIntersectTest::testTriangle10()
{
-
+ LOG(1, "+++++++ Testing triangle 10" );
typedef TransformedTriangle TT;
double coords[9] =
void TransformedTriangleIntersectTest::testTriangle11()
{
-
+ LOG(1, "+++++++ Testing triangle 11" );
typedef TransformedTriangle TT;
double coords[9] =
void TransformedTriangleIntersectTest::testTriangle12()
{
-
+ LOG(1, "+++++++ Testing triangle 12" );
typedef TransformedTriangle TT;
double coords[9] =
void TransformedTriangleIntersectTest::testTriangle13()
{
-
+ LOG(1, "+++++++ Testing triangle 13" );
typedef TransformedTriangle TT;
double coords[9] =
void TransformedTriangleIntersectTest::testTriangleX()
{
-
+ LOG(1, "+++++++ Testing triangle X" );
typedef TransformedTriangle TT;
double coords[9] =
// if(consistency != 0.0) {
// if(num_zeros == 2 || num_neg == 0 || num_neg == 3)
- if((num_zero == 1 && num_neg != 1) || num_zero == 2 || num_neg == 0 || num_neg == 3 )
+ if((num_zero == 1 && num_neg != 1) || num_zero == 2 || num_neg == 0 && num_zero !=3 || num_neg == 3 )
{
++num_cases;
CPPUNIT_ASSERT_DOUBLES_EQUAL(correct_c_vals[i], c_vals[i], ERR_TOL);
}
}
+
+#if 0
+void TransformedTriangleTest::inconsistent2()
+{
+
+ typedef TransformedTriangle::TriSegment TriSegment;
+ typedef TransformedTriangle::TetraCorner TetraCorner;
+
+ const double x = 54946.35168415842;
+ const double y = 84351.32165113264;
+ const double z = 14845.65498715654;
+
+ double coords[9] =
+ {
+ x+0.000001, y, z,
+ -x, -y, -z,
+ 0.5*x, 0.25*y, -10.0
+ };
+
+ std::cout << "test : " << coords[0]*coords[4] - coords[1]*coords[3] << std::endl;
+ std::cout << "test : " << coords[0]*coords[5] - coords[2]*coords[3] << std::endl;
+ std::cout << "test : " << coords[1]*coords[5] - coords[2]*coords[4] << std::endl;
+
+ TransformedTriangle* tri = new TransformedTriangle(&coords[0], &coords[3], &coords[6]);
+
+ for(TriSegment seg = TransformedTriangle::PQ ; seg <= TransformedTriangle::RP ; seg = TriSegment(seg + 1))
+ {
+ std::cout << "seg " << seg << " consistent? = " << tri->areDoubleProductsConsistent(seg) << std::endl;
+ }
+
+ tri->calculateIntersectionVolume();
+
+ delete tri;
+
+}
+#endif
void test_calcUnstableT();
void test_calcStableC_Consistency();
-
+#if 0
+ void inconsistent2();
+#endif
+
double p1[3], q1[3], r1[3];
double hp1, hq1, hr1;
double Hp1, Hq1, Hr1;
