1 // Copyright (C) 2007-2019 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 // File : UnitTetraIntersectionBaryTest.cxx
21 // Created : Thu Dec 11 15:54:41 2008
22 // Author : Edward AGAPOV (eap)
24 #include "UnitTetraIntersectionBaryTest.hxx"
26 #include "UnitTetraIntersectionBary.hxx"
27 #include "TetraAffineTransform.hxx"
28 #include "InterpolationUtils.hxx"
29 #include "SplitterTetra.txx"
30 #include "MCIdType.hxx"
34 using namespace INTERP_KERNEL;
38 void fill_UnitTetraIntersectionBary(UnitTetraIntersectionBary& bary, double nodes[][3])
40 int faceConn[4][3] = { { 0, 1, 2 },// inverse order
44 // int faceConn[4][3] = { { 0, 2, 1 },
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();
56 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_1()
58 // cutting tetra coincides with the unit one
59 double nodes[4][3] = { { 0.0, 0.0, 0.0 },
63 UnitTetraIntersectionBary bary;
64 fill_UnitTetraIntersectionBary(bary,nodes);
66 bool ok = bary.getBary( baryCenter );
67 double vol = bary.getVolume();
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);
74 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_2()
76 // cutting tetra fully include the unit one
77 double nodes[4][3] = { {-0.1,-0.1,-0.1 },
81 UnitTetraIntersectionBary bary;
82 fill_UnitTetraIntersectionBary(bary,nodes);
84 bool ok = bary.getBary( baryCenter );
85 double vol = bary.getVolume();
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);
92 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_3()
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 },
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);
110 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_4()
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 },
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);
128 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_5()
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 },
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);
146 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_6()
148 // cutting tetra is deeped into the UT with one corner
149 double nodes[4][3] = { { 0.2, 0.2, 0.2 },
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 );
164 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_7()
166 // cutting tetra passes through the UT with one corner
167 double nodes[4][3] = { {-0.2, 0.2, 0.2 },
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 );
182 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_8()
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 );
200 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_9()
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 );
218 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_10()
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);
236 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_11()
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);
254 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_12()
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);
275 typedef mcIdType MyConnType;
278 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_13()
281 66.6666666666666714,133.333333333333343,66.6666666666666714,
287 100,166.666666666666657,66.6666666666666714,
292 mcIdType conn[4] = { 0,1,2,3 };
294 const double* tnodes[4]={ T, T+3, T+6, T+9 };
295 const double* snodes[4]={ S, S+3, S+6, S+9 };
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);
303 void UnitTetraIntersectionBaryTest::test_TetraAffineTransform_reverseApply()
305 double nodes[12] = { -4.0, 9.0, 3.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);
320 void UnitTetraIntersectionBaryTest::test_barycentric_coords()
322 // compute barycentric coordinates
323 double nodes[4][3] = { {11.0, 0.0, 2.0 },
327 std::vector<const double*> n (4);
332 double p [3] = { 2, 2, 5 }, bc[4];
333 barycentric_coords(n, p, bc);
335 double p2 [3] = { 0,0,0 };
336 for ( int i = 0; i < 4; ++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];
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);
349 * - for HEXA8, point 5 is taken to be the origin (see med file ref connec):
360 void UnitTetraIntersectionBaryTest::test_cuboid_mapped_coords_3D()
362 double nodes[8][3] = { { 0.0, 2.0, 4.0 }, //0
366 { 0.0, 2.0, 0.0 }, // 4
372 for (int i=0; i < 8; ++i)
373 for (int j=0; j < 3; ++j)
376 std::vector<const double*> n (8);
377 for (int i=0; i<8; i++)
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);
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);
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);
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);
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);
421 - for QUAD4, point 0 is taken to be the origin (again see med file ref connec):
428 void UnitTetraIntersectionBaryTest::test_quad_mapped_coords_2D()
431 double nodes[4][2] = { { 0.0, 0.0 },
437 for (int i=0; i < 4; ++i)
438 for (int j=0; j < 2; ++j)
441 std::vector<const double*> n (4);
442 for (int i=0; i<4; i++)
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);
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);
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);
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);
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);
484 double p[2] = { 18.0, 18.0 }, bc[2];
485 CPPUNIT_ASSERT_THROW(quad_mapped_coords(n, p, bc), INTERP_KERNEL::Exception);