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 __TRIANGULATIONINTERSECTOR_TXX__
20 #define __TRIANGULATIONINTERSECTOR_TXX__
22 #include "TriangulationIntersector.hxx"
23 #include "PlanarIntersectorP0P0.txx"
24 #include "PlanarIntersectorP0P1.txx"
25 #include "PlanarIntersectorP1P0.txx"
27 #include "InterpolationUtils.hxx"
28 #include "PlanarIntersector.hxx"
32 namespace INTERP_KERNEL
34 template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
35 TriangulationIntersector<MyMeshType,MyMatrix,InterpType>::TriangulationIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
36 double DimCaracteristic, double Precision,
37 double MedianPlane, int orientation, int PrintLevel)
38 :InterpType<MyMeshType,MyMatrix,TriangulationIntersector<MyMeshType,MyMatrix,InterpType> >(meshT,meshS,DimCaracteristic, Precision, MedianPlane, true, orientation, PrintLevel)
40 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 1)
42 std::cout << " - intersection type = triangles " << std::endl;
43 if(SPACEDIM==3) std::cout << "_do_rotate = true"<< std::endl;
47 template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
48 double TriangulationIntersector<MyMeshType,MyMatrix,InterpType>::intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS)
53 //Obtain the coordinates of T and S
54 std::vector<double> CoordsT;
55 std::vector<double> CoordsS;
56 PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(icellT,icellS,nbNodesT,nbNodesS,CoordsT,CoordsS,orientation);
57 //Compute the intersection area
58 double area[SPACEDIM];
59 for(ConnType iT = 1; iT<nbNodesT-1; iT++)
61 for(ConnType iS = 1; iS<nbNodesS-1; iS++)
63 std::vector<double> inter;
64 INTERP_KERNEL::intersec_de_triangle(&CoordsT[0],&CoordsT[SPACEDIM*iT],&CoordsT[SPACEDIM*(iT+1)],
65 &CoordsS[0],&CoordsS[SPACEDIM*iS],&CoordsS[SPACEDIM*(iS+1)],
66 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
67 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
68 ConnType nb_inter=((ConnType)inter.size())/2;
69 if(nb_inter >3) inter=reconstruct_polygon(inter);
70 for(ConnType i = 1; i<nb_inter-1; i++)
72 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
73 result +=0.5*fabs(area[0]);
76 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
78 std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl;
79 for(ConnType i=0; i< nb_inter; i++)
80 {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; }
86 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
87 std::cout << std::endl <<"Intersection area = " << result << std::endl;
89 return orientation*result;
92 template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
93 double TriangulationIntersector<MyMeshType,MyMatrix,InterpType>::intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector<double>& sourceCoords, bool isSourceQuad)
96 ConnType nbNodesS=sourceCoords.size()/SPACEDIM;
97 //Compute the intersection area
98 double area[SPACEDIM];
99 for(ConnType iT = 1; iT<3; iT++)
101 for(ConnType iS = 1; iS<nbNodesS-1; iS++)
103 std::vector<double> inter;
104 INTERP_KERNEL::intersec_de_triangle(quadrangle,&quadrangle[SPACEDIM*iT],&quadrangle[SPACEDIM*(iT+1)],
105 &sourceCoords[0],&sourceCoords[SPACEDIM*iS],&sourceCoords[SPACEDIM*(iS+1)],
106 inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
107 PlanarIntersector<MyMeshType,MyMatrix>::_precision);
108 ConnType nb_inter=((ConnType)inter.size())/2;
109 if(nb_inter >3) inter=reconstruct_polygon(inter);
110 for(ConnType i = 1; i<nb_inter-1; i++)
112 INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
113 result +=0.5*fabs(area[0]);
116 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
118 std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl;
119 for(ConnType i=0; i< nb_inter; i++)
120 {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; }
126 if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 3)
127 std::cout << std::endl <<"Intersection area = " << result << std::endl;