From 3926b9ab64c157a152953e28bf269a832ffc5def Mon Sep 17 00:00:00 2001 From: ndjinga Date: Fri, 26 Sep 2008 14:36:44 +0000 Subject: [PATCH] Added orientation parameter for filtering --- src/INTERP_KERNEL/ConvexIntersector.txx | 9 +- src/INTERP_KERNEL/Interpolation3DSurf.hxx | 2 +- src/INTERP_KERNEL/Interpolation3DSurf.txx | 4 +- src/INTERP_KERNEL/InterpolationPlanar.hxx | 4 +- src/INTERP_KERNEL/InterpolationPlanar.txx | 18 +- src/INTERP_KERNEL/InterpolationUtils.hxx | 193 ++++++++++++++---- .../TriangulationIntersector.txx | 11 +- 7 files changed, 183 insertions(+), 58 deletions(-) diff --git a/src/INTERP_KERNEL/ConvexIntersector.txx b/src/INTERP_KERNEL/ConvexIntersector.txx index 3252e5c06..53af77e48 100644 --- a/src/INTERP_KERNEL/ConvexIntersector.txx +++ b/src/INTERP_KERNEL/ConvexIntersector.txx @@ -36,7 +36,9 @@ namespace INTERP_KERNEL double ConvexIntersector::intersectCells(ConnType icell_A, ConnType icell_B, int nb_NodesA, int nb_NodesB) { - double result=0; + double result = 0; + int orientation = 1; + /*** Obtain the coordinates of A and B ***/ std::vector Coords_A (SPACEDIM*nb_NodesA); std::vector Coords_B (SPACEDIM*nb_NodesB); @@ -74,7 +76,8 @@ namespace INTERP_KERNEL /*** project cells A and B on the median plane ***/ /*** and rotate the median plane ***/ - if(SPACEDIM==3) projection(Coords_A, Coords_B, nb_dist_NodesA, nb_dist_NodesB,_epsilon, + if(SPACEDIM==3) + orientation = projection(Coords_A, Coords_B, nb_dist_NodesA, nb_dist_NodesB,_epsilon, PlanarIntersector::_medianPlane, PlanarIntersector::_doRotate); @@ -111,7 +114,7 @@ namespace INTERP_KERNEL std::cout << std::endl <<"Intersection area = " << result << std::endl; } - return result; + return orientation*result; } } diff --git a/src/INTERP_KERNEL/Interpolation3DSurf.hxx b/src/INTERP_KERNEL/Interpolation3DSurf.hxx index f6c884c72..7aeb2f69a 100644 --- a/src/INTERP_KERNEL/Interpolation3DSurf.hxx +++ b/src/INTERP_KERNEL/Interpolation3DSurf.hxx @@ -10,7 +10,7 @@ namespace INTERP_KERNEL public: Interpolation3DSurf(); void setOptions(double precision, int printLevel, double medianPlane, - IntersectionType intersectionType, bool doRotate); + IntersectionType intersectionType, bool doRotate, int orientation); public: bool doRotate() const { return _doRotate; } double medianPlane() const { return _medianPlane; } diff --git a/src/INTERP_KERNEL/Interpolation3DSurf.txx b/src/INTERP_KERNEL/Interpolation3DSurf.txx index b9e14f57f..cd0b7189b 100644 --- a/src/INTERP_KERNEL/Interpolation3DSurf.txx +++ b/src/INTERP_KERNEL/Interpolation3DSurf.txx @@ -38,9 +38,9 @@ namespace INTERP_KERNEL - Default: 0. */ void Interpolation3DSurf::setOptions(double precision, int printLevel, double medianPlane, - IntersectionType intersectionType, bool doRotate) + IntersectionType intersectionType, bool doRotate, int orientation) { - InterpolationPlanar::setOptions(precision,printLevel,intersectionType); + InterpolationPlanar::setOptions(precision,printLevel,intersectionType, orientation); _doRotate=doRotate; _medianPlane=medianPlane; } diff --git a/src/INTERP_KERNEL/InterpolationPlanar.hxx b/src/INTERP_KERNEL/InterpolationPlanar.hxx index 22307e8b3..94223f375 100755 --- a/src/INTERP_KERNEL/InterpolationPlanar.hxx +++ b/src/INTERP_KERNEL/InterpolationPlanar.hxx @@ -12,6 +12,7 @@ namespace INTERP_KERNEL template class InterpolationPlanar : public Interpolation< InterpolationPlanar >, public OptionManager + { private: int _printLevel; @@ -19,12 +20,13 @@ namespace INTERP_KERNEL double _dimCaracteristic; IntersectionType _intersectionType; static const double DEFAULT_PRECISION; + int _orientation; public: InterpolationPlanar(); // geometric precision, debug print level, coice of the median plane, intersection etc ... void setOptions(double precision, int printLevel, - IntersectionType intersectionType); + IntersectionType intersectionType, int orientation); // Main function to interpolate triangular and quadratic meshes template diff --git a/src/INTERP_KERNEL/InterpolationPlanar.txx b/src/INTERP_KERNEL/InterpolationPlanar.txx index 3aa78b5c3..4f3464794 100644 --- a/src/INTERP_KERNEL/InterpolationPlanar.txx +++ b/src/INTERP_KERNEL/InterpolationPlanar.txx @@ -2,6 +2,7 @@ #define __INTERPOLATIONPLANAR_TXX__ #include "InterpolationPlanar.hxx" +#include "OptionManager.hxx" #include "PlanarIntersector.hxx" #include "PlanarIntersector.txx" #include "TriangulationIntersector.hxx" @@ -29,11 +30,13 @@ namespace INTERP_KERNEL */ template InterpolationPlanar::InterpolationPlanar():_printLevel(0),_precision(DEFAULT_PRECISION),_dimCaracteristic(1), - _intersectionType(Triangulation) + _intersectionType(Triangulation),_orientation(1) { registerOption(&_precision, "Precision",DEFAULT_PRECISION); registerOption(&_intersectionType, "IntersectionType",Triangulation); registerOption(&_printLevel, "PrintLevel",0); + registerOption(&_orientation, "IntersectionSign",1); + } /** @@ -50,11 +53,12 @@ namespace INTERP_KERNEL * - Default: 0. */ template - void InterpolationPlanar::setOptions(double precision, int printLevel, IntersectionType intersectionType) + void InterpolationPlanar::setOptions(double precision, int printLevel, IntersectionType intersectionType, int orientation) { _precision=precision; _printLevel=printLevel; _intersectionType=intersectionType; + _orientation = orientation; } @@ -177,9 +181,15 @@ namespace INTERP_KERNEL int nb_nodesS=connIndxS[i_S+1]-connIndxS[i_S]; double surf=intersector->intersectCells(OTT::indFC(i_P), OTT::indFC(i_S),nb_nodesP,nb_nodesS); - result[i_P].insert(std::make_pair(OTT::indFC(i_S),surf)); + + //filtering out zero surfaces and badly oriented surfaces + // orientation = -1,0,1 + // -1 : the intersection is taken into account if target and cells have different orientation + // 0 : the intersection is always taken into account + // 1 : the intersection is taken into account if target and cells have the same orientation + if (( surf > 0.0 && (double)_orientation >=0 ) || ( surf < 0.0 && (double)_orientation <=0 )) + result[i_P].insert(std::make_pair(OTT::indFC(i_S),surf)); counter++; - //rempli_vect_d_intersect(i_P+1,i_S+1,MaP,MaS,result_areas); } intersecting_elems.clear(); } diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index d9544d1e1..71b59d16c 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -9,6 +9,7 @@ #include #include #include +#include namespace INTERP_KERNEL { @@ -252,12 +253,12 @@ namespace INTERP_KERNEL double k_1=-((P_3[1]-P_4[1])*(P_3[0]-P_1[0])+(P_4[0]-P_3[0])*(P_3[1]-P_1[1]))/det; - if( k_1>=0 && k_1<=1) + if( k_1>=-precision && k_1<=1+precision) { double k_2= ((P_1[1]-P_2[1])*(P_1[0]-P_3[0])+(P_2[0]-P_1[0])*(P_1[1]-P_3[1]))/det; - if( k_2>=0 && k_2<=1) + if( k_2>=-precision && k_2<=1+precision) { double P_0[2]; P_0[0]=P_1[0]+k_1*(P_2[0]-P_1[0]); @@ -315,8 +316,25 @@ namespace INTERP_KERNEL {V.push_back(Num); } } - - + /*! Function that compares two angles from the values of the pairs (sin,cos)*/ + /*! Angles are considered in [0, 2Pi] bt are not computed explicitely */ + class AngleLess + { + public: + bool operator()(std::pairtheta1, std::pair theta2) + { + double norm1 = sqrt(theta1.first*theta1.first +theta1.second*theta1.second); + double norm2 = sqrt(theta2.first*theta2.first +theta2.second*theta2.second); + + double epsilon = 1.e-12; + + if( norm1 < epsilon || norm2 < epsilon ) + std::cout << "Warning InterpolationUtils.hxx: AngleLess : Vector with zero norm, cannot define the angle !!!! " << std::endl; + + return theta1.second*(norm2 + theta2.first) < theta2.second*(norm1 + theta1.first); + +} + }; /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ @@ -337,38 +355,49 @@ namespace INTERP_KERNEL { double *COS=new double[taille/2]; double *SIN=new double[taille/2]; - double *angle=new double[taille/2]; - std::vector Bary=bary_poly(V); + //double *angle=new double[taille/2]; + std::vector Bary=bary_poly(V); COS[0]=1.0; SIN[0]=0.0; - angle[0]=0.0; + //angle[0]=0.0; for(int i=0; i Trigo=calcul_cos_et_sin(&Bary[0],&V[0],&V[2*(i+1)]); COS[i+1]=Trigo[0]; SIN[i+1]=Trigo[1]; - if(SIN[i+1]>=0) - {angle[i+1]=acos(COS[i+1]);} - else - {angle[i+1]=-acos(COS[i+1]);} + //if(SIN[i+1]>=0) + // {angle[i+1]=atan2(SIN[i+1],COS[i+1]);} +// else +// {angle[i+1]=-atan2(SIN[i+1],COS[i+1]);} } //ensuite on ordonne les angles. std::vector Pt_ordonne; Pt_ordonne.reserve(taille); - std::multimap Ordre; + // std::multimap Ordre; + std::multimap,int, AngleLess> CosSin; for(int i=0;i::iterator mi; - for(mi=Ordre.begin();mi!=Ordre.end();mi++) { - int j=(*mi).second; + // Ordre.insert(std::make_pair(angle[i],i)); + CosSin.insert(std::make_pair(std::make_pair(SIN[i],COS[i]),i)); + } + // std::multimap ::iterator mi; + std::multimap,int, AngleLess>::iterator micossin; +// for(mi=Ordre.begin();mi!=Ordre.end();mi++) +// { +// int j=(*mi).second; +// Pt_ordonne.push_back(V[2*j]); +// Pt_ordonne.push_back(V[2*j+1]); +// } + for(micossin=CosSin.begin();micossin!=CosSin.end();micossin++) + { + int j=(*micossin).second; Pt_ordonne.push_back(V[2*j]); Pt_ordonne.push_back(V[2*j+1]); } delete [] COS; delete [] SIN; - delete [] angle; + // delete [] angle; return Pt_ordonne; } } @@ -400,7 +429,8 @@ namespace INTERP_KERNEL /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ /* Computes the dot product of a and b */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ - template inline double dotprod( double * a, double * b) + template + inline double dotprod( double * a, double * b) { double result=0; for(int idim = 0; idim < dim ; idim++) result += a[idim]*b[idim]; @@ -409,20 +439,29 @@ namespace INTERP_KERNEL /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ /* Computes the norm of vector v */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ - template inline double norm( double * v) + template + inline double norm( double * v) { double result =0; for(int idim =0; idim inline double distance2( const double * a, const double * b) + template + inline double distance2( const double * a, const double * b) { double result =0; for(int idim =0; idim + inline double distance2( T * a, int inda, T * b, int indb) + { + double result =0; + for(int idim =0; idim + bool checkEqualPolygonsOneDirection(T * L1, T * L2, int size1, int size2, int istart1, int istart2, double epsilon, int sign) + { + int i1 = istart1; + int i2 = istart2; + int i1next = ( i1 + 1 ) % size1; + int i2next = ( i2 + sign +size2) % size2; + + while(true) + { + while( i1next!=istart1 && distance2(L1,i1*dim, L1,i1next*dim) < epsilon ) i1next = ( i1next + 1 ) % size1; + while( i2next!=istart2 && distance2(L2,i2*dim, L2,i2next*dim) < epsilon ) i2next = ( i2next + sign +size2 ) % size2; + + if(i1next == istart1) + { + if(i2next == istart2) + return true; + else return false; + } + else + if(i2next == istart2) + return false; + else + { + if(distance2(L1,i1next*dim, L2,i2next*dim) > epsilon ) + return false; + else + { + i1 = i1next; + i2 = i2next; + i1next = ( i1 + 1 ) % size1; + i2next = ( i2 + sign + size2 ) % size2; + } + } + } + } + + /*! Tests if two list of nodes (not necessarily distincts) describe the same polygon.*/ + /*! Existence of multiple points in the list is considered.*/ + template + bool checkEqualPolygons(T * L1, T * L2, double epsilon) + { + if(L1==NULL || L2==NULL) + { + std::cout << "Warning InterpolationUtils.hxx:checkEqualPolygonsPointer: Null pointer " << std::endl; + throw(Exception("big error: not closed polygon...")); + } + + int size1 = (*L1).size()/dim; + int size2 = (*L2).size()/dim; + int istart1 = 0; + int istart2 = 0; + + while( istart2 < size2 && distance2(L1,istart1*dim, L2,istart2*dim) > epsilon ) istart2++; + + if(istart2 == size2) + { + return (size1 == 0) && (size2 == 0); + } + else + return checkEqualPolygonsOneDirection( L1, L2, size1, size2, istart1, istart2, epsilon, 1) + || checkEqualPolygonsOneDirection( L1, L2, size1, size2, istart1, istart2, epsilon, -1); + + } } -/*! Computes the triple product (XA^XB).XC (in 3D)*/ -inline double triple_product(const double* A, const double*B, const double*C, const double*X) -{ - double XA[3]; - XA[0]=A[0]-X[0]; - XA[1]=A[1]-X[1]; - XA[2]=A[2]-X[2]; - double XB[3]; - XB[0]=B[0]-X[0]; - XB[1]=B[1]-X[1]; - XB[2]=B[2]-X[2]; - double XC[3]; - XC[0]=C[0]-X[0]; - XC[1]=C[1]-X[1]; - XC[2]=C[2]-X[2]; - - return - (XA[1]*XB[2]-XA[2]*XB[1])*XC[0]+ - (XA[2]*XB[0]-XA[0]*XB[2])*XC[1]+ - (XA[0]*XB[1]-XA[1]*XB[0])*XC[2]; -} + #endif diff --git a/src/INTERP_KERNEL/TriangulationIntersector.txx b/src/INTERP_KERNEL/TriangulationIntersector.txx index fc4dd7eab..ac4114189 100644 --- a/src/INTERP_KERNEL/TriangulationIntersector.txx +++ b/src/INTERP_KERNEL/TriangulationIntersector.txx @@ -33,7 +33,8 @@ namespace INTERP_KERNEL double TriangulationIntersector::intersectCells(ConnType icell_A, ConnType icell_B, int nb_NodesA, int nb_NodesB) { - double result=0.; + double result = 0.; + int orientation = 1; //Obtain the coordinates of A and B std::vector Coords_A (SPACEDIM*nb_NodesA); @@ -48,9 +49,9 @@ namespace INTERP_KERNEL //project cells A and B on the median plane and rotate the median plane if(SPACEDIM==3) - projection(Coords_A, Coords_B, nb_NodesA, nb_NodesB, - PlanarIntersector::_dimCaracteristic * PlanarIntersector::_precision, - PlanarIntersector::_medianPlane, PlanarIntersector::_doRotate); + orientation = projection(Coords_A, Coords_B, nb_NodesA, nb_NodesB, + PlanarIntersector::_dimCaracteristic * PlanarIntersector::_precision, + PlanarIntersector::_medianPlane, PlanarIntersector::_doRotate); //DEBUG PRINTS if(PlanarIntersector::_printLevel >= 3) @@ -96,7 +97,7 @@ namespace INTERP_KERNEL if(PlanarIntersector::_printLevel >= 3) std::cout << std::endl <<"Intersection area = " << result << std::endl; - return result; + return orientation*result; } } -- 2.39.2