src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF
   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 // File      : UnitTetraIntersectionBaryTest.cxx
  21 // Created   : Thu Dec 11 15:54:41 2008
  22 // Author    : Edward AGAPOV (eap)
  23 //
  24 #include "UnitTetraIntersectionBaryTest.hxx"
  25
  26 #include "UnitTetraIntersectionBary.hxx"
  27 #include "TetraAffineTransform.hxx"
  28 #include "InterpolationUtils.hxx"
  29 #include "SplitterTetra.txx"
  30 #include "MCIdType.hxx"
  31
  32 #include <iostream>
  33
  34 using namespace INTERP_KERNEL;
  35
  36 namespace INTERP_TEST
  37 {
  38   void fill_UnitTetraIntersectionBary(UnitTetraIntersectionBary& bary, double nodes[][3])
  39   {
  40     int faceConn[4][3] = { { 0, 1, 2 },// inverse order
  41                            { 0, 3, 1 },
  42                            { 1, 3, 2 },
  43                            { 3, 0, 2 } };
  44 //     int faceConn[4][3] = { { 0, 2, 1 },
  45 //                            { 0, 1, 3 },
  46 //                            { 1, 2, 3 },
  47 //                            { 3, 2, 0 } };
  48     bary.init(true);
  49     for ( int i = 0; i < 4; ++i ) {
  50       int* faceNodes = faceConn[ i ];
  51       TransformedTriangle tri(nodes[faceNodes[0]], nodes[faceNodes[1]], nodes[faceNodes[2]]);
  52       tri.calculateIntersectionVolume();
  53       bary.addSide( tri );
  54     }
  55   }
  56   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_1()
  57   {
  58     // cutting tetra coincides with the unit one
  59     double nodes[4][3] = { { 0.0, 0.0, 0.0 },
  60                            { 1.0, 0.0, 0.0 },
  61                            { 0.0, 1.0, 0.0 },
  62                            { 0.0, 0.0, 1.0 } };
  63     UnitTetraIntersectionBary bary;
  64     fill_UnitTetraIntersectionBary(bary,nodes);
  65     double baryCenter[3];
  66     bool ok    = bary.getBary( baryCenter );
  67     double vol = bary.getVolume();
  68     CPPUNIT_ASSERT( ok );
  69     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.166667, vol, 1e-5);
  70     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[0], 1e-5);
  71     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[1], 1e-5);
  72     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[2], 1e-5);
  73   }
  74   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_2()
  75   {
  76     // cutting tetra fully include the unit one
  77     double nodes[4][3] = { {-0.1,-0.1,-0.1 },
  78                            { 1.5,-0.1,-0.1 },
  79                            {-0.1, 1.5,-0.1 },
  80                            {-0.1,-0.1, 1.5 } };
  81     UnitTetraIntersectionBary bary;
  82     fill_UnitTetraIntersectionBary(bary,nodes);
  83     double baryCenter[3];
  84     bool ok    = bary.getBary( baryCenter );
  85     double vol = bary.getVolume();
  86     CPPUNIT_ASSERT( ok );
  87     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.166667, vol, 1e-5);
  88     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[0], 1e-5);
  89     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[1], 1e-5);
  90     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[2], 1e-5);
  91   }
  92   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_3()
  93   {
  94     // cutting tetra is same as the unit one but moved up by 0.5
  95     double nodes[4][3] = { { 0.0, 0.0, 0.5 },
  96                            { 1.0, 0.0, 0.5 },
  97                            { 0.0, 1.0, 0.5 },
  98                            { 0.0, 0.0, 1.5 } };
  99     UnitTetraIntersectionBary bary;
 100     fill_UnitTetraIntersectionBary(bary,nodes);
 101     double baryCenter[3];
 102     bool ok    = bary.getBary( baryCenter );
 103     double vol = bary.getVolume();
 104     CPPUNIT_ASSERT( ok );
 105     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.020833333333333332, vol, 1e-5);
 106     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125, baryCenter[0], 1e-5);
 107     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125, baryCenter[1], 1e-5);
 108     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.625, baryCenter[2], 1e-5);
 109   }
 110   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_4()
 111   {
 112     // same as previous but no cutting sides lay on the sides of unit tetra
 113     double nodes[4][3] = { {-0.2,-0.2, 0.5 },
 114                            { 1.0, 0.0, 0.5 },
 115                            { 0.0, 1.0, 0.5 },
 116                            { 0.0, 0.0, 2.0 } };
 117     UnitTetraIntersectionBary bary;
 118     fill_UnitTetraIntersectionBary(bary,nodes);
 119     double baryCenter[3];
 120     bool ok    = bary.getBary( baryCenter );
 121     double vol = bary.getVolume();
 122     CPPUNIT_ASSERT( ok );
 123     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.020833333333333332, vol, 1e-5);
 124     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125, baryCenter[0], 1e-5);
 125     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125, baryCenter[1], 1e-5);
 126     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.625, baryCenter[2], 1e-5);
 127   }
 128   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_5()
 129   {
 130     // cutting tetra is similar and parallel to the UT but moved (-0.1,-0.1,-0.1)
 131     double nodes[4][3] = { {-0.1,-0.1,-0.1 },
 132                            { 1.1,-0.1,-0.1 },
 133                            {-0.1, 1.1,-0.1 },
 134                            {-0.1,-0.1, 1.1 } };
 135     UnitTetraIntersectionBary bary;
 136     fill_UnitTetraIntersectionBary(bary,nodes);
 137     double baryCenter[3];
 138     bool ok    = bary.getBary( baryCenter );
 139     double vol = bary.getVolume();
 140     CPPUNIT_ASSERT( ok );
 141     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.1215, vol, 1e-5);
 142     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.225, baryCenter[0], 1e-5);
 143     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.225, baryCenter[1], 1e-5);
 144     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.225, baryCenter[2], 1e-5);
 145   }
 146   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_6()
 147   {
 148     // cutting tetra is deeped into the UT with one corner
 149     double nodes[4][3] = { { 0.2, 0.2, 0.2 },
 150                            { 1.0, 0.2, 0.2 },
 151                            { 0.9, 1.0, 0.2 },
 152                            { 0.9, 9.0, 1.0 } };
 153     UnitTetraIntersectionBary bary;
 154     fill_UnitTetraIntersectionBary(bary,nodes);
 155     double baryCenter[3];
 156     bool ok    = bary.getBary( baryCenter );
 157     double vol = bary.getVolume();
 158     CPPUNIT_ASSERT( ok );
 159     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.000441855, vol, 1e-5);
 160     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.353463 , baryCenter[0], 1e-5 );
 161     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.33877  , baryCenter[1], 1e-5 );
 162     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.207767 , baryCenter[2], 1e-5 );
 163   }
 164   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_7()
 165   {
 166     // cutting tetra passes through the UT with one corner
 167     double nodes[4][3] = { {-0.2, 0.2, 0.2 },
 168                            { 1.0, 0.2, 0.2 },
 169                            { 0.9, 1.0, 0.2 },
 170                            { 0.9, 0.9, 1.0 } };
 171     UnitTetraIntersectionBary bary;
 172     fill_UnitTetraIntersectionBary(bary,nodes);
 173     double baryCenter[3];
 174     bool ok    = bary.getBary( baryCenter );
 175     double vol = bary.getVolume();
 176     CPPUNIT_ASSERT( ok );
 177     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0103501, vol, 1e-5);
 178     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.215578 , baryCenter[0], 1e-5 );
 179     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.341363 , baryCenter[1], 1e-5 );
 180     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.263903 , baryCenter[2], 1e-5 );
 181   }
 182   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_8()
 183   {
 184     // cutting tetra passes through the UT with one edge
 185     double nodes[4][3] = { { 0.5, 0.2, -0.2 }, // O
 186                            {-0.5,-0.2, -0.2 }, // OX
 187                            { 1.0,-0.5, -0.2 }, // OY
 188                            { 0.5, 0.2,  1.5 } };//OZ
 189     UnitTetraIntersectionBary bary;
 190     fill_UnitTetraIntersectionBary(bary,nodes);
 191     double baryCenter[3];
 192     bool ok    = bary.getBary( baryCenter );
 193     double vol = bary.getVolume();
 194     CPPUNIT_ASSERT( ok );
 195     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0349217, vol, 1e-5);
 196     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.332275  , baryCenter[0], 1e-2 );
 197     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0565892 , baryCenter[1], 1e-3 );
 198     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.308713  , baryCenter[2], 1e-2 );
 199   }
 200   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_9()
 201   {
 202     // cutting tetra touches the UT at an edge, intersection volume == 0
 203     double nodes[4][3] = { { 1.0, 0.0, 0.0 }, // 0
 204                            {-1.0, 2.0, 2.0 }, // OX
 205                            {-1.0,-2.0, 2.0 }, // OY
 206                            { 1.0, 0.0, 2.0 } };//OZ
 207     UnitTetraIntersectionBary bary;
 208     fill_UnitTetraIntersectionBary(bary,nodes);
 209     double baryCenter[3];
 210     bool ok    = bary.getBary( baryCenter );
 211     double vol = bary.getVolume();
 212     CPPUNIT_ASSERT( !ok );
 213     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, vol, 1e-15);
 214     CPPUNIT_ASSERT_DOUBLES_EQUAL( -1. , baryCenter[0], 1e-5 );
 215     CPPUNIT_ASSERT_DOUBLES_EQUAL( -1. , baryCenter[1], 1e-5 );
 216     CPPUNIT_ASSERT_DOUBLES_EQUAL( -1. , baryCenter[2], 1e-5 );
 217   }
 218   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_10()
 219   {
 220     // cutting tetra fully includes theUT and touches it at an edge
 221     double nodes[4][3] = { { 1.0, 0.0, 0.0 }, // 0
 222                            {-1.0,-4.0, 2.0 }, // OX
 223                            {-1.0, 4.0, 2.0 }, // OY
 224                            { 1.0, 0.0,-2.0 } };//OZ
 225     UnitTetraIntersectionBary bary;
 226     fill_UnitTetraIntersectionBary(bary,nodes);
 227     double baryCenter[3];
 228     bool ok    = bary.getBary( baryCenter );
 229     double vol = bary.getVolume();
 230     CPPUNIT_ASSERT( ok );
 231     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.166667, vol, 1e-5);
 232     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[0], 1e-5);
 233     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[1], 1e-5);
 234     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[2], 1e-5);
 235   }
 236   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_11()
 237   {
 238     // cutting tetra intersects the UT and touches it at an edge
 239     double nodes[4][3] = { { 1.0, 0.0, 0.0 }, // 0
 240                            {-1.0,-4.0, 2.0 }, // OX
 241                            {-1.0, 4.0, 2.0 }, // OY
 242                            {-1.0, 0.0,-1.0 } };//OZ
 243     UnitTetraIntersectionBary bary;
 244     fill_UnitTetraIntersectionBary(bary,nodes);
 245     double baryCenter[3];
 246     bool ok    = bary.getBary( baryCenter );
 247     double vol = bary.getVolume();
 248     CPPUNIT_ASSERT( ok );
 249     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.15873 , vol, 1e-5);
 250     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.250000, baryCenter[0], 1e-5);
 251     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.230952, baryCenter[1], 1e-5);
 252     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.260714, baryCenter[2], 1e-5);
 253   }
 254   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_12()
 255   {
 256     // cutting tetra has one corner inside the UT and one its side passes through an UT edge
 257     double nodes[4][3] = { { 0.25, 0.25, 0.25 }, // 0
 258                            { 1.75,-0.25,-0.25 }, // OX
 259                            { 0.5 , 0.25, 0.25 }, // OY
 260                            { 0.5 , 0   , 0.5  } };//OZ
 261     UnitTetraIntersectionBary bary;
 262     fill_UnitTetraIntersectionBary(bary,nodes);
 263     double baryCenter[3];
 264     bool ok    = bary.getBary( baryCenter );
 265     double vol = bary.getVolume();
 266     CPPUNIT_ASSERT( ok );
 267     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.005208 , vol, 1e-5);
 268     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.562500, baryCenter[0], 1e-5);
 269     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125000, baryCenter[1], 1e-5);
 270     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.250000, baryCenter[2], 1e-5);
 271   }
 272
 273   struct __MESH_DUMMY
 274   {
 275     typedef mcIdType MyConnType;
 276   };
 277
 278   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_13()
 279   {
 280     double T[] = {
 281       66.6666666666666714,133.333333333333343,66.6666666666666714,
 282       100,200,100,
 283       100,100,100,
 284       200,200,0 };
 285
 286     double S[] = {
 287       100,166.666666666666657,66.6666666666666714,
 288       100,150,50,
 289       75,150,75,
 290       100,100,100};
 291
 292     mcIdType conn[4] = { 0,1,2,3 };
 293
 294     const double* tnodes[4]={ T, T+3, T+6, T+9 };
 295     const double* snodes[4]={ S, S+3, S+6, S+9 };
 296
 297     __MESH_DUMMY dummyMesh;
 298     SplitterTetra<__MESH_DUMMY> src( dummyMesh, snodes, conn );
 299     double volume = src.intersectTetra( tnodes );
 300     CPPUNIT_ASSERT_DOUBLES_EQUAL(6944.4444444444443,volume,1e-9);
 301   }
 302
 303   /**
 304    * Extract from a single tetra from BoxHexa1.med and BoxHexa2.med.
 305    * Symetry test was failing.
 306    */
 307   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_14()
 308   {
 309     double S[] = {
 310       25.0, 96.0, 0.0,
 311       25.0, 120.0, 0.0,
 312       37.5, 120.0, 0.0,
 313       25.0, 120.0, 11.333333333333339255};
 314
 315     double T[] = {
 316       25.0, 90.0, 6.333333333333339255,
 317       25.0, 120.0, 6.3333333333333357018,
 318       41.6666666666666714036, 120.0, 6.3333333333333348136,
 319       25.0, 120.0, 17.6666666666666714036};
 320
 321     mcIdType conn[4] = { 0,1,2,3 };
 322     const double* tnodes[4]={ T, T+3, T+6, T+9 };
 323     const double* snodes[4]={ S, S+3, S+6, S+9 };
 324     const double refVol = 48.6591695501729;
 325
 326     __MESH_DUMMY dummyMesh;
 327     SplitterTetra<__MESH_DUMMY> src( dummyMesh, snodes, conn );
 328     double volume = src.intersectTetra( tnodes );
 329     CPPUNIT_ASSERT_DOUBLES_EQUAL(refVol,volume,1e-9);
 330
 331     // Now the other way round:
 332     SplitterTetra<__MESH_DUMMY> tgt( dummyMesh, tnodes, conn );
 333     double volume2 = tgt.intersectTetra( snodes );
 334     CPPUNIT_ASSERT_DOUBLES_EQUAL(refVol,volume2,1e-9);
 335   }
 336
 337
 338   void UnitTetraIntersectionBaryTest::test_TetraAffineTransform_reverseApply()
 339   {
 340     double nodes[12] = { -4.0, 9.0, 3.0,
 341                          11.0, 0.0, 2.0,
 342                          0.0, 0.0, 0.0,
 343                          2.0, 1.0,10.0 };
 344     //    double pSrc[3] = { -4.0, 9.0, 3.0 };
 345     double pSrc[3] = { 40., -20., 100. };
 346     double pDest[] = {1,1,1};
 347     TetraAffineTransform a(nodes);
 348     a.apply( pDest, pSrc );
 349     a.reverseApply( pDest, pDest );
 350     CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[0], pDest[0], 1e-12);
 351     CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[1], pDest[1], 1e-12);
 352     CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[2], pDest[2], 1e-12);
 353   }
 354
 355   void UnitTetraIntersectionBaryTest::test_barycentric_coords()
 356   {
 357     // compute barycentric coordinates
 358     double nodes[4][3] = { {11.0, 0.0, 2.0 },
 359                            {-4.0, 9.0, 3.0 },
 360                            { 0.0, 0.0, 0.0 },
 361                            { 6.0, 1.0,10.0 }};
 362     std::vector<const double*> n (4);
 363     n[0] = &nodes[0][0];
 364     n[1] = &nodes[1][0];
 365     n[2] = &nodes[2][0];
 366     n[3] = &nodes[3][0];
 367     double p  [3] = { 2, 2, 5 }, bc[4];
 368     barycentric_coords(n, p, bc);
 369     double bcSum = 0;
 370     double p2 [3] = { 0,0,0 };
 371     for ( int i = 0; i < 4; ++i ) {
 372       bcSum += bc[i];
 373       p2[0] += bc[i]*n[i][0];
 374       p2[1] += bc[i]*n[i][1];
 375       p2[2] += bc[i]*n[i][2];
 376     }
 377     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., bcSum, 1e-12);
 378     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[0], p2[0], 1e-12);
 379     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[1], p2[1], 1e-12);
 380     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[2], p2[2], 1e-12);
 381   }
 382
 383   /* Conventions:
 384   *   - for HEXA8, point 5 is taken to be the origin (see med file ref connec):
 385   *          0 ------ 3
 386             /|       /|
 387            / |      / |
 388           1 ------ 2  |
 389           |  |     |  |
 390           |  |     |  |
 391           |  4-----|- 7
 392           | /      | /
 393           5 ------ 6
 394    */
 395   void UnitTetraIntersectionBaryTest::test_cuboid_mapped_coords_3D()
 396   {
 397     double nodes[8][3] = { { 0.0, 2.0, 4.0 }, //0
 398                            { 0.0, 0.0, 4.0 },
 399                            { 1.0, 0.0, 4.0 },
 400                            { 1.0, 2.0, 4.0 },
 401                            { 0.0, 2.0, 0.0 }, // 4
 402                            { 0.0, 0.0, 0.0 },
 403                            { 1.0, 0.0, 0.0 },
 404                            { 1.0, 2.0, 0.0 }
 405     };
 406     // Translate cube:
 407     for (int i=0; i < 8; ++i)
 408       for (int j=0; j < 3; ++j)
 409         nodes[i][j] += 15.0;
 410
 411     std::vector<const double*> n (8);
 412     for (int i=0; i<8; i++)
 413       n[i] = &nodes[i][0];
 414
 415     {
 416         // middle point
 417         double p[3] = { 15.5, 16.0, 17.0 }, bc[3];
 418         cuboid_mapped_coords(n, p, bc);
 419         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[0], 1e-12);
 420         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[1], 1e-12);
 421         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[2], 1e-12);
 422     }
 423     {
 424       // point 1
 425       double p[3] = { 15.0, 15.0, 19.0 }, bc[3];
 426       cuboid_mapped_coords(n, p, bc);
 427       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[0], 1e-12);
 428       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[1], 1e-12);
 429       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[2], 1e-12);
 430     }
 431     {
 432       // point 7
 433       double p[3] = { 16.0, 17.0, 15.0 }, bc[3];
 434       cuboid_mapped_coords(n, p, bc);
 435       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
 436       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
 437       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[2], 1e-12);
 438     }
 439     {
 440       // point 3
 441       double p[3] = { 16.0, 17.0, 19.0 }, bc[3];
 442       cuboid_mapped_coords(n, p, bc);
 443       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
 444       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
 445       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[2], 1e-12);
 446     }
 447     {
 448       // point outside
 449       double p[3] = { 2.0, 16.0, 18.0 }, bc[3];
 450       CPPUNIT_ASSERT_THROW(cuboid_mapped_coords(n, p, bc), INTERP_KERNEL::Exception);
 451     }
 452
 453   }
 454
 455   /* Convention
 456       - for QUAD4, point 0 is taken to be the origin (again see med file ref connec):
 457
 458          1------2
 459          |      |
 460          |      |
 461          0------3
 462   */
 463   void UnitTetraIntersectionBaryTest::test_quad_mapped_coords_2D()
 464   {
 465
 466     double nodes[4][2] = { { 0.0, 0.0 },
 467                            { 0.0, 1.0 },
 468                            { 2.0, 3.0 },
 469                            { 1.0, 0.0 } };
 470
 471     // Translate quad4:
 472     for (int i=0; i < 4; ++i)
 473       for (int j=0; j < 2; ++j)
 474         nodes[i][j] += 15.0;
 475
 476     std::vector<const double*> n (4);
 477     for (int i=0; i<4; i++)
 478       n[i] = &nodes[i][0];
 479
 480     {
 481       // middle point
 482       double p[2] = { 15.75, 16.0 }, bc[2];
 483       quad_mapped_coords(n, p, bc);
 484       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[0], 1e-12);
 485       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[1], 1e-12);
 486     }
 487
 488     {
 489       // middle point of seg
 490       double p[2] = { 15.5, 15.0 }, bc[2];
 491       quad_mapped_coords(n, p, bc);
 492       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[0], 1e-12);
 493       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[1], 1e-12);
 494     }
 495
 496     {
 497       // point 1
 498       double p[2] = { 15.0, 16.0 }, bc[2];
 499       quad_mapped_coords(n, p, bc);
 500       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[0], 1e-12);
 501       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
 502     }
 503     {
 504       // point 2
 505       double p[2] = { 17.0, 18.0 }, bc[2];
 506       quad_mapped_coords(n, p, bc);
 507       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
 508       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
 509     }
 510     {
 511       // point 3
 512       double p[2] = { 16.0, 15.0 }, bc[2];
 513       quad_mapped_coords(n, p, bc);
 514       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
 515       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[1], 1e-12);
 516     }
 517     {
 518       // point outside
 519       double p[2] = { 18.0, 18.0 }, bc[2];
 520       CPPUNIT_ASSERT_THROW(quad_mapped_coords(n, p, bc), INTERP_KERNEL::Exception);
 521     }
 522   }
 523
 524
 525 }