src/INTERP_KERNEL/PlanarIntersector.txx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 #ifndef __PLANARINTERSECTOR_TXX__
  20 #define __PLANARINTERSECTOR_TXX__
  21
  22 #include "PlanarIntersector.hxx"
  23 #include "InterpolationUtils.hxx"
  24 #include "TranslationRotationMatrix.hxx"
  25
  26 #include <iostream>
  27 #include <limits>
  28
  29 namespace INTERP_KERNEL
  30 {
  31   template<class MyMeshType, class MyMatrix>
  32   PlanarIntersector<MyMeshType,MyMatrix>::PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel):
  33     _meshT(meshT),_meshS(meshS),
  34     _dim_caracteristic(dimCaracteristic),_max_distance_3Dsurf_intersect(md3DSurf),_precision(precision),_median_plane(medianPlane),
  35     _do_rotate(doRotate),_orientation(orientation),_print_level(printLevel)
  36   {
  37     _connectT=meshT.getConnectivityPtr();
  38     _connectS=meshS.getConnectivityPtr();
  39     _connIndexT=meshT.getConnectivityIndexPtr();
  40     _connIndexS=meshS.getConnectivityIndexPtr();
  41     _coordsT=meshT.getCoordinatesPtr();
  42     _coordsS=meshS.getCoordinatesPtr();
  43   }
  44
  45   template<class MyMeshType, class MyMatrix>
  46   PlanarIntersector<MyMeshType,MyMatrix>::~PlanarIntersector()
  47   {
  48   }
  49
  50   /*!
  51     \brief creates the bounding boxes for all the cells of mesh \a mesh
  52
  53     The method accepts mixed meshes (containing triangles and quadrangles).
  54     The vector returned is of dimension 6*nb_elems with bounding boxes stored as xmin1, xmax1, ymin1, ymax1, zmin1, zmax1, xmin2, xmax2, ymin2,...
  55     The returned pointer must be deleted by the calling code.
  56
  57     \param mesh structure pointing to the mesh
  58     \param bbox vector containing the bounding boxes
  59   */
  60   template<class MyMeshType, class MyMatrix>
  61   void PlanarIntersector<MyMeshType,MyMatrix>::createBoundingBoxes(const MyMeshType& mesh, std::vector<double>& bbox)
  62   {
  63     /* We build the segment tree for locating possible matching intersections*/
  64     long nbelems = mesh.getNumberOfElements();
  65     bbox.resize(2*SPACEDIM* nbelems);
  66     const double* coords = mesh.getCoordinatesPtr();
  67     const ConnType* conn = mesh.getConnectivityPtr();
  68     const ConnType* conn_index = mesh.getConnectivityIndexPtr();
  69     int ibox=0;
  70     for(long icell=0; icell<nbelems; icell++)
  71       {
  72         int nb_nodes_per_elem =conn_index[icell+1]-conn_index[icell];
  73         //initializing bounding box limits
  74         for(int idim=0; idim<SPACEDIM; idim++)
  75           {
  76             bbox[2*SPACEDIM*ibox+2*idim]   =  std::numeric_limits<double>::max();
  77             bbox[2*SPACEDIM*ibox+2*idim+1] = -std::numeric_limits<double>::max();
  78           }
  79         //updating the bounding box with each node of the element
  80         for (int j=0; j<nb_nodes_per_elem; j++)
  81           {
  82             const double* coord_node=coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(conn[OTT<ConnType,numPol>::conn2C(conn_index[icell]+j)]);
  83             for(int idim=0; idim<SPACEDIM; idim++)
  84               {
  85                 double x=*(coord_node+idim);
  86                 bbox[ibox*2*SPACEDIM + 2*idim]   = (bbox[ibox*2*SPACEDIM + 2*idim]  <x)?bbox[ibox*2*SPACEDIM + 2*idim  ]:x;
  87                 bbox[ibox*2*SPACEDIM + 2*idim+1] = (bbox[ibox*2*SPACEDIM + 2*idim+1]>x)?bbox[ibox*2*SPACEDIM + 2*idim+1]:x;
  88               }
  89           }
  90         ibox++;
  91       }
  92   }
  93
  94   /*!
  95     Computes the bouding box of a given element. iP in numPol mode.
  96   */
  97   template<class MyMeshType, class MyMatrix>
  98   void PlanarIntersector<MyMeshType,MyMatrix>::getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes)
  99   {
 100     const double* coords = mesh.getCoordinatesPtr();
 101     const ConnType* conn_index = mesh.getConnectivityIndexPtr();
 102     const ConnType* conn = mesh.getConnectivityPtr();
 103     //initializing bounding box limits
 104     for(int idim=0; idim<SPACEDIM; idim++)
 105       {
 106         bb[2*idim  ] =  std::numeric_limits<double>::max();
 107         bb[2*idim+1] = -std::numeric_limits<double>::max();
 108       }
 109
 110     for (ConnType i=0; i<nb_nodes; i++)
 111       {
 112         //MN: iP= cell index, not node index, use of connectivity array ?
 113         const double* coord_node=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(conn[OTT<ConnType,numPol>::conn2C(conn_index[OTT<ConnType,numPol>::ind2C(iP)]+i)]));
 114         for(int idim=0; idim<SPACEDIM; idim++)
 115           {
 116             double x = *(coord_node+idim);
 117             //double y = *(mesh.getCoordinates(MED_FULL_INTERLACE)+3*(iP+i)+1);
 118             bb[2*idim  ] = (x<bb[2*idim  ])?x:bb[2*idim  ];
 119             bb[2*idim+1] = (x>bb[2*idim+1])?x:bb[2*idim+1];
 120           }
 121       }
 122   }
 123
 124   /*! Readjusts a set of bounding boxes so that they are extended
 125     in all dimensions for avoiding missing interesting intersections
 126
 127     \param bbox vector containing the bounding boxes
 128   */
 129   template<class MyMeshType, class MyMatrix>
 130   void PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double surf3DAdjustmentEps, double surf3DAdjustmentEpsAbs)
 131   {
 132     /* We build the segment tree for locating possible matching intersections*/
 133
 134     long size = bbox.size()/(2*SPACEDIM);
 135     for (int i=0; i<size; i++)
 136       {
 137         double max=- std::numeric_limits<double>::max();
 138         for(int idim=0; idim<SPACEDIM; idim++)
 139           {
 140             double Dx=bbox[i*2*SPACEDIM+1+2*idim]-bbox[i*2*SPACEDIM+2*idim];
 141             max=(max<Dx)?Dx:max;
 142           }
 143         for(int idim=0; idim<SPACEDIM; idim++)
 144           {
 145             bbox[i*2*SPACEDIM+2*idim  ] -= surf3DAdjustmentEps*max+surf3DAdjustmentEpsAbs;
 146             bbox[i*2*SPACEDIM+2*idim+1] += surf3DAdjustmentEps*max+surf3DAdjustmentEpsAbs;
 147           }
 148       }
 149   }
 150
 151   /*!
 152    * @param icellT id in target mesh in format of MyMeshType.
 153    * @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
 154    */
 155   template<class MyMeshType, class MyMatrix>
 156   void PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT)
 157   {
 158     int nbNodesT=_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1]-_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
 159     coordsT.resize(SPACEDIM*nbNodesT);
 160     for (ConnType iT=0; iT<nbNodesT; iT++)
 161       for(int idim=0; idim<SPACEDIM; idim++)
 162         coordsT[SPACEDIM*iT+idim]=_coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
 163   }
 164
 165   /*!
 166    * @param icellS id in source mesh in format of MyMeshType.
 167    * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
 168    */
 169   template<class MyMeshType, class MyMatrix>
 170   void PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS)
 171   {
 172     int nbNodesS=_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1]-_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
 173     coordsS.resize(SPACEDIM*nbNodesS);
 174     for (ConnType iS=0; iS<nbNodesS; iS++)
 175       for(int idim=0; idim<SPACEDIM; idim++)
 176         coordsS[SPACEDIM*iS+idim]=_coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
 177   }
 178
 179   /*!
 180    * @param icellT id in target mesh in format of MyMeshType.
 181    * @param offset is a value in C format that indicates the number of circular permutation.
 182    * @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
 183    */
 184   template<class MyMeshType, class MyMatrix>
 185   void PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinatesPermute(ConnType icellT, int offset, std::vector<double>& coordsT)
 186   {
 187     int nbNodesT=_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1]-_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
 188     coordsT.resize(SPACEDIM*nbNodesT);
 189     for (ConnType iTTmp=0; iTTmp<nbNodesT; iTTmp++)
 190       {
 191         ConnType iT=(iTTmp+offset)%nbNodesT;
 192         for(int idim=0; idim<SPACEDIM; idim++)
 193           coordsT[SPACEDIM*iTTmp+idim]=_coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
 194       }
 195   }
 196
 197   /*!
 198    * @param icellS id in source mesh in format of MyMeshType.
 199    * @param offset is a value in C format that indicates the number of circular permutation.
 200    * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
 201    */
 202   template<class MyMeshType, class MyMatrix>
 203   void PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinatesPermute(ConnType icellS, int offset, std::vector<double>& coordsS)
 204   {
 205     int nbNodesS=_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1]-_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
 206     coordsS.resize(SPACEDIM*nbNodesS);
 207     for (ConnType iSTmp=0; iSTmp<nbNodesS; iSTmp++)
 208       {
 209         ConnType iS=(iSTmp+offset)%nbNodesS;
 210         for(int idim=0; idim<SPACEDIM; idim++)
 211           coordsS[SPACEDIM*iSTmp+idim]=_coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
 212       }
 213   }
 214
 215   /*!
 216    * @param icellT id in target mesh in format of MyMeshType.
 217    * @param icellS id in source mesh in format of MyMeshType.
 218    * @param nbNodesT nb of nodes of the target cell.
 219    * @param nbNodesS nb of nodes of the source cell.
 220    * @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
 221    * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
 222    * @param orientation is an output value too, only set if SPACEDIM==3.
 223    */
 224   template<class MyMeshType, class MyMatrix>
 225   void PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, std::vector<double>& coordsT, std::vector<double>& coordsS, int& orientation)
 226   {
 227     coordsT.resize(SPACEDIM*nbNodesT);
 228     coordsS.resize(SPACEDIM*nbNodesS);
 229     for (int idim=0; idim<SPACEDIM; idim++)
 230       {
 231         for (ConnType iT=0; iT<nbNodesT; iT++)
 232           coordsT[SPACEDIM*iT+idim] = _coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
 233         for (ConnType iS=0; iS<nbNodesS; iS++)
 234           coordsS[SPACEDIM*iS+idim] = _coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
 235       }
 236
 237     //project cells T and S on the median plane and rotate the median plane
 238     if(SPACEDIM==3)
 239       orientation = projectionThis(&coordsT[0], &coordsS[0], nbNodesT, nbNodesS);
 240
 241     //DEBUG PRINTS
 242     if(_print_level >= 3)
 243       {
 244         std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl;
 245         std::cout << std::endl << "icellT= " << icellT << ", nb nodes T= " <<  nbNodesT << std::endl;
 246         for(int iT =0; iT< nbNodesT; iT++)
 247           {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsT[SPACEDIM*iT+idim] << " "; std::cout << std::endl;}
 248         std::cout << std::endl << "icellS= " << icellS << ", nb nodes S= " <<  nbNodesS << std::endl;
 249         for(int iS =0; iS< nbNodesS; iS++)
 250           {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsS[SPACEDIM*iS+idim]<< " "; std::cout << std::endl; }
 251       }
 252   }
 253
 254   /*!
 255    * Filtering out zero surfaces and badly oriented surfaces
 256    * _orientation = -1,0,1,2
 257    * -1 : the intersection is taken into account if target and cells have different orientation
 258    * 0 : the intersection is always taken into account
 259    * 1 : the intersection is taken into account if target and cells have the same orientation
 260    * 2 : the absolute value of intersection is always taken into account
 261    */
 262   template<class MyMeshType, class MyMatrix>
 263   double PlanarIntersector<MyMeshType,MyMatrix>::getValueRegardingOption(double val) const
 264   {
 265     if(_orientation==0)
 266       return val;
 267     if(_orientation==2)
 268       return fabs(val);
 269     if (( val > 0.0 && _orientation==1) || ( val < 0.0 && _orientation==-1 ))
 270       return _orientation*val;
 271     return 0.;
 272   }
 273
 274   template<class MyMeshType, class MyMatrix>
 275   int PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB)
 276   {
 277     return projection(Coords_A,Coords_B,nb_NodesA,nb_NodesB,_dim_caracteristic*_precision,_max_distance_3Dsurf_intersect,_median_plane,_do_rotate);
 278   }
 279
 280   template<class MyMeshType, class MyMatrix>
 281   int PlanarIntersector<MyMeshType,MyMatrix>::projection(double *Coords_A, double *Coords_B,
 282                                                          int nb_NodesA, int nb_NodesB, double epsilon, double md3DSurf, double median_plane, bool do_rotate)
 283   {
 284     double normal_A[3]={0,0,0};
 285     double normal_B[3]={0,0,0};
 286     double linear_comb[3];
 287     double proj;
 288     bool same_orientation;
 289
 290     //Find the normal to cells A and B
 291     int i_A1=1;
 292     while(i_A1<nb_NodesA && distance2<SPACEDIM>(Coords_A,&Coords_A[SPACEDIM*i_A1])< epsilon) i_A1++;
 293     int i_A2=i_A1+1;
 294     crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
 295     double normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
 296     while(i_A2<nb_NodesA && normA < epsilon)
 297       {
 298         crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
 299         i_A2++;
 300         normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
 301
 302       }
 303     int i_B1=1;
 304     while(i_B1<nb_NodesB && distance2<SPACEDIM>(Coords_B,Coords_B+SPACEDIM*i_B1)< epsilon) i_B1++;
 305     int i_B2=i_B1+1;
 306     crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
 307     double normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
 308     while(i_B2<nb_NodesB && normB < epsilon)
 309       {
 310         crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
 311         i_B2++;
 312         normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
 313       }
 314
 315     //fabien option
 316     if(md3DSurf>0.)
 317       {
 318         double coords_GA[3];
 319         for (int i=0;i<3;i++)
 320           {
 321             coords_GA[i]=0.;
 322             for (int j=0;j<nb_NodesA;j++)
 323               coords_GA[i]+=Coords_A[3*j+i];
 324             coords_GA[i]/=nb_NodesA;
 325           }
 326         double G1[3],G2[3],G3[3];
 327       for (int i=0;i<3;i++)
 328         {
 329           G1[i]=Coords_B[i]-coords_GA[i];
 330           G2[i]=Coords_B[i+3]-coords_GA[i];
 331           G3[i]=Coords_B[i+6]-coords_GA[i];
 332         }
 333       double prodvect[3];
 334       prodvect[0]=G1[1]*G2[2]-G1[2]*G2[1];
 335       prodvect[1]=G1[2]*G2[0]-G1[0]*G2[2];
 336       prodvect[2]=G1[0]*G2[1]-G1[1]*G2[0];
 337       double prodscal=prodvect[0]*G3[0]+prodvect[1]*G3[1]+prodvect[2]*G3[2];
 338       if(fabs(prodscal)>md3DSurf)
 339         return 0;
 340       }
 341     if(i_A2<nb_NodesA && i_B2<nb_NodesB)
 342       {
 343         //Build the normal of the median plane
 344         same_orientation = dotprod<SPACEDIM>(normal_A,normal_B)>=0;
 345
 346         if(!same_orientation)
 347           for(int idim =0; idim< SPACEDIM; idim++) normal_A[idim] *=-1;
 348
 349         double normBB= sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
 350
 351         for(int idim =0; idim< SPACEDIM; idim++)
 352           linear_comb[idim] = median_plane*normal_A[idim]/normA + (1-median_plane)*normal_B[idim]/normBB;
 353         double norm= sqrt(dotprod<SPACEDIM>(linear_comb,linear_comb));
 354
 355         //Necessarily: norm>epsilon, no need to check
 356         for(int idim =0; idim< SPACEDIM; idim++) linear_comb[idim]/=norm;
 357
 358         //Project the nodes of A and B on the median plane
 359         for(int i_A=0; i_A<nb_NodesA; i_A++)
 360           {
 361             proj = dotprod<SPACEDIM>(&Coords_A[SPACEDIM*i_A],linear_comb);
 362             for(int idim =0; idim< SPACEDIM; idim++)
 363               Coords_A[SPACEDIM*i_A+idim] -=  proj*linear_comb[idim];
 364           }
 365         for(int i_B=0; i_B<nb_NodesB; i_B++)
 366           {
 367             proj = dotprod<SPACEDIM>(Coords_B+SPACEDIM*i_B,linear_comb);
 368             for(int idim =0; idim< SPACEDIM; idim++)
 369               Coords_B[SPACEDIM*i_B+idim] -=  proj*linear_comb[idim];
 370           }
 371
 372         //Buid the matrix sending  A into the Oxy plane and apply it to A and B
 373         if(do_rotate)
 374           {
 375             TranslationRotationMatrix rotation;
 376             //rotate3DTriangle(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2], rotation);
 377             rotate3DTriangle(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2, rotation);
 378             for (int i=0; i<nb_NodesA; i++)    rotation.transform_vector(Coords_A+SPACEDIM*i);
 379             for (int i=0; i<nb_NodesB; i++)    rotation.transform_vector(Coords_B+SPACEDIM*i);
 380           }
 381         if(same_orientation)
 382           return 1;
 383         else return -1;
 384       }
 385     else
 386       {
 387         std::cout << " Degenerated cell " << "epsilon = " << epsilon << std::endl;
 388         std::cout << " i_A1= " << i_A1 << " i_A2= " << i_A2 << std::endl;
 389         std::cout << " distance2<SPACEDIM>(Coords_A,&Coords_A[i_A1])= " <<  distance2<SPACEDIM>(Coords_A,&Coords_A[i_A1]) << std::endl;
 390         std::cout << "abs(normal_A) = " << fabs(normal_A[0]) << " ; " <<fabs( normal_A[1]) << " ; " << fabs(normal_A[2]) << std::endl;
 391         std::cout << " i_B1= " << i_B1 << " i_B2= " << i_B2 << std::endl;
 392         std::cout << " distance2<SPACEDIM>(&Coords_B[0],&Coords_B[i_B1])= " <<  distance2<SPACEDIM>(Coords_B,Coords_B+i_B1) << std::endl;
 393         std::cout << "normal_B = " << normal_B[0] << " ; " << normal_B[1] << " ; " << normal_B[2] << std::endl;
 394
 395         return 1;
 396       }
 397   }
 398
 399   template<class MyMeshType, class MyMatrix>
 400   void PlanarIntersector<MyMeshType,MyMatrix>::rotate3DTriangle(double* PP1, double*PP2, double*PP3,
 401                                                                 TranslationRotationMatrix& rotation_matrix)
 402   {
 403     //initializes
 404     rotation_matrix.translate(PP1);
 405
 406     double P2w[3];
 407     double P3w[3];
 408     P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2];
 409     P3w[0]=PP3[0]; P3w[1]=PP3[1];P3w[2]=PP3[2];
 410
 411     // translating to set P1 at the origin
 412     for (int i=0; i<3; i++)
 413       {
 414         P2w[i]-=PP1[i];
 415         P3w[i]-=PP1[i];
 416       }
 417
 418     // rotating to set P2 on the Oxy plane
 419     TranslationRotationMatrix A;
 420     A.rotate_x(P2w);
 421     A.rotate_vector(P3w);
 422     rotation_matrix.multiply(A);
 423
 424     //rotating to set P2 on the Ox axis
 425     TranslationRotationMatrix B;
 426     B.rotate_z(P2w);
 427     B.rotate_vector(P3w);
 428     rotation_matrix.multiply(B);
 429
 430     //rotating to set P3 on the Oxy plane
 431     TranslationRotationMatrix C;
 432     C.rotate_x(P3w);
 433     rotation_matrix.multiply(C);
 434   }
 435 }
 436
 437 #endif