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 __PLANARINTERSECTOR_TXX__
20 #define __PLANARINTERSECTOR_TXX__
22 #include "PlanarIntersector.hxx"
23 #include "InterpolationUtils.hxx"
24 #include "TranslationRotationMatrix.hxx"
29 namespace INTERP_KERNEL
31 template<class MyMeshType, class MyMatrix>
32 PlanarIntersector<MyMeshType,MyMatrix>::PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel):
33 _meshT(meshT),_meshS(meshS),
34 _dim_caracteristic(dimCaracteristic),_precision(precision),_median_plane(medianPlane),
35 _do_rotate(doRotate),_orientation(orientation),_print_level(printLevel)
37 _connectT=meshT.getConnectivityPtr();
38 _connectS=meshS.getConnectivityPtr();
39 _connIndexT=meshT.getConnectivityIndexPtr();
40 _connIndexS=meshS.getConnectivityIndexPtr();
41 _coordsT=meshT.getCoordinatesPtr();
42 _coordsS=meshS.getCoordinatesPtr();
45 template<class MyMeshType, class MyMatrix>
46 PlanarIntersector<MyMeshType,MyMatrix>::~PlanarIntersector()
51 \brief creates the bounding boxes for all the cells of mesh \a mesh
53 The method accepts mixed meshes (containing triangles and quadrangles).
54 The vector returned is of dimension 6*nb_elems with bounding boxes stored as xmin1, xmax1, ymin1, ymax1, zmin1, zmax1, xmin2, xmax2, ymin2,...
55 The returned pointer must be deleted by the calling code.
57 \param mesh structure pointing to the mesh
58 \param bbox vector containing the bounding boxes
60 template<class MyMeshType, class MyMatrix>
61 void PlanarIntersector<MyMeshType,MyMatrix>::createBoundingBoxes(const MyMeshType& mesh, std::vector<double>& bbox)
63 /* We build the segment tree for locating possible matching intersections*/
64 long nbelems = mesh.getNumberOfElements();
65 bbox.resize(2*SPACEDIM* nbelems);
66 const double* coords = mesh.getCoordinatesPtr();
67 const ConnType* conn = mesh.getConnectivityPtr();
68 const ConnType* conn_index = mesh.getConnectivityIndexPtr();
70 for(long icell=0; icell<nbelems; icell++)
72 int nb_nodes_per_elem =conn_index[icell+1]-conn_index[icell];
73 //initializing bounding box limits
74 for(int idim=0; idim<SPACEDIM; idim++)
76 bbox[2*SPACEDIM*ibox+2*idim] = std::numeric_limits<double>::max();
77 bbox[2*SPACEDIM*ibox+2*idim+1] = -std::numeric_limits<double>::max();
79 //updating the bounding box with each node of the element
80 for (int j=0; j<nb_nodes_per_elem; j++)
82 const double* coord_node=coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(conn[OTT<ConnType,numPol>::conn2C(conn_index[icell]+j)]);
83 for(int idim=0; idim<SPACEDIM; idim++)
85 double x=*(coord_node+idim);
86 bbox[ibox*2*SPACEDIM + 2*idim] = (bbox[ibox*2*SPACEDIM + 2*idim] <x)?bbox[ibox*2*SPACEDIM + 2*idim ]:x;
87 bbox[ibox*2*SPACEDIM + 2*idim+1] = (bbox[ibox*2*SPACEDIM + 2*idim+1]>x)?bbox[ibox*2*SPACEDIM + 2*idim+1]:x;
95 Computes the bouding box of a given element. iP in numPol mode.
97 template<class MyMeshType, class MyMatrix>
98 void PlanarIntersector<MyMeshType,MyMatrix>::getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes)
100 const double* coords = mesh.getCoordinatesPtr();
101 const ConnType* conn_index = mesh.getConnectivityIndexPtr();
102 const ConnType* conn = mesh.getConnectivityPtr();
103 //initializing bounding box limits
104 for(int idim=0; idim<SPACEDIM; idim++)
106 bb[2*idim ] = std::numeric_limits<double>::max();
107 bb[2*idim+1] = -std::numeric_limits<double>::max();
110 for (ConnType i=0; i<nb_nodes; i++)
112 //MN: iP= cell index, not node index, use of connectivity array ?
113 const double* coord_node=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(conn[OTT<ConnType,numPol>::conn2C(conn_index[OTT<ConnType,numPol>::ind2C(iP)]+i)]));
114 for(int idim=0; idim<SPACEDIM; idim++)
116 double x = *(coord_node+idim);
117 //double y = *(mesh.getCoordinates(MED_FULL_INTERLACE)+3*(iP+i)+1);
118 bb[2*idim ] = (x<bb[2*idim ])?x:bb[2*idim ];
119 bb[2*idim+1] = (x>bb[2*idim+1])?x:bb[2*idim+1];
124 /*! Readjusts a set of bounding boxes so that they are extended
125 in all dimensions for avoiding missing interesting intersections
127 \param bbox vector containing the bounding boxes
129 template<class MyMeshType, class MyMatrix>
130 void PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double Surf3DAdjustmentEps)
132 /* We build the segment tree for locating possible matching intersections*/
134 long size = bbox.size()/(2*SPACEDIM);
135 for (int i=0; i<size; i++)
137 double max=- std::numeric_limits<double>::max();
138 for(int idim=0; idim<SPACEDIM; idim++)
140 double Dx=bbox[i*2*SPACEDIM+1+2*idim]-bbox[i*2*SPACEDIM+2*idim];
143 for(int idim=0; idim<SPACEDIM; idim++)
145 bbox[i*2*SPACEDIM+2*idim ] -= Surf3DAdjustmentEps*max;
146 bbox[i*2*SPACEDIM+2*idim+1] += Surf3DAdjustmentEps*max;
152 * @param icellT id in target mesh in format of MyMeshType.
153 * @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
155 template<class MyMeshType, class MyMatrix>
156 void PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT)
158 int nbNodesT=_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1]-_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
159 coordsT.resize(SPACEDIM*nbNodesT);
160 for (ConnType iT=0; iT<nbNodesT; iT++)
161 for(int idim=0; idim<SPACEDIM; idim++)
162 coordsT[SPACEDIM*iT+idim]=_coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
166 * @param icellS id in source mesh in format of MyMeshType.
167 * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
169 template<class MyMeshType, class MyMatrix>
170 void PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS)
172 int nbNodesS=_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1]-_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
173 coordsS.resize(SPACEDIM*nbNodesS);
174 for (ConnType iS=0; iS<nbNodesS; iS++)
175 for(int idim=0; idim<SPACEDIM; idim++)
176 coordsS[SPACEDIM*iS+idim]=_coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
180 * @param icellT id in target mesh in format of MyMeshType.
181 * @param icellS id in source mesh in format of MyMeshType.
182 * @param nbNodesT nb of nodes of the target cell.
183 * @param nbNodesS nb of nodes of the source cell.
184 * @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
185 * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
186 * @param orientation is an output value too, only set if SPACEDIM==3.
188 template<class MyMeshType, class MyMatrix>
189 void PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, std::vector<double>& coordsT, std::vector<double>& coordsS, int& orientation)
191 /*double epsilon=_precision*_dim_caracteristic;
192 coordsT.resize(SPACEDIM*nbNodesT);
193 coordsS.resize(SPACEDIM*nbNodesS);
194 int nb_dist_NodesT=nbNodesT;
195 int nb_dist_NodesS=nbNodesS;
196 int i_last = nbNodesT - 1;
197 const double * Pi_last=_coordsT +_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+i_last)];
199 for (int iT=0; iT<nbNodesT; iT++)
201 const double * PiT = _coordsT + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)]);
202 if(distance2<SPACEDIM>(Pi_last, PiT)>epsilon)
204 for (int idim=0; idim<SPACEDIM; idim++)
205 coordsT[SPACEDIM*iT+idim]=PiT[idim];
206 i_last=iT; Pi_last = PiT;
211 coordsT.resize(nb_dist_NodesT*SPACEDIM);
212 i_last = nbNodesS - 1;
213 Pi_last=_coordsS + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+i_last)]);
214 for (int iS=0; iS<nbNodesS; iS++)
216 const double * PiS=_coordsS+SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)]);
217 if(distance2<SPACEDIM>(Pi_last, PiS)>epsilon)
219 for (int idim=0; idim<SPACEDIM; idim++)
220 coordsS[SPACEDIM*iS+idim]=PiS[idim];
221 i_last=iS; Pi_last = PiS;
226 coordsS.resize(nb_dist_NodesS*SPACEDIM);
227 //project cells T and S on the median plane
228 // and rotate the median plane
230 orientation = projectionThis(&coordsT[0], &coordsS[0], nb_dist_NodesT, nb_dist_NodesS);
233 if(_print_level >= 3)
235 std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl;
236 std::cout << std::endl << "icellT= " << icellT << ", nb nodes T= " << nbNodesT << std::endl;
237 for(int iT =0; iT< nbNodesT; iT++)
238 {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsT[SPACEDIM*iT+idim] << " "; std::cout << std::endl;}
239 std::cout << std::endl << "icellS= " << icellS << ", nb nodes S= " << nbNodesS << std::endl;
240 for(int iS =0; iS< nbNodesS; iS++)
241 {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsS[SPACEDIM*iS+idim]<< " "; std::cout << std::endl; }
243 coordsT.resize(SPACEDIM*nbNodesT);
244 coordsS.resize(SPACEDIM*nbNodesS);
245 for (int idim=0; idim<SPACEDIM; idim++)
247 for (ConnType iT=0; iT<nbNodesT; iT++)
248 coordsT[SPACEDIM*iT+idim] = _coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
249 for (ConnType iS=0; iS<nbNodesS; iS++)
250 coordsS[SPACEDIM*iS+idim] = _coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
253 //project cells T and S on the median plane and rotate the median plane
255 orientation = projectionThis(&coordsT[0], &coordsS[0], nbNodesT, nbNodesS);
258 if(_print_level >= 3)
260 std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl;
261 std::cout << std::endl << "icellT= " << icellT << ", nb nodes T= " << nbNodesT << std::endl;
262 for(int iT =0; iT< nbNodesT; iT++)
263 {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsT[SPACEDIM*iT+idim] << " "; std::cout << std::endl;}
264 std::cout << std::endl << "icellS= " << icellS << ", nb nodes S= " << nbNodesS << std::endl;
265 for(int iS =0; iS< nbNodesS; iS++)
266 {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsS[SPACEDIM*iS+idim]<< " "; std::cout << std::endl; }
270 template<class MyMeshType, class MyMatrix>
271 int PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB)
273 return projection(Coords_A,Coords_B,nb_NodesA,nb_NodesB,_dim_caracteristic*_precision,_median_plane,_do_rotate);
276 template<class MyMeshType, class MyMatrix>
277 int PlanarIntersector<MyMeshType,MyMatrix>::projection(double *Coords_A, double *Coords_B,
278 int nb_NodesA, int nb_NodesB, double epsilon, double median_plane, bool do_rotate)
280 double normal_A[3]={0,0,0};
281 double normal_B[3]={0,0,0};
282 double linear_comb[3];
284 bool same_orientation;
286 //Find the normal to cells A and B
288 while(i_A1<nb_NodesA && distance2<SPACEDIM>(Coords_A,&Coords_A[SPACEDIM*i_A1])< epsilon) i_A1++;
290 crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
291 double normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
292 while(i_A2<nb_NodesA && normA < epsilon)
294 crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
296 normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
300 while(i_B1<nb_NodesB && distance2<SPACEDIM>(Coords_B,Coords_B+SPACEDIM*i_B1)< epsilon) i_B1++;
302 crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
303 double normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
304 while(i_B2<nb_NodesB && normB < epsilon)
306 crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
308 normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
311 if(i_A2<nb_NodesA && i_B2<nb_NodesB)
313 //Build the normal of the median plane
314 same_orientation = dotprod<SPACEDIM>(normal_A,normal_B)>=0;
316 if(!same_orientation)
317 for(int idim =0; idim< SPACEDIM; idim++) normal_A[idim] *=-1;
319 double normB= sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
321 for(int idim =0; idim< SPACEDIM; idim++)
322 linear_comb[idim] = median_plane*normal_A[idim]/normA + (1-median_plane)*normal_B[idim]/normB;
323 double norm= sqrt(dotprod<SPACEDIM>(linear_comb,linear_comb));
325 //Necessarily: norm>epsilon, no need to check
326 for(int idim =0; idim< SPACEDIM; idim++) linear_comb[idim]/=norm;
328 //Project the nodes of A and B on the median plane
329 for(int i_A=0; i_A<nb_NodesA; i_A++)
331 proj = dotprod<SPACEDIM>(&Coords_A[SPACEDIM*i_A],linear_comb);
332 for(int idim =0; idim< SPACEDIM; idim++)
333 Coords_A[SPACEDIM*i_A+idim] -= proj*linear_comb[idim];
335 for(int i_B=0; i_B<nb_NodesB; i_B++)
337 proj = dotprod<SPACEDIM>(Coords_B+SPACEDIM*i_B,linear_comb);
338 for(int idim =0; idim< SPACEDIM; idim++)
339 Coords_B[SPACEDIM*i_B+idim] -= proj*linear_comb[idim];
342 //Buid the matrix sending A into the Oxy plane and apply it to A and B
345 TranslationRotationMatrix rotation;
346 //rotate3DTriangle(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2], rotation);
347 rotate3DTriangle(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2, rotation);
348 for (int i=0; i<nb_NodesA; i++) rotation.transform_vector(Coords_A+SPACEDIM*i);
349 for (int i=0; i<nb_NodesB; i++) rotation.transform_vector(Coords_B+SPACEDIM*i);
357 std::cout << " Maille dégénérée " << "epsilon = " << epsilon << std::endl;
358 std::cout << " i_A1= " << i_A1 << " i_A2= " << i_A2 << std::endl;
359 std::cout << " distance2<SPACEDIM>(Coords_A,&Coords_A[i_A1])= " << distance2<SPACEDIM>(Coords_A,&Coords_A[i_A1]) << std::endl;
360 std::cout << "abs(normal_A) = " << fabs(normal_A[0]) << " ; " <<fabs( normal_A[1]) << " ; " << fabs(normal_A[2]) << std::endl;
361 std::cout << " i_B1= " << i_B1 << " i_B2= " << i_B2 << std::endl;
362 std::cout << " distance2<SPACEDIM>(&Coords_B[0],&Coords_B[i_B1])= " << distance2<SPACEDIM>(Coords_B,Coords_B+i_B1) << std::endl;
363 std::cout << "normal_B = " << normal_B[0] << " ; " << normal_B[1] << " ; " << normal_B[2] << std::endl;
369 template<class MyMeshType, class MyMatrix>
370 void PlanarIntersector<MyMeshType,MyMatrix>::rotate3DTriangle(double* PP1, double*PP2, double*PP3,
371 TranslationRotationMatrix& rotation_matrix)
374 rotation_matrix.translate(PP1);
378 P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2];
379 P3w[0]=PP3[0]; P3w[1]=PP3[1];P3w[2]=PP3[2];
381 // translating to set P1 at the origin
382 for (int i=0; i<3; i++)
388 // rotating to set P2 on the Oxy plane
389 TranslationRotationMatrix A;
391 A.rotate_vector(P3w);
392 rotation_matrix.multiply(A);
394 //rotating to set P2 on the Ox axis
395 TranslationRotationMatrix B;
397 B.rotate_vector(P3w);
398 rotation_matrix.multiply(B);
400 //rotating to set P3 on the Oxy plane
401 TranslationRotationMatrix C;
403 rotation_matrix.multiply(C);