src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.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 // File      : UnitTetraIntersectionBaryTest.cxx
  20 // Created   : Thu Dec 11 15:54:41 2008
  21 // Author    : Edward AGAPOV (eap)
  22 #include "UnitTetraIntersectionBaryTest.hxx"
  23
  24 #include "UnitTetraIntersectionBary.hxx"
  25 #include "TetraAffineTransform.hxx"
  26 #include "InterpolationUtils.hxx"
  27
  28 #include <iostream>
  29
  30 using namespace INTERP_KERNEL;
  31 using namespace std;
  32
  33 namespace INTERP_TEST
  34 {
  35   void fill_UnitTetraIntersectionBary(UnitTetraIntersectionBary& bary, double nodes[][3])
  36   {
  37     int faceConn[4][3] = { { 0, 2, 1 },
  38                            { 0, 1, 3 },
  39                            { 1, 2, 3 },
  40                            { 3, 2, 0 } };
  41     bary.init();
  42     for ( int i = 0; i < 4; ++i ) {
  43       int* faceNodes = faceConn[ i ];
  44       TransformedTriangle tri(nodes[faceNodes[0]], nodes[faceNodes[1]], nodes[faceNodes[2]]);
  45       tri.calculateIntersectionVolume();
  46       bary.addSide( tri );
  47     }
  48   }
  49   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_1()
  50   {
  51     // cutting tetra coincides with the unit one
  52     double nodes[4][3] = { { 0.0, 0.0, 0.0 },
  53                            { 1.0, 0.0, 0.0 },
  54                            { 0.0, 1.0, 0.0 },
  55                            { 0.0, 0.0, 1.0 } };
  56     UnitTetraIntersectionBary bary;
  57     fill_UnitTetraIntersectionBary(bary,nodes);
  58     double baryCenter[3];
  59     bool ok    = bary.getBary( baryCenter );
  60     double vol = bary.getVolume();
  61     CPPUNIT_ASSERT( ok );
  62     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.166667, vol, 1e-5);
  63     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[0], 1e-5);
  64     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[1], 1e-5);
  65     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[2], 1e-5);
  66   }
  67   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_2()
  68   {
  69     // cutting tetra fully include the unit one
  70     double nodes[4][3] = { {-0.1,-0.1,-0.1 },
  71                            { 1.5,-0.1,-0.1 },
  72                            {-0.1, 1.5,-0.1 },
  73                            {-0.1,-0.1, 1.5 } };
  74     UnitTetraIntersectionBary bary;
  75     fill_UnitTetraIntersectionBary(bary,nodes);
  76     double baryCenter[3];
  77     bool ok    = bary.getBary( baryCenter );
  78     double vol = bary.getVolume();
  79     CPPUNIT_ASSERT( ok );
  80     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.166667, vol, 1e-5);
  81     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[0], 1e-5);
  82     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[1], 1e-5);
  83     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[2], 1e-5);
  84   }
  85   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_3()
  86   {
  87     // cutting tetra is same as the unit one but moved up by 0.5
  88     double nodes[4][3] = { { 0.0, 0.0, 0.5 },
  89                            { 1.0, 0.0, 0.5 },
  90                            { 0.0, 1.0, 0.5 },
  91                            { 0.0, 0.0, 1.5 } };
  92     UnitTetraIntersectionBary bary;
  93     fill_UnitTetraIntersectionBary(bary,nodes);
  94     double baryCenter[3];
  95     bool ok    = bary.getBary( baryCenter );
  96     double vol = bary.getVolume();
  97     CPPUNIT_ASSERT( ok );
  98     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.020833333333333332, vol, 1e-5);
  99     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125, baryCenter[0], 1e-5);
 100     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125, baryCenter[1], 1e-5);
 101     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.625, baryCenter[2], 1e-5);
 102   }
 103   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_4()
 104   {
 105     // same as previous but no cutting sides lay on the sides of unit tetra
 106     double nodes[4][3] = { {-0.2,-0.2, 0.5 },
 107                            { 1.0, 0.0, 0.5 },
 108                            { 0.0, 1.0, 0.5 },
 109                            { 0.0, 0.0, 2.0 } };
 110     UnitTetraIntersectionBary bary;
 111     fill_UnitTetraIntersectionBary(bary,nodes);
 112     double baryCenter[3];
 113     bool ok    = bary.getBary( baryCenter );
 114     double vol = bary.getVolume();
 115     CPPUNIT_ASSERT( ok );
 116     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.020833333333333332, vol, 1e-5);
 117     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125, baryCenter[0], 1e-5);
 118     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.125, baryCenter[1], 1e-5);
 119     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.625, baryCenter[2], 1e-5);
 120   }
 121   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_5()
 122   {
 123     // cutting tetra is similar and parallel to the UT but moved (-0.1,-0.1,-0.1)
 124     double nodes[4][3] = { {-0.1,-0.1,-0.1 },
 125                            { 1.1,-0.1,-0.1 },
 126                            {-0.1, 1.1,-0.1 },
 127                            {-0.1,-0.1, 1.1 } };
 128     UnitTetraIntersectionBary bary;
 129     fill_UnitTetraIntersectionBary(bary,nodes);
 130     double baryCenter[3];
 131     bool ok    = bary.getBary( baryCenter );
 132     double vol = bary.getVolume();
 133     CPPUNIT_ASSERT( ok );
 134     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.1215, vol, 1e-5);
 135     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.225, baryCenter[0], 1e-5);
 136     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.225, baryCenter[1], 1e-5);
 137     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.225, baryCenter[2], 1e-5);
 138   }
 139   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_6()
 140   {
 141     // cutting tetra is deeped into the UT with one corner
 142     double nodes[4][3] = { { 0.2, 0.2, 0.2 },
 143                            { 1.0, 0.2, 0.2 },
 144                            { 0.9, 1.0, 0.2 },
 145                            { 0.9, 9.0, 1.0 } };
 146     UnitTetraIntersectionBary bary;
 147     fill_UnitTetraIntersectionBary(bary,nodes);
 148     double baryCenter[3];
 149     bool ok    = bary.getBary( baryCenter );
 150     double vol = bary.getVolume();
 151     CPPUNIT_ASSERT( ok );
 152     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.000441855, vol, 1e-5);
 153     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.353463 , baryCenter[0], 1e-5 );
 154     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.33877  , baryCenter[1], 1e-5 );
 155     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.207767 , baryCenter[2], 1e-5 );
 156   }
 157   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_7()
 158   {
 159     // cutting tetra passes through the UT with one corner
 160     double nodes[4][3] = { {-0.2, 0.2, 0.2 },
 161                            { 1.0, 0.2, 0.2 },
 162                            { 0.9, 1.0, 0.2 },
 163                            { 0.9, 0.9, 1.0 } };
 164     UnitTetraIntersectionBary bary;
 165     fill_UnitTetraIntersectionBary(bary,nodes);
 166     double baryCenter[3];
 167     bool ok    = bary.getBary( baryCenter );
 168     double vol = bary.getVolume();
 169     CPPUNIT_ASSERT( ok );
 170     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0103501, vol, 1e-5);
 171     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.215578 , baryCenter[0], 1e-5 );
 172     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.341363 , baryCenter[1], 1e-5 );
 173     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.263903 , baryCenter[2], 1e-5 );
 174   }
 175   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_8()
 176   {
 177     // cutting tetra passes through the UT with one edge
 178     double nodes[4][3] = { { 0.5, 0.2, -0.2 }, // O
 179                            {-0.5,-0.2, -0.2 }, // OX
 180                            { 1.0,-0.5, -0.2 }, // OY
 181                            { 0.5, 0.2,  1.5 } };//OZ
 182     UnitTetraIntersectionBary bary;
 183     fill_UnitTetraIntersectionBary(bary,nodes);
 184     double baryCenter[3];
 185     bool ok    = bary.getBary( baryCenter );
 186     double vol = bary.getVolume();
 187     CPPUNIT_ASSERT( ok );
 188     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0349217, vol, 1e-5);
 189     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.332275  , baryCenter[0], 1e-2 );
 190     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0565892 , baryCenter[1], 1e-3 );
 191     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.308713  , baryCenter[2], 1e-2 );
 192   }
 193   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_9()
 194   {
 195     // cutting tetra touches the UT at an edge, intersection volume == 0
 196     double nodes[4][3] = { { 1.0, 0.0, 0.0 }, // 0
 197                            {-1.0, 2.0, 2.0 }, // OX
 198                            {-1.0,-2.0, 2.0 }, // OY
 199                            { 1.0, 0.0, 2.0 } };//OZ
 200     UnitTetraIntersectionBary bary;
 201     fill_UnitTetraIntersectionBary(bary,nodes);
 202     double baryCenter[3];
 203     bool ok    = bary.getBary( baryCenter );
 204     double vol = bary.getVolume();
 205     CPPUNIT_ASSERT( !ok );
 206     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, vol, 1e-15);
 207     CPPUNIT_ASSERT_DOUBLES_EQUAL( -1. , baryCenter[0], 1e-5 );
 208     CPPUNIT_ASSERT_DOUBLES_EQUAL( -1. , baryCenter[1], 1e-5 );
 209     CPPUNIT_ASSERT_DOUBLES_EQUAL( -1. , baryCenter[2], 1e-5 );
 210   }
 211   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_10()
 212   {
 213     // cutting tetra fully includes theUT and touches it at an edge
 214     double nodes[4][3] = { { 1.0, 0.0, 0.0 }, // 0
 215                            {-1.0,-4.0, 2.0 }, // OX
 216                            {-1.0, 4.0, 2.0 }, // OY
 217                            { 1.0, 0.0,-2.0 } };//OZ
 218     UnitTetraIntersectionBary bary;
 219     fill_UnitTetraIntersectionBary(bary,nodes);
 220     double baryCenter[3];
 221     bool ok    = bary.getBary( baryCenter );
 222     double vol = bary.getVolume();
 223     CPPUNIT_ASSERT( ok );
 224     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.166667, vol, 1e-5);
 225     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[0], 1e-5);
 226     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[1], 1e-5);
 227     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, baryCenter[2], 1e-5);
 228   }
 229   void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_11()
 230   {
 231     // cutting tetra intersects the UT and touches it at an edge
 232     double nodes[4][3] = { { 1.0, 0.0, 0.0 }, // 0
 233                            {-1.0,-4.0, 2.0 }, // OX
 234                            {-1.0, 4.0, 2.0 }, // OY
 235                            {-1.0, 0.0,-1.0 } };//OZ
 236     UnitTetraIntersectionBary bary;
 237     fill_UnitTetraIntersectionBary(bary,nodes);
 238     double baryCenter[3];
 239     bool ok    = bary.getBary( baryCenter );
 240     double vol = bary.getVolume();
 241     CPPUNIT_ASSERT( ok );
 242     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.15873 , vol, 1e-5);
 243     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.250000, baryCenter[0], 1e-5);
 244     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.230952, baryCenter[1], 1e-5);
 245     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.260714, baryCenter[2], 1e-5);
 246   }
 247
 248   void UnitTetraIntersectionBaryTest::test_TetraAffineTransform_reverseApply()
 249   {
 250     double nodes[4][3] = { {-4.0, 9.0, 3.0 },
 251                            {11.0, 0.0, 2.0 },
 252                            { 0.0, 0.0, 0.0 },
 253                            { 2.0, 1.0,10.0 }};
 254     //    double pSrc[3] = { -4.0, 9.0, 3.0 };
 255     double pSrc[3] = { 40., -20., 100. };
 256     double pDest[] = {1,1,1};
 257     const double* n[4] = { &nodes[0][0], &nodes[1][0], &nodes[2][0], &nodes[3][0] };
 258     TetraAffineTransform a(&n[0]);
 259     a.apply( pDest, pSrc );
 260     a.reverseApply( pDest, pDest );
 261     CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[0], pDest[0], 1e-12);
 262     CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[1], pDest[1], 1e-12);
 263     CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[2], pDest[2], 1e-12);
 264   }
 265
 266   void UnitTetraIntersectionBaryTest::test_barycentric_coords()
 267   {
 268     // compute barycentric coordinates
 269     double nodes[4][3] = { {11.0, 0.0, 2.0 },
 270                            {-4.0, 9.0, 3.0 },
 271                            { 0.0, 0.0, 0.0 },
 272                            { 6.0, 1.0,10.0 }};
 273     vector<const double*> n (4);
 274     n[0] = &nodes[0][0];
 275     n[1] = &nodes[1][0];
 276     n[2] = &nodes[2][0];
 277     n[3] = &nodes[3][0];
 278     double p  [3] = { 2, 2, 5 }, bc[4];
 279     barycentric_coords(n, p, bc);
 280     double bcSum = 0;
 281     double p2 [3] = { 0,0,0 };
 282     for ( int i = 0; i < 4; ++i ) {
 283       bcSum += bc[i];
 284       p2[0] += bc[i]*n[i][0];
 285       p2[1] += bc[i]*n[i][1];
 286       p2[2] += bc[i]*n[i][2];
 287     }
 288     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., bcSum, 1e-12);
 289     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[0], p2[0], 1e-12);
 290     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[1], p2[1], 1e-12);
 291     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[2], p2[2], 1e-12);
 292   }
 293 }