#define __INTERPOLATIONPLANAR_TXX__
#include "InterpolationPlanar.hxx"
+#include "OptionManager.hxx"
#include "PlanarIntersector.hxx"
#include "PlanarIntersector.txx"
#include "TriangulationIntersector.hxx"
*/
template<class RealPlanar>
InterpolationPlanar<RealPlanar>::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);
+
}
/**
* - Default: 0.
*/
template<class RealPlanar>
- void InterpolationPlanar<RealPlanar>::setOptions(double precision, int printLevel, IntersectionType intersectionType)
+ void InterpolationPlanar<RealPlanar>::setOptions(double precision, int printLevel, IntersectionType intersectionType, int orientation)
{
_precision=precision;
_printLevel=printLevel;
_intersectionType=intersectionType;
+ _orientation = orientation;
}
int nb_nodesS=connIndxS[i_S+1]-connIndxS[i_S];
double surf=intersector->intersectCells(OTT<ConnType,numPol>::indFC(i_P),
OTT<ConnType,numPol>::indFC(i_S),nb_nodesP,nb_nodesS);
- result[i_P].insert(std::make_pair(OTT<ConnType,numPol>::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<ConnType,numPol>::indFC(i_S),surf));
counter++;
- //rempli_vect_d_intersect(i_P+1,i_S+1,MaP,MaS,result_areas);
}
intersecting_elems.clear();
}
#include <string>
#include <vector>
#include <algorithm>
+#include <iostream>
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]);
{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::pair<double,double>theta1, std::pair<double,double> 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);
+
+}
+ };
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
{
double *COS=new double[taille/2];
double *SIN=new double[taille/2];
- double *angle=new double[taille/2];
- std::vector<double> Bary=bary_poly(V);
+ //double *angle=new double[taille/2];
+ std::vector<double> 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<taille/2-1;i++)
{
std::vector<double> 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<double> Pt_ordonne;
Pt_ordonne.reserve(taille);
- std::multimap<double,int> Ordre;
+ // std::multimap<double,int> Ordre;
+ std::multimap<std::pair<double,double>,int, AngleLess> CosSin;
for(int i=0;i<taille/2;i++)
- {Ordre.insert(std::make_pair(angle[i],i));}
- std::multimap <double,int>::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 <double,int>::iterator mi;
+ std::multimap<std::pair<double,double>,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;
}
}
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
/* Computes the dot product of a and b */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
- template<int dim> inline double dotprod( double * a, double * b)
+ template<int dim>
+ inline double dotprod( double * a, double * b)
{
double result=0;
for(int idim = 0; idim < dim ; idim++) result += a[idim]*b[idim];
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
/* Computes the norm of vector v */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
- template<int dim> inline double norm( double * v)
+ template<int dim>
+ inline double norm( double * v)
{
double result =0;
for(int idim =0; idim<dim; idim++) result+=v[idim]*v[idim];
return sqrt(result);
}
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
- /* Computes the norm of vector a-b */
+ /* Computes the square norm of vector a-b */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
- template<int dim> inline double distance2( const double * a, const double * b)
+ template<int dim>
+ inline double distance2( const double * a, const double * b)
{
double result =0;
for(int idim =0; idim<dim; idim++) result+=(a[idim]-b[idim])*(a[idim]-b[idim]);
return result;
+ }
+ template<class T, int dim>
+ inline double distance2( T * a, int inda, T * b, int indb)
+ {
+ double result =0;
+ for(int idim =0; idim<dim; idim++) result += ((*a)[inda+idim] - (*b)[indb+idim])* ((*a)[inda+idim] - (*b)[indb+idim]);
+ return result;
}
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
/* Computes the determinant of a and b */
}
return result;
}
+
+ /*! 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];
+ }
+
+ /*! Subroutine of checkEqualPolygins that tests if two list of nodes (not necessarily distincts) describe the same polygon, assuming they share a comon point.*/
+ /*! Indexes istart1 and istart2 designate two points P1 in L1 and P2 in L2 that have identical coordinates. Generally called with istart1=0.*/
+ /*! Integer sign ( 1 or -1) indicate the direction used in going all over L2. */
+ template<class T, int dim>
+ 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<T,dim>(L1,i1*dim, L1,i1next*dim) < epsilon ) i1next = ( i1next + 1 ) % size1;
+ while( i2next!=istart2 && distance2<T,dim>(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<T,dim>(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<class T, int dim>
+ 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<T,dim>(L1,istart1*dim, L2,istart2*dim) > epsilon ) istart2++;
+
+ if(istart2 == size2)
+ {
+ return (size1 == 0) && (size2 == 0);
+ }
+ else
+ return checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon, 1)
+ || checkEqualPolygonsOneDirection<T,dim>( 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
