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 }