]> SALOME platform Git repositories - modules/med.git/blob - medtool/src/INTERP_KERNELTest/MeshTestToolkit.txx
Salome HOME
Move medtool folder to MED base repository
[modules/med.git] / medtool / src / INTERP_KERNELTest / MeshTestToolkit.txx
1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D
2 //
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, or (at your option) any later version.
7 //
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.
12 //
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
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 #include "TestInterpKernelUtils.hxx"
20
21 #include "MeshTestToolkit.hxx"
22
23 #include "MEDFileMesh.hxx"
24
25 #include "MEDCouplingNormalizedUnstructuredMesh.hxx"
26 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
27 #include "MEDCouplingFieldDouble.hxx"
28
29 #include "Interpolation3DSurf.hxx"
30 #include "Interpolation2D.txx"
31 #include "Interpolation3D.txx"
32
33 #include <map>
34 #include <cmath>
35 #include <vector>
36 #include <cstring>
37 #include <iostream>
38 #include <algorithm>
39
40
41 // levels :
42 // 1 - titles and volume results
43 // 2 - symmetry / diagonal results and intersection matrix output
44 // 3 - empty
45 // 4 - empty
46 // 5 - misc
47 #include "Log.hxx"
48
49 #include <cppunit/extensions/HelperMacros.h>
50
51 //#define VOL_PREC 1.0e-6
52 using namespace ParaMEDMEM;
53 using namespace INTERP_KERNEL;
54
55 namespace INTERP_TEST
56 {
57   /**
58    * Calculates the sum of a row of an intersection matrix
59    *
60    * @param m  an intersection matrix
61    * @param i  the index of the row (1 <= i <= no. rows)
62    * @return the sum of the values of row i
63    *
64    */
65   template <int SPACEDIM, int MESHDIM>
66   double MeshTestToolkit<SPACEDIM,MESHDIM>::sumRow(const IntersectionMatrix& m, int i) const
67   {
68     double vol = 0.0;
69     for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
70       {
71         if(iter->count(i) != 0.0)
72           {
73             std::map<int, double>::const_iterator iter2 = iter->find(i);
74             vol += fabs(iter2->second);
75           }
76       }
77     return vol;
78   }
79
80   /**
81    * Calculates the sum of a column of an intersection matrix
82    *
83    * @param m  an intersection matrix
84    * @param i  the index of the column (0 <= i <= no. rows - 1)
85    * @return the sum of the values of column i
86    *
87    */
88   template <int SPACEDIM, int MESHDIM>
89   double MeshTestToolkit<SPACEDIM,MESHDIM>::sumCol(const IntersectionMatrix& m, int i) const
90   {
91     double vol = 0.0;
92     const std::map<int, double>& col = m[i];
93     for(std::map<int, double>::const_iterator iter = col.begin() ; iter != col.end() ; ++iter)
94       {
95         vol += fabs(iter->second);
96       }
97     return vol;
98   }
99
100   /**
101    * Gets the volumes of the elements in a mesh.
102    *
103    * @param mesh   the mesh
104    * @param tab    pointer to double[no. elements of mesh] array in which to store the volumes
105    */
106   template <int SPACEDIM, int MESHDIM>
107   void MeshTestToolkit<SPACEDIM,MESHDIM>::getVolumes(ParaMEDMEM::MEDCouplingUMesh& mesh, double *tab) const
108   {
109     MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=mesh.getMeasureField(true);
110     std::copy(vol->getArray()->begin(),vol->getArray()->end(),tab);
111   }
112
113   /**
114    * Sums all the elements (volumes) of an intersection matrix
115    *
116    * @param m  the intersection matrix
117    * @return   the sum of the elements of m
118    */
119
120   template <int SPACEDIM, int MESHDIM>
121   double MeshTestToolkit<SPACEDIM,MESHDIM>::sumVolume(const IntersectionMatrix& m) const
122   {
123     std::vector<double> volumes;
124     for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
125       {
126         for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
127           {
128             volumes.push_back(fabs(iter2->second));
129           }
130       }
131
132     // sum in ascending order to avoid rounding errors
133
134     sort(volumes.begin(), volumes.end());
135     const double vol = accumulate(volumes.begin(), volumes.end(), 0.0);
136
137     return vol;
138   }
139
140   /**
141    * Verifies if for a given intersection matrix the sum of each row is equal to the volumes
142    * of the corresponding source elements and the sum of each column is equal to the volumes
143    * of the corresponding target elements. This will be true as long as the meshes correspond
144    * to the same geometry. The equalities are in the "epsilon-sense", making sure the relative
145    * error is small enough.
146    *
147    * @param  m       the intersection matrix
148    * @param  sMesh   the source mesh
149    * @param  tMesh   the target mesh
150    * @return true if the condition is verified, false if not.
151    */
152   template <int SPACEDIM, int MESHDIM>
153   bool MeshTestToolkit<SPACEDIM,MESHDIM>::testVolumes(const IntersectionMatrix& m,  ParaMEDMEM::MEDCouplingUMesh& sMesh,  ParaMEDMEM::MEDCouplingUMesh& tMesh) const
154   {
155     bool ok = true;
156
157     // source elements
158     double* sVol = new double[sMesh.getNumberOfCells()];
159     getVolumes(sMesh, sVol);
160
161     for(int i = 0; i < sMesh.getNumberOfCells(); ++i)
162       {
163         const double sum_row = sumRow(m, i);
164         if(!epsilonEqualRelative(sum_row, fabs(sVol[i]), _precision))
165           {
166             LOG(1, "Source volume inconsistent : vol of cell " << i << " = " << sVol[i] << " but the row sum is " << sum_row );
167             ok = false;
168           }
169         LOG(1, "diff = " <<sum_row - sVol[i] );
170       }
171
172     // target elements
173     double* tVol = new double[tMesh.getNumberOfCells()];
174     getVolumes(tMesh, tVol);
175     for(int i = 0; i < tMesh.getNumberOfCells(); ++i)
176       {
177         const double sum_col = sumCol(m, i);
178         if(!epsilonEqualRelative(sum_col,fabs(tVol[i]), _precision))
179           {
180             LOG(1, "Target volume inconsistent : vol of cell " << i << " = " << tVol[i] << " but the col sum is " << sum_col);
181             ok = false;
182           }
183         LOG(1, "diff = " <<sum_col - tVol[i] );
184       }
185     delete[] sVol;
186     delete[] tVol;
187
188     return ok;
189   }
190
191   /**
192    * Verifies that two intersection matrices have the necessary elements to be able to be each others' transposes.
193    *
194    * @param m1  the first intersection matrix
195    * @param m2  the second intersection matrix
196    *
197    * @return true if for each element (i,j) of m1, the element (j,i) exists in m2, false if not.
198    */
199   template <int SPACEDIM, int MESHDIM>
200   bool MeshTestToolkit<SPACEDIM,MESHDIM>::areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
201   {
202     bool compatitable = true;
203     int i = 0;
204     for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
205       {
206         for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
207           {
208             int j = iter2->first;
209             if(m2.at(j).count(i) == 0)
210               {
211                 if(!epsilonEqual(iter2->second, 0.0, _precision))
212                   {
213                     LOG(2, "V1( " << i << ", " << j << ") exists, but V2( " << j << ", " << i << ") " << " does not " );
214                     LOG(2, "(" << i << ", " << j << ") fails");
215                     compatitable = false;
216                   }
217               }
218           }
219         ++i;
220       }
221     if(!compatitable)
222       {
223         LOG(1, "*** matrices are not compatitable");
224       }
225     return compatitable;
226   }
227
228   /**
229    * Tests if two intersection matrices are each others' transposes.
230    *
231    * @param m1  the first intersection matrix
232    * @param m2  the second intersection matrix
233    * @return true if m1 = m2^T, false if not.
234    */
235   template <int SPACEDIM, int MESHDIM>
236   bool MeshTestToolkit<SPACEDIM,MESHDIM>::testTranspose(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
237   {
238     int i = 0;
239     bool isSymmetric = true;
240
241     LOG(1, "Checking symmetry src - target" );
242     isSymmetric = isSymmetric & areCompatitable(m1, m2) ;
243     LOG(1, "Checking symmetry target - src" );
244     isSymmetric = isSymmetric & areCompatitable(m2, m1);
245
246     for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
247       {
248         for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
249           {
250             int j = iter2->first;
251             const double v1 = fabs(iter2->second);
252             //if(m2[j - 1].count(i+1) > 0)
253             //  {
254             std::map<int, double> theMap =  m2.at(j);
255             const double v2 = fabs(theMap[i]);
256             if(v1 != v2)
257               {
258                 LOG(2, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j << ", " << i << ") = " << v2 << " | diff = " << v1 - v2 );
259                 if(!epsilonEqualRelative(v1, v2, _precision))
260                   {
261                     LOG(2, "(" << i << ", " << j << ") fails");
262                     isSymmetric = false;
263                   }
264               }
265           }
266         ++i;
267       }
268     if(!isSymmetric)
269       {
270         LOG(1, "*** matrices are not symmetric");
271       }
272     return isSymmetric;
273   }
274
275   /**
276    * Tests if an intersection matrix is diagonal.
277    *
278    * @param m   the intersection matrix
279    * @return true if m is diagonal; false if not
280    *
281    */
282   template <int SPACEDIM, int MESHDIM>
283   bool MeshTestToolkit<SPACEDIM,MESHDIM>::testDiagonal(const IntersectionMatrix& m) const
284   {
285     LOG(1, "Checking if matrix is diagonal" );
286     int i = 0;
287     bool isDiagonal = true;
288     for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
289       {
290         for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
291           {
292             int j = iter2->first;
293             const double vol = iter2->second;
294             if(vol != 0.0 && (i != j))
295               {
296                 LOG(2, "V( " << i << ", " << j << ") = " << vol << " which is not zero" );
297                 if(!epsilonEqual(vol, 0.0, _precision))
298                   {
299                     LOG(2, "(" << i << ", " << j << ") fails");
300                     isDiagonal = false;
301                   }
302               }
303           }
304         ++i;
305       }
306     if(!isDiagonal)
307       {
308         LOG(1, "*** matrix is not diagonal");
309       }
310     return isDiagonal;
311   }
312
313   /**
314    * Outputs the intersection matrix as a list of all its elements to std::cout.
315    *
316    * @param m the intersection matrix to output
317    */
318   template <int SPACEDIM, int MESHDIM>
319   void MeshTestToolkit<SPACEDIM,MESHDIM>::dumpIntersectionMatrix(const IntersectionMatrix& m) const
320   {
321     int i = 0;
322     std::cout << "Intersection matrix is " << std::endl;
323     for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
324       {
325         for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
326           {
327             std::cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << std::endl;
328           }
329         ++i;
330       }
331     std::cout << "Sum of volumes = " << sumVolume(m) << std::endl;
332   }
333
334   /**
335    * Calculates the intersection matrix for two meshes.
336    * If the source and target meshes are the same, a CppUnit assertion raised if testVolumes() returns false.
337    *
338    * @param  mesh1path   the path to the file containing the source mesh, relative to {$MEDTOOL_ROOT_DIR}/share/resources/med/
339    * @param  mesh1       the name of the source mesh
340    * @param  mesh2path   the path to the file containing the target mesh, relative to {$MEDTOOL_ROOT_DIR}/share/resources/med/
341    * @param  mesh2       the name of the target mesh
342    * @param  m           intersection matrix in which to store the result of the intersection
343    */
344   template <int SPACEDIM, int MESHDIM>
345   void MeshTestToolkit<SPACEDIM,MESHDIM>::calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const
346   {
347     LOG(1, std::endl << "=== -> intersecting src = " << mesh1path << ", target = " << mesh2path );
348
349     LOG(5, "Loading " << mesh1 << " from " << mesh1path);
350     MEDCouplingAutoRefCountObjectPtr<MEDFileUMesh> sMeshML=MEDFileUMesh::New(INTERP_TEST::getResourceFile(mesh1path).c_str(),mesh1);
351     MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> sMesh=sMeshML->getMeshAtLevel(0);
352
353     LOG(5, "Loading " << mesh2 << " from " << mesh2path);
354     MEDCouplingAutoRefCountObjectPtr<MEDFileUMesh> tMeshML=MEDFileUMesh::New(INTERP_TEST::getResourceFile(mesh2path).c_str(),mesh2);
355     MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> tMesh=tMeshML->getMeshAtLevel(0);
356
357     MEDCouplingNormalizedUnstructuredMesh<SPACEDIM,MESHDIM> sMesh_wrapper(sMesh);
358     MEDCouplingNormalizedUnstructuredMesh<SPACEDIM,MESHDIM> tMesh_wrapper(tMesh);
359
360     if (SPACEDIM==2 && MESHDIM==2)
361       {
362         Interpolation2D interpolator;
363         interpolator.setOptions(_precision, LOG_LEVEL, _intersectionType,1);
364         interpolator.interpolateMeshes(sMesh_wrapper, tMesh_wrapper,m,"P0P0");
365       }
366     else if (SPACEDIM==3 && MESHDIM==2)
367       {
368         Interpolation3DSurf interpolator;
369         interpolator.setOptions(_precision,LOG_LEVEL, 0.5,_intersectionType,false,1);
370         interpolator.interpolateMeshes(sMesh_wrapper, tMesh_wrapper,m,"P0P0");
371       }
372     else if (SPACEDIM==3 && MESHDIM==3)
373       {
374         Interpolation3D interpolator;
375         interpolator.interpolateMeshes(sMesh_wrapper, tMesh_wrapper,m,"P0P0");
376       }
377     else
378       {
379         throw INTERP_KERNEL::Exception("Wrong dimensions");
380       }
381     // if reflexive, check volumes
382     if(strcmp(mesh1path,mesh2path) == 0)
383       {
384         const bool row_and_col_sums_ok = testVolumes(m, *sMesh, *tMesh);
385         CPPUNIT_ASSERT_EQUAL_MESSAGE("Row or column sums incorrect", true, row_and_col_sums_ok);
386         const bool is_diagonal =testDiagonal(m);
387         CPPUNIT_ASSERT_EQUAL_MESSAGE("Self intersection matrix is not diagonal", true, is_diagonal);
388       }
389
390     LOG(1, "Intersection calculation done. " << std::endl );
391   }
392
393   /**
394    * Tests the intersection algorithm for two meshes.
395    * Depending on the nature of the meshes, different tests will be performed. The sum of the elements will
396    * be compared to the given total volume of the intersection in all cases. If the two meshes are the same, then
397    * it will be confirmed that the intersection matrix is diagonal, otherwise the intersection matrices will be
398    * calculated once which each mesh as source mesh, and it will be verified that the they are each others' transpose.
399    *
400    * @param  mesh1path   the path to the file containing the source mesh, relative to {$MEDTOOL_ROOT_DIR}/share/resources/med/
401    * @param  mesh1       the name of the source mesh
402    * @param  mesh2path   the path to the file containing the target mesh, relative to {$MEDTOOL_ROOT_DIR}/share/resources/med/
403    * @param  mesh2       the name of the target mesh
404    * @param  correctVol  the total volume of the intersection of the two meshes
405    * @param  prec        maximum relative error to be tolerated in volume comparisions
406    * @param  doubleTest  if false, only the test with mesh 1 as the source mesh and mesh 2 as the target mesh will be performed
407    *
408    */
409   template <int SPACEDIM, int MESHDIM>
410   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
411   {
412     LOG(1, std::endl << std::endl << "=============================" );
413
414     using std::string;
415     const string path1 = string(mesh1path) + string(mesh1);
416     const string path2 = string(mesh2path) + string(mesh2);
417
418     const bool isTestReflexive = (path1.compare(path2) == 0);
419
420     IntersectionMatrix matrix1;
421     calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
422
423 #if LOG_LEVEL >= 2
424     dumpIntersectionMatrix(matrix1);
425 #endif
426
427     std::cout.precision(16);
428
429     const double vol1 = sumVolume(matrix1);
430
431     if(!doubleTest)
432       {
433         LOG(1, "vol =  " << vol1 <<"  correctVol = " << correctVol );
434         CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1));
435
436         if(isTestReflexive)
437           {
438             CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
439           }
440       }
441     else
442       {
443         IntersectionMatrix matrix2;
444         calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);
445
446 #if LOG_LEVEL >= 2
447         dumpIntersectionMatrix(matrix2);
448 #endif
449
450         const double vol2 = sumVolume(matrix2);
451
452         LOG(1, "vol1 =  " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol );
453
454         CPPUNIT_ASSERT_EQUAL_MESSAGE("Symmetry test failed", true, testTranspose(matrix1, matrix2));
455         CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(vol1, correctVol));
456         CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec * std::max(vol2, correctVol));
457         CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec * std::max(vol1, vol2));
458       }
459   }
460
461   /**
462    * Utility method used to facilitate the call to intersect meshes.
463    * It calls intersectMeshes, using "mesh1.med" as file name for the mesh with name "mesh1" and
464    * "mesh2.med" as file name for the mesh with name "mesh2". The rest of the arguments are passed
465    * along as they are.
466    *
467    * @param  mesh1       the name of the source mesh
468    * @param  mesh2       the name of the target mesh
469    * @param  correctVol  the total volume of the intersection of the two meshes
470    * @param  prec        maximum relative error to be tolerated in volume comparisions
471    * @param  doubleTest  if false, only the test with mesh 1 as the source mesh and mesh 2 as the target mesh will be performed
472    *
473    */
474   template <int SPACEDIM, int MESHDIM>
475   void MeshTestToolkit<SPACEDIM,MESHDIM>::intersectMeshes(const char* mesh1, const char* mesh2, const double correctVol, const double prec, bool doubleTest) const
476   {
477     const std::string path1 = std::string(mesh1) + std::string(".med");
478     std::cout << "here :" << path1 << std::endl;
479     const std::string path2 = std::string(mesh2) + std::string(".med");
480
481     intersectMeshes(path1.c_str(), mesh1, path2.c_str(), mesh2, correctVol, prec, doubleTest);
482   }
483 }