1 // Copyright (C) 2007-2014 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, or (at your option) any later version.
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 // Author : Anthony Geay (CEA/DEN)
20 #ifndef __CONVEXINTERSECTOR_TXX__
21 #define __CONVEXINTERSECTOR_TXX__
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"
31 #include "PolygonAlgorithms.txx"
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>
38 namespace INTERP_KERNEL
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)
47 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 1)
49 std::cout << " - intersection type = convex " << std::endl;
51 if(PlanarIntersector<MyMeshType,MyMatrix>::_do_rotate) std::cout << " _do_rotate = true" << std::endl;
52 else std::cout << " _do_rotate = false" << std::endl;
57 CONVINTERSECTOR_TEMPLATE
58 double CONVEX_INTERSECTOR_::intersectGeometry(ConnType icellT, ConnType icellS,
59 ConnType nbNodesT, ConnType nbNodesS)
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 CoordsT.size()/SPACEDIM, 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++)
76 INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
77 result +=0.5*norm<SPACEDIM>(area);
81 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
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;
89 return orientation*result;
92 CONVINTERSECTOR_TEMPLATE
93 double CONVEX_INTERSECTOR_::intersectGeometryWithQuadrangle(const double * quadrangle,
94 const std::vector<double>& sourceCoords,
98 int nbOfNodesS=sourceCoords.size()/SPACEDIM;
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],
104 double area[SPACEDIM];
105 int nb_inter =((int)inter.size())/SPACEDIM;
106 for(int i = 1; i<nb_inter-1; i++)
108 INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
109 result +=0.5*norm<SPACEDIM>(area);
113 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
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;
124 CONVINTERSECTOR_TEMPLATE
125 double CONVEX_INTERSECTOR_::intersectGeometryGeneral(const std::vector<double>& targetCoords,
126 const std::vector<double>& sourceCoords)
129 int nbOfNodesS=sourceCoords.size()/SPACEDIM;
130 int nbOfNodesT=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++)
139 INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
140 result +=0.5*norm<SPACEDIM>(area);
145 //================================================================================
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
153 //================================================================================
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)
162 double barycenter[SPACEDIM] = {0., 0.};
163 int nbOfNodesT=targetCell.size()/SPACEDIM;
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++)
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];
177 if ( area > std::numeric_limits<double>::min() )
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;
182 barycentric_coords<2>( sourceTria, &barycenter[0], &res[0]);