src/INTERP_KERNEL/TriangulationIntersector.txx

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