--- /dev/null
+#include <iostream>
+#include <vector>
+#include <algorithm>
+
+const int MIN_NB_ELEMS =10;
+const int MAX_LEVEL=10;
+using namespace std;
+
+template <int dim> class BBTree
+{
+
+private:
+ BBTree* _left;
+ BBTree* _right;
+ BBTree* _center;
+ int _level;
+ double _median;
+ double* _bb;
+ vector <int> _elems;
+ bool _terminal;
+
+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)
+{
+ if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL)
+ {
+ _terminal=true;
+
+ }
+ vector<double> nodes (nbelems*2);
+ for (int i=0; i<nbelems; i++)
+ {
+ int elem;
+ if (elems!=0)
+ elem= elems[i];
+ 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];
+ }
+ if (_terminal) return;
+
+ //elem list is not useful if the node is not terminal
+ _elems.clear();
+
+ vector<double>::iterator median = nodes.begin()+nbelems;
+ nth_element< vector<double>::iterator>(nodes.begin(), nodes.begin()+nbelems, nodes.end());
+ _median=*median;
+ // std:: cout << *median <<std::endl;
+
+ vector<int> new_elems_left;
+ vector<int> new_elems_right;
+ vector<int> new_elems_center;
+ for (int i=0; i<nbelems;i++)
+ {
+ int elem;
+ if (elems!=0)
+ elem= elems[i];
+ 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)
+ {
+ new_elems_right.push_back(elem);
+ continue;
+ }
+ new_elems_center.push_back(elem);
+ }
+
+ _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());
+
+
+}
+
+~BBTree()
+{
+ if (_left!=0) delete _left;
+ if (_right!=0) delete _right;
+ if (_center!=0) delete _center;
+}
+
+void getIntersectingElems(double* bb, vector<int>& elems)
+{
+ // terminal node : return list of elements intersecting bb
+ if (_terminal)
+ {
+ for (int i=0; i<_elems.size(); i++)
+ {
+ double* bb_ptr=_bb+_elems[i]*2*dim;
+ 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])
+ intersects=false;
+ }
+ if (intersects)
+ 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)
+ {
+ _left->getIntersectingElems(bb, elems);
+ _center->getIntersectingElems(bb,elems);
+ return;
+ }
+ if (min > _median)
+ {
+ _right->getIntersectingElems(bb,elems);
+ _center->getIntersectingElems(bb,elems);
+ return;
+ }
+ _left->getIntersectingElems(bb,elems);
+ _center->getIntersectingElems(bb,elems);
+ _right->getIntersectingElems(bb,elems);
+}
+};
ElementLocator::ElementLocator(const ParaMESH& mesh, const ProcessorGroup& distant_group)
:_local_mesh(mesh.getMesh()),
_local_group(*mesh.getBlockTopology()->getProcGroup()),
-_distant_group(distant_group)
+ _distant_group(distant_group),
+ _adjustment_eps(0.1)
{
_union_group = _local_group.fuse(distant_group);
_computeBoundingBoxes();
bool ElementLocator::_intersectsBoundingBox(double* bb1, double* bb2, int dim)
{
+ double bbtemp[2*dim];
+ double deltamax=0.0;
+ for (int i=0; i< dim; i++)
+ {
+ double delta = bb1[2*i+1]-bb1[2*i];
+ deltamax = (delta>deltamax)?delta:deltamax;
+ }
+ for (int i=0; i<dim; i++)
+ {
+ bbtemp[i*2]=bb1[i*2]-deltamax*_adjustment_eps;
+ bbtemp[i*2+1]=bb1[i*2+1]+deltamax*_adjustment_eps;
+ }
+
for (int idim=0; idim < dim; idim++)
{
- bool intersects = bb1[idim*2]<bb2[idim*2+1] && bb2[idim*2]<bb1[idim*2+1];
+ bool intersects = bbtemp[idim*2]<bb2[idim*2+1] && bb2[idim*2]<bbtemp[idim*2+1];
if (!intersects) return false;
}
return true;
MEDMEM::MESHING* meshing = new MEDMEM::MESHING ();
int* recv_buffer_ptr = recv_buffer;
meshing->setCoordinates(distant_space_dim, distant_nnodes, recv_coords, "CARTESIAN", MED_EN::MED_FULL_INTERLACE);
+
meshing->setNumberOfTypes(distant_nb_types,MED_EN::MED_CELL);
// converting the types from int to medGeometryElement
{
distant_ids_recv[i]=*recv_buffer_ptr++;
}
-
+ meshing->setMeshDimension(distant_mesh_dim);
+
distant_mesh=meshing;
delete[] recv_buffer;
delete[] nbtypes;
meshing->setNumberOfTypes(nbelems_per_type.size(),MED_EN::MED_CELL);
meshing->setTypes(new_types,MED_EN::MED_CELL);
meshing->setNumberOfElements(nbtypes,MED_EN::MED_CELL);
-
+
+ int dimmax=0;
int* conn_ptr= conn;
for (int i=0; i<nbelems_per_type.size(); i++)
{
meshing->setConnectivity(conn_ptr, MED_EN::MED_CELL,new_types[i]);
conn_ptr+=nbtypes[i]*(new_types[i]%100);
+ if (new_types[i]/100>dimmax) dimmax=new_types[i]/100;
}
+ meshing->setMeshDimension(dimmax);
delete [] small2big;
delete [] nbtypes;
delete [] conn;
const ProcessorGroup& _local_group;
ProcessorGroup* _union_group;
std::vector<int> _distant_proc_ids;
+ double _adjustment_eps;
//InterpolationMatrix _matrix;
//MxN_Mapping _mapping;
--- /dev/null
+#ifndef _INTERPOLATION_HXX_
+#define _INTERPOLATION_HXX_
+
+
+#include <vector>
+#include <map>
+
+#include "MEDMEM_Mesh.hxx"
+
+
+namespace ParaMEDMEM
+{
+
+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;
+
+};
+
+};
+#endif
+#include "INTERPOLATION_Utils.hxx"
#include "INTERPOLATION_2D.hxx"
using namespace std;
using namespace MED_EN;
using namespace MEDMEM;
-
-struct monMaillageS
-{
- int _NbMaille;
- int _NbNoeud;
- const int* _Connect;
- const double* _Coord;
- const int* _ReversNConnect;
- const int* _ReversNConnectIndex;
- double _DimCaracteristic;
-};
-
-struct monMaillageP
+namespace ParaMEDMEM
{
- int _NbMaille;
- int _NbNoeud;
- const int* _Connect;
- const double* _Coord;
- double _DimCaracteristic;
-};
-
-
INTERPOLATION_2D::INTERPOLATION_2D()
{
_PrintLevel=PrintLevel;
}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/* calcul la surface d'un triangle */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-double INTERPOLATION_2D::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);
- return Surface;
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/* fonction qui calcul le déterminant */
-/* de deux vecteur(cf doc CGAL). */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-
-//fonction qui calcul le déterminant des vecteurs: P3P1 et P3P2
-//(cf doc CGAL).
-
-double INTERPOLATION_2D::mon_determinant(const double* P_1,const double* P_2,
- const double* P_3)
-{
- double mon_det=(P_1[0]-P_3[0])*(P_2[1]-P_3[1])-(P_2[0]-P_3[0])*(P_1[1]-P_3[1]);
- return mon_det;
-}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-//calcul la norme du vecteur P1P2
-
-double INTERPOLATION_2D::norme_vecteur(const double* P_1,const double* P_2)
-{
- double X=P_1[0]-P_2[0];
- double Y=P_1[1]-P_2[1];
- double norme=sqrt(X*X+Y*Y);
- return norme;
-}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/* calcul le cos et le sin de l'angle P1P2,P1P3 */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-vector<double> INTERPOLATION_2D::calcul_cos_et_sin(const double* P_1,const double* P_2,
- const double* P_3)
-{
-
- vector<double> Vect;
- double P1_P2=norme_vecteur(P_1,P_2);
- double P2_P3=norme_vecteur(P_2,P_3);
- double P3_P1=norme_vecteur(P_3,P_1);
-
- double N=P1_P2*P1_P2+P3_P1*P3_P1-P2_P3*P2_P3;
- double D=2.0*P1_P2*P3_P1;
- double COS=N/D;
- Vect.push_back(COS);
- double V=mon_determinant(P_2,P_3,P_1);
- double D_1=P1_P2*P3_P1;
- double SIN=V/D_1;
- Vect.push_back(SIN);
-
- return Vect;
-
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/*fonction pour vérifier qu'un point n'a pas déja été considérer dans */
-/* le vecteur et le rajouter au vecteur sinon. */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-void INTERPOLATION_2D::verif_point_dans_vect(const double* P, vector<double>& V)
-{
-
- int taille=V.size();
- // bool isPresent=false;
- for(int i=0;i<taille/2;i++)
- {
- double dx=P[0]-V[2*i];
- dx=dx>0.0?dx:-dx;
- if (dx>_Precision)
- continue;
- double dy=P[1]-V[2*i+1];
- dy=dy>0.0?dy:-dy;
- if (dy<_Precision)
-
- return;
- // if( sqrt((P[0]-V[2*i])*(P[0]-V[2*i])+(P[1]-V[2*i+1])*(P[1]-V[2*i+1])))/_DimCaracteristic<_Precision)
- // isPresent=true;
-
- }
- // if(!isPresent)
- //{
-
- V.push_back(P[0]);
- V.push_back(P[1]);
-
- //}
-}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/* calcul les coordonnées du barycentre d'un polygone */
-/* le vecteur en entrée est constitué des coordonnées */
-/* des sommets du polygone */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-vector<double> INTERPOLATION_2D::bary_poly(const vector<double>& V)
-{
- vector<double> Bary;
- int taille=V.size();
- double x=0;
- double y=0;
-
- for (int i=0;i<taille/2;i++)
- {
- x=x+V[2*i];
- y=y+V[2*i+1];
- }
- double A=2*x/taille;
- double B=2*y/taille;
- Bary.push_back(A);//taille vecteur=2*nb de points.
- Bary.push_back(B);
-
-
- return Bary;
-}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/* calcul la surface d'un polygone. */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-double INTERPOLATION_2D::Surf_Poly(const vector<double>& Poly)
-{
-
- double Surface=0;
- for (int i=0; i<(Poly.size())/2-2; i++)
- {
- double Surf=Surf_Tri( &Poly[0],&Poly[2*(i+1)],&Poly[2*(i+2)] );
- Surface=Surface + Surf ;
- }
- return Surface ;
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-/* calcul de l'intersection de deux segments: segments P1P2 avec P3P4 */
-/* . Si l'intersection est non nulle et si celle-ci n'est */
-/* n'est pas déjà contenue dans Vect on la rajoute à Vect */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-
-void INTERPOLATION_2D::inters_de_segment(const double * P_1,const double * P_2,
- const double * P_3,const double * P_4,vector<double>& Vect)
-{
-
-
- // calcul du déterminant de P1_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)/_DimCaracteristic>_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]));
-
- 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]));
-
-
- 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);
-
- }
- }
- }
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/* fonction qui teste si un point est dans une maille */
-/* point: P_0 */
-/* P_1, P_2, P_3 sommet des mailles */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-int INTERPOLATION_2D::point_dans_triangle(const double* P_0,const double* P_1,
- const double* P_2,const double* P_3)
-{
-
- int A=0;
- 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) )
- {
- A=1;
- }
-
- return A;
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/* fonction qui rajoute les sommet du triangle P dans le vecteur V */
-/* si ceux-ci sont compris dans le triangle S et ne sont pas déjà dans */
-/* V. */
-/*sommets de P: P_1, P_2, P_3 */
-/*sommets de S: P_4, P_5, P_6 */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-void INTERPOLATION_2D::rajou_sommet_triangl(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>& V)
-{
-
- int A_1=point_dans_triangle(P_1,P_4,P_5,P_6);
- if(A_1==1)
- {verif_point_dans_vect(P_1,V);
- // if (V.size()==1)
-// {
-// // all nodes are necessarily in the triangle
-// verif_point_dans_vect(P_2,V);
-// verif_point_dans_vect(P_3,V);
-// return;
-// }
- }
- int A_2=point_dans_triangle(P_2,P_4,P_5,P_6);
- if(A_2==1)
- {verif_point_dans_vect(P_2,V);}
-
- int A_3=point_dans_triangle(P_3,P_4,P_5,P_6);
- if(A_3==1)
- {verif_point_dans_vect(P_3,V);}
-
-
-
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-/* calcul l'intersection de deux triangles */
-/* P_1, P_2, P_3: sommets du premier triangle */
-/* P_4, P_5, P_6: sommets du deuxième triangle */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-
-vector<double> INTERPOLATION_2D::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)
-{
-
- //debug cout << " T1 : " << *P_1 << " , " << *(P_1+1) << " , " << *P_2 << " , " << *(P_2+1) << " , " << *P_3 << " , " << *(P_3+1)<< endl;
- //debug cout << " T2 : " << *P_4 << " , " << *(P_4+1) << " , " << *P_5 << " , " << *(P_5+1) << " , " << *P_6 << " , " << *(P_6+1)<< endl;
- vector<double> Vect;
-
- inters_de_segment(P_1,P_2,P_4,P_5,Vect);
- inters_de_segment(P_1,P_2,P_5,P_6,Vect);
- inters_de_segment(P_1,P_2,P_6,P_4,Vect);
- inters_de_segment(P_2,P_3,P_4,P_5,Vect);
- inters_de_segment(P_2,P_3,P_5,P_6,Vect);
- inters_de_segment(P_2,P_3,P_6,P_4,Vect);
- inters_de_segment(P_3,P_1,P_4,P_5,Vect);
- inters_de_segment(P_3,P_1,P_5,P_6,Vect);
- inters_de_segment(P_3,P_1,P_6,P_4,Vect);
- rajou_sommet_triangl(P_1,P_2,P_3,P_4,P_5,P_6,Vect);
- rajou_sommet_triangl(P_4,P_5,P_6,P_1,P_2,P_3,Vect);
- //debug cout << " Inter : ";
- //debug for (int i=0; i<Vect.size(); ++i)
- //debug cout << Vect[i] << " ";
- //debug cout << endl << endl;
-
- return Vect;
-}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/* fonction pour reconstituer un polygone convexe à partir */
-/* d'un nuage de point. */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-vector<double> INTERPOLATION_2D::reconstruct_polygon(const vector<double>& V)
-{
-
- int taille=V.size();
- if(taille<=6)
- {return V;}
- else
- {
- double COS[taille/2];
- double SIN[taille/2];
- double angle[taille/2];
- vector<double> Bary=bary_poly(V);
- COS[0]=1.0;
- SIN[0]=0.0;
- angle[0]=0.0;
- for(int i=0; i<taille/2-1;i++)
- {
- vector<double> Trigo=calcul_cos_et_sin(&Bary[0],&V[0],&V[2*(i+1)]);
- COS[i+1]=Trigo[0];
- SIN[i+1]=Trigo[1];
- if(SIN[i+1]>=0)
- {angle[i+1]=acos(COS[i+1]);}
- else
- {angle[i+1]=-acos(COS[i+1]);}
- }
-
- //ensuite on ordonne les angles.
- vector<double> Pt_ordonne;
- Pt_ordonne.reserve(taille);
- map<double,int> Ordre;
- for(int i=0;i<taille/2;i++)
- {Ordre[angle[i]]=i;}
- map <double,int>::iterator mi;
- for(mi=Ordre.begin();mi!=Ordre.end();mi++)
- {
- int j=(*mi).second;
- Pt_ordonne.push_back(V[2*j]);
- Pt_ordonne.push_back(V[2*j+1]);
- }
- return Pt_ordonne;
- }
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ */
-/*fonction pour trouver les mailles voisines d'une maille triangle.*/
-/* (mailles voisines par les faces) */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ */
-
-vector<int> INTERPOLATION_2D::touv_maill_vois(int i_S,const monMaillageS& MaS)
-{
-
- int N_1=MaS._Connect[3*(i_S-1)];
- int N_2=MaS._Connect[3*(i_S-1)+1];
- int N_3=MaS._Connect[3*(i_S-1)+2];
- vector<int> tab;
- int I_1=MaS._ReversNConnectIndex[N_1-1]-1;// Dans MED le n°des noeuds commencent à 1.
- int I_2=MaS._ReversNConnectIndex[N_2-1]-1;
- int I_3=MaS._ReversNConnectIndex[N_3-1]-1;
- int F_1=MaS._ReversNConnectIndex[N_1]-1;
- int F_2=MaS._ReversNConnectIndex[N_2]-1;
- int F_3=MaS._ReversNConnectIndex[N_3]-1;
-
-
- vector<int> V_12;
- V_12.reserve(2);
- insert_iterator<vector<int> > ib_1(V_12,V_12.begin());
- vector<int> V_23;
- V_23.reserve(2);
- insert_iterator<vector<int> > ib_2(V_23,V_23.begin());
- vector<int> V_13;
- V_13.reserve(2);
- insert_iterator<vector<int> > ib_3(V_13,V_13.begin());
-
-
- set_intersection(MaS._ReversNConnect+I_1,MaS._ReversNConnect+F_1,MaS._ReversNConnect+I_2,
- MaS._ReversNConnect+F_2,ib_1);
- set_intersection(MaS._ReversNConnect+I_2,MaS._ReversNConnect+F_2,MaS._ReversNConnect+I_3,
- MaS._ReversNConnect+F_3,ib_2);
- set_intersection(MaS._ReversNConnect+I_1,MaS._ReversNConnect+F_1,MaS._ReversNConnect+I_3,
- MaS._ReversNConnect+F_3,ib_3);
-
-
- for(int i=0;i<(V_12).size();i++)
- {
- if(V_12[i]!=i_S)
- {
- tab.push_back(V_12[i]);
- break;
- }
- }
-
- for(int i=0;i<V_23.size();i++)
- {
- if(V_23[i]!=i_S)
- {
- tab.push_back(V_23[i]);
- break;
- }
- }
-
- for(int i=0;i<V_13.size();i++)
- {
- if(V_13[i]!=i_S)
- {
- tab.push_back(V_13[i]);
- break;
- }
- }
- return tab;
-}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-/* fonction pour vérifier qu'un n°de maille n'a pas déja été considérer */
-/* dans le vecteur et le rajouter au vecteur sinon. */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-
-void INTERPOLATION_2D::verif_maill_dans_vect(int Num, vector<int>& V)
-{
- int taille=V.size();
- int A=0;
- for(int i=0;i<taille;i++)
- {
- if(Num==V[i])
- {
- A=1;
- break;
- }
- }
- if(A==0)
- {V.push_back(Num); }
-}
-
-
-
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
/*localise le barycentre des mailles de P dans le maillage S*/
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-void INTERPOLATION_2D::local_iP_dans_S(MESH& myMesh_S,MESH& myMesh_P,monMaillageP& MaP,
- monMaillageS& MaS,vector<int>& localis )
+void INTERPOLATION_2D::local_iP_dans_S(MESH& myMesh_S,MESH& myMesh_P,ParaMEDMEM_utils::monMaillageP& MaP,
+ ParaMEDMEM_utils::monMaillageS& MaS,vector<int>& localis )
{
//calcul des coordonnées des barycentre des mailles de P
int NP_3=MaP._Connect[3*i_P+2];
double Coord_bary_i_P[2]={Bary[2*i_P],Bary[2*i_P+1]};
- int Vois_N=monInterpol_S.getNearestNode(Coord_bary_i_P);
+ int Vois_N=monInterpol_S.getNearestNode(Coord_bary_i_P)+1;
int Ni=MaS._ReversNConnectIndex[Vois_N-1];
int Nf=MaS._ReversNConnectIndex[Vois_N];
int diff=Nf-Ni;
//debug cout << "mailles sources testées : " << N_Maille << endl;
- vector<double> Inter=intersec_de_triangle(MaS._Coord+2*(NS_1-1),MaS._Coord+2*(NS_2-1),
+ vector<double> Inter=ParaMEDMEM_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));
+ 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);
if(Inter.size()>0)
{
int N_3=MaS._Connect[3*(N_Maille-1)+2];
//debug cout << "mailles sources testées : " << N_Maille << endl;
- vector<double> Inter=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) );
+ vector<double> Inter=ParaMEDMEM_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);
if(Inter.size()>0)
{
/* _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-inline void INTERPOLATION_2D::rempli_vect_d_intersect(int i_P,int i_S1,const monMaillageP& MaP,const monMaillageS& MaS,
+inline void INTERPOLATION_2D::rempli_vect_d_intersect(int i_P,int i_S1,const ParaMEDMEM_utils::monMaillageP& MaP,const ParaMEDMEM_utils::monMaillageS& MaS,
vector<map<int,double> >& Surface_d_intersect,FIELD<double>& myField_P)
{
//on calcule la surface de la maille i_P
- double Surf_i_P =Surf_Tri(MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),
+ double Surf_i_P =ParaMEDMEM_utils::Surf_Tri(MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),
MaP._Coord+2*(NP_3-1));
double Surface = 0 ;
int NS_3=MaS._Connect[3*(i_S-1)+2];
- vector<double> Inter=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));
-
-
-
+ vector<double> Inter=ParaMEDMEM_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);
//on teste l'intersection.
int taille=Inter.size()/2;
//debug cout << " -> maille source : " << i_S << " , nb intersection=" << taille << endl;
/* CAS1 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- if(taille==0)
- {int Rien=0;}
+ // taille ==0, on ne fait rien
/* CAS 2 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- else if(taille==1 || taille==2)
+ if(taille==1 || taille==2)
{
//on ne remplit pas le vecteur des correspondances n° maille-surfaces
//d'intersections mais on récupère le numéro des mailles voisines de la maille i_S.
- vector<int> M_Vois=touv_maill_vois(i_S,MaS);
+ vector<int> M_Vois=ParaMEDMEM_utils::touv_maill_vois(i_S,MaS);
for(int i=0;i< M_Vois.size();i++)
- {verif_maill_dans_vect(M_Vois[i],maill_deja_vu_S);}
+ {ParaMEDMEM_utils::verif_maill_dans_vect(M_Vois[i],maill_deja_vu_S);}
}
/* CAS 3 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
else if(taille==3)
{
- double Surf_inter = Surf_Poly(Inter) ;
+ double Surf_inter = ParaMEDMEM_utils::Surf_Poly(Inter) ;
//on remplit le vecteur des correspondances n° maille-surfaces d'intersections.
- Surface_d_intersect[i_P-1][i_S]=Surf_inter;
-
- // on récupère le numéro des mailles voisines de la maille i_S.
- vector<int> M_Vois=touv_maill_vois(i_S,MaS);
- for(int i_M=0;i_M<(M_Vois).size();i_M++)
- {verif_maill_dans_vect(M_Vois[i_M],maill_deja_vu_S);}
-
- Surface = Surface + Surf_inter ;
+ Surface_d_intersect[i_P-1][i_S]=Surf_inter;
+
+ // on récupère le numéro des mailles voisines de la maille i_S.
+ vector<int> M_Vois=ParaMEDMEM_utils::touv_maill_vois(i_S,MaS);
+ for(int i_M=0;i_M<(M_Vois).size();i_M++)
+ {ParaMEDMEM_utils::verif_maill_dans_vect(M_Vois[i_M],maill_deja_vu_S);}
+
+ Surface = Surface + Surf_inter ;
}
-
+
/* CAS 4 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- else //taille>3
+ else if (taille > 3) //taille>3
{
- vector<double> Poly_Inter=reconstruct_polygon(Inter);
+ vector<double> Poly_Inter=ParaMEDMEM_utils::reconstruct_polygon(Inter);
- double Surf_inter =Surf_Poly(Poly_Inter) ;
+ double Surf_inter =ParaMEDMEM_utils::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;
// on récupère le numéro des mailles voisines de la maille i_S.
- vector<int> M_Vois=touv_maill_vois(i_S,MaS);
+ vector<int> M_Vois=ParaMEDMEM_utils::touv_maill_vois(i_S,MaS);
for(int i_M=0;i_M<(M_Vois).size();i_M++)
- {verif_maill_dans_vect(M_Vois[i_M],maill_deja_vu_S);}
+ {ParaMEDMEM_utils::verif_maill_dans_vect(M_Vois[i_M],maill_deja_vu_S);}
Surface = Surface + Surf_inter ;
}
delete[] Type_P;
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- monMaillageS MaS;
- monMaillageP MaP;
+ ParaMEDMEM_utils::monMaillageS MaS;
+ ParaMEDMEM_utils::monMaillageP MaP;
MaS._NbNoeud = myMesh_S.getNumberOfNodes() ;
MaP._NbNoeud = myMesh_P.getNumberOfNodes() ;
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
- int SpaceDimension_S = myMesh_S.getNumberOfNodes() ;
- int SpaceDimension_P = myMesh_P.getNumberOfNodes() ;
-
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
//on cherche la dimension caractéristique des maillages
vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
// }
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ // for (int i=0; i< Surface_d_intersect.size();i++)
+// {
+// double size=0;
+// for (map<int,double>::const_iterator iter = Surface_d_intersect[i].begin();
+// iter !=Surface_d_intersect[i].end();
+// iter++)
+// size+= iter->second;
+// //debug cout << "\nintersect maille " << i_P << endl;
+// int NP_1=MaP._Connect[3*(i)];
+// int NP_2=MaP._Connect[3*(i)+1];
+// int NP_3=MaP._Connect[3*(i)+2];
+
+
+// //on calcule la surface de la maille i_P
+
+// double Surf_i_P =Surf_Tri(MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),
+// MaP._Coord+2*(NP_3-1));
+// cout <<size<<" "<<Surf_i_P<<endl;
+// }
return Surface_d_intersect;
}
-
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/* fonction qui écrit les résultats d'annelise dans un fichier: */
-/* pour chaque maille du maillage 1 on stoque la fraction volumique */
-/* considéré lors de l'algorithme */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-void WriteFractSurfInFile(MEDMEM::FIELD<double>& FractSurf,MEDMEM::MESH* myMesh_P,
- MEDMEM::MESH* myMesh_S,string NameMesh_P,string NameMesh_S)
-{
- //On sauvegarde le champ crée sur la maillage P.
- string FileName="FracSurf"+NameMesh_S+"_to_"+NameMesh_P;
- int id1=(*myMesh_P).addDriver(MED_DRIVER,FileName,"MP");
- int id2=(*myMesh_S).addDriver(MED_DRIVER,FileName,"MS");
- int id3=FractSurf.addDriver(MED_DRIVER,FileName,"% MP");
- (*myMesh_P).write(id1);
- (*myMesh_S).write(id2);
- FractSurf.write(id3);
-}
+};
#ifndef _INTERPOLATION_2D_HXX_
#define _INTERPOLATION_2D_HXX_
-
#include "MEDMEM_InterpolationHighLevelObjects.hxx"
#include "MEDMEM_Mesh.hxx"
#include "MEDMEM_Field.hxx"
#include "MEDMEM_DriverFactory.hxx"
#include "MEDMEM_Meshing.hxx"
#include "MEDMEM_define.hxx"
+#include "INTERPOLATION.hxx"
#include "MEDMEM_Interpolation.hxx"
#include <complex>
-struct monMaillageS;
-
+namespace ParaMEDMEM_utils
+{
+ struct monMaillageS;
+ struct monMaillageP;
+};
-struct monMaillageP;
+namespace ParaMEDMEM
+{
+ class INTERPOLATION;
-class INTERPOLATION_2D
+class INTERPOLATION_2D : public INTERPOLATION
{
-
private:
double _Precision;
double _DimCaracteristic;
int _SecondMethodOfLocalisation;
int _PrintLevel;
+
// Méthodes publiques
public:
INTERPOLATION_2D();
-
// precision geometrique, choix de la methode comlementaire de localisation, niveau d'impressions
void setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel);
- //pour calculer la surface d'un triangle.
- double Surf_Tri(const double* P_1,const double* P_2,
- const double* P_3);
-
- //calcul déterminant de trois points (cf doc CGAL)
- double mon_determinant(const double* P_1,const double* P_2,
- const double* P_3);
-
- //calcul la norme d'un vecteur.
- double norme_vecteur(const double* P_1,const double* P_2);
-
- //calcul le cos et le sin.
- vector<double> calcul_cos_et_sin(const double* P_1,
- const double* P_2,const double* P_3);
-
- //fonction pour vérifier qu'un point n'a pas déja été considérer dans
- //le vecteur et le rajouter au vecteur sinon.
- void verif_point_dans_vect(const double* P, vector<double>& V);
-
- //calcul les coordonnées du barycentre
- //d'un polygone
- vector<double> bary_poly(const vector<double>& V);
-
- //fonction qui calcul la surface d'un polygone
- double Surf_Poly(const vector<double>& Poly);
-
- //calcul de l'intersection de deux segments.
- void inters_de_segment(const double* P_1,const double* P_2,
- const double* P_3,const double* P_4,vector<double>& Vect);
-
- //fonction qui teste si un point est dans une maille
- int point_dans_triangle(const double* P_0,const double* P_1,
- const double* P_2,const double* P_3);
-
- //fonction qui rajouter les sommets du triangle P2P2P3 si ceux-ci sont compris
- // dans le triangle P4P5P6 et n'ont pas encore été considérer.
- //Rq:point du triangle donné dans ordre quelconque.
- void rajou_sommet_triangl(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>& V);
-
- //calcul l'intersection de deux triangles (sommets du triangle donné dans
- //ordre quelconque)
- 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);
-
- //fonction pour reconstituer un polygon convexe à partir
- //d'un nuage de point.
- vector<double> reconstruct_polygon(const vector<double>& V);
-
- //fonction pour trouver les mailles voisines d'une maille triangle.
- vector<int> touv_maill_vois(int i_S,const monMaillageS& MaS);
-
- //fonction pour vérifier qu'un n°de maille n'a pas déja été considérer
- // dans le vecteur et le rajouter au vecteur sinon.
- void INTERPOLATION_2D::verif_maill_dans_vect(int Num, vector<int>& V);
-
//localise le barycentre des mailles de P dans le maillage S
void local_iP_dans_S(MEDMEM::MESH& myMesh_S,MEDMEM::MESH& myMesh_P,
- monMaillageP& MaP,monMaillageS& MaS,vector<int>& localis );
+ ParaMEDMEM_utils::monMaillageP& MaP,ParaMEDMEM_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 monMaillageP& MaP,const monMaillageS& MaS,
+ inline void rempli_vect_d_intersect(int i_P,int i_S,const ParaMEDMEM_utils::monMaillageP& MaP,const ParaMEDMEM_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);
-
- // fonction qui écrit les résultats d'annelise dans un fichier:
- // pour chaque maille du maillage 1 on stoque la fraction volumique
- // considéré lors de l'algorithme
-
-
-
-
-
- void WriteFractSurfInFile(MEDMEM::FIELD<double>& FractSurf,MEDMEM::MESH* myMesh_P,
- MEDMEM::MESH* myMesh_S,string FileName,string NameMesh_P,string NameMesh_S);
-
-
-
- // MEDMEM::FIELD<double>* createField();
- // ...
-
private:
+};
+
};
#endif
--- /dev/null
+#include "MEDMEM_Field.hxx"
+#include "MEDMEM_Mesh.hxx"
+#include "TranslationRotationMatrix.hxx"
+#include "INTERPOLATION_Utils.hxx"
+#include "INTERPOLATION_3D_surf.hxx"
+#include "BBTree.H"
+
+using namespace std;
+using namespace MED_EN;
+using namespace MEDMEM;
+using namespace ParaMEDMEM_utils;
+
+namespace ParaMEDMEM
+{
+
+ INTERPOLATION_3D_surf::INTERPOLATION_3D_surf()
+ {
+ _Precision=1.0E-12;
+ _DimCaracteristic=1;
+ _SecondMethodOfLocalisation=1;
+ _PrintLevel=0;
+ _Surf3DAdjustmentEps=0.1;
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* Options : */
+ /* Precision : for geometric computation */
+ /* SecondMethodOfLocalisation : 0/1 */
+ /* PrintLevel : between 0 and 3 */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ void INTERPOLATION_3D_surf::setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel)
+ {
+ _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 INTERPOLATION_3D_surf::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 INTERPOLATION_3D_surf::rempli_vect_d_intersect(int i_P,int i_S,const monMaillageP& MaP,const monMaillageS& MaS,
+ vector<map<int,double> >& Surface_d_intersect
+ )
+
+ {
+
+
+ //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 =ParaMEDMEM_utils::Surf_Tri(MaP._Coord+3*(NP_1-1),MaP._Coord+3*(NP_2-1),
+ // MaP._Coord+3*(NP_3-1));
+
+ // double Surface = 0 ;
+
+
+
+ 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=ParaMEDMEM_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;
+
+
+ /* CAS 3 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ if(taille==3)
+ {
+ 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 ;
+ }
+
+ /* CAS 4 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+ else if (taille>3) //taille>3
+ {
+ 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 ;
+ }
+
+ //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> > INTERPOLATION_3D_surf::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;
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ ParaMEDMEM_utils::monMaillageS MaS;
+ ParaMEDMEM_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
+
+ BBTree<3> my_tree(&bbox[0], 0, 0,MaS._NbMaille );//creating the search structure
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ //on crée un tableau où l'on rentrera la valeurs des intersections.
+ cout << "INTERPOLATION_3D_surf::calcul intersec"<<endl;
+ map<int,double> MAP_init;
+ vector<map<int,double> > Surface_d_intersect(MaP._NbMaille,MAP_init);//on initialise.
+
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ //boucle sur les mailles de P.
+ //Coeur de l'algorithme
+
+ for(int i_P=0; i_P<MaP._NbMaille ; i_P++)
+ {
+ vector<int> intersecting_elems;
+ double bb[6];
+ ParaMEDMEM_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++)
+ {
+ //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);
+ }
+ intersecting_elems.clear();
+ }
+
+
+ //DEBUG prints
+ if (_PrintLevel >= 2)
+ {
+ 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)
+ {
+ map<int,double>::iterator surface;
+ for( surface=Surface_d_intersect[i].begin();surface!=Surface_d_intersect[i].end();surface++)
+ {
+ cout<<" ("<<i+1<<" , "<<(*surface).first<<")"<<" : "<<(*surface).second<<endl;
+ total+=(*surface).second;
+ }
+ }
+ }
+ 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 INTERPOLATION_3D_surf::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;
+ }
+ }
+ }
+
+
+ /*! 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 INTERPOLATION_3D_surf::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;
+
+ }
+ }
+
+};
--- /dev/null
+#ifndef _INTERPOLATION_3D_SURF_HXX_
+#define _INTERPOLATION_3D_SURF_HXX_
+
+
+
+#include "MEDMEM_Mesh.hxx"
+#include "INTERPOLATION_Utils.hxx"
+#include "TranslationRotationMatrix.hxx"
+#include "INTERPOLATION.hxx"
+
+#include <cmath>
+#include <map>
+#include <sstream>
+#include <string>
+#include <set>
+#include <algorithm>
+#include <vector>
+
+namespace ParaMEDMEM
+{
+
+ class INTERPOLATION_3D_surf : public INTERPOLATION
+{
+
+
+ private:
+ double _Precision;
+ double _DimCaracteristic;
+ int _SecondMethodOfLocalisation;
+ int _PrintLevel;
+ double _Surf3DAdjustmentEps;
+ // Méthodes publiques
+ public:
+
+ INTERPOLATION_3D_surf();
+
+
+ // 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);
+
+
+ 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 ParaMEDMEM_utils::monMaillageP& MaP,
+ const ParaMEDMEM_utils::monMaillageS& MaS,
+ vector<map<int,double> >& Surface_d_intersect);
+
+
+ void createBoundingBoxes(const struct ParaMEDMEM_utils::monMaillageS mesh, vector<double>& bbox);
+
+ void rotate3DTriangle(const double* PP1,
+ const double* PP2,
+ const double* PP3,
+ TranslationRotationMatrix& matrix);
+
+ void adjustBoundingBoxes(vector<double>& bbox);
+
+
+};
+
+};
+
+#endif
--- /dev/null
+#ifndef _INTERPOLATION_UTILS_HXX_
+#define _INTERPOLATION_UTILS_HXX_
+
+
+#include "MEDMEM_Mesh.hxx"
+#include "MEDMEM_Field.hxx"
+
+
+#include <cmath>
+#include <map>
+#include <algorithm>
+#include <vector>
+
+
+namespace ParaMEDMEM_utils
+{
+ // const double HUGE=1e300;
+ struct monMaillageS
+ {
+ int _NbMaille;
+ int _NbNoeud;
+ const int* _Connect;
+ const double* _Coord;
+ const int* _ReversNConnect;
+ const int* _ReversNConnectIndex;
+ double _DimCaracteristic;
+ };
+
+ struct monMaillageP
+ {
+ int _NbMaille;
+ int _NbNoeud;
+ const int* _Connect;
+ const double* _Coord;
+ double _DimCaracteristic;
+ };
+
+// void verif_point_dans_vect(const double* P, vector<double>& V,double dim_caracteristic, double precision );
+// inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3);
+// inline double norme_vecteur(const double* P_1,const double* P_2);
+// inline double mon_determinant(const double* P_1,
+// const double* P_2,
+// const double* P_3);
+// vector<double> calcul_cos_et_sin(const double* P_1,
+// const double* P_2,
+// const double* P_3);
+// inline vector<double> bary_poly(const vector<double>& V);
+// bool point_dans_triangle(const double* P_0,const double* P_1,
+// const double* P_2,const double* P_3);
+// double Surf_Poly(const vector<double>& Poly);
+// vector<double> reconstruct_polygon(const vector<double>& V);
+
+// void rajou_sommet_triangl(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>& V,
+// double dim_caracteristic,
+// double precision);
+
+
+// inline void inters_de_segment(const double * P_1,
+// const double * P_2,
+// const double * P_3,
+// const double * P_4,
+// vector<double>& Vect,
+// double dim_caracteristic,
+// double precision);
+// 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);
+// vector<int> touv_maill_vois(int i_S,const monMaillageS& MaS);
+// void verif_point_dans_vect(const double* P, vector<double>& V,double dim_caracteristic, double precision );
+// void WriteFractSurfInFile(MEDMEM::FIELD<double>& FractSurf,MEDMEM::MESH* myMesh_P,
+// MEDMEM::MESH* myMesh_S,string NameMesh_P,string NameMesh_S);
+
+
+
+
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* calcul la surface d'un triangle */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+ 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);
+ return Surface;
+ }
+
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* fonction qui calcul le déterminant */
+ /* de deux vecteur(cf doc CGAL). */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+
+ //fonction qui calcul le déterminant des vecteurs: P3P1 et P3P2
+ //(cf doc CGAL).
+
+ inline double mon_determinant(const double* P_1,
+ const double* P_2,
+ const double* P_3)
+ {
+ double mon_det=(P_1[0]-P_3[0])*(P_2[1]-P_3[1])-(P_2[0]-P_3[0])*(P_1[1]-P_3[1]);
+ return mon_det;
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ //calcul la norme du vecteur P1P2
+
+ inline double norme_vecteur(const double* P_1,const double* P_2)
+ {
+ double X=P_1[0]-P_2[0];
+ double Y=P_1[1]-P_2[1];
+ double norme=sqrt(X*X+Y*Y);
+ return norme;
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* calcul le cos et le sin de l'angle P1P2,P1P3 */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+ inline vector<double> calcul_cos_et_sin(const double* P_1,
+ const double* P_2,
+ const double* P_3)
+ {
+
+ vector<double> Vect;
+ double P1_P2=norme_vecteur(P_1,P_2);
+ double P2_P3=norme_vecteur(P_2,P_3);
+ double P3_P1=norme_vecteur(P_3,P_1);
+
+ double N=P1_P2*P1_P2+P3_P1*P3_P1-P2_P3*P2_P3;
+ double D=2.0*P1_P2*P3_P1;
+ double COS=N/D;
+ Vect.push_back(COS);
+ double V=mon_determinant(P_2,P_3,P_1);
+ double D_1=P1_P2*P3_P1;
+ double SIN=V/D_1;
+ Vect.push_back(SIN);
+
+ return Vect;
+
+ }
+
+
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* calcul les coordonnées du barycentre d'un polygone */
+ /* le vecteur en entrée est constitué des coordonnées */
+ /* des sommets du polygone */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+ inline vector<double> bary_poly(const vector<double>& V)
+ {
+ vector<double> Bary;
+ int taille=V.size();
+ double x=0;
+ double y=0;
+
+ for (int i=0;i<taille/2;i++)
+ {
+ x=x+V[2*i];
+ y=y+V[2*i+1];
+ }
+ double A=2*x/taille;
+ double B=2*y/taille;
+ Bary.push_back(A);//taille vecteur=2*nb de points.
+ Bary.push_back(B);
+
+
+ return Bary;
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* calcul la surface d'un polygone. */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+ inline double Surf_Poly(const vector<double>& Poly)
+ {
+
+ double Surface=0;
+ for (int i=0; i<(Poly.size())/2-2; i++)
+ {
+ double Surf=Surf_Tri( &Poly[0],&Poly[2*(i+1)],&Poly[2*(i+2)] );
+ Surface=Surface + Surf ;
+ }
+ return Surface ;
+ }
+
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* fonction qui teste si un point est dans une maille */
+ /* point: P_0 */
+ /* P_1, P_2, P_3 sommet des mailles */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+ inline bool point_dans_triangle(const double* P_0,const double* P_1,
+ const double* P_2,const double* P_3)
+ {
+
+ 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) )
+ {
+ A=true;
+ }
+
+ return A;
+ }
+
+
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /*fonction pour vérifier qu'un point n'a pas déja été considérer dans */
+ /* 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 )
+ {
+ // int taille=V.size();
+ // for(int i=0;i<taille/2;i++)
+ // {
+ // double dx=P[0]-V[2*i];
+ // dx=dx>0.0?dx:-dx;
+ // if (dx>_Precision)
+ // continue;
+ // double dy=P[1]-V[2*i+1];
+ // dy=dy>0.0?dy:-dy;
+ // if (dy<_Precision)
+
+ // return;
+
+ // }
+ // V.push_back(P[0]);
+ // V.push_back(P[1]);
+
+ int taille=V.size();
+ 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)
+ isPresent=true;
+
+ }
+ if(!isPresent)
+ {
+
+ V.push_back(P[0]);
+ V.push_back(P[1]);
+ }
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* fonction qui rajoute les sommet du triangle P dans le vecteur V */
+ /* si ceux-ci sont compris dans le triangle S et ne sont pas déjà dans */
+ /* V. */
+ /*sommets de P: P_1, P_2, P_3 */
+ /*sommets de S: P_4, P_5, P_6 */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+ inline void rajou_sommet_triangl(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>& V, double dim_caracteristic, double precision)
+ {
+
+ bool A_1=ParaMEDMEM_utils::point_dans_triangle(P_1,P_4,P_5,P_6);
+ if(A_1)
+ {verif_point_dans_vect(P_1,V,dim_caracteristic, precision);
+ // if (V.size()==1)
+ // {
+ // // all nodes are necessarily in the triangle
+ // verif_point_dans_vect(P_2,V);
+ // verif_point_dans_vect(P_3,V);
+ // return;
+ // }
+ }
+ bool A_2=ParaMEDMEM_utils::point_dans_triangle(P_2,P_4,P_5,P_6);
+ if(A_2)
+ {verif_point_dans_vect(P_2,V,dim_caracteristic,precision);}
+
+ bool A_3=ParaMEDMEM_utils::point_dans_triangle(P_3,P_4,P_5,P_6);
+ if(A_3)
+ {verif_point_dans_vect(P_3,V,dim_caracteristic, precision);}
+
+ }
+
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /* calcul de l'intersection de deux segments: segments P1P2 avec P3P4 */
+ /* . Si l'intersection est non nulle et si celle-ci n'est */
+ /* n'est pas déjà contenue dans Vect on la rajoute à Vect */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+
+ inline void inters_de_segment(const double * P_1,const double * P_2,
+ const double * P_3,const double * P_4,vector<double>& Vect,
+ double dim_caracteristic, double precision)
+ {
+
+
+ // calcul du déterminant de P1_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 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]));
+
+ 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]));
+
+
+ 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);
+
+ }
+ }
+ }
+ }
+
+
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /* calcul l'intersection de deux triangles */
+ /* P_1, P_2, P_3: sommets du premier triangle */
+ /* 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)
+ {
+
+ // 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;
+
+ 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);
+ inters_de_segment(P_1,P_2,P_6,P_4,Vect, dim_caracteristic, precision);
+ inters_de_segment(P_2,P_3,P_4,P_5,Vect, dim_caracteristic, precision);
+ inters_de_segment(P_2,P_3,P_5,P_6,Vect, dim_caracteristic, precision);
+ inters_de_segment(P_2,P_3,P_6,P_4,Vect, dim_caracteristic, precision);
+ inters_de_segment(P_3,P_1,P_4,P_5,Vect, dim_caracteristic, precision);
+ inters_de_segment(P_3,P_1,P_5,P_6,Vect, dim_caracteristic, precision);
+ inters_de_segment(P_3,P_1,P_6,P_4,Vect, dim_caracteristic, precision);
+ rajou_sommet_triangl(P_1,P_2,P_3,P_4,P_5,P_6,Vect, dim_caracteristic, precision);
+ rajou_sommet_triangl(P_4,P_5,P_6,P_1,P_2,P_3,Vect, dim_caracteristic, precision);
+ //debug cout << " Inter : ";
+ //debug for (int i=0; i<Vect.size(); ++i)
+ //debug cout << Vect[i] << " ";
+ //debug cout << endl << endl;
+
+ return Vect;
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ */
+ /*fonction pour trouver les mailles voisines d'une maille triangle.*/
+ /* (mailles voisines par les faces) */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ */
+
+ inline vector<int> touv_maill_vois(int i_S,const monMaillageS& MaS)
+ {
+
+ int N_1=MaS._Connect[3*(i_S-1)];
+ int N_2=MaS._Connect[3*(i_S-1)+1];
+ int N_3=MaS._Connect[3*(i_S-1)+2];
+ vector<int> tab;
+ int I_1=MaS._ReversNConnectIndex[N_1-1]-1;// Dans MED le n°des noeuds commencent à 1.
+ int I_2=MaS._ReversNConnectIndex[N_2-1]-1;
+ int I_3=MaS._ReversNConnectIndex[N_3-1]-1;
+ int F_1=MaS._ReversNConnectIndex[N_1]-1;
+ int F_2=MaS._ReversNConnectIndex[N_2]-1;
+ int F_3=MaS._ReversNConnectIndex[N_3]-1;
+
+
+ vector<int> V_12;
+ V_12.reserve(2);
+ insert_iterator<vector<int> > ib_1(V_12,V_12.begin());
+ vector<int> V_23;
+ V_23.reserve(2);
+ insert_iterator<vector<int> > ib_2(V_23,V_23.begin());
+ vector<int> V_13;
+ V_13.reserve(2);
+ insert_iterator<vector<int> > ib_3(V_13,V_13.begin());
+
+
+ set_intersection(MaS._ReversNConnect+I_1,MaS._ReversNConnect+F_1,MaS._ReversNConnect+I_2,
+ MaS._ReversNConnect+F_2,ib_1);
+ set_intersection(MaS._ReversNConnect+I_2,MaS._ReversNConnect+F_2,MaS._ReversNConnect+I_3,
+ MaS._ReversNConnect+F_3,ib_2);
+ set_intersection(MaS._ReversNConnect+I_1,MaS._ReversNConnect+F_1,MaS._ReversNConnect+I_3,
+ MaS._ReversNConnect+F_3,ib_3);
+
+
+ for(int i=0;i<(V_12).size();i++)
+ {
+ if(V_12[i]!=i_S)
+ {
+ tab.push_back(V_12[i]);
+ break;
+ }
+ }
+
+ for(int i=0;i<V_23.size();i++)
+ {
+ if(V_23[i]!=i_S)
+ {
+ tab.push_back(V_23[i]);
+ break;
+ }
+ }
+
+ for(int i=0;i<V_13.size();i++)
+ {
+ if(V_13[i]!=i_S)
+ {
+ tab.push_back(V_13[i]);
+ break;
+ }
+ }
+ return tab;
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /* fonction pour vérifier qu'un n°de maille n'a pas déja été considérer */
+ /* dans le vecteur et le rajouter au vecteur sinon. */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+
+ inline void verif_maill_dans_vect(int Num, vector<int>& V)
+ {
+ int taille=V.size();
+ int A=0;
+ for(int i=0;i<taille;i++)
+ {
+ if(Num==V[i])
+ {
+ A=1;
+ break;
+ }
+ }
+ if(A==0)
+ {V.push_back(Num); }
+ }
+
+
+
+
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* fonction pour reconstituer un polygone convexe à partir */
+ /* d'un nuage de point. */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+ inline vector<double> reconstruct_polygon(const vector<double>& V)
+ {
+
+ int taille=V.size();
+
+ //VB : why 6 ?
+
+ if(taille<=6)
+ {return V;}
+ else
+ {
+ double COS[taille/2];
+ double SIN[taille/2];
+ double angle[taille/2];
+ vector<double> Bary=bary_poly(V);
+ COS[0]=1.0;
+ SIN[0]=0.0;
+ angle[0]=0.0;
+ for(int i=0; i<taille/2-1;i++)
+ {
+ vector<double> Trigo=calcul_cos_et_sin(&Bary[0],&V[0],&V[2*(i+1)]);
+ COS[i+1]=Trigo[0];
+ SIN[i+1]=Trigo[1];
+ if(SIN[i+1]>=0)
+ {angle[i+1]=acos(COS[i+1]);}
+ else
+ {angle[i+1]=-acos(COS[i+1]);}
+ }
+
+ //ensuite on ordonne les angles.
+ vector<double> Pt_ordonne;
+ Pt_ordonne.reserve(taille);
+ multimap<double,int> Ordre;
+ for(int i=0;i<taille/2;i++)
+ {Ordre.insert(make_pair(angle[i],i));}
+ multimap <double,int>::iterator mi;
+ for(mi=Ordre.begin();mi!=Ordre.end();mi++)
+ {
+ int j=(*mi).second;
+ Pt_ordonne.push_back(V[2*j]);
+ Pt_ordonne.push_back(V[2*j+1]);
+ }
+ return Pt_ordonne;
+ }
+ }
+
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* fonction qui écrit les résultats d'annelise dans un fichier: */
+ /* pour chaque maille du maillage 1 on stoque la fraction volumique */
+ /* considéré lors de l'algorithme */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+ inline void WriteFractSurfInFile(MEDMEM::FIELD<double>& FractSurf,MEDMEM::MESH* myMesh_P,
+ MEDMEM::MESH* myMesh_S,string NameMesh_P,string NameMesh_S)
+ {
+ //On sauvegarde le champ crée sur la maillage P.
+ string FileName="FracSurf"+NameMesh_S+"_to_"+NameMesh_P;
+ int id1=(*myMesh_P).addDriver(MED_DRIVER,FileName,"MP");
+ int id2=(*myMesh_S).addDriver(MED_DRIVER,FileName,"MS");
+ int id3=FractSurf.addDriver(MED_DRIVER,FileName,"% MP");
+ (*myMesh_P).write(id1);
+ (*myMesh_S).write(id2);
+ FractSurf.write(id3);
+ }
+
+ inline void getElemBB(double* bb, const monMaillageP& MaP, int iP)
+ {
+ 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++)
+ {
+ 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];
+ 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];
+ }
+
+ }
+
+};
+#endif
#include "MxN_Mapping.hxx"
#include "InterpolationMatrix.hxx"
#include "INTERPOLATION_2D.hxx"
+#include "INTERPOLATION_3D_surf.hxx"
/*! \class InterpolationMatrix
This class enables the storage of an interpolation matrix Wij mapping
*/
void InterpolationMatrix::addContribution(MEDMEM::MESH& distant_support, int iproc_distant, int* distant_elems)
{
- INTERPOLATION_2D interpolator;
- int nbelems= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
- //_row_offsets.resize(nbelems);
- vector<map<int,double> > surfaces = interpolator.interpol_maillages(distant_support,_source_support);
+
+ // computing the intersections between distant and local elements
+ INTERPOLATION* interpolator;
+ if (distant_support.getMeshDimension()==2 && distant_support.getSpaceDimension()==3)
+ interpolator=new INTERPOLATION_3D_surf();
+ else if (distant_support.getMeshDimension()==2 && distant_support.getSpaceDimension()==2)
+ interpolator=new INTERPOLATION_2D();
+
+ int nbelems= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ vector<map<int,double> > surfaces = interpolator->interpol_maillages(distant_support,_source_support);
+ delete interpolator;
+
+
if (surfaces.size() != nbelems)
throw MEDEXCEPTION("uncoherent number of rows in interpolation matrix");
iter!=surfaces[ielem].end();
iter++)
{
+ //surface of the source triangle
double surf = triangle_surf->getValueIJ(iter->first,1);
+ //oriented surfaces are not considered
surf = (surf>0)?surf:-surf;
+ //locating the (iproc, itriangle) pair in the list of columns
vector<pair<int,int> >::iterator iter2 = find (_col_offsets.begin(), _col_offsets.end(),make_pair(iproc_distant,iter->first));
int col_id;
+
+ //(iproc, itriangle) is not registered in the list
+ //of distant elements
if (iter2==_col_offsets.end())
{
_col_offsets.push_back(make_pair(iproc_distant,iter->first));
col_id =_col_offsets.size();
- _mapping.addElementFromSource(col_id,iproc_distant,distant_elems[iter->first-1]);
+ _mapping.addElementFromSource(iproc_distant,distant_elems[iter->first-1]);
}
- else
+ else
{
col_id=iter2-_col_offsets.begin()+1;
}
+ //the column indirection number is registered
+ //with the interpolation coefficient
+ //actual column number is obtained with _col_offsets[col_id]
_col_numbers.push_back(col_id);
_coeffs.push_back(iter->second/surf);
-
}
}
delete triangle_surf;
{
target_value[i]=0.0;
}
- // int nbrows = _source_indices.size();
int nbrows= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
for (int irow=0; irow<nbrows; irow++)
{
int* distant_ids=0;
for (int i=0; i<_target_group->size(); i++)
{
- int idistant_proc = (i+_source_group->myRank())%_target_group->size();
- // idistant_proc=i;
+ // int idistant_proc = (i+_source_group->myRank())%_target_group->size();
+ int idistant_proc=i;
//gathers pieces of the target meshes that can intersect the local mesh
locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
int* distant_ids=0;
for (int i=0; i<_source_group->size(); i++)
{
- int idistant_proc = (i+_target_group->myRank())%_source_group->size();
-
+ // int idistant_proc = (i+_target_group->myRank())%_source_group->size();
+ int idistant_proc=i;
//gathers pieces of the target meshes that can intersect the local mesh
locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
cout << " Data sent from "<<_union_group->myRank()<<" to source proc "<< idistant_proc<<endl;
DEC.hxx\
MxN_Mapping.hxx\
INTERPOLATION_2D.hxx\
+INTERPOLATION_3D_surf.hxx\
+INTERPOLATION.hxx\
+INTERPOLATION_Utils.hxx\
StructuredCoincidentDEC.hxx\
InterpolationMatrix.hxx\
IntersectionDEC.hxx\
UnstructuredParaSUPPORT.hxx\
ExplicitCoincidentDEC.hxx\
ElementLocator.hxx\
-ExplicitMapping.hxx
+ExplicitMapping.hxx\
+TranslationRotationMatrix.hxx\
+BBTree.H
# Libraries targets
ParaFIELD.cxx\
DEC.cxx\
INTERPOLATION_2D.cxx\
+INTERPOLATION_3D_surf.cxx\
MxN_Mapping.cxx\
InterpolationMatrix.cxx\
StructuredCoincidentDEC.cxx\
BIN_CLIENT_IDL =
TEST_PROGS = test_ProcessorGroup test_BlockTopology test_ParaStructuredSupport \
-test_ParaField test_DEC test_UnstructuredDEC test_ExplicitDEC test_IntersectionDEC
+test_ParaField test_DEC test_UnstructuredDEC test_ExplicitDEC test_IntersectionDEC test_SeqIntersectionDEC test_Seq3DsurfIntersection test_MEDMEMConstructor make_cylinder make_plane test_INTERPOL_2D
LDFLAGS+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome
LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome
delete _union_group;
}
-void MxN_Mapping::addElementFromSource(int local_element, int distant_proc, int distant_element)
+
+/*!
+Method registering a new element for correspondence with a distant element
+\param distant_proc proc rank of the distant processor (in terms of the union group)
+\param distant_element id of the element on the distant processor
+ */
+void MxN_Mapping::addElementFromSource(int distant_proc, int distant_element)
{
_sending_ids.push_back(make_pair(distant_proc,distant_element));
- // _local_ids.push_back(local_element);
for (int i=distant_proc; i<_union_group->size(); i++)
_send_proc_offsets[i+1]++;
}
MxN_Mapping();
MxN_Mapping(const ProcessorGroup& local_group, const ProcessorGroup& distant_group);
virtual ~MxN_Mapping();
- void addElementFromSource(int local_elem, int distant_proc, int distant_elem);
+ void addElementFromSource(int distant_proc, int distant_elem);
void prepareSendRecv();
void sendRecv(MEDMEM::FIELD<double>& field);
void sendRecv(double* field, MEDMEM::FIELD<double>& field) const ;
_field->setNumberOfValues(para_support->getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS));
string* compnames=new string[nb_components];
string* compdesc=new string[nb_components];
+ string* compunit=new string[nb_components];
for (int i=0; i<nb_components; i++)
{
ostringstream stream(compnames[i]);
ostringstream stream2(compdesc[i]);
+ ostringstream stream3(compunit[i]);
stream<<"component "<<i;
stream2<<"component description "<<i;
+ stream3<<"component unit "<<i;
+
+ compnames[i]=stream.str();
+ compdesc[i]=stream2.str();
+ compunit[i]=stream3.str();
+
}
_field->setComponentsNames(compnames);
_field->setComponentsDescriptions(compdesc);
+ _field->setMEDComponentsUnits(compunit);
_field->setIterationNumber(0);
_field->setOrderNumber(0);
_field->setTime(0.0);
ParaMESH::ParaMESH(MEDMEM::MESH& subdomain_mesh, const ProcessorGroup& proc_group, const string& name):
_mesh(&subdomain_mesh),
-_name (name),
_my_domain_id(proc_group.myRank()),
_block_topology (new BlockTopology(proc_group, subdomain_mesh.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS)))
{
+ ostringstream stream;
+ stream<<name<<"_"<<_my_domain_id+1;
+ _name=stream.str();
+
_cellglobal = new int[subdomain_mesh.getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS)];
int offset = _block_topology->localToGlobal(make_pair(_my_domain_id,0));
for (int i=0; i<subdomain_mesh.getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS); i++)
--- /dev/null
+#ifndef _TRANSLATIONROTATIONMATRIX_HXX_
+#define _TRANSLATIONROTATIONMATRIX_HXX_
+
+namespace ParaMEDMEM
+{
+
+ const double EPS=1e-12;
+
+class TranslationRotationMatrix
+{
+
+public:
+
+ TranslationRotationMatrix():_rotation_coeffs(9,0.0), _translation_coeffs(3,0.0)
+ {
+ _rotation_coeffs[0]=1.0;
+ _rotation_coeffs[4]=1.0;
+ _rotation_coeffs[8]=1.0;
+ }
+
+ void multiply(const TranslationRotationMatrix& A)
+ {
+ TranslationRotationMatrix result;
+ //setting output matrix to zero
+ for (int i=0; i<3; i++)
+ result._rotation_coeffs[i*4]=0.0;
+ //multiplying
+ for (int i=0; i<3;i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<3; k++)
+ result._rotation_coeffs[j+i*3]+=A._rotation_coeffs[3*i+k]*_rotation_coeffs[j+k*3];
+
+ for (int i=0;i<9; i++)
+ _rotation_coeffs[i]=result._rotation_coeffs[i];
+
+ // cout << "Matrix "<<endl;
+// for (int i=0; i<3; i++)
+// {
+// cout<<"| ";
+// for (int j=0; j<3;j++)
+// cout <<_rotation_coeffs[i*3+j]<<" ";
+// cout<<" "<<endl;
+// }
+ }
+
+ void rotate_vector(double* P)
+ {
+ double temp[3]={0.0, 0.0, 0.0};
+
+ for (int i=0; i<3;i++)
+ for (int j=0; j<3; j++)
+ temp[i] +=_rotation_coeffs[3*i+j]*P[j];
+
+ P[0]=temp[0];P[1]=temp[1];P[2]=temp[2];
+ }
+
+ void transform_vector(double*P)
+ {
+ P[0]+=_translation_coeffs[0];
+ P[1]+=_translation_coeffs[1];
+ P[2]+=_translation_coeffs[2];
+ rotate_vector(P);
+ }
+ void translate(const double* P)
+ {
+ _translation_coeffs[0]=P[0];
+ _translation_coeffs[1]=P[1];
+ _translation_coeffs[2]=P[2];
+ }
+
+ void rotate_x (double* P)
+ {
+ _rotation_coeffs[0]=1.0;
+ double r_sqr = P[1]*P[1]+P[2]*P[2];
+ if (r_sqr < EPS)
+ {_rotation_coeffs[4]=1.0; _rotation_coeffs[8]=1.0; return;}
+ double r = sqrt(r_sqr);
+ double cos =P[1]/r;
+ double sin =P[2]/r;
+
+ _rotation_coeffs[4]=cos;
+ _rotation_coeffs[5]=sin;
+ _rotation_coeffs[7]=-sin;
+ _rotation_coeffs[8]=cos;
+
+
+ rotate_vector(P);
+ }
+ void rotate_z (double* P)
+ {
+ _rotation_coeffs[8]=1.0;
+ double r_sqr = P[0]*P[0]+P[1]*P[1];
+ if (r_sqr < EPS)
+ {_rotation_coeffs[4]=1.0; _rotation_coeffs[0]=1.0; return;}
+ double r = sqrt(r_sqr);
+ double cos =P[0]/r;
+ double sin =P[1]/r;
+
+ _rotation_coeffs[0]=cos;
+ _rotation_coeffs[1]=sin;
+ _rotation_coeffs[3]=-sin;
+ _rotation_coeffs[4]=cos;
+
+ rotate_vector(P);
+ }
+
+
+private:
+ vector<double> _rotation_coeffs;
+ vector<double> _translation_coeffs;
+
+};
+
+};
+#endif
