src/INTERP_KERNELTest/TransformedTriangleTest.cxx

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