Back to index

salome-med  6.5.0
TestingUtils.hxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #ifndef _TESTING_UTILS_HXX_
00021 #define _TESTING_UTILS_HXX_
00022 
00023 #include "Interpolation3D.hxx"
00024 #include "MEDMEM_Mesh.hxx"
00025 
00026 #include <iostream>
00027 #include <map>
00028 #include <vector>
00029 #include <cmath>
00030 #include <algorithm>
00031 
00032 #include "VectorUtils.hxx"
00033 
00034 // levels : 
00035 // 1 - titles and volume results
00036 // 2 - symmetry / diagonal results and intersection matrix output
00037 // 3 - empty
00038 // 4 - empty
00039 // 5 - misc
00040 #include "Log.hxx"
00041 
00042 using namespace MEDMEM;
00043 using namespace INTERP_KERNEL;
00044 using namespace MED_EN;
00045 
00046 
00047 double sumVolume(const IntersectionMatrix& m) 
00048 {
00049   
00050   vector<double> volumes;
00051   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
00052     {
00053       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
00054         {
00055           volumes.push_back(iter2->second);
00056           //    vol += std::fabs(iter2->second);
00057         }
00058     }
00059   
00060   // sum in ascending order to avoid rounding errors
00061 
00062   sort(volumes.begin(), volumes.end());
00063   const double vol = accumulate(volumes.begin(), volumes.end(), 0.0);
00064 
00065   return vol;
00066 }
00067 
00068 #if 0
00069 
00070 bool areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2)
00071 {
00072   bool compatitable = true;
00073   int i = 0;
00074   for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
00075     {
00076       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
00077         {
00078           int j = iter2->first;
00079           if(m2.at(j-1).count(i+1) == 0)
00080             {
00081               if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
00082                 {
00083                   LOG(2, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
00084                   LOG(2, "(" << i << ", " << j << ") fails");
00085                   compatitable = false;
00086                 }
00087             }
00088         }
00089       ++i;
00090     }
00091   if(!compatitable)
00092     {
00093       LOG(1, "*** matrices are not compatitable");
00094     }
00095   return compatitable;
00096 }
00097       
00098 bool testSymmetric(const IntersectionMatrix& m1, const IntersectionMatrix& m2)
00099 {
00100 
00101   int i = 0;
00102   bool isSymmetric = true;
00103 
00104   LOG(1, "Checking symmetry src - target" );
00105   isSymmetric = isSymmetric & areCompatitable(m1, m2) ;
00106   LOG(1, "Checking symmetry target - src" );
00107   isSymmetric = isSymmetric & areCompatitable(m2, m1);
00108 
00109   for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
00110     {
00111       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
00112         {
00113           int j = iter2->first;
00114           const double v1 = iter2->second;
00115           //if(m2[j - 1].count(i+1) > 0)
00116           //  {
00117           std::map<int, double> theMap =  m2.at(j-1);
00118           const double v2 = theMap[i + 1];
00119           if(v1 != v2)
00120             {
00121               LOG(2, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 );
00122               if(!epsilonEqualRelative(v1, v2, VOL_PREC))
00123                 {
00124                   LOG(2, "(" << i << ", " << j << ") fails");
00125                   isSymmetric = false;
00126                 }
00127             }
00128         }
00129       ++i;
00130     }
00131   if(!isSymmetric)
00132     {
00133       LOG(1, "*** matrices are not symmetric");
00134     }
00135   return isSymmetric;
00136 }
00137 
00138 bool testDiagonal(const IntersectionMatrix& m)
00139 {
00140   LOG(1, "Checking if matrix is diagonal" );
00141   int i = 1;
00142   bool isDiagonal = true;
00143   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
00144     {
00145       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
00146         {
00147           int j = iter2->first;
00148           const double vol = iter2->second;
00149           if(vol != 0.0 && (i != j))
00150             {
00151               LOG(2, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" );
00152               if(!epsilonEqual(vol, 0.0, VOL_PREC))
00153                 {
00154                   LOG(2, "(" << i << ", " << j << ") fails");
00155                   isDiagonal = false;
00156                 }
00157             }
00158         }
00159       ++i;
00160     }
00161   if(!isDiagonal)
00162     {
00163       LOG(1, "*** matrix is not diagonal");
00164     }
00165   return isDiagonal;
00166 }
00167 
00168 #endif
00169 
00170 void dumpIntersectionMatrix(const IntersectionMatrix& m) 
00171 {
00172   int i = 0;
00173   std::cout << "Intersection matrix is " << std::endl;
00174   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
00175     {
00176       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
00177         {
00178     
00179           std::cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << std::endl;
00180     
00181         }
00182       ++i;
00183     }
00184   std::cout << "Sum of volumes = " << sumVolume(m) << std::endl;
00185 }
00186 
00187 std::pair<int,int> countNumberOfMatrixEntries(const IntersectionMatrix& m)
00188 {
00189   
00190   int numElems = 0;
00191   int numNonZero = 0;
00192   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
00193     {
00194       numElems += iter->size();
00195       for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
00196         {
00197           if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
00198             {
00199               ++numNonZero;
00200             }
00201         }
00202     }
00203   return std::make_pair(numElems, numNonZero);
00204 }
00205 
00206 
00207 void calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) 
00208 {
00209   const std::string dataBaseDir = getenv("MED_ROOT_DIR");
00210   const std::string dataDir = dataBaseDir + "/share/salome/resources/med/";
00211 
00212   LOG(1, std::endl << "=== -> intersecting src = " << mesh1 << ", target = " << mesh2 );
00213 
00214   LOG(5, "Loading " << mesh1 << " from " << mesh1path);
00215   const MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1);
00216   const int numSrcElems = sMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
00217   LOG(1, "Source mesh has " << numSrcElems << " elements");
00218 
00219 
00220   LOG(5, "Loading " << mesh2 << " from " << mesh2path);
00221   const MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
00222   const int numTargetElems = tMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
00223 
00224   LOG(1, "Target mesh has " << numTargetElems << " elements");
00225 
00226   Interpolation3D* interpolator = new Interpolation3D();
00227 
00228   m = interpolator->interpolateMeshes(sMesh, tMesh);
00229 
00230   std::pair<int, int> eff = countNumberOfMatrixEntries(m);
00231   LOG(1, eff.first << " of " << numTargetElems * numSrcElems << " intersections calculated : ratio = " 
00232       << double(eff.first) / double(numTargetElems * numSrcElems));
00233   LOG(1, eff.second << " non-zero elements of " << eff.first << " total : filter efficiency = " 
00234       << double(eff.second) / double(eff.first));
00235 
00236   delete interpolator;
00237 
00238   LOG(1, "Intersection calculation done. " << std::endl );
00239   
00240 }
00241 
00242 
00243 
00244 
00245 
00246 
00247 
00248 
00249 #if 0
00250 void intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec, bool doubleTest) 
00251 {
00252   LOG(1, std::endl << std::endl << "=============================" );
00253 
00254   using std::string;
00255   const string path1 = string(mesh1path) + string(mesh1);
00256   const string path2 = string(mesh2path) + string(mesh2);
00257 
00258   const bool isTestReflexive = (path1.compare(path2) == 0);
00259 
00260   IntersectionMatrix matrix1;
00261   calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
00262 
00263 #if LOG_LEVEL >= 2 
00264   dumpIntersectionMatrix(matrix1);
00265 #endif
00266 
00267   std::cout.precision(16);
00268 
00269   const double vol1 = sumVolume(matrix1);
00270 
00271   if(!doubleTest)
00272     {
00273       LOG(1, "vol =  " << vol1 <<"  correctVol = " << correctVol );
00274       // CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1));
00275      
00276       if(isTestReflexive)
00277         {
00278           // CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
00279         }
00280     }
00281   else
00282     {
00283       
00284       IntersectionMatrix matrix2;
00285       calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);    
00286 
00287 #if LOG_LEVEL >= 2
00288       dumpIntersectionMatrix(matrix2);
00289 #endif
00290       
00291       const double vol2 = sumVolume(matrix2);
00292 
00293       LOG(1, "vol1 =  " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol );
00294 
00295       // CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(vol1, correctVol));
00296       // CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec * std::max(vol2, correctVol));
00297       // CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec * std::max(vol1, vol2));
00298       // CPPUNIT_ASSERT_EQUAL_MESSAGE("Symmetry test failed", true, testSymmetric(matrix1, matrix2));
00299     }
00300 
00301 }
00302 
00303 
00304 
00305 #endif
00306 
00307 
00308 #endif