1 // Copyright (C) 2007-2022 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 #include "ThreeDSurfProjectionTest.hxx"
21 #include "PlanarIntersector.txx"
22 #include "MCIdType.hxx"
27 static const int MY_SPACEDIM=3;
28 static const int MY_MESHDIM=3;
29 static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
30 typedef mcIdType MyConnType;
37 void INTERP_TEST::ThreeDSurfProjectionTest::test1()
39 // Two triangles coo and coo2 are perfectly // each others with a distance equal to 1e-6.
40 // A little rotation to make it more funny.
41 //coo=DataArrayDouble([0.,0.,0.,1.,0.,0.,0.,1.,0.],3,3)
43 //coo2=DataArrayDouble([0.,0.,eps,1.,0.,eps,0.,1.,eps],3,3)
44 //MEDCouplingPointSet.Rotate3DAlg([0.,0.,0.],[2.,1.,3.],0.3,coo)
45 //MEDCouplingPointSet.Rotate3DAlg([0.,0.,0.],[2.,1.,3.],0.3,coo2)
46 const double coo[9]={0.,0.,0.,0.96809749223257568,0.24332379388106262,-0.059839592782071335,-0.23056279077409292,0.95852673990234838,0.16753294721527912};
47 const double coo2[9]={9.8122602102980502e-08,-1.4839144255482456e-7,9.8404874611628791e-7,0.96809759035517784,0.24332364548962007,-0.059838608733325221,-0.23056269265149082,0.9585265915109058,0.16753393126402524};
48 double *tmp0(new double[9]),*tmp1(new double[9]);
50 //eps=1e-2. eps is a tolerance to detect that two points are the same or not in a same polygon.
51 // here the max 3D distance is 1e-5 > 1e-6 so 1 is expected
52 std::copy(coo,coo+9,tmp0);
53 std::copy(coo2,coo2+9,tmp1);
54 ret=INTERP_KERNEL::PlanarIntersector<MyMeshType,MyMatrixType>::Projection(tmp0,tmp1,3,3,1e-2,1e-5/* <- */,-1.,0.5,true);
55 CPPUNIT_ASSERT_EQUAL(1,ret);
56 const double expected0[9]={0.,0.,0.,1.,0.,0.,0.,1.,0.};
59 CPPUNIT_ASSERT_DOUBLES_EQUAL(expected0[i],tmp0[i],1e-15);
60 CPPUNIT_ASSERT_DOUBLES_EQUAL(expected0[i],tmp1[i],1e-15);
62 // here the max 3D distance is 1e-8 < 1e-6 so 0 is expected
63 std::copy(coo,coo+9,tmp0);
64 std::copy(coo2,coo2+9,tmp1);
65 ret=INTERP_KERNEL::PlanarIntersector<MyMeshType,MyMatrixType>::Projection(tmp0,tmp1,3,3,1e-2,1e-8/* <- */,-1.,0.5,true);
66 CPPUNIT_ASSERT_EQUAL(0,ret);
67 // here testing when max 3D distance is 1e-5 > 1e-6 with inverted cells
68 std::copy(coo,coo+9,tmp0);
69 std::copy(coo2,coo2+3,tmp1+6); std::copy(coo2+3,coo2+6,tmp1+3); std::copy(coo2+6,coo2+9,tmp1);
70 ret=INTERP_KERNEL::PlanarIntersector<MyMeshType,MyMatrixType>::Projection(tmp0,tmp1,3,3,1e-2,1e-5/* <- */,-1.,0.5,true);
71 CPPUNIT_ASSERT_EQUAL(-1,ret);
72 const double expected1[9]={-0.7071067811865476,-0.7071067811865476,0.,0.,-1.4142135623730951,0.,-1.4142135623730951,-1.4142135623730951,0.};
73 const double expected2[9]={-1.4142135623730951,-1.4142135623730951,0.,0.,-1.4142135623730951,0.,-0.7071067811865476,-0.7071067811865476,0.};
76 CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],tmp0[i],1e-14);
77 CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],tmp1[i],1e-14);
84 void INTERP_TEST::ThreeDSurfProjectionTest::test2()
85 {// here the two triangles have their center of inertia very close (eps) but the angle between the two planes is "big"
86 //coo=DataArrayDouble([0.,0.,0.,1.,0.,0.,0.,1.,0.],3,3)
87 //coocpy=coo.deepCopy()
88 //MEDCouplingPointSet.Rotate3DAlg([0.,0.,0.],[-1,-1.,0.],pi/3,coocpy)
89 //coocpy+=[eps*sqrt(3)/2,eps/2,eps*0.]
91 const double coo[9]={0.,0.,0.,0.96809749223257568,0.24332379388106262,-0.059839592782071335,-0.23056279077409292,0.95852673990234838,0.16753294721527912};
92 const double coo2[9]={7.2311562622637225e-07,6.8998795679738294e-07,3.1943866106249849e-08,0.72852072144314628,0.33125439126063028,0.5996079016637561,0.0090154262465889021,0.87059752249869415,-0.49191448334281612};
93 double *tmp0(new double[9]),*tmp1(new double[9]);
95 //eps=1e-2. eps is a tolerance to detect that two points are the same or not in a same polygon.
96 // here the max 3D distance is 1e-5 > 1e-6 so 1 is expected
97 std::copy(coo,coo+9,tmp0);
98 std::copy(coo2,coo2+9,tmp1);
99 ret=INTERP_KERNEL::PlanarIntersector<MyMeshType,MyMatrixType>::Projection(tmp0,tmp1,3,3,1e-2,1e-5/* <- */,-1.,0.5,true);
100 CPPUNIT_ASSERT_EQUAL(1,ret);
101 // here the max 3D distance is 1e-8 < 1e-6 so 0 is expected
102 std::copy(coo,coo+9,tmp0);
103 std::copy(coo2,coo2+9,tmp1);
104 ret=INTERP_KERNEL::PlanarIntersector<MyMeshType,MyMatrixType>::Projection(tmp0,tmp1,3,3,1e-2,1e-8/* <- */,-1.,0.5,true);
105 CPPUNIT_ASSERT_EQUAL(0,ret);
106 // again max 3D distance is 1e-5 > 1e-6 so 1 is expected
107 std::copy(coo,coo+9,tmp0);
108 std::copy(coo2,coo2+9,tmp1);
109 ret=INTERP_KERNEL::PlanarIntersector<MyMeshType,MyMatrixType>::Projection(tmp0,tmp1,3,3,1e-2,1e-5/* <- */,-1.,0.5,true);
110 CPPUNIT_ASSERT_EQUAL(1,ret);
111 // again max 3D distance is 1e-5 > 1e-6 but minDot set to 0.8. 0 expected. because the angle is pi/4 so cos(pi/3) > 0.8
112 std::copy(coo,coo+9,tmp0);
113 std::copy(coo2,coo2+9,tmp1);
114 ret=INTERP_KERNEL::PlanarIntersector<MyMeshType,MyMatrixType>::Projection(tmp0,tmp1,3,3,1e-2,1e-5/* <- */,0.8/* <- */,0.5,true);
115 CPPUNIT_ASSERT_EQUAL(0,ret);
116 // again max 3D distance is 1e-5 > 1e-6 but minDot set to 0.7. 1 expected. because the angle is pi/4 so cos(pi/3) < 0.49
117 std::copy(coo,coo+9,tmp0);
118 std::copy(coo2,coo2+9,tmp1);
119 ret=INTERP_KERNEL::PlanarIntersector<MyMeshType,MyMatrixType>::Projection(tmp0,tmp1,3,3,1e-2,1e-5/* <- */,0.49/* <- */,0.5,true);
120 CPPUNIT_ASSERT_EQUAL(1,ret);
121 // again max 3D distance is 1e-5 > 1e-6 but minDot set to 0.7. 0 expected. because the angle is pi/4 so cos(pi/3) > 0.51
122 std::copy(coo,coo+9,tmp0);
123 std::copy(coo2,coo2+9,tmp1);
124 ret=INTERP_KERNEL::PlanarIntersector<MyMeshType,MyMatrixType>::Projection(tmp0,tmp1,3,3,1e-2,1e-5/* <- */,0.51/* <- */,0.5,true);
125 CPPUNIT_ASSERT_EQUAL(0,ret);