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 #include "Interpolation3DTest.hxx"
21 #include "MEDMEM_Mesh.hxx"
29 #include "VectorUtils.hxx"
31 #include "MEDMEM_Field.hxx"
32 #include "MEDMEM_Support.hxx"
35 // 1 - titles and volume results
36 // 2 - symmetry / diagonal results and intersection matrix output
43 //#define VOL_PREC 1.0e-6
45 using namespace MEDMEM;
46 using namespace INTERP_KERNEL;
47 using namespace MED_EN;
49 double Interpolation3DTest::sumRow(const IntersectionMatrix& m, int i) const
52 for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
54 if(iter->count(i) != 0.0)
56 std::map<int, double>::const_iterator iter2 = iter->find(i);
63 double Interpolation3DTest::sumCol(const IntersectionMatrix& m, int i) const
66 const std::map<int, double>& col = m[i];
67 for(std::map<int, double>::const_iterator iter = col.begin() ; iter != col.end() ; ++iter)
69 vol += std::fabs(iter->second);
75 void Interpolation3DTest::getVolumes(MESH& mesh, double* tab) const
77 SUPPORT *sup=new SUPPORT(&mesh,"dummy",MED_CELL);
78 FIELD<double>* f=mesh.getVolume(sup);
79 const double *tabS=f->getValue();
80 std::copy(tabS,tabS+mesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS),tab)
84 double Interpolation3DTest::sumVolume(const IntersectionMatrix& m) const
87 std::vector<double> volumes;
88 for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
90 for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
92 volumes.push_back(iter2->second);
93 // vol += std::abs(iter2->second);
97 // sum in ascending order to avoid rounding errors
99 sort(volumes.begin(), volumes.end());
100 const double vol = accumulate(volumes.begin(), volumes.end(), 0.0);
105 bool Interpolation3DTest::testVolumes(const IntersectionMatrix& m, MESH& sMesh, MESH& tMesh) const
110 double* sVol = new double[sMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)];
111 getVolumes(sMesh, sVol);
113 for(int i = 0; i < sMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS); ++i)
115 const double sum_row = sumRow(m, i+1);
116 if(!epsilonEqualRelative(sum_row, sVol[i], VOL_PREC))
118 LOG(1, "Source volume inconsistent : vol of cell " << i << " = " << sVol[i] << " but the row sum is " << sum_row );
121 LOG(1, "diff = " <<sum_row - sVol[i] );
125 double* tVol = new double[tMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)];
126 getVolumes(tMesh, tVol);
127 for(int i = 0; i < tMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS); ++i)
129 const double sum_col = sumCol(m, i);
130 if(!epsilonEqualRelative(sum_col, tVol[i], VOL_PREC))
132 LOG(1, "Target volume inconsistent : vol of cell " << i << " = " << tVol[i] << " but the col sum is " << sum_col);
135 LOG(1, "diff = " <<sum_col - tVol[i] );
143 bool Interpolation3DTest::areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
145 bool compatitable = true;
147 for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
149 for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
151 int j = iter2->first;
152 if(m2.at(j-1).count(i+1) == 0)
154 if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
156 LOG(2, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
157 LOG(2, "(" << i << ", " << j << ") fails");
158 compatitable = false;
166 LOG(1, "*** matrices are not compatitable");
171 bool Interpolation3DTest::testSymmetric(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
175 bool isSymmetric = true;
177 LOG(1, "Checking symmetry src - target" );
178 isSymmetric = isSymmetric & areCompatitable(m1, m2) ;
179 LOG(1, "Checking symmetry target - src" );
180 isSymmetric = isSymmetric & areCompatitable(m2, m1);
182 for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
184 for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
186 int j = iter2->first;
187 const double v1 = iter2->second;
188 //if(m2[j - 1].count(i+1) > 0)
190 std::map<int, double> theMap = m2.at(j-1);
191 const double v2 = theMap[i + 1];
194 LOG(2, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 );
195 if(!epsilonEqualRelative(v1, v2, VOL_PREC))
197 LOG(2, "(" << i << ", " << j << ") fails");
206 LOG(1, "*** matrices are not symmetric");
211 bool Interpolation3DTest::testDiagonal(const IntersectionMatrix& m) const
213 LOG(1, "Checking if matrix is diagonal" );
215 bool isDiagonal = true;
216 for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
218 for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
220 int j = iter2->first;
221 const double vol = iter2->second;
222 if(vol != 0.0 && (i != j))
224 LOG(2, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" );
225 if(!epsilonEqual(vol, 0.0, VOL_PREC))
227 LOG(2, "(" << i << ", " << j << ") fails");
236 LOG(1, "*** matrix is not diagonal");
241 void Interpolation3DTest::dumpIntersectionMatrix(const IntersectionMatrix& m) const
244 std::cout << "Intersection matrix is " << std::endl;
245 for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
247 for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
250 std::cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << std::endl;
255 std::cout << "Sum of volumes = " << sumVolume(m) << std::endl;
258 void Interpolation3DTest::setUp()
260 interpolator = new Interpolation3D();
263 void Interpolation3DTest::tearDown()
268 void Interpolation3DTest::calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const
270 const string dataBaseDir = getenv("MED_ROOT_DIR");
271 const string dataDir = dataBaseDir + "/share/salome/resources/med/";
273 LOG(1, std::endl << "=== -> intersecting src = " << mesh1 << ", target = " << mesh2 );
275 LOG(5, "Loading " << mesh1 << " from " << mesh1path);
276 MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1);
278 LOG(5, "Loading " << mesh2 << " from " << mesh2path);
279 MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
281 m = interpolator->interpolateMeshes(sMesh, tMesh);
283 // if reflexive, check volumes
284 if(strcmp(mesh1path,mesh2path) == 0)
286 const bool row_and_col_sums_ok = testVolumes(m, sMesh, tMesh);
287 CPPUNIT_ASSERT_EQUAL_MESSAGE("Row or column sums incorrect", true, row_and_col_sums_ok);
290 LOG(1, "Intersection calculation done. " << std::endl );
294 void Interpolation3DTest::intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec, bool doubleTest) const
296 LOG(1, std::endl << std::endl << "=============================" );
299 const string path1 = string(mesh1path) + string(mesh1);
300 const string path2 = string(mesh2path) + string(mesh2);
302 const bool isTestReflexive = (path1.compare(path2) == 0);
304 IntersectionMatrix matrix1;
305 calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
308 dumpIntersectionMatrix(matrix1);
311 std::cout.precision(16);
313 const double vol1 = sumVolume(matrix1);
317 LOG(1, "vol = " << vol1 <<" correctVol = " << correctVol );
318 CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1));
322 CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
328 IntersectionMatrix matrix2;
329 calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);
332 dumpIntersectionMatrix(matrix2);
335 const double vol2 = sumVolume(matrix2);
337 LOG(1, "vol1 = " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol );
339 CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(vol1, correctVol));
340 CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec * std::max(vol2, correctVol));
341 CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec * std::max(vol1, vol2));
342 CPPUNIT_ASSERT_EQUAL_MESSAGE("Symmetry test failed", true, testSymmetric(matrix1, matrix2));