From 0d7c1dbf0aa40e4d5d4f1164814cd3adce375a89 Mon Sep 17 00:00:00 2001 From: geay Date: Fri, 21 Feb 2014 16:14:13 +0100 Subject: [PATCH] On the road of 3D surf interpolation taking into account of non coincident 3D surf cells --- src/INTERP_KERNEL/PlanarIntersector.txx | 67 +++++------ src/INTERP_KERNELTest/CMakeLists.txt | 1 + src/INTERP_KERNELTest/TestInterpKernel.cxx | 3 +- .../ThreeDSurfProjectionTest.cxx | 106 ++++++++++++++++++ .../ThreeDSurfProjectionTest.hxx | 42 +++++++ 5 files changed, 186 insertions(+), 33 deletions(-) create mode 100644 src/INTERP_KERNELTest/ThreeDSurfProjectionTest.cxx create mode 100644 src/INTERP_KERNELTest/ThreeDSurfProjectionTest.hxx diff --git a/src/INTERP_KERNEL/PlanarIntersector.txx b/src/INTERP_KERNEL/PlanarIntersector.txx index ffc88f196..8c1df161a 100644 --- a/src/INTERP_KERNEL/PlanarIntersector.txx +++ b/src/INTERP_KERNEL/PlanarIntersector.txx @@ -289,11 +289,12 @@ namespace INTERP_KERNEL bool same_orientation; //Find the normal to cells A and B - int i_A1=1; - while(i_A1(Coords_A,&Coords_A[SPACEDIM*i_A1])< epsilon) i_A1++; - int i_A2=i_A1+1; + int i_A1(1); + while(i_A1(Coords_A,&Coords_A[SPACEDIM*i_A1])< epsilon) + i_A1++; + int i_A2(i_A1+1); crossprod(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A); - double normA = sqrt(dotprod(normal_A,normal_A)); + double normA(sqrt(dotprod(normal_A,normal_A))); while(i_A2(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A); @@ -301,11 +302,12 @@ namespace INTERP_KERNEL normA = sqrt(dotprod(normal_A,normal_A)); } - int i_B1=1; - while(i_B1(Coords_B,Coords_B+SPACEDIM*i_B1)< epsilon) i_B1++; - int i_B2=i_B1+1; + int i_B1(1); + while(i_B1(Coords_B,Coords_B+SPACEDIM*i_B1)< epsilon) + i_B1++; + int i_B2(i_B1+1); crossprod(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B); - double normB = sqrt(dotprod(normal_B,normal_B)); + double normB(sqrt(dotprod(normal_B,normal_B))); while(i_B2(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B); @@ -317,7 +319,7 @@ namespace INTERP_KERNEL if(md3DSurf>0.) { double coords_GA[3]; - for (int i=0;i<3;i++) + for(int i=0;i<3;i++) { coords_GA[i]=0.; for (int j=0;jmd3DSurf) - return 0; + for(int i=0;i<3;i++) + { + G1[i]=Coords_B[i]-coords_GA[i]; + G2[i]=Coords_B[i+3]-coords_GA[i]; + G3[i]=Coords_B[i+6]-coords_GA[i]; + } + double prodvect[3]; + prodvect[0]=G1[1]*G2[2]-G1[2]*G2[1]; + prodvect[1]=G1[2]*G2[0]-G1[0]*G2[2]; + prodvect[2]=G1[0]*G2[1]-G1[1]*G2[0]; + double prodscal=prodvect[0]*G3[0]+prodvect[1]*G3[1]+prodvect[2]*G3[2]; + if(fabs(prodscal)>md3DSurf) + return 0; } if(i_A2(normal_A,normal_B)>=0; + same_orientation=dotprod(normal_A,normal_B)>=0; if(!same_orientation) - for(int idim =0; idim< SPACEDIM; idim++) normal_A[idim] *=-1; + for(int idim =0; idim< SPACEDIM; idim++) + normal_A[idim] *=-1; - double normBB= sqrt(dotprod(normal_B,normal_B)); + double normBB(sqrt(dotprod(normal_B,normal_B))); for(int idim =0; idim< SPACEDIM; idim++) linear_comb[idim] = median_plane*normal_A[idim]/normA + (1-median_plane)*normal_B[idim]/normBB; double norm= sqrt(dotprod(linear_comb,linear_comb)); //Necessarily: norm>epsilon, no need to check - for(int idim =0; idim< SPACEDIM; idim++) linear_comb[idim]/=norm; + for(int idim =0; idim< SPACEDIM; idim++) + linear_comb[idim]/=norm; //Project the nodes of A and B on the median plane for(int i_A=0; i_A(&Coords_B[0],&Coords_B[i_B1])= " << distance2(Coords_B,Coords_B+i_B1) << std::endl; std::cout << "normal_B = " << normal_B[0] << " ; " << normal_B[1] << " ; " << normal_B[2] << std::endl; - return 1; } } diff --git a/src/INTERP_KERNELTest/CMakeLists.txt b/src/INTERP_KERNELTest/CMakeLists.txt index f011e163e..2a6855850 100644 --- a/src/INTERP_KERNELTest/CMakeLists.txt +++ b/src/INTERP_KERNELTest/CMakeLists.txt @@ -49,6 +49,7 @@ SET(InterpKernelTest_SOURCES UnitTetra3D2DIntersectionTest.cxx UnitTetraIntersectionBaryTest.cxx TestInterpKernelUtils.cxx + ThreeDSurfProjectionTest.cxx ) SET(TestINTERP_KERNEL_SOURCES diff --git a/src/INTERP_KERNELTest/TestInterpKernel.cxx b/src/INTERP_KERNELTest/TestInterpKernel.cxx index f157962cc..618b3b8c1 100644 --- a/src/INTERP_KERNELTest/TestInterpKernel.cxx +++ b/src/INTERP_KERNELTest/TestInterpKernel.cxx @@ -32,6 +32,7 @@ #include "MultiElement2DTests.hxx" #include "MultiElementTetraTests.hxx" #include "SingleElementTetraTests.hxx" +#include "ThreeDSurfProjectionTest.hxx" using namespace INTERP_TEST; @@ -51,7 +52,7 @@ CPPUNIT_TEST_SUITE_REGISTRATION( HexaTests ); CPPUNIT_TEST_SUITE_REGISTRATION( MultiElement2DTests ); CPPUNIT_TEST_SUITE_REGISTRATION( MultiElementTetraTests ); CPPUNIT_TEST_SUITE_REGISTRATION( SingleElementTetraTests ); - +CPPUNIT_TEST_SUITE_REGISTRATION( ThreeDSurfProjectionTest ); // --- generic Main program from KERNEL_SRC/src/Basics/Test #include "BasicMainTest.hxx" diff --git a/src/INTERP_KERNELTest/ThreeDSurfProjectionTest.cxx b/src/INTERP_KERNELTest/ThreeDSurfProjectionTest.cxx new file mode 100644 index 000000000..f052edbb4 --- /dev/null +++ b/src/INTERP_KERNELTest/ThreeDSurfProjectionTest.cxx @@ -0,0 +1,106 @@ +// Copyright (C) 2007-2014 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#include "ThreeDSurfProjectionTest.hxx" +#include "PlanarIntersector.txx" + +class MyMeshType +{ +public: + static const int MY_SPACEDIM=3; + typedef int MyConnType; +}; + +class MyMatrixType +{ +}; + +void INTERP_TEST::ThreeDSurfProjectionTest::test1() +{ + // Two triangles coo and coo2 are perfectly // each others with a distance equal to 1e-6. + // A little rotation to make it more funny. + //coo=DataArrayDouble([0.,0.,0.,1.,0.,0.,0.,1.,0.],3,3) + //eps=1e-6 + //coo2=DataArrayDouble([0.,0.,eps,1.,0.,eps,0.,1.,eps],3,3) + //MEDCouplingPointSet.Rotate3DAlg([0.,0.,0.],[2.,1.,3.],0.3,coo) + //MEDCouplingPointSet.Rotate3DAlg([0.,0.,0.],[2.,1.,3.],0.3,coo2) + const double coo[9]={0.,0.,0.,0.96809749223257568,0.24332379388106262,-0.059839592782071335,-0.23056279077409292,0.95852673990234838,0.16753294721527912}; + 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}; + double *tmp0(new double[9]),*tmp1(new double[9]); + int ret; + //eps=1e-2. eps is a tolerance to detect that two points are the same or not in a same polygon. + // here the max 3D distance is 1e-5 > 1e-6 so 1 is expected + std::copy(coo,coo+9,tmp0); + std::copy(coo2,coo2+9,tmp1); + ret=INTERP_KERNEL::PlanarIntersector::projection(tmp0,tmp1,3,3,1e-2,1e-5/* <- */,0.5,true); + CPPUNIT_ASSERT_EQUAL(1,ret); + const double expected0[9]={0.,0.,0.,1.,0.,0.,0.,1.,0.}; + for(int i=0;i<9;i++) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected0[i],tmp0[i],1e-15); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected0[i],tmp1[i],1e-15); + } + // here the max 3D distance is 1e-8 < 1e-6 so 0 is expected + std::copy(coo,coo+9,tmp0); + std::copy(coo2,coo2+9,tmp1); + ret=INTERP_KERNEL::PlanarIntersector::projection(tmp0,tmp1,3,3,1e-2,1e-8/* <- */,0.5,true); + CPPUNIT_ASSERT_EQUAL(0,ret); + // here testing when max 3D distance is 1e-5 > 1e-6 with inverted cells + std::copy(coo,coo+9,tmp0); + std::copy(coo2,coo2+3,tmp1+6); std::copy(coo2+3,coo2+6,tmp1+3); std::copy(coo2+6,coo2+9,tmp1); + ret=INTERP_KERNEL::PlanarIntersector::projection(tmp0,tmp1,3,3,1e-2,1e-5/* <- */,0.5,true); + CPPUNIT_ASSERT_EQUAL(-1,ret); + const double expected1[9]={-0.7071067811865476,-0.7071067811865476,0.,0.,-1.4142135623730951,0.,-1.4142135623730951,-1.4142135623730951,0.}; + const double expected2[9]={-1.4142135623730951,-1.4142135623730951,0.,0.,-1.4142135623730951,0.,-0.7071067811865476,-0.7071067811865476,0.}; + for(int i=0;i<9;i++) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],tmp0[i],1e-14); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],tmp1[i],1e-14); + } + // + delete [] tmp0; + delete [] tmp1; +} + +void INTERP_TEST::ThreeDSurfProjectionTest::test2() +{// here the two triangles have their center of inertia very close (eps) but the angle between the two planes is "big" + //coo=DataArrayDouble([0.,0.,0.,1.,0.,0.,0.,1.,0.],3,3) + //coocpy=coo.deepCpy() + //MEDCouplingPointSet.Rotate3DAlg([0.,0.,0.],[-1,-1.,0.],pi/3,coocpy) + //coocpy+=[eps*sqrt(3)/2,eps/2,eps*0.] + // + const double coo[9]={0.,0.,0.,0.96809749223257568,0.24332379388106262,-0.059839592782071335,-0.23056279077409292,0.95852673990234838,0.16753294721527912}; + 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}; + double *tmp0(new double[9]),*tmp1(new double[9]); + int ret; + //eps=1e-2. eps is a tolerance to detect that two points are the same or not in a same polygon. + // here the max 3D distance is 1e-5 > 1e-6 so 1 is expected + std::copy(coo,coo+9,tmp0); + std::copy(coo2,coo2+9,tmp1); + ret=INTERP_KERNEL::PlanarIntersector::projection(tmp0,tmp1,3,3,1e-2,1e-5/* <- */,0.5,true); + CPPUNIT_ASSERT_EQUAL(1,ret); + // here the max 3D distance is 1e-8 < 1e-6 so 0 is expected + std::copy(coo,coo+9,tmp0); + std::copy(coo2,coo2+9,tmp1); + ret=INTERP_KERNEL::PlanarIntersector::projection(tmp0,tmp1,3,3,1e-2,1e-8/* <- */,0.5,true); + CPPUNIT_ASSERT_EQUAL(0,ret); + // + delete [] tmp0; + delete [] tmp1; +} diff --git a/src/INTERP_KERNELTest/ThreeDSurfProjectionTest.hxx b/src/INTERP_KERNELTest/ThreeDSurfProjectionTest.hxx new file mode 100644 index 000000000..f344cd92d --- /dev/null +++ b/src/INTERP_KERNELTest/ThreeDSurfProjectionTest.hxx @@ -0,0 +1,42 @@ +// Copyright (C) 2007-2014 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#ifndef __THREEDSURFPROJECTIONTEST_HXX__ +#define __THREEDSURFPROJECTIONTEST_HXX__ + +#include + +namespace INTERP_TEST +{ + /** + * \brief Class dedicated of the test of the preprocessing of 3D surf cells before performing invoking 2D algorithms. + */ + class ThreeDSurfProjectionTest : public CppUnit::TestFixture + { + CPPUNIT_TEST_SUITE( ThreeDSurfProjectionTest ); + CPPUNIT_TEST ( test1 ); + CPPUNIT_TEST ( test2 ); + CPPUNIT_TEST_SUITE_END(); + public: + void test1(); + void test2(); + }; +} + +#endif -- 2.39.2