Salome HOME
Indices are stored as mcIdType type instead of int to support switch to 64bits indexing
[tools/medcoupling.git] / src / INTERP_KERNELTest / ThreeDSurfProjectionTest.cxx
1 // Copyright (C) 2007-2019  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 #include "ThreeDSurfProjectionTest.hxx"
21 #include "PlanarIntersector.txx"
22 #include "MCIdType.hxx"
23
24 class MyMeshType
25 {
26 public:
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;
31 };
32
33 class MyMatrixType
34 {
35 };
36
37 void INTERP_TEST::ThreeDSurfProjectionTest::test1()
38 {
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)
42   //eps=1e-6
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]);
49   int ret;
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.};
57   for(int i=0;i<9;i++)
58     {
59       CPPUNIT_ASSERT_DOUBLES_EQUAL(expected0[i],tmp0[i],1e-15);
60       CPPUNIT_ASSERT_DOUBLES_EQUAL(expected0[i],tmp1[i],1e-15);
61     }
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.};
74   for(int i=0;i<9;i++)
75     {
76       CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],tmp0[i],1e-14);
77       CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],tmp1[i],1e-14);
78     }
79   //
80   delete [] tmp0;
81   delete [] tmp1;
82 }
83
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.]
90   //
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]);
94   int ret;
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);
126   //
127   delete [] tmp0;
128   delete [] tmp1;
129 }