Back to index

salome-med  6.5.0
TransformedTriangleTest.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 "TransformedTriangleTest.hxx"
00021 
00022 #include <iostream>
00023 
00024 using namespace INTERP_KERNEL;
00025 
00026 namespace INTERP_TEST
00027 {
00028 
00033   void TransformedTriangleTest::setUp() 
00034   {
00035     // tri1 -> no unstable double products - no changes brought about by preCalculateDoubleProducts
00036     //         this allows the testing of calcUnstableT
00037     // tri2 -> unstable double products - for testing calcStableC / preCalculateDoubleProducts
00038 
00039     // triangle to test unstable C and T calculations
00040     p1[0] = -1.5 ; p1[1] = 0.5; p1[2] = 0.5;
00041     q1[0] = 2.0 ; q1[1] = 0.4; q1[2] = 0.6;
00042     r1[0] = 1.0 ; r1[1] = 2.4; r1[2] = 1.2;
00043     hp1 = 1 - p1[0] - p1[1] - p1[2];
00044     hq1 = 1 - q1[0] - q1[1] - q1[2];
00045     hr1 = 1 - r1[0] - r1[1] - r1[2]; 
00046     Hp1 = 1 - p1[0] - p1[1];
00047     Hq1 = 1 - q1[0] - q1[1];
00048     Hr1 = 1 - r1[0] - r1[1];
00049 
00050     //  std::cout <<std::endl<< "constructing tri1..." << std::endl;
00051     tri1 = new TransformedTriangle(p1, q1, r1);
00052  
00053 
00054     // triangle to test stable C calculation
00055     const double err = 1.5e-3;
00056     
00057     p2[0] = 0.000000000084654984189118; p2[1] = -0.000000000000000027536546231654231688873; p2[2] = 0.0000000000000001649875466831349431;
00058     q2[0] = -p2[0] +err; q2[1] = -p2[1] + err; q2[2] = -p2[2] +err;
00059     r2[0] = 2.01 ; r2[1] = 1.8; r2[2] = 0.92;
00060     
00061     hp2 = 1 - p2[0] - p2[1] - p2[2];
00062     hq2 = 1 - q2[0] - q2[1] - q2[2];
00063     hr2 = 1 - r2[0] - r2[1] - r2[2]; 
00064     Hp2 = 1 - p2[0] - p2[1];
00065     Hq2 = 1 - q2[0] - q2[1];
00066     Hr2 = 1 - r2[0] - r2[1];
00067     
00068     tri2 = new TransformedTriangle(p2, q2, r2);
00069   
00070   
00071 
00072   }
00073 
00078   void TransformedTriangleTest::tearDown() 
00079   {
00080     delete tri1;
00081     delete tri2;
00082   }
00083 
00086   void TransformedTriangleTest::test_constructor() {
00087     // test that _coords has correct values after constructor is called
00088 
00089     double good_values1[15] = 
00090       {
00091         p1[0], p1[1], p1[2], hp1, Hp1,
00092         q1[0], q1[1], q1[2], hq1, Hq1,
00093         r1[0], r1[1], r1[2], hr1, Hr1
00094       };
00095 
00096     double good_values2[15] = 
00097       {
00098         p2[0], p2[1], p2[2], hp2, Hp2,
00099         q2[0], q2[1], q2[2], hq2, Hq2,
00100         r2[0], r2[1], r2[2], hr2, Hr2
00101       };
00102 
00103   
00104     for(int i = 0 ; i < 15 ; ++i)
00105       {
00106         CPPUNIT_ASSERT_DOUBLES_EQUAL(good_values1[i], tri1->_coords[i], ERR_TOL);
00107         CPPUNIT_ASSERT_DOUBLES_EQUAL(good_values2[i], tri2->_coords[i], ERR_TOL);
00108       }
00109 
00110     CPPUNIT_ASSERT_EQUAL(true, tri1->_is_double_products_calculated);
00111     CPPUNIT_ASSERT_EQUAL(true, tri2->_is_double_products_calculated);
00112   }
00113 
00116   void TransformedTriangleTest::test_calcUnstableC() {
00117     typedef TransformedTriangle::TriSegment TriSegment;
00118 
00119     // test that the correct c-values are calculated
00120   
00121     double correct_c_vals[24] = 
00122       { 
00123         p1[0] * q1[1] - p1[1] * q1[0], 
00124         p1[1] * q1[2] - p1[2] * q1[1], 
00125         p1[2] * q1[0] - p1[0] * q1[2],
00126         p1[0] * hq1 - hp1 * q1[0],
00127         p1[1] * hq1 - hp1 * q1[1],
00128         p1[2] * hq1 - hp1 * q1[2],
00129         Hp1 * q1[0] - p1[0] * Hq1,
00130         p1[1] * Hq1 - Hp1 * q1[1],
00131         q1[0] * r1[1] - q1[1] * r1[0], 
00132         q1[1] * r1[2] - q1[2] * r1[1], 
00133         q1[2] * r1[0] - q1[0] * r1[2],
00134         q1[0] * hr1 - hq1 * r1[0],
00135         q1[1] * hr1 - hq1 * r1[1],
00136         q1[2] * hr1 - hq1 * r1[2],
00137         Hq1 * r1[0] - q1[0] * Hr1,
00138         q1[1] * Hr1 - Hq1 * r1[1],
00139         r1[0]*p1[1]-r1[1]*p1[0], 
00140         r1[1]*p1[2]-r1[2]*p1[1], 
00141         r1[2]*p1[0]-r1[0]*p1[2],
00142         r1[0] * hp1 - hr1 * p1[0],
00143         r1[1] * hp1 - hr1 * p1[1],
00144         r1[2] * hp1 - hr1 * p1[2],
00145         Hr1 * p1[0] - r1[0] * Hp1,
00146         r1[1] * Hp1 - Hr1 * p1[1]
00147       };
00148 
00149     double c_vals[3 * 8];
00150     for(TriSegment seg = TransformedTriangle::PQ ; seg <= TransformedTriangle::RP ; seg = TriSegment(seg + 1)) {
00151       
00152       c_vals[seg*8 + 0] = tri1->calcUnstableC(seg, TransformedTriangle::C_XY);
00153       c_vals[seg*8 + 1] = tri1->calcUnstableC(seg, TransformedTriangle::C_YZ);
00154       c_vals[seg*8 + 2] = tri1->calcUnstableC(seg, TransformedTriangle::C_ZX);
00155       c_vals[seg*8 + 3] = tri1->calcUnstableC(seg, TransformedTriangle::C_XH);
00156       c_vals[seg*8 + 4] = tri1->calcUnstableC(seg, TransformedTriangle::C_YH);
00157       c_vals[seg*8 + 5] = tri1->calcUnstableC(seg, TransformedTriangle::C_ZH);
00158       c_vals[seg*8 + 6] = tri1->calcUnstableC(seg, TransformedTriangle::C_01);
00159       c_vals[seg*8 + 7] = tri1->calcUnstableC(seg, TransformedTriangle::C_10);
00160 
00161     }
00162     
00163     for(int i = 0 ; i < 3*8 ; ++i) {
00164       CPPUNIT_ASSERT_DOUBLES_EQUAL( correct_c_vals[i], c_vals[i], ERR_TOL );
00165     }
00166 
00167 
00168   }
00169 
00172   void TransformedTriangleTest::test_calcUnstableT()
00173   {
00174     typedef TransformedTriangle::TetraCorner TetraCorner;
00175 
00176     // correct values calculated by determinants (Grandy, [15])
00177     const double correct_t_vals[4] = 
00178       {
00179         p1[0]*(q1[1]*r1[2] - q1[2]*r1[1]) -
00180         q1[0]*(p1[1]*r1[2] - p1[2]*r1[1]) +
00181         r1[0]*(p1[1]*q1[2] - p1[2]*q1[1]),
00182 
00183         -(hp1*(q1[1]*r1[2] - q1[2]*r1[1]) -
00184           hq1*(p1[1]*r1[2] - p1[2]*r1[1]) +
00185           hr1*(p1[1]*q1[2] - p1[2]*q1[1])),
00186 
00187         -(p1[0]*(hq1*r1[2] - q1[2]*hr1) -
00188           q1[0]*(hp1*r1[2] - p1[2]*hr1) +
00189           r1[0]*(hp1*q1[2] - p1[2]*hq1)),
00190     
00191         -(p1[0]*(q1[1]*hr1 - r1[1]*hq1) -
00192           q1[0]*(p1[1]*hr1 - r1[1]*hp1) +
00193           r1[0]*(p1[1]*hq1 - q1[1]*hp1))
00194       };
00195     
00196 
00197     // test that triple products are correctly calculated
00198     for(TetraCorner corner = TransformedTriangle::O ; corner <= TransformedTriangle::Z ; corner = TetraCorner(corner + 1)) 
00199       {
00200       
00201         for(int row = 1 ; row < 4 ; ++row)
00202           {
00203             const double t = tri1->calcTByDevelopingRow(corner, row, false);
00204             //    std::cout << std::endl  << " Corner = " << corner  << " Row = " << row << " got: " << t << 
00205             //  " expected: " << correct_t_vals[corner]<< std::endl;
00206             CPPUNIT_ASSERT_DOUBLES_EQUAL(correct_t_vals[corner], t, ERR_TOL);    
00207           }
00208       }
00209   }
00210 
00213   void TransformedTriangleTest::test_calcStableC_Consistency()
00214   {
00215 
00216     typedef TransformedTriangle::TriSegment TriSegment;
00217     typedef TransformedTriangle::TetraCorner TetraCorner;
00218 
00219     // grandy, eq 14
00220     double correct_c_vals[24] = 
00221       { 
00222         p2[0] * q2[1] - p2[1] * q2[0], 
00223         p2[1] * q2[2] - p2[2] * q2[1], 
00224         p2[2] * q2[0] - p2[0] * q2[2],
00225         p2[0] * hq2 - hp2 * q2[0],
00226         p2[1] * hq2 - hp2 * q2[1],
00227         p2[2] * hq2 - hp2 * q2[2],
00228         Hp2 * q2[0] - p2[0] * Hq2,
00229         p2[1] * Hq2 - Hp2 * q2[1],
00230         q2[0] * r2[1] - q2[1] * r2[0], 
00231         q2[1] * r2[2] - q2[2] * r2[1], 
00232         q2[2] * r2[0] - q2[0] * r2[2],
00233         q2[0] * hr2 - hq2 * r2[0],
00234         q2[1] * hr2 - hq2 * r2[1],
00235         q2[2] * hr2 - hq2 * r2[2],
00236         Hq2 * r2[0] - q2[0] * Hr2,
00237         q2[1] * Hr2 - Hq2 * r2[1],
00238         r2[0]*p2[1]-r2[1]*p2[0], 
00239         r2[1]*p2[2]-r2[2]*p2[1], 
00240         r2[2]*p2[0]-r2[0]*p2[2],
00241         r2[0] * hp2 - hr2 * p2[0],
00242         r2[1] * hp2 - hr2 * p2[1],
00243         r2[2] * hp2 - hr2 * p2[2],
00244         Hr2 * p2[0] - r2[0] * Hp2,
00245         r2[1] * Hp2 - Hr2 * p2[1]
00246       };
00247 
00248 
00249     // number of inconsistent cases found : 
00250     // should be (at least) 1 for the test to be meaningful
00251     int num_cases = 0; 
00252 
00253     // find unstable products to check for consistency (Grandy [46])  
00254     for(TriSegment seg = TransformedTriangle::PQ ; seg <= TransformedTriangle::RP ; seg = TriSegment(seg + 1)) 
00255       {
00256         const double c_xy = tri2->calcUnstableC(seg, TransformedTriangle::C_XY);
00257         const double c_yz = tri2->calcUnstableC(seg, TransformedTriangle::C_YZ);
00258         const double c_zx = tri2->calcUnstableC(seg, TransformedTriangle::C_ZX);
00259         const double c_xh = tri2->calcUnstableC(seg, TransformedTriangle::C_XH);
00260         const double c_yh = tri2->calcUnstableC(seg, TransformedTriangle::C_YH);
00261         const double c_zh = tri2->calcUnstableC(seg, TransformedTriangle::C_ZH);
00262       
00263         const int num_zero = (c_yz*c_xh == 0.0 ? 1 : 0) + (c_zx*c_yh == 0.0 ? 1 : 0) + (c_xy*c_zh == 0.0 ? 1 : 0);
00264         const int num_neg = (c_yz*c_xh < 0.0 ? 1 : 0) + (c_zx*c_yh < 0.0 ? 1 : 0) + (c_xy*c_zh < 0.0 ? 1 : 0);
00265       
00266         if((num_zero == 1 && num_neg != 1) || num_zero == 2 || num_neg == 0 && num_zero !=3 || num_neg == 3 )
00267           {
00268             ++num_cases;
00269   
00270             double min_dist = -1.0; // initialised first time through loop
00271             TetraCorner min_corner = TransformedTriangle::O;
00272   
00273             for(TetraCorner corner = TransformedTriangle::O ; corner <= TransformedTriangle::Z ; corner = TetraCorner(corner + 1))
00274               {
00275                 // calculate distance from each corner of tetraeder to the segment
00276                 // formula : ( (Q-P) x (P - corner) )^2 / norm(Q-P)^2
00277       
00278                 const double ptP[3] = { tri2->_coords[5*seg], tri2->_coords[5*seg + 1], tri2->_coords[5*seg + 2] };
00279                 const double ptQ[3] = { tri2->_coords[5*( (seg+1) % 3)], tri2->_coords[5*( (seg+1) % 3) + 1], tri2->_coords[5*( (seg+1) % 3) + 2] };
00280                 const double ptCorner[3] = { 
00281                   corner == TransformedTriangle::X ? 1.0 : 0.0,
00282                   corner == TransformedTriangle::Y ? 1.0 : 0.0,
00283                   corner == TransformedTriangle::Z ? 1.0 : 0.0,
00284                 };
00285 
00286                 const double diff_21[3] = { ptQ[0] - ptP[0], ptQ[1] - ptP[1], ptQ[2] - ptP[2] };
00287                 const double diff_1_corner[3] = { ptP[0] - ptCorner[0], ptP[1] - ptCorner[1], ptP[2] - ptCorner[2] };
00288       
00289                 const double cross[3] = { 
00290                   diff_21[1]*diff_1_corner[2] - diff_21[2]*diff_1_corner[1],  
00291                   diff_21[2]*diff_1_corner[0] - diff_21[0]*diff_1_corner[2],
00292                   diff_21[0]*diff_1_corner[1] - diff_21[1]*diff_1_corner[0]
00293                 };
00294 
00295                 const double cross_sq = cross[0]*cross[0] + cross[1]*cross[1] + cross[2]*cross[2];
00296 
00297                 const double norm_pq = diff_21[0]*diff_21[0] + diff_21[1]*diff_21[1] + diff_21[2]*diff_21[2];
00298 
00299                 if(corner == TransformedTriangle::O || (cross_sq / norm_pq) < min_dist)
00300                   {
00301                     min_dist = cross_sq / norm_pq;
00302                     min_corner = corner;
00303                   }
00304               }
00305   
00306             // now check if the corresponding double products are zero
00307             static const DoubleProduct DOUBLE_PRODUCTS[12] = 
00308               {
00309                 TransformedTriangle::C_YZ, TransformedTriangle::C_XY, TransformedTriangle::C_ZX, // O
00310                 TransformedTriangle::C_ZH, TransformedTriangle::C_YZ, TransformedTriangle::C_YH, // X
00311                 TransformedTriangle::C_ZH, TransformedTriangle::C_ZX, TransformedTriangle::C_XH, // Y
00312                 TransformedTriangle::C_XY, TransformedTriangle::C_YH, TransformedTriangle::C_XH  // Z
00313               };
00314 
00315             for(int i = 0; i < 3 ; ++i) 
00316               {
00317                 DoubleProduct dp = DOUBLE_PRODUCTS[3*min_corner + i];
00318                 //        std::cout << std::endl << "in test inconsistent (seg,dp) :(" << seg <<", " << dp << ")" << std::endl;
00319                 CPPUNIT_ASSERT_EQUAL(0.0, tri2->calcStableC(seg, dp));
00320                 correct_c_vals[8*seg + dp] = 0.0;
00321               }
00322           }
00323 
00324       }
00325   
00326     if(num_cases < 1)
00327       {
00328         CPPUNIT_FAIL("Consistency test not pertinent");
00329       }
00330 
00331     //  std::cout << std::endl << "Number of geometric inconsistencies : " << num_cases << std::endl; 
00332     
00333     // check that all other double products have right value too
00334     double c_vals[8*3];
00335 
00336     for(TriSegment seg = TransformedTriangle::PQ ; seg <= TransformedTriangle::RP ; seg = TriSegment(seg + 1)) {
00337       
00338       c_vals[seg*8 + 0] = tri2->calcStableC(seg, TransformedTriangle::C_XY);
00339       c_vals[seg*8 + 1] = tri2->calcStableC(seg, TransformedTriangle::C_YZ);
00340       c_vals[seg*8 + 2] = tri2->calcStableC(seg, TransformedTriangle::C_ZX);
00341       c_vals[seg*8 + 3] = tri2->calcStableC(seg, TransformedTriangle::C_XH);
00342       c_vals[seg*8 + 4] = tri2->calcStableC(seg, TransformedTriangle::C_YH);
00343       c_vals[seg*8 + 5] = tri2->calcStableC(seg, TransformedTriangle::C_ZH);
00344       c_vals[seg*8 + 6] = tri2->calcStableC(seg, TransformedTriangle::C_01);
00345       c_vals[seg*8 + 7] = tri2->calcStableC(seg, TransformedTriangle::C_10);
00346 
00347     }
00348 
00349     for(int i = 0 ; i < 24 ; ++i)
00350       {
00351         CPPUNIT_ASSERT_DOUBLES_EQUAL(correct_c_vals[i], c_vals[i], ERR_TOL);
00352       }
00353   }
00354 
00355 }