#include "VectorUtils.hxx"
+#define VOL_PREC 1.0e-11
+
using namespace MEDMEM;
using namespace std;
using namespace INTERP_UTILS;
-double Interpolation3DTest::sumVolume(IntersectionMatrix m)
+double Interpolation3DTest::sumVolume(const IntersectionMatrix& m) const
{
double vol = 0.0;
for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
return vol;
}
-bool Interpolation3DTest::isReflexive(IntersectionMatrix m1, IntersectionMatrix m2)
+bool Interpolation3DTest::areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
{
+ bool compatitable = true;
int i = 0;
- bool isReflexive = true;
for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
{
for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
{
int j = iter2->first;
- const double v1 = iter2->second;
- if(m2[j - 1].count(i+1) > 0)
+ if(m2.at(j-1).count(i+1) == 0)
{
- const double v2 = m2[j - 1][i + 1];
- if(!epsilonEqual(v1, v2))
+ if(!epsilonEqual(iter2->second, 0.0))
{
- std::cout << "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 << std::endl;
- isReflexive = false;
+ std::cout << "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " << std::endl;
+ compatitable = false;
}
}
- else
+ }
+ ++i;
+ }
+ return compatitable;
+}
+
+bool Interpolation3DTest::testSymmetric(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
+{
+
+ int i = 0;
+ bool isSymmetric = true;
+
+ std::cout << "Checking symmetry src - target" << std::endl;
+ isSymmetric = isSymmetric & areCompatitable(m1, m2) ;
+ std::cout << "Checking symmetry target - src" << std::endl;
+ isSymmetric = isSymmetric & areCompatitable(m2, m1);
+
+ for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
+ {
+ for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
+ {
+ int j = iter2->first;
+ const double v1 = iter2->second;
+ //if(m2[j - 1].count(i+1) > 0)
+ // {
+ map<int, double> theMap = m2.at(j-1);
+ const double v2 = theMap[i + 1];
+ if(v1 != v2)
{
- if(!epsilonEqual(v1, 0.0))
+ 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))
{
- std::cout << "V2( " << iter2->first - 1 << ", " << i + 1 << ") " << " does not exist" << std::endl;
- isReflexive = false;
+ isSymmetric = false;
}
}
}
++i;
}
- return isReflexive;
+ return isSymmetric;
}
- //? this is not a good test
-bool Interpolation3DTest::isIntersectionConsistent(IntersectionMatrix m)
+bool Interpolation3DTest::testDiagonal(const IntersectionMatrix& m) const
{
- bool res = true;
- int i = 0;
+ std::cout << "Checking if matrix is diagonal" << std::endl;
+ int i = 1;
+ bool isDiagonal = true;
for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
{
for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
{
- // volume should be between 0 and 1 / 6
- if(iter2->second < 0.0 - ERR_TOL || iter2->second > 1.0 / 6.0 + ERR_TOL)
+ int j = iter2->first;
+ const double vol = iter2->second;
+ if(vol != 0.0 && (i != j))
{
- cout << "Inconsistent intersection volume : V(" << i << ", " << iter2->first << ") = " << iter2->second << endl;
- res = false;
+ std::cout << "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" << std::endl;
+ if(!epsilonEqual(vol, 0.0, VOL_PREC))
+ {
+ isDiagonal = false;
+ }
}
}
++i;
}
- return res;
+ return isDiagonal;
}
-void Interpolation3DTest::dumpIntersectionMatrix(IntersectionMatrix m)
+void Interpolation3DTest::dumpIntersectionMatrix(const IntersectionMatrix& m) const
{
int i = 0;
cout << "Intersection matrix is " << endl;
delete interpolator;
}
-void Interpolation3DTest::calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m)
+void Interpolation3DTest::calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const
{
const string dataDir = getenv("DATA_DIR");
std::cout << std::endl << "=== -> intersecting src = " << mesh1 << ", target = " << mesh2 << std::endl;
std::cout << "Loading " << mesh1 << " from " << mesh1path << endl;
- MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1);
+ const MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1);
std::cout << "Loading " << mesh2 << " from " << mesh2path << endl;
- MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
+ const MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
m = interpolator->interpol_maillages(sMesh, tMesh);
std::cout << "Intersection calculation done. " << std::endl << std::endl;
}
-void Interpolation3DTest::intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec)
+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;
+ using std::string;
+ const string path1 = string(mesh1path) + string(mesh1);
+ const string path2 = string(mesh2path) + string(mesh2);
+
+ const bool isTestReflexive = (path1.compare(path2) == 0);
+
IntersectionMatrix matrix1;
calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
-
dumpIntersectionMatrix(matrix1);
-
- IntersectionMatrix matrix2;
- calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);
-
- dumpIntersectionMatrix(matrix2);
-
- bool reflexive = isReflexive(matrix1, matrix2);
+ std::cout.precision(16);
const double vol1 = sumVolume(matrix1);
- const double vol2 = sumVolume(matrix2);
- std::cout.precision(8);
- std::cout << "vol1 = " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol << std::endl;
+ if(isTestReflexive)
+ {
+ std::cout << "vol = " << vol1 <<" correctVol = " << correctVol << std::endl;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec);
+ CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
+ }
+ else
+ {
+ IntersectionMatrix matrix2;
+ calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);
+ dumpIntersectionMatrix(matrix2);
+
+ const double vol2 = sumVolume(matrix2);
+
+ std::cout << "vol1 = " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol << std::endl;
+
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec);
+ CPPUNIT_ASSERT_EQUAL_MESSAGE("Symmetry test failed", true, testSymmetric(matrix1, matrix2));
+ }
- //CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, 1.0e-6);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec);
- CPPUNIT_ASSERT_EQUAL(true, reflexive);
}
CPPUNIT_TEST_SUITE( Interpolation3DTest );
// single - element
+ //#if 0
+
CPPUNIT_TEST( tetraReflexiveUnit );
+
CPPUNIT_TEST( tetraReflexiveGeneral );
+
CPPUNIT_TEST( tetraNudgedSimpler );
CPPUNIT_TEST( tetraNudged );
CPPUNIT_TEST( tetraCorner );
CPPUNIT_TEST( tetraSimpleIncluded );
CPPUNIT_TEST( tetraDegenEdge );
CPPUNIT_TEST( tetraDegenFace );
+ CPPUNIT_TEST( tetraDegenTranslatedInPlane );
CPPUNIT_TEST( tetraHalfstripOnly );
CPPUNIT_TEST( tetraHalfstripOnly2 );
CPPUNIT_TEST( tetraSimpleHalfstripOnly );
CPPUNIT_TEST( generalTetra );
+ CPPUNIT_TEST( trickyTetra1 );
+
// multi - element
CPPUNIT_TEST( tetraComplexIncluded );
CPPUNIT_TEST( dividedUnitTetraSimplerReflexive );
CPPUNIT_TEST( nudgedDividedUnitTetra );
CPPUNIT_TEST( nudgedDividedUnitTetraSimpler );
CPPUNIT_TEST( dividedGenTetra );
- //CPPUNIT_TEST( boxReflexive );
+ //#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();
intersectMeshes("meshes/HalfstripOnly.med", "HalfstripOnly", "meshes/UnitTetra.med", "UnitTetra", 0.0);
}
- void tetraHalfstripOnly2()
+ void tetraHalfstripOnly2()
{
// NB this test is not completely significant : we should also verify that
// there are triangles on the element that give a non-zero volume
- intersectMeshes("meshes/HalfstripOnly2.med", "HalfstripOnly2", "meshes/UnitTetra.med", "UnitTetra", 0.0);
+ intersectMeshes("meshes/HalfstripOnly2.med", "HalfstripOnly2", "meshes/UnitTetra.med", "UnitTetra", 0.0);
}
void tetraSimpleHalfstripOnly()
intersectMeshes("meshes/GenTetra1.med", "GenTetra1", "meshes/GenTetra2.med", "GenTetra2", 4.91393, 1.0e-5);
}
+ void trickyTetra1()
+ {
+ intersectMeshes("meshes/UnitTetra.med", "UnitTetra", "meshes/TrickyTetra1.med", "TrickyTetra1", 0.0, 1.0e-8);
+ }
+
void tetraDegenEdge()
{
intersectMeshes("meshes/UnitTetraDegenT.med", "UnitTetraDegenT", "meshes/DegenEdgeXY.med", "DegenEdgeXY", 0.0);
intersectMeshes("meshes/UnitTetraDegenT.med", "UnitTetraDegenT", "meshes/DegenFaceXYZ.med", "DegenFaceXYZ", 0.0);
}
+ void tetraDegenTranslatedInPlane()
+ {
+ intersectMeshes("meshes/UnitTetraDegenT.med", "UnitTetraDegenT", "meshes/DegenTranslatedInPlane.med", "DegenTranslatedInPlane", 0.0571667);
+ }
+
void dividedUnitTetraReflexive()
{
intersectMeshes("meshes/DividedUnitTetra.med", "DividedUnitTetra", "meshes/DividedUnitTetra.med", "DividedUnitTetra", 0.1666667);
intersectMeshes("meshes/Box3.med", "Box3", "meshes/Box3.med", "Box3", 13.9954, 1.0e-4);
}
+ void boxReflexiveModerate()
+ {
+ intersectMeshes("meshes/Box1Moderate.med", "Box1Moderate", "meshes/Box1Moderate.med", "Box1Moderate", 1.0e6, 1.0);
+ }
+
void tetraBoxes()
{
intersectMeshes("meshes/Box1.med", "Box1", "meshes/Box2.med", "Box2", 124.197, 1.0e-3);
}
-
-
-#if 0
- void tetraHalfstripOnly()
+ void moderateBoxes()
{
- std::cout << std::endl << std::endl << "=============================" << std::endl;
- IntersectionMatrix matrix;
- // we want no transformation - unit tetra is target
- calcIntersectionMatrix("meshes/HalfstripOnly.med", "HalfstripOnly", "meshes/UnitTetra.med", "UnitTetra", matrix);
-
- // check that the total volume is zero
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, sumVolume(matrix), 1.0e-6);
-
- // check that we have non-zero volumes that cancel
- bool allVolumesZero = true;
-
- for(IntersectionMatrix::const_iterator iter = matrix.begin() ; iter != matrix.end() ; ++iter)
- {
- for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
- {
- allVolumesZero = allVolumesZero && (iter2->second == 0.0);
- }
- }
-
- CPPUNIT_ASSERT_EQUAL(false, allVolumesZero);
-
+ intersectMeshes("meshes/Box1Moderate.med", "Box1Moderate", "meshes/Box2Moderate.med", "Box2Moderate", 376856, 1.0);
}
- void tetraHalfstripOnly2()
+ void moderateBoxesSmaller()
{
- std::cout << std::endl << std::endl << "=============================" << std::endl;
- IntersectionMatrix matrix;
- // we want no transformation - unit tetra is target
- calcIntersectionMatrix("meshes/HalfstripOnly2.med", "HalfstripOnly2", "meshes/UnitTetra.med", "UnitTetra", matrix);
-
- // check that the total volume is zero
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, sumVolume(matrix), 1.0e-6);
-
- // check that we have non-zero volumes that cancel
- bool allVolumesZero = true;
-
- for(IntersectionMatrix::const_iterator iter = matrix.begin() ; iter != matrix.end() ; ++iter)
- {
- for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
- {
- allVolumesZero = allVolumesZero && (iter2->second == 0.0);
- }
- }
-
- CPPUNIT_ASSERT_EQUAL(false, allVolumesZero);
-
+ intersectMeshes("meshes/BoxModSmall1.med", "BoxModSmall1", "meshes/BoxModSmall2.med", "BoxModSmall2", 321853, 1.0);
}
- void tetraSimpleHalfstripOnly()
+ void moderateBoxSmallReflexive()
{
- std::cout << std::endl << std::endl << "=============================" << std::endl;
- IntersectionMatrix matrix;
- // we want no transformation - unit tetra is target
- calcIntersectionMatrix("meshes/SimpleHalfstripOnly.med", "SimpleHalfstripOnly", "meshes/UnitTetra.med", "UnitTetra", matrix);
-
- // check that the total volume is zero
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, sumVolume(matrix), 1.0e-6);
-
- // check that we have non-zero volumes that cancel
- bool allVolumesZero = true;
-
- for(IntersectionMatrix::const_iterator iter = matrix.begin() ; iter != matrix.end() ; ++iter)
- {
- for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
- {
- allVolumesZero = allVolumesZero && (iter2->second == 0.0);
- }
- }
- CPPUNIT_ASSERT_EQUAL(false, allVolumesZero);
-
+ intersectMeshes("meshes/BoxModSmall1.med", "BoxModSmall1", "meshes/BoxModSmall1.med", "BoxModSmall1", 1.44018e6, 1.0);
}
-#endif
- void reflexiveTetra();
- void tetraTetraScale();
-
- void cyl1();
-
- void box1();
-
- void tetra1();
+ void moderateBoxEvenSmallerReflexive()
+ {
+ intersectMeshes("meshes/BoxEvenSmaller1.med", "BoxEvenSmaller1", "meshes/BoxEvenSmaller1.med", "BoxEvenSmaller1", 1.44018e6, 1.0);
+ }
- void tetra3();
+ void tinyBoxReflexive()
+ {
+ intersectMeshes("meshes/TinyBox.med", "TinyBox", "meshes/TinyBox.med", "TinyBox", 979200, 1.0);
+ }
private:
Interpolation3D* interpolator;
- double sumVolume(IntersectionMatrix m);
+ double sumVolume(const IntersectionMatrix& m) const;
- bool Interpolation3DTest::isReflexive(IntersectionMatrix m1, IntersectionMatrix m2);
-
- bool isIntersectionConsistent(IntersectionMatrix m);
+ bool areCompatitable( const IntersectionMatrix& m1, const IntersectionMatrix& m2) const;
+
+ bool testSymmetric(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const;
+
+ bool testDiagonal(const IntersectionMatrix& m) const;
- void dumpIntersectionMatrix(IntersectionMatrix m);
+ void dumpIntersectionMatrix(const IntersectionMatrix& m) const;
- void intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec = 1.0e-6);
+ void intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec = 1.0e-6) const;
- void calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m);
+ void calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const;
};
