1 // Copyright (C) 2007-2013 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 // Author : Anthony Geay (CEA/DEN)
20 #ifndef __TRIANGULATIONINTERSECTOR_TXX__
21 #define __TRIANGULATIONINTERSECTOR_TXX__
23 #include "TriangulationIntersector.hxx"
24 #include "PlanarIntersectorP0P0.txx"
25 #include "PlanarIntersectorP0P1.txx"
26 #include "PlanarIntersectorP1P0.txx"
27 #include "PlanarIntersectorP1P0Bary.txx"
29 #include "InterpolationUtils.hxx"
30 #include "PlanarIntersector.hxx"
34 #define TRI_INTERSECTOR TriangulationIntersector<MyMeshType,MyMatrix,InterpType>
35 #define TRI_INTER_TEMPLATE template<class MyMeshType, class MyMatrix, \
36 template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
38 namespace INTERP_KERNEL
41 TRI_INTERSECTOR::TriangulationIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
42 double DimCaracteristic, double Precision, double md3DSurf, double minDot3DSurf,
43 double MedianPlane, int orientation, int PrintLevel)
44 :InterpType<MyMeshType,MyMatrix,TRI_INTERSECTOR >(meshT,meshS,DimCaracteristic, Precision, md3DSurf, minDot3DSurf,
45 MedianPlane, true, orientation, PrintLevel)
47 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 1)
49 std::cout << " - intersection type = triangles " << std::endl;
50 if(SPACEDIM==3) std::cout << "_do_rotate = true"<< std::endl;
55 double TRI_INTERSECTOR::intersectGeometry(ConnType icellT, ConnType icellS,
56 ConnType nbNodesT, ConnType nbNodesS)
61 //Obtain the coordinates of T and S
62 std::vector<double> CoordsT;
63 std::vector<double> CoordsS;
64 PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(icellT,icellS,nbNodesT,nbNodesS,CoordsT,CoordsS,orientation);
65 //Compute the intersection area
66 double area[SPACEDIM];
67 for(ConnType iT = 1; iT<nbNodesT-1; iT++)
69 for(ConnType iS = 1; iS<nbNodesS-1; iS++)
71 std::vector<double> inter;
72 INTERP_KERNEL::intersec_de_triangle(&CoordsT[0],&CoordsT[SPACEDIM*iT],&CoordsT[SPACEDIM*(iT+1)],
73 &CoordsS[0],&CoordsS[SPACEDIM*iS],&CoordsS[SPACEDIM*(iS+1)],
74 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
75 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
76 ConnType nb_inter=((ConnType)inter.size())/2;
77 if(nb_inter >3) inter=reconstruct_polygon(inter);
78 for(ConnType i = 1; i<nb_inter-1; i++)
80 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
81 result +=0.5*fabs(area[0]);
84 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
86 std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl;
87 for(ConnType i=0; i< nb_inter; i++)
88 {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; }
94 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
95 std::cout << std::endl <<"Intersection area = " << result << std::endl;
97 return orientation*result;
101 double TRI_INTERSECTOR::intersectGeometryWithQuadrangle(const double * quadrangle,
102 const std::vector<double>& sourceCoords,
106 ConnType nbNodesS=sourceCoords.size()/SPACEDIM;
107 //Compute the intersection area
108 double area[SPACEDIM];
109 for(ConnType iT = 1; iT<3; iT++)
111 for(ConnType iS = 1; iS<nbNodesS-1; iS++)
113 std::vector<double> inter;
114 INTERP_KERNEL::intersec_de_triangle(quadrangle,&quadrangle[SPACEDIM*iT],&quadrangle[SPACEDIM*(iT+1)],
115 &sourceCoords[0],&sourceCoords[SPACEDIM*iS],&sourceCoords[SPACEDIM*(iS+1)],
116 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
117 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
118 ConnType nb_inter=((ConnType)inter.size())/2;
119 if(nb_inter >3) inter=reconstruct_polygon(inter);
120 for(ConnType i = 1; i<nb_inter-1; i++)
122 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
123 result +=0.5*fabs(area[0]);
126 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
128 std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl;
129 for(ConnType i=0; i< nb_inter; i++)
130 {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; }
136 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
137 std::cout << std::endl <<"Intersection area = " << result << std::endl;
143 double TRI_INTERSECTOR::intersectGeometryGeneral(const std::vector<double>& targetCoords,
144 const std::vector<double>& sourceCoords)
147 ConnType nbNodesS=sourceCoords.size()/SPACEDIM;
148 ConnType nbNodesT=targetCoords.size()/SPACEDIM;
149 //Compute the intersection area
150 double area[SPACEDIM];
151 for(ConnType iT = 1; iT<nbNodesT-1; iT++)
153 for(ConnType iS = 1; iS<nbNodesS-1; iS++)
155 std::vector<double> inter;
156 INTERP_KERNEL::intersec_de_triangle(&targetCoords[0],&targetCoords[SPACEDIM*iT],&targetCoords[SPACEDIM*(iT+1)],
157 &sourceCoords[0],&sourceCoords[SPACEDIM*iS],&sourceCoords[SPACEDIM*(iS+1)],
158 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
159 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
160 ConnType nb_inter=((ConnType)inter.size())/2;
161 if(nb_inter >3) inter=reconstruct_polygon(inter);
162 for(ConnType i = 1; i<nb_inter-1; i++)
164 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
165 result +=0.5*fabs(area[0]);
172 //================================================================================
174 * \brief Intersect a triangle and a polygon for P1P0 barycentric algorithm
175 * \param targetCell - list of coordinates of target polygon in full interlace
176 * \param targetCellQuadratic - specifies if target polygon is quadratic or not
177 * \param sourceTria - list of coordinates of source triangle
178 * \param res - coefficients a,b and c associated to nodes of sourceTria
180 //================================================================================
183 double TRI_INTERSECTOR::intersectGeoBary(const std::vector<double>& targetCell,
184 bool targetCellQuadratic,
185 const double * sourceTria,
186 std::vector<double>& res)
188 std::vector<const double* > sourceCell(3);
189 sourceCell[0] = &sourceTria[0];
190 sourceCell[1] = &sourceTria[SPACEDIM];
191 sourceCell[2] = &sourceTria[SPACEDIM*2];
193 //Compute the intersection area
194 double inter_area[SPACEDIM], total_area = 0.;
195 double total_barycenter[SPACEDIM]={0.,0.};
197 const ConnType nbNodesT=targetCell.size()/SPACEDIM;
198 for(ConnType iT = 1; iT<nbNodesT-1; iT++)
200 std::vector<double> inter;
201 INTERP_KERNEL::intersec_de_triangle(&targetCell[0],&targetCell[SPACEDIM*iT],&targetCell[SPACEDIM*(iT+1)],
202 sourceCell[0], sourceCell[1], sourceCell[2],
203 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
204 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
205 ConnType nb_inter=((ConnType)inter.size())/2;
206 if(nb_inter >3) inter=reconstruct_polygon(inter);
207 for(ConnType i = 1; i<nb_inter-1; i++)
209 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],inter_area);
210 inter_area[0] = 0.5 * fabs( inter_area[0] );
211 total_area += inter_area[0];
212 std::vector<double> inter_bary=INTERP_KERNEL::bary_poly(inter);
213 total_barycenter[0] += inter_area[0] * inter_bary[0];
214 total_barycenter[1] += inter_area[0] * inter_bary[1];
217 if ( total_area > std::numeric_limits<double>::min() )
219 total_barycenter[0] /= total_area;
220 total_barycenter[1] /= total_area;
222 barycentric_coords( sourceCell, &total_barycenter[0], &res[0]);
223 res[0] *= total_area;
224 res[1] *= total_area;
225 res[2] *= total_area;