]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/INTERP_KERNEL/Geometric2DIntersector.txx
Salome HOME
786f41c5631d1be652ad4c65a62f9db333d2e197
[tools/medcoupling.git] / src / INTERP_KERNEL / Geometric2DIntersector.txx
1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D
2 //
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.
7 //
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.
12 //
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
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 // Author : Anthony Geay (CEA/DEN)
20 #ifndef __GEOMETRIC2DINTERSECTOR_TXX__
21 #define __GEOMETRIC2DINTERSECTOR_TXX__
22
23 #include "Geometric2DIntersector.hxx"
24 #include "PlanarIntersectorP0P0.txx"
25 #include "Planar2D1DIntersectorP0P0.txx"
26 #include "PlanarIntersectorP0P1.txx"
27 #include "PlanarIntersectorP1P0.txx"
28 #include "PlanarIntersectorP1P1.txx"
29 #include "PlanarIntersectorP1P0Bary.txx"
30 #include "PlanarIntersectorP0P1Bary.txx"
31 #include "CellModel.hxx"
32
33 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
34 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
35 #include "InterpKernelGeo2DEdgeLin.hxx"
36 #include "InterpKernelGeo2DNode.hxx"
37
38 #define GEO2D_INTERSECTOR    Geometric2DIntersector<MyMeshType,MyMatrix,InterpType>
39 #define INTERSECTOR_TEMPLATE template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
40
41 namespace INTERP_KERNEL
42 {
43   INTERSECTOR_TEMPLATE
44   GEO2D_INTERSECTOR::Geometric2DIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
45                                             double dimCaracteristic, double md3DSurf, double minDot3DSurf, double medianPlane,
46                                             double precision, int orientation):
47     InterpType<MyMeshType,MyMatrix,GEO2D_INTERSECTOR >(meshT,meshS,dimCaracteristic, precision, md3DSurf, minDot3DSurf, medianPlane, true, orientation, 0)
48   {
49     QUADRATIC_PLANAR::_precision=precision;
50   }
51   
52   INTERSECTOR_TEMPLATE
53   double GEO2D_INTERSECTOR::intersectGeometry(ConnType icellT,   ConnType icellS,
54                                               ConnType nbNodesT, ConnType nbNodesS)
55   {
56     int orientation = 1;
57     std::vector<double> CoordsT;
58     std::vector<double> CoordsS;
59     PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(icellT,icellS,nbNodesT,nbNodesS,CoordsT,CoordsS,orientation);
60     NormalizedCellType tT=PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getTypeOfElement(icellT);
61     NormalizedCellType tS=PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getTypeOfElement(icellS);
62     QuadraticPolygon *p1=buildPolygonFrom(CoordsT,tT);
63     QuadraticPolygon *p2=buildPolygonFrom(CoordsS,tS);
64     double ret=p1->intersectWithAbs(*p2);
65     delete p1; delete p2;
66     return orientation*ret;
67   }
68
69   INTERSECTOR_TEMPLATE
70   double GEO2D_INTERSECTOR::intersectGeometry1D(ConnType icellT,   ConnType icellS,
71                                                 ConnType nbNodesT, ConnType nbNodesS,
72                                                 bool& isColinear)
73   {
74     int orientation = 1;
75     std::vector<double> CoordsT;
76     std::vector<double> CoordsS;
77     PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(icellT,icellS,nbNodesT,nbNodesS,CoordsT,CoordsS,orientation);
78     NormalizedCellType tT=PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getTypeOfElement(icellT);
79     NormalizedCellType tS=PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getTypeOfElement(icellS);
80     QuadraticPolygon *p1=buildPolygonFrom(CoordsT,tT);
81     QuadraticPolygon *p2=buildPolygonOfOneEdgeFrom(CoordsS,tS);
82     double ret=p1->intersectWithAbs1D(*p2, isColinear);
83     delete p1; delete p2;
84     return orientation*ret;
85   }
86
87   INTERSECTOR_TEMPLATE
88   double GEO2D_INTERSECTOR::intersectGeometryWithQuadrangle(const double             * quadrangle,
89                                                             const std::vector<double>& sourceCoords,
90                                                             bool                       isSourceQuad)
91   {
92     std::vector<Node *> nodes(4);
93     nodes[0]=new Node(quadrangle[0],quadrangle[1]);
94     nodes[1]=new Node(quadrangle[SPACEDIM],quadrangle[SPACEDIM+1]);
95     nodes[2]=new Node(quadrangle[2*SPACEDIM],quadrangle[2*SPACEDIM+1]);
96     nodes[3]=new Node(quadrangle[3*SPACEDIM],quadrangle[3*SPACEDIM+1]);
97     int nbOfSourceNodes=sourceCoords.size()/SPACEDIM;
98     std::vector<Node *> nodes2(nbOfSourceNodes);
99     for(int i=0;i<nbOfSourceNodes;i++)
100       nodes2[i]=new Node(sourceCoords[i*SPACEDIM],sourceCoords[i*SPACEDIM+1]);
101     QuadraticPolygon *p1=QuadraticPolygon::BuildLinearPolygon(nodes);
102     QuadraticPolygon *p2;
103     if(!isSourceQuad)
104       p2=QuadraticPolygon::BuildLinearPolygon(nodes2);
105     else
106       p2=QuadraticPolygon::BuildArcCirclePolygon(nodes2);
107     double ret=p1->intersectWithAbs(*p2);
108     delete p1; delete p2;
109     return ret;
110   }
111
112   INTERSECTOR_TEMPLATE
113   double GEO2D_INTERSECTOR::intersectGeometryGeneral(const std::vector<double>& targetCoords,
114                                                      const std::vector<double>& sourceCoords)
115   {
116     int nbOfTargetNodes=targetCoords.size()/SPACEDIM;
117     std::vector<Node *> nodes(nbOfTargetNodes);
118     for(int i=0;i<nbOfTargetNodes;i++)
119       nodes[i]=new Node(targetCoords[i*SPACEDIM],targetCoords[i*SPACEDIM+1]);
120     int nbOfSourceNodes=sourceCoords.size()/SPACEDIM;
121     std::vector<Node *> nodes2(nbOfSourceNodes);
122     for(int i=0;i<nbOfSourceNodes;i++)
123       nodes2[i]=new Node(sourceCoords[i*SPACEDIM],sourceCoords[i*SPACEDIM+1]);
124     QuadraticPolygon *p1=QuadraticPolygon::BuildLinearPolygon(nodes);
125     QuadraticPolygon *p2=QuadraticPolygon::BuildLinearPolygon(nodes2);
126     double ret=p1->intersectWithAbs(*p2);
127     delete p1; delete p2;
128     return ret;
129   }
130
131   //================================================================================
132   /*!
133    * \brief Intersect a triangle and a polygon for P1P0 barycentric algorithm
134    *  \param targetCell - list of coordinates of target polygon in full interlace
135    *  \param targetCellQuadratic - specifies if target polygon is quadratic or not
136    *  \param sourceTria - list of coordinates of source triangle
137    *  \param res - coefficients a,b and c associated to nodes of sourceTria
138    */
139   //================================================================================
140
141   INTERSECTOR_TEMPLATE
142   double GEO2D_INTERSECTOR::intersectGeoBary(const std::vector<double>& targetCell,
143                                              bool                       targetCellQuadratic,
144                                              const double *             sourceTria,
145                                              std::vector<double>&       res)
146   {
147     std::vector<Node *> nodes(3);
148     nodes[0]=new Node(sourceTria[0*SPACEDIM],sourceTria[0*SPACEDIM+1]);
149     nodes[1]=new Node(sourceTria[1*SPACEDIM],sourceTria[1*SPACEDIM+1]);
150     nodes[2]=new Node(sourceTria[2*SPACEDIM],sourceTria[2*SPACEDIM+1]);
151     int nbOfTargetNodes=targetCell.size()/SPACEDIM;
152     std::vector<Node *> nodes2(nbOfTargetNodes);
153     for(int i=0;i<nbOfTargetNodes;i++)
154       nodes2[i]=new Node(targetCell[i*SPACEDIM],targetCell[i*SPACEDIM+1]);
155     QuadraticPolygon *p1=QuadraticPolygon::BuildLinearPolygon(nodes);
156     QuadraticPolygon *p2;
157     if(!targetCellQuadratic)
158       p2=QuadraticPolygon::BuildLinearPolygon(nodes2);
159     else
160       p2=QuadraticPolygon::BuildArcCirclePolygon(nodes2);
161     double barycenter[2];
162     double ret=p1->intersectWithAbs(*p2,barycenter);
163     delete p1; delete p2;
164     if ( ret > std::numeric_limits<double>::min() )
165     {
166       std::vector<const double* > sourceCell(3);
167       sourceCell[0] = &sourceTria[0];
168       sourceCell[1] = &sourceTria[SPACEDIM];
169       sourceCell[2] = &sourceTria[SPACEDIM*2];
170       res.resize(3);
171       barycentric_coords( sourceCell, barycenter, &res[0]);
172       res[0] *= ret;
173       res[1] *= ret;
174       res[2] *= ret;
175     }
176     else
177     {
178       ret = 0;
179     }
180     return ret;
181   }
182
183   INTERSECTOR_TEMPLATE
184   QuadraticPolygon *GEO2D_INTERSECTOR::buildPolygonFrom(const std::vector<double>& coords, NormalizedCellType type)
185   {
186     int nbNodes=coords.size()/SPACEDIM;
187     std::vector<Node *> nodes(nbNodes);
188     for(int i=0;i<nbNodes;i++)
189       nodes[i]=new Node(coords[i*SPACEDIM],coords[i*SPACEDIM+1]);
190     if(!CellModel::GetCellModel(type).isQuadratic())
191       return QuadraticPolygon::BuildLinearPolygon(nodes);
192     else
193       return QuadraticPolygon::BuildArcCirclePolygon(nodes);
194   }
195
196   INTERSECTOR_TEMPLATE
197   QuadraticPolygon *GEO2D_INTERSECTOR::buildPolygonOfOneEdgeFrom(const std::vector<double>& coords, NormalizedCellType type)
198   {
199     if(type==NORM_SEG2)
200       {
201         Node *node0=new Node(coords[0],coords[1]);
202         Node *node1=new Node(coords[SPACEDIM],coords[SPACEDIM+1]);
203         QuadraticPolygon *ret=new QuadraticPolygon;
204         ret->pushBack(new EdgeLin(node0,node1));
205         node0->decrRef(); node1->decrRef();
206         return ret;
207       }
208     else if(type==NORM_SEG3)
209       {
210         Node *nodeBg=new Node(coords[0],coords[1]);
211         Node *nodeEnd=new Node(coords[SPACEDIM],coords[SPACEDIM+1]);
212         Node *nodeMiddle=new Node(coords[2*SPACEDIM],coords[2*SPACEDIM+1]);
213         QuadraticPolygon *ret=new QuadraticPolygon;
214         ret->pushBack(new EdgeArcCircle(nodeBg,nodeMiddle,nodeEnd));
215         nodeBg->decrRef(); nodeEnd->decrRef(); nodeMiddle->decrRef();
216         return ret;
217       }
218     else
219       throw INTERP_KERNEL::Exception("buildPolygonOfOneEdgeFrom : trying to build such non close QuadraticPolygon with 1D type !");
220   }
221
222   INTERSECTOR_TEMPLATE
223   QuadraticPolygon *GEO2D_INTERSECTOR::buildPolygonAFrom(ConnType cell, int nbOfPoints, NormalizedCellType type)
224   {
225     const ConnType *startOfCellNodeConn=PlanarIntersector<MyMeshType,MyMatrix>::_connectT+OTT<ConnType,numPol>::conn2C(PlanarIntersector<MyMeshType,MyMatrix>::_connIndexT[OTT<ConnType,numPol>::ind2C(cell)]);
226     std::vector<Node *> nodes(nbOfPoints);
227     for(int i=0;i<nbOfPoints;i++)
228       nodes[i]=new Node(PlanarIntersector<MyMeshType,MyMatrix>::_coordsT+OTT<ConnType,numPol>::coo2C(startOfCellNodeConn[i])*SPACEDIM);
229     if(CellModel::GetCellModel(type).isQuadratic())
230       return QuadraticPolygon::BuildLinearPolygon(nodes);
231     else
232       return QuadraticPolygon::BuildArcCirclePolygon(nodes);
233   }
234
235   INTERSECTOR_TEMPLATE
236   QuadraticPolygon *GEO2D_INTERSECTOR::buildPolygonBFrom(ConnType cell, int nbOfPoints, NormalizedCellType type)
237   {
238     const ConnType *startOfCellNodeConn=PlanarIntersector<MyMeshType,MyMatrix>::_connectS+OTT<ConnType,numPol>::conn2C(PlanarIntersector<MyMeshType,MyMatrix>::_connIndexS[OTT<ConnType,numPol>::ind2C(cell)]);
239     std::vector<Node *> nodes(nbOfPoints);
240     for(int i=0;i<nbOfPoints;i++)
241       nodes[i]=new Node(PlanarIntersector<MyMeshType,MyMatrix>::_coordsS+OTT<ConnType,numPol>::coo2C(startOfCellNodeConn[i])*SPACEDIM);
242     const CellModel& cm=CellModel::GetCellModel(type);
243     if(!cm.isQuadratic())
244       return QuadraticPolygon::BuildLinearPolygon(nodes);
245     else
246       return QuadraticPolygon::BuildArcCirclePolygon(nodes);
247   }
248 }
249
250 #endif