src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx

   1 // Copyright (C) 2007-2021  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 // 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   void UnitTetraIntersectionBaryTest::test_TetraAffineTransform_reverseApply()
 304   {
 305     double nodes[12] = { -4.0, 9.0, 3.0,
 306                          11.0, 0.0, 2.0,
 307                          0.0, 0.0, 0.0,
 308                          2.0, 1.0,10.0 };
 309     //    double pSrc[3] = { -4.0, 9.0, 3.0 };
 310     double pSrc[3] = { 40., -20., 100. };
 311     double pDest[] = {1,1,1};
 312     TetraAffineTransform a(nodes);
 313     a.apply( pDest, pSrc );
 314     a.reverseApply( pDest, pDest );
 315     CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[0], pDest[0], 1e-12);
 316     CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[1], pDest[1], 1e-12);
 317     CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[2], pDest[2], 1e-12);
 318   }
 319
 320   void UnitTetraIntersectionBaryTest::test_barycentric_coords()
 321   {
 322     // compute barycentric coordinates
 323     double nodes[4][3] = { {11.0, 0.0, 2.0 },
 324                            {-4.0, 9.0, 3.0 },
 325                            { 0.0, 0.0, 0.0 },
 326                            { 6.0, 1.0,10.0 }};
 327     std::vector<const double*> n (4);
 328     n[0] = &nodes[0][0];
 329     n[1] = &nodes[1][0];
 330     n[2] = &nodes[2][0];
 331     n[3] = &nodes[3][0];
 332     double p  [3] = { 2, 2, 5 }, bc[4];
 333     barycentric_coords(n, p, bc);
 334     double bcSum = 0;
 335     double p2 [3] = { 0,0,0 };
 336     for ( int i = 0; i < 4; ++i ) {
 337       bcSum += bc[i];
 338       p2[0] += bc[i]*n[i][0];
 339       p2[1] += bc[i]*n[i][1];
 340       p2[2] += bc[i]*n[i][2];
 341     }
 342     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., bcSum, 1e-12);
 343     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[0], p2[0], 1e-12);
 344     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[1], p2[1], 1e-12);
 345     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[2], p2[2], 1e-12);
 346   }
 347
 348   /* Conventions:
 349   *   - for HEXA8, point 5 is taken to be the origin (see med file ref connec):
 350   *          0 ------ 3
 351             /|       /|
 352            / |      / |
 353           1 ------ 2  |
 354           |  |     |  |
 355           |  |     |  |
 356           |  4-----|- 7
 357           | /      | /
 358           5 ------ 6
 359    */
 360   void UnitTetraIntersectionBaryTest::test_cuboid_mapped_coords_3D()
 361   {
 362     double nodes[8][3] = { { 0.0, 2.0, 4.0 }, //0
 363                            { 0.0, 0.0, 4.0 },
 364                            { 1.0, 0.0, 4.0 },
 365                            { 1.0, 2.0, 4.0 },
 366                            { 0.0, 2.0, 0.0 }, // 4
 367                            { 0.0, 0.0, 0.0 },
 368                            { 1.0, 0.0, 0.0 },
 369                            { 1.0, 2.0, 0.0 }
 370     };
 371     // Translate cube:
 372     for (int i=0; i < 8; ++i)
 373       for (int j=0; j < 3; ++j)
 374         nodes[i][j] += 15.0;
 375
 376     std::vector<const double*> n (8);
 377     for (int i=0; i<8; i++)
 378       n[i] = &nodes[i][0];
 379
 380     {
 381         // middle point
 382         double p[3] = { 15.5, 16.0, 17.0 }, bc[3];
 383         cuboid_mapped_coords(n, p, bc);
 384         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[0], 1e-12);
 385         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[1], 1e-12);
 386         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[2], 1e-12);
 387     }
 388     {
 389       // point 1
 390       double p[3] = { 15.0, 15.0, 19.0 }, bc[3];
 391       cuboid_mapped_coords(n, p, bc);
 392       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[0], 1e-12);
 393       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[1], 1e-12);
 394       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[2], 1e-12);
 395     }
 396     {
 397       // point 7
 398       double p[3] = { 16.0, 17.0, 15.0 }, bc[3];
 399       cuboid_mapped_coords(n, p, bc);
 400       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
 401       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
 402       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[2], 1e-12);
 403     }
 404     {
 405       // point 3
 406       double p[3] = { 16.0, 17.0, 19.0 }, bc[3];
 407       cuboid_mapped_coords(n, p, bc);
 408       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
 409       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
 410       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[2], 1e-12);
 411     }
 412     {
 413       // point outside
 414       double p[3] = { 2.0, 16.0, 18.0 }, bc[3];
 415       CPPUNIT_ASSERT_THROW(cuboid_mapped_coords(n, p, bc), INTERP_KERNEL::Exception);
 416     }
 417
 418   }
 419
 420   /* Convention
 421       - for QUAD4, point 0 is taken to be the origin (again see med file ref connec):
 422
 423          1------2
 424          |      |
 425          |      |
 426          0------3
 427   */
 428   void UnitTetraIntersectionBaryTest::test_quad_mapped_coords_2D()
 429   {
 430
 431     double nodes[4][2] = { { 0.0, 0.0 },
 432                            { 0.0, 1.0 },
 433                            { 2.0, 3.0 },
 434                            { 1.0, 0.0 } };
 435
 436     // Translate quad4:
 437     for (int i=0; i < 4; ++i)
 438       for (int j=0; j < 2; ++j)
 439         nodes[i][j] += 15.0;
 440
 441     std::vector<const double*> n (4);
 442     for (int i=0; i<4; i++)
 443       n[i] = &nodes[i][0];
 444
 445     {
 446       // middle point
 447       double p[2] = { 15.75, 16.0 }, bc[2];
 448       quad_mapped_coords(n, p, bc);
 449       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[0], 1e-12);
 450       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[1], 1e-12);
 451     }
 452
 453     {
 454       // middle point of seg
 455       double p[2] = { 15.5, 15.0 }, bc[2];
 456       quad_mapped_coords(n, p, bc);
 457       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[0], 1e-12);
 458       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[1], 1e-12);
 459     }
 460
 461     {
 462       // point 1
 463       double p[2] = { 15.0, 16.0 }, bc[2];
 464       quad_mapped_coords(n, p, bc);
 465       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[0], 1e-12);
 466       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
 467     }
 468     {
 469       // point 2
 470       double p[2] = { 17.0, 18.0 }, bc[2];
 471       quad_mapped_coords(n, p, bc);
 472       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
 473       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
 474     }
 475     {
 476       // point 3
 477       double p[2] = { 16.0, 15.0 }, bc[2];
 478       quad_mapped_coords(n, p, bc);
 479       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
 480       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[1], 1e-12);
 481     }
 482     {
 483       // point outside
 484       double p[2] = { 18.0, 18.0 }, bc[2];
 485       CPPUNIT_ASSERT_THROW(quad_mapped_coords(n, p, bc), INTERP_KERNEL::Exception);
 486     }
 487   }
 488
 489
 490 }