Back to index

salome-med  6.5.0
MEDMEMTest_Formulae.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
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 "MEDMEMTest.hxx"
00021 #include <cppunit/TestAssert.h>
00022 
00023 #include "VolSurfFormulae.hxx"
00024 #include "MEDMEM_STRING.hxx"
00025 #include "MEDMEM_Exception.hxx"
00026 
00027 #include <iostream>
00028 #include <sstream>
00029 #include <cmath>
00030 #include <cfloat>
00031 
00032 // use this define to enable lines, execution of which leads to Segmentation Fault
00033 //#define ENABLE_FAULTS
00034 
00035 // use this define to enable CPPUNIT asserts and fails, showing bugs
00036 //#define ENABLE_FORCED_FAILURES
00037 
00038 using namespace std;
00039 using namespace MEDMEM;
00040 
00041 // #17: MEDMEM_Formulae.hxx  }  MEDMEMTest_Formulae.cxx
00042 
00084 void MEDMEMTest::testFormulae()
00085 {
00086   double val;
00087 
00088   //      ^y
00089   //      |
00090   //      *3
00091   //  4   |
00092   //  .*. . . . . . *6,7
00093   //      |
00094   //  . .8* . . . . .
00095   //      |  5   2,9
00096   //  . . . .*. * . .
00097   //     1|
00098   // -.-.-*-.-.-.-.-.-->x
00099   //      |
00100 
00101   // S_12634 = (3 + (3+4.8)*2 + 4.8 + 1.5*4)/2 = 1.5 + 7.8 + 2.4 + 3 = 14.7
00102   // S_143652 = S_14362 - S_625 = 14.7 - 2 * 1.5 / 2 = 13.2
00103 
00104   double xy1[2]  = { 0.0, 0.0};
00105   double xy2[2]  = { 3.0, 1.0};
00106   double xy3[2]  = { 0.0, 4.0};
00107   double xy4[2]  = {-1.5, 3.0};
00108   double xy5[2]  = { 1.5, 1.0};
00109   double xy6[2]  = { 4.8, 3.0};
00110 
00111   double xyz1[3] = { 0.0, 0.0,  0.0};
00112   double xyz2[3] = { 3.0, 1.0,  4.0}; // cos(alpha) = 3/5
00113   double xyz3[3] = { 0.0, 4.0,  0.0};
00114   double xyz4[3] = {-1.5, 3.0, -2.0}; // z4 = z2 * x4 / x2 = - 4 * 1.5 / 3
00115   double xyz5[3] = { 1.5, 1.0,  2.0}; // z5 = z2 * x5 / x2 = 4 * 1.5 / 3
00116   double xyz6[3] = { 4.8, 3.0,  6.4}; // z6 = z2 * x6 / x2 = 4 * 4.8 / 3
00117   double xyz7[3] = { 4.8, 3.0,  0.0};
00118   double xyz8[3] = { 0.0, 2.0,  0.0};
00119   double xyz9[3] = { 3.0, 1.0,  0.0};
00120 
00121   // S_3d = S_2d * 5.0 / 3.0
00122 
00124   // CalculateAreaForPolyg //
00126   {
00127     // 2D: Convex polygon
00128     const double * poly_2d_cc[5] = {xy1, xy2, xy6, xy3, xy4};
00129     const double * poly_2d_cw[5] = {xy1, xy4, xy3, xy6, xy2};
00130 
00131     // counter-clockwise
00132     val = INTERP_KERNEL::calculateAreaForPolyg(poly_2d_cc, /*nbOfPtsInPolygs*/5, /*dim*/2);
00133     CPPUNIT_ASSERT_DOUBLES_EQUAL(-14.7, val, 0.000001);
00134 
00135     // clockwise
00136     val = INTERP_KERNEL::calculateAreaForPolyg(poly_2d_cw, /*nbOfPtsInPolygs*/5, /*dim*/2);
00137     CPPUNIT_ASSERT_DOUBLES_EQUAL(14.7, val, 0.000001);
00138 
00139     // 2D: Non-convex polygon
00140     const double * poly_2d_nc[6] = {xy1, xy4, xy3, xy6, xy5, xy2};
00141     val = INTERP_KERNEL::calculateAreaForPolyg(poly_2d_nc, /*nbOfPtsInPolygs*/6, /*dim*/2);
00142     CPPUNIT_ASSERT_DOUBLES_EQUAL(13.2, val, 0.000001);
00143 
00144     // 3D: Convex polygon
00145     const double * poly_3d_cc[5] = {xyz1, xyz2, xyz6, xyz3, xyz4};
00146 
00147     val = INTERP_KERNEL::calculateAreaForPolyg(poly_3d_cc, /*nbOfPtsInPolygs*/5, /*dim*/3);
00148     CPPUNIT_ASSERT_DOUBLES_EQUAL(24.5, val, 0.000001);
00149 
00150     // 3D: Non-convex polygon
00151     const double * poly_3d_nc[6] = {xyz1, xyz4, xyz3, xyz6, xyz5, xyz2};
00152     val = INTERP_KERNEL::calculateAreaForPolyg(poly_3d_nc, /*nbOfPtsInPolygs*/6, /*dim*/3);
00153     //#ifdef ENABLE_FORCED_FAILURES
00154     // (BUG) Wrong area calculation for non-convex polygons in 3D space,
00155     //       because area of triangle is always positive
00156     // It's a spec
00157     //CPPUNIT_ASSERT_DOUBLES_EQUAL(22.0, val, 0.000001);
00158     CPPUNIT_ASSERT( abs(22.0-val)>0.000001 );
00159     //#endif
00160   }
00161 
00163   // CalculateAreaForTria //
00165   {
00166     // 2D: counter-clockwise
00167     val = INTERP_KERNEL::calculateAreaForTria(xy1, xy2, xy3, /*dim*/2);
00168     CPPUNIT_ASSERT_DOUBLES_EQUAL(-6.0, val, 0.000001);
00169 
00170     // 2D: clockwise
00171     val = INTERP_KERNEL::calculateAreaForTria(xy2, xy1, xy3, /*dim*/2);
00172     CPPUNIT_ASSERT_DOUBLES_EQUAL(6.0, val, 0.000001);
00173 
00174     // 3D
00175     val = INTERP_KERNEL::calculateAreaForTria(xyz1, xyz2, xyz3, /*dim*/3);
00176     CPPUNIT_ASSERT_DOUBLES_EQUAL(10.0, val, 0.000001);
00177 
00178     // Invalid (three points on one line)
00179     val = INTERP_KERNEL::calculateAreaForTria(xyz1, xyz8, xyz3, /*dim*/3);
00180     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, val, 0.000001);
00181   }
00182 
00184   // CalculateAreaForQuad //
00186   {
00187     // 2D: Convex quadrangle
00188 
00189     // counter-clockwise
00190     val = INTERP_KERNEL::calculateAreaForQuad(xy1, xy2, xy3, xy4, /*dim*/2);
00191     CPPUNIT_ASSERT_DOUBLES_EQUAL(-9.0, val, 0.000001);
00192 
00193     // clockwise
00194     val = INTERP_KERNEL::calculateAreaForQuad(xy2, xy1, xy4, xy3, /*dim*/2);
00195     CPPUNIT_ASSERT_DOUBLES_EQUAL(9.0, val, 0.000001);
00196 
00197     // wrong order
00198     CPPUNIT_ASSERT_NO_THROW(INTERP_KERNEL::calculateAreaForQuad(xy2, xy1, xy3, xy4, /*dim*/2));
00199 
00200     // 2D: Non-convex quadrangle
00201 
00202     // counter-clockwise
00203     val = INTERP_KERNEL::calculateAreaForQuad(xy1, xy2, xy3, xy5, /*dim*/2);
00204     CPPUNIT_ASSERT_DOUBLES_EQUAL(-3.0, val, 0.000001);
00205 
00206     // clockwise
00207     val = INTERP_KERNEL::calculateAreaForQuad(xy1, xy5, xy3, xy2, /*dim*/2);
00208     CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0, val, 0.000001);
00209 
00210     // 3D: Convex quadrangle
00211 
00212     // good order
00213     val = INTERP_KERNEL::calculateAreaForQuad(xyz1, xyz2, xyz3, xyz4, /*dim*/3);
00214     CPPUNIT_ASSERT_DOUBLES_EQUAL(15.0, val, 0.000001);
00215 
00216     // wrong order
00217     CPPUNIT_ASSERT_NO_THROW(INTERP_KERNEL::calculateAreaForQuad(xyz1, xyz4, xyz2, xyz3, /*dim*/3));
00218 
00219     // 3D: Non-convex quadrangle
00220     val = INTERP_KERNEL::calculateAreaForQuad(xyz1, xyz2, xyz3, xyz5, /*dim*/3);
00221     CPPUNIT_ASSERT_DOUBLES_EQUAL(5.0, val, 0.000001);
00222 
00223     // 3D: Non-planar quadrangle
00224     double xyz7[3] = {-1.5, 3.0,  2.0};
00225     CPPUNIT_ASSERT_NO_THROW(INTERP_KERNEL::calculateAreaForQuad(xyz1, xyz2, xyz3, xyz7, /*dim*/3));
00226   }
00227 
00229   // CalculateNormalForTria //
00231   {
00232     double tria_normal [3];
00233 
00234     // Triangle in plane XOY, normal is opposit to axis Z
00235     INTERP_KERNEL::calculateNormalForTria(xyz1, xyz3, xyz7, tria_normal);
00236     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, tria_normal[0], 0.0000001); // Nx
00237     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, tria_normal[1], 0.0000001); // Ny
00238     CPPUNIT_ASSERT(tria_normal[2] < -0.0000001); // Nz
00239 
00240     // Triangle in plane XOY, normal co-directed with axis Z
00241     INTERP_KERNEL::calculateNormalForTria(xyz1, xyz7, xyz3, tria_normal);
00242     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, tria_normal[0], 0.0000001); // Nx
00243     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, tria_normal[1], 0.0000001); // Ny
00244     CPPUNIT_ASSERT(tria_normal[2] > 0.0000001); // Nz
00245 
00246     // Triangle in 3D
00247     INTERP_KERNEL::calculateNormalForTria(xyz1, xyz3, xyz6, tria_normal);
00248     double koeff = tria_normal[0]/12.8;
00249     //CPPUNIT_ASSERT_DOUBLES_EQUAL(12.8, tria_normal[0], 0.0000001); // Nx
00250     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0        , tria_normal[1], 0.0000001); // Ny
00251     CPPUNIT_ASSERT_DOUBLES_EQUAL(-9.6 * koeff, tria_normal[2], 0.0000001); // Nz
00252 
00253     // Invalid Triangle (three points on one line)
00254     CPPUNIT_ASSERT_NO_THROW(INTERP_KERNEL::calculateNormalForTria(xyz1, xyz8, xyz3, tria_normal));
00255     //MEDMEMTest_DumpArray<double>(cout, tria_normal, 3, "Invalid Triangle normal");
00256     //Invalid Triangle normal: {0, 0, 0}
00257   }
00258 
00260   // CalculateNormalForQuad //
00262   {
00263     double quad_normal [3];
00264 
00265     // Quadrangle in plane XOY, normal is opposit to axis Z
00266     INTERP_KERNEL::calculateNormalForQuad(xyz1, xyz3, xyz7, xyz9, quad_normal);
00267     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, quad_normal[0], 0.0000001); // Nx
00268     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, quad_normal[1], 0.0000001); // Ny
00269     CPPUNIT_ASSERT(quad_normal[2] < -0.0000001); // Nz
00270 
00271     // Quadrangle in plane XOY, normal co-directed with axis Z
00272     INTERP_KERNEL::calculateNormalForQuad(xyz1, xyz9, xyz7, xyz3, quad_normal);
00273     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, quad_normal[0], 0.0000001); // Nx
00274     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, quad_normal[1], 0.0000001); // Ny
00275     CPPUNIT_ASSERT(quad_normal[2] > 0.0000001); // Nz
00276 
00277     // Quadrangle in 3D
00278     INTERP_KERNEL::calculateNormalForQuad(xyz1, xyz3, xyz6, xyz2, quad_normal);
00279     //MEDMEMTest_DumpArray<double>(cout, quad_normal, 3, "Quadrangle in 3D normal");
00280     double koeff = quad_normal[0]/15.6;
00281     //CPPUNIT_ASSERT_DOUBLES_EQUAL(15.6, quad_normal[0], 0.0000001); // Nx
00282     CPPUNIT_ASSERT_DOUBLES_EQUAL(  0.0        , quad_normal[1], 0.0000001); // Ny
00283     CPPUNIT_ASSERT_DOUBLES_EQUAL(-11.7 * koeff, quad_normal[2], 0.0000001); // Nz
00284 
00285     // Invalid Quadrangle (four points on one line)
00286     CPPUNIT_ASSERT_NO_THROW(INTERP_KERNEL::calculateNormalForQuad(xyz1, xyz8, xyz3, xyz3, quad_normal));
00287     //MEDMEMTest_DumpArray<double>(cout, quad_normal, 3, "Invalid Quadrangle normal");
00288     //#ifdef ENABLE_FORCED_FAILURES
00289     // (BUG) division on zero in CalculateNormalForQuad(), if quadrangle is singular
00290     //Invalid Quadrangle normal: {nan, nan, nan}
00291     // Spec. Wait for an improvement
00292     // CPPUNIT_ASSERT(quad_normal[0] < DBL_MAX);
00293     // CPPUNIT_ASSERT(quad_normal[1] < DBL_MAX);
00294     // CPPUNIT_ASSERT(quad_normal[2] < DBL_MAX);
00295     //#endif
00296   }
00297 
00299   // CalculateNormalForPolyg //
00301   {
00302     double poly_normal [3];
00303     const double * polygon_cc[4] = {xyz1, xyz3, xyz7, xyz9};
00304     const double * polygon_er[4] = {xyz1, xyz8, xyz3, xyz7};
00305     const double * polygon_cw[4] = {xyz1, xyz9, xyz7, xyz3};
00306     const double * polygon_3d[4] = {xyz1, xyz3, xyz6, xyz2};
00307     const double * polygon_si[4] = {xyz1, xyz8, xyz3, xyz3};
00308 
00309     // Polygon in plane XOY, normal is opposit to axis Z
00310     INTERP_KERNEL::calculateNormalForPolyg(polygon_cc, /*nbOfPtsInPolygs*/4, poly_normal);
00311     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, poly_normal[0], 0.0000001); // Nx
00312     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, poly_normal[1], 0.0000001); // Ny
00313     //#ifdef ENABLE_FORCED_FAILURES
00314     // (BUG) Normal for polygon is wrong. Confired to be a bug
00315     CPPUNIT_ASSERT_DOUBLES_EQUAL(-11.7, poly_normal[2], 0.0000001); // Nz
00316     //#endif
00317 
00318     // Polygon in plane XOY, normal co-directed with axis Z
00319     INTERP_KERNEL::calculateNormalForPolyg(polygon_cw, /*nbOfPtsInPolygs*/4, poly_normal);
00320     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, poly_normal[0], 0.0000001); // Nx
00321     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, poly_normal[1], 0.0000001); // Ny
00322     //#ifdef ENABLE_FORCED_FAILURES
00323     // (BUG) Normal for polygon is wrong. Confired to be a bug
00324     CPPUNIT_ASSERT_DOUBLES_EQUAL(11.7, poly_normal[2], 0.0000001); // Nz
00325     //#endif
00326 
00327     // Polygon in plane XOY, normal is opposit to axis Z, first three points lay on one line
00328     //CPPUNIT_ASSERT_THROW(CalculateNormalForPolyg(polygon_er, /*nbOfPtsInPolygs*/4,
00329     //                                             poly_normal),MEDEXCEPTION);
00330     CPPUNIT_ASSERT_NO_THROW(INTERP_KERNEL::calculateNormalForPolyg(polygon_er, /*nbOfPtsInPolygs*/4,
00331                                                     poly_normal));
00332     //CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, poly_normal[0], 0.0000001); // Nx
00333     //CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, poly_normal[1], 0.0000001); // Ny
00334     //#ifdef ENABLE_FORCED_FAILURES
00335     // (BUG) Normal for polygon is wrong if first (three) points are on one line
00336     // Spec, can be improved
00337     // CPPUNIT_ASSERT(poly_normal[2] < -0.0000001); // Nz
00338     //#endif
00339 
00340     // Polygon in 3D
00341     INTERP_KERNEL::calculateNormalForPolyg(polygon_3d, /*nbOfPtsInPolygs*/4, poly_normal);
00342     double koeff = poly_normal[0]/15.6;
00343     //CPPUNIT_ASSERT_DOUBLES_EQUAL(15.6, poly_normal[0], 0.0000001); // Nx
00344     CPPUNIT_ASSERT_DOUBLES_EQUAL(  0.0        , poly_normal[1], 0.0000001); // Ny
00345     CPPUNIT_ASSERT_DOUBLES_EQUAL(-11.7 * koeff, poly_normal[2], 0.0000001); // Nz
00346 
00347     // Invalid Polygon (four points on one line)
00348     bool isException=false;
00349     try
00350       {
00351         INTERP_KERNEL::calculateNormalForPolyg(polygon_si, /*nbOfPtsInPolygs*/4,
00352                                                poly_normal);
00353       }
00354     catch(std::exception& e)
00355       {
00356         isException=true;
00357       }
00358     CPPUNIT_ASSERT(isException);
00359     //MEDMEMTest_DumpArray<double>(cout, poly_normal, 3, "Invalid Polygon normal");
00360     //#ifdef ENABLE_FORCED_FAILURES
00361     // (BUG) division on zero in CalculateNormalForPolyg(), if polygon is singular
00362     //Invalid Polygon normal: {nan, nan, nan}
00363     // Cinfirmed to be a bug
00364     //CPPUNIT_ASSERT(poly_normal[0] < DBL_MAX);
00365     //CPPUNIT_ASSERT(poly_normal[1] < DBL_MAX);
00366     //CPPUNIT_ASSERT(poly_normal[2] < DBL_MAX);
00367     //#endif
00368   }
00369 
00371   // CalculateVolumeForTetra //
00373   {
00374     // good
00375     val = INTERP_KERNEL::calculateVolumeForTetra(xyz1, xyz3, xyz7, xyz5);
00376     CPPUNIT_ASSERT_DOUBLES_EQUAL(6.4, val, 0.000001);
00377 
00378     // reversed
00379     val = INTERP_KERNEL::calculateVolumeForTetra(xyz1, xyz7, xyz3, xyz5);
00380     CPPUNIT_ASSERT_DOUBLES_EQUAL(-6.4, val, 0.000001);
00381 
00382     // good
00383     val = INTERP_KERNEL::calculateVolumeForTetra(xyz1, xyz7, xyz3, xyz4);
00384     CPPUNIT_ASSERT_DOUBLES_EQUAL(6.4, val, 0.000001);
00385 
00386     // reversed
00387     val = INTERP_KERNEL::calculateVolumeForTetra(xyz1, xyz3, xyz7, xyz4);
00388     CPPUNIT_ASSERT_DOUBLES_EQUAL(-6.4, val, 0.000001);
00389 
00390     // singular (in plane)
00391     val = INTERP_KERNEL::calculateVolumeForTetra(xyz1, xyz3, xyz7, xyz9);
00392     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, val, 0.000001);
00393   }
00394 
00396   // CalculateVolumeForPyra //
00398   {
00399     // good
00400     val = INTERP_KERNEL::calculateVolumeForPyra(xyz1, xyz3, xyz7, xyz9, xyz5);
00401     CPPUNIT_ASSERT_DOUBLES_EQUAL(7.8, val, 0.000001);
00402 
00403     // reversed
00404     val = INTERP_KERNEL::calculateVolumeForPyra(xyz1, xyz9, xyz7, xyz3, xyz5);
00405     CPPUNIT_ASSERT_DOUBLES_EQUAL(-7.8, val, 0.000001);
00406 
00407     // good
00408     val = INTERP_KERNEL::calculateVolumeForPyra(xyz1, xyz9, xyz7, xyz3, xyz4);
00409     CPPUNIT_ASSERT_DOUBLES_EQUAL(7.8, val, 0.000001);
00410 
00411     // reversed
00412     val = INTERP_KERNEL::calculateVolumeForPyra(xyz1, xyz3, xyz7, xyz9, xyz4);
00413     CPPUNIT_ASSERT_DOUBLES_EQUAL(-7.8, val, 0.000001);
00414 
00415     // singular (in plane)
00416     val = INTERP_KERNEL::calculateVolumeForPyra(xyz1, xyz3, xyz7, xyz9, xyz8);
00417     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, val, 0.000001);
00418   }
00419 
00421   // CalculateVolumeForPenta //
00423   {
00424     double top1[3] = {xyz1[0], xyz1[1], xyz1[2] + 10.0};
00425     double top3[3] = {xyz3[0], xyz3[1], xyz3[2] + 10.0};
00426     double top7[3] = {xyz7[0], xyz7[1], xyz7[2] + 10.0};
00427 
00428     // good
00429     val = INTERP_KERNEL::calculateVolumeForPenta(xyz1, xyz3, xyz7, top1, top3, top7);
00430     CPPUNIT_ASSERT_DOUBLES_EQUAL(96.0, val, 0.000001);
00431 
00432     // reversed
00433     val = INTERP_KERNEL::calculateVolumeForPenta(xyz1, xyz7, xyz3, top1, top7, top3);
00434     CPPUNIT_ASSERT_DOUBLES_EQUAL(-96.0, val, 0.000001);
00435 
00436     // good
00437     top1[0] = top1[0] + 7.0;
00438     top3[0] = top3[0] + 7.0;
00439     top7[0] = top7[0] + 7.0;
00440 
00441     val = INTERP_KERNEL::calculateVolumeForPenta(xyz1, xyz3, xyz7, top1, top3, top7);
00442     CPPUNIT_ASSERT_DOUBLES_EQUAL(96.0, val, 0.000001);
00443 
00444     // reversed
00445     val = INTERP_KERNEL::calculateVolumeForPenta(xyz1, xyz7, xyz3, top1, top7, top3);
00446     CPPUNIT_ASSERT_DOUBLES_EQUAL(-96.0, val, 0.000001);
00447 
00448     // singular (in plane)
00449     top1[2] = top1[2] - 10.0;
00450     top3[2] = top3[2] - 10.0;
00451     top7[2] = top7[2] - 10.0;
00452 
00453     val = INTERP_KERNEL::calculateVolumeForPenta(xyz1, xyz3, xyz7, top1, top3, top7);
00454     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, val, 0.000001);
00455   }
00456 
00458   // CalculateVolumeForHexa //
00460   {
00461     double top1[3] = {xyz1[0], xyz1[1], xyz1[2] + 10.0};
00462     double top3[3] = {xyz3[0], xyz3[1], xyz3[2] + 10.0};
00463     double top7[3] = {xyz7[0], xyz7[1], xyz7[2] + 10.0};
00464     double top9[3] = {xyz9[0], xyz9[1], xyz9[2] + 10.0};
00465 
00466     // good
00467     val = INTERP_KERNEL::calculateVolumeForHexa(xyz1, xyz3, xyz7, xyz9, top1, top3, top7, top9);
00468     CPPUNIT_ASSERT_DOUBLES_EQUAL(117.0, val, 0.000001);
00469 
00470     // reversed
00471     val = INTERP_KERNEL::calculateVolumeForHexa(xyz1, xyz9, xyz7, xyz3, top1, top9, top7, top3);
00472     CPPUNIT_ASSERT_DOUBLES_EQUAL(-117.0, val, 0.000001);
00473 
00474     // good
00475     top1[0] = top1[0] + 7.0;
00476     top3[0] = top3[0] + 7.0;
00477     top7[0] = top7[0] + 7.0;
00478     top9[0] = top9[0] + 7.0;
00479 
00480     val = INTERP_KERNEL::calculateVolumeForHexa(xyz1, xyz3, xyz7, xyz9, top1, top3, top7, top9);
00481     CPPUNIT_ASSERT_DOUBLES_EQUAL(117.0, val, 0.000001);
00482 
00483     // reversed
00484     val = INTERP_KERNEL::calculateVolumeForHexa(xyz1, xyz9, xyz7, xyz3, top1, top9, top7, top3);
00485     CPPUNIT_ASSERT_DOUBLES_EQUAL(-117.0, val, 0.000001);
00486 
00487     // singular (in plane)
00488     top1[2] = top1[2] - 10.0;
00489     top3[2] = top3[2] - 10.0;
00490     top7[2] = top7[2] - 10.0;
00491     top9[2] = top9[2] - 10.0;
00492 
00493     val = INTERP_KERNEL::calculateVolumeForHexa(xyz1, xyz3, xyz7, xyz9, top1, top3, top7, top9);
00494     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, val, 0.000001);
00495   }
00496 
00498   // CalculateVolumeForPolyh //
00500   // inline double CalculateVolumeForPolyh(const double ***pts, const int *nbOfNodesPerFaces,
00501   //                                       int nbOfFaces, const double *bary);
00502   {
00503     int nbFaces = 4;
00504     int nbOfNodesPerFaces[4] = {3,3,3,3};
00505 
00506     const double * nodes[4] = {xyz1, xyz7, xyz3, xyz5};
00507     double bary[3];
00508     INTERP_KERNEL::calculateBarycenterDyn(nodes, 4, /*dim*/3, bary);
00509 
00510     // good
00511     const double * fa1[3] = {xyz1, xyz7, xyz3};
00512     const double * fa2[3] = {xyz1, xyz3, xyz5};
00513     const double * fa3[3] = {xyz3, xyz7, xyz5};
00514     const double * fa4[3] = {xyz7, xyz1, xyz5};
00515     const double ** polyh[4] = {fa1, fa2, fa3, fa4};
00516 
00517     val = INTERP_KERNEL::calculateVolumeForPolyh(polyh, nbOfNodesPerFaces, nbFaces, bary);
00518     CPPUNIT_ASSERT_DOUBLES_EQUAL(-6.4, val, 0.000001);
00519 
00520     // reversed
00521     //const double * fa1_r[3] = {xyz1, xyz3, xyz7};
00522     //const double * fa2_r[3] = {xyz3, xyz1, xyz5};
00523     //const double * fa3_r[3] = {xyz7, xyz3, xyz5};
00524     //const double * fa4_r[3] = {xyz1, xyz7, xyz5};
00525     //const double ** polyh_r[4] = {fa1_r, fa2_r, fa3_r, fa4_r};
00526     //
00527     //val = CalculateVolumeForPolyh(polyh_r, nbOfNodesPerFaces, nbFaces, bary);
00528     //CPPUNIT_ASSERT_DOUBLES_EQUAL(-6.4, val, 0.000001);
00529 
00530     // singular (in plane)
00531     const double * nodes_si[4] = {xyz1, xyz7, xyz3, xyz9};
00532     double bary_si[3];
00533     INTERP_KERNEL::calculateBarycenterDyn(nodes_si, 4, /*dim*/3, bary_si);
00534 
00535     const double * fa1_si[3] = {xyz1, xyz7, xyz3};
00536     const double * fa2_si[3] = {xyz1, xyz3, xyz9};
00537     const double * fa3_si[3] = {xyz3, xyz7, xyz9};
00538     const double * fa4_si[3] = {xyz7, xyz1, xyz9};
00539     const double ** polyh_si[4] = {fa1_si, fa2_si, fa3_si, fa4_si};
00540 
00541     val = INTERP_KERNEL::calculateVolumeForPolyh(polyh_si, nbOfNodesPerFaces, nbFaces, bary_si);
00542     CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, val, 0.000001);
00543   }
00544 
00546   // addComponentsOfVec //
00548   {
00549     // five points
00550     const double * pts[5] = {xyz2, xyz1, xyz3, xyz4, xyz5};
00551 
00552     val = INTERP_KERNEL::addComponentsOfVec<5>(pts, 0); // x
00553     CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0, val, 0.000001);
00554 
00555     val = INTERP_KERNEL::addComponentsOfVec<5>(pts, 1); // y
00556     CPPUNIT_ASSERT_DOUBLES_EQUAL(9.0, val, 0.000001);
00557 
00558     val = INTERP_KERNEL::addComponentsOfVec<5>(pts, 2); // z
00559     CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0, val, 0.000001);
00560 
00561     // one point: xyz2
00562     val = INTERP_KERNEL::addComponentsOfVec<1>(pts, 0); // x
00563     CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0, val, 0.000001);
00564 
00565     val = INTERP_KERNEL::addComponentsOfVec<1>(pts, 1); // y
00566     CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, val, 0.000001);
00567 
00568     val = INTERP_KERNEL::addComponentsOfVec<1>(pts, 2); // z
00569     CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0, val, 0.000001);
00570   }
00571 
00573   // CalculateBarycenterDyn //
00575   {
00576     // five points
00577     const double * pts[5] = {xyz2, xyz1, xyz3, xyz4, xyz5};
00578     double bary_3d[3];
00579     double bary_2d[2];
00580 
00581     INTERP_KERNEL::calculateBarycenterDyn(pts, /*nbPts*/5, /*dim*/3, bary_3d);
00582     CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0/5.0, bary_3d[0], 0.000001); // x
00583     CPPUNIT_ASSERT_DOUBLES_EQUAL(9.0/5.0, bary_3d[1], 0.000001); // y
00584     CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0/5.0, bary_3d[2], 0.000001); // z
00585 
00586     INTERP_KERNEL::calculateBarycenterDyn(pts, /*nbPts*/5, /*dim*/2, bary_2d);
00587     CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0/5.0, bary_2d[0], 0.000001); // x
00588     CPPUNIT_ASSERT_DOUBLES_EQUAL(9.0/5.0, bary_2d[1], 0.000001); // y
00589 
00590     // one point: xyz2
00591     INTERP_KERNEL::calculateBarycenterDyn(pts, /*nbPts*/1, /*dim*/3, bary_3d);
00592     CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0, bary_3d[0], 0.000001); // x
00593     CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, bary_3d[1], 0.000001); // y
00594     CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0, bary_3d[2], 0.000001); // z
00595   }
00596 
00598   // CalculateBarycenter //
00600   // template<int N, int DIM> inline void CalculateBarycenter(const double **pts, double *bary);
00601   // template<> inline void CalculateBarycenter<2,0>(const double **pts, double *bary);
00602   // template<> inline void CalculateBarycenter<3,0>(const double **pts, double *bary);
00603   // template<> inline void CalculateBarycenter<4,0>(const double **pts, double *bary);
00604   // template<> inline void CalculateBarycenter<5,0>(const double **pts, double *bary);
00605   // template<> inline void CalculateBarycenter<6,0>(const double **pts, double *bary);
00606   // template<> inline void CalculateBarycenter<7,0>(const double **pts, double *bary);
00607   // template<> inline void CalculateBarycenter<8,0>(const double **pts, double *bary);
00608   {
00609     // five points
00610     const double * pts[5] = {xyz2, xyz1, xyz3, xyz4, xyz5};
00611     double bary_3d[3];
00612     double bary_2d[2];
00613 
00614     INTERP_KERNEL::calculateBarycenter<5,3>(pts, bary_3d);
00615     CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0/5.0, bary_3d[0], 0.000001); // x
00616     CPPUNIT_ASSERT_DOUBLES_EQUAL(9.0/5.0, bary_3d[1], 0.000001); // y
00617     CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0/5.0, bary_3d[2], 0.000001); // z
00618 
00619     INTERP_KERNEL::calculateBarycenter<5,2>(pts, bary_2d);
00620     CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0/5.0, bary_2d[0], 0.000001); // x
00621     CPPUNIT_ASSERT_DOUBLES_EQUAL(9.0/5.0, bary_2d[1], 0.000001); // y
00622 
00623     // one point: xyz2 : NOT IMPLEMENTED!!!
00624     //CalculateBarycenter<1,3>(pts, bary_3d);
00625     //CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0, bary_3d[0], 0.000001);
00626     //CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, bary_3d[1], 0.000001);
00627     //CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0, bary_3d[2], 0.000001);
00628   }
00629 }