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 __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 md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel):
33 _meshT(meshT),_meshS(meshS),
34 _dim_caracteristic(dimCaracteristic),_max_distance_3Dsurf_intersect(md3DSurf),_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, double surf3DAdjustmentEpsAbs)
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+surf3DAdjustmentEpsAbs;
146 bbox[i*2*SPACEDIM+2*idim+1] += surf3DAdjustmentEps*max+surf3DAdjustmentEpsAbs;
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 offset is a value in C format that indicates the number of circular permutation.
182 * @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
184 template<class MyMeshType, class MyMatrix>
185 void PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinatesPermute(ConnType icellT, int offset, std::vector<double>& coordsT)
187 int nbNodesT=_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1]-_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
188 coordsT.resize(SPACEDIM*nbNodesT);
189 for (ConnType iTTmp=0; iTTmp<nbNodesT; iTTmp++)
191 ConnType iT=(iTTmp+offset)%nbNodesT;
192 for(int idim=0; idim<SPACEDIM; idim++)
193 coordsT[SPACEDIM*iTTmp+idim]=_coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
198 * @param icellS id in source mesh in format of MyMeshType.
199 * @param offset is a value in C format that indicates the number of circular permutation.
200 * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
202 template<class MyMeshType, class MyMatrix>
203 void PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinatesPermute(ConnType icellS, int offset, std::vector<double>& coordsS)
205 int nbNodesS=_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1]-_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
206 coordsS.resize(SPACEDIM*nbNodesS);
207 for (ConnType iSTmp=0; iSTmp<nbNodesS; iSTmp++)
209 ConnType iS=(iSTmp+offset)%nbNodesS;
210 for(int idim=0; idim<SPACEDIM; idim++)
211 coordsS[SPACEDIM*iSTmp+idim]=_coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
216 * @param icellT id in target mesh in format of MyMeshType.
217 * @param icellS id in source mesh in format of MyMeshType.
218 * @param nbNodesT nb of nodes of the target cell.
219 * @param nbNodesS nb of nodes of the source cell.
220 * @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
221 * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
222 * @param orientation is an output value too, only set if SPACEDIM==3.
224 template<class MyMeshType, class MyMatrix>
225 void PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, std::vector<double>& coordsT, std::vector<double>& coordsS, int& orientation)
227 coordsT.resize(SPACEDIM*nbNodesT);
228 coordsS.resize(SPACEDIM*nbNodesS);
229 for (int idim=0; idim<SPACEDIM; idim++)
231 for (ConnType iT=0; iT<nbNodesT; iT++)
232 coordsT[SPACEDIM*iT+idim] = _coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
233 for (ConnType iS=0; iS<nbNodesS; iS++)
234 coordsS[SPACEDIM*iS+idim] = _coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
237 //project cells T and S on the median plane and rotate the median plane
239 orientation = projectionThis(&coordsT[0], &coordsS[0], nbNodesT, nbNodesS);
242 if(_print_level >= 3)
244 std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl;
245 std::cout << std::endl << "icellT= " << icellT << ", nb nodes T= " << nbNodesT << std::endl;
246 for(int iT =0; iT< nbNodesT; iT++)
247 {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsT[SPACEDIM*iT+idim] << " "; std::cout << std::endl;}
248 std::cout << std::endl << "icellS= " << icellS << ", nb nodes S= " << nbNodesS << std::endl;
249 for(int iS =0; iS< nbNodesS; iS++)
250 {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsS[SPACEDIM*iS+idim]<< " "; std::cout << std::endl; }
255 * Filtering out zero surfaces and badly oriented surfaces
256 * _orientation = -1,0,1,2
257 * -1 : the intersection is taken into account if target and cells have different orientation
258 * 0 : the intersection is always taken into account
259 * 1 : the intersection is taken into account if target and cells have the same orientation
260 * 2 : the absolute value of intersection is always taken into account
262 template<class MyMeshType, class MyMatrix>
263 double PlanarIntersector<MyMeshType,MyMatrix>::getValueRegardingOption(double val) const
269 if (( val > 0.0 && _orientation==1) || ( val < 0.0 && _orientation==-1 ))
270 return _orientation*val;
274 template<class MyMeshType, class MyMatrix>
275 int PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB)
277 return projection(Coords_A,Coords_B,nb_NodesA,nb_NodesB,_dim_caracteristic*_precision,_max_distance_3Dsurf_intersect,_median_plane,_do_rotate);
280 template<class MyMeshType, class MyMatrix>
281 int PlanarIntersector<MyMeshType,MyMatrix>::projection(double *Coords_A, double *Coords_B,
282 int nb_NodesA, int nb_NodesB, double epsilon, double md3DSurf, double median_plane, bool do_rotate)
284 double normal_A[3]={0,0,0};
285 double normal_B[3]={0,0,0};
286 double linear_comb[3];
288 bool same_orientation;
290 //Find the normal to cells A and B
292 while(i_A1<nb_NodesA && distance2<SPACEDIM>(Coords_A,&Coords_A[SPACEDIM*i_A1])< epsilon) i_A1++;
294 crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
295 double normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
296 while(i_A2<nb_NodesA && normA < epsilon)
298 crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
300 normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
304 while(i_B1<nb_NodesB && distance2<SPACEDIM>(Coords_B,Coords_B+SPACEDIM*i_B1)< epsilon) i_B1++;
306 crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
307 double normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
308 while(i_B2<nb_NodesB && normB < epsilon)
310 crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
312 normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
319 for (int i=0;i<3;i++)
322 for (int j=0;j<nb_NodesA;j++)
323 coords_GA[i]+=Coords_A[3*j+i];
324 coords_GA[i]/=nb_NodesA;
326 double G1[3],G2[3],G3[3];
327 for (int i=0;i<3;i++)
329 G1[i]=Coords_B[i]-coords_GA[i];
330 G2[i]=Coords_B[i+3]-coords_GA[i];
331 G3[i]=Coords_B[i+6]-coords_GA[i];
334 prodvect[0]=G1[1]*G2[2]-G1[2]*G2[1];
335 prodvect[1]=G1[2]*G2[0]-G1[0]*G2[2];
336 prodvect[2]=G1[0]*G2[1]-G1[1]*G2[0];
337 double prodscal=prodvect[0]*G3[0]+prodvect[1]*G3[1]+prodvect[2]*G3[2];
338 if(fabs(prodscal)>md3DSurf)
341 if(i_A2<nb_NodesA && i_B2<nb_NodesB)
343 //Build the normal of the median plane
344 same_orientation = dotprod<SPACEDIM>(normal_A,normal_B)>=0;
346 if(!same_orientation)
347 for(int idim =0; idim< SPACEDIM; idim++) normal_A[idim] *=-1;
349 double normBB= sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
351 for(int idim =0; idim< SPACEDIM; idim++)
352 linear_comb[idim] = median_plane*normal_A[idim]/normA + (1-median_plane)*normal_B[idim]/normBB;
353 double norm= sqrt(dotprod<SPACEDIM>(linear_comb,linear_comb));
355 //Necessarily: norm>epsilon, no need to check
356 for(int idim =0; idim< SPACEDIM; idim++) linear_comb[idim]/=norm;
358 //Project the nodes of A and B on the median plane
359 for(int i_A=0; i_A<nb_NodesA; i_A++)
361 proj = dotprod<SPACEDIM>(&Coords_A[SPACEDIM*i_A],linear_comb);
362 for(int idim =0; idim< SPACEDIM; idim++)
363 Coords_A[SPACEDIM*i_A+idim] -= proj*linear_comb[idim];
365 for(int i_B=0; i_B<nb_NodesB; i_B++)
367 proj = dotprod<SPACEDIM>(Coords_B+SPACEDIM*i_B,linear_comb);
368 for(int idim =0; idim< SPACEDIM; idim++)
369 Coords_B[SPACEDIM*i_B+idim] -= proj*linear_comb[idim];
372 //Buid the matrix sending A into the Oxy plane and apply it to A and B
375 TranslationRotationMatrix rotation;
376 //rotate3DTriangle(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2], rotation);
377 rotate3DTriangle(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2, rotation);
378 for (int i=0; i<nb_NodesA; i++) rotation.transform_vector(Coords_A+SPACEDIM*i);
379 for (int i=0; i<nb_NodesB; i++) rotation.transform_vector(Coords_B+SPACEDIM*i);
387 std::cout << " Degenerated cell " << "epsilon = " << epsilon << std::endl;
388 std::cout << " i_A1= " << i_A1 << " i_A2= " << i_A2 << std::endl;
389 std::cout << " distance2<SPACEDIM>(Coords_A,&Coords_A[i_A1])= " << distance2<SPACEDIM>(Coords_A,&Coords_A[i_A1]) << std::endl;
390 std::cout << "abs(normal_A) = " << fabs(normal_A[0]) << " ; " <<fabs( normal_A[1]) << " ; " << fabs(normal_A[2]) << std::endl;
391 std::cout << " i_B1= " << i_B1 << " i_B2= " << i_B2 << std::endl;
392 std::cout << " distance2<SPACEDIM>(&Coords_B[0],&Coords_B[i_B1])= " << distance2<SPACEDIM>(Coords_B,Coords_B+i_B1) << std::endl;
393 std::cout << "normal_B = " << normal_B[0] << " ; " << normal_B[1] << " ; " << normal_B[2] << std::endl;
399 template<class MyMeshType, class MyMatrix>
400 void PlanarIntersector<MyMeshType,MyMatrix>::rotate3DTriangle(double* PP1, double*PP2, double*PP3,
401 TranslationRotationMatrix& rotation_matrix)
404 rotation_matrix.translate(PP1);
408 P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2];
409 P3w[0]=PP3[0]; P3w[1]=PP3[1];P3w[2]=PP3[2];
411 // translating to set P1 at the origin
412 for (int i=0; i<3; i++)
418 // rotating to set P2 on the Oxy plane
419 TranslationRotationMatrix A;
421 A.rotate_vector(P3w);
422 rotation_matrix.multiply(A);
424 //rotating to set P2 on the Ox axis
425 TranslationRotationMatrix B;
427 B.rotate_vector(P3w);
428 rotation_matrix.multiply(B);
430 //rotating to set P3 on the Oxy plane
431 TranslationRotationMatrix C;
433 rotation_matrix.multiply(C);