src/INTERP_KERNELTest/Interpolation3DTest.cxx

   1 // Copyright (C) 2007-2021  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
  20 #include "Interpolation3DTest.hxx"
  21
  22 #include "MEDFileMesh.hxx"
  23 #include "MEDCouplingUMesh.hxx"
  24 #include "MEDCouplingFieldDouble.hxx"
  25
  26 #include <map>
  27 #include <cmath>
  28 #include <vector>
  29 #include <iostream>
  30 #include <algorithm>
  31
  32 #include "VectorUtils.hxx"
  33
  34 // levels :
  35 // 1 - titles and volume results
  36 // 2 - symmetry / diagonal results and intersection matrix output
  37 // 3 - empty
  38 // 4 - empty
  39 // 5 - misc
  40 #include "Log.hxx"
  41
  42
  43 //#define VOL_PREC 1.0e-6
  44
  45 using namespace MEDCoupling;
  46 using namespace INTERP_KERNEL;
  47
  48 double Interpolation3DTest::sumRow(const IntersectionMatrix& m, int i) const
  49 {
  50   double vol = 0.0;
  51   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
  52     {
  53       if(iter->count(i) != 0.0)
  54         {
  55           std::map<int, double>::const_iterator iter2 = iter->find(i);
  56           vol += iter2->second;
  57         }
  58     }
  59   return vol;
  60 }
  61
  62 double Interpolation3DTest::sumCol(const IntersectionMatrix& m, int i) const
  63 {
  64   double vol = 0.0;
  65   const std::map<int, double>& col = m[i];
  66   for(std::map<int, double>::const_iterator iter = col.begin() ; iter != col.end() ; ++iter)
  67     {
  68       vol += std::fabs(iter->second);
  69     }
  70   return vol;
  71 }
  72
  73
  74 void Interpolation3DTest::getVolumes(MEDCoupling::MEDCouplingUMesh& mesh, double *tab) const
  75 {
  76   MCAuto<MEDCouplingFieldDouble> vol=mesh->getMeasureField(true);
  77   std::copy(vol->getArray()->begin(),vol->getArray()->end(),tab);
  78 }
  79
  80 double Interpolation3DTest::sumVolume(const IntersectionMatrix& m) const
  81 {
  82
  83   std::vector<double> volumes;
  84   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
  85     {
  86       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
  87         {
  88           volumes.push_back(iter2->second);
  89           //    vol += std::abs(iter2->second);
  90         }
  91     }
  92
  93   // sum in ascending order to avoid rounding errors
  94
  95   sort(volumes.begin(), volumes.end());
  96   const double vol = accumulate(volumes.begin(), volumes.end(), 0.0);
  97
  98   return vol;
  99 }
 100
 101 bool Interpolation3DTest::testVolumes(const IntersectionMatrix& m,  MEDCouplingUMesh& sMesh,  MEDCouplingUMesh& tMesh) const
 102 {
 103   bool ok = true;
 104
 105   // source elements
 106   double* sVol = new double[sMesh.getNumberOfCells()];
 107   getVolumes(sMesh, sVol);
 108
 109   for(int i = 0; i < sMesh.getNumberOfCells(); ++i)
 110     {
 111       const double sum_row = sumRow(m, i+1);
 112       if(!epsilonEqualRelative(sum_row, sVol[i], VOL_PREC))
 113         {
 114           LOG(1, "Source volume inconsistent : vol of cell " << i << " = " << sVol[i] << " but the row sum is " << sum_row );
 115           ok = false;
 116         }
 117       LOG(1, "diff = " <<sum_row - sVol[i] );
 118     }
 119
 120   // target elements
 121   double* tVol = new double[tMesh.getNumberOfCells()];
 122   getVolumes(tMesh, tVol);
 123   for(int i = 0; i < tMesh.getNumberOfCells(); ++i)
 124     {
 125       const double sum_col = sumCol(m, i);
 126       if(!epsilonEqualRelative(sum_col, tVol[i], VOL_PREC))
 127         {
 128           LOG(1, "Target volume inconsistent : vol of cell " << i << " = " << tVol[i] << " but the col sum is " << sum_col);
 129           ok = false;
 130         }
 131       LOG(1, "diff = " <<sum_col - tVol[i] );
 132     }
 133   delete[] sVol;
 134   delete[] tVol;
 135
 136   return ok;
 137 }
 138
 139 bool Interpolation3DTest::areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
 140 {
 141   bool compatitable = true;
 142   int i = 0;
 143   for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
 144     {
 145       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
 146         {
 147           int j = iter2->first;
 148           if(m2.at(j-1).count(i+1) == 0)
 149             {
 150               if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
 151                 {
 152                   LOG(2, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
 153                   LOG(2, "(" << i << ", " << j << ") fails");
 154                   compatitable = false;
 155                 }
 156             }
 157         }
 158       ++i;
 159     }
 160   if(!compatitable)
 161     {
 162       LOG(1, "*** matrices are not compatitable");
 163     }
 164   return compatitable;
 165 }
 166
 167 bool Interpolation3DTest::testSymmetric(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
 168 {
 169
 170   int i = 0;
 171   bool isSymmetric = true;
 172
 173   LOG(1, "Checking symmetry src - target" );
 174   isSymmetric = isSymmetric & areCompatitable(m1, m2) ;
 175   LOG(1, "Checking symmetry target - src" );
 176   isSymmetric = isSymmetric & areCompatitable(m2, m1);
 177
 178   for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
 179     {
 180       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
 181         {
 182           int j = iter2->first;
 183           const double v1 = iter2->second;
 184           //if(m2[j - 1].count(i+1) > 0)
 185           //  {
 186           std::map<int, double> theMap =  m2.at(j-1);
 187           const double v2 = theMap[i + 1];
 188           if(v1 != v2)
 189             {
 190               LOG(2, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 );
 191               if(!epsilonEqualRelative(v1, v2, VOL_PREC))
 192                 {
 193                   LOG(2, "(" << i << ", " << j << ") fails");
 194                   isSymmetric = false;
 195                 }
 196             }
 197         }
 198       ++i;
 199     }
 200   if(!isSymmetric)
 201     {
 202       LOG(1, "*** matrices are not symmetric");
 203     }
 204   return isSymmetric;
 205 }
 206
 207 bool Interpolation3DTest::testDiagonal(const IntersectionMatrix& m) const
 208 {
 209   LOG(1, "Checking if matrix is diagonal" );
 210   int i = 1;
 211   bool isDiagonal = true;
 212   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
 213     {
 214       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
 215         {
 216           int j = iter2->first;
 217           const double vol = iter2->second;
 218           if(vol != 0.0 && (i != j))
 219             {
 220               LOG(2, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" );
 221               if(!epsilonEqual(vol, 0.0, VOL_PREC))
 222                 {
 223                   LOG(2, "(" << i << ", " << j << ") fails");
 224                   isDiagonal = false;
 225                 }
 226             }
 227         }
 228       ++i;
 229     }
 230   if(!isDiagonal)
 231     {
 232       LOG(1, "*** matrix is not diagonal");
 233     }
 234   return isDiagonal;
 235 }
 236
 237 void Interpolation3DTest::dumpIntersectionMatrix(const IntersectionMatrix& m) const
 238 {
 239   int i = 0;
 240   std::cout << "Intersection matrix is " << std::endl;
 241   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
 242     {
 243       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
 244         {
 245
 246           std::cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << std::endl;
 247
 248         }
 249       ++i;
 250     }
 251   std::cout << "Sum of volumes = " << sumVolume(m) << std::endl;
 252 }
 253
 254 void Interpolation3DTest::setUp()
 255 {
 256   interpolator = new Interpolation3D();
 257 }
 258
 259 void Interpolation3DTest::tearDown()
 260 {
 261   delete interpolator;
 262 }
 263
 264 void Interpolation3DTest::calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const
 265 {
 266   string dataDir = "";
 267   if ( getenv("MEDCOUPLING_ROOT_DIR") ) {
 268     dataDir = getenv("MEDCOUPLING_ROOT_DIR");
 269     dataDir += "/share/resources/med/";
 270   }
 271   else {
 272     dataDir = get_current_dir_name();
 273     dataDir += "/../../resources/";
 274   }
 275
 276   LOG(1, std::endl << "=== -> intersecting src = " << mesh1 << ", target = " << mesh2 );
 277
 278   LOG(5, "Loading " << mesh1 << " from " << mesh1path);
 279   MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1);
 280
 281   LOG(5, "Loading " << mesh2 << " from " << mesh2path);
 282   MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
 283
 284   m = interpolator->interpolateMeshes(sMesh, tMesh);
 285
 286   // if reflexive, check volumes
 287   if(strcmp(mesh1path,mesh2path) == 0)
 288     {
 289       const bool row_and_col_sums_ok = testVolumes(m, sMesh, tMesh);
 290       CPPUNIT_ASSERT_EQUAL_MESSAGE("Row or column sums incorrect", true, row_and_col_sums_ok);
 291     }
 292
 293   LOG(1, "Intersection calculation done. " << std::endl );
 294
 295 }
 296
 297 void Interpolation3DTest::intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec, bool doubleTest) const
 298 {
 299   LOG(1, std::endl << std::endl << "=============================" );
 300
 301   using std::string;
 302   const string path1 = string(mesh1path) + string(mesh1);
 303   const string path2 = string(mesh2path) + string(mesh2);
 304
 305   const bool isTestReflexive = (path1.compare(path2) == 0);
 306
 307   IntersectionMatrix matrix1;
 308   calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
 309
 310 #if LOG_LEVEL >= 2
 311   dumpIntersectionMatrix(matrix1);
 312 #endif
 313
 314   std::cout.precision(16);
 315
 316   const double vol1 = sumVolume(matrix1);
 317
 318   if(!doubleTest)
 319     {
 320       LOG(1, "vol =  " << vol1 <<"  correctVol = " << correctVol );
 321       CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1));
 322
 323       if(isTestReflexive)
 324         {
 325           CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
 326         }
 327     }
 328   else
 329     {
 330
 331       IntersectionMatrix matrix2;
 332       calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);
 333
 334 #if LOG_LEVEL >= 2
 335       dumpIntersectionMatrix(matrix2);
 336 #endif
 337
 338       const double vol2 = sumVolume(matrix2);
 339
 340       LOG(1, "vol1 =  " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol );
 341
 342       CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(vol1, correctVol));
 343       CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec * std::max(vol2, correctVol));
 344       CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec * std::max(vol1, vol2));
 345       CPPUNIT_ASSERT_EQUAL_MESSAGE("Symmetry test failed", true, testSymmetric(matrix1, matrix2));
 346     }
 347
 348 }
 349
 350
 351