src/INTERP_KERNELTest/MeshTestToolkit.txx

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