src/INTERP_KERNELTest/Interpolation3DTest.cxx

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