src/INTERP_KERNEL/PlanarIntersector.txx

   1 //  Copyright (C) 2007-2008  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 medianPlane, bool doRotate, int orientation, int printLevel):
  33     _meshT(meshT),_meshS(meshS),
  34     _dim_caracteristic(dimCaracteristic),_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)
 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;
 146             bbox[i*2*SPACEDIM+2*idim+1] += Surf3DAdjustmentEps*max;
 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 icellS id in source mesh in format of MyMeshType.
 182    * @param nbNodesT nb of nodes of the target cell.
 183    * @param nbNodesS nb of nodes of the source cell.
 184    * @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
 185    * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
 186    * @param orientation is an output value too, only set if SPACEDIM==3.
 187    */
 188   template<class MyMeshType, class MyMatrix>
 189   void PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, std::vector<double>& coordsT, std::vector<double>& coordsS, int& orientation)
 190   {
 191     /*double epsilon=_precision*_dim_caracteristic;
 192       coordsT.resize(SPACEDIM*nbNodesT);
 193       coordsS.resize(SPACEDIM*nbNodesS);
 194       int nb_dist_NodesT=nbNodesT;
 195       int nb_dist_NodesS=nbNodesS;
 196       int i_last = nbNodesT - 1;
 197       const double * Pi_last=_coordsT +_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+i_last)];
 198
 199       for (int iT=0; iT<nbNodesT; iT++)
 200       {
 201       const double * PiT = _coordsT + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)]);
 202       if(distance2<SPACEDIM>(Pi_last, PiT)>epsilon)
 203       {
 204       for (int idim=0; idim<SPACEDIM; idim++)
 205       coordsT[SPACEDIM*iT+idim]=PiT[idim];
 206       i_last=iT; Pi_last = PiT;
 207       }
 208       else
 209       nb_dist_NodesT--;
 210       }
 211       coordsT.resize(nb_dist_NodesT*SPACEDIM);
 212       i_last = nbNodesS - 1;
 213       Pi_last=_coordsS + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+i_last)]);
 214       for (int iS=0; iS<nbNodesS; iS++)
 215       {
 216       const double * PiS=_coordsS+SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)]);
 217       if(distance2<SPACEDIM>(Pi_last, PiS)>epsilon)
 218       {
 219       for (int idim=0; idim<SPACEDIM; idim++)
 220       coordsS[SPACEDIM*iS+idim]=PiS[idim];
 221       i_last=iS; Pi_last = PiS;
 222       }
 223       else
 224       nb_dist_NodesS--;
 225       }
 226       coordsS.resize(nb_dist_NodesS*SPACEDIM);
 227       //project cells T and S on the median plane
 228       // and rotate the median plane
 229       if(SPACEDIM==3)
 230       orientation = projectionThis(&coordsT[0], &coordsS[0], nb_dist_NodesT, nb_dist_NodesS);
 231
 232       //DEBUG PRINTS
 233       if(_print_level >= 3)
 234       {
 235       std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl;
 236       std::cout << std::endl << "icellT= " << icellT << ", nb nodes T= " <<  nbNodesT << std::endl;
 237       for(int iT =0; iT< nbNodesT; iT++)
 238       {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsT[SPACEDIM*iT+idim] << " "; std::cout << std::endl;}
 239       std::cout << std::endl << "icellS= " << icellS << ", nb nodes S= " <<  nbNodesS << std::endl;
 240       for(int iS =0; iS< nbNodesS; iS++)
 241       {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsS[SPACEDIM*iS+idim]<< " "; std::cout << std::endl; }
 242       }*/
 243     coordsT.resize(SPACEDIM*nbNodesT);
 244     coordsS.resize(SPACEDIM*nbNodesS);
 245     for (int idim=0; idim<SPACEDIM; idim++)
 246       {
 247         for (ConnType iT=0; iT<nbNodesT; iT++)
 248           coordsT[SPACEDIM*iT+idim] = _coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
 249         for (ConnType iS=0; iS<nbNodesS; iS++)
 250           coordsS[SPACEDIM*iS+idim] = _coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
 251       }
 252
 253     //project cells T and S on the median plane and rotate the median plane
 254     if(SPACEDIM==3)
 255       orientation = projectionThis(&coordsT[0], &coordsS[0], nbNodesT, nbNodesS);
 256
 257     //DEBUG PRINTS
 258     if(_print_level >= 3)
 259       {
 260         std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl;
 261         std::cout << std::endl << "icellT= " << icellT << ", nb nodes T= " <<  nbNodesT << std::endl;
 262         for(int iT =0; iT< nbNodesT; iT++)
 263           {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsT[SPACEDIM*iT+idim] << " "; std::cout << std::endl;}
 264         std::cout << std::endl << "icellS= " << icellS << ", nb nodes S= " <<  nbNodesS << std::endl;
 265         for(int iS =0; iS< nbNodesS; iS++)
 266           {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsS[SPACEDIM*iS+idim]<< " "; std::cout << std::endl; }
 267       }
 268   }
 269
 270   template<class MyMeshType, class MyMatrix>
 271   int PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB)
 272   {
 273     return projection(Coords_A,Coords_B,nb_NodesA,nb_NodesB,_dim_caracteristic*_precision,_median_plane,_do_rotate);
 274   }
 275
 276   template<class MyMeshType, class MyMatrix>
 277   int PlanarIntersector<MyMeshType,MyMatrix>::projection(double *Coords_A, double *Coords_B,
 278                                                          int nb_NodesA, int nb_NodesB, double epsilon, double median_plane, bool do_rotate)
 279   {
 280     double normal_A[3]={0,0,0};
 281     double normal_B[3]={0,0,0};
 282     double linear_comb[3];
 283     double proj;
 284     bool same_orientation;
 285
 286     //Find the normal to cells A and B
 287     int i_A1=1;
 288     while(i_A1<nb_NodesA && distance2<SPACEDIM>(Coords_A,&Coords_A[SPACEDIM*i_A1])< epsilon) i_A1++;
 289     int i_A2=i_A1+1;
 290     crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
 291     double normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
 292     while(i_A2<nb_NodesA && normA < epsilon)
 293       {
 294         crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
 295         i_A2++;
 296         normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
 297
 298       }
 299     int i_B1=1;
 300     while(i_B1<nb_NodesB && distance2<SPACEDIM>(Coords_B,Coords_B+SPACEDIM*i_B1)< epsilon) i_B1++;
 301     int i_B2=i_B1+1;
 302     crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
 303     double normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
 304     while(i_B2<nb_NodesB && normB < epsilon)
 305       {
 306         crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
 307         i_B2++;
 308         normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
 309       }
 310
 311     if(i_A2<nb_NodesA && i_B2<nb_NodesB)
 312       {
 313         //Build the normal of the median plane
 314         same_orientation = dotprod<SPACEDIM>(normal_A,normal_B)>=0;
 315
 316         if(!same_orientation)
 317           for(int idim =0; idim< SPACEDIM; idim++) normal_A[idim] *=-1;
 318
 319         double normB= sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
 320
 321         for(int idim =0; idim< SPACEDIM; idim++)
 322           linear_comb[idim] = median_plane*normal_A[idim]/normA + (1-median_plane)*normal_B[idim]/normB;
 323         double norm= sqrt(dotprod<SPACEDIM>(linear_comb,linear_comb));
 324
 325         //Necessarily: norm>epsilon, no need to check
 326         for(int idim =0; idim< SPACEDIM; idim++) linear_comb[idim]/=norm;
 327
 328         //Project the nodes of A and B on the median plane
 329         for(int i_A=0; i_A<nb_NodesA; i_A++)
 330           {
 331             proj = dotprod<SPACEDIM>(&Coords_A[SPACEDIM*i_A],linear_comb);
 332             for(int idim =0; idim< SPACEDIM; idim++)
 333               Coords_A[SPACEDIM*i_A+idim] -=  proj*linear_comb[idim];
 334           }
 335         for(int i_B=0; i_B<nb_NodesB; i_B++)
 336           {
 337             proj = dotprod<SPACEDIM>(Coords_B+SPACEDIM*i_B,linear_comb);
 338             for(int idim =0; idim< SPACEDIM; idim++)
 339               Coords_B[SPACEDIM*i_B+idim] -=  proj*linear_comb[idim];
 340           }
 341
 342         //Buid the matrix sending  A into the Oxy plane and apply it to A and B
 343         if(do_rotate)
 344           {
 345             TranslationRotationMatrix rotation;
 346             //rotate3DTriangle(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2], rotation);
 347             rotate3DTriangle(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2, rotation);
 348             for (int i=0; i<nb_NodesA; i++)    rotation.transform_vector(Coords_A+SPACEDIM*i);
 349             for (int i=0; i<nb_NodesB; i++)    rotation.transform_vector(Coords_B+SPACEDIM*i);
 350           }
 351         if(same_orientation)
 352           return 1;
 353         else return -1;
 354       }
 355     else
 356       {
 357         std::cout << " Maille dégénérée " << "epsilon = " << epsilon << std::endl;
 358         std::cout << " i_A1= " << i_A1 << " i_A2= " << i_A2 << std::endl;
 359         std::cout << " distance2<SPACEDIM>(Coords_A,&Coords_A[i_A1])= " <<  distance2<SPACEDIM>(Coords_A,&Coords_A[i_A1]) << std::endl;
 360         std::cout << "abs(normal_A) = " << fabs(normal_A[0]) << " ; " <<fabs( normal_A[1]) << " ; " << fabs(normal_A[2]) << std::endl;
 361         std::cout << " i_B1= " << i_B1 << " i_B2= " << i_B2 << std::endl;
 362         std::cout << " distance2<SPACEDIM>(&Coords_B[0],&Coords_B[i_B1])= " <<  distance2<SPACEDIM>(Coords_B,Coords_B+i_B1) << std::endl;
 363         std::cout << "normal_B = " << normal_B[0] << " ; " << normal_B[1] << " ; " << normal_B[2] << std::endl;
 364
 365         return 1;
 366       }
 367   }
 368
 369   template<class MyMeshType, class MyMatrix>
 370   void PlanarIntersector<MyMeshType,MyMatrix>::rotate3DTriangle(double* PP1, double*PP2, double*PP3,
 371                                                                 TranslationRotationMatrix& rotation_matrix)
 372   {
 373     //initializes
 374     rotation_matrix.translate(PP1);
 375
 376     double P2w[3];
 377     double P3w[3];
 378     P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2];
 379     P3w[0]=PP3[0]; P3w[1]=PP3[1];P3w[2]=PP3[2];
 380
 381     // translating to set P1 at the origin
 382     for (int i=0; i<3; i++)
 383       {
 384         P2w[i]-=PP1[i];
 385         P3w[i]-=PP1[i];
 386       }
 387
 388     // rotating to set P2 on the Oxy plane
 389     TranslationRotationMatrix A;
 390     A.rotate_x(P2w);
 391     A.rotate_vector(P3w);
 392     rotation_matrix.multiply(A);
 393
 394     //rotating to set P2 on the Ox axis
 395     TranslationRotationMatrix B;
 396     B.rotate_z(P2w);
 397     B.rotate_vector(P3w);
 398     rotation_matrix.multiply(B);
 399
 400     //rotating to set P3 on the Oxy plane
 401     TranslationRotationMatrix C;
 402     C.rotate_x(P3w);
 403     rotation_matrix.multiply(C);
 404   }
 405 }
 406
 407 #endif