src/INTERP_KERNELTest/MeshTestToolkit.txx

   1 //  Copyright (C) 2007-2008  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 "MEDNormalizedUnstructuredMesh.hxx"
  20 #include "MEDNormalizedUnstructuredMesh.txx"
  21
  22 #include "MeshTestToolkit.hxx"
  23 #include "MEDMEM_Mesh.hxx"
  24
  25 #include "Interpolation3DSurf.txx"
  26 #include "Interpolation2D.txx"
  27 #include "Interpolation3D.txx"
  28
  29 #include <iostream>
  30 #include <map>
  31 #include <vector>
  32 #include <cmath>
  33 #include <algorithm>
  34
  35 //#include "VectorUtils.hxx"
  36
  37 #include "MEDMEM_Field.hxx"
  38 #include "MEDMEM_Support.hxx"
  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 std;
  54 using namespace MED_EN;
  55 using namespace INTERP_KERNEL;
  56
  57 namespace INTERP_TEST
  58 {
  59
  60   /**
  61    * Calculates the sum of a row of an intersection matrix
  62    *
  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
  66    *
  67    */
  68   template <int SPACEDIM, int MESHDIM>
  69   double MeshTestToolkit<SPACEDIM,MESHDIM>::sumRow(const IntersectionMatrix& m, int i) const
  70   {
  71     double vol = 0.0;
  72     for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
  73       {
  74         if(iter->count(i) != 0.0)
  75           {
  76             map<int, double>::const_iterator iter2 = iter->find(i);
  77             vol += fabs(iter2->second);
  78           }
  79       }
  80     return vol;
  81   }
  82
  83   /**
  84    * Calculates the sum of a column of an intersection matrix
  85    *
  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
  89    *
  90    */
  91   template <int SPACEDIM, int MESHDIM>
  92   double MeshTestToolkit<SPACEDIM,MESHDIM>::sumCol(const IntersectionMatrix& m, int i) const
  93   {
  94     double vol = 0.0;
  95     const std::map<int, double>& col = m[i];
  96     for(map<int, double>::const_iterator iter = col.begin() ; iter != col.end() ; ++iter)
  97       {
  98         vol += fabs(iter->second);
  99       }
 100     return vol;
 101   }
 102
 103   /**
 104    * Gets the volumes of the elements in a mesh.
 105    *
 106    * @param mesh   the mesh
 107    * @param tab    pointer to double[no. elements of mesh] array in which to store the volumes
 108    */
 109   template <int SPACEDIM, int MESHDIM>
 110   void MeshTestToolkit<SPACEDIM,MESHDIM>::getVolumes(MEDMEM::MESH& mesh, double* tab) const
 111   {
 112     SUPPORT *sup=new SUPPORT(&mesh,"dummy",MED_CELL);
 113     FIELD<double>* f;
 114     switch (MESHDIM)
 115       {
 116       case 2:
 117         f=mesh.getArea(sup);
 118         break;
 119       case 3:
 120         f=mesh.getVolume(sup);
 121         break;
 122       }
 123     const double *tabS = f->getValue();
 124     std::copy(tabS,tabS+mesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS),tab);
 125     delete sup;
 126     delete f;
 127   }
 128
 129   /**
 130    * Sums all the elements (volumes) of an intersection matrix
 131    *
 132    * @param m  the intersection matrix
 133    * @return   the sum of the elements of m
 134    */
 135
 136   template <int SPACEDIM, int MESHDIM>
 137   double MeshTestToolkit<SPACEDIM,MESHDIM>::sumVolume(const IntersectionMatrix& m) const
 138   {
 139
 140     vector<double> volumes;
 141     for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
 142       {
 143         for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
 144           {
 145             volumes.push_back(fabs(iter2->second));
 146           }
 147       }
 148
 149     // sum in ascending order to avoid rounding errors
 150
 151     sort(volumes.begin(), volumes.end());
 152     const double vol = accumulate(volumes.begin(), volumes.end(), 0.0);
 153
 154     return vol;
 155   }
 156
 157   /**
 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.
 163    *
 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.
 168    */
 169   template <int SPACEDIM, int MESHDIM>
 170   bool MeshTestToolkit<SPACEDIM,MESHDIM>::testVolumes(const IntersectionMatrix& m,  MEDMEM::MESH& sMesh,  MEDMEM::MESH& tMesh) const
 171   {
 172     bool ok = true;
 173
 174     // source elements
 175     double* sVol = new double[sMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)];
 176     getVolumes(sMesh, sVol);
 177
 178     for(int i = 0; i < sMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS); ++i)
 179       {
 180         const double sum_row = sumRow(m, i+1);
 181         if(!epsilonEqualRelative(sum_row, fabs(sVol[i]), _precision))
 182           {
 183             LOG(1, "Source volume inconsistent : vol of cell " << i << " = " << sVol[i] << " but the row sum is " << sum_row );
 184             ok = false;
 185           }
 186         LOG(1, "diff = " <<sum_row - sVol[i] );
 187       }
 188
 189     // target elements
 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)
 193       {
 194         const double sum_col = sumCol(m, i);
 195         if(!epsilonEqualRelative(sum_col,fabs(tVol[i]), _precision))
 196           {
 197             LOG(1, "Target volume inconsistent : vol of cell " << i << " = " << tVol[i] << " but the col sum is " << sum_col);
 198             ok = false;
 199           }
 200         LOG(1, "diff = " <<sum_col - tVol[i] );
 201       }
 202     delete[] sVol;
 203     delete[] tVol;
 204
 205     return ok;
 206   }
 207
 208   /**
 209    * Verifies that two intersection matrices have the necessary elements to be able to be each others' transposes.
 210    *
 211    * @param m1  the first intersection matrix
 212    * @param m2  the second intersection matrix
 213    *
 214    * @return true if for each element (i,j) of m1, the element (j,i) exists in m2, false if not.
 215    */
 216   template <int SPACEDIM, int MESHDIM>
 217   bool MeshTestToolkit<SPACEDIM,MESHDIM>::areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
 218   {
 219     bool compatitable = true;
 220     int i = 0;
 221     for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
 222       {
 223         for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
 224           {
 225             int j = iter2->first;
 226             if(m2.at(j-1).count(i+1) == 0)
 227               {
 228                 if(!epsilonEqual(iter2->second, 0.0, _precision))
 229                   {
 230                     LOG(2, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
 231                     LOG(2, "(" << i << ", " << j << ") fails");
 232                     compatitable = false;
 233                   }
 234               }
 235           }
 236         ++i;
 237       }
 238     if(!compatitable)
 239       {
 240         LOG(1, "*** matrices are not compatitable");
 241       }
 242     return compatitable;
 243   }
 244
 245   /**
 246    * Tests if two intersection matrices are each others' transposes.
 247    *
 248    * @param m1  the first intersection matrix
 249    * @param m2  the second intersection matrix
 250    * @return true if m1 = m2^T, false if not.
 251    */
 252   template <int SPACEDIM, int MESHDIM>
 253   bool MeshTestToolkit<SPACEDIM,MESHDIM>::testTranspose(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
 254   {
 255
 256     int i = 0;
 257     bool isSymmetric = true;
 258
 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);
 263
 264     for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
 265       {
 266         for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
 267           {
 268             int j = iter2->first;
 269             const double v1 = fabs(iter2->second);
 270             //if(m2[j - 1].count(i+1) > 0)
 271             //  {
 272             map<int, double> theMap =  m2.at(j-1);
 273             const double v2 = fabs(theMap[i + 1]);
 274             if(v1 != v2)
 275               {
 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))
 278                   {
 279                     LOG(2, "(" << i << ", " << j << ") fails");
 280                     isSymmetric = false;
 281                   }
 282               }
 283           }
 284         ++i;
 285       }
 286     if(!isSymmetric)
 287       {
 288         LOG(1, "*** matrices are not symmetric");
 289       }
 290     return isSymmetric;
 291   }
 292
 293   /**
 294    * Tests if an intersection matrix is diagonal.
 295    *
 296    * @param m   the intersection matrix
 297    * @return true if m is diagonal; false if not
 298    *
 299    */
 300   template <int SPACEDIM, int MESHDIM>
 301   bool MeshTestToolkit<SPACEDIM,MESHDIM>::testDiagonal(const IntersectionMatrix& m) const
 302   {
 303     LOG(1, "Checking if matrix is diagonal" );
 304     int i = 1;
 305     bool isDiagonal = true;
 306     for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
 307       {
 308         for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
 309           {
 310             int j = iter2->first;
 311             const double vol = iter2->second;
 312             if(vol != 0.0 && (i != j))
 313               {
 314                 LOG(2, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" );
 315                 if(!epsilonEqual(vol, 0.0, _precision))
 316                   {
 317                     LOG(2, "(" << i << ", " << j << ") fails");
 318                     isDiagonal = false;
 319                   }
 320               }
 321           }
 322         ++i;
 323       }
 324     if(!isDiagonal)
 325       {
 326         LOG(1, "*** matrix is not diagonal");
 327       }
 328     return isDiagonal;
 329   }
 330
 331   /**
 332    * Outputs the intersection matrix as a list of all its elements to std::cout.
 333    *
 334    * @param m the intersection matrix to output
 335    */
 336   template <int SPACEDIM, int MESHDIM>
 337   void MeshTestToolkit<SPACEDIM,MESHDIM>::dumpIntersectionMatrix(const IntersectionMatrix& m) const
 338   {
 339     int i = 0;
 340     std::cout << "Intersection matrix is " << endl;
 341     for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
 342       {
 343         for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
 344           {
 345
 346             std::cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << endl;
 347
 348           }
 349         ++i;
 350       }
 351     std::cout << "Sum of volumes = " << sumVolume(m) << std::endl;
 352   }
 353
 354   /**
 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.
 357    *
 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
 363    */
 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
 366   {
 367     const string dataBaseDir = getenv("MED_ROOT_DIR");
 368     const string dataDir = dataBaseDir + string("/share/salome/resources/med/");
 369
 370     LOG(1, std::endl << "=== -> intersecting src = " << mesh1path << ", target = " << mesh2path );
 371
 372     LOG(5, "Loading " << mesh1 << " from " << mesh1path);
 373     MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1);
 374
 375     LOG(5, "Loading " << mesh2 << " from " << mesh2path);
 376     MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
 377
 378     MEDNormalizedUnstructuredMesh<SPACEDIM,MESHDIM> sMesh_wrapper(&sMesh);
 379     MEDNormalizedUnstructuredMesh<SPACEDIM,MESHDIM> tMesh_wrapper(&tMesh);
 380
 381     if (SPACEDIM==2 && MESHDIM==2)
 382       {
 383         Interpolation2D interpolator;
 384         interpolator.setOptions(_precision, LOG_LEVEL, _intersectionType,1);
 385         interpolator.interpolateMeshes(sMesh_wrapper, tMesh_wrapper,m,"P0P0");
 386       }
 387     else if (SPACEDIM==3 && MESHDIM==2)
 388       {
 389         Interpolation3DSurf interpolator;
 390         interpolator.setOptions(_precision,LOG_LEVEL, 0.5,_intersectionType,false,1);
 391         interpolator.interpolateMeshes(sMesh_wrapper, tMesh_wrapper,m,"P0P0");
 392       }
 393     else if (SPACEDIM==3 && MESHDIM==3)
 394       {
 395         Interpolation3D interpolator;
 396         interpolator.interpolateMeshes(sMesh_wrapper, tMesh_wrapper,m,"P0P0");
 397       }
 398     else
 399       {
 400         throw MEDEXCEPTION("Wrong dimensions");
 401       }
 402     // if reflexive, check volumes
 403     if(strcmp(mesh1path,mesh2path) == 0)
 404       {
 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);
 409       }
 410
 411     LOG(1, "Intersection calculation done. " << std::endl );
 412
 413   }
 414
 415   /**
 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.
 421    *
 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
 429    *
 430    */
 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
 433   {
 434     LOG(1, std::endl << std::endl << "=============================" );
 435
 436     using std::string;
 437     const string path1 = string(mesh1path) + string(mesh1);
 438     const string path2 = string(mesh2path) + string(mesh2);
 439
 440     const bool isTestReflexive = (path1.compare(path2) == 0);
 441
 442     IntersectionMatrix matrix1;
 443     calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
 444
 445 #if LOG_LEVEL >= 2
 446     dumpIntersectionMatrix(matrix1);
 447 #endif
 448
 449     std::cout.precision(16);
 450
 451     const double vol1 = sumVolume(matrix1);
 452
 453     if(!doubleTest)
 454       {
 455         LOG(1, "vol =  " << vol1 <<"  correctVol = " << correctVol );
 456         CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1));
 457
 458         if(isTestReflexive)
 459           {
 460             CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
 461           }
 462       }
 463     else
 464       {
 465
 466         IntersectionMatrix matrix2;
 467         calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);
 468
 469 #if LOG_LEVEL >= 2
 470         dumpIntersectionMatrix(matrix2);
 471 #endif
 472
 473         const double vol2 = sumVolume(matrix2);
 474
 475         LOG(1, "vol1 =  " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol );
 476
 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));
 481       }
 482
 483   }
 484
 485   /**
 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
 489    * along as they are.
 490    *
 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
 496    *
 497    */
 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
 500   {
 501     const string path1 = string(mesh1) + string(".med");
 502     std::cout << "here :" << path1 << std::endl;
 503     const string path2 = string(mesh2) + string(".med");
 504
 505     intersectMeshes(path1.c_str(), mesh1, path2.c_str(), mesh2, correctVol, prec, doubleTest);
 506   }
 507
 508
 509 }