medtool/src/INTERP_KERNELTest/TransformedTriangleTest.cxx

   1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "TransformedTriangleTest.hxx"
  21
  22 #include <iostream>
  23
  24 using namespace INTERP_KERNEL;
  25
  26 namespace INTERP_TEST
  27 {
  28
  29   /**
  30    * Creates the TransformedTriangle objects used by the tests.
  31    *
  32    */
  33   void TransformedTriangleTest::setUp()
  34   {
  35     // tri1 -> no unstable double products - no changes brought about by preCalculateDoubleProducts
  36     //         this allows the testing of calcUnstableT
  37     // tri2 -> unstable double products - for testing calcStableC / preCalculateDoubleProducts
  38
  39     // triangle to test unstable C and T calculations
  40     p1[0] = -1.5 ; p1[1] = 0.5; p1[2] = 0.5;
  41     q1[0] = 2.0 ; q1[1] = 0.4; q1[2] = 0.6;
  42     r1[0] = 1.0 ; r1[1] = 2.4; r1[2] = 1.2;
  43     hp1 = 1 - p1[0] - p1[1] - p1[2];
  44     hq1 = 1 - q1[0] - q1[1] - q1[2];
  45     hr1 = 1 - r1[0] - r1[1] - r1[2];
  46     Hp1 = 1 - p1[0] - p1[1];
  47     Hq1 = 1 - q1[0] - q1[1];
  48     Hr1 = 1 - r1[0] - r1[1];
  49
  50     //  std::cout <<std::endl<< "constructing tri1..." << std::endl;
  51     tri1 = new TransformedTriangle(p1, q1, r1);
  52
  53
  54     // triangle to test stable C calculation
  55     const double err = 1.5e-3;
  56
  57     p2[0] = 0.000000000084654984189118; p2[1] = -0.000000000000000027536546231654231688873; p2[2] = 0.0000000000000001649875466831349431;
  58     q2[0] = -p2[0] +err; q2[1] = -p2[1] + err; q2[2] = -p2[2] +err;
  59     r2[0] = 2.01 ; r2[1] = 1.8; r2[2] = 0.92;
  60
  61     hp2 = 1 - p2[0] - p2[1] - p2[2];
  62     hq2 = 1 - q2[0] - q2[1] - q2[2];
  63     hr2 = 1 - r2[0] - r2[1] - r2[2];
  64     Hp2 = 1 - p2[0] - p2[1];
  65     Hq2 = 1 - q2[0] - q2[1];
  66     Hr2 = 1 - r2[0] - r2[1];
  67
  68     tri2 = new TransformedTriangle(p2, q2, r2);
  69
  70
  71
  72   }
  73
  74   /**
  75    * Liberates the transformed triangle objects used by the test suite
  76    *
  77    */
  78   void TransformedTriangleTest::tearDown()
  79   {
  80     delete tri1;
  81     delete tri2;
  82   }
  83
  84   /// Tests that _coords has correct values after construction of object is finished
  85   /// \brief Status : pass
  86   void TransformedTriangleTest::test_constructor() {
  87     // test that _coords has correct values after constructor is called
  88
  89     double good_values1[15] =
  90       {
  91         p1[0], p1[1], p1[2], hp1, Hp1,
  92         q1[0], q1[1], q1[2], hq1, Hq1,
  93         r1[0], r1[1], r1[2], hr1, Hr1
  94       };
  95
  96     double good_values2[15] =
  97       {
  98         p2[0], p2[1], p2[2], hp2, Hp2,
  99         q2[0], q2[1], q2[2], hq2, Hq2,
 100         r2[0], r2[1], r2[2], hr2, Hr2
 101       };
 102
 103
 104     for(int i = 0 ; i < 15 ; ++i)
 105       {
 106         CPPUNIT_ASSERT_DOUBLES_EQUAL(good_values1[i], tri1->_coords[i], ERR_TOL);
 107         CPPUNIT_ASSERT_DOUBLES_EQUAL(good_values2[i], tri2->_coords[i], ERR_TOL);
 108       }
 109
 110     CPPUNIT_ASSERT_EQUAL(true, tri1->_is_double_products_calculated);
 111     CPPUNIT_ASSERT_EQUAL(true, tri2->_is_double_products_calculated);
 112   }
 113
 114   /// Tests the calculation of double products (without the corrections)
 115   /// \brief Status : pass
 116   void TransformedTriangleTest::test_calcUnstableC() {
 117     typedef TransformedTriangle::TriSegment TriSegment;
 118
 119     // test that the correct c-values are calculated
 120
 121     double correct_c_vals[24] =
 122       {
 123         p1[0] * q1[1] - p1[1] * q1[0],
 124         p1[1] * q1[2] - p1[2] * q1[1],
 125         p1[2] * q1[0] - p1[0] * q1[2],
 126         p1[0] * hq1 - hp1 * q1[0],
 127         p1[1] * hq1 - hp1 * q1[1],
 128         p1[2] * hq1 - hp1 * q1[2],
 129         Hp1 * q1[0] - p1[0] * Hq1,
 130         p1[1] * Hq1 - Hp1 * q1[1],
 131         q1[0] * r1[1] - q1[1] * r1[0],
 132         q1[1] * r1[2] - q1[2] * r1[1],
 133         q1[2] * r1[0] - q1[0] * r1[2],
 134         q1[0] * hr1 - hq1 * r1[0],
 135         q1[1] * hr1 - hq1 * r1[1],
 136         q1[2] * hr1 - hq1 * r1[2],
 137         Hq1 * r1[0] - q1[0] * Hr1,
 138         q1[1] * Hr1 - Hq1 * r1[1],
 139         r1[0]*p1[1]-r1[1]*p1[0],
 140         r1[1]*p1[2]-r1[2]*p1[1],
 141         r1[2]*p1[0]-r1[0]*p1[2],
 142         r1[0] * hp1 - hr1 * p1[0],
 143         r1[1] * hp1 - hr1 * p1[1],
 144         r1[2] * hp1 - hr1 * p1[2],
 145         Hr1 * p1[0] - r1[0] * Hp1,
 146         r1[1] * Hp1 - Hr1 * p1[1]
 147       };
 148
 149     double c_vals[3 * 8];
 150     for(TriSegment seg = TransformedTriangle::PQ ; seg <= TransformedTriangle::RP ; seg = TriSegment(seg + 1)) {
 151
 152       c_vals[seg*8 + 0] = tri1->calcUnstableC(seg, TransformedTriangle::C_XY);
 153       c_vals[seg*8 + 1] = tri1->calcUnstableC(seg, TransformedTriangle::C_YZ);
 154       c_vals[seg*8 + 2] = tri1->calcUnstableC(seg, TransformedTriangle::C_ZX);
 155       c_vals[seg*8 + 3] = tri1->calcUnstableC(seg, TransformedTriangle::C_XH);
 156       c_vals[seg*8 + 4] = tri1->calcUnstableC(seg, TransformedTriangle::C_YH);
 157       c_vals[seg*8 + 5] = tri1->calcUnstableC(seg, TransformedTriangle::C_ZH);
 158       c_vals[seg*8 + 6] = tri1->calcUnstableC(seg, TransformedTriangle::C_01);
 159       c_vals[seg*8 + 7] = tri1->calcUnstableC(seg, TransformedTriangle::C_10);
 160
 161     }
 162
 163     for(int i = 0 ; i < 3*8 ; ++i) {
 164       CPPUNIT_ASSERT_DOUBLES_EQUAL( correct_c_vals[i], c_vals[i], ERR_TOL );
 165     }
 166
 167
 168   }
 169
 170   /// Tests the calculation of triple products (without corrections)
 171   /// \brief Status : pass
 172   void TransformedTriangleTest::test_calcUnstableT()
 173   {
 174     typedef TransformedTriangle::TetraCorner TetraCorner;
 175
 176     // correct values calculated by determinants (Grandy, [15])
 177     const double correct_t_vals[4] =
 178       {
 179         p1[0]*(q1[1]*r1[2] - q1[2]*r1[1]) -
 180         q1[0]*(p1[1]*r1[2] - p1[2]*r1[1]) +
 181         r1[0]*(p1[1]*q1[2] - p1[2]*q1[1]),
 182
 183         -(hp1*(q1[1]*r1[2] - q1[2]*r1[1]) -
 184           hq1*(p1[1]*r1[2] - p1[2]*r1[1]) +
 185           hr1*(p1[1]*q1[2] - p1[2]*q1[1])),
 186
 187         -(p1[0]*(hq1*r1[2] - q1[2]*hr1) -
 188           q1[0]*(hp1*r1[2] - p1[2]*hr1) +
 189           r1[0]*(hp1*q1[2] - p1[2]*hq1)),
 190
 191         -(p1[0]*(q1[1]*hr1 - r1[1]*hq1) -
 192           q1[0]*(p1[1]*hr1 - r1[1]*hp1) +
 193           r1[0]*(p1[1]*hq1 - q1[1]*hp1))
 194       };
 195
 196
 197     // test that triple products are correctly calculated
 198     for(TetraCorner corner = TransformedTriangle::O ; corner <= TransformedTriangle::Z ; corner = TetraCorner(corner + 1))
 199       {
 200
 201         for(int row = 1 ; row < 4 ; ++row)
 202           {
 203             const double t = tri1->calcTByDevelopingRow(corner, row, false);
 204             //    std::cout << std::endl  << " Corner = " << corner  << " Row = " << row << " got: " << t <<
 205             //  " expected: " << correct_t_vals[corner]<< std::endl;
 206             CPPUNIT_ASSERT_DOUBLES_EQUAL(correct_t_vals[corner], t, ERR_TOL);
 207           }
 208       }
 209   }
 210
 211   /// Tests the consistency correction
 212   /// \brief Status : fails because it is not significant - the consistency correction is not brought into play
 213   void TransformedTriangleTest::test_calcStableC_Consistency()
 214   {
 215
 216     typedef TransformedTriangle::TriSegment TriSegment;
 217     typedef TransformedTriangle::TetraCorner TetraCorner;
 218
 219     // grandy, eq 14
 220     double correct_c_vals[24] =
 221       {
 222         p2[0] * q2[1] - p2[1] * q2[0],
 223         p2[1] * q2[2] - p2[2] * q2[1],
 224         p2[2] * q2[0] - p2[0] * q2[2],
 225         p2[0] * hq2 - hp2 * q2[0],
 226         p2[1] * hq2 - hp2 * q2[1],
 227         p2[2] * hq2 - hp2 * q2[2],
 228         Hp2 * q2[0] - p2[0] * Hq2,
 229         p2[1] * Hq2 - Hp2 * q2[1],
 230         q2[0] * r2[1] - q2[1] * r2[0],
 231         q2[1] * r2[2] - q2[2] * r2[1],
 232         q2[2] * r2[0] - q2[0] * r2[2],
 233         q2[0] * hr2 - hq2 * r2[0],
 234         q2[1] * hr2 - hq2 * r2[1],
 235         q2[2] * hr2 - hq2 * r2[2],
 236         Hq2 * r2[0] - q2[0] * Hr2,
 237         q2[1] * Hr2 - Hq2 * r2[1],
 238         r2[0]*p2[1]-r2[1]*p2[0],
 239         r2[1]*p2[2]-r2[2]*p2[1],
 240         r2[2]*p2[0]-r2[0]*p2[2],
 241         r2[0] * hp2 - hr2 * p2[0],
 242         r2[1] * hp2 - hr2 * p2[1],
 243         r2[2] * hp2 - hr2 * p2[2],
 244         Hr2 * p2[0] - r2[0] * Hp2,
 245         r2[1] * Hp2 - Hr2 * p2[1]
 246       };
 247
 248
 249     // number of inconsistent cases found :
 250     // should be (at least) 1 for the test to be meaningful
 251     int num_cases = 0;
 252
 253     // find unstable products to check for consistency (Grandy [46])
 254     for(TriSegment seg = TransformedTriangle::PQ ; seg <= TransformedTriangle::RP ; seg = TriSegment(seg + 1))
 255       {
 256         const double c_xy = tri2->calcUnstableC(seg, TransformedTriangle::C_XY);
 257         const double c_yz = tri2->calcUnstableC(seg, TransformedTriangle::C_YZ);
 258         const double c_zx = tri2->calcUnstableC(seg, TransformedTriangle::C_ZX);
 259         const double c_xh = tri2->calcUnstableC(seg, TransformedTriangle::C_XH);
 260         const double c_yh = tri2->calcUnstableC(seg, TransformedTriangle::C_YH);
 261         const double c_zh = tri2->calcUnstableC(seg, TransformedTriangle::C_ZH);
 262
 263         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);
 264         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);
 265
 266         if((num_zero == 1 && num_neg != 1) || num_zero == 2 || num_neg == 0 && num_zero !=3 || num_neg == 3 )
 267           {
 268             ++num_cases;
 269
 270             double min_dist = -1.0; // initialised first time through loop
 271             TetraCorner min_corner = TransformedTriangle::O;
 272
 273             for(TetraCorner corner = TransformedTriangle::O ; corner <= TransformedTriangle::Z ; corner = TetraCorner(corner + 1))
 274               {
 275                 // calculate distance from each corner of tetraeder to the segment
 276                 // formula : ( (Q-P) x (P - corner) )^2 / norm(Q-P)^2
 277
 278                 const double ptP[3] = { tri2->_coords[5*seg], tri2->_coords[5*seg + 1], tri2->_coords[5*seg + 2] };
 279                 const double ptQ[3] = { tri2->_coords[5*( (seg+1) % 3)], tri2->_coords[5*( (seg+1) % 3) + 1], tri2->_coords[5*( (seg+1) % 3) + 2] };
 280                 const double ptCorner[3] = {
 281                   corner == TransformedTriangle::X ? 1.0 : 0.0,
 282                   corner == TransformedTriangle::Y ? 1.0 : 0.0,
 283                   corner == TransformedTriangle::Z ? 1.0 : 0.0,
 284                 };
 285
 286                 const double diff_21[3] = { ptQ[0] - ptP[0], ptQ[1] - ptP[1], ptQ[2] - ptP[2] };
 287                 const double diff_1_corner[3] = { ptP[0] - ptCorner[0], ptP[1] - ptCorner[1], ptP[2] - ptCorner[2] };
 288
 289                 const double cross[3] = {
 290                   diff_21[1]*diff_1_corner[2] - diff_21[2]*diff_1_corner[1],
 291                   diff_21[2]*diff_1_corner[0] - diff_21[0]*diff_1_corner[2],
 292                   diff_21[0]*diff_1_corner[1] - diff_21[1]*diff_1_corner[0]
 293                 };
 294
 295                 const double cross_sq = cross[0]*cross[0] + cross[1]*cross[1] + cross[2]*cross[2];
 296
 297                 const double norm_pq = diff_21[0]*diff_21[0] + diff_21[1]*diff_21[1] + diff_21[2]*diff_21[2];
 298
 299                 if(corner == TransformedTriangle::O || (cross_sq / norm_pq) < min_dist)
 300                   {
 301                     min_dist = cross_sq / norm_pq;
 302                     min_corner = corner;
 303                   }
 304               }
 305
 306             // now check if the corresponding double products are zero
 307             static const DoubleProduct DOUBLE_PRODUCTS[12] =
 308               {
 309                 TransformedTriangle::C_YZ, TransformedTriangle::C_XY, TransformedTriangle::C_ZX, // O
 310                 TransformedTriangle::C_ZH, TransformedTriangle::C_YZ, TransformedTriangle::C_YH, // X
 311                 TransformedTriangle::C_ZH, TransformedTriangle::C_ZX, TransformedTriangle::C_XH, // Y
 312                 TransformedTriangle::C_XY, TransformedTriangle::C_YH, TransformedTriangle::C_XH  // Z
 313               };
 314
 315             for(int i = 0; i < 3 ; ++i)
 316               {
 317                 DoubleProduct dp = DOUBLE_PRODUCTS[3*min_corner + i];
 318                 //        std::cout << std::endl << "in test inconsistent (seg,dp) :(" << seg <<", " << dp << ")" << std::endl;
 319                 CPPUNIT_ASSERT_EQUAL(0.0, tri2->calcStableC(seg, dp));
 320                 correct_c_vals[8*seg + dp] = 0.0;
 321               }
 322           }
 323
 324       }
 325
 326     if(num_cases < 1)
 327       {
 328         CPPUNIT_FAIL("Consistency test not pertinent");
 329       }
 330
 331     //  std::cout << std::endl << "Number of geometric inconsistencies : " << num_cases << std::endl;
 332
 333     // check that all other double products have right value too
 334     double c_vals[8*3];
 335
 336     for(TriSegment seg = TransformedTriangle::PQ ; seg <= TransformedTriangle::RP ; seg = TriSegment(seg + 1)) {
 337
 338       c_vals[seg*8 + 0] = tri2->calcStableC(seg, TransformedTriangle::C_XY);
 339       c_vals[seg*8 + 1] = tri2->calcStableC(seg, TransformedTriangle::C_YZ);
 340       c_vals[seg*8 + 2] = tri2->calcStableC(seg, TransformedTriangle::C_ZX);
 341       c_vals[seg*8 + 3] = tri2->calcStableC(seg, TransformedTriangle::C_XH);
 342       c_vals[seg*8 + 4] = tri2->calcStableC(seg, TransformedTriangle::C_YH);
 343       c_vals[seg*8 + 5] = tri2->calcStableC(seg, TransformedTriangle::C_ZH);
 344       c_vals[seg*8 + 6] = tri2->calcStableC(seg, TransformedTriangle::C_01);
 345       c_vals[seg*8 + 7] = tri2->calcStableC(seg, TransformedTriangle::C_10);
 346
 347     }
 348
 349     for(int i = 0 ; i < 24 ; ++i)
 350       {
 351         CPPUNIT_ASSERT_DOUBLES_EQUAL(correct_c_vals[i], c_vals[i], ERR_TOL);
 352       }
 353   }
 354
 355 }