#include <vector>
#include <algorithm>
-const int MIN_NB_ELEMS =10;
-const int MAX_LEVEL=10;
+const int MIN_NB_ELEMS =15;
+const int MAX_LEVEL=20;
using namespace std;
+
template <int dim> class BBTree
{
private:
BBTree* _left;
BBTree* _right;
- BBTree* _center;
int _level;
- double _median;
+ double _max_left;
+ double _min_right;
double* _bb;
vector <int> _elems;
bool _terminal;
+ int _nbelems;
public:
BBTree(double* bbs, int* elems, int level, int nbelems):
- _left(0), _right(0), _center(0), _level(level), _median(0.0), _bb(bbs), _terminal(false)
+ _left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems)
{
if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL)
{
_terminal=true;
-
+
}
- vector<double> nodes (nbelems*2);
+ double* nodes=new double [nbelems];
+ _elems.resize(nbelems);
for (int i=0; i<nbelems; i++)
{
int elem;
else
elem=i;
- _elems.push_back(elem);
- nodes[i*2]=bbs[elem*dim*2+(level%dim)*2];
- nodes[i*2+1]=bbs[elem*dim*2+(level%dim)*2+1];
+ _elems[i]=elem;
+ nodes[i]=bbs[elem*dim*2+(level%dim)*2];
}
- if (_terminal) return;
-
- //elem list is not useful if the node is not terminal
- _elems.clear();
+ if (_terminal) { delete[] nodes; return;}
- vector<double>::iterator median = nodes.begin()+nbelems;
- nth_element< vector<double>::iterator>(nodes.begin(), nodes.begin()+nbelems, nodes.end());
- _median=*median;
+ nth_element< double*>(nodes, nodes+nbelems/2, nodes+nbelems);
+ double median = *(nodes+nbelems/2);
+ delete[] nodes;
// std:: cout << *median <<std::endl;
vector<int> new_elems_left;
vector<int> new_elems_right;
- vector<int> new_elems_center;
+
+ new_elems_left.reserve(nbelems/2+1);
+ new_elems_right.reserve(nbelems/2+1);
+ double max_left=-HUGE;
+ double min_right=HUGE;
for (int i=0; i<nbelems;i++)
{
int elem;
else
elem=i;
double max=bbs[elem*dim*2+(level%dim)*2+1];
- if (max < *median)
- {
- new_elems_left.push_back(elem);
- continue;
- }
double min = bbs[elem*dim*2+(level%dim)*2];
- if (min>*median)
+
+ if (min>median)
{
new_elems_right.push_back(elem);
- continue;
+ if (min<min_right) min_right = min;
+ }
+ else
+ {
+ new_elems_left.push_back(elem);
+ if (max>max_left) max_left = max;
}
- new_elems_center.push_back(elem);
- }
+ }
+ _max_left=max_left;
+ _min_right=min_right;
_left=new BBTree(bbs, &(new_elems_left[0]), level+1, new_elems_left.size());
_right=new BBTree(bbs, &(new_elems_right[0]), level+1, new_elems_right.size());
- _center=new BBTree(bbs, &(new_elems_center[0]), level+1, new_elems_center.size());
-
}
{
if (_left!=0) delete _left;
if (_right!=0) delete _right;
- if (_center!=0) delete _center;
}
-void getIntersectingElems(double* bb, vector<int>& elems)
+void getIntersectingElems(const double* bb, vector<int>& elems) const
{
// terminal node : return list of elements intersecting bb
if (_terminal)
{
- for (int i=0; i<_elems.size(); i++)
+ for (int i=0; i<_nbelems; i++)
{
- double* bb_ptr=_bb+_elems[i]*2*dim;
+ const double* const bb_ptr=_bb+_elems[i]*2*dim;
+// cout << "trying "<<_elems[i]<<endl;
bool intersects = true;
for (int idim=0; idim<dim; idim++)
{
- if (bb_ptr[idim*2]>bb[idim*2+1] || bb_ptr[idim*2+1]<bb[idim*2])
+ if (bb_ptr[idim*2]-bb[idim*2+1]>1e-12|| bb_ptr[idim*2+1]-bb[idim*2]<-1e-12)
intersects=false;
}
if (intersects)
- elems.push_back(_elems[i]);
+ {
+ // std::c4ut << "in BB "<<_elems[i]<<std::endl;
+ elems.push_back(_elems[i]);
+ }
}
return;
}
//non terminal node
double min = bb[(_level%dim)*2];
double max = bb[(_level%dim)*2+1];
- if (max < _median)
+ if (max < _min_right)
{
_left->getIntersectingElems(bb, elems);
- _center->getIntersectingElems(bb,elems);
return;
}
- if (min > _median)
+ if (min> _max_left)
{
_right->getIntersectingElems(bb,elems);
- _center->getIntersectingElems(bb,elems);
return;
}
_left->getIntersectingElems(bb,elems);
- _center->getIntersectingElems(bb,elems);
_right->getIntersectingElems(bb,elems);
}
-};
+
+int size()
+{
+ if (_terminal) return _nbelems;
+ return _left->size()+_right->size();
+}};
--- /dev/null
+#include "MEDMEM_Mesh.hxx"
+#include "InterpolationUtils.hxx"
+#include "PlanarIntersector.hxx"
+#include "PolygonAlgorithms.hxx"
+#include "ConvexIntersector.hxx"
+
+using namespace MED_EN;
+using namespace MEDMEM;
+
+namespace MEDMEM
+{
+ template <int DIM>
+ ConvexIntersector<DIM>::ConvexIntersector(const MEDMEM::MESH& mesh_A, const MEDMEM::MESH& mesh_B,
+ double DimCaracteristic, double Precision,
+ double MedianPlane, bool do_rotate , int PrintLevel)
+ :PlanarIntersector(mesh_A,mesh_B, DimCaracteristic, Precision, MedianPlane)
+ {
+ _Connect_A= mesh_A.getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,
+ MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ _Connect_B= mesh_B.getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,
+ MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ _Conn_index_A= mesh_A.getConnectivityIndex(MED_EN::MED_NODAL, MED_EN::MED_CELL);
+ _Conn_index_B= mesh_B.getConnectivityIndex(MED_EN::MED_NODAL, MED_EN::MED_CELL);
+ _Coords_A = mesh_A.getCoordinates(MED_EN::MED_FULL_INTERLACE);
+ _Coords_B = mesh_B.getCoordinates(MED_EN::MED_FULL_INTERLACE);
+ _Precision = Precision;
+ _Epsilon = Precision*DimCaracteristic;
+ _MedianPlane = MedianPlane;
+ _Do_rotate = do_rotate;
+ _PrintLevel= PrintLevel;
+
+// vector< vector< const double *, int > _Ordered_nodes_A;
+// vector< vector< const double *, int > _Ordered_nodes_B;
+// vector< vector< double > > _Coords_A;
+// vector< vector< double > > _Coords_B;
+// vector< int > _Nb_dist_nodes_A;
+// vector< int > _Nb_dist_nodes_B;
+
+ if( _PrintLevel >= 1)
+ {
+ cout<< " _intersection_type = convex " <<endl;
+ if(DIM==3){
+ if(_Do_rotate) cout<< " _do_rotate = true"<<endl;
+ else cout<< " _do_rotate = false"<<endl;
+ }
+ }
+ }
+
+ template <int DIM>
+ double ConvexIntersector<DIM>::intersectCells(int icell_A, int icell_B,
+ int nb_NodesA, int nb_NodesB)
+ {
+ double result=0;
+
+ /*** Obtain the coordinates of A and B ***/
+ vector< double> Coords_A (DIM*nb_NodesA);
+ vector< double> Coords_B (DIM*nb_NodesB);
+ int nb_dist_NodesA=nb_NodesA;
+ int nb_dist_NodesB=nb_NodesB;
+ int i_last = nb_NodesA - 1;
+ const double * Pi_last= _Coords_A + DIM*(_Connect_A[_Conn_index_A[icell_A-1]-1+i_last]-1);
+
+ for (int i_A=0; i_A<nb_NodesA; i_A++)
+ {
+ const double * Pi_A = _Coords_A + DIM*(_Connect_A[_Conn_index_A[icell_A-1]-1+i_A]-1);
+ if(distance2<DIM>(Pi_last, Pi_A)> _Epsilon)
+ {
+ for (int idim=0; idim<DIM; idim++) Coords_A[DIM*i_A+idim] = *(Pi_A+idim);
+ i_last=i_A; Pi_last = Pi_A;
+ }
+ else nb_dist_NodesA--;
+ }
+ i_last = nb_NodesB - 1;
+ Pi_last= _Coords_B + DIM*(_Connect_B[_Conn_index_B[icell_B-1]-1+i_last]-1);
+ for (int i_B=0; i_B<nb_NodesB; i_B++)
+ {
+ const double * Pi_B = _Coords_B + DIM*(_Connect_B[_Conn_index_B[icell_B-1]-1+i_B]-1);
+ if(distance2<DIM>(Pi_last, Pi_B)> _Epsilon)
+ {
+ for (int idim=0; idim<DIM; idim++) Coords_B[DIM*i_B+idim] = *(Pi_B + idim);
+ i_last=i_B; Pi_last = Pi_B;
+ }
+ else nb_dist_NodesB--;
+ }
+
+ /*** project cells A and B on the median plane ***/
+ /*** and rotate the median plane ***/
+ if(DIM==3) projection(Coords_A, Coords_B, nb_dist_NodesA, nb_dist_NodesB, _Epsilon,_MedianPlane,_Do_rotate);
+
+ /*** Compute the intersection area ***/
+ INTERP_UTILS::PolygonAlgorithms<DIM> P( _Epsilon, _Precision);
+ deque<double> inter = P.PolygonAlgorithms<DIM>::intersect_convex_polygons(&Coords_A[0], &Coords_B[0],
+ nb_dist_NodesA, nb_dist_NodesB);
+ double area[DIM];
+ int nb_inter =((int)inter.size())/DIM;
+ for(int i = 1; i<nb_inter-1; i++)
+ {
+ INTERP_UTILS::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
+ result +=0.5*norm<DIM>(area);
+ }
+
+ //DEBUG prints
+ if( _PrintLevel >= 3)
+ {
+ cout<<endl<< "Cell coordinates after projection"<<endl;
+ cout<<endl<< "icell_A= " << icell_A << ", nb nodes A= " << nb_NodesA << endl;
+ for(int i_A =0; i_A< nb_dist_NodesA; i_A++)
+ {for (int idim=0; idim<DIM; idim++) cout<< Coords_A[DIM*i_A+idim]<< " "; cout << endl;}
+ cout<<endl<< "icell_B= " << icell_B << ", nb nodes B= " << nb_NodesB << endl;
+ for(int i_B =0; i_B< nb_dist_NodesB; i_B++)
+ {for (int idim=0; idim<DIM; idim++) cout<< Coords_B[DIM*i_B+idim]<< " "; cout << endl;}
+ cout<<endl<<"Number of nodes of the intersection = "<< nb_inter <<endl;
+ for(int i=0; i< nb_inter; i++)
+ {for (int idim=0; idim<DIM; idim++) cout<< inter[DIM*i+idim]<< " "; cout << endl;}
+ cout<<endl<<"Intersection area = " << result << endl;
+ exit(1);
+ }
+
+ return result;
+ }
+
+}
--- /dev/null
+#ifndef _CONVEX_INTERSECTOR_HXX_
+#define _CONVEX_INTERSECTOR_HXX_
+
+#include "MEDMEM_Mesh.hxx"
+#include "InterpolationUtils.hxx"
+#include "PolygonAlgorithms.hxx"
+
+namespace MEDMEM
+{
+ template <int DIM>
+ class ConvexIntersector: public PlanarIntersector
+ {
+ public:
+ ConvexIntersector(const MEDMEM::MESH& mesh_A, const MEDMEM::MESH& mesh_B,
+ const double DimCaracteristic, double Precision, double MedianPlane,
+ bool do_rotate, int PrintLevel);
+ double intersectCells(int icell_A, int icell_B, int nb_NodesA, int nb_NodesB);
+
+ private :
+ const int * _Connect_A;
+ const int * _Connect_B;
+ const int *_Conn_index_A;
+ const int *_Conn_index_B;
+ const double * _Coords_A;
+ const double * _Coords_B;
+// vector< vector< const double *, int > _Ordered_nodes_A;
+// vector< vector< const double *, int > _Ordered_nodes_B;
+// vector< vector< double > > _Coords_A;
+// vector< vector< double > > _Coords_B;
+// vector< int > _Nb_dist_nodes_A;
+// vector< int > _Nb_dist_nodes_B;
+ double _Epsilon;
+ double _Precision;
+ double _MedianPlane;
+ bool _Do_rotate;
+ int _PrintLevel;
+ };
+}
+
+#endif
--- /dev/null
+#include "MEDMEM_Mesh.hxx"
+#include "InterpolationUtils.hxx"
+#include "PlanarIntersector.hxx"
+#include "PolygonAlgorithms.hxx"
+#include "GenericIntersector.hxx"
+
+using namespace MED_EN;
+using namespace MEDMEM;
+
+namespace MEDMEM
+{
+ template <int DIM>
+ GenericIntersector<DIM>::GenericIntersector(const MEDMEM::MESH& mesh_A, const MEDMEM::MESH& mesh_B,
+ double DimCaracteristic, double Precision,
+ double MedianPlane, bool do_rotate, int PrintLevel )
+ :PlanarIntersector(mesh_A,mesh_B, DimCaracteristic, Precision, MedianPlane)
+ {
+ _Connect_A= mesh_A.getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,
+ MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ _Connect_B= mesh_B.getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,
+ MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ _Conn_index_A= mesh_A.getConnectivityIndex(MED_EN::MED_NODAL, MED_EN::MED_CELL);
+ _Conn_index_B= mesh_B.getConnectivityIndex(MED_EN::MED_NODAL, MED_EN::MED_CELL);
+ _Coords_A = mesh_A.getCoordinates(MED_EN::MED_FULL_INTERLACE);
+ _Coords_B = mesh_B.getCoordinates(MED_EN::MED_FULL_INTERLACE);
+ _Precision = Precision;
+ _Epsilon = Precision*DimCaracteristic;
+ _MedianPlane = MedianPlane;
+ _Do_rotate = do_rotate;
+ _PrintLevel= PrintLevel;
+
+ if( _PrintLevel >= 1)
+ {
+ cout<< " _intersection_type = Generic " <<endl;
+ if(_Do_rotate) cout<< " _do_rotate = true"<<endl;
+ else cout<< " _do_rotate = false"<<endl;
+ }
+ }
+
+ template <int DIM>
+ double GenericIntersector<DIM>::intersectCells(int icell_A, int icell_B,
+ int nb_NodesA, int nb_NodesB)
+ {
+ double result=0;
+
+ /******* Obtain the coordinates of A and B ********/
+ vector< double> Coords_A (DIM*nb_NodesA);
+ vector< double> Coords_B (DIM*nb_NodesB);
+ for (int idim=0; idim<DIM; idim++)
+ {
+ for (int i_A=0; i_A<nb_NodesA; i_A++)
+ Coords_A[DIM*i_A+idim] = _Coords_A[DIM*(_Connect_A[_Conn_index_A[icell_A-1]-1+i_A]-1)+idim];
+ for (int i_B=0; i_B<nb_NodesB; i_B++)
+ Coords_B[DIM*i_B+idim] = _Coords_B[DIM*(_Connect_B[_Conn_index_B[icell_B-1]-1+i_B]-1)+idim];
+ }
+
+ /*** project cells A and B on the median plane ***/
+ /*** and rotate the median plane ***/
+ if(DIM==3) projection(Coords_A, Coords_B, nb_NodesA, nb_NodesB, _Epsilon,_MedianPlane,_Do_rotate);
+
+ /**********Decompose each polygon into convex component***********/
+ INTERP_UTILS::PolygonAlgorithms<DIM> P(_Epsilon, _Precision);
+ vector< map< int,int > > components_A (0), components_B(0);
+ vector< int > components_index_A(0), components_index_B(0);
+ int nb_comp_A =P.convex_decomposition( &Coords_A[0], nb_NodesA,
+ components_A, components_index_A, _Epsilon);
+ int nb_comp_B =P.convex_decomposition( &Coords_B[0], nb_NodesB,
+ components_B, components_index_B, _Epsilon);
+
+ if( _PrintLevel >= 3)
+ {
+ cout<<endl<< "Cell coordinates after projection"<<endl;
+ cout<< "icell_A= " << icell_A << ", nb nodes A= " << nb_NodesA << endl;
+ for(int i_A =0; i_A< nb_NodesA; i_A++)
+ for (int idim=0; idim<DIM; idim++) cout<< Coords_A[DIM*i_A+idim]<< " "; cout << endl;
+ cout<< "icell_B= " << icell_B << ", nb nodes B= " << nb_NodesB << endl;
+ for(int i_B =0; i_B< nb_NodesB; i_B++)
+ for (int idim=0; idim<DIM; idim++) cout<< Coords_B[DIM*i_B+idim]<< " "; cout << endl;
+ }
+
+ /*** Compute the intersection area ***/
+ for(int i_comp_A = 0; i_comp_A<nb_comp_A; i_comp_A++)
+ {
+ int sign_A = components_index_A[i_comp_A] <0 ? -1:1;
+ int nb_nodes_A = sign_A * components_index_A[i_comp_A];
+ double * Coords_subP_A;
+ map< int,int >::iterator mi_A = (components_A[i_comp_A]).begin();
+
+ for(int i_comp_B = 0; i_comp_B<nb_comp_B; i_comp_B++)
+ {
+ int sign_B = components_index_B[i_comp_B] <0 ? -1:1;
+ int nb_nodes_B = sign_B * components_index_B[i_comp_B];
+ double * Coords_subP_B;
+ map< int,int >::iterator mi_B = (components_B[i_comp_B]).begin();
+
+ for (int i_A=0; i_A< nb_nodes_A; i_A++)
+ {
+ for (int idim=0; idim<DIM; idim++)
+ Coords_subP_A[DIM*i_A+idim] = Coords_A[DIM*(* mi_A).second+idim];
+ mi_A++;
+ }
+ for (int i_B=0; i_B< nb_nodes_B; i_B++)
+ {
+ for (int idim=0; idim<DIM; idim++)
+ Coords_subP_B[DIM*i_B+idim] = Coords_B[DIM*(* mi_B).second+idim];
+ mi_B++;
+ }
+
+ deque<double> inter =
+ P.PolygonAlgorithms<DIM>::intersect_convex_polygons(&Coords_subP_A[0], &Coords_subP_B[0],
+ nb_nodes_A, nb_nodes_B);
+ double area[DIM];
+ int nb_inter =((int)inter.size())/DIM;
+ for(int i = 1; i<nb_inter-1; i++)
+ {
+ INTERP_UTILS::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
+ result += 0.5 * norm<DIM>(area)*sign_A*sign_B;
+ }
+ //DEBUG prints
+ if(_PrintLevel >= 3)
+ {
+ cout<<endl<<"Number of nodes of the intersection = "<< nb_inter <<endl<<endl;
+ for(int i=0; i< nb_inter; i++)
+ for (int idim=0; idim<2; idim++) cout<< inter[2*i+idim]<< " "; cout << endl;
+ }
+ }
+ }
+
+ return result;
+ }
+
+}
--- /dev/null
+#ifndef _GENERIC_INTERSECTOR_HXX_
+#define _GENERIC_INTERSECTOR_HXX_
+
+#include "MEDMEM_Mesh.hxx"
+#include "InterpolationUtils.hxx"
+
+namespace MEDMEM
+{
+ class MESH;
+ template <int DIM>
+ class GenericIntersector: public PlanarIntersector
+ {
+ public:
+ GenericIntersector (const MEDMEM::MESH& mesh_A, const MEDMEM::MESH& mesh_B,
+ const double DimCaracteristic, double Precision, double MedianPlane,
+ bool do_rotate, int PrintLevel);
+ double intersectCells(int icell_A, int icell_B, int nb_NodesA, int nb_NodesB);
+
+ private :
+ const int * _Connect_A;
+ const int * _Connect_B;
+ const int *_Conn_index_A;
+ const int *_Conn_index_B;
+ const double * _Coords_A;
+ const double * _Coords_B;
+ double _Epsilon;
+ double _Precision;
+ double _MedianPlane;
+ bool _Do_rotate;
+ int _PrintLevel;
+ };
+}
+
+#endif
namespace MEDMEM
{
-
-class Interpolation
-{
-public:
- Interpolation(){}
- virtual ~Interpolation(){}
-
- //interpolation of two triangular meshes.
- virtual std::vector<std::map<int, double> > interpol_maillages(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2)=0;
-
+ typedef enum {Triangulation, Convex, Generic} IntersectionType;
+
+ class Interpolation
+ {
+ public:
+ Interpolation(){}
+ virtual ~Interpolation(){}
+
+ //interpolation of two triangular meshes.
+ virtual std::vector<std::map<int, double> > interpolateMeshes(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2)=0;
+
};
};
#include "InterpolationUtils.hxx"
+#include "TranslationRotationMatrix.hxx"
#include "Interpolation.hxx"
#include "Interpolation2D.hxx"
namespace MEDMEM
{
- Interpolation2D::Interpolation2D()
+ Interpolation2D::Interpolation2D():_counter(0)
{
_Precision=1.0E-12;
_DimCaracteristic=1;
//debug cout << "mailles sources testées : " << N_Maille << endl;
- vector<double> Inter=INTERP_UTILS::intersec_de_triangle(MaS._Coord+2*(NS_1-1),MaS._Coord+2*(NS_2-1),
-
- MaS._Coord+2*(NS_3-1),MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1), _DimCaracteristic, _Precision);
-
+ vector<double> Inter;
+ INTERP_UTILS::intersec_de_triangle(MaS._Coord+2*(NS_1-1),MaS._Coord+2*(NS_2-1),
+ MaS._Coord+2*(NS_3-1),MaP._Coord+2*(NP_1-1),
+ MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1),
+ Inter,_DimCaracteristic, _Precision);
+ _counter++;
if(Inter.size()>0)
{
//debug cout << "Localise : maille proche barycentre trouvée : " << N_Maille << endl;
int N_3=MaS._Connect[3*(N_Maille-1)+2];
//debug cout << "mailles sources testées : " << N_Maille << endl;
- vector<double> Inter=INTERP_UTILS::intersec_de_triangle(MaS._Coord+2*(N_1-1),MaS._Coord+2*(N_2-1),
- MaS._Coord+2*(N_3-1),MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1) , _DimCaracteristic, _Precision);
-
+ vector<double> Inter;
+ INTERP_UTILS::intersec_de_triangle(MaS._Coord+2*(N_1-1),MaS._Coord+2*(N_2-1),
+ MaS._Coord+2*(N_3-1),MaP._Coord+2*(NP_1-1),
+ MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1),
+ Inter, _DimCaracteristic, _Precision);
+ _counter++;
if(Inter.size()>0)
{
//debug cout << "Localise : maille proche sommet " << Num+1 << " trouvée : " << N_Maille << endl;
vector<map<int,double> >& Surface_d_intersect,FIELD<double>& myField_P)
{
-
+
//on récupere le n° des noeuds.
//debug cout << "\nintersect maille " << i_P << endl;
int NS_3=MaS._Connect[3*(i_S-1)+2];
- vector<double> Inter=INTERP_UTILS::intersec_de_triangle(MaS._Coord+2*(NS_1-1),MaS._Coord+2*(NS_2-1),MaS._Coord+2*(NS_3-1),MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1), _DimCaracteristic, _Precision);
+ vector<double> Inter;
+ INTERP_UTILS::intersec_de_triangle(MaS._Coord+2*(NS_1-1),MaS._Coord+2*(NS_2-1),
+ MaS._Coord+2*(NS_3-1),MaP._Coord+2*(NP_1-1),
+ MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1),
+ Inter,_DimCaracteristic, _Precision);
+ _counter++;
//on teste l'intersection.
int taille=Inter.size()/2;
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
//fonction principale pour interpoler deux maillages triangulaires.
- vector<map<int,double> > Interpolation2D::interpol_maillages(const MEDMEM::MESH& myMesh_S,
+ vector<map<int,double> > Interpolation2D::interpolateMeshes(const MEDMEM::MESH& myMesh_S,
const MEDMEM::MESH& myMesh_P)
{
-
+ long global_start =clock();
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
// on vérifie que les deux maillages sont formés de triangles.
int NumberOfCellsTypes_S = myMesh_S.getNumberOfTypes(MED_CELL);
//on cherche pour chaque maille i_P du maillage projetté
// une maille i_S du maillage S qui chevauche la maille i_P.
+ long start_filtering=clock();
vector<int> localis;
localis.reserve(MaP._NbMaille);
MEDMEM::MESH* ptr_S = const_cast<MEDMEM::MESH*>(&myMesh_S);
cout << "Interpolation2D::local_iP_dansS"<<endl;
local_iP_dans_S(*ptr_S,*ptr_P,MaP,MaS,localis);
+ long end_filtering=clock();
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
//on crée un tableau où l'on rentrera la valeurs des intersections.
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
//boucle sur les mailles de P.
//Coeur de l'algorithme
-
+ long start_intersection=clock();
for(int i_P=0; i_P<MaP._NbMaille ; i_P++)
{
int i_S_init=localis[i_P];
rempli_vect_d_intersect(i_P+1,i_S_init,MaP,MaS,Surface_d_intersect,myField_P);
}
-
- if (_PrintLevel >= 2)
+ if (_PrintLevel >= 1)
{
- cout<<endl<<"Impression des surfaces d'intersection:"<<endl<<endl;
- cout<<"(maille source, maille cible):surface d'intersection"<<endl;
- for(int i=0; i< Surface_d_intersect.size();i++)
- {
- if(Surface_d_intersect[i].size()!=0)
- {
+ long end_intersection=clock();
+ if (_PrintLevel >= 2)
+ {
+ cout<<endl<<"Printing intersection areas:"<<endl<<endl;
+ cout<<"(source cell, target cell): intersection areas"<<endl;
+ double total=0.0;
+ double total_interm=0.0;
+ int nb_result_areas = Surface_d_intersect.size();
+ for(int i=0; i< nb_result_areas;i++)
+ {
map<int,double>::iterator surface;
+ total_interm=0.0;
for( surface=Surface_d_intersect[i].begin();surface!=Surface_d_intersect[i].end();surface++)
- cout<<" ("<<i+1<<" , "<<(*surface).first<<")"<<" : "<<(*surface).second<<endl;
+ {
+ cout<<" ("<<i+1<<" , "<<(*surface).first<<")"<<" : "<<(*surface).second<<endl;
+ total_interm +=(*surface).second;
+ }
+ cout<< " elem " << i+1 << " area= " << total_interm << endl;
+ total+=total_interm;
}
- }
+ cout << "total area "<<total<<endl;
+ }
+ cout<< "Barycenter localisation time= " << end_filtering-start_filtering<< endl;
+ cout<< "Intersection time= " << end_intersection-start_intersection<< endl;
+ cout<< "counter= " << _counter << endl;
+ long global_end =clock();
+ cout<< "Global time= " << global_end - global_start << endl;
}
-
+
// /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
// //On sauvegarde le champ crée sur la maillage P.
// if (_PrintLevel>=3)
#include "MEDMEM_Mesh.hxx"
#include "MEDMEM_Field.hxx"
-
+#include "Interpolation.hxx"
#include <cmath>
#include <map>
namespace INTERP_UTILS
{
- struct monMaillageS;
- struct monMaillageP;
+ struct monMaillageS;
+ struct monMaillageP;
};
namespace MEDMEM
{
-
- class Interpolation;
-
-class Interpolation2D : public Interpolation
-{
-
- private:
- double _Precision;
- double _DimCaracteristic;
- int _SecondMethodOfLocalisation;
- int _PrintLevel;
-
- // Méthodes publiques
- public:
-
- Interpolation2D();
-
- // precision geometrique, choix de la methode comlementaire de localisation, niveau d'impressions
- void setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel);
-
-
- //localise le barycentre des mailles de P dans le maillage S
- void local_iP_dans_S(MEDMEM::MESH& myMesh_S,MEDMEM::MESH& myMesh_P,
- INTERP_UTILS::monMaillageP& MaP,INTERP_UTILS::monMaillageS& MaS,vector<int>& localis );
-
-
- //fonction qui permet de remplir le vecteur donnant la surface d'intersection
- //de la mailles i_P du maillage projetté avec la maille i_S du maillage support.
- inline void rempli_vect_d_intersect(int i_P,int i_S,const INTERP_UTILS::monMaillageP& MaP,const INTERP_UTILS::monMaillageS& MaS,
- vector<map<int,double> >& Surface_d_intersect,
- FIELD<double>& myField_P);
-
-
- //fonction principale pour interpoler deux maillages triangulaires.
- std::vector<map<int, double> > interpol_maillages(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
-
- private:
-
-
-};
-
+
+
+ class Interpolation2D : public Interpolation
+ {
+
+ private:
+ double _Precision;
+ double _DimCaracteristic;
+ int _SecondMethodOfLocalisation;
+ int _PrintLevel;
+ IntersectionType _intersection_type;
+ int _counter;
+ // Méthodes publiques
+ public:
+
+ Interpolation2D();
+
+ // precision geometrique, choix de la methode comlementaire de localisation, niveau d'impressions
+ void setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel);
+
+ //declare Intersection polygons as convex
+ void setConvex(bool flag)
+ {
+ if (flag) _intersection_type=Convex;
+ }
+
+ //localise le barycentre des mailles de P dans le maillage S
+ void local_iP_dans_S(MEDMEM::MESH& myMesh_S,MEDMEM::MESH& myMesh_P,
+ INTERP_UTILS::monMaillageP& MaP,INTERP_UTILS::monMaillageS& MaS,vector<int>& localis );
+
+
+ //fonction qui permet de remplir le vecteur donnant la surface d'intersection
+ //de la mailles i_P du maillage projetté avec la maille i_S du maillage support.
+ inline void rempli_vect_d_intersect(int i_P,int i_S,
+ const INTERP_UTILS::monMaillageP& MaP,const INTERP_UTILS::monMaillageS& MaS,
+ vector<map<int,double> >& Surface_d_intersect,FIELD<double>& myField_P);
+
+
+ //fonction principale pour interpoler deux maillages triangulaires.
+ std::vector<map<int, double> >interpolateMeshes(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
+
+ };
+
};
#endif
#include "Interpolation3D.hxx"
#include "MeshElement.hxx"
#include "TransformedTriangle.hxx"
-//#include "VectorUtils.hxx"
#include "IntersectorTetra.hxx"
#include "IntersectorHexa.hxx"
#include "TargetIntersector.hxx"
* @return vector containing for each element i of the source mesh, a map giving for each element j
* of the target mesh which i intersects, the volume of the intersection
*/
- IntersectionMatrix Interpolation3D::interpol_maillages(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh)
+ IntersectionMatrix Interpolation3D::interpolateMeshes(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh)
{
// it seems wasteful to make a copy here
IntersectionMatrix matrix;
Interpolation3D();
virtual ~Interpolation3D();
- virtual IntersectionMatrix interpol_maillages(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh);
+ virtual IntersectionMatrix interpolateMeshes(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh);
private:
-#include "MEDMEM_Field.hxx"
#include "MEDMEM_Mesh.hxx"
-#include "TranslationRotationMatrix.hxx"
-#include "InterpolationUtils.hxx"
-#include "Interpolation.hxx"
#include "Interpolation3DSurf.hxx"
+#include "PlanarIntersector.hxx"
+#include "PlanarIntersector.H"
+#include "TriangulationIntersector.hxx"
+#include "TriangulationIntersector.cxx"
+#include "ConvexIntersector.hxx"
+#include "ConvexIntersector.cxx"
+#include "GenericIntersector.hxx"
+#include "GenericIntersector.cxx"
#include "BBTree.H"
+#include<time.h>
-using namespace std;
using namespace MED_EN;
-using namespace MEDMEM;
-using namespace INTERP_UTILS;
namespace MEDMEM
{
Interpolation3DSurf::Interpolation3DSurf()
- {
+ {
+ _Intersection_type= MEDMEM::Triangulation;
+ _MedianPlane=0.5;
+ _Do_rotate=true;
_Precision=1.0E-12;
_DimCaracteristic=1;
- _SecondMethodOfLocalisation=1;
+ _Surf3DAdjustmentEps=0.0001;
_PrintLevel=0;
- _Surf3DAdjustmentEps=0.1;
}
-
+
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
/* Options : */
/* Precision : for geometric computation */
- /* SecondMethodOfLocalisation : 0/1 */
/* PrintLevel : between 0 and 3 */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- void Interpolation3DSurf::setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel)
+ void Interpolation3DSurf::setOptions(double Precision, int PrintLevel, double MedianPlane, IntersectionType intersection_type, bool do_rotate)
{
+ _Intersection_type=intersection_type;
+ _MedianPlane=MedianPlane;
+ _Do_rotate=do_rotate;
_Precision=Precision;
- _SecondMethodOfLocalisation=SecondMethodOfLocalisation;
_PrintLevel=PrintLevel;
- }
-
-
-
- /*! method computing the translation/rotation matrix
- necessary to transform a triangle into a triangle located inside the Oxy plane. The method fails only if two of the three points are coincident
- \param P1 pointer to three coordinates defining the first point
- \param P2 pointer to three coordinates defining the second point
- \param P3 pointer to three coordinates defining the third point
- \return TranslationRotationMatrix object containing the transform that must be applied to the triangle for putting it inside the Oxy plane.
- */
- void Interpolation3DSurf::rotate3DTriangle(const double* PP1, const double*PP2, const double*PP3, TranslationRotationMatrix& rotation_matrix)
- {
- //initializes
- rotation_matrix.translate(PP1);
-
- double P2w[3];
- double P3w[3];
- P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2];
- P3w[0]=PP3[0]; P3w[1]=PP3[1];P3w[2]=PP3[2];
-
- // translating to set P1 at the origin
- for (int i=0; i<2; i++)
- {
- P2w[i]-=PP1[i];
- P3w[i]-=PP1[i];
- }
-
- // rotating to set P2 on the Oxy plane
- TranslationRotationMatrix A;
- A.rotate_x(P2w);
- A.rotate_vector(P3w);
- rotation_matrix.multiply(A);
-
- //rotating to set P2 on the Ox axis
- TranslationRotationMatrix B;
- B.rotate_z(P2w);
- B.rotate_vector(P3w);
- rotation_matrix.multiply(B);
+ }
- //rotating to set P3 on the Oxy plane
- TranslationRotationMatrix C;
- C.rotate_x(P3w);
- rotation_matrix.multiply(C);
- }
-
-
-
- /* _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- /* fonction qui permet de remplir le vecteur donnant la surface d'intersection */
- /* de la maille i_P du maillage P (projeté) avec la maille i_S du maillage S (support) */
- /* _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-
- inline void Interpolation3DSurf::rempli_vect_d_intersect(int i_P,int i_S,const monMaillageP& MaP,const monMaillageS& MaS,
- vector<map<int,double> >& Surface_d_intersect
- )
-
+ /***************************************************************/
+ /* Main function to interpolate triangular or quadratic meshes */
+ /***************************************************************/
+ vector<map<int,double> > Interpolation3DSurf::interpolateMeshes(const MEDMEM::MESH& myMesh_S,
+ const MEDMEM::MESH& myMesh_P)
{
-
-
- //on récupere le n° des noeuds.
-
- //debug cout << "\nintersect maille " << i_P << endl;
- int NP_1=MaP._Connect[3*(i_P-1)];
- int NP_2=MaP._Connect[3*(i_P-1)+1];
- int NP_3=MaP._Connect[3*(i_P-1)+2];
-
-
- TranslationRotationMatrix rotation;
- rotate3DTriangle(MaP._Coord+3*(NP_1-1),MaP._Coord+3*(NP_2-1), MaP._Coord+3*(NP_3-1), rotation);
-
- //on calcule la surface de la maille i_P
-
- // double Surf_i_P =INTERP_UTILS::Surf_Tri(MaP._Coord+3*(NP_1-1),MaP._Coord+3*(NP_2-1),
- // MaP._Coord+3*(NP_3-1));
-
- // double Surface = 0 ;
+ long global_start =clock();
+ /***********************************************************/
+ /* Check both meshes are made of triangles and quadrangles */
+ /***********************************************************/
+ int NumberOfCellsTypes_S = myMesh_S.getNumberOfTypes(MED_EN::MED_CELL);
+ int NumberOfCellsTypes_P = myMesh_P.getNumberOfTypes(MED_EN::MED_CELL);
+ if ( NumberOfCellsTypes_S > 2 || NumberOfCellsTypes_P >2)
+ cout<<"Interpolation3DSurf::interpolateMeshes: one of the two meshes contains more than two types of cells"<<endl;
+ string* Type_S = myMesh_S.getCellTypeNames(MED_EN::MED_CELL);
+ string* Type_P = myMesh_P.getCellTypeNames(MED_EN::MED_CELL);
+ if((Type_S[0] != "MED_TRIA3" && Type_S[0] != "MED_QUAD4") || (Type_P[0] != "MED_TRIA3" && Type_P[0] != "MED_QUAD4"))
+ cout<<"Interpolation3DSurf::interpolateMeshes: one of the two meshes is neither linear triangular nor linear quadratic"<<endl;
+ //VB
+ delete[] Type_S;
+ delete[] Type_P;
+ int nbMaille_S =myMesh_S.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ int nbMaille_P =myMesh_P.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
- int NS_1=MaS._Connect[3*(i_S-1)];
- int NS_2=MaS._Connect[3*(i_S-1)+1];
- int NS_3=MaS._Connect[3*(i_S-1)+2];
-
- double PS1[3];
- double PS2[3];
- double PS3[3];
- double PP1[3];
- double PP2[3];
- double PP3[3];
-
- for (int i=0; i<3; i++)
- {
- PS1[i] = *(MaS._Coord+3*(NS_1-1)+i);
- PS2[i] = *(MaS._Coord+3*(NS_2-1)+i);
- PS3[i] = *(MaS._Coord+3*(NS_3-1)+i);
- PP1[i] = *(MaP._Coord+3*(NP_1-1)+i);
- PP2[i] = *(MaP._Coord+3*(NP_2-1)+i);
- PP3[i] = *(MaP._Coord+3*(NP_3-1)+i);
- }
-
- rotation.transform_vector(PS1);
- rotation.transform_vector(PS2);
- rotation.transform_vector(PS3);
- rotation.transform_vector(PP1);
- rotation.transform_vector(PP2);
- rotation.transform_vector(PP3);
-
-
-
- //intersects 3D triangle only considering the two first coordinates
- //
- //> MaP is located in the Oxy plane
- //> PS is not and is therefore projected on MaP orthogonally along the z
- // direction
- vector<double> Inter=INTERP_UTILS::
- intersec_de_triangle(PS1,PS2,PS3,
- PP1,PP2,PP3,
- _DimCaracteristic,
- _Precision);
-
-
- //on teste l'intersection.
- int taille=Inter.size()/2;
- //debug cout << " -> maille source : " << i_S << " , nb intersection=" << taille << endl;
+ /**************************************************/
+ /* Search the characteristic size of the meshes */
+ /**************************************************/
+ vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
+ vector<vector<double> > BoxP=myMesh_P.getBoundingBox();
+ double diagonalS=sqrt((BoxS[1][0]-BoxS[0][0])*(BoxS[1][0]-BoxS[0][0])
+ +(BoxS[1][1]-BoxS[0][1])*(BoxS[1][1]-BoxS[0][1])
+ +(BoxS[1][2]-BoxS[0][2])*(BoxS[1][2]-BoxS[0][2]));
+ double DimCaracteristic_S=diagonalS/nbMaille_S;
+ double diagonalP=sqrt((BoxP[1][0]-BoxP[0][0])*(BoxP[1][0]-BoxP[0][0])
+ +(BoxP[1][1]-BoxP[0][1])*(BoxP[1][1]-BoxP[0][1])
+ +(BoxP[1][2]-BoxP[0][2])*(BoxP[1][2]-BoxP[0][2]));
+ double DimCaracteristic_P=diagonalP/nbMaille_P;
- /* CAS 3 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
- if(taille==3)
+ _DimCaracteristic=min(DimCaracteristic_S, DimCaracteristic_P);
+ if (_PrintLevel>=1)
{
- double Surf_inter = Surf_Poly(Inter) ;
-
- //on remplit le vecteur des correspondances n° maille-surfaces d'intersections.
- Surface_d_intersect[i_P-1][i_S]=Surf_inter;
-
-
- // Surface = Surface + Surf_inter ;
+ cout<<" - Characteristic size of the source mesh : "<< DimCaracteristic_S<<endl;
+ cout<<" - Characteristic size of the target mesh: "<< DimCaracteristic_P<<endl;
+ cout << "Interpolation3DSurf::computation of the intersections"<<endl;
}
- /* CAS 4 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ PlanarIntersector* intersector;
- else if (taille>3) //taille>3
+ switch (_Intersection_type)
{
- vector<double> Poly_Inter=reconstruct_polygon(Inter);
-
- double Surf_inter =Surf_Poly(Poly_Inter) ;
-
- //on remplit le vecteur des correspondances n° maille-surfaces d'intersection
- (Surface_d_intersect[i_P-1])[i_S]=Surf_inter;
-
- // Surface = Surface + Surf_inter ;
+ case MEDMEM::Triangulation:
+ intersector=new TriangulationIntersector<3>(myMesh_P,myMesh_S, _DimCaracteristic,_Precision,
+ _MedianPlane, _PrintLevel);
+ break;
+ case MEDMEM::Convex:
+ intersector=new ConvexIntersector<3>(myMesh_P,myMesh_S, _DimCaracteristic, _Precision,
+ _MedianPlane, _Do_rotate, _PrintLevel);
+ break;
+ case MEDMEM::Generic:
+ intersector=new GenericIntersector<3>(myMesh_P,myMesh_S, _DimCaracteristic,_Precision,
+ _MedianPlane, _Do_rotate, _PrintLevel);
+ break;
}
- //debug cout << " surface = " << Surface << endl << flush;
-
-
- /* _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- //on rentre la fraction totale de la maille i_P qui a été considérer lors de l'interpolation.
- // double Fract=Surface/Surf_i_P;
- // myField_P.setValueIJ(i_P,1,Fract);
-
- }
-
-
- /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- //fonction principale pour interpoler deux maillages triangulaires.
- vector<map<int,double> > Interpolation3DSurf::interpol_maillages(const MEDMEM::MESH& myMesh_S,
- const MEDMEM::MESH& myMesh_P)
- {
-
-
- /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- // on vérifie que les deux maillages sont formés de triangles.
- int NumberOfCellsTypes_S = myMesh_S.getNumberOfTypes(MED_CELL);
- int NumberOfCellsTypes_P = myMesh_P.getNumberOfTypes(MED_CELL);
- if ( NumberOfCellsTypes_S != 1 || NumberOfCellsTypes_P != 1)
- { cout<<"l'un au moins des deux maillages n'est pas triangulaires"<<endl;}
- string* Type_S = myMesh_S.getCellTypeNames(MED_CELL);
- string* Type_P = myMesh_P.getCellTypeNames(MED_CELL);
- if ( Type_S[0] != "MED_TRIA3" || Type_P[0] != "MED_TRIA3")
- { cout<<"l'un au moins des deux maillages n'est pas triangulaires"<<endl;}
- //VB
- delete[] Type_S;
- delete[] Type_P;
-
- /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- INTERP_UTILS::monMaillageS MaS;
- INTERP_UTILS::monMaillageP MaP;
-
- MaS._NbNoeud = myMesh_S.getNumberOfNodes() ;
- MaP._NbNoeud = myMesh_P.getNumberOfNodes() ;
- MaS._NbMaille =myMesh_S.getNumberOfElements(MED_CELL,MED_TRIA3);
- MaP._NbMaille =myMesh_P.getNumberOfElements(MED_CELL,MED_TRIA3);
-
- //on charge les connectivités des deux maillages.
- MaS._Connect =myMesh_S.getConnectivity(MED_FULL_INTERLACE,
- MED_NODAL, MED_CELL, MED_TRIA3) ;
- MaP._Connect =myMesh_P.getConnectivity(MED_FULL_INTERLACE,
- MED_NODAL, MED_CELL, MED_TRIA3) ;
-
- //on charge les coordonnées des noeuds.
- MaS._Coord = myMesh_S.getCoordinates(MED_FULL_INTERLACE);
- MaP._Coord = myMesh_P.getCoordinates(MED_FULL_INTERLACE);
-
- /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- //on cherche la dimension caractéristique des maillages
- /*______________________________________________________________*/
-
- vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
- vector<vector<double> > BoxP=myMesh_P.getBoundingBox();
- double VolumeS=sqrt((BoxS[1][0]-BoxS[0][0])*(BoxS[1][0]-BoxS[0][0])+(BoxS[1][1]-BoxS[0][1])*(BoxS[1][1]-BoxS[0][1]))+(BoxS[1][2]-BoxS[0][2])*(BoxS[1][2]-BoxS[0][2]);
- MaS._DimCaracteristic=sqrt(VolumeS/MaS._NbMaille);
- double VolumeP=sqrt((BoxP[1][0]-BoxP[0][0])*(BoxP[1][0]-BoxP[0][0])+(BoxP[1][1]-BoxP[0][1])*(BoxP[1][1]-BoxP[0][1]))+(BoxP[1][2]-BoxP[0][2])*(BoxP[1][2]-BoxP[0][2]);
- MaP._DimCaracteristic=sqrt(VolumeP/MaP._NbMaille);
-
- _DimCaracteristic=min(MaS._DimCaracteristic,MaP._DimCaracteristic);
- if (_PrintLevel)
- {
- cout<<"recherche de la dimension caractéristique des maillages 3D_surf :"<<endl;
- cout<<" - dimension caractéristique du maillage source : "<<MaS._DimCaracteristic<<endl;
- cout<<" - dimension caractéristique du maillage projeté: "<<MaS._DimCaracteristic<<endl;
- cout<<" - d'où la dimension caractéristique: "<<_DimCaracteristic<<endl;
- }
-
- //creating a search structure based on the bounding boxes
- //of the elements to know
-
- vector<double> bbox;
- createBoundingBoxes(MaS,bbox); // create the bounding boxes
- adjustBoundingBoxes(bbox); // expand them so that no intersection is missed
+ /****************************************************************/
+ /* Create a search tree based on the bounding boxes */
+ /* Instanciate the intersector and initialise the result vector */
+ /****************************************************************/
+
+ long start_filtering=clock();
- BBTree<3> my_tree(&bbox[0], 0, 0,MaS._NbMaille );//creating the search structure
+ vector<double> bbox;
+ intersector->createBoundingBoxes<3>(myMesh_S,bbox); // create the bounding boxes
+ intersector->adjustBoundingBoxes<3>(bbox, _Surf3DAdjustmentEps); // expand them so that no intersection is missed
+ BBTree<3> my_tree(&bbox[0], 0, 0,nbMaille_S );//creating the search structure
+
+ long end_filtering=clock();
- /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- //on crée un tableau où l'on rentrera la valeurs des intersections.
- cout << "Interpolation3DSurf::calcul intersec"<<endl;
map<int,double> MAP_init;
- vector<map<int,double> > Surface_d_intersect(MaP._NbMaille,MAP_init);//on initialise.
+ vector<map<int,double> > result_areas(nbMaille_P,MAP_init);//on initialise.
+ /****************************************************/
+ /* Loop on the target cells - core of the algorithm */
+ /****************************************************/
+ const MED_EN::medGeometryElement* types = myMesh_P.getTypes(MED_EN::MED_CELL);
+ int type_max_index=0;//maximum cell number for a given type
+ int type_min_index=0;//minimum cell number for a given type
+ int i_P=0;//global index of cell
- /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- //boucle sur les mailles de P.
- //Coeur de l'algorithme
+ long start_intersection=clock();
- for(int i_P=0; i_P<MaP._NbMaille ; i_P++)
+ for (int itype = 0; itype<NumberOfCellsTypes_P; itype++)
{
- vector<int> intersecting_elems;
- double bb[6];
- INTERP_UTILS::getElemBB(bb,MaP,i_P);
- my_tree.getIntersectingElems(bb, intersecting_elems);
-
- //browsing all the i_S (from mesh S) elems that can
- //intersect elem i_P (from mesh P)
- for (int ielem=0; ielem<intersecting_elems.size();ielem++)
+ int nb_nodesP=types[itype]%100;
+ int nbelem_type = myMesh_P.getNumberOfElements(MED_EN::MED_CELL, types[itype]);
+ type_max_index += nbelem_type;
+ for( i_P=type_min_index; i_P<type_max_index ; i_P++)
{
- //BBTree structure returns numbers between 0 and n-1
- int i_S=intersecting_elems[ielem]+1;
-
- rempli_vect_d_intersect(i_P+1,i_S,MaP,MaS,Surface_d_intersect);
+ vector<int> intersecting_elems;
+ double bb[6];
+ intersector->getElemBB<3>(bb,myMesh_P,i_P+1,nb_nodesP);
+ my_tree.getIntersectingElems(bb, intersecting_elems);
+ int nb_intersecting_elems = intersecting_elems.size();
+ //browsing all the i_S (from mesh S) elems that can
+ //intersect elem i_P (from mesh P)
+ for (int ielem=0; ielem<nb_intersecting_elems;ielem++)
+ {
+ //BBTree structure returns numbers between 0 and n-1
+ int i_S=intersecting_elems[ielem]; //MN: Global number of cell ?
+ int nb_nodesS=myMesh_S.getElementType(MED_EN::MED_CELL,i_S+1)%100;
+ double surf=intersector->intersectCells(i_P+1,i_S+1,nb_nodesP,nb_nodesS);
+ (result_areas[i_P]).insert(make_pair(i_S+1,surf));
+ //rempli_vect_d_intersect(i_P+1,i_S+1,MaP,MaS,result_areas);
+ }
+ intersecting_elems.clear();
}
- intersecting_elems.clear();
+ type_min_index = type_max_index;
}
+ delete intersector;
+ /***********************************/
+ /* DEBUG prints */
+ /***********************************/
- //DEBUG prints
- if (_PrintLevel >= 2)
+ if (_PrintLevel >=1)
{
- cout<<endl<<"Impression des surfaces d'intersection:"<<endl<<endl;
- cout<<"(maille source, maille cible):surface d'intersection"<<endl;
- double total=0.0;
- for(int i=0; i< Surface_d_intersect.size();i++)
- {
- if(Surface_d_intersect[i].size()!=0)
- {
+ long end_intersection=clock();
+ if (_PrintLevel >=2)
+ {
+ cout<<endl<<"Printing intersection areas:"<<endl<<endl;
+ cout<<"(source cell, target cell): intersection areas"<<endl;
+ double total=0.0;
+ double total_interm=0.0;
+ int nb_result_areas = result_areas.size();
+ for(int i=0; i< nb_result_areas;i++)
+ {
map<int,double>::iterator surface;
- for( surface=Surface_d_intersect[i].begin();surface!=Surface_d_intersect[i].end();surface++)
+ total_interm=0.0;
+ for( surface=result_areas[i].begin();surface!=result_areas[i].end();surface++)
{
cout<<" ("<<i+1<<" , "<<(*surface).first<<")"<<" : "<<(*surface).second<<endl;
- total+=(*surface).second;
+ total_interm +=(*surface).second;
}
+ cout<< " elem " << i+1 << " area= " << total_interm << endl;
+ total+=total_interm;
}
- }
- cout << "surface totale "<<total<<endl;
- }
-
- return Surface_d_intersect;
-
- }
-
-
-
-
-
- /*! Creates a set of bounding boxes that correspond to mesh \a mesh. The vector returned is of dimension 6*nb_elems with bounding boxes stored as xmin1, xmax1, ymin1, ymax1, zmin1, zmax1, xmin2, xmax2, ymin2,...
- The returned pointer must be deleted by the calling code.
-
- \param mesh structure pointing to the mesh
- \param bbox vector containing the bounding boxes
- */
-
- void Interpolation3DSurf::createBoundingBoxes(const struct monMaillageS Ma_S, vector<double>& bbox)
- {
- /* We build the segment tree for locating possible matching intersections*/
- bbox.resize(6* Ma_S._NbMaille);
- for (int i=0; i< Ma_S._NbMaille; i++)
- {
- bbox[i*6]=HUGE;
- bbox[i*6+1]=-HUGE;
- bbox[i*6+2]=HUGE;
- bbox[i*6+3]=-HUGE;
- bbox[i*6+4]=HUGE;
- bbox[i*6+5]=-HUGE;
-
- for (int j=0; j<3; j++)
- {
- double x=Ma_S._Coord[(Ma_S._Connect[i*3+j]-1)*3];
- bbox[i*6]=(bbox[i*6]<x)?bbox[i*6]:x;
- bbox[i*6+1]=(bbox[i*6+1]>x)?bbox[i*6+1]:x;
- double y=Ma_S._Coord[(Ma_S._Connect[i*3+j]-1)*3+1];
- bbox[i*6+2]=(bbox[i*6+2]<y)?bbox[i*6+2]:y;
- bbox[i*6+3]=(bbox[i*6+3]>y)?bbox[i*6+3]:y;
- double z=Ma_S._Coord[(Ma_S._Connect[i*3+j]-1)*3+2];
- bbox[i*6+4]=(bbox[i*6+4]<z)?bbox[i*6+4]:z;
- bbox[i*6+5]=(bbox[i*6+5]>z)?bbox[i*6+5]:z;
+ cout << "total area "<<total<<endl;
}
+ cout<< "Filtering time= " << end_filtering-start_filtering<< endl;
+ cout<< "Intersection time= " << end_intersection-start_intersection<< endl;
+ long global_end =clock();
+ cout<< "Global time= " << global_end - global_start << endl;
}
+
+ return result_areas;
}
-
-
- /*! Readjusts a set of bounding boxes so that they are extended
- in all dimensions for avoiding missing interesting intersections
-
- \param bbox vector containing the bounding boxes
- */
-
- void Interpolation3DSurf::adjustBoundingBoxes(vector<double>& bbox)
- {
- /* We build the segment tree for locating possible matching intersections*/
-
- int size = bbox.size()/6;
- for (int i=0; i<size; i++)
- {
- double Dx=bbox[i*6+1]-bbox[i*6];
- double Dy=bbox[i*6+3]-bbox[i*6+2];
- double Dz=bbox[i*6+5]-bbox[i*6+4];
- double max=(Dx<Dy)?Dy:Dx;
- max=(max<Dz)?Dz:max;
- bbox[i*6]-=_Surf3DAdjustmentEps*max;
- bbox[i*6+1]+=_Surf3DAdjustmentEps*max;
- bbox[i*6+2]-=_Surf3DAdjustmentEps*max;
- bbox[i*6+3]+=_Surf3DAdjustmentEps*max;
- bbox[i*6+4]-=_Surf3DAdjustmentEps*max;
- bbox[i*6+5]+=_Surf3DAdjustmentEps*max;
-
- }
- }
-
};
#ifndef _INTERPOLATION_3D_SURF_HXX_
#define _INTERPOLATION_3D_SURF_HXX_
-
-
#include "MEDMEM_Mesh.hxx"
-#include "InterpolationUtils.hxx"
-#include "TranslationRotationMatrix.hxx"
#include "Interpolation.hxx"
-#include <cmath>
-#include <map>
-#include <sstream>
-#include <string>
-#include <set>
-#include <algorithm>
-#include <vector>
-
namespace MEDMEM
{
-
class Interpolation3DSurf : public Interpolation
{
-
-
private:
+ IntersectionType _Intersection_type;
+ double _MedianPlane;
+ bool _Do_rotate;
double _Precision;
double _DimCaracteristic;
- int _SecondMethodOfLocalisation;
- int _PrintLevel;
double _Surf3DAdjustmentEps;
- // Méthodes publiques
- public:
+ int _PrintLevel;
+ public:
Interpolation3DSurf();
-
-
- // precision geometrique, choix de la methode comlementaire de localisation, niveau d'impressions
- void setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel);
-
-
-
- //fonction principale pour interpoler deux maillages triangulaires.
- std::vector<map<int, double> > interpol_maillages(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
-
+
+ // geometric precision, debug print level, coice of the median plane, intersection etc ...
+ void setOptions(double Precision, int PrintLevel, double MedianPlane,
+ IntersectionType intersectionType, bool do_rotate);
+
+ // Main function to interpolate triangular and qudadratic meshes
+ std::vector<map<int, double> > interpolateMeshes(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
private:
- //fonction qui permet de remplir le vecteur donnant la surface d'intersection
- //de la mailles i_P du maillage projetté avec la maille i_S du maillage support.
- inline void rempli_vect_d_intersect(int i_P,int i_S,
- const INTERP_UTILS::monMaillageP& MaP,
- const INTERP_UTILS::monMaillageS& MaS,
- vector<map<int,double> >& Surface_d_intersect);
-
-
- void createBoundingBoxes(const struct INTERP_UTILS::monMaillageS mesh, vector<double>& bbox);
-
- void rotate3DTriangle(const double* PP1,
- const double* PP2,
- const double* PP3,
- TranslationRotationMatrix& matrix);
-
- void adjustBoundingBoxes(vector<double>& bbox);
-
-
};
};
#include "MEDMEM_Mesh.hxx"
#include "MEDMEM_Field.hxx"
-
#include <cmath>
#include <map>
#include <algorithm>
inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3)
{
double A=(P_3[1]-P_1[1])*(P_2[0]-P_1[0])-(P_2[1]-P_1[1])*(P_3[0]-P_1[0]);
- double Surface = 1/2.0*abs(A);
+ double Surface = 0.5*abs(A);
return Surface;
}
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
inline bool point_dans_triangle(const double* P_0,const double* P_1,
- const double* P_2,const double* P_3)
+ const double* P_2,const double* P_3,
+ double eps)
{
bool A=false;
double det_1=mon_determinant(P_1,P_3,P_0);
double det_2=mon_determinant(P_3,P_2,P_0);
double det_3=mon_determinant(P_2,P_1,P_0);
- if( (det_1>=0 && det_2>=0 && det_3>=0) || (det_1<=0 && det_2<=0 && det_3<=0) )
+ if( (det_1>=-eps && det_2>=-eps && det_3>=-eps) || (det_1<=eps && det_2<=eps && det_3<=eps) )
{
A=true;
}
/* le vecteur et le rajouter au vecteur sinon. */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- inline void verif_point_dans_vect(const double* P, vector<double>& V,double dim_caracteristic, double precision )
+ inline void verif_point_dans_vect(const double* P, vector<double>& V, double absolute_precision )
{
// int taille=V.size();
// for(int i=0;i<taille/2;i++)
bool isPresent=false;
for(int i=0;i<taille/2;i++)
{
- if (sqrt(((P[0]-V[2*i])*(P[0]-V[2*i])+(P[1]-V[2*i+1])*(P[1]-V[2*i+1])))/dim_caracteristic<precision)
+ if (sqrt(((P[0]-V[2*i])*(P[0]-V[2*i])+(P[1]-V[2*i+1])*(P[1]-V[2*i+1])))<absolute_precision)
isPresent=true;
}
const double* P_4,const double* P_5,const double* P_6,vector<double>& V, double dim_caracteristic, double precision)
{
- bool A_1=INTERP_UTILS::point_dans_triangle(P_1,P_4,P_5,P_6);
+ double absolute_precision = precision*dim_caracteristic;
+ bool A_1=INTERP_UTILS::point_dans_triangle(P_1,P_4,P_5,P_6,absolute_precision);
if(A_1)
- {verif_point_dans_vect(P_1,V,dim_caracteristic, precision);
+ {verif_point_dans_vect(P_1,V,absolute_precision);
// if (V.size()==1)
// {
// // all nodes are necessarily in the triangle
// return;
// }
}
- bool A_2=INTERP_UTILS::point_dans_triangle(P_2,P_4,P_5,P_6);
+ bool A_2=INTERP_UTILS::point_dans_triangle(P_2,P_4,P_5,P_6,absolute_precision);
if(A_2)
- {verif_point_dans_vect(P_2,V,dim_caracteristic,precision);}
+ {verif_point_dans_vect(P_2,V,absolute_precision);}
- bool A_3=INTERP_UTILS::point_dans_triangle(P_3,P_4,P_5,P_6);
+ bool A_3=INTERP_UTILS::point_dans_triangle(P_3,P_4,P_5,P_6,absolute_precision);
if(A_3)
- {verif_point_dans_vect(P_3,V,dim_caracteristic, precision);}
+ {verif_point_dans_vect(P_3,V,absolute_precision);}
}
{
- // calcul du déterminant de P1_1P_2 et P_3P_4.
+ // calcul du déterminant de P_1P_2 et P_3P_4.
double det=(P_2[0]-P_1[0])*(P_4[1]-P_3[1])-(P_4[0]-P_3[0])*(P_2[1]-P_1[1]);
-
- if(abs(det)/dim_caracteristic>precision)
+ double absolute_precision = dim_caracteristic*precision;
+ if(abs(det)>absolute_precision)
{
- double k_1=1/((P_2[0]-P_1[0])*(P_3[1]-P_4[1])+(P_4[0]-P_3[0])*(P_2[1]-P_1[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]));
+ 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)
{
- double k_2=1/((P_4[0]-P_3[0])*(P_1[1]-P_2[1])+(P_2[0]-P_1[0])*(P_4[1]-P_3[1]))
- *((P_1[1]-P_2[1])*(P_1[0]-P_3[0])+(P_2[0]-P_1[0])*(P_1[1]-P_3[1]));
+ 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)
double P_0[2];
P_0[0]=P_1[0]+k_1*(P_2[0]-P_1[0]);
P_0[1]=P_1[1]+k_1*(P_2[1]-P_1[1]);
- verif_point_dans_vect(P_0,Vect, dim_caracteristic, precision);
+ verif_point_dans_vect(P_0,Vect,absolute_precision);
}
}
/* P_4, P_5, P_6: sommets du deuxième triangle */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
- inline vector<double> intersec_de_triangle(const double* P_1,const double* P_2,
- const double* P_3,const double* P_4,const double* P_5,const double* P_6, double dim_caracteristic, double precision)
+ inline void intersec_de_triangle(const double* P_1,const double* P_2, const double* P_3,
+ const double* P_4,const double* P_5,const double* P_6,
+ vector<double>& Vect, double dim_caracteristic, double precision)
{
// cout << " T1 : " << *P_1 << " , " << *(P_1+1) << " , " << *P_2 << " , " << *(P_2+1) << " , " << *P_3 << " , " << *(P_3+1)<< endl;
// cout << " T2 : " << *P_4 << " , " << *(P_4+1) << " , " << *P_5 << " , " << *(P_5+1) << " , " << *P_6 << " , " << *(P_6+1)<< endl;
- vector<double> Vect;
+ //MN: vector<double> Vect;
inters_de_segment(P_1,P_2,P_4,P_5,Vect, dim_caracteristic, precision);
inters_de_segment(P_1,P_2,P_5,P_6,Vect, dim_caracteristic, precision);
//debug cout << Vect[i] << " ";
//debug cout << endl << endl;
- return Vect;
+ //MN: return Vect;
}
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ */
FractSurf.write(id3);
}
- inline void getElemBB(double* bb, const monMaillageP& MaP, int iP)
+ inline void getElemBB(double* bb, const MEDMEM::MESH& mesh, int iP, int nb_nodes)
{
bb[0]=HUGE;bb[1]=-HUGE;
bb[2]=HUGE;bb[3]=-HUGE;
bb[4]=HUGE;bb[5]=-HUGE;
- for (int i=0; i<3; i++)
+ for (int i=0; i<nb_nodes; i++)
{
- double x = MaP._Coord[(MaP._Connect[3*iP+i]-1)*3];
- double y = MaP._Coord[(MaP._Connect[3*iP+i]-1)*3+1];
- double z = MaP._Coord[(MaP._Connect[3*iP+i]-1)*3+2];
+ double x = *(mesh.getCoordinates(MED_FULL_INTERLACE)+3*(iP+i));
+ double y = *(mesh.getCoordinates(MED_FULL_INTERLACE)+3*(iP+i)+1);
+ double z = *(mesh.getCoordinates(MED_FULL_INTERLACE)+3*(iP+i)+2);
bb[0]=(x<bb[0])?x:bb[0];
bb[1]=(x>bb[1])?x:bb[1];
bb[2]=(y<bb[2])?y:bb[2];
bb[3]=(y>bb[3])?y:bb[3];
bb[4]=(z<bb[4])?z:bb[4];
bb[5]=(z>bb[5])?z:bb[5];
- }
-
+ }
}
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /* Computes the dot product of a and 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];
+ return result;
+ }
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /* Computes the norm of vector 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 */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ 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;
+ }
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /* Computes the determinant of a and b */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ inline double determinant ( double * a, double * b)
+ {
+ return a[0]*b[1]-a[1]*b[0];
+ }
+ inline double determinant ( double * a, double * b, double * c)
+ {
+ return a[0]*determinant(b+1,c+1)-b[0]*determinant(a+1,c+1)+c[0]*determinant(a+1,b+1);
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /* Computes the cross product of AB and AC */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-
+ template<int dim> inline void crossprod( const double * A, const double * B, const double * C, double * V);
+
+ template<> inline
+ void crossprod<2>( const double * A, const double * B, const double * C, double * V)
+ {
+ double AB[2];
+ double AC[2];
+ for(int idim =0; idim<2; idim++) AB[idim] = B[idim]-A[idim];//B-A
+ for(int idim =0; idim<2; idim++) AC[idim] = C[idim]-A[idim];//C-A;
+
+ V[0]=determinant(AB,AC);
+ V[1]=0;
+ }
+ template<> inline
+ void crossprod<3>( const double * A, const double * B, const double * C, double * V)
+ {
+ double AB[3];
+ double AC[3];
+ for(int idim =0; idim<3; idim++) AB[idim] = B[idim]-A[idim];//B-A
+ for(int idim =0; idim<3; idim++) AC[idim] = C[idim]-A[idim];//C-A;
+
+ V[0]=AB[1]*AC[2]-AB[2]*AC[1];
+ V[1]=-AB[0]*AC[2]+AB[2]*AC[0];
+ V[2]=AB[0]*AC[1]-AB[1]*AC[0];
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /* Checks wether point A is inside the quadrangle BCDE */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+
+ template<int dim> inline double check_inside(const double* A,const double* B,const double* C,const double* D,
+ const double* E,double* ABC, double* ADE)
+ {
+ crossprod<dim>(A,B,C,ABC);
+ crossprod<dim>(A,D,E,ADE);
+ return dotprod<dim>(ABC,ADE);
+ }
-};
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /* Computes the geometric angle (in [0,Pi]) between two non zero vectors AB and AC */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ template<int dim> inline double angle(const double * A, const double * B, const double * C, double * n)
+ {
+ double AB[dim];
+ double AC[dim];
+ double orthAB[dim];
+
+ for(int idim =0; idim<dim; idim++) AB[idim] = B[idim]-A[idim];//B-A;
+ for(int idim =0; idim<dim; idim++) AC[idim] = C[idim]-A[idim];//C-A;
+
+ double normAB= norm<dim>(AB);
+ for(int idim =0; idim<dim; idim++) AB[idim]/=normAB;
+
+ double normAC= norm<dim>(AC);
+ double AB_dot_AC=dotprod<dim>(AB,AC);
+ for(int idim =0; idim<dim; idim++) orthAB[idim] = AC[idim]-AB_dot_AC*AB[idim];
+
+ double denom= normAC+AB_dot_AC;
+ double numer=norm<dim>(orthAB);
+
+ return 2*atan2(numer,denom);
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /* Tells whether the frame constituted of vectors AB, AC and n is direct */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ template<int dim> inline double direct_frame(const double * A, const double * B, const double * C, double * n);
+ template<> inline
+ double direct_frame<2>(const double * A, const double * B, const double * C, double * n)
+ {
+ double AB[2];
+ double AC[2];
+ for(int idim =0; idim<2; idim++) AB[idim] = B[idim]-A[idim];//B-A;
+ for(int idim =0; idim<2; idim++) AC[idim] = C[idim]-A[idim];//C-A;
+
+ return determinant(AB,AC)*n[0];
+ }
+ template<> inline
+ double direct_frame<3>(const double * A, const double * B, const double * C, double * n)
+ {
+ double AB[3];
+ double AC[3];
+ for(int idim =0; idim<3; idim++) AB[idim] = B[idim]-A[idim];//B-A;
+ for(int idim =0; idim<3; idim++) AC[idim] = C[idim]-A[idim];//C-A;
+
+ return determinant(AB,AC,n)>0;
+ }
+};
#endif
TetraAffineTransform.hxx\
TranslationRotationMatrix.hxx\
Interpolation2D.hxx\
+Interpolation2Dbis.hxx\
Interpolation3DSurf.hxx\
InterpolationUtils.hxx\
BBTree.H\
TransformedTriangle_inline.hxx\
IntersectorTetra.hxx\
TargetIntersector.hxx\
-IntersectorHexa.hxx
+IntersectorHexa.hxx\
+PolygonAlgorithms.hxx\
+PolygonAlgorithms.cxx\
+PlanarIntersector.hxx\
+TriangulationIntersector.hxx\
+TriangulationIntersector.cxx\
+ConvexIntersector.hxx\
+ConvexIntersector.cxx\
+GenericIntersector.hxx\
+GenericIntersector.cxx
# Libraries targets
Interpolation3D.cxx\
Interpolation3DSurf.cxx\
Interpolation2D.cxx\
+Interpolation2Dbis.cxx\
IntersectorTetra.cxx\
-IntersectorHexa.cxx
+IntersectorHexa.cxx\
+PlanarIntersector.cxx
# Executables targets
BIN =
+# test_InterpolationUtils
BIN_SRC =
BIN_SERVER_IDL =
BIN_CLIENT_IDL =
-TEST_PROGS = testUnitTetra
+TEST_PROGS = testUnitTetra test_PolygonAlgorithms test_Interpol_2D test_Interpol_2Dbis test_Interpol_3DSurf
LDFLAGS+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome
LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome
-CXXFLAGS+=@CXXTMPDPTHFLAGS@ $(MED2_INCLUDES)
+CXXFLAGS+=@CXXTMPDPTHFLAGS@ $(MED2_INCLUDES) -I$(srcdir)
CPPFLAGS+=$(BOOST_CPPFLAGS)
# optimization
-CXXFLAGS+= -O2 -DOPTIMIZE -Wall
-CPPFLAGS+= -O2 -DOPTIMIZE -Wall
+CXXFLAGS+= -DOPTIMIZE -Wall
+CPPFLAGS+= -DOPTIMIZE -Wall
# for log
CXXFLAGS+= -DLOG_LEVEL=0
LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 $(STDLIB) -lmedmem -lutil
#LDFLAGS+= $(HDF5_LIBS) $(STDLIB) -lutil
# for gcov
-#LDFLAGS+= -lgcov
-#LDFLAGS+= -pg
+LDFLAGS+= -lgcov
+LDFLAGS+= -pg
#LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS)
# change motivated by the bug KERNEL4778.
--- /dev/null
+#include "PlanarIntersector.hxx"
+#include "MEDMEM_Mesh.hxx"
+
+namespace MEDMEM {
+/*!
+ \brief creates the bounding boxes for all the cells of mesh \a mesh
+
+ The method accepts mixed meshes (containing triangles and quadrangles).
+ The vector returned is of dimension 6*nb_elems with bounding boxes stored as xmin1, xmax1, ymin1, ymax1, zmin1, zmax1, xmin2, xmax2, ymin2,...
+ The returned pointer must be deleted by the calling code.
+
+ \param mesh structure pointing to the mesh
+ \param bbox vector containing the bounding boxes
+ */
+ template <int DIM>
+ void PlanarIntersector::createBoundingBoxes(const MEDMEM::MESH& mesh, vector<double>& bbox)
+ {
+ /* We build the segment tree for locating possible matching intersections*/
+ int nbelems = mesh.getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
+ int nbtypes = mesh.getNumberOfTypes(MED_EN::MED_CELL);
+ const MED_EN::medGeometryElement* types = mesh.getTypes(MED_EN::MED_CELL);
+ bbox.resize(2*DIM* nbelems);
+ const double* coords = mesh.getCoordinates(MED_EN::MED_FULL_INTERLACE);
+ const int* conn = mesh.getConnectivity(MED_EN::MED_FULL_INTERLACE,
+ MED_EN::MED_NODAL,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ const int* conn_index = mesh.getConnectivityIndex(MED_EN::MED_NODAL, MED_EN::MED_CELL);
+ int ibox=0;
+ int type_max_index;//maximum cell number for a given type
+ int type_min_index=1;//minimum cell number for a given type
+ int icell;
+ for (int itype = 0; itype<nbtypes; itype++)
+ {
+ int nbelem_type = mesh.getNumberOfElements(MED_EN::MED_CELL, types[itype]);
+ int nb_nodes_per_elem = types[itype]%100;
+ type_max_index = type_min_index + nbelem_type;
+
+ for (icell=type_min_index; icell<type_max_index; icell++)
+ {
+ //initializing bounding box limits
+ for(int idim=0; idim<DIM; idim++)
+ {
+ bbox[2*DIM*ibox+2*idim] = HUGE;
+ bbox[2*DIM*ibox+2*idim+1] = -HUGE;
+ }
+ //updating the bounding box with each node of the element
+ for (int j=0; j<nb_nodes_per_elem; j++)
+ {
+ const double* coord_node=coords + DIM*(conn[conn_index[icell-1]-1+j]-1);
+ for(int idim=0; idim<DIM; idim++)
+ {
+ double x=*(coord_node+idim);
+ bbox[ibox*2*DIM + 2*idim] = (bbox[ibox*2*DIM + 2*idim] <x)?bbox[ibox*2*DIM + 2*idim ]:x;
+ bbox[ibox*2*DIM + 2*idim+1] = (bbox[ibox*2*DIM + 2*idim+1]>x)?bbox[ibox*2*DIM + 2*idim+1]:x;
+ }
+ }
+ ibox++;
+ }
+ type_min_index = type_max_index;
+ }
+ }
+ /*
+ Computes the bouding box of a given element
+ */
+ template <int DIM>
+ inline void PlanarIntersector::getElemBB(double* bb, const MEDMEM::MESH& mesh, int iP, int nb_nodes)
+ {
+ const double* coords = mesh.getCoordinates(MED_EN::MED_FULL_INTERLACE);
+ const int* conn_index = mesh.getConnectivityIndex(MED_EN::MED_NODAL, MED_EN::MED_CELL);
+ const int* conn = mesh.getConnectivity(MED_EN::MED_FULL_INTERLACE,
+ MED_EN::MED_NODAL,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ //initializing bounding box limits
+ for(int idim=0; idim<DIM; idim++)
+ {
+ bb[2*idim ] = HUGE;
+ bb[2*idim+1] = -HUGE;
+ }
+
+ for (int i=0; i<nb_nodes; i++)
+ {
+ //MN: iP= cell index, not node index, use of connectivity array ?
+ const double* coord_node=coords + DIM*(conn[conn_index[iP-1]-1+i]-1);
+ for(int idim=0; idim<DIM; idim++)
+ {
+ double x = *(coord_node+idim);
+ //double y = *(mesh.getCoordinates(MED_FULL_INTERLACE)+3*(iP+i)+1);
+ bb[2*idim ] = (x<bb[2*idim ])?x:bb[2*idim ];
+ bb[2*idim+1] = (x>bb[2*idim+1])?x:bb[2*idim+1];
+ }
+ }
+ }
+
+ /*! Readjusts a set of bounding boxes so that they are extended
+ in all dimensions for avoiding missing interesting intersections
+
+ \param bbox vector containing the bounding boxes
+ */
+
+ template <int DIM>
+ void PlanarIntersector::adjustBoundingBoxes(vector<double>& bbox, double Surf3DAdjustmentEps)
+ {
+ /* We build the segment tree for locating possible matching intersections*/
+
+ int size = bbox.size()/(2*DIM);
+ for (int i=0; i<size; i++)
+ {
+ double max=-HUGE;
+ for(int idim=0; idim<DIM; idim++)
+ {
+ double Dx=bbox[i*2*DIM+1+2*idim]-bbox[i*2*DIM+2*idim];
+ max=(max<Dx)?Dx:max;
+ }
+ for(int idim=0; idim<DIM; idim++)
+ {
+ bbox[i*2*DIM+2*idim ] -= Surf3DAdjustmentEps*max;
+ bbox[i*2*DIM+2*idim+1] += Surf3DAdjustmentEps*max;
+ }
+ }
+ }
+}
--- /dev/null
+#include "TranslationRotationMatrix.hxx"
+#include "InterpolationUtils.hxx"
+#include "PlanarIntersector.hxx"
+#include "MEDMEM_Mesh.hxx"
+
+using namespace MED_EN;
+using namespace INTERP_UTILS;
+namespace MEDMEM
+{
+
+ void PlanarIntersector::projection( vector< double>& Coords_A, vector< double>& Coords_B,
+ int nb_NodesA, int nb_NodesB, double epsilon,
+ double median_plane, bool do_rotate)
+ {
+ double normal_A[3]={0,0,0};
+ double normal_B[3]={0,0,0};
+ double linear_comb[3];
+ double proj;
+ const int DIM=3;
+
+ //Find the normal to cells A and B
+ int i_A1=1;
+ while(i_A1<nb_NodesA && distance2<DIM>(&Coords_A[0],&Coords_A[DIM*i_A1])< epsilon) i_A1++;
+ int i_A2=i_A1+1;
+ crossprod<DIM>(&Coords_A[0], &Coords_A[DIM*i_A1], &Coords_A[DIM*i_A2],normal_A);
+ while(i_A2<nb_NodesA && dotprod<DIM>(normal_A,normal_A)<epsilon)
+ {
+ crossprod<DIM>(&Coords_A[0], &Coords_A[DIM*i_A1], &Coords_A[DIM*i_A2],normal_A);
+ i_A2++;
+ }
+ int i_B1=1;
+ while(i_B1<nb_NodesB && distance2<DIM>(&Coords_B[0],&Coords_B[DIM*i_B1])< epsilon) i_B1++;
+ int i_B2=i_B1+1;
+ crossprod<DIM>(&Coords_B[0], &Coords_B[DIM*i_B1], &Coords_B[DIM*i_B2],normal_B);
+ while(i_B2<nb_NodesB && dotprod<DIM>(normal_B,normal_B)< epsilon)
+ {
+ crossprod<DIM>(&Coords_B[0], &Coords_B[DIM*i_B1], &Coords_B[DIM*i_B2],normal_B);
+ i_B2++;
+ }
+
+ if(i_A2<nb_NodesA && i_B2<nb_NodesB)
+ {
+ //Build the normal of the median plane
+ if(dotprod<DIM>(normal_A,normal_B)<0)
+ for(int idim =0; idim< DIM; idim++) normal_A[idim] *=-1;
+ for(int idim =0; idim< DIM; idim++)
+ linear_comb[idim] = median_plane*normal_A[idim] + (1-median_plane)*normal_B[idim];
+ double norm= sqrt(dotprod<DIM>(linear_comb,linear_comb));
+
+ if(norm>epsilon)
+ {
+ for(int idim =0; idim< DIM; idim++) linear_comb[idim]/=norm;
+
+ //Project the nodes of A and B on the median plane
+ for(int i_A=0; i_A<nb_NodesA; i_A++)
+ {
+ proj = dotprod<DIM>(&Coords_A[DIM*i_A],linear_comb);
+ for(int idim =0; idim< DIM; idim++)
+ Coords_A[DIM*i_A+idim] -= proj*linear_comb[idim];
+ }
+ for(int i_B=0; i_B<nb_NodesB; i_B++)
+ {
+ proj = dotprod<DIM>(&Coords_B[DIM*i_B],linear_comb);
+ for(int idim =0; idim< DIM; idim++)
+ Coords_B[DIM*i_B+idim] -= proj*linear_comb[idim];
+ }
+
+ //Buid the matrix sending A into the Oxy plane and apply it to A and B
+ if(do_rotate)
+ {
+ TranslationRotationMatrix rotation;
+ //rotate3DTriangle(&Coords_A[0], &Coords_A[DIM*i_A1], &Coords_A[DIM*i_A2], rotation);
+ rotate3DTriangle(&Coords_B[0], &Coords_B[DIM*i_B1], &Coords_B[DIM*i_B2], rotation);
+ for (int i=0; i<nb_NodesA; i++) rotation.transform_vector(&Coords_A[DIM*i]);
+ for (int i=0; i<nb_NodesB; i++) rotation.transform_vector(&Coords_B[DIM*i]);
+ }
+ }
+ }
+ else
+ {
+ cout<< " Maille dégénérée " << "epsilon = " << epsilon << endl;
+ cout<< " i_A1= " << i_A1 << " i_A2= " << i_A2 << endl;
+ cout<< " distance2<DIM>(&Coords_A[0],&Coords_A[i_A1])= " << distance2<DIM>(&Coords_A[0],&Coords_A[i_A1]) << endl;
+ cout<< "abs(normal_A) = " << abs(normal_A[0]) << " ; " <<abs( normal_A[1]) << " ; " << abs(normal_A[2]) << endl;
+ cout<< " i_B1= " << i_B1 << " i_B2= " << i_B2 << endl;
+ cout<< " distance2<DIM>(&Coords_B[0],&Coords_B[i_B1])= " << distance2<DIM>(&Coords_B[0],&Coords_B[i_B1]) << endl;
+ cout<< "normal_B = " << normal_B[0] << " ; " << normal_B[1] << " ; " << normal_B[2] << endl;
+ }
+ }
+
+ void PlanarIntersector::rotate3DTriangle( double* PP1, double*PP2, double*PP3,
+ TranslationRotationMatrix& rotation_matrix)
+ {
+ //initializes
+ rotation_matrix.translate(PP1);
+
+ double P2w[3];
+ double P3w[3];
+ P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2];
+ P3w[0]=PP3[0]; P3w[1]=PP3[1];P3w[2]=PP3[2];
+
+ // translating to set P1 at the origin
+ for (int i=0; i<2; i++)
+ {
+ P2w[i]-=PP1[i];
+ P3w[i]-=PP1[i];
+ }
+
+ // rotating to set P2 on the Oxy plane
+ TranslationRotationMatrix A;
+ A.rotate_x(P2w);
+ A.rotate_vector(P3w);
+ rotation_matrix.multiply(A);
+
+ //rotating to set P2 on the Ox axis
+ TranslationRotationMatrix B;
+ B.rotate_z(P2w);
+ B.rotate_vector(P3w);
+ rotation_matrix.multiply(B);
+
+ //rotating to set P3 on the Oxy plane
+ TranslationRotationMatrix C;
+ C.rotate_x(P3w);
+ rotation_matrix.multiply(C);
+ }
+
+}
--- /dev/null
+#ifndef _PLANAR_INTERSECTOR_HXX_
+#define _PLANAR_INTERSECTOR_HXX_
+
+#include "TranslationRotationMatrix.hxx"
+
+namespace MEDMEM
+{
+ class MESH;
+ class PlanarIntersector
+ {
+ public:
+ PlanarIntersector (const MEDMEM::MESH& mesh_A, const MEDMEM::MESH& mesh_B,
+ double DimCaracteristic, double Precision, double MedianPlane){};
+ virtual ~ PlanarIntersector(){};
+
+ //Tool for cell intersection
+ virtual double intersectCells(int icell_A, int icell_B, int nb_NodesA, int nb_NodesB)=0;
+
+ //Tool for cell filtering
+ template <int DIM>
+ void createBoundingBoxes(const MEDMEM::MESH& mesh, vector<double>& bbox);
+ template <int DIM>
+ void adjustBoundingBoxes(vector<double>& bbox, double Surf3DAdjustmentEps );
+ template <int DIM>
+ void getElemBB(double* bb, const MEDMEM::MESH& mesh, int iP, int nb_nodes);
+
+ protected :
+ void projection( vector< double>& Coords_A, vector< double>& Coords_B,
+ int nb_NodesA, int nb_NodesB, double epsilon, double median_plane, bool do_rotate);
+ void rotate3DTriangle( double* PP1, double*PP2, double*PP3,
+ TranslationRotationMatrix& rotation_matrix);
+
+ };
+}
+
+#endif
--- /dev/null
+#ifndef _POLYGON_ALGORITHMS_CXX_
+#define _POLYGON_ALGORITHMS_CXX_
+
+#include "PolygonAlgorithms.hxx"
+#include "InterpolationUtils.hxx"
+#include "MEDMEM_Exception.hxx"
+
+#include <deque>
+
+namespace INTERP_UTILS
+{
+ template< int DIM>
+ PolygonAlgorithms<DIM>::PolygonAlgorithms(double epsilon, double precision)//: (0)
+ {
+ _Inter_size = 0;//To do remove des push_back and update of _Inter_size
+ _Is_in_intersection = false;
+ _Epsilon = epsilon;
+ _Precision = precision;
+ }
+ /*************************************************************/
+ /* Computes the 3D intersection between two COPLANAR */
+ /* Segments [A,B] and [C,D], stores the result in V. */
+ /* If A belongs to [CD] then the vertex E (preceeding A) */
+ /* is used to decide if the crossing is real. If A coincides */
+ /* with C or D, a special treatment is performed */
+ /*************************************************************/
+ template< int DIM>
+ bool PolygonAlgorithms<DIM>::intersect_segment_segment(const double * A, const double * B, const double * C,
+ const double * D, const double * E, double * V)
+ {
+ double AB[DIM], DC[DIM], AC[DIM];
+ for(int idim=0;idim<DIM;idim++)
+ {
+ AB[idim] = B[idim]-A[idim];//B-A
+ DC[idim] = C[idim]-D[idim];//C-D
+ AC[idim] = C[idim]-A[idim];//C-A
+ }
+ double det = determinant(AB,DC);//determinant of the first two coefficients
+ double t1, t2;
+ int z_plane=0;
+ if(abs(det) < _Epsilon)
+ {
+ switch(DIM)
+ {
+ case 2:
+ return false;
+ case 3:
+ det = determinant(&AB[1],&DC[1]);//determinant of the two last coefficients
+ if(abs(det) < _Epsilon) return false;//case of paralell segments
+ else z_plane=1;
+ }
+ }
+
+ t1 = determinant(&AC[z_plane],&DC[z_plane])/det;//solves the linear system t1*AB+t2*DC=AC
+ t2 = determinant(&AB[z_plane],&AC[z_plane])/det;//solves the linear system t1*AB+t2*DC=AC
+
+ if(t1>_Precision && t1<1-_Precision)
+ {
+ if( t2>_Precision && t2<1-_Precision)
+ {
+ for(int idim=0;idim<DIM;idim++) V[idim]=t1*AB[idim] + A[idim];
+ return true;
+ }
+ }
+ else if(abs(t1) <= _Precision)
+ {
+ if( t2>_Precision && t2<1-_Precision)//vertex on an edge
+ {
+ double V12[DIM];
+ double V34[DIM];
+ crossprod<DIM>(A,D,B,V12);
+ crossprod<DIM>(A,D,E,V34);
+ double same_side =dotprod<DIM>(V12, V34);
+ if( same_side < -_Epsilon ) // <= epsilon or 0 ?//crossing
+ {
+ for(int idim=0;idim<DIM;idim++) V[idim]=A[idim];
+ return true;
+ }
+ else if( same_side > _Epsilon ) _Terminus= !_Is_in_intersection;//reflexion
+ else //separation of overlaping edges
+ {
+ if(_Inter.size() ==0) _Terminus=true;
+ else if(!_Is_in_intersection)
+ {
+ for(int idim=0;idim<DIM;idim++) V[idim]=A[idim];
+ return true;
+ }
+ }
+ }
+ else if(abs(t2-1) <= _Precision)//vertex on a vertex (A=D), first run
+ crossprod<DIM>(A,C,E,_Vdouble);//angle between vectors AC and AE (E=vertex preceding A)
+ else if(abs(t2) <= _Precision)//vertex on a vertex (A=C), second run
+ {
+ double Vdoublebis[DIM];
+ crossprod<DIM>(A,B,D,Vdoublebis);
+ double in_between =dotprod<DIM>(Vdoublebis,_Vdouble);
+ if(in_between>_Epsilon)//crossing
+ {
+ for(int idim=0;idim<DIM;idim++) V[idim]=A[idim];
+ return true;
+ }
+ else if(abs(in_between)<=_Epsilon && dotprod<DIM>(Vdoublebis,Vdoublebis) > _Epsilon)
+ //ie _Vdouble=0, separation of overlaping edges at a double point
+ {
+ crossprod<DIM>(A,E,B,_Vdouble);
+ if(dotprod<DIM>(_Vdouble,Vdoublebis) >=_Epsilon )//crossing
+ {
+ if(_Inter.size()==0) _Terminus=true;
+ else if(!_Is_in_intersection)
+ {
+ for(int idim=0;idim<DIM;idim++) V[idim]=A[idim];
+ return true;
+ }
+ }
+ }
+ }
+ }
+ return false;
+ }
+
+ /*************************************************************/
+ /* adds vertex i to the list inter and updates _End_segments */
+ /* i is the local index of the current vertex */
+ /*************************************************************/
+ template< int DIM>
+ inline void PolygonAlgorithms<DIM>::add_new_vertex( int i, int i_glob, int i_next_glob, int i_prev_glob,
+ const double * P)
+ {
+ /* Question:Should we add vertex i to the front or back ? */
+ if( _End_segments[1].second == i_glob)
+ {
+ for(int idim=0;idim<DIM;idim++) _Inter.push_back(P[DIM*i+idim]);
+ _End_segments[1] = make_pair(i_glob, i_next_glob);
+ }
+ else
+ {
+ for(int idim=DIM-1;idim>-1;idim--) _Inter.push_front(P[DIM*i+idim]);
+ _End_segments[0] = make_pair(i_glob, i_next_glob);
+ }
+ }
+
+ /************************************************************/
+ /* adds a crossing between two segments starting at i and j */
+ /* to the double ended list inter in the correct order */
+ /* according to endsegments, updates _End_segments */
+ /************************************************************/
+ template< int DIM>
+ inline void PolygonAlgorithms<DIM>::add_crossing( double * ABCD, pair< int,int > i_i_next,
+ pair< int,int > j_j_next)
+ {
+ if( _Inter.size()>0 )
+ {
+ if(_End_segments[1] ==i_i_next)
+ {
+ for(int idim=0;idim<DIM;idim++) _Inter.push_back(ABCD[idim]);
+ _Terminus= (_End_segments[0]== j_j_next);
+ _End_segments[1] = j_j_next;
+ }
+ else
+ {
+ if( _End_segments[1]== j_j_next)
+ {
+ for(int idim=0;idim<DIM;idim++) _Inter.push_back(ABCD[idim]);
+ _Terminus= (_End_segments[0]== i_i_next);
+ _End_segments[1] = i_i_next;
+ }
+ else
+ {
+ for(int idim=DIM-1;idim>-1;idim--) _Inter.push_front(ABCD[idim]);
+ _End_segments[0] = (_End_segments[0]== i_i_next) ? j_j_next : i_i_next;
+ }
+ }
+ }
+ else
+ {
+ for(int i=0;i<DIM;i++) _Inter.push_back(ABCD[i]);
+ _End_segments.push_back(i_i_next);
+ _End_segments.push_back(j_j_next);
+ }
+ }
+
+ /*******************************************************/
+ /* checks the possible crossing between segments [A,B] */
+ /* (with end-point global indices i and i_next) */
+ /* and [C,D] (end-point global indices j and j_next). */
+ /* If no intersection is detected, checks whether B is */
+ /* inside the quadrangle AEDC. */
+ /* Updates the lists inter and _End_segments */
+ /*******************************************************/
+
+ template< int DIM>
+ void PolygonAlgorithms<DIM>::add_crossing0(const double * A, const double * B, int i, int i_next,
+ const double * C, const double * D, int j, int j_next)
+ {
+ double ABCD[DIM];
+ if(intersect_segment_segment(A,B,C,D,ABCD, ABCD))
+ //fifth and sixth arguments are useless here
+ {
+ add_crossing(ABCD, make_pair(i, i_next), make_pair(j, j_next));
+ _Status.insert(make_pair(i_next,make_pair(i, false)));
+ multimap< int, pair< int,bool> >::iterator mi;
+ mi=_Status.find(j_next);
+ ((* mi).second).second= !((* mi).second).second;
+ }
+ else _Status.insert(make_pair(i_next,make_pair(i,true)));
+ }
+
+ /*******************************************************/
+ /* adds the possible crossings between segments [A,B] (with end-point global indices i and i_next) */
+ /*and segments [C,D] and [E,F] to the list inter and updates _End_segments */
+ /* In cases of ambiguity, the vertex G is used to decide wether the crossing should be accepted */
+ /*******************************************************/
+ template< int DIM>
+ inline void PolygonAlgorithms<DIM>::add_crossings( const double * A, const double * B, int i , int i_next,
+ const double * C, const double * D, int j1, int j2,
+ const double * E, const double * F, int j3, int j4,
+ const double * G)
+ {
+ double ABCD[DIM];
+ double ABEF[DIM];
+ multimap< int, pair< int,bool> >::iterator mi;
+
+ if(intersect_segment_segment(A,B,C,D,G,ABCD))
+ {
+ if(intersect_segment_segment(A,B,E,F,G,ABEF))
+ {
+ VertexLess<DIM> vl;
+ if (vl(ABCD,ABEF))
+ {
+ add_crossing(ABCD, make_pair(i, i_next), make_pair(j1, j2));
+ add_crossing(ABEF, make_pair(i, i_next), make_pair(j3, j4));
+ }
+ else
+ {
+ add_crossing(ABEF, make_pair(i, i_next), make_pair(j3, j4));
+ add_crossing(ABCD, make_pair(i, i_next), make_pair(j1, j2));
+ }
+ _Status.insert(make_pair(i_next,make_pair(i, _Is_in_intersection)));
+ mi=_Status.find(j2);
+ ((* mi).second).second= !((* mi).second).second;
+ mi=_Status.find(j4);
+ ((* mi).second).second= !((* mi).second).second;
+ }
+ else
+ {
+ add_crossing(ABCD, make_pair( i, i_next), make_pair(j1,j2));
+ _Status.insert(make_pair(i_next,make_pair(i, !_Is_in_intersection)));
+ mi=_Status.find(j2);
+ ((* mi).second).second= !((* mi).second).second;
+ }
+ }
+ else
+ {
+ if(intersect_segment_segment(A,B,E,F,G, ABEF))
+ {
+ add_crossing(ABEF, make_pair( i, i_next), make_pair( j3, j4));
+ _Status.insert(make_pair(i_next,make_pair(i, !_Is_in_intersection)));
+ mi=_Status.find(j4);
+ ((* mi).second).second= !((* mi).second).second;
+ }
+ else _Status.insert(make_pair(i_next,make_pair(i, _Is_in_intersection)));
+ }
+ }
+
+
+ /* define various indices required in the function intersect_conv_polygon */
+ /* vertices from the both polygons are supposed to be present in the status */
+ template< int DIM>
+ inline void PolygonAlgorithms<DIM>::define_indices(int& i_loc, int& i_next, int& i_prev,
+ const double *& Poly1, const double *& Poly2,
+ int& j1, int& j1_glob, int& j2, int& j2_glob,
+ int& j3, int& j3_glob, int& j4, int& j4_glob,
+ int& i_glob, int& i_next_glob, int& i_prev_glob,
+ const double * P_1, const double * P_2,
+ int N1, int N2, int sign)
+ {
+ int N0, shift;
+ if(i_glob < N1)
+ {
+ N0 = N1;
+ shift = 0;
+ Poly1 = P_1;
+ Poly2 = P_2;
+
+ multimap< int, pair< int,bool> >::reverse_iterator mi1=_Status.rbegin();
+ j1_glob=((*mi1).second).first;
+ j1=j1_glob-N1;
+ j2_glob=(*mi1).first;
+ j2=j2_glob-N1;
+ mi1++;
+ j3_glob=((*mi1).second).first;
+ j3=j3_glob-N1;
+ j4_glob=(*mi1).first;
+ j4=j4_glob-N1;
+ }
+ else
+ {
+ N0 = N2;
+ shift = N1;
+ Poly1 = P_2;
+ Poly2 = P_1;
+
+ multimap< int, pair< int,bool> >::iterator mi2= _Status.begin();
+ j1_glob=((*mi2).second).first;
+ j1=j1_glob;
+ j2_glob=(*mi2).first;
+ j2=j2_glob;
+ mi2++;
+ j3_glob=((*mi2).second).first;
+ j3=j3_glob;
+ j4_glob=(*mi2).first;
+ j4=j4_glob;
+ }
+ i_loc = i_glob-shift;
+ i_next = (i_next_glob-shift+N0)%N0;//end-point of segment starting at i
+ i_prev = (i_prev_glob-shift+N0)%N0;
+ i_next_glob = i_next+shift;
+ i_prev_glob = i_prev+shift;
+ //warning: sign is either 1 or -1;
+ //To do test and remove from Convex_intersecor.cxx
+ while(distance2<DIM>(&Poly1[DIM*i_loc],&Poly1[DIM*i_next])< _Epsilon && i_next != i_loc)
+ i_next =(i_next+sign+N0)%N0;
+ while(distance2<DIM>(&Poly1[DIM*i_loc],&Poly1[DIM*i_prev])< _Epsilon && i_prev != i_loc)
+ i_prev =(i_prev+sign+N0)%N0;
+ }
+ /*******************************************************/
+ /* computes the vertices of the intersection of two COPLANAR */
+ /* simple (no dble points)convex polygons using line sweep algorithm */
+ /* P1 and P2 contain the 3D coordinates of the successive vertices */
+ /*******************************************************/
+ template< int DIM>
+ deque< double > PolygonAlgorithms<DIM>::intersect_convex_polygons(const double* P_1,const double* P_2,
+ int N1, int N2)
+ {
+ int i_loc, i_glob, j1, j1_glob, j2,j2_glob, j3, j3_glob, j4,j4_glob,
+ i_prev, i_prev_glob, i_next, i_next_glob, nb_prev, sign, idim;
+ const double * Poly1, * Poly2;
+ bool four_neighbours=false;
+
+ int minN1N2 = N1 < N2 ? N1 : N2;
+// _Inter.resize(2*minN1N2);
+ _Terminus = minN1N2 < 3;
+
+ /* list of future events ordered according to their coordinates (x,y,z) (lexicographical order) */
+ multimap< const double *, int, VertexLess<DIM> > events;
+ typename multimap< const double *, int, VertexLess<DIM> >::iterator mi1,mi2;
+
+ multimap< int, pair< int,bool> >::iterator mi;
+
+ /********** Initalisation of events with P1 and P2 vertices ************/
+ for(i_loc=0;i_loc<N1;i_loc++)
+ events.insert(make_pair(&P_1[DIM*i_loc],i_loc));
+ for(i_loc=0;i_loc<N2;i_loc++)
+ events.insert(make_pair(&P_2[DIM*i_loc],i_loc+N1));
+
+ if(!_Terminus)
+ {
+ /******** Treatment of the first vertex ********/
+ mi1=events.begin();
+ i_glob = (* mi1).second;
+ bool which_start = i_glob < N1;
+ if(i_glob < N1){ i_next_glob = (i_glob +1)%N1; i_prev_glob = (i_glob -1+N1)%N1;}
+ else{ i_next_glob = (i_glob-N1+1)%N2 + N1;i_prev_glob = (i_glob-N1-1+N2)%N2 + N1;}
+ _Status.insert(make_pair(i_next_glob,make_pair(i_glob, false)));
+ _Status.insert(make_pair(i_prev_glob,make_pair(i_glob, false)));
+ mi1++;
+
+ /******* Loop until the second polygon is reached *******/
+ while( !four_neighbours)
+ {
+ i_glob=(* mi1).second;//global index of vertex i
+ nb_prev = _Status.count(i_glob);//counts the number of segments ending at i
+ switch (nb_prev)
+ {
+ case 1 :
+ mi=_Status.find(i_glob);// pointer to the segment ending at i
+ i_prev_glob = ((*mi).second).first;//starting point of the segment ending at i
+ i_next= (i_prev_glob - i_glob > 0) == (abs(i_prev_glob - i_glob) == 1) ? i_glob - 1 : i_glob + 1;
+ if(i_glob < N1) i_next_glob = (i_next +N1)%N1;
+ else i_next_glob = (i_next-N1+N2)%N2 + N1;
+ _Status.erase(mi);
+ _Status.insert(make_pair(i_next_glob,make_pair(i_glob, false)));
+ mi1++;
+ break;
+ case 2 :
+ return _Inter;
+ case 0 :
+ if( (i_glob < N1) != which_start)
+ {
+ four_neighbours = true;
+ /* detection of double points */
+ const double * Pi1=(* mi1).first;
+ mi2=mi1++;
+ const double * Pi2=(* mi2).first;
+ int i_double= (*mi2).second;
+ if(distance2<DIM>(Pi1, Pi2) < _Epsilon && _Status.count(i_double)==1)
+ {
+ i_glob=i_double;
+ nb_prev=1;
+ events.erase(mi2);
+ }
+ mi1--;
+ }
+ break;
+ default:
+ throw MED_EXCEPTION
+ ("intersect_convex_polygons: sequence of nodes does not describe a simple polygon (1)");
+ }
+ }
+ /******** Loop until a terminal point or crossing is reached ************/
+ while( !_Terminus)
+ {
+ switch (nb_prev)
+ {
+ case 1 :
+ mi=_Status.find(i_glob);// pointer to the segment ending at i
+ i_prev_glob = ((*mi).second).first;//starting point of the segment ending at i
+ sign = (i_prev_glob - i_glob > 0) == (abs(i_prev_glob - i_glob) == 1) ? - 1 : + 1;
+ i_next_glob = i_glob+sign;
+ _Is_in_intersection = ((*mi).second).second;//boolean that tells if i is in the intersection
+ _Status.erase(mi);
+ define_indices(i_loc,i_next,i_prev, Poly1,Poly2,
+ j1,j1_glob,j2,j2_glob,j3,j3_glob,j4,j4_glob,
+ i_glob,i_next_glob,i_prev_glob, P_1,P_2, N1, N2, sign);
+ if( _Is_in_intersection ) add_new_vertex(i_loc, i_glob, i_next_glob, i_prev_glob, Poly1);
+ add_crossings(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob,
+ &Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob,
+ &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob,
+ &Poly1[DIM*i_prev]);
+ break;
+ case 2 :
+ if(_Inter.size()>0)
+ if(i_glob < N1) for(idim=0;idim<DIM;idim++) _Inter.push_back(P_1[DIM*i_glob+idim]);
+ else for(idim=0;idim<DIM;idim++) _Inter.push_back(P_2[DIM*(i_glob-N1)+idim]);
+ return _Inter;
+ case 0 ://MN: à virer d'ici
+ if(_Inter.size()==0){
+ i_next_glob = i_glob+1;
+ i_prev_glob = i_glob-1;
+ define_indices(i_loc,i_next,i_prev, Poly1,Poly2,
+ j1,j1_glob,j2,j2_glob,j3,j3_glob,j4,j4_glob,
+ i_glob,i_next_glob,i_prev_glob, P_1,P_2, N1, N2, 1);
+ double V12[DIM], V34[DIM];
+ double inside = check_inside<DIM>(&Poly1[DIM*i_loc],&Poly2[DIM*j1],&Poly2[DIM*j2],
+ &Poly2[DIM*j3], &Poly2[DIM*j4],V12, V34);
+ _Is_in_intersection=( inside < _Epsilon ); // <= epsilon or 0 ?
+
+ if(abs(inside) > _Epsilon)//vertex clearly inside or outside
+ {
+ if( _Is_in_intersection)
+ {
+ for(int idim=0;idim<DIM;idim++) _Inter.push_back(Poly1[DIM*i_loc+idim]);
+ _End_segments.push_back(make_pair(i_glob,i_next_glob));
+ _End_segments.push_back(make_pair(i_glob,i_prev_glob));
+ }
+ add_crossings(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob,
+ &Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob,
+ &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob,
+ &Poly1[DIM*i_prev]);
+ add_crossings(&Poly1[DIM*i_loc], &Poly1[DIM*i_prev], i_glob, i_prev_glob,
+ &Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob,
+ &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob,
+ &Poly1[DIM*i_next]);
+ }
+ else //vertex on an edge
+ {
+ bool inside_next, inside_prev;
+ double Vnext[DIM], Vprev[DIM];
+ for(idim=0;idim<DIM;idim++) _Inter.push_back(Poly1[DIM*i_loc+idim]);
+
+ if(dotprod<DIM>(V34,V34) > _Epsilon)//vertex i on edge (j1,j2), not on (j3,j4)
+ {
+ crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_next],Vnext);
+ crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_prev],Vprev);
+ inside_next= (dotprod<DIM>(Vnext,V34)<0);
+ inside_prev= (dotprod<DIM>(Vprev,V34)<0);
+
+ if(!(inside_next || inside_prev)) return deque< double >();
+
+ if(inside_next)
+ {
+ _End_segments.push_back(make_pair(i_glob,i_next_glob));
+ add_crossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob,
+ &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob);
+ }
+ else
+ {
+ _End_segments.push_back(make_pair(j1_glob,j2_glob));
+ _Status.insert(make_pair(i_next_glob,make_pair(i_glob, false)));
+ mi=_Status.find(j2_glob);
+ ((* mi).second).second= !((* mi).second).second;
+ }
+ if(inside_prev)
+ {
+ _End_segments.push_back(make_pair(i_glob,i_prev_glob));
+ add_crossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_prev], i_glob, i_prev_glob,
+ &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob);
+ }
+ else
+ {
+ _End_segments.push_back(make_pair(j1_glob,j2_glob));
+ _Status.insert(make_pair(i_prev_glob,make_pair(i_glob, false)));
+ mi=_Status.find(j2_glob);
+ ((* mi).second).second= !((* mi).second).second;
+ }
+ }
+ else if(dotprod<DIM>(V12,V12) > _Epsilon)//vertex i on a edge (j3,j4), not on (j1,j2)
+ {
+ crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_next],Vnext);
+ crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_prev],Vprev);
+ inside_next= dotprod<DIM>(Vnext,V12)<0;
+ inside_prev= dotprod<DIM>(Vprev,V12)<0;
+
+ if(!(inside_next || inside_prev)) return deque< double >();
+
+ if(inside_next)
+ {
+ _End_segments.push_back(make_pair(i_glob,i_next_glob));
+ add_crossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob,
+ &Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob);
+ }
+ else
+ {
+ _End_segments.push_back(make_pair(j3_glob,j4_glob));
+ _Status.insert(make_pair(i_next_glob,make_pair(i_glob, false)));
+ mi=_Status.find(j4_glob);
+ ((* mi).second).second= ! ((* mi).second).second;
+ }
+ if(inside_prev)
+ {
+ _End_segments.push_back(make_pair(i_glob,i_prev_glob));
+ add_crossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_prev], i_glob, i_prev_glob,
+ &Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob);
+ }
+ else
+ {
+ _End_segments.push_back(make_pair(j3_glob,j4_glob));
+ _Status.insert(make_pair(i_prev_glob,make_pair(i_glob, false)));
+ mi=_Status.find(j4_glob);
+ ((* mi).second).second= !((* mi).second).second;
+ }
+ }
+ else //vertices i, j1 and j3 share the same coordinates
+ {
+ crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_next],Vnext);
+ crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_prev],Vprev);
+ crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_next],V12);
+ crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_prev],V34);
+
+ inside_next= dotprod<DIM>(Vnext,V12)<=_Epsilon;
+ inside_prev= dotprod<DIM>(Vprev,V34)<=_Epsilon;
+ bool inside_j2= dotprod<DIM>(Vnext,Vprev)<-_Epsilon;
+
+ if(!( inside_next || inside_prev || inside_j2 )) return deque< double >();
+ if(inside_next) _End_segments.push_back(make_pair(i_glob,i_next_glob));
+ if(inside_prev) _End_segments.push_back(make_pair(i_glob,i_prev_glob));
+ if(inside_j2)
+ {
+ _End_segments.push_back(make_pair(j1_glob,j2_glob));
+ mi=_Status.find(j2_glob);
+ ((* mi).second).second= !((* mi).second).second;
+ }
+ if(!( (inside_next && (inside_prev || inside_j2)) || (inside_prev && inside_j2)))
+ {
+ _End_segments.push_back(make_pair(j3_glob,j4_glob));
+ mi=_Status.find(j4_glob);
+ ((* mi).second).second= !((* mi).second).second;
+ }
+ _Status.insert(make_pair(i_prev_glob,make_pair(i_glob,inside_prev)));
+ _Status.insert(make_pair(i_next_glob,make_pair(i_glob,inside_next)));
+ }
+ }
+ }
+ break;
+ default:
+ cout<<"Problem: nbprev= "<< nb_prev << " ; i_glob = "<< i_glob << endl;
+ throw MED_EXCEPTION
+ ("intersect_convex_polygons: sequence of nodes does not describe a simple polygon (2)");
+ }
+ mi1++;
+ i_glob=(* mi1).second;//global index of vertex i
+ nb_prev = _Status.count(i_glob);
+ }
+ }
+ return _Inter;
+ }
+
+ /**************************************************************************/
+ /* computes the convex hull of a polygon subP which is a sub polygon of P */
+ /* P is the array of coordinates, subP is a map containing initially the indices of a subpolygon of P */
+ /* in the end, subP contains only the elements belonging to the convex hull, and not_in_hull the others */
+ /**************************************************************************/
+ template< int DIM>
+ inline void PolygonAlgorithms<DIM>::conv_hull(const double *P, int N, double * normal,
+ map< int,int >& subP, map< int,int >& not_in_hull,
+ int& NsubP, const double epsilon)
+ {
+ if(NsubP>3)
+ {
+ map< int,int >::iterator mi_prev = subP.begin();
+ map< int,int >::iterator mi = mi_prev++;
+ map< int,int >::iterator mi_next = mi++;
+ double directframe=0.;
+
+ /* Check if the polygon subP is positively oriented */
+ map< int,int >::iterator mi1=mi;
+ while(mi1 != subP.end() && distance2<DIM>(&P[DIM*(*subP.begin()).second],&P[DIM*(*mi1).second])< epsilon)
+ mi1++;
+ map< int,int >::iterator mi2=mi1;
+ while(mi2 != subP.end() && abs(directframe)<epsilon)
+ {
+ directframe =direct_frame<DIM>(&P[DIM* (*mi1).second],
+ &P[DIM* (*subP.begin()).second],
+ &P[DIM* (*mi2).second], normal);
+ mi2++;
+ }
+ if(directframe < 0) for(int idim=0; idim< DIM; idim++) normal[idim] *= -1;
+
+ /* Core of the algorithm */
+ while(mi_next != subP.end())
+ {
+ directframe = direct_frame<DIM>(&P[DIM* (*mi).second],
+ &P[DIM* (*mi_prev).second],
+ &P[DIM* (*mi_next).second], normal);
+ if(directframe > -epsilon){
+ mi = mi++;
+ mi_prev=mi_prev++;
+ mi_next=mi_next++;
+ }
+ else
+ {
+ not_in_hull.insert(*mi);
+ subP.erase(mi);
+ NsubP--;
+ mi--;
+ }
+ }
+ directframe = direct_frame<DIM>(&P[DIM*(*mi).second],
+ &P[DIM*(*mi_prev).second],
+ &P[DIM*(*subP.begin()).second], normal);
+ if(directframe < -epsilon)
+ {
+ not_in_hull.insert(*mi);
+ subP.erase(mi);
+ NsubP--;
+ }
+ }
+ }
+
+ template< int DIM>
+ void PolygonAlgorithms<DIM>::convex_decomposition(const double * P, int N, double *normal, vector< int > subP, int NsubP,
+ vector< map< int,int > >& components, vector< int >& components_index,
+ int& Ncomp, int sign, const double epsilon)
+ {
+ int i;
+ map< int, int > hull;
+ map< int, int > not_in_hull;
+ map< int, int >::iterator mi, mj;
+ vector< int > reflex_region;
+ int Nreflex;
+ int i_xmax=0;
+ const double * xmax=&P[DIM*subP[0]];
+ /* checking an extremal point of subP */
+ for(i=0; i<NsubP; i++)
+ {
+ if(&P[DIM*subP[i]]> xmax)
+ {
+ i_xmax=i;
+ xmax=&P[DIM*subP[i]];
+ }
+ }
+ /* renumbering of SubP elements for the convex hull*/
+ for(i=0; i<NsubP; i++) hull.insert(hull.end(),make_pair(i,subP[(i+i_xmax)%NsubP]));
+ /* compute the convex hull */
+ conv_hull(P, N, normal, hull, not_in_hull, NsubP,epsilon);
+ /* convex hull is the next component */
+ components.push_back(hull);
+ components_index.push_back(sign*NsubP);
+ Ncomp++;
+ /* searching for reflex regions */
+ for(mi=not_in_hull.begin(); mi!=not_in_hull.end(); mi++)
+ {
+ reflex_region.clear();
+ reflex_region.push_back(hull[(*mi).first-1]);
+ reflex_region.push_back( (*mi).second );
+ Nreflex=2;
+ mj=mi;
+ mj++;
+ while((mj != not_in_hull.end()) && ((*mj).first == (*mi).first+1))
+ {
+ reflex_region.push_back((*mj).second);
+ Nreflex++;
+ mi++;
+ mj++;
+ }
+ reflex_region.push_back(hull[(*mi).first+1]);
+ Nreflex++;
+ convex_decomposition( P, N,normal, reflex_region, Nreflex, components, components_index, Ncomp, -sign, epsilon);
+ }
+ }
+
+ /**************************************************************************/
+ /* decomposes a non convex polygon P with N vertices contained in a plane */
+ /* into a sequence of convex polygons */
+ /* the input vectors 'components' and 'components_index' should be empty */
+ /* returns the number of convex components */
+ /* if P is composed of a single point, then an empty polygon is returned */
+ /**************************************************************************/
+ template< int DIM>
+ int PolygonAlgorithms<DIM>::convex_decomposition(const double * P, int N, vector< map< int,int > >& components,
+ vector< int >& components_index, const double epsilon)
+ {
+ int Ncomp=0;
+ vector< int > subP(N);
+ double normal[3]={0,0,0};
+
+ for(int i = 0; i<N; i++) subP[i]=i;
+
+ //Build the normal of polygon P
+ int i1=1;
+ while(i1<N && distance2<DIM>(&P[0],&P[i1])< epsilon) i1++;
+ int i2=i1+1;
+ while(i2<N && abs(dotprod<DIM>(normal,normal))<epsilon)
+ {
+ crossprod<DIM>(&P[i1], &P[0], &P[i2],normal);
+ i2++;
+ }
+
+ convex_decomposition(P, N, normal, subP, N, components, components_index, Ncomp, 1, epsilon);
+ return Ncomp;
+ }
+};
+
+#endif
--- /dev/null
+#ifndef _PolygonAlgorithms_HXX_
+#define _PolygonAlgorithms_HXX_
+
+#include <deque>
+#include <map>
+
+namespace INTERP_UTILS
+{
+ template< int DIM>
+ class VertexLess
+ {
+ public:
+ bool operator()(const double * P_1, const double * P_2)
+ {
+ for(int idim=0; idim<DIM; idim++)
+ {
+ if(P_1[idim] < P_2[idim] ) return true;
+ else if( P_1[idim] > P_2[idim]) return false;
+ }
+ return false;
+ }
+ };
+
+ template< int DIM>
+ class PolygonAlgorithms
+ {
+ public:
+ int _Inter_size;//size of the intersection
+
+ PolygonAlgorithms(double epsilon, double precision);
+ deque<double> intersect_convex_polygons(const double* P_1,const double* P_2, int N1, int N2);
+ //Not yet tested
+ int convex_decomposition(const double * P, int N, vector< map< int,int > >& components,
+ vector< int >& components_index, const double epsilon);
+
+ private:
+ deque< double > _Inter;/* vertices of the intersection P1^P2 */
+ vector< pair< int,int > > _End_segments; /* segments containing inter final edges */
+ /* status list of segments (ending point, starting point) intersected by the sweeping line */
+ /* and a boolean true if the ending point is in the intersection */
+ multimap< int, pair< int,bool> > _Status;
+ bool _Is_in_intersection;
+ bool _Terminus;
+ double _Vdouble[DIM];
+ double _Epsilon;
+ double _Precision;
+
+
+ void define_indices(int& i_loc, int& i_next, int& i_prev,
+ const double *& Poly1, const double *& Poly2,
+ int& j1, int& j1_glob, int& j2, int& j2_glob,
+ int& j3, int& j3_glob, int& j4, int& j4_glob,
+ int& i_glob, int& i_next_glob, int& i_prev_glob,
+ const double * P_1, const double * P_2,
+ int N1, int N2, int sign);
+ void add_crossings( const double * A, const double * B, int i , int i_next,
+ const double * C, const double * D, int j1, int j2,
+ const double * E, const double * F, int j3, int j4,
+ const double * G);
+ void add_crossing0(const double * A, const double * B, int i, int i_next,
+ const double * C, const double * D, int j, int j_next);
+ void add_crossing( double * ABCD, pair< int,int > i_i_next, pair< int,int > j_j_next);
+ void add_new_vertex( int i, int i_glob, int i_next_glob, int i_prev_glob, const double * P);
+ bool intersect_segment_segment(const double * A, const double * B, const double * C, const double * D,
+ const double * E, double * V);
+
+
+ //Not yet tested
+ void convex_decomposition(const double* P, int N, double* n, vector< int > subP, int NsubP,
+ vector< map< int,int > >& components, vector< int >& components_index,
+ int& Ncomp, int sign, const double epsilon);
+ void conv_hull(const double *P, int N, double * n, map< int,int >& subP,
+ map< int,int >& not_in_hull, int& NsubP, const double epsilon);
+ };
+};
+
+#include "PolygonAlgorithms.cxx"
+
+#endif
LOG(5, "Loading " << mesh2 << " from " << mesh2path);
MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
- m = interpolator->interpol_maillages(sMesh, tMesh);
+ m = interpolator->interpolateMeshes(sMesh, tMesh);
// if reflexive, check volumes
if(strcmp(mesh1path,mesh2path) == 0)
CPPFLAGS+=$(BOOST_CPPFLAGS) $(MED2_INCLUDES) $(HDF5_INCLUDES) -I$(top_srcdir)/src/INTERP_KERNEL
#include CppUnit
-CXXFLAGS+= -I/usr/include/cppunit
-CPPFLAGS+= -I/usr/include/cppunit
+CXXFLAGS+= ${CPPUNIT_INCLUDES}
+CPPFLAGS+= ${CPPUNIT_INCLUDES}
# for log
CXXFLAGS+= -DLOG_LEVEL=3 -DOPTIMIZE -O2 -Wall
#CPPFLAGS+=-fprofile-arcs -ftest-coverage
#for gprof
-#CXXFLAGS+=-pg
-#CPPFLAGS+=-pg
+CXXFLAGS+=-pg
+CPPFLAGS+=-pg
#LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS)
# change motivated by the bug KERNEL4778.
#LDFLAGS+= $(HDF5_LIBS) $(STDLIB) -lmedmem -lutil -lmed
# for gcov
-#LDFLAGS+= -lgcov
-#LDFLAGS+= -pg
+LDFLAGS+= -lgcov
+LDFLAGS+= -pg
#LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS)
# change motivated by the bug KERNEL4778.
LDFLAGSFORBIN += $(LDFLAGS) -lm $(HDF5_LIBS) \
-L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome \
- -lcppunit -linterpkernel
+ $(CPPUNIT_LIBS) -linterpkernel
#LDFLAGSFORBIN += -pg
void MeshTestToolkit::calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const
{
const string dataBaseDir = getenv("MED_ROOT_DIR");
- const string dataDir = dataBaseDir + string("share/salome/resources/med/");
+ const string dataDir = dataBaseDir + string("/share/salome/resources/med/");
LOG(1, std::endl << "=== -> intersecting src = " << mesh1path << ", target = " << mesh2path );
LOG(5, "Loading " << mesh2 << " from " << mesh2path);
MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2);
- m = _interpolator->interpol_maillages(sMesh, tMesh);
+ m = _interpolator->interpolateMeshes(sMesh, tMesh);
// if reflexive, check volumes
if(strcmp(mesh1path,mesh2path) == 0)
LOG(1, "Target mesh has " << numTargetElems << " elements");
- m = _interpolator->interpol_maillages(sMesh, tMesh);
+ m = _interpolator->interpolateMeshes(sMesh, tMesh);
std::pair<int, int> eff = countNumberOfMatrixEntries(m);
LOG(1, eff.first << " of " << numTargetElems * numSrcElems << " intersections calculated : ratio = "
Interpolation3D* interpolator = new Interpolation3D();
- m = interpolator->interpol_maillages(sMesh, tMesh);
+ m = interpolator->interpolateMeshes(sMesh, tMesh);
std::pair<int, int> eff = countNumberOfMatrixEntries(m);
LOG(1, eff.first << " of " << numTargetElems * numSrcElems << " intersections calculated : ratio = "
#ifndef _TRANSLATIONROTATIONMATRIX_HXX_
#define _TRANSLATIONROTATIONMATRIX_HXX_
+#include "InterpolationUtils.hxx"
+#include <cmath>
+//#include <vector>
+
namespace MEDMEM
{
--- /dev/null
+#include "MEDMEM_Mesh.hxx"
+#include "InterpolationUtils.hxx"
+#include "PlanarIntersector.hxx"
+#include "TriangulationIntersector.hxx"
+
+//using namespace std;
+using namespace MED_EN;
+using namespace INTERP_UTILS;
+//using namespace MEDMEM;
+
+namespace MEDMEM
+{
+ template <int DIM>
+ TriangulationIntersector<DIM>::TriangulationIntersector(const MEDMEM::MESH& mesh_A, const MEDMEM::MESH& mesh_B,
+ double DimCaracteristic, double Precision,
+ double MedianPlane, int PrintLevel)
+ :PlanarIntersector(mesh_A,mesh_B, DimCaracteristic, Precision, MedianPlane)
+ {
+ _Connect_A= mesh_A.getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,
+ MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ _Connect_B= mesh_B.getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,
+ MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ _Conn_index_A= mesh_A.getConnectivityIndex(MED_EN::MED_NODAL, MED_EN::MED_CELL);
+ _Conn_index_B= mesh_B.getConnectivityIndex(MED_EN::MED_NODAL, MED_EN::MED_CELL);
+ _Coords_A = mesh_A.getCoordinates(MED_EN::MED_FULL_INTERLACE);
+ _Coords_B = mesh_B.getCoordinates(MED_EN::MED_FULL_INTERLACE);
+ _DimCaracteristic = DimCaracteristic;
+ _Precision = Precision;
+ _MedianPlane = MedianPlane;
+ _do_rotate = true;
+ _PrintLevel= PrintLevel;
+
+ if( _PrintLevel >= 1)
+ {
+ cout<< "_intersection_type = triangles " <<endl;
+ if(DIM==3) cout<< "_do_rotate = true"<<endl;
+ }
+ }
+
+ template <int DIM>
+ double TriangulationIntersector<DIM>::intersectCells(int icell_A, int icell_B,
+ int nb_NodesA, int nb_NodesB)
+ {
+ double result=0;
+
+ //Obtain the coordinates of A and B
+ vector< double> Coords_A (DIM*nb_NodesA);
+ vector< double> Coords_B (DIM*nb_NodesB);
+ for (int idim=0; idim<DIM; idim++)
+ {
+ for (int i_A=0; i_A<nb_NodesA; i_A++)
+ Coords_A[DIM*i_A+idim] = _Coords_A[DIM*(_Connect_A[_Conn_index_A[icell_A-1]-1+i_A]-1)+idim];
+ for (int i_B=0; i_B<nb_NodesB; i_B++)
+ Coords_B[DIM*i_B+idim] = _Coords_B[DIM*(_Connect_B[_Conn_index_B[icell_B-1]-1+i_B]-1)+idim];
+ }
+
+ //project cells A and B on the median plane and rotate the median plane
+ if(DIM==3) projection(Coords_A, Coords_B, nb_NodesA, nb_NodesB,
+ _DimCaracteristic * _Precision, _MedianPlane, _do_rotate);
+
+ //Compute the intersection area
+ double area[DIM];
+ for(int i_A = 1; i_A<nb_NodesA-1; i_A++)
+ {
+ for(int i_B = 1; i_B<nb_NodesB-1; i_B++)
+ {
+ vector<double> inter;
+ INTERP_UTILS::intersec_de_triangle(&Coords_A[0],&Coords_A[DIM*i_A],&Coords_A[DIM*(i_A+1)],
+ &Coords_B[0],&Coords_B[DIM*i_B],&Coords_B[DIM*(i_B+1)],
+ inter, _DimCaracteristic, _Precision);
+ int nb_inter=((int)inter.size())/2;
+ if(nb_inter >3) inter=reconstruct_polygon(inter);
+ for(int i = 1; i<nb_inter-1; i++)
+ {
+ INTERP_UTILS::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
+ result +=0.5*norm<2>(area);
+ }
+ //DEBUG prints
+ if(_PrintLevel >= 3)
+ {
+ cout<<endl<<"Number of nodes of the intersection = "<< nb_inter <<endl;
+ for(int i=0; i< nb_inter; i++)
+ {for (int idim=0; idim<2; idim++) cout<< inter[2*i+idim]<< " "; cout << endl;}
+ }
+ }
+ }
+
+ if(_PrintLevel >= 3)
+ {
+ cout<<endl<< "Cell coordinates after projection"<<endl;
+ cout<<endl<< "icell_A= " << icell_A << ", nb nodes A= " << nb_NodesA << endl;
+ for(int i_A =0; i_A< nb_NodesA; i_A++)
+ {for (int idim=0; idim<2; idim++) cout<< Coords_A[DIM*i_A+idim]<< " "; cout << endl;}
+ cout<<endl<< "icell_B= " << icell_B << ", nb nodes B= " << nb_NodesB << endl;
+ for(int i_B =0; i_B< nb_NodesB; i_B++)
+ {for (int idim=0; idim<2; idim++) cout<< Coords_B[DIM*i_B+idim]<< " "; cout << endl;}
+ cout<<endl<<"Intersection area= " << result << endl;
+ }
+
+ return result;
+ }
+
+}
--- /dev/null
+#ifndef _TRIANGULATION_INTERSECTOR_HXX_
+#define _TRIANGULATION_INTERSECTOR_HXX_
+
+#include "MEDMEM_Mesh.hxx"
+#include "InterpolationUtils.hxx"
+
+namespace MEDMEM
+{
+ class MESH;
+ template <int DIM>
+ class TriangulationIntersector: public PlanarIntersector
+ {
+ public:
+ TriangulationIntersector (const MEDMEM::MESH& mesh_A, const MEDMEM::MESH& mesh_B,
+ double DimCaracteristic, double Precision, double MedianPlane, int PrintLevel);
+ // virtual ~Intersector();// Warum virtual ?
+ double intersectCells(int icell_A, int icell_B, int nb_NodesA, int nb_NodesB);
+
+ private :
+ const int * _Connect_A;
+ const int * _Connect_B;
+ const double * _Coords_A;
+ const double * _Coords_B;
+ const int *_Conn_index_A;
+ const int *_Conn_index_B;
+ double _DimCaracteristic;
+ double _Precision;
+ double _MedianPlane;
+ bool _do_rotate;
+ int _PrintLevel;
+ };
+}
+
+#endif
MEDMEM::MESH mesh1(MED_DRIVER,"/home/vb144235/resources/DividedUnitTetraSimpler.med","DividedUnitTetraSimpler");
MEDMEM::MESH mesh2(MED_DRIVER,"/home/vb144235/resources/DividedUnitTetraSimpler.med","DividedUnitTetraSimpler");
Interpolation3D interpolator;
- vector<map<int,double> > matrix = interpolator.interpol_maillages(mesh1,mesh2);
+ vector<map<int,double> > matrix = interpolator.interpolateMeshes(mesh1,mesh2);
// dump
dumpIntersectionMatrix(matrix);
--- /dev/null
+#include "Interpolation2D.hxx"
+#include "MEDMEM_Mesh.hxx"
+
+int main()
+{
+ MEDMEM::MESH source(MED_DRIVER,"/home/vb144235/resources/square128000.med","Mesh_1");
+ MEDMEM::MESH target(MED_DRIVER,"/home/vb144235/resources/square30000.med","Mesh_1");
+
+ MEDMEM::Interpolation2D interp;
+ interp.setOptions(1e-6,1,2);
+ interp.interpolateMeshes(source,target);
+}
--- /dev/null
+// Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// 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.
+//
+// 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
+//
+//
+//
+// File : testUUnit.cxx
+// Module : MED
+
+#include "InterpolationUtils.hxx"
+#include "PolygonAlgorithms.hxx"
+#include "PolygonAlgorithms.cxx"
+
+using namespace std;
+
+int main ()
+{
+ cerr<< endl;
+ cerr<< " ************ Test 1 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_1(1e-6,1e-6);
+
+ const double Losange1[12]=
+ {
+ 1,0,0,
+ 0,1,0,
+ -1,0,0,
+ 0,-1,0
+ };
+
+ const double Losange2[12]=
+ {
+ 2,0,0,
+ 1,1,0,
+ 0,0,0,
+ 1,-1,0
+ };
+
+ deque< double > resultat_test1 = P_1.intersect_convex_polygons(Losange1,Losange2,4,4);
+
+ cerr<< " Test1: Résultat théorique" << endl;
+ cerr<< 0.5 << " ,"<< -0.5 << " ," << 0 << endl;
+ cerr<< 0 << " ,"<< 0 << " ,"<< 0 << endl;
+ cerr<< 0.5 << " ,"<< 0.5 << " ,"<< 0 << endl;
+ cerr<< 1 << " ," << 0 << " ," << 0 << endl;
+
+ cerr<< " Test1: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test1.size()/3; i++)
+ {
+ cerr << resultat_test1[3*i] << " ";
+ cerr << resultat_test1[3*i+1] << " ";
+ cerr << resultat_test1[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ************ Test 2 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_2(1e-6,1e-6);
+
+ const double Losange3[12]=
+ {
+ 2.5,0.5,0,
+ 1.5,1.5,0,
+ 0.5,0.5,0,
+ 1.5,-0.5,0
+ };
+ deque< double > resultat_test2 = P_2.intersect_convex_polygons(Losange1,Losange3,4,4);
+
+ cerr<< " Test2: Résultat théorique" << endl;
+ // cerr<< 0.5 << " ,"<< 0.5 << " ," << 0 << endl;
+ // cerr<< 1 << " ,"<< 0 << " ," << 0 << endl;
+
+ cerr<< " Test2: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test2.size()/3; i++)
+ {
+ cerr << resultat_test2[3*i] << " ";
+ cerr << resultat_test2[3*i+1] << " ";
+ cerr << resultat_test2[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ************ Test 3 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_3(1e-6,1e-6);
+
+ const double Carre1[12]=
+ {
+ -1,-1,0,
+ -1,1,0,
+ 1,1,0,
+ 1,-1,0
+ };
+ const double Carre2[12]=
+ {
+ 1,-0.25,0,
+ 0,-0.25,0,
+ 0,0.25,0,
+ 1,0.25,0
+ };
+ deque< double > resultat_test3 = P_3.intersect_convex_polygons( Carre1, Carre2,4,4);
+
+ cerr<< " Test3: Résultat théorique" << endl;
+ cerr<< 0 << " ,"<< 0.25 << " ," << 0 << endl;
+ cerr<< 0 << " ,"<< -0.25 << " ," << 0 << endl;
+ cerr<< 1 << " ,"<< -0.25 << " ,"<< 0 << endl;
+ cerr<< 1 << " ,"<< 0.25 << " ," << 0 << endl;
+
+ cerr<< " Test3: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test3.size()/3; i++)
+ {
+ cerr << resultat_test3[3*i] << " ";
+ cerr << resultat_test3[3*i+1] << " ";
+ cerr << resultat_test3[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ***************** Test 4 ***************** " << endl;
+
+ INTERP_UTILS::PolygonAlgorithms<3> P_4(1e-6,1e-6);
+ const double Losange4[12]=
+ {
+ 3,0,0,
+ 2,1,0,
+ 1,0,0,
+ 2,-1,0
+ };
+ deque< double > resultat_test4 = P_4.intersect_convex_polygons( Losange1, Losange4,4,4);
+
+ cerr<< " Test4: Résultat théorique" << endl;
+ cerr<< 1 << " ,"<< 0 << " ,"<< 0 << endl;
+
+ cerr<< " Test4: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test4.size()/3; i++)
+ {
+ cerr << resultat_test4[3*i] << " ";
+ cerr << resultat_test4[3*i+1] << " ";
+ cerr << resultat_test4[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ***************** Test 5 ***************** " << endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_5(1e-6,1e-6);
+ deque< double > resultat_test5 = P_5.intersect_convex_polygons( Carre1, Carre1,4,4);
+
+ cerr<< " Test5: Résultat théorique" << endl;
+ cerr<< -1 << " ,"<< 1 << " ," << 0 << endl;
+ cerr<< -1 << " ,"<< -1 << " ," << 0 << endl;
+ cerr<< 1 << " ,"<< -1 << " ," << 0 << endl;
+ cerr<< 1 << " ,"<< 1 << " ," << 0 << endl;
+
+ cerr<< " Test5: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test5.size()/3; i++)
+ {
+ cerr << resultat_test5[3*i] << " ";
+ cerr << resultat_test5[3*i+1] << " ";
+ cerr << resultat_test5[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ***************** Test 6 ***************** " << endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_6(1e-6,1e-6);
+
+ const double Losange5[12]=
+ {
+ 1.5,0,0,
+ 0,1.5,0,
+ -1.5,0,0,
+ 0,-1.5,0
+ };
+ deque< double > resultat_test6 = P_6.intersect_convex_polygons( Carre1, Losange5,4,4);
+
+ cerr<< " Test6: Résultat théorique" << endl;
+ cerr<< 1 << " ,"<< -0.5 << " ," << 0 << endl;
+ cerr<< 0.5 << " ,"<< -1 << " ," << 0 << endl;
+ cerr<< -0.5 << " ,"<< -1 << " ," << 0 << endl;
+ cerr<< -1 << " ,"<< -0.5 << " ," << 0 << endl;
+ cerr<< -1 << " ,"<< 0.5 << " ," << 0 << endl;
+ cerr<< -0.5 << " ,"<< 1 << " ," << 0 << endl;
+ cerr<< 0.5 << " ,"<< 1 << " ,"<< 0 << endl;
+ cerr<< 1 << " ,"<< 0.5 << " ," << 0 << endl;
+
+ cerr<< " Test6: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test6.size()/3; i++)
+ {
+ cerr << resultat_test6[3*i] << " ";
+ cerr << resultat_test6[3*i+1] << " ";
+ cerr << resultat_test6[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ***************** Test 7 ***************** " << endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_7(1e-6,1e-6);
+
+ deque< double > resultat_test7 = P_7.intersect_convex_polygons( Losange1, Carre1,4,4);
+
+ cerr<< " Test7: Résultat théorique" << endl;
+ cerr<< 0 << " ,"<< -1 << " ," << 0 << endl;
+ cerr<< -1 << " ,"<< 0 << " ," << 0 << endl;
+ cerr<< 0 << " ,"<< 1 << " ,"<< 0 << endl;
+ cerr<< 1 << " ,"<< 0 << " ," << 0 << endl;
+
+ cerr<< " Test7: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test7.size()/3; i++)
+ {
+ cerr << resultat_test7[3*i] << " ";
+ cerr << resultat_test7[3*i+1] << " ";
+ cerr << resultat_test7[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ************ Test 8 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_8(1e-6,1e-6);
+
+ const double Losange6[18]=
+ {
+ 2,0,0,
+ 1,1,0,
+ 0.5,0.5,0,
+ 0,0,0,
+ 0.5,-0.5,0,
+ 1,-1,0
+ };
+ const double Losange7[15]=
+ {
+ 1,0,0,
+ 0,1,0,
+ -1,0,0,
+ 0,-1,0,
+ 0.5,-0.5,0
+ };
+
+ deque< double > resultat_test8 = P_8.intersect_convex_polygons(Losange6,Losange7,6,5);
+
+ cerr<< " Test8: Résultat théorique" << endl;
+ cerr<< 0.5 << " ,"<< -0.5 << " ," << 0 << endl;
+ cerr<< 0.5 << " ,"<< -0.5 << " ," << 0 << endl;
+ cerr<< 0 << " ,"<< 0 << " ,"<< 0 << endl;
+ cerr<< 0.5 << " ,"<< 0.5 << " ,"<< 0 << endl;
+ cerr<< 0.5 << " ,"<< 0.5 << " ,"<< 0 << endl;
+ cerr<< 1 << " ," << 0 << " ," << 0 << endl;
+
+ cerr<< " Test8: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test8.size()/3; i++)
+ {
+ cerr << resultat_test8[3*i] << " ";
+ cerr << resultat_test8[3*i+1] << " ";
+ cerr << resultat_test8[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ************ Test 9 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_9(1e-6,1e-6);
+
+ const double Carre3[15]=
+ {
+ -1,-1,0,
+ -1,1,0,
+ 0.5,1,0,
+ 1,1,0,
+ 1,-1,0
+ };
+ const double Carre4[12]=
+ {
+ -0.5,-1,0,
+ -0.5,1,0,
+ 1.5,1,0,
+ 1.5,-1,0
+ };
+
+ deque< double > resultat_test9 = P_9.intersect_convex_polygons(Carre4,Carre3,4,5);
+
+ cerr<< " Test9: Résultat théorique" << endl;
+ cerr<< -0.5 << " ,"<< 1 << " ," << 0 << endl;
+ cerr<< -0.5 << " ,"<< -1 << " ," << 0 << endl;
+ cerr<< 1 << " ,"<< -1 << " ," << 0 << endl;
+ cerr<< 1 << " ,"<< 1 << " ," << 0 << endl;
+
+ cerr<< " Test9: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test9.size()/3; i++)
+ {
+ cerr << resultat_test9[3*i] << " ";
+ cerr << resultat_test9[3*i+1] << " ";
+ cerr << resultat_test9[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ************ Test 10 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_10(1e-6,1e-6);
+
+ const double Carre5[15]=
+ {
+ -1,-1,0,
+ -1,1,0,
+ 0,1,0,
+ 1,1,0,
+ 1,-1,0
+ };
+ const double Losange8[12]=
+ {
+ 0,1,0,
+ 1,-1,0,
+ 0,-1.5,0,
+ -0.5,-1,0
+ };
+ deque< double > resultat_test10 = P_10.intersect_convex_polygons(Losange8,Carre5,4,5);
+ cerr<< " Test10: Résultat théorique" << endl;
+ cerr<< 0 << " ,"<< 1 << " ," << 0 << endl;
+ cerr<< -0.5 << " ,"<< -1 << " ," << 0 << endl;
+ cerr<< 1 << " ,"<< -1 << " ," << 0 << endl;
+
+ cerr<< " Test10: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test10.size()/3; i++)
+ {
+ cerr << resultat_test10[3*i] << " ";
+ cerr << resultat_test10[3*i+1] << " ";
+ cerr << resultat_test10[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ************ Test 11 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_11(1e-6,1e-6);
+
+ const double Losange9[12]=
+ {
+ 0.5,0,0,
+ 0,1,0,
+ -1.5,0,0,
+ 0,-1,0
+ };
+
+ deque< double > resultat_test11 = P_11.intersect_convex_polygons(Losange1,Losange9,4,4);
+
+ cerr<< " Test11: Résultat théorique" << endl;
+ cerr<< 0 << " ,"<< -1 << " ," << 0 << endl;
+ cerr<< 0 << " ,"<< -1 << " ," << 0 << endl;
+ cerr<< -1 << " ,"<< 0 << " ,"<< 0 << endl;
+ cerr<< 0 << " ,"<< 1 << " ,"<< 0 << endl;
+ cerr<< 0 << " ,"<< 1 << " ,"<< 0 << endl;
+ cerr<< 0.5 << " ," << 0 << " ," << 0 << endl;
+
+ cerr<< " Test11: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test11.size()/3; i++)
+ {
+ cerr << resultat_test11[3*i] << " ";
+ cerr << resultat_test11[3*i+1] << " ";
+ cerr << resultat_test11[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ************ Test 12 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_12(1e-6,1e-6);
+
+ const double hexagone1[18]=
+ {
+ -2,0,0,
+ -1,-1,0,
+ 1,-1,0,
+ 2,0,0,
+ 1,1,0,
+ -1,1,0
+ };
+ const double hexagone2[18]=
+ {
+ -1.5,0.5,0,
+ -1,-1,0,
+ 1,-1,0,
+ 2,1,0,
+ 1,1,0,
+ -1,1,0
+ };
+ deque< double > resultat_test12 = P_12.intersect_convex_polygons(hexagone1,hexagone2,6,6);
+
+ cerr<< " Test12: Résultat théorique" << endl;
+ cerr<< 1 << " ,"<< -1 << " ," << 0 << endl;
+ cerr<< -1 << " ,"<< -1 << " ,"<< 0 << endl;
+ cerr<< -1.5 << " ,"<< 0.5 << " ,"<< 0 << endl;
+ cerr<< -1 << " ," << 1 << " ," << 0 << endl;
+ cerr<< 1 << " ,"<< 1 << " ," << 0 << endl;
+ cerr<< 5./3 << " ,"<< 1/3. << " ," << 0 << endl;
+
+ cerr<< " Test12: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test12.size()/3; i++)
+ {
+ cerr << resultat_test12[3*i] << " ";
+ cerr << resultat_test12[3*i+1] << " ";
+ cerr << resultat_test12[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ************ Test 13 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_13(1e-6,1e-6);
+
+ const double hexagone3[18]=
+ {
+ -2,2,0,
+ -1,1,0,
+ 1,1,0,
+ 2,2,0,
+ 1,3,0,
+ -1,3,0
+ };
+ deque< double > resultat_test13 = P_13.intersect_convex_polygons(hexagone1,hexagone3,6,6);
+
+ cerr<< " Test13: Résultat théorique" << endl;
+
+ cerr<< " Test13: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test13.size()/3; i++)
+ {
+ cerr << resultat_test13[3*i] << " ";
+ cerr << resultat_test13[3*i+1] << " ";
+ cerr << resultat_test13[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ************ Test 14 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_14(1e-6,1e-6);
+
+ const double Carre6[12]=
+ {
+ -1,1,0,
+ -1,3,0,
+ 0.5,3,0,
+ 0.5,1,0
+ };
+ deque< double > resultat_test14 = P_14.intersect_convex_polygons( Carre1, Carre6,4,4);
+
+ cerr<< " Test14: Résultat théorique" << endl;
+
+ cerr<< " Test14: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test14.size()/3; i++)
+ {
+ cerr << resultat_test14[3*i] << " ";
+ cerr << resultat_test14[3*i+1] << " ";
+ cerr << resultat_test14[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ***************** Test 15 ***************** " << endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_15(1e-6,1e-6);
+
+ const double Losange10[12]=
+ {
+ 0,-1,0,
+ 1,-2,0,
+ 0,-3,0,
+ -1,-2,0
+ };
+ deque< double > resultat_test15 = P_15.intersect_convex_polygons( Losange1, Losange10,4,4);
+
+ cerr<< " Test15: Résultat théorique" << endl;
+
+ cerr<< " Test15: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test15.size()/3; i++)
+ {
+ cerr << resultat_test15[3*i] << " ";
+ cerr << resultat_test15[3*i+1] << " ";
+ cerr << resultat_test15[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ************ Test 16 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_16(1e-6,1e-6);
+
+
+ const double triangle1[9]=
+ {
+ 0.5,0,0,
+ 1,1,0,
+ 0,1,0
+ };
+
+ deque< double > resultat_test16 = P_16.intersect_convex_polygons(Losange1,triangle1,4,3);
+
+ cerr<< " Test16: Résultat théorique" << endl;
+ cerr<< 0.5 << " ,"<< 0 << " ,"<< 0 << endl;
+ cerr<< 0 << " ,"<< 1 << " ," << 0 << endl;
+ cerr<< 2/3 << " ," << 1/3 << " ," << 0 << endl;
+
+ cerr<< " Test16: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test16.size()/3; i++)
+ {
+ cerr << resultat_test16[3*i] << " ";
+ cerr << resultat_test16[3*i+1] << " ";
+ cerr << resultat_test16[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+ cerr<< " ************ Test 17 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_17(1e-6,1e-6);
+
+
+ const double triangle2[9]=
+ {
+ 0,0.5,0,
+ 0,-0.5,0,
+ 1.5,0,0
+ };
+
+ deque< double > resultat_test17 = P_17.intersect_convex_polygons(Losange1,triangle1,4,3);
+
+ cerr<< " Test17: Résultat théorique" << endl;
+ cerr<< 0 << " ,"<< 0.5 << " ,"<< 0 << endl;
+ cerr<< 0 << " ,"<< -0.5 << " ," << 0 << endl;
+ cerr<< "??" << " ," << 1/3 << " ," << 0 << endl;
+
+ cerr<< " Test17: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test17.size()/3; i++)
+ {
+ cerr << resultat_test17[3*i] << " ";
+ cerr << resultat_test17[3*i+1] << " ";
+ cerr << resultat_test17[3*i+2] << " ";
+ cerr << endl ;
+ }
+
+};
