src/INTERP_KERNEL/TriangulationIntersector.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 __TRIANGULATIONINTERSECTOR_TXX__
  20 #define __TRIANGULATIONINTERSECTOR_TXX__
  21
  22 #include "TriangulationIntersector.hxx"
  23 #include "PlanarIntersectorP0P0.txx"
  24 #include "PlanarIntersectorP0P1.txx"
  25 #include "PlanarIntersectorP1P0.txx"
  26 #include "PlanarIntersectorP1P0Bary.txx"
  27
  28 #include "InterpolationUtils.hxx"
  29 #include "PlanarIntersector.hxx"
  30
  31 #include <iostream>
  32
  33 #define TRI_INTERSECTOR TriangulationIntersector<MyMeshType,MyMatrix,InterpType>
  34 #define TRI_INTER_TEMPLATE template<class MyMeshType, class MyMatrix, \
  35                     template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
  36
  37 namespace INTERP_KERNEL
  38 {
  39   TRI_INTER_TEMPLATE
  40   TRI_INTERSECTOR::TriangulationIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
  41                                             double DimCaracteristic, double Precision, double md3DSurf,
  42                                             double MedianPlane, int orientation, int PrintLevel)
  43     :InterpType<MyMeshType,MyMatrix,TRI_INTERSECTOR >(meshT,meshS,DimCaracteristic, Precision, md3DSurf,
  44                                                       MedianPlane, true, orientation, PrintLevel)
  45   {
  46     if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 1)
  47       {
  48         std::cout << "  - intersection type = triangles " << std::endl;
  49         if(SPACEDIM==3) std::cout << "_do_rotate = true"<< std::endl;
  50       }
  51   }
  52
  53   TRI_INTER_TEMPLATE
  54   double TRI_INTERSECTOR::intersectGeometry(ConnType icellT,   ConnType icellS,
  55                                             ConnType nbNodesT, ConnType nbNodesS)
  56   {
  57     double result = 0.;
  58     int orientation = 1;
  59
  60     //Obtain the coordinates of T and S
  61     std::vector<double> CoordsT;
  62     std::vector<double> CoordsS;
  63     PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(icellT,icellS,nbNodesT,nbNodesS,CoordsT,CoordsS,orientation);
  64     //Compute the intersection area
  65     double area[SPACEDIM];
  66     for(ConnType iT = 1; iT<nbNodesT-1; iT++)
  67       {
  68         for(ConnType iS = 1; iS<nbNodesS-1; iS++)
  69           {
  70             std::vector<double> inter;
  71             INTERP_KERNEL::intersec_de_triangle(&CoordsT[0],&CoordsT[SPACEDIM*iT],&CoordsT[SPACEDIM*(iT+1)],
  72                                                 &CoordsS[0],&CoordsS[SPACEDIM*iS],&CoordsS[SPACEDIM*(iS+1)],
  73                                                 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
  74                                                 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
  75             ConnType nb_inter=((ConnType)inter.size())/2;
  76             if(nb_inter >3) inter=reconstruct_polygon(inter);
  77             for(ConnType i = 1; i<nb_inter-1; i++)
  78               {
  79                 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
  80                 result +=0.5*fabs(area[0]);
  81               }
  82             //DEBUG prints
  83             if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
  84               {
  85                 std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl;
  86                 for(ConnType i=0; i< nb_inter; i++)
  87                   {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; }
  88               }
  89           }
  90       }
  91
  92     //DEBUG PRINTS
  93     if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
  94       std::cout << std::endl <<"Intersection area = " << result << std::endl;
  95
  96     return orientation*result;
  97   }
  98
  99   TRI_INTER_TEMPLATE
 100   double TRI_INTERSECTOR::intersectGeometryWithQuadrangle(const double             * quadrangle,
 101                                                           const std::vector<double>& sourceCoords,
 102                                                           bool                       isSourceQuad)
 103   {
 104     double result = 0.;
 105     ConnType nbNodesS=sourceCoords.size()/SPACEDIM;
 106     //Compute the intersection area
 107     double area[SPACEDIM];
 108     for(ConnType iT = 1; iT<3; iT++)
 109       {
 110         for(ConnType iS = 1; iS<nbNodesS-1; iS++)
 111           {
 112             std::vector<double> inter;
 113             INTERP_KERNEL::intersec_de_triangle(quadrangle,&quadrangle[SPACEDIM*iT],&quadrangle[SPACEDIM*(iT+1)],
 114                                                 &sourceCoords[0],&sourceCoords[SPACEDIM*iS],&sourceCoords[SPACEDIM*(iS+1)],
 115                                                 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
 116                                                 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
 117             ConnType nb_inter=((ConnType)inter.size())/2;
 118             if(nb_inter >3) inter=reconstruct_polygon(inter);
 119             for(ConnType i = 1; i<nb_inter-1; i++)
 120               {
 121                 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
 122                 result +=0.5*fabs(area[0]);
 123               }
 124             //DEBUG prints
 125             if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
 126               {
 127                 std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl;
 128                 for(ConnType i=0; i< nb_inter; i++)
 129                   {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; }
 130               }
 131           }
 132       }
 133
 134     //DEBUG PRINTS
 135     if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
 136       std::cout << std::endl <<"Intersection area = " << result << std::endl;
 137
 138     return result;
 139   }
 140
 141   TRI_INTER_TEMPLATE
 142   double TRI_INTERSECTOR::intersectGeometryGeneral(const std::vector<double>& targetCoords,
 143                                                    const std::vector<double>& sourceCoords)
 144   {
 145     double result = 0.;
 146     ConnType nbNodesS=sourceCoords.size()/SPACEDIM;
 147     ConnType nbNodesT=targetCoords.size()/SPACEDIM;
 148     //Compute the intersection area
 149     double area[SPACEDIM];
 150     for(ConnType iT = 1; iT<nbNodesT-1; iT++)
 151       {
 152         for(ConnType iS = 1; iS<nbNodesS-1; iS++)
 153           {
 154             std::vector<double> inter;
 155             INTERP_KERNEL::intersec_de_triangle(&targetCoords[0],&targetCoords[SPACEDIM*iT],&targetCoords[SPACEDIM*(iT+1)],
 156                                                 &sourceCoords[0],&sourceCoords[SPACEDIM*iS],&sourceCoords[SPACEDIM*(iS+1)],
 157                                                 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
 158                                                 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
 159             ConnType nb_inter=((ConnType)inter.size())/2;
 160             if(nb_inter >3) inter=reconstruct_polygon(inter);
 161             for(ConnType i = 1; i<nb_inter-1; i++)
 162             {
 163               INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
 164               result +=0.5*fabs(area[0]);
 165             }
 166           }
 167       }
 168     return result;
 169   }
 170
 171   //================================================================================
 172   /*!
 173    * \brief Intersect a triangle and a polygon for P1P0 barycentric algorithm
 174    *  \param targetCell - list of coordinates of target polygon in full interlace
 175    *  \param targetCellQuadratic - specifies if target polygon is quadratic or not
 176    *  \param sourceTria - list of coordinates of source triangle
 177    *  \param res - coefficients a,b and c associated to nodes of sourceTria
 178    */
 179   //================================================================================
 180
 181   TRI_INTER_TEMPLATE
 182   double TRI_INTERSECTOR::intersectGeoBary(const std::vector<double>& targetCell,
 183                                            bool                       targetCellQuadratic,
 184                                            const double *             sourceTria,
 185                                            std::vector<double>&       res)
 186   {
 187     std::vector<const double* > sourceCell(3);
 188     sourceCell[0] = &sourceTria[0];
 189     sourceCell[1] = &sourceTria[SPACEDIM];
 190     sourceCell[2] = &sourceTria[SPACEDIM*2];
 191
 192     //Compute the intersection area
 193     double inter_area[SPACEDIM], total_area = 0.;
 194     double total_barycenter[SPACEDIM]={0.,0.};
 195
 196     const ConnType nbNodesT=targetCell.size()/SPACEDIM;
 197     for(ConnType iT = 1; iT<nbNodesT-1; iT++)
 198     {
 199       std::vector<double> inter;
 200       INTERP_KERNEL::intersec_de_triangle(&targetCell[0],&targetCell[SPACEDIM*iT],&targetCell[SPACEDIM*(iT+1)],
 201                                           sourceCell[0], sourceCell[1], sourceCell[2],
 202                                           inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
 203                                           PlanarIntersector<MyMeshType,MyMatrix>::_precision);
 204       ConnType nb_inter=((ConnType)inter.size())/2;
 205       if(nb_inter >3) inter=reconstruct_polygon(inter);
 206       for(ConnType i = 1; i<nb_inter-1; i++)
 207       {
 208         INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],inter_area);
 209         inter_area[0] = 0.5 * fabs( inter_area[0] );
 210         total_area += inter_area[0];
 211         std::vector<double> inter_bary=INTERP_KERNEL::bary_poly(inter);
 212         total_barycenter[0] += inter_area[0] * inter_bary[0];
 213         total_barycenter[1] += inter_area[0] * inter_bary[1];
 214       }
 215     }
 216     if ( total_area > std::numeric_limits<double>::min() )
 217     {
 218       total_barycenter[0] /= total_area;
 219       total_barycenter[1] /= total_area;
 220       res.resize(3);
 221       barycentric_coords( sourceCell, &total_barycenter[0], &res[0]);
 222       res[0] *= total_area;
 223       res[1] *= total_area;
 224       res[2] *= total_area;
 225     }
 226     else
 227     {
 228       total_area = 0;
 229     }
 230     return total_area;
 231   }
 232
 233 }
 234
 235 #endif