1 // Copyright (C) 2007-2024 CEA, EDF
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);
304 * Extract from a single tetra from BoxHexa1.med and BoxHexa2.med.
305 * Symetry test was failing.
307 void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_14()
313 25.0, 120.0, 11.333333333333339255};
316 25.0, 90.0, 6.333333333333339255,
317 25.0, 120.0, 6.3333333333333357018,
318 41.6666666666666714036, 120.0, 6.3333333333333348136,
319 25.0, 120.0, 17.6666666666666714036};
321 mcIdType conn[4] = { 0,1,2,3 };
322 const double* tnodes[4]={ T, T+3, T+6, T+9 };
323 const double* snodes[4]={ S, S+3, S+6, S+9 };
324 const double refVol = 48.6591695501729;
326 __MESH_DUMMY dummyMesh;
327 SplitterTetra<__MESH_DUMMY> src( dummyMesh, snodes, conn );
328 double volume = src.intersectTetra( tnodes );
329 CPPUNIT_ASSERT_DOUBLES_EQUAL(refVol,volume,1e-9);
331 // Now the other way round:
332 SplitterTetra<__MESH_DUMMY> tgt( dummyMesh, tnodes, conn );
333 double volume2 = tgt.intersectTetra( snodes );
334 CPPUNIT_ASSERT_DOUBLES_EQUAL(refVol,volume2,1e-9);
338 void UnitTetraIntersectionBaryTest::test_TetraAffineTransform_reverseApply()
340 double nodes[12] = { -4.0, 9.0, 3.0,
344 // double pSrc[3] = { -4.0, 9.0, 3.0 };
345 double pSrc[3] = { 40., -20., 100. };
346 double pDest[] = {1,1,1};
347 TetraAffineTransform a(nodes);
348 a.apply( pDest, pSrc );
349 a.reverseApply( pDest, pDest );
350 CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[0], pDest[0], 1e-12);
351 CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[1], pDest[1], 1e-12);
352 CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[2], pDest[2], 1e-12);
355 void UnitTetraIntersectionBaryTest::test_barycentric_coords()
357 // compute barycentric coordinates
358 double nodes[4][3] = { {11.0, 0.0, 2.0 },
362 std::vector<const double*> n (4);
367 double p [3] = { 2, 2, 5 }, bc[4];
368 barycentric_coords(n, p, bc);
370 double p2 [3] = { 0,0,0 };
371 for ( int i = 0; i < 4; ++i ) {
373 p2[0] += bc[i]*n[i][0];
374 p2[1] += bc[i]*n[i][1];
375 p2[2] += bc[i]*n[i][2];
377 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., bcSum, 1e-12);
378 CPPUNIT_ASSERT_DOUBLES_EQUAL( p[0], p2[0], 1e-12);
379 CPPUNIT_ASSERT_DOUBLES_EQUAL( p[1], p2[1], 1e-12);
380 CPPUNIT_ASSERT_DOUBLES_EQUAL( p[2], p2[2], 1e-12);
384 * - for HEXA8, point 5 is taken to be the origin (see med file ref connec):
395 void UnitTetraIntersectionBaryTest::test_cuboid_mapped_coords_3D()
397 double nodes[8][3] = { { 0.0, 2.0, 4.0 }, //0
401 { 0.0, 2.0, 0.0 }, // 4
407 for (int i=0; i < 8; ++i)
408 for (int j=0; j < 3; ++j)
411 std::vector<const double*> n (8);
412 for (int i=0; i<8; i++)
417 double p[3] = { 15.5, 16.0, 17.0 }, bc[3];
418 cuboid_mapped_coords(n, p, bc);
419 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[0], 1e-12);
420 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[1], 1e-12);
421 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[2], 1e-12);
425 double p[3] = { 15.0, 15.0, 19.0 }, bc[3];
426 cuboid_mapped_coords(n, p, bc);
427 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[0], 1e-12);
428 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[1], 1e-12);
429 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[2], 1e-12);
433 double p[3] = { 16.0, 17.0, 15.0 }, bc[3];
434 cuboid_mapped_coords(n, p, bc);
435 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
436 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
437 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[2], 1e-12);
441 double p[3] = { 16.0, 17.0, 19.0 }, bc[3];
442 cuboid_mapped_coords(n, p, bc);
443 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
444 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
445 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[2], 1e-12);
449 double p[3] = { 2.0, 16.0, 18.0 }, bc[3];
450 CPPUNIT_ASSERT_THROW(cuboid_mapped_coords(n, p, bc), INTERP_KERNEL::Exception);
456 - for QUAD4, point 0 is taken to be the origin (again see med file ref connec):
463 void UnitTetraIntersectionBaryTest::test_quad_mapped_coords_2D()
466 double nodes[4][2] = { { 0.0, 0.0 },
472 for (int i=0; i < 4; ++i)
473 for (int j=0; j < 2; ++j)
476 std::vector<const double*> n (4);
477 for (int i=0; i<4; i++)
482 double p[2] = { 15.75, 16.0 }, bc[2];
483 quad_mapped_coords(n, p, bc);
484 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[0], 1e-12);
485 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[1], 1e-12);
489 // middle point of seg
490 double p[2] = { 15.5, 15.0 }, bc[2];
491 quad_mapped_coords(n, p, bc);
492 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[0], 1e-12);
493 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[1], 1e-12);
498 double p[2] = { 15.0, 16.0 }, bc[2];
499 quad_mapped_coords(n, p, bc);
500 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[0], 1e-12);
501 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
505 double p[2] = { 17.0, 18.0 }, bc[2];
506 quad_mapped_coords(n, p, bc);
507 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
508 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
512 double p[2] = { 16.0, 15.0 }, bc[2];
513 quad_mapped_coords(n, p, bc);
514 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
515 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[1], 1e-12);
519 double p[2] = { 18.0, 18.0 }, bc[2];
520 CPPUNIT_ASSERT_THROW(quad_mapped_coords(n, p, bc), INTERP_KERNEL::Exception);