1 // Copyright (C) 2007-2008 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 #ifndef __CONVEXINTERSECTOR_TXX__
20 #define __CONVEXINTERSECTOR_TXX__
22 #include "ConvexIntersector.hxx"
23 #include "PlanarIntersectorP0P0.txx"
24 #include "PlanarIntersectorP0P1.txx"
25 #include "PlanarIntersectorP1P0.txx"
27 #include "PolygonAlgorithms.txx"
31 namespace INTERP_KERNEL
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)
40 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 1)
42 std::cout << " - intersection type = convex " << std::endl;
44 if(PlanarIntersector<MyMeshType,MyMatrix>::_do_rotate) std::cout << " _do_rotate = true" << std::endl;
45 else std::cout << " _do_rotate = false" << std::endl;
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)
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++)
68 INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
69 result +=0.5*norm<SPACEDIM>(area);
73 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
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;
81 return orientation*result;
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)
88 int nbOfNodesS=sourceCoords.size()/SPACEDIM;
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],
94 double area[SPACEDIM];
95 int nb_inter =((int)inter.size())/SPACEDIM;
96 for(int i = 1; i<nb_inter-1; i++)
98 INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
99 result +=0.5*norm<SPACEDIM>(area);
103 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
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;