src/INTERP_KERNELTest/TestingUtils.hxx

   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 #ifndef _TESTING_UTILS_HXX_
  21 #define _TESTING_UTILS_HXX_
  22
  23 #include "Interpolation3D.hxx"
  24 #include "MEDMEM_Mesh.hxx"
  25
  26 #include <iostream>
  27 #include <map>
  28 #include <vector>
  29 #include <cmath>
  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 using namespace MEDMEM;
  43 using namespace INTERP_KERNEL;
  44 using namespace MED_EN;
  45
  46
  47 double sumVolume(const IntersectionMatrix& m)
  48 {
  49
  50   vector<double> volumes;
  51   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
  52     {
  53       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
  54         {
  55           volumes.push_back(iter2->second);
  56           //    vol += std::fabs(iter2->second);
  57         }
  58     }
  59
  60   // sum in ascending order to avoid rounding errors
  61
  62   sort(volumes.begin(), volumes.end());
  63   const double vol = accumulate(volumes.begin(), volumes.end(), 0.0);
  64
  65   return vol;
  66 }
  67
  68 #if 0
  69
  70 bool areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2)
  71 {
  72   bool compatitable = true;
  73   int i = 0;
  74   for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
  75     {
  76       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
  77         {
  78           int j = iter2->first;
  79           if(m2.at(j-1).count(i+1) == 0)
  80             {
  81               if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
  82                 {
  83                   LOG(2, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
  84                   LOG(2, "(" << i << ", " << j << ") fails");
  85                   compatitable = false;
  86                 }
  87             }
  88         }
  89       ++i;
  90     }
  91   if(!compatitable)
  92     {
  93       LOG(1, "*** matrices are not compatitable");
  94     }
  95   return compatitable;
  96 }
  97
  98 bool testSymmetric(const IntersectionMatrix& m1, const IntersectionMatrix& m2)
  99 {
 100
 101   int i = 0;
 102   bool isSymmetric = true;
 103
 104   LOG(1, "Checking symmetry src - target" );
 105   isSymmetric = isSymmetric & areCompatitable(m1, m2) ;
 106   LOG(1, "Checking symmetry target - src" );
 107   isSymmetric = isSymmetric & areCompatitable(m2, m1);
 108
 109   for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
 110     {
 111       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
 112         {
 113           int j = iter2->first;
 114           const double v1 = iter2->second;
 115           //if(m2[j - 1].count(i+1) > 0)
 116           //  {
 117           std::map<int, double> theMap =  m2.at(j-1);
 118           const double v2 = theMap[i + 1];
 119           if(v1 != v2)
 120             {
 121               LOG(2, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 );
 122               if(!epsilonEqualRelative(v1, v2, VOL_PREC))
 123                 {
 124                   LOG(2, "(" << i << ", " << j << ") fails");
 125                   isSymmetric = false;
 126                 }
 127             }
 128         }
 129       ++i;
 130     }
 131   if(!isSymmetric)
 132     {
 133       LOG(1, "*** matrices are not symmetric");
 134     }
 135   return isSymmetric;
 136 }
 137
 138 bool testDiagonal(const IntersectionMatrix& m)
 139 {
 140   LOG(1, "Checking if matrix is diagonal" );
 141   int i = 1;
 142   bool isDiagonal = true;
 143   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.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           const double vol = iter2->second;
 149           if(vol != 0.0 && (i != j))
 150             {
 151               LOG(2, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" );
 152               if(!epsilonEqual(vol, 0.0, VOL_PREC))
 153                 {
 154                   LOG(2, "(" << i << ", " << j << ") fails");
 155                   isDiagonal = false;
 156                 }
 157             }
 158         }
 159       ++i;
 160     }
 161   if(!isDiagonal)
 162     {
 163       LOG(1, "*** matrix is not diagonal");
 164     }
 165   return isDiagonal;
 166 }
 167
 168 #endif
 169
 170 void dumpIntersectionMatrix(const IntersectionMatrix& m)
 171 {
 172   int i = 0;
 173   std::cout << "Intersection matrix is " << std::endl;
 174   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
 175     {
 176       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
 177         {
 178
 179           std::cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << std::endl;
 180
 181         }
 182       ++i;
 183     }
 184   std::cout << "Sum of volumes = " << sumVolume(m) << std::endl;
 185 }
 186
 187 std::pair<int,int> countNumberOfMatrixEntries(const IntersectionMatrix& m)
 188 {
 189
 190   int numElems = 0;
 191   int numNonZero = 0;
 192   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
 193     {
 194       numElems += iter->size();
 195       for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
 196         {
 197           if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
 198             {
 199               ++numNonZero;
 200             }
 201         }
 202     }
 203   return std::make_pair(numElems, numNonZero);
 204 }
 205
 206
 207 void calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m)
 208 {
 209   const std::string dataBaseDir = getenv("MED_ROOT_DIR");
 210   const std::string dataDir = dataBaseDir + "/share/salome/resources/med/";
 211
 212   LOG(1, std::endl << "=== -> intersecting src = " << mesh1 << ", target = " << mesh2 );
 213
 214   LOG(5, "Loading " << mesh1 << " from " << mesh1path);
 215   const MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1);
 216   const int numSrcElems = sMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
 217   LOG(1, "Source mesh has " << numSrcElems << " elements");
 218
 219
 220   LOG(5, "Loading " << mesh2 << " from " << mesh2path);
 221   const MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
 222   const int numTargetElems = tMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
 223
 224   LOG(1, "Target mesh has " << numTargetElems << " elements");
 225
 226   Interpolation3D* interpolator = new Interpolation3D();
 227
 228   m = interpolator->interpolateMeshes(sMesh, tMesh);
 229
 230   std::pair<int, int> eff = countNumberOfMatrixEntries(m);
 231   LOG(1, eff.first << " of " << numTargetElems * numSrcElems << " intersections calculated : ratio = "
 232       << double(eff.first) / double(numTargetElems * numSrcElems));
 233   LOG(1, eff.second << " non-zero elements of " << eff.first << " total : filter efficiency = "
 234       << double(eff.second) / double(eff.first));
 235
 236   delete interpolator;
 237
 238   LOG(1, "Intersection calculation done. " << std::endl );
 239
 240 }
 241
 242
 243
 244
 245
 246
 247
 248
 249 #if 0
 250 void intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec, bool doubleTest)
 251 {
 252   LOG(1, std::endl << std::endl << "=============================" );
 253
 254   using std::string;
 255   const string path1 = string(mesh1path) + string(mesh1);
 256   const string path2 = string(mesh2path) + string(mesh2);
 257
 258   const bool isTestReflexive = (path1.compare(path2) == 0);
 259
 260   IntersectionMatrix matrix1;
 261   calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
 262
 263 #if LOG_LEVEL >= 2
 264   dumpIntersectionMatrix(matrix1);
 265 #endif
 266
 267   std::cout.precision(16);
 268
 269   const double vol1 = sumVolume(matrix1);
 270
 271   if(!doubleTest)
 272     {
 273       LOG(1, "vol =  " << vol1 <<"  correctVol = " << correctVol );
 274       // CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1));
 275
 276       if(isTestReflexive)
 277         {
 278           // CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
 279         }
 280     }
 281   else
 282     {
 283
 284       IntersectionMatrix matrix2;
 285       calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);
 286
 287 #if LOG_LEVEL >= 2
 288       dumpIntersectionMatrix(matrix2);
 289 #endif
 290
 291       const double vol2 = sumVolume(matrix2);
 292
 293       LOG(1, "vol1 =  " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol );
 294
 295       // CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(vol1, correctVol));
 296       // CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec * std::max(vol2, correctVol));
 297       // CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec * std::max(vol1, vol2));
 298       // CPPUNIT_ASSERT_EQUAL_MESSAGE("Symmetry test failed", true, testSymmetric(matrix1, matrix2));
 299     }
 300
 301 }
 302
 303
 304
 305 #endif
 306
 307
 308 #endif