From: ageay <ageay>
Date: Mon, 19 Apr 2010 14:22:11 +0000 (+0000)
Subject: Correction of bug in extrudedmesh interpolation.
X-Git-Tag: V5_1_main_FINAL~121
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=88c35e80a28df39d8a26888c012dd8616ff6c0c5;p=tools%2Fmedcoupling.git

Correction of bug in extrudedmesh interpolation.
---

diff --git a/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx b/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx
index 29ef6031b..9a73e6a72 100644
--- a/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx
+++ b/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx
@@ -26,6 +26,7 @@
 #include "MEDCouplingBasicsTest.hxx"
 
 #include <cmath>
+#include <numeric>
 
 using namespace ParaMEDMEM;
 
@@ -696,6 +697,174 @@ void MEDCouplingRemapperTest::testExtruded()
   extT->decrRef();
 }
 
+void MEDCouplingRemapperTest::testExtruded2()
+{
+  MEDCouplingUMesh *meshN,*meshTT,*meshTF;
+  MEDCouplingBasicsTest::build3DExtrudedUMesh_2(meshN,meshTT,meshTF);
+  std::vector<int> n;
+  double pt[3]={300.,300.,0.};
+  double v[3]={0.,0.,2.};
+  meshN->findNodesOnPlane(pt,v,1e-12,n);
+  MEDCouplingUMesh *meshN2D=(MEDCouplingUMesh *)meshN->buildFacePartOfMySelfNode(&n[0],&n[0]+n.size(),true);
+  n.clear();
+  bool b=false;
+  DataArrayInt *da=meshTT->mergeNodes(1e-12,b);
+  CPPUNIT_ASSERT(b);
+  da->decrRef();
+  meshTT->findNodesOnPlane(pt,v,1e-12,n);
+  MEDCouplingUMesh *meshTT2D=(MEDCouplingUMesh *)meshTT->buildFacePartOfMySelfNode(&n[0],&n[0]+n.size(),true);
+  n.clear();
+  meshTF->findNodesOnPlane(pt,v,1e-12,n);
+  MEDCouplingUMesh *meshTF2D=(MEDCouplingUMesh *)meshTF->buildFacePartOfMySelfNode(&n[0],&n[0]+n.size(),true);
+  n.clear();
+  //
+  MEDCouplingExtrudedMesh *meshNE=MEDCouplingExtrudedMesh::New(meshN,meshN2D,0);
+  MEDCouplingExtrudedMesh *meshTTE=MEDCouplingExtrudedMesh::New(meshTT,meshTT2D,0);
+  MEDCouplingExtrudedMesh *meshTFE=MEDCouplingExtrudedMesh::New(meshTF,meshTF2D,0);
+  //
+  MEDCouplingRemapper remapper;
+  remapper.setPrecision(1e-12);
+  remapper.setIntersectionType(INTERP_KERNEL::Geometric2D);
+  CPPUNIT_ASSERT_EQUAL(1,remapper.prepare(meshNE,meshTTE,"P0P0"));
+  MEDCouplingFieldDouble *srcField=MEDCouplingFieldDouble::New(ON_CELLS);
+  srcField->setNature(IntegralGlobConstraint);
+  srcField->setMesh(meshNE);
+  DataArrayDouble *array=DataArrayDouble::New();
+  array->alloc(meshNE->getNumberOfCells(),1);
+  srcField->setArray(array);
+  double vals1[40]={
+    1000.,1000.,1020.,1030.,1040.,1000.,1000.,1070.,1080.,1090.,1000.,1000.,1120.,1130.,1140.,1000.,1000.,1170.,1180.,1190.,
+    2000.,2000.,2020.,2030.,2040.,2000.,2000.,2070.,2080.,2090.,2000.,2000.,2120.,2130.,2140.,2000.,2000.,2170.,2180.,2190.,
+  };
+  CPPUNIT_ASSERT_EQUAL((int)(sizeof(vals1)/sizeof(double)),meshNE->getNumberOfCells());
+  std::copy(vals1,vals1+meshNE->getNumberOfCells(),array->getPointer());
+  array->decrRef();
+  MEDCouplingFieldDouble *trgField=remapper.transferField(srcField,4.220173);
+  double expected1[200]={
+    800.,800.,800.,800.,800.,800.,800.,800.,800.,800.,1600.,1600.,1600.,1600.,1600.,1600.,1600.,1600.,1600.,1600.,
+    102.,102.,102.,102.,102.,102.,102.,102.,102.,102.,202.,202.,202.,202.,202.,202.,202.,202.,202.,202.,
+    103.,103.,103.,103.,103.,103.,103.,103.,103.,103.,203.,203.,203.,203.,203.,203.,203.,203.,203.,203.,
+    104.,104.,104.,104.,104.,104.,104.,104.,104.,104.,204.,204.,204.,204.,204.,204.,204.,204.,204.,204.,
+    219.,219.,219.,219.,219.,219.,219.,219.,219.,219.,419.,419.,419.,419.,419.,419.,419.,419.,419.,419.,
+    221.,221.,221.,221.,221.,221.,221.,221.,221.,221.,421.,421.,421.,421.,421.,421.,421.,421.,421.,421.,
+    223.,223.,223.,223.,223.,223.,223.,223.,223.,223.,423.,423.,423.,423.,423.,423.,423.,423.,423.,423.,
+    117.,117.,117.,117.,117.,117.,117.,117.,117.,117.,217.,217.,217.,217.,217.,217.,217.,217.,217.,217.,
+    118.,118.,118.,118.,118.,118.,118.,118.,118.,118.,218.,218.,218.,218.,218.,218.,218.,218.,218.,218.,
+    119.,119.,119.,119.,119.,119.,119.,119.,119.,119.,219.,219.,219.,219.,219.,219.,219.,219.,219.,219.
+    };
+  for(int i=0;i<200;i++)
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],trgField->getArray()->getConstPointer()[i],1e-3);//1e-3 precision due to non coincidence in 1D mesh
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(std::accumulate(expected1,expected1+200,0.),std::accumulate(trgField->getArray()->getConstPointer(),trgField->getArray()->getConstPointer()+200,0.),1e-10);
+  trgField->decrRef();
+  //
+  CPPUNIT_ASSERT_EQUAL(1,remapper.prepare(meshNE,meshTFE,"P0P0"));
+  trgField=remapper.transferField(srcField,4.220173);
+  double expected2[340]={25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75,
+                         160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75,
+                         29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75,
+                         26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25,
+                         172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5,
+                         51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5,
+                         85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25,
+                         29.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154.,
+                         308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75,
+                         161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5,
+                         101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25,
+                         155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25,
+                         108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75,
+                         51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25,
+                         322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 800., 800., 800., 800., 800., 800., 800., 800., 800., 800., 1600., 1600., 1600., 1600., 1600., 1600., 1600.,
+                         1600., 1600., 1600.};
+  for(int i=0;i<340;i++)
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],trgField->getArray()->getConstPointer()[i],1e-3);//1e-3 precision due to non coincidence in 1D mesh
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(std::accumulate(expected2,expected2+340,0.),std::accumulate(trgField->getArray()->getConstPointer(),trgField->getArray()->getConstPointer()+340,0.),1e-10);
+  trgField->decrRef();
+  srcField->decrRef();
+  //
+  double vals2[200]={
+    100., 200., 300., 400., 500., 600., 700., 800., 900., 1000., 1100., 1200., 1300., 1400., 1500., 1600., 1700., 1800., 1900., 2000,
+    101., 201., 301., 401., 501., 601., 701., 801., 901., 1001., 1101., 1201., 1301., 1401., 1501., 1601., 1701., 1801., 1901., 2001,
+    102., 202., 302., 402., 502., 602., 702., 802., 902., 1002., 1102., 1202., 1302., 1402., 1502., 1602., 1702., 1802., 1902., 2002,
+    103., 203., 303., 403., 503., 603., 703., 803., 903., 1003., 1103., 1203., 1303., 1403., 1503., 1603., 1703., 1803., 1903., 2003,
+    104., 204., 304., 404., 504., 604., 704., 804., 904., 1004., 1104., 1204., 1304., 1404., 1504., 1604., 1704., 1804., 1904., 2004,
+    105., 205., 305., 405., 505., 605., 705., 805., 905., 1005., 1105., 1205., 1305., 1405., 1505., 1605., 1705., 1805., 1905., 2005,
+    106., 206., 306., 406., 506., 606., 706., 806., 906., 1006., 1106., 1206., 1306., 1406., 1506., 1606., 1706., 1806., 1906., 2006,
+    107., 207., 307., 407., 507., 607., 707., 807., 907., 1007., 1107., 1207., 1307., 1407., 1507., 1607., 1707., 1807., 1907., 2007,
+    108., 208., 308., 408., 508., 608., 708., 808., 908., 1008., 1108., 1208., 1308., 1408., 1508., 1608., 1708., 1808., 1908., 2008,
+    109., 209., 309., 409., 509., 609., 709., 809., 909., 1009., 1109., 1209., 1309., 1409., 1509., 1609., 1709., 1809., 1909., 2009.
+  };
+  CPPUNIT_ASSERT_EQUAL(1,remapper.prepare(meshNE,meshTTE,"P0P0"));
+  trgField=MEDCouplingFieldDouble::New(ON_CELLS);
+  trgField->setNature(ConservativeVolumic);
+  trgField->setMesh(meshTTE);
+  array=DataArrayDouble::New();
+  array->alloc(meshTTE->getNumberOfCells(),1);
+  trgField->setArray(array);
+  std::copy(vals2,vals2+meshTTE->getNumberOfCells(),array->getPointer());
+  array->decrRef();
+  srcField=remapper.reverseTransferField(trgField,4.220173);
+  double expected3[40]={
+    550.,550.,551.,552.,553.,550.,550.,554.,555.,556.,550.,550.,554.,555.,556.,550.,550.,557.,558.,559.,
+    1550.,1550.,1551.,1552.,1553.,1550.,1550.,1554.,1555.,1556.,1550.,1550.,1554.,1555.,1556.,1550.,1550.,1557.,1558.,1559.
+  };
+  for(int i=0;i<40;i++)
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(expected3[i],srcField->getArray()->getConstPointer()[i],1e-3);//1e-3 precision due to non coincidence in 1D mesh
+  srcField->decrRef();
+  trgField->decrRef();
+  //
+  double vals3[340]={
+    100., 101., 102., 103., 104., 105., 106., 107., 108., 109., 110., 111., 112., 113., 114., 115.,
+    200., 201., 202., 203., 204., 205., 206., 207., 208., 209., 210., 211., 212., 213., 214., 215.,
+    300., 301., 302., 303., 304., 305., 306., 307., 308., 309., 310., 311., 312., 313., 314., 315.,
+    400., 401., 402., 403., 404., 405., 406., 407., 408., 409., 410., 411., 412., 413., 414., 415.,
+    500., 501., 502., 503., 504., 505., 506., 507., 508., 509., 510., 511., 512., 513., 514., 515.,
+    600., 601., 602., 603., 604., 605., 606., 607., 608., 609., 610., 611., 612., 613., 614., 615.,
+    700., 701., 702., 703., 704., 705., 706., 707., 708., 709., 710., 711., 712., 713., 714., 715.,
+    800., 801., 802., 803., 804., 805., 806., 807., 808., 809., 810., 811., 812., 813., 814., 815.,
+    900., 901., 902., 903., 904., 905., 906., 907., 908., 909., 910., 911., 912., 913., 914., 915.,
+    1000., 1001., 1002., 1003., 1004., 1005., 1006., 1007., 1008., 1009., 1010., 1011., 1012., 1013., 1014., 1015.,
+    1100., 1101., 1102., 1103., 1104., 1105., 1106., 1107., 1108., 1109., 1110., 1111., 1112., 1113., 1114., 1115.,
+    1200., 1201., 1202., 1203., 1204., 1205., 1206., 1207., 1208., 1209., 1210., 1211., 1212., 1213., 1214., 1215.,
+    1300., 1301., 1302., 1303., 1304., 1305., 1306., 1307., 1308., 1309., 1310., 1311., 1312., 1313., 1314., 1315.,
+    1400., 1401., 1402., 1403., 1404., 1405., 1406., 1407., 1408., 1409., 1410., 1411., 1412., 1413., 1414., 1415.,
+    1500., 1501., 1502., 1503., 1504., 1505., 1506., 1507., 1508., 1509., 1510., 1511., 1512., 1513., 1514., 1515.,
+    1600., 1601., 1602., 1603., 1604., 1605., 1606., 1607., 1608., 1609., 1610., 1611., 1612., 1613., 1614., 1615.,
+    1700., 1701., 1702., 1703., 1704., 1705., 1706., 1707., 1708., 1709., 1710., 1711., 1712., 1713., 1714., 1715.,
+    1800., 1801., 1802., 1803., 1804., 1805., 1806., 1807., 1808., 1809., 1810., 1811., 1812., 1813., 1814., 1815.,
+    1900., 1901., 1902., 1903., 1904., 1905., 1906., 1907., 1908., 1909., 1910., 1911., 1912., 1913., 1914., 1915.,
+    2000., 2001., 2002., 2003., 2004., 2005., 2006., 2007., 2008., 2009., 2010., 2011., 2012., 2013., 2014., 2015.,
+    116.,216.,316.,416.,516.,616.,716.,816.,916.,1016.,1116.,1216.,1316.,1416.,1516.,1616.,1716.,1816.,1916.,2016.
+  };
+  CPPUNIT_ASSERT_EQUAL(1,remapper.prepare(meshNE,meshTFE,"P0P0"));
+  trgField=MEDCouplingFieldDouble::New(ON_CELLS);
+  trgField->setNature(ConservativeVolumic);
+  trgField->setMesh(meshTFE);
+  array=DataArrayDouble::New();
+  array->alloc(meshTFE->getNumberOfCells(),1);
+  trgField->setArray(array);
+  std::copy(vals3,vals3+meshTFE->getNumberOfCells(),array->getPointer());
+  array->decrRef();
+  srcField=remapper.reverseTransferField(trgField,4.220173);
+  double expected4[40]={
+    566.,566.,552.5,553.5,554.5,566.,566.,554.5,555.5,556.5,566.,566.,558.5,559.5,560.5,566.,566.,560.5,561.5,562.5,
+    1566.,1566.,1552.5,1553.5,1554.5,1566.,1566.,1554.5,1555.5,1556.5,1566.,1566.,1558.5,1559.5,1560.5,1566.,1566.,1560.5,1561.5,1562.5
+  };
+  for(int i=0;i<40;i++)
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(expected4[i],srcField->getArray()->getConstPointer()[i],1e-3);//1e-3 precision due to non coincidence in 1D mesh
+  srcField->decrRef();
+  trgField->decrRef();
+  //
+  meshN2D->decrRef();
+  meshTT2D->decrRef();
+  meshTF2D->decrRef();
+  meshNE->decrRef();
+  meshTTE->decrRef();
+  meshTFE->decrRef();
+  meshN->decrRef();
+  meshTT->decrRef();
+  meshTF->decrRef();
+}
+
 MEDCouplingUMesh *MEDCouplingRemapperTest::build1DTargetMesh_2()
 {
   double targetCoords[20]={
diff --git a/src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx b/src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx
index 5bcebfb89..29cf7cbcd 100644
--- a/src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx
+++ b/src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx
@@ -37,6 +37,7 @@ namespace ParaMEDMEM
     CPPUNIT_TEST( testMultiDimCombi );
     CPPUNIT_TEST( testNatureOfField );
     CPPUNIT_TEST( testExtruded );
+    CPPUNIT_TEST( testExtruded2 );
     CPPUNIT_TEST_SUITE_END();
   public:
     void test2DInterpP0P0_1();
@@ -45,6 +46,7 @@ namespace ParaMEDMEM
     void testMultiDimCombi();
     void testNatureOfField();
     void testExtruded();
+    void testExtruded2();
   private:
     static MEDCouplingUMesh *build1DTargetMesh_2();
     static MEDCouplingUMesh *build2DTargetMesh_3();