src/INTERP_KERNEL/ConvexIntersector.txx

   1 // Copyright (C) 2007-2021  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, or (at your option) any later version.
   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 __CONVEXINTERSECTOR_TXX__
  21 #define __CONVEXINTERSECTOR_TXX__
  22
  23 #include "ConvexIntersector.hxx"
  24 #include "PlanarIntersectorP0P0.txx"
  25 #include "PlanarIntersectorP0P1.txx"
  26 #include "PlanarIntersectorP1P0.txx"
  27 #include "PlanarIntersectorP1P1.txx"
  28 #include "PlanarIntersectorP1P0Bary.txx"
  29 #include "PlanarIntersectorP0P1Bary.txx"
  30
  31 #include "PolygonAlgorithms.txx"
  32
  33 #include <iostream>
  34
  35 #define CONVINTERSECTOR_TEMPLATE template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
  36 #define CONVEX_INTERSECTOR_ ConvexIntersector<MyMeshType,MyMatrix,InterpType>
  37
  38 namespace INTERP_KERNEL
  39 {
  40   CONVINTERSECTOR_TEMPLATE
  41   CONVEX_INTERSECTOR_::ConvexIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
  42                                          double dimCaracteristic, double precision, double md3DSurf, double minDot3DSurf,
  43                                          double medianPlane, bool doRotate , int oriantation, int printLevel)
  44     :InterpType<MyMeshType,MyMatrix,CONVEX_INTERSECTOR_ >(meshT,meshS,dimCaracteristic, precision, md3DSurf, minDot3DSurf, medianPlane, doRotate, oriantation, printLevel),
  45      _epsilon(precision*dimCaracteristic)
  46   {
  47     if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 1)
  48       {
  49         std::cout << " - intersection type = convex " << std::endl;
  50         if(SPACEDIM==3){
  51           if(PlanarIntersector<MyMeshType,MyMatrix>::_do_rotate) std::cout << "  _do_rotate = true" << std::endl;
  52           else std::cout << "  _do_rotate = false" << std::endl;
  53         }
  54       }
  55   }
  56
  57   CONVINTERSECTOR_TEMPLATE
  58   double CONVEX_INTERSECTOR_::intersectGeometry(ConnType icellT,   ConnType icellS,
  59                                                 ConnType nbNodesT, ConnType nbNodesS)
  60   {
  61     double result = 0;
  62     int orientation = 1;
  63
  64     /*** Obtain the coordinates of T and S ***/
  65     std::vector<double> CoordsT;
  66     std::vector<double> CoordsS;
  67     PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(icellT,icellS,nbNodesT,nbNodesS,CoordsT,CoordsS,orientation);
  68     /*** Compute the intersection area ***/
  69     INTERP_KERNEL::PolygonAlgorithms<SPACEDIM> P(_epsilon, PlanarIntersector<MyMeshType,MyMatrix>::_precision);
  70     std::deque<double> inter =  P.intersectConvexPolygons(&CoordsT[0], &CoordsS[0],
  71                                                           int(CoordsT.size())/SPACEDIM, int(CoordsS.size())/SPACEDIM);
  72     double area[SPACEDIM];
  73     int nb_inter =((int)inter.size())/SPACEDIM;
  74     for(int i = 1; i<nb_inter-1; i++)
  75       {
  76         INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
  77         result +=0.5*norm<SPACEDIM>(area);
  78       }
  79
  80     //DEBUG prints
  81     if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
  82       {
  83         std::cout << std::endl << "Number of nodes of the intersection = "<<  nb_inter << std::endl;
  84         for(int i=0; i<  nb_inter; i++)
  85           {for (int idim=0; idim<SPACEDIM; idim++) std::cout << inter[SPACEDIM*i+idim]<< " "; std::cout << std::endl;}
  86         std::cout << std::endl <<"Intersection area = " << result << std::endl;
  87       }
  88
  89     return orientation*result;
  90   }
  91
  92   CONVINTERSECTOR_TEMPLATE
  93   double CONVEX_INTERSECTOR_::intersectGeometryWithQuadrangle(const double             * quadrangle,
  94                                                               const std::vector<double>& sourceCoords,
  95                                                               bool                       isSourceQuad)
  96   {
  97     double result = 0;
  98     int nbOfNodesS=int(sourceCoords.size())/SPACEDIM;
  99
 100     /*** Compute the intersection area ***/
 101     INTERP_KERNEL::PolygonAlgorithms<SPACEDIM> P(_epsilon, PlanarIntersector<MyMeshType,MyMatrix>::_precision);
 102     std::deque<double> inter =  P.intersectConvexPolygons(quadrangle, &sourceCoords[0],
 103                                                           4, nbOfNodesS);
 104     double area[SPACEDIM];
 105     int nb_inter =((int)inter.size())/SPACEDIM;
 106     for(int i = 1; i<nb_inter-1; i++)
 107       {
 108         INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
 109         result +=0.5*norm<SPACEDIM>(area);
 110       }
 111
 112     //DEBUG prints
 113     if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
 114       {
 115         std::cout << std::endl << "Number of nodes of the intersection = "<<  nb_inter << std::endl;
 116         for(int i=0; i<  nb_inter; i++)
 117           {for (int idim=0; idim<SPACEDIM; idim++) std::cout << inter[SPACEDIM*i+idim]<< " "; std::cout << std::endl;}
 118         std::cout << std::endl <<"Intersection area = " << result << std::endl;
 119       }
 120
 121     return result;
 122   }
 123
 124   CONVINTERSECTOR_TEMPLATE
 125   double CONVEX_INTERSECTOR_::intersectGeometryGeneral(const std::vector<double>& targetCoords,
 126                                                        const std::vector<double>& sourceCoords)
 127   {
 128     double result = 0;
 129     int nbOfNodesS=int(sourceCoords.size())/SPACEDIM;
 130     int nbOfNodesT=int(targetCoords.size())/SPACEDIM;
 131     /*** Compute the intersection area ***/
 132     INTERP_KERNEL::PolygonAlgorithms<SPACEDIM> P(_epsilon, PlanarIntersector<MyMeshType,MyMatrix>::_precision);
 133     std::deque<double> inter =  P.intersectConvexPolygons(&targetCoords[0], &sourceCoords[0],
 134                                                           nbOfNodesT, nbOfNodesS);
 135     double area[SPACEDIM];
 136     int nb_inter =((int)inter.size())/SPACEDIM;
 137     for(int i = 1; i<nb_inter-1; i++)
 138       {
 139         INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
 140         result +=0.5*norm<SPACEDIM>(area);
 141       }
 142     return result;
 143   }
 144
 145   //================================================================================
 146   /*!
 147    * \brief Intersect a triangle and a polygon for P1P0 barycentric algorithm
 148    *  \param targetCell - list of coordinates of target polygon in full interlace
 149    *  \param targetCellQuadratic - specifies if target polygon is quadratic or not
 150    *  \param sourceTria - list of coordinates of source triangle
 151    *  \param res - coefficients a,b and c associated to nodes of sourceTria
 152    */
 153   //================================================================================
 154
 155   CONVINTERSECTOR_TEMPLATE
 156   double CONVEX_INTERSECTOR_::intersectGeoBary(const std::vector<double>& targetCell,
 157                                                bool                       targetCellQuadratic,
 158                                                const double *             sourceTria,
 159                                                std::vector<double>&       res)
 160   {
 161     double area = 0;
 162     double barycenter[SPACEDIM] = {0., 0.};
 163     int nbOfNodesT=targetCell.size()/SPACEDIM;
 164
 165     /*** Compute the intersection area ***/
 166     INTERP_KERNEL::PolygonAlgorithms<SPACEDIM> P(_epsilon, PlanarIntersector<MyMeshType,MyMatrix>::_precision);
 167     std::deque<double> inter =  P.intersectConvexPolygons(sourceTria, &targetCell[0], 3, nbOfNodesT);
 168     double cross[SPACEDIM];
 169     int nb_inter =((int)inter.size())/SPACEDIM;
 170     for(int i = 1; i<nb_inter-1; i++)
 171     {
 172       INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],cross);
 173       area += 0.5*norm<SPACEDIM>(cross);
 174       barycenter[0] += inter[SPACEDIM*i];
 175       barycenter[1] += inter[SPACEDIM*i+1];
 176     }
 177     if ( area > std::numeric_limits<double>::min() )
 178     {
 179       barycenter[0] = ( barycenter[0] + inter[0] + inter[SPACEDIM*(nb_inter-1)]  ) / nb_inter;
 180       barycenter[1] = ( barycenter[1] + inter[1] + inter[SPACEDIM*(nb_inter-1)+1]) / nb_inter;
 181       res.resize(3);
 182       barycentric_coords<2>( sourceTria, &barycenter[0], &res[0]);
 183       res[0] *= area;
 184       res[1] *= area;
 185       res[2] *= area;
 186     }
 187     else
 188     {
 189       area = 0;
 190     }
 191     return area;
 192   }
 193 }
 194
 195 #endif