Back to index

salome-med  6.5.0
Interpolation3DTest.cxx
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 #include "Interpolation3DTest.hxx"
00021 #include "MEDMEM_Mesh.hxx"
00022 
00023 #include <iostream>
00024 #include <map>
00025 #include <vector>
00026 #include <cmath>
00027 #include <algorithm>
00028 
00029 #include "VectorUtils.hxx"
00030 
00031 #include "MEDMEM_Field.hxx"
00032 #include "MEDMEM_Support.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 
00043 //#define VOL_PREC 1.0e-6
00044 
00045 using namespace MEDMEM;
00046 using namespace INTERP_KERNEL;
00047 using namespace MED_EN;
00048 
00049 double Interpolation3DTest::sumRow(const IntersectionMatrix& m, int i) const
00050 {
00051   double vol = 0.0;
00052   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
00053     {
00054       if(iter->count(i) != 0.0)
00055         {
00056           std::map<int, double>::const_iterator iter2 = iter->find(i);
00057           vol += iter2->second;
00058         }
00059     }
00060   return vol;
00061 }
00062 
00063 double Interpolation3DTest::sumCol(const IntersectionMatrix& m, int i) const
00064 {
00065   double vol = 0.0;
00066   const std::map<int, double>& col = m[i];
00067   for(std::map<int, double>::const_iterator iter = col.begin() ; iter != col.end() ; ++iter)
00068     {
00069       vol += std::fabs(iter->second);
00070     }
00071   return vol;
00072 }
00073 
00074 
00075 void Interpolation3DTest::getVolumes(MESH& mesh, double* tab) const
00076 {
00077   SUPPORT *sup=new SUPPORT(&mesh,"dummy",MED_CELL);
00078   FIELD<double>* f=mesh.getVolume(sup);
00079   const double *tabS=f->getValue();
00080   std::copy(tabS,tabS+mesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS),tab)
00081     delete sup;
00082 }
00083 
00084 double Interpolation3DTest::sumVolume(const IntersectionMatrix& m) const
00085 {
00086   
00087   std::vector<double> volumes;
00088   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
00089     {
00090       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
00091         {
00092           volumes.push_back(iter2->second);
00093           //    vol += std::abs(iter2->second);
00094         }
00095     }
00096   
00097   // sum in ascending order to avoid rounding errors
00098 
00099   sort(volumes.begin(), volumes.end());
00100   const double vol = accumulate(volumes.begin(), volumes.end(), 0.0);
00101 
00102   return vol;
00103 }
00104 
00105 bool Interpolation3DTest::testVolumes(const IntersectionMatrix& m,  MESH& sMesh,  MESH& tMesh) const
00106 {
00107   bool ok = true;
00108 
00109   // source elements
00110   double* sVol = new double[sMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)];
00111   getVolumes(sMesh, sVol);
00112 
00113   for(int i = 0; i < sMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS); ++i)
00114     {
00115       const double sum_row = sumRow(m, i+1);
00116       if(!epsilonEqualRelative(sum_row, sVol[i], VOL_PREC))
00117         {
00118           LOG(1, "Source volume inconsistent : vol of cell " << i << " = " << sVol[i] << " but the row sum is " << sum_row );
00119           ok = false;
00120         }
00121       LOG(1, "diff = " <<sum_row - sVol[i] );
00122     }
00123 
00124   // target elements
00125   double* tVol = new double[tMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)];
00126   getVolumes(tMesh, tVol);
00127   for(int i = 0; i < tMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS); ++i)
00128     {
00129       const double sum_col = sumCol(m, i);
00130       if(!epsilonEqualRelative(sum_col, tVol[i], VOL_PREC))
00131         {
00132           LOG(1, "Target volume inconsistent : vol of cell " << i << " = " << tVol[i] << " but the col sum is " << sum_col);
00133           ok = false;
00134         }
00135       LOG(1, "diff = " <<sum_col - tVol[i] );
00136     }
00137   delete[] sVol;
00138   delete[] tVol;
00139 
00140   return ok;
00141 }
00142 
00143 bool Interpolation3DTest::areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
00144 {
00145   bool compatitable = true;
00146   int i = 0;
00147   for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
00148     {
00149       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
00150         {
00151           int j = iter2->first;
00152           if(m2.at(j-1).count(i+1) == 0)
00153             {
00154               if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
00155                 {
00156                   LOG(2, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
00157                   LOG(2, "(" << i << ", " << j << ") fails");
00158                   compatitable = false;
00159                 }
00160             }
00161         }
00162       ++i;
00163     }
00164   if(!compatitable)
00165     {
00166       LOG(1, "*** matrices are not compatitable");
00167     }
00168   return compatitable;
00169 }
00170       
00171 bool Interpolation3DTest::testSymmetric(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
00172 {
00173 
00174   int i = 0;
00175   bool isSymmetric = true;
00176 
00177   LOG(1, "Checking symmetry src - target" );
00178   isSymmetric = isSymmetric & areCompatitable(m1, m2) ;
00179   LOG(1, "Checking symmetry target - src" );
00180   isSymmetric = isSymmetric & areCompatitable(m2, m1);
00181 
00182   for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
00183     {
00184       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
00185         {
00186           int j = iter2->first;
00187           const double v1 = iter2->second;
00188           //if(m2[j - 1].count(i+1) > 0)
00189           //  {
00190           std::map<int, double> theMap =  m2.at(j-1);
00191           const double v2 = theMap[i + 1];
00192           if(v1 != v2)
00193             {
00194               LOG(2, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 );
00195               if(!epsilonEqualRelative(v1, v2, VOL_PREC))
00196                 {
00197                   LOG(2, "(" << i << ", " << j << ") fails");
00198                   isSymmetric = false;
00199                 }
00200             }
00201         }
00202       ++i;
00203     }
00204   if(!isSymmetric)
00205     {
00206       LOG(1, "*** matrices are not symmetric");
00207     }
00208   return isSymmetric;
00209 }
00210 
00211 bool Interpolation3DTest::testDiagonal(const IntersectionMatrix& m) const
00212 {
00213   LOG(1, "Checking if matrix is diagonal" );
00214   int i = 1;
00215   bool isDiagonal = true;
00216   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
00217     {
00218       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
00219         {
00220           int j = iter2->first;
00221           const double vol = iter2->second;
00222           if(vol != 0.0 && (i != j))
00223             {
00224               LOG(2, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" );
00225               if(!epsilonEqual(vol, 0.0, VOL_PREC))
00226                 {
00227                   LOG(2, "(" << i << ", " << j << ") fails");
00228                   isDiagonal = false;
00229                 }
00230             }
00231         }
00232       ++i;
00233     }
00234   if(!isDiagonal)
00235     {
00236       LOG(1, "*** matrix is not diagonal");
00237     }
00238   return isDiagonal;
00239 }
00240 
00241 void Interpolation3DTest::dumpIntersectionMatrix(const IntersectionMatrix& m) const
00242 {
00243   int i = 0;
00244   std::cout << "Intersection matrix is " << std::endl;
00245   for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
00246     {
00247       for(std::map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
00248         {
00249     
00250           std::cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << std::endl;
00251     
00252         }
00253       ++i;
00254     }
00255   std::cout << "Sum of volumes = " << sumVolume(m) << std::endl;
00256 }
00257 
00258 void Interpolation3DTest::setUp()
00259 {
00260   interpolator = new Interpolation3D();
00261 }
00262 
00263 void Interpolation3DTest::tearDown()
00264 {
00265   delete interpolator;
00266 } 
00267 
00268 void Interpolation3DTest::calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const
00269 {
00270   const string dataBaseDir = getenv("MED_ROOT_DIR");
00271   const string dataDir = dataBaseDir + "/share/salome/resources/med/";
00272 
00273   LOG(1, std::endl << "=== -> intersecting src = " << mesh1 << ", target = " << mesh2 );
00274 
00275   LOG(5, "Loading " << mesh1 << " from " << mesh1path);
00276   MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1);
00277   
00278   LOG(5, "Loading " << mesh2 << " from " << mesh2path);
00279   MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
00280 
00281   m = interpolator->interpolateMeshes(sMesh, tMesh);
00282 
00283   // if reflexive, check volumes
00284   if(strcmp(mesh1path,mesh2path) == 0)
00285     {
00286       const bool row_and_col_sums_ok = testVolumes(m, sMesh, tMesh);
00287       CPPUNIT_ASSERT_EQUAL_MESSAGE("Row or column sums incorrect", true, row_and_col_sums_ok);
00288     }
00289 
00290   LOG(1, "Intersection calculation done. " << std::endl );
00291   
00292 }
00293 
00294 void Interpolation3DTest::intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec, bool doubleTest) const
00295 {
00296   LOG(1, std::endl << std::endl << "=============================" );
00297 
00298   using std::string;
00299   const string path1 = string(mesh1path) + string(mesh1);
00300   const string path2 = string(mesh2path) + string(mesh2);
00301 
00302   const bool isTestReflexive = (path1.compare(path2) == 0);
00303 
00304   IntersectionMatrix matrix1;
00305   calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1);
00306 
00307 #if LOG_LEVEL >= 2 
00308   dumpIntersectionMatrix(matrix1);
00309 #endif
00310 
00311   std::cout.precision(16);
00312 
00313   const double vol1 = sumVolume(matrix1);
00314 
00315   if(!doubleTest)
00316     {
00317       LOG(1, "vol =  " << vol1 <<"  correctVol = " << correctVol );
00318       CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1));
00319 
00320       if(isTestReflexive)
00321         {
00322           CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1));
00323         }
00324     }
00325   else
00326     {
00327       
00328       IntersectionMatrix matrix2;
00329       calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2);    
00330 
00331 #if LOG_LEVEL >= 2
00332       dumpIntersectionMatrix(matrix2);
00333 #endif
00334       
00335       const double vol2 = sumVolume(matrix2);
00336 
00337       LOG(1, "vol1 =  " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol );
00338 
00339       CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(vol1, correctVol));
00340       CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec * std::max(vol2, correctVol));
00341       CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec * std::max(vol1, vol2));
00342       CPPUNIT_ASSERT_EQUAL_MESSAGE("Symmetry test failed", true, testSymmetric(matrix1, matrix2));
00343     }
00344 
00345 }
00346 
00347 
00348