src/INTERP_KERNEL/TriangulationIntersector.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 __TRIANGULATIONINTERSECTOR_TXX__
  20 #define __TRIANGULATIONINTERSECTOR_TXX__
  21
  22 #include "TriangulationIntersector.hxx"
  23 #include "PlanarIntersectorP0P0.txx"
  24 #include "PlanarIntersectorP0P1.txx"
  25 #include "PlanarIntersectorP1P0.txx"
  26
  27 #include "InterpolationUtils.hxx"
  28 #include "PlanarIntersector.hxx"
  29
  30 #include <iostream>
  31
  32 namespace INTERP_KERNEL
  33 {
  34   template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
  35   TriangulationIntersector<MyMeshType,MyMatrix,InterpType>::TriangulationIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
  36                                                                               double DimCaracteristic, double Precision,
  37                                                                               double MedianPlane, int orientation, int PrintLevel)
  38     :InterpType<MyMeshType,MyMatrix,TriangulationIntersector<MyMeshType,MyMatrix,InterpType> >(meshT,meshS,DimCaracteristic, Precision, MedianPlane, true, orientation, PrintLevel)
  39   {
  40     if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 1)
  41       {
  42         std::cout << "  - intersection type = triangles " << std::endl;
  43         if(SPACEDIM==3) std::cout << "_do_rotate = true"<< std::endl;
  44       }
  45   }
  46
  47   template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
  48   double TriangulationIntersector<MyMeshType,MyMatrix,InterpType>::intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS)
  49   {
  50     double result = 0.;
  51     int orientation = 1;
  52
  53     //Obtain the coordinates of T and S
  54     std::vector<double> CoordsT;
  55     std::vector<double> CoordsS;
  56     PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(icellT,icellS,nbNodesT,nbNodesS,CoordsT,CoordsS,orientation);
  57     //Compute the intersection area
  58     double area[SPACEDIM];
  59     for(ConnType iT = 1; iT<nbNodesT-1; iT++)
  60       {
  61         for(ConnType iS = 1; iS<nbNodesS-1; iS++)
  62           {
  63             std::vector<double> inter;
  64             INTERP_KERNEL::intersec_de_triangle(&CoordsT[0],&CoordsT[SPACEDIM*iT],&CoordsT[SPACEDIM*(iT+1)],
  65                                                 &CoordsS[0],&CoordsS[SPACEDIM*iS],&CoordsS[SPACEDIM*(iS+1)],
  66                                                 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
  67                                                 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
  68             ConnType nb_inter=((ConnType)inter.size())/2;
  69             if(nb_inter >3) inter=reconstruct_polygon(inter);
  70             for(ConnType i = 1; i<nb_inter-1; i++)
  71               {
  72                 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
  73                 result +=0.5*fabs(area[0]);
  74               }
  75             //DEBUG prints
  76             if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
  77               {
  78                 std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl;
  79                 for(ConnType i=0; i< nb_inter; i++)
  80                   {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; }
  81               }
  82           }
  83       }
  84
  85     //DEBUG PRINTS
  86     if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
  87       std::cout << std::endl <<"Intersection area = " << result << std::endl;
  88
  89     return orientation*result;
  90   }
  91
  92   template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
  93   double TriangulationIntersector<MyMeshType,MyMatrix,InterpType>::intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector<double>& sourceCoords, bool isSourceQuad)
  94   {
  95     double result = 0.;
  96     ConnType nbNodesS=sourceCoords.size()/SPACEDIM;
  97     //Compute the intersection area
  98     double area[SPACEDIM];
  99     for(ConnType iT = 1; iT<3; iT++)
 100       {
 101         for(ConnType iS = 1; iS<nbNodesS-1; iS++)
 102           {
 103             std::vector<double> inter;
 104             INTERP_KERNEL::intersec_de_triangle(quadrangle,&quadrangle[SPACEDIM*iT],&quadrangle[SPACEDIM*(iT+1)],
 105                                                 &sourceCoords[0],&sourceCoords[SPACEDIM*iS],&sourceCoords[SPACEDIM*(iS+1)],
 106                                                 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
 107                                                 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
 108             ConnType nb_inter=((ConnType)inter.size())/2;
 109             if(nb_inter >3) inter=reconstruct_polygon(inter);
 110             for(ConnType i = 1; i<nb_inter-1; i++)
 111               {
 112                 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
 113                 result +=0.5*fabs(area[0]);
 114               }
 115             //DEBUG prints
 116             if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
 117               {
 118                 std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl;
 119                 for(ConnType i=0; i< nb_inter; i++)
 120                   {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; }
 121               }
 122           }
 123       }
 124
 125     //DEBUG PRINTS
 126     if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
 127       std::cout << std::endl <<"Intersection area = " << result << std::endl;
 128
 129     return result;
 130   }
 131 }
 132
 133 #endif