1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #ifndef _TESTING_UTILS_HXX_
21 #define _TESTING_UTILS_HXX_
23 #include "Interpolation3D.hxx"
24 #include "MEDMEM_Mesh.hxx"
32 #include "VectorUtils.hxx"
35 // 1 - titles and volume results
36 // 2 - symmetry / diagonal results and intersection matrix output
42 using namespace MEDMEM;
43 using namespace INTERP_KERNEL;
44 using namespace MED_EN;
47 double sumVolume(const IntersectionMatrix& m)
50 vector<double> volumes;
51 for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
53 for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
55 volumes.push_back(iter2->second);
56 // vol += std::fabs(iter2->second);
60 // sum in ascending order to avoid rounding errors
62 sort(volumes.begin(), volumes.end());
63 const double vol = accumulate(volumes.begin(), volumes.end(), 0.0);
70 bool areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2)
72 bool compatitable = true;
74 for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
76 for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
79 if(m2.at(j-1).count(i+1) == 0)
81 if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
83 LOG(2, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
84 LOG(2, "(" << i << ", " << j << ") fails");
93 LOG(1, "*** matrices are not compatitable");
98 bool testSymmetric(const IntersectionMatrix& m1, const IntersectionMatrix& m2)
102 bool isSymmetric = true;
104 LOG(1, "Checking symmetry src - target" );
105 isSymmetric = isSymmetric & areCompatitable(m1, m2) ;
106 LOG(1, "Checking symmetry target - src" );
107 isSymmetric = isSymmetric & areCompatitable(m2, m1);
109 for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
111 for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
113 int j = iter2->first;
114 const double v1 = iter2->second;
115 //if(m2[j - 1].count(i+1) > 0)
117 std::map<int, double> theMap = m2.at(j-1);
118 const double v2 = theMap[i + 1];
121 LOG(2, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 );
122 if(!epsilonEqualRelative(v1, v2, VOL_PREC))
124 LOG(2, "(" << i << ", " << j << ") fails");
133 LOG(1, "*** matrices are not symmetric");
138 bool testDiagonal(const IntersectionMatrix& m)
140 LOG(1, "Checking if matrix is diagonal" );
142 bool isDiagonal = true;
143 for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
145 for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
147 int j = iter2->first;
148 const double vol = iter2->second;
149 if(vol != 0.0 && (i != j))
151 LOG(2, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" );
152 if(!epsilonEqual(vol, 0.0, VOL_PREC))
154 LOG(2, "(" << i << ", " << j << ") fails");
163 LOG(1, "*** matrix is not diagonal");
170 void dumpIntersectionMatrix(const IntersectionMatrix& m)
173 std::cout << "Intersection matrix is " << std::endl;
174 for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
176 for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
179 std::cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << std::endl;
184 std::cout << "Sum of volumes = " << sumVolume(m) << std::endl;
187 std::pair<int,int> countNumberOfMatrixEntries(const IntersectionMatrix& m)
192 for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
194 numElems += iter->size();
195 for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
197 if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
203 return std::make_pair(numElems, numNonZero);
207 void calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m)
209 const std::string dataBaseDir = getenv("MED_ROOT_DIR");
210 const std::string dataDir = dataBaseDir + "/share/salome/resources/med/";
212 LOG(1, std::endl << "=== -> intersecting src = " << mesh1 << ", target = " << mesh2 );
214 LOG(5, "Loading " << mesh1 << " from " << mesh1path);
215 const MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1);
216 const int numSrcElems = sMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
217 LOG(1, "Source mesh has " << numSrcElems << " elements");
220 LOG(5, "Loading " << mesh2 << " from " << mesh2path);
221 const MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
222 const int numTargetElems = tMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
224 LOG(1, "Target mesh has " << numTargetElems << " elements");
226 Interpolation3D* interpolator = new Interpolation3D();
228 m = interpolator->interpolateMeshes(sMesh, tMesh);
230 std::pair<int, int> eff = countNumberOfMatrixEntries(m);
231 LOG(1, eff.first << " of " << numTargetElems * numSrcElems << " intersections calculated : ratio = "
232 << double(eff.first) / double(numTargetElems * numSrcElems));
233 LOG(1, eff.second << " non-zero elements of " << eff.first << " total : filter efficiency = "
234 << double(eff.second) / double(eff.first));
238 LOG(1, "Intersection calculation done. " << std::endl );
250 void intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec, bool doubleTest)
252 LOG(1, std::endl << std::endl << "=============================" );
255 const string path1 = string(mesh1path) + string(mesh1);
256 const string path2 = string(mesh2path) + string(mesh2);
258 const bool isTestReflexive = (path1.compare(path2) == 0);
260 IntersectionMatrix matrix1;
261 calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
264 dumpIntersectionMatrix(matrix1);
267 std::cout.precision(16);
269 const double vol1 = sumVolume(matrix1);
273 LOG(1, "vol = " << vol1 <<" correctVol = " << correctVol );
274 // CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1));
278 // CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
284 IntersectionMatrix matrix2;
285 calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);
288 dumpIntersectionMatrix(matrix2);
291 const double vol2 = sumVolume(matrix2);
293 LOG(1, "vol1 = " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol );
295 // CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(vol1, correctVol));
296 // CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec * std::max(vol2, correctVol));
297 // CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec * std::max(vol1, vol2));
298 // CPPUNIT_ASSERT_EQUAL_MESSAGE("Symmetry test failed", true, testSymmetric(matrix1, matrix2));