src/INTERP_KERNEL/ConvexIntersector.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 __CONVEXINTERSECTOR_TXX__
  20 #define __CONVEXINTERSECTOR_TXX__
  21
  22 #include "ConvexIntersector.hxx"
  23 #include "PlanarIntersectorP0P0.txx"
  24 #include "PlanarIntersectorP0P1.txx"
  25 #include "PlanarIntersectorP1P0.txx"
  26
  27 #include "PolygonAlgorithms.txx"
  28
  29 #include <iostream>
  30
  31 namespace INTERP_KERNEL
  32 {
  33   template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
  34   ConvexIntersector<MyMeshType,MyMatrix,InterpType>::ConvexIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
  35                                                                 double dimCaracteristic, double precision,
  36                                                                 double medianPlane, bool doRotate , int oriantation, int printLevel)
  37     :InterpType<MyMeshType,MyMatrix,ConvexIntersector<MyMeshType,MyMatrix,InterpType> >(meshT,meshS,dimCaracteristic, precision, medianPlane, doRotate, oriantation, printLevel),
  38      _epsilon(precision*dimCaracteristic)
  39   {
  40     if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 1)
  41       {
  42         std::cout << " - intersection type = convex " << std::endl;
  43         if(SPACEDIM==3){
  44           if(PlanarIntersector<MyMeshType,MyMatrix>::_do_rotate) std::cout << "  _do_rotate = true" << std::endl;
  45           else std::cout << "  _do_rotate = false" << std::endl;
  46         }
  47       }
  48   }
  49
  50   template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
  51   double ConvexIntersector<MyMeshType,MyMatrix,InterpType>::intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS)
  52   {
  53     double result = 0;
  54     int orientation = 1;
  55
  56     /*** Obtain the coordinates of T and S ***/
  57     std::vector<double> CoordsT;
  58     std::vector<double> CoordsS;
  59     PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(icellT,icellS,nbNodesT,nbNodesS,CoordsT,CoordsS,orientation);
  60     /*** Compute the intersection area ***/
  61     INTERP_KERNEL::PolygonAlgorithms<SPACEDIM> P(_epsilon, PlanarIntersector<MyMeshType,MyMatrix>::_precision);
  62     std::deque<double> inter =  P.intersectConvexPolygons(&CoordsT[0], &CoordsS[0],
  63                                                             CoordsT.size()/SPACEDIM, CoordsS.size()/SPACEDIM);
  64     double area[SPACEDIM];
  65     int nb_inter =((int)inter.size())/SPACEDIM;
  66     for(int i = 1; i<nb_inter-1; i++)
  67       {
  68         INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
  69         result +=0.5*norm<SPACEDIM>(area);
  70       }
  71
  72     //DEBUG prints
  73     if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
  74       {
  75         std::cout << std::endl << "Number of nodes of the intersection = "<<  nb_inter << std::endl;
  76         for(int i=0; i<  nb_inter; i++)
  77           {for (int idim=0; idim<SPACEDIM; idim++) std::cout << inter[SPACEDIM*i+idim]<< " "; std::cout << std::endl;}
  78         std::cout << std::endl <<"Intersection area = " << result << std::endl;
  79       }
  80
  81     return orientation*result;
  82   }
  83
  84   template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
  85   double ConvexIntersector<MyMeshType,MyMatrix,InterpType>::intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector<double>& sourceCoords, bool isSourceQuad)
  86   {
  87     double result = 0;
  88     int nbOfNodesS=sourceCoords.size()/SPACEDIM;
  89
  90     /*** Compute the intersection area ***/
  91     INTERP_KERNEL::PolygonAlgorithms<SPACEDIM> P(_epsilon, PlanarIntersector<MyMeshType,MyMatrix>::_precision);
  92     std::deque<double> inter =  P.intersectConvexPolygons(quadrangle, &sourceCoords[0],
  93                                                             4, nbOfNodesS);
  94     double area[SPACEDIM];
  95     int nb_inter =((int)inter.size())/SPACEDIM;
  96     for(int i = 1; i<nb_inter-1; i++)
  97       {
  98         INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
  99         result +=0.5*norm<SPACEDIM>(area);
 100       }
 101
 102     //DEBUG prints
 103     if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
 104       {
 105         std::cout << std::endl << "Number of nodes of the intersection = "<<  nb_inter << std::endl;
 106         for(int i=0; i<  nb_inter; i++)
 107           {for (int idim=0; idim<SPACEDIM; idim++) std::cout << inter[SPACEDIM*i+idim]<< " "; std::cout << std::endl;}
 108         std::cout << std::endl <<"Intersection area = " << result << std::endl;
 109       }
 110
 111     return result;
 112   }
 113 }
 114
 115 #endif