]> SALOME platform Git repositories - modules/med.git/blob - medtool/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx
Salome HOME
Move medtool folder to MED base repository
[modules/med.git] / medtool / src / INTERP_KERNELTest / UnitTetraIntersectionBaryTest.cxx
1 // Copyright (C) 2007-2015  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 }