1 // Copyright (C) 2007-2012 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 __TRIANGULATIONINTERSECTOR_TXX__
20 #define __TRIANGULATIONINTERSECTOR_TXX__
22 #include "TriangulationIntersector.hxx"
23 #include "PlanarIntersectorP0P0.txx"
24 #include "PlanarIntersectorP0P1.txx"
25 #include "PlanarIntersectorP1P0.txx"
26 #include "PlanarIntersectorP1P0Bary.txx"
28 #include "InterpolationUtils.hxx"
29 #include "PlanarIntersector.hxx"
33 #define TRI_INTERSECTOR TriangulationIntersector<MyMeshType,MyMatrix,InterpType>
34 #define TRI_INTER_TEMPLATE template<class MyMeshType, class MyMatrix, \
35 template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
37 namespace INTERP_KERNEL
40 TRI_INTERSECTOR::TriangulationIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
41 double DimCaracteristic, double Precision, double md3DSurf,
42 double MedianPlane, int orientation, int PrintLevel)
43 :InterpType<MyMeshType,MyMatrix,TRI_INTERSECTOR >(meshT,meshS,DimCaracteristic, Precision, md3DSurf,
44 MedianPlane, true, orientation, PrintLevel)
46 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 1)
48 std::cout << " - intersection type = triangles " << std::endl;
49 if(SPACEDIM==3) std::cout << "_do_rotate = true"<< std::endl;
54 double TRI_INTERSECTOR::intersectGeometry(ConnType icellT, ConnType icellS,
55 ConnType nbNodesT, ConnType nbNodesS)
60 //Obtain the coordinates of T and S
61 std::vector<double> CoordsT;
62 std::vector<double> CoordsS;
63 PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(icellT,icellS,nbNodesT,nbNodesS,CoordsT,CoordsS,orientation);
64 //Compute the intersection area
65 double area[SPACEDIM];
66 for(ConnType iT = 1; iT<nbNodesT-1; iT++)
68 for(ConnType iS = 1; iS<nbNodesS-1; iS++)
70 std::vector<double> inter;
71 INTERP_KERNEL::intersec_de_triangle(&CoordsT[0],&CoordsT[SPACEDIM*iT],&CoordsT[SPACEDIM*(iT+1)],
72 &CoordsS[0],&CoordsS[SPACEDIM*iS],&CoordsS[SPACEDIM*(iS+1)],
73 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
74 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
75 ConnType nb_inter=((ConnType)inter.size())/2;
76 if(nb_inter >3) inter=reconstruct_polygon(inter);
77 for(ConnType i = 1; i<nb_inter-1; i++)
79 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
80 result +=0.5*fabs(area[0]);
83 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
85 std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl;
86 for(ConnType i=0; i< nb_inter; i++)
87 {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; }
93 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
94 std::cout << std::endl <<"Intersection area = " << result << std::endl;
96 return orientation*result;
100 double TRI_INTERSECTOR::intersectGeometryWithQuadrangle(const double * quadrangle,
101 const std::vector<double>& sourceCoords,
105 ConnType nbNodesS=sourceCoords.size()/SPACEDIM;
106 //Compute the intersection area
107 double area[SPACEDIM];
108 for(ConnType iT = 1; iT<3; iT++)
110 for(ConnType iS = 1; iS<nbNodesS-1; iS++)
112 std::vector<double> inter;
113 INTERP_KERNEL::intersec_de_triangle(quadrangle,&quadrangle[SPACEDIM*iT],&quadrangle[SPACEDIM*(iT+1)],
114 &sourceCoords[0],&sourceCoords[SPACEDIM*iS],&sourceCoords[SPACEDIM*(iS+1)],
115 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
116 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
117 ConnType nb_inter=((ConnType)inter.size())/2;
118 if(nb_inter >3) inter=reconstruct_polygon(inter);
119 for(ConnType i = 1; i<nb_inter-1; i++)
121 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
122 result +=0.5*fabs(area[0]);
125 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
127 std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl;
128 for(ConnType i=0; i< nb_inter; i++)
129 {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; }
135 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
136 std::cout << std::endl <<"Intersection area = " << result << std::endl;
142 double TRI_INTERSECTOR::intersectGeometryGeneral(const std::vector<double>& targetCoords,
143 const std::vector<double>& sourceCoords)
146 ConnType nbNodesS=sourceCoords.size()/SPACEDIM;
147 ConnType nbNodesT=targetCoords.size()/SPACEDIM;
148 //Compute the intersection area
149 double area[SPACEDIM];
150 for(ConnType iT = 1; iT<nbNodesT-1; iT++)
152 for(ConnType iS = 1; iS<nbNodesS-1; iS++)
154 std::vector<double> inter;
155 INTERP_KERNEL::intersec_de_triangle(&targetCoords[0],&targetCoords[SPACEDIM*iT],&targetCoords[SPACEDIM*(iT+1)],
156 &sourceCoords[0],&sourceCoords[SPACEDIM*iS],&sourceCoords[SPACEDIM*(iS+1)],
157 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
158 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
159 ConnType nb_inter=((ConnType)inter.size())/2;
160 if(nb_inter >3) inter=reconstruct_polygon(inter);
161 for(ConnType i = 1; i<nb_inter-1; i++)
163 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
164 result +=0.5*fabs(area[0]);
171 //================================================================================
173 * \brief Intersect a triangle and a polygon for P1P0 barycentric algorithm
174 * \param targetCell - list of coordinates of target polygon in full interlace
175 * \param targetCellQuadratic - specifies if target polygon is quadratic or not
176 * \param sourceTria - list of coordinates of source triangle
177 * \param res - coefficients a,b and c associated to nodes of sourceTria
179 //================================================================================
182 double TRI_INTERSECTOR::intersectGeoBary(const std::vector<double>& targetCell,
183 bool targetCellQuadratic,
184 const double * sourceTria,
185 std::vector<double>& res)
187 std::vector<const double* > sourceCell(3);
188 sourceCell[0] = &sourceTria[0];
189 sourceCell[1] = &sourceTria[SPACEDIM];
190 sourceCell[2] = &sourceTria[SPACEDIM*2];
192 //Compute the intersection area
193 double inter_area[SPACEDIM], total_area = 0.;
194 double total_barycenter[SPACEDIM]={0.,0.};
196 const ConnType nbNodesT=targetCell.size()/SPACEDIM;
197 for(ConnType iT = 1; iT<nbNodesT-1; iT++)
199 std::vector<double> inter;
200 INTERP_KERNEL::intersec_de_triangle(&targetCell[0],&targetCell[SPACEDIM*iT],&targetCell[SPACEDIM*(iT+1)],
201 sourceCell[0], sourceCell[1], sourceCell[2],
202 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
203 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
204 ConnType nb_inter=((ConnType)inter.size())/2;
205 if(nb_inter >3) inter=reconstruct_polygon(inter);
206 for(ConnType i = 1; i<nb_inter-1; i++)
208 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],inter_area);
209 inter_area[0] = 0.5 * fabs( inter_area[0] );
210 total_area += inter_area[0];
211 std::vector<double> inter_bary=INTERP_KERNEL::bary_poly(inter);
212 total_barycenter[0] += inter_area[0] * inter_bary[0];
213 total_barycenter[1] += inter_area[0] * inter_bary[1];
216 if ( total_area > std::numeric_limits<double>::min() )
218 total_barycenter[0] /= total_area;
219 total_barycenter[1] /= total_area;
221 barycentric_coords( sourceCell, &total_barycenter[0], &res[0]);
222 res[0] *= total_area;
223 res[1] *= total_area;
224 res[2] *= total_area;