1 // Copyright (C) 2007-2008 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
19 #include "MEDNormalizedUnstructuredMesh.hxx"
20 #include "MEDNormalizedUnstructuredMesh.txx"
22 #include "MeshTestToolkit.hxx"
23 #include "MEDMEM_Mesh.hxx"
25 #include "Interpolation3DSurf.txx"
26 #include "Interpolation2D.txx"
27 #include "Interpolation3D.txx"
35 //#include "VectorUtils.hxx"
37 #include "MEDMEM_Field.hxx"
38 #include "MEDMEM_Support.hxx"
41 // 1 - titles and volume results
42 // 2 - symmetry / diagonal results and intersection matrix output
48 #include <cppunit/extensions/HelperMacros.h>
50 //#define VOL_PREC 1.0e-6
52 using namespace MEDMEM;
54 using namespace MED_EN;
55 using namespace INTERP_KERNEL;
61 * Calculates the sum of a row of an intersection matrix
63 * @param m an intersection matrix
64 * @param i the index of the row (1 <= i <= no. rows)
65 * @return the sum of the values of row i
68 template <int SPACEDIM, int MESHDIM>
69 double MeshTestToolkit<SPACEDIM,MESHDIM>::sumRow(const IntersectionMatrix& m, int i) const
72 for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
74 if(iter->count(i) != 0.0)
76 map<int, double>::const_iterator iter2 = iter->find(i);
77 vol += fabs(iter2->second);
84 * Calculates the sum of a column of an intersection matrix
86 * @param m an intersection matrix
87 * @param i the index of the column (0 <= i <= no. rows - 1)
88 * @return the sum of the values of column i
91 template <int SPACEDIM, int MESHDIM>
92 double MeshTestToolkit<SPACEDIM,MESHDIM>::sumCol(const IntersectionMatrix& m, int i) const
95 const std::map<int, double>& col = m[i];
96 for(map<int, double>::const_iterator iter = col.begin() ; iter != col.end() ; ++iter)
98 vol += fabs(iter->second);
104 * Gets the volumes of the elements in a mesh.
106 * @param mesh the mesh
107 * @param tab pointer to double[no. elements of mesh] array in which to store the volumes
109 template <int SPACEDIM, int MESHDIM>
110 void MeshTestToolkit<SPACEDIM,MESHDIM>::getVolumes(MEDMEM::MESH& mesh, double* tab) const
112 SUPPORT *sup=new SUPPORT(&mesh,"dummy",MED_CELL);
120 f=mesh.getVolume(sup);
123 const double *tabS = f->getValue();
124 std::copy(tabS,tabS+mesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS),tab);
130 * Sums all the elements (volumes) of an intersection matrix
132 * @param m the intersection matrix
133 * @return the sum of the elements of m
136 template <int SPACEDIM, int MESHDIM>
137 double MeshTestToolkit<SPACEDIM,MESHDIM>::sumVolume(const IntersectionMatrix& m) const
140 vector<double> volumes;
141 for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
143 for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
145 volumes.push_back(fabs(iter2->second));
149 // sum in ascending order to avoid rounding errors
151 sort(volumes.begin(), volumes.end());
152 const double vol = accumulate(volumes.begin(), volumes.end(), 0.0);
158 * Verifies if for a given intersection matrix the sum of each row is equal to the volumes
159 * of the corresponding source elements and the sum of each column is equal to the volumes
160 * of the corresponding target elements. This will be true as long as the meshes correspond
161 * to the same geometry. The equalities are in the "epsilon-sense", making sure the relative
162 * error is small enough.
164 * @param m the intersection matrix
165 * @param sMesh the source mesh
166 * @param tMesh the target mesh
167 * @return true if the condition is verified, false if not.
169 template <int SPACEDIM, int MESHDIM>
170 bool MeshTestToolkit<SPACEDIM,MESHDIM>::testVolumes(const IntersectionMatrix& m, MEDMEM::MESH& sMesh, MEDMEM::MESH& tMesh) const
175 double* sVol = new double[sMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)];
176 getVolumes(sMesh, sVol);
178 for(int i = 0; i < sMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS); ++i)
180 const double sum_row = sumRow(m, i+1);
181 if(!epsilonEqualRelative(sum_row, fabs(sVol[i]), _precision))
183 LOG(1, "Source volume inconsistent : vol of cell " << i << " = " << sVol[i] << " but the row sum is " << sum_row );
186 LOG(1, "diff = " <<sum_row - sVol[i] );
190 double* tVol = new double[tMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)];
191 getVolumes(tMesh, tVol);
192 for(int i = 0; i < tMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS); ++i)
194 const double sum_col = sumCol(m, i);
195 if(!epsilonEqualRelative(sum_col,fabs(tVol[i]), _precision))
197 LOG(1, "Target volume inconsistent : vol of cell " << i << " = " << tVol[i] << " but the col sum is " << sum_col);
200 LOG(1, "diff = " <<sum_col - tVol[i] );
209 * Verifies that two intersection matrices have the necessary elements to be able to be each others' transposes.
211 * @param m1 the first intersection matrix
212 * @param m2 the second intersection matrix
214 * @return true if for each element (i,j) of m1, the element (j,i) exists in m2, false if not.
216 template <int SPACEDIM, int MESHDIM>
217 bool MeshTestToolkit<SPACEDIM,MESHDIM>::areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
219 bool compatitable = true;
221 for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
223 for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
225 int j = iter2->first;
226 if(m2.at(j-1).count(i+1) == 0)
228 if(!epsilonEqual(iter2->second, 0.0, _precision))
230 LOG(2, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
231 LOG(2, "(" << i << ", " << j << ") fails");
232 compatitable = false;
240 LOG(1, "*** matrices are not compatitable");
246 * Tests if two intersection matrices are each others' transposes.
248 * @param m1 the first intersection matrix
249 * @param m2 the second intersection matrix
250 * @return true if m1 = m2^T, false if not.
252 template <int SPACEDIM, int MESHDIM>
253 bool MeshTestToolkit<SPACEDIM,MESHDIM>::testTranspose(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
257 bool isSymmetric = true;
259 LOG(1, "Checking symmetry src - target" );
260 isSymmetric = isSymmetric & areCompatitable(m1, m2) ;
261 LOG(1, "Checking symmetry target - src" );
262 isSymmetric = isSymmetric & areCompatitable(m2, m1);
264 for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
266 for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
268 int j = iter2->first;
269 const double v1 = fabs(iter2->second);
270 //if(m2[j - 1].count(i+1) > 0)
272 map<int, double> theMap = m2.at(j-1);
273 const double v2 = fabs(theMap[i + 1]);
276 LOG(2, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 );
277 if(!epsilonEqualRelative(v1, v2, _precision))
279 LOG(2, "(" << i << ", " << j << ") fails");
288 LOG(1, "*** matrices are not symmetric");
294 * Tests if an intersection matrix is diagonal.
296 * @param m the intersection matrix
297 * @return true if m is diagonal; false if not
300 template <int SPACEDIM, int MESHDIM>
301 bool MeshTestToolkit<SPACEDIM,MESHDIM>::testDiagonal(const IntersectionMatrix& m) const
303 LOG(1, "Checking if matrix is diagonal" );
305 bool isDiagonal = true;
306 for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
308 for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
310 int j = iter2->first;
311 const double vol = iter2->second;
312 if(vol != 0.0 && (i != j))
314 LOG(2, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" );
315 if(!epsilonEqual(vol, 0.0, _precision))
317 LOG(2, "(" << i << ", " << j << ") fails");
326 LOG(1, "*** matrix is not diagonal");
332 * Outputs the intersection matrix as a list of all its elements to std::cout.
334 * @param m the intersection matrix to output
336 template <int SPACEDIM, int MESHDIM>
337 void MeshTestToolkit<SPACEDIM,MESHDIM>::dumpIntersectionMatrix(const IntersectionMatrix& m) const
340 std::cout << "Intersection matrix is " << endl;
341 for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
343 for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
346 std::cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << endl;
351 std::cout << "Sum of volumes = " << sumVolume(m) << std::endl;
355 * Calculates the intersection matrix for two meshes.
356 * If the source and target meshes are the same, a CppUnit assertion raised if testVolumes() returns false.
358 * @param mesh1path the path to the file containing the source mesh, relative to {$MED_ROOT_DIR}/share/salome/resources/med/
359 * @param mesh1 the name of the source mesh
360 * @param mesh2path the path to the file containing the target mesh, relative to {$MED_ROOT_DIR}/share/salome/resources/med/
361 * @param mesh2 the name of the target mesh
362 * @param m intersection matrix in which to store the result of the intersection
364 template <int SPACEDIM, int MESHDIM>
365 void MeshTestToolkit<SPACEDIM,MESHDIM>::calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const
367 const string dataBaseDir = getenv("MED_ROOT_DIR");
368 const string dataDir = dataBaseDir + string("/share/salome/resources/med/");
370 LOG(1, std::endl << "=== -> intersecting src = " << mesh1path << ", target = " << mesh2path );
372 LOG(5, "Loading " << mesh1 << " from " << mesh1path);
373 MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1);
375 LOG(5, "Loading " << mesh2 << " from " << mesh2path);
376 MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
378 MEDNormalizedUnstructuredMesh<SPACEDIM,MESHDIM> sMesh_wrapper(&sMesh);
379 MEDNormalizedUnstructuredMesh<SPACEDIM,MESHDIM> tMesh_wrapper(&tMesh);
381 if (SPACEDIM==2 && MESHDIM==2)
383 Interpolation2D interpolator;
384 interpolator.setOptions(_precision, LOG_LEVEL, _intersectionType,1);
385 interpolator.interpolateMeshes(sMesh_wrapper, tMesh_wrapper,m,"P0P0");
387 else if (SPACEDIM==3 && MESHDIM==2)
389 Interpolation3DSurf interpolator;
390 interpolator.setOptions(_precision,LOG_LEVEL, 0.5,_intersectionType,false,1);
391 interpolator.interpolateMeshes(sMesh_wrapper, tMesh_wrapper,m,"P0P0");
393 else if (SPACEDIM==3 && MESHDIM==3)
395 Interpolation3D interpolator;
396 interpolator.interpolateMeshes(sMesh_wrapper, tMesh_wrapper,m,"P0P0");
400 throw MEDEXCEPTION("Wrong dimensions");
402 // if reflexive, check volumes
403 if(strcmp(mesh1path,mesh2path) == 0)
405 const bool row_and_col_sums_ok = testVolumes(m, sMesh, tMesh);
406 CPPUNIT_ASSERT_EQUAL_MESSAGE("Row or column sums incorrect", true, row_and_col_sums_ok);
407 const bool is_diagonal =testDiagonal(m);
408 CPPUNIT_ASSERT_EQUAL_MESSAGE("Self intersection matrix is not diagonal", true, is_diagonal);
411 LOG(1, "Intersection calculation done. " << std::endl );
416 * Tests the intersection algorithm for two meshes.
417 * Depending on the nature of the meshes, different tests will be performed. The sum of the elements will
418 * be compared to the given total volume of the intersection in all cases. If the two meshes are the same, then
419 * it will be confirmed that the intersection matrix is diagonal, otherwise the intersection matrices will be
420 * calculated once which each mesh as source mesh, and it will be verified that the they are each others' transpose.
422 * @param mesh1path the path to the file containing the source mesh, relative to {$MED_ROOT_DIR}/share/salome/resources/med/
423 * @param mesh1 the name of the source mesh
424 * @param mesh2path the path to the file containing the target mesh, relative to {$MED_ROOT_DIR}/share/salome/resources/med/
425 * @param mesh2 the name of the target mesh
426 * @param correctVol the total volume of the intersection of the two meshes
427 * @param prec maximum relative error to be tolerated in volume comparisions
428 * @param doubleTest if false, only the test with mesh 1 as the source mesh and mesh 2 as the target mesh will be performed
431 template <int SPACEDIM, int MESHDIM>
432 void MeshTestToolkit<SPACEDIM,MESHDIM>::intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec, bool doubleTest) const
434 LOG(1, std::endl << std::endl << "=============================" );
437 const string path1 = string(mesh1path) + string(mesh1);
438 const string path2 = string(mesh2path) + string(mesh2);
440 const bool isTestReflexive = (path1.compare(path2) == 0);
442 IntersectionMatrix matrix1;
443 calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
446 dumpIntersectionMatrix(matrix1);
449 std::cout.precision(16);
451 const double vol1 = sumVolume(matrix1);
455 LOG(1, "vol = " << vol1 <<" correctVol = " << correctVol );
456 CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1));
460 CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
466 IntersectionMatrix matrix2;
467 calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);
470 dumpIntersectionMatrix(matrix2);
473 const double vol2 = sumVolume(matrix2);
475 LOG(1, "vol1 = " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol );
477 CPPUNIT_ASSERT_EQUAL_MESSAGE("Symmetry test failed", true, testTranspose(matrix1, matrix2));
478 CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(vol1, correctVol));
479 CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec * std::max(vol2, correctVol));
480 CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec * std::max(vol1, vol2));
486 * Utility method used to facilitate the call to intersect meshes.
487 * It calls intersectMeshes, using "mesh1.med" as file name for the mesh with name "mesh1" and
488 * "mesh2.med" as file name for the mesh with name "mesh2". The rest of the arguments are passed
491 * @param mesh1 the name of the source mesh
492 * @param mesh2 the name of the target mesh
493 * @param correctVol the total volume of the intersection of the two meshes
494 * @param prec maximum relative error to be tolerated in volume comparisions
495 * @param doubleTest if false, only the test with mesh 1 as the source mesh and mesh 2 as the target mesh will be performed
498 template <int SPACEDIM, int MESHDIM>
499 void MeshTestToolkit<SPACEDIM,MESHDIM>::intersectMeshes(const char* mesh1, const char* mesh2, const double correctVol, const double prec, bool doubleTest) const
501 const string path1 = string(mesh1) + string(".med");
502 std::cout << "here :" << path1 << std::endl;
503 const string path2 = string(mesh2) + string(".med");
505 intersectMeshes(path1.c_str(), mesh1, path2.c_str(), mesh2, correctVol, prec, doubleTest);