src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx

   1 // Copyright (C) 2007-2016  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
  31 #include <iostream>
  32
  33 using namespace INTERP_KERNEL;
  34
  35 namespace INTERP_TEST
  36 {
  37   void fill_UnitTetraIntersectionBary(UnitTetraIntersectionBary& bary, double nodes[][3])
  38   {
  39     int faceConn[4][3] = { { 0, 1, 2 },// inverse order
  40                            { 0, 3, 1 },
  41                            { 1, 3, 2 },
  42                            { 3, 0, 2 } };
  43 //     int faceConn[4][3] = { { 0, 2, 1 },
  44 //                            { 0, 1, 3 },
  45 //                            { 1, 2, 3 },
  46 //                            { 3, 2, 0 } };
  47     bary.init(true);
  48     for ( int i = 0; i < 4; ++i ) {
  49       int* faceNodes = faceConn[ i ];
  50       TransformedTriangle tri(nodes[faceNodes[0]], nodes[faceNodes[1]], nodes[faceNodes[2]]);
  51       tri.calculateIntersectionVolume();
  52       bary.addSide( tri );
  53     }
  54   }
  55   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_1()
  56   {
  57     // cutting tetra coincides with the unit one
  58     double nodes[4][3] = { { 0.0, 0.0, 0.0 },
  59                            { 1.0, 0.0, 0.0 },
  60                            { 0.0, 1.0, 0.0 },
  61                            { 0.0, 0.0, 1.0 } };
  62     UnitTetraIntersectionBary bary;
  63     fill_UnitTetraIntersectionBary(bary,nodes);
  64     double baryCenter[3];
  65     bool ok    = bary.getBary( baryCenter );
  66     double vol = bary.getVolume();
  67     CPPUNIT_ASSERT( ok );
  68     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.166667, vol, 1e-5);
  69     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[0], 1e-5);
  70     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[1], 1e-5);
  71     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[2], 1e-5);
  72   }
  73   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_2()
  74   {
  75     // cutting tetra fully include the unit one
  76     double nodes[4][3] = { {-0.1,-0.1,-0.1 },
  77                            { 1.5,-0.1,-0.1 },
  78                            {-0.1, 1.5,-0.1 },
  79                            {-0.1,-0.1, 1.5 } };
  80     UnitTetraIntersectionBary bary;
  81     fill_UnitTetraIntersectionBary(bary,nodes);
  82     double baryCenter[3];
  83     bool ok    = bary.getBary( baryCenter );
  84     double vol = bary.getVolume();
  85     CPPUNIT_ASSERT( ok );
  86     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.166667, vol, 1e-5);
  87     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[0], 1e-5);
  88     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[1], 1e-5);
  89     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[2], 1e-5);
  90   }
  91   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_3()
  92   {
  93     // cutting tetra is same as the unit one but moved up by 0.5
  94     double nodes[4][3] = { { 0.0, 0.0, 0.5 },
  95                            { 1.0, 0.0, 0.5 },
  96                            { 0.0, 1.0, 0.5 },
  97                            { 0.0, 0.0, 1.5 } };
  98     UnitTetraIntersectionBary bary;
  99     fill_UnitTetraIntersectionBary(bary,nodes);
 100     double baryCenter[3];
 101     bool ok    = bary.getBary( baryCenter );
 102     double vol = bary.getVolume();
 103     CPPUNIT_ASSERT( ok );
 104     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.020833333333333332, vol, 1e-5);
 105     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125, baryCenter[0], 1e-5);
 106     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125, baryCenter[1], 1e-5);
 107     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.625, baryCenter[2], 1e-5);
 108   }
 109   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_4()
 110   {
 111     // same as previous but no cutting sides lay on the sides of unit tetra
 112     double nodes[4][3] = { {-0.2,-0.2, 0.5 },
 113                            { 1.0, 0.0, 0.5 },
 114                            { 0.0, 1.0, 0.5 },
 115                            { 0.0, 0.0, 2.0 } };
 116     UnitTetraIntersectionBary bary;
 117     fill_UnitTetraIntersectionBary(bary,nodes);
 118     double baryCenter[3];
 119     bool ok    = bary.getBary( baryCenter );
 120     double vol = bary.getVolume();
 121     CPPUNIT_ASSERT( ok );
 122     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.020833333333333332, vol, 1e-5);
 123     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125, baryCenter[0], 1e-5);
 124     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125, baryCenter[1], 1e-5);
 125     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.625, baryCenter[2], 1e-5);
 126   }
 127   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_5()
 128   {
 129     // cutting tetra is similar and parallel to the UT but moved (-0.1,-0.1,-0.1)
 130     double nodes[4][3] = { {-0.1,-0.1,-0.1 },
 131                            { 1.1,-0.1,-0.1 },
 132                            {-0.1, 1.1,-0.1 },
 133                            {-0.1,-0.1, 1.1 } };
 134     UnitTetraIntersectionBary bary;
 135     fill_UnitTetraIntersectionBary(bary,nodes);
 136     double baryCenter[3];
 137     bool ok    = bary.getBary( baryCenter );
 138     double vol = bary.getVolume();
 139     CPPUNIT_ASSERT( ok );
 140     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.1215, vol, 1e-5);
 141     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.225, baryCenter[0], 1e-5);
 142     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.225, baryCenter[1], 1e-5);
 143     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.225, baryCenter[2], 1e-5);
 144   }
 145   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_6()
 146   {
 147     // cutting tetra is deeped into the UT with one corner
 148     double nodes[4][3] = { { 0.2, 0.2, 0.2 },
 149                            { 1.0, 0.2, 0.2 },
 150                            { 0.9, 1.0, 0.2 },
 151                            { 0.9, 9.0, 1.0 } };
 152     UnitTetraIntersectionBary bary;
 153     fill_UnitTetraIntersectionBary(bary,nodes);
 154     double baryCenter[3];
 155     bool ok    = bary.getBary( baryCenter );
 156     double vol = bary.getVolume();
 157     CPPUNIT_ASSERT( ok );
 158     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.000441855, vol, 1e-5);
 159     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.353463 , baryCenter[0], 1e-5 );
 160     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.33877  , baryCenter[1], 1e-5 );
 161     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.207767 , baryCenter[2], 1e-5 );
 162   }
 163   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_7()
 164   {
 165     // cutting tetra passes through the UT with one corner
 166     double nodes[4][3] = { {-0.2, 0.2, 0.2 },
 167                            { 1.0, 0.2, 0.2 },
 168                            { 0.9, 1.0, 0.2 },
 169                            { 0.9, 0.9, 1.0 } };
 170     UnitTetraIntersectionBary bary;
 171     fill_UnitTetraIntersectionBary(bary,nodes);
 172     double baryCenter[3];
 173     bool ok    = bary.getBary( baryCenter );
 174     double vol = bary.getVolume();
 175     CPPUNIT_ASSERT( ok );
 176     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0103501, vol, 1e-5);
 177     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.215578 , baryCenter[0], 1e-5 );
 178     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.341363 , baryCenter[1], 1e-5 );
 179     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.263903 , baryCenter[2], 1e-5 );
 180   }
 181   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_8()
 182   {
 183     // cutting tetra passes through the UT with one edge
 184     double nodes[4][3] = { { 0.5, 0.2, -0.2 }, // O
 185                            {-0.5,-0.2, -0.2 }, // OX
 186                            { 1.0,-0.5, -0.2 }, // OY
 187                            { 0.5, 0.2,  1.5 } };//OZ
 188     UnitTetraIntersectionBary bary;
 189     fill_UnitTetraIntersectionBary(bary,nodes);
 190     double baryCenter[3];
 191     bool ok    = bary.getBary( baryCenter );
 192     double vol = bary.getVolume();
 193     CPPUNIT_ASSERT( ok );
 194     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0349217, vol, 1e-5);
 195     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.332275  , baryCenter[0], 1e-2 );
 196     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0565892 , baryCenter[1], 1e-3 );
 197     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.308713  , baryCenter[2], 1e-2 );
 198   }
 199   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_9()
 200   {
 201     // cutting tetra touches the UT at an edge, intersection volume == 0
 202     double nodes[4][3] = { { 1.0, 0.0, 0.0 }, // 0
 203                            {-1.0, 2.0, 2.0 }, // OX
 204                            {-1.0,-2.0, 2.0 }, // OY
 205                            { 1.0, 0.0, 2.0 } };//OZ
 206     UnitTetraIntersectionBary bary;
 207     fill_UnitTetraIntersectionBary(bary,nodes);
 208     double baryCenter[3];
 209     bool ok    = bary.getBary( baryCenter );
 210     double vol = bary.getVolume();
 211     CPPUNIT_ASSERT( !ok );
 212     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, vol, 1e-15);
 213     CPPUNIT_ASSERT_DOUBLES_EQUAL( -1. , baryCenter[0], 1e-5 );
 214     CPPUNIT_ASSERT_DOUBLES_EQUAL( -1. , baryCenter[1], 1e-5 );
 215     CPPUNIT_ASSERT_DOUBLES_EQUAL( -1. , baryCenter[2], 1e-5 );
 216   }
 217   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_10()
 218   {
 219     // cutting tetra fully includes theUT and touches it at an edge
 220     double nodes[4][3] = { { 1.0, 0.0, 0.0 }, // 0
 221                            {-1.0,-4.0, 2.0 }, // OX
 222                            {-1.0, 4.0, 2.0 }, // OY
 223                            { 1.0, 0.0,-2.0 } };//OZ
 224     UnitTetraIntersectionBary bary;
 225     fill_UnitTetraIntersectionBary(bary,nodes);
 226     double baryCenter[3];
 227     bool ok    = bary.getBary( baryCenter );
 228     double vol = bary.getVolume();
 229     CPPUNIT_ASSERT( ok );
 230     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.166667, vol, 1e-5);
 231     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[0], 1e-5);
 232     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[1], 1e-5);
 233     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[2], 1e-5);
 234   }
 235   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_11()
 236   {
 237     // cutting tetra intersects the UT and touches it at an edge
 238     double nodes[4][3] = { { 1.0, 0.0, 0.0 }, // 0
 239                            {-1.0,-4.0, 2.0 }, // OX
 240                            {-1.0, 4.0, 2.0 }, // OY
 241                            {-1.0, 0.0,-1.0 } };//OZ
 242     UnitTetraIntersectionBary bary;
 243     fill_UnitTetraIntersectionBary(bary,nodes);
 244     double baryCenter[3];
 245     bool ok    = bary.getBary( baryCenter );
 246     double vol = bary.getVolume();
 247     CPPUNIT_ASSERT( ok );
 248     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.15873 , vol, 1e-5);
 249     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.250000, baryCenter[0], 1e-5);
 250     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.230952, baryCenter[1], 1e-5);
 251     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.260714, baryCenter[2], 1e-5);
 252   }
 253   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_12()
 254   {
 255     // cutting tetra has one corner inside the UT and one its side passes through an UT edge
 256     double nodes[4][3] = { { 0.25, 0.25, 0.25 }, // 0
 257                            { 1.75,-0.25,-0.25 }, // OX
 258                            { 0.5 , 0.25, 0.25 }, // OY
 259                            { 0.5 , 0   , 0.5  } };//OZ
 260     UnitTetraIntersectionBary bary;
 261     fill_UnitTetraIntersectionBary(bary,nodes);
 262     double baryCenter[3];
 263     bool ok    = bary.getBary( baryCenter );
 264     double vol = bary.getVolume();
 265     CPPUNIT_ASSERT( ok );
 266     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.005208 , vol, 1e-5);
 267     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.562500, baryCenter[0], 1e-5);
 268     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125000, baryCenter[1], 1e-5);
 269     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.250000, baryCenter[2], 1e-5);
 270   }
 271
 272   struct __MESH_DUMMY
 273   {
 274     typedef int MyConnType;
 275   };
 276
 277   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_13()
 278   {
 279     double T[] = {
 280       66.6666666666666714,133.333333333333343,66.6666666666666714,
 281       100,200,100,
 282       100,100,100,
 283       200,200,0 };
 284
 285     double S[] = {
 286       100,166.666666666666657,66.6666666666666714,
 287       100,150,50,
 288       75,150,75,
 289       100,100,100};
 290
 291     int conn[4] = { 0,1,2,3 };
 292
 293     const double* tnodes[4]={ T, T+3, T+6, T+9 };
 294     const double* snodes[4]={ S, S+3, S+6, S+9 };
 295
 296     __MESH_DUMMY dummyMesh;
 297     SplitterTetra<__MESH_DUMMY> src( dummyMesh, snodes, conn );
 298     double volume = src.intersectTetra( tnodes );
 299     CPPUNIT_ASSERT_DOUBLES_EQUAL(6944.4444444444443,volume,1e-9);
 300   }
 301
 302   void UnitTetraIntersectionBaryTest::test_TetraAffineTransform_reverseApply()
 303   {
 304     double nodes[12] = { -4.0, 9.0, 3.0,
 305                          11.0, 0.0, 2.0,
 306                          0.0, 0.0, 0.0,
 307                          2.0, 1.0,10.0 };
 308     //    double pSrc[3] = { -4.0, 9.0, 3.0 };
 309     double pSrc[3] = { 40., -20., 100. };
 310     double pDest[] = {1,1,1};
 311     TetraAffineTransform a(nodes);
 312     a.apply( pDest, pSrc );
 313     a.reverseApply( pDest, pDest );
 314     CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[0], pDest[0], 1e-12);
 315     CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[1], pDest[1], 1e-12);
 316     CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[2], pDest[2], 1e-12);
 317   }
 318
 319   void UnitTetraIntersectionBaryTest::test_barycentric_coords()
 320   {
 321     // compute barycentric coordinates
 322     double nodes[4][3] = { {11.0, 0.0, 2.0 },
 323                            {-4.0, 9.0, 3.0 },
 324                            { 0.0, 0.0, 0.0 },
 325                            { 6.0, 1.0,10.0 }};
 326     std::vector<const double*> n (4);
 327     n[0] = &nodes[0][0];
 328     n[1] = &nodes[1][0];
 329     n[2] = &nodes[2][0];
 330     n[3] = &nodes[3][0];
 331     double p  [3] = { 2, 2, 5 }, bc[4];
 332     barycentric_coords(n, p, bc);
 333     double bcSum = 0;
 334     double p2 [3] = { 0,0,0 };
 335     for ( int i = 0; i < 4; ++i ) {
 336       bcSum += bc[i];
 337       p2[0] += bc[i]*n[i][0];
 338       p2[1] += bc[i]*n[i][1];
 339       p2[2] += bc[i]*n[i][2];
 340     }
 341     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., bcSum, 1e-12);
 342     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[0], p2[0], 1e-12);
 343     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[1], p2[1], 1e-12);
 344     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[2], p2[2], 1e-12);
 345   }
 346 }