-// Copyright (C) 2007-2008 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2019 CEA/DEN, EDF R&D
//
-// This library is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
//
-// This library is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// Lesser General Public License for more details.
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
//
-// You should have received a copy of the GNU Lesser General Public
-// License along with this library; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
-// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
+// Author : Anthony Geay (CEA/DEN)
#ifndef __PLANARINTERSECTOR_TXX__
#define __PLANARINTERSECTOR_TXX__
namespace INTERP_KERNEL
{
template<class MyMeshType, class MyMatrix>
- PlanarIntersector<MyMeshType,MyMatrix>::PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel):
+ PlanarIntersector<MyMeshType,MyMatrix>::PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double minDot3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel):
_meshT(meshT),_meshS(meshS),
- _dim_caracteristic(dimCaracteristic),_precision(precision),_median_plane(medianPlane),
+ _dim_caracteristic(dimCaracteristic),_max_distance_3Dsurf_intersect(md3DSurf),_min_dot_btw_3Dsurf_intersect(minDot3DSurf),_precision(precision),_median_plane(medianPlane),
_do_rotate(doRotate),_orientation(orientation),_print_level(printLevel)
{
_connectT=meshT.getConnectivityPtr();
int ibox=0;
for(long icell=0; icell<nbelems; icell++)
{
- int nb_nodes_per_elem =conn_index[icell+1]-conn_index[icell];
+ ConnType nb_nodes_per_elem =conn_index[icell+1]-conn_index[icell];
//initializing bounding box limits
for(int idim=0; idim<SPACEDIM; idim++)
{
bbox[2*SPACEDIM*ibox+2*idim+1] = -std::numeric_limits<double>::max();
}
//updating the bounding box with each node of the element
- for (int j=0; j<nb_nodes_per_elem; j++)
+ for (ConnType j=0; j<nb_nodes_per_elem; j++)
{
const double* coord_node=coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(conn[OTT<ConnType,numPol>::conn2C(conn_index[icell]+j)]);
for(int idim=0; idim<SPACEDIM; idim++)
}
/*!
- Computes the bouding box of a given element. iP in numPol mode.
+ Computes the bounding box of a given element. iP in numPol mode.
*/
template<class MyMeshType, class MyMatrix>
void PlanarIntersector<MyMeshType,MyMatrix>::getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes)
\param bbox vector containing the bounding boxes
*/
template<class MyMeshType, class MyMatrix>
- void PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double Surf3DAdjustmentEps)
+ void PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double surf3DAdjustmentEps, double surf3DAdjustmentEpsAbs)
{
/* We build the segment tree for locating possible matching intersections*/
- long size = bbox.size()/(2*SPACEDIM);
- for (int i=0; i<size; i++)
+ std::size_t size = bbox.size()/(2*SPACEDIM);
+ for (std::size_t i=0; i<size; i++)
{
double max=- std::numeric_limits<double>::max();
for(int idim=0; idim<SPACEDIM; idim++)
}
for(int idim=0; idim<SPACEDIM; idim++)
{
- bbox[i*2*SPACEDIM+2*idim ] -= Surf3DAdjustmentEps*max;
- bbox[i*2*SPACEDIM+2*idim+1] += Surf3DAdjustmentEps*max;
+ bbox[i*2*SPACEDIM+2*idim ] -= surf3DAdjustmentEps*max+surf3DAdjustmentEpsAbs;
+ bbox[i*2*SPACEDIM+2*idim+1] += surf3DAdjustmentEps*max+surf3DAdjustmentEpsAbs;
}
}
}
template<class MyMeshType, class MyMatrix>
void PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT)
{
- int nbNodesT=_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1]-_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
+ ConnType nbNodesT=_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1]-_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
coordsT.resize(SPACEDIM*nbNodesT);
for (ConnType iT=0; iT<nbNodesT; iT++)
for(int idim=0; idim<SPACEDIM; idim++)
template<class MyMeshType, class MyMatrix>
void PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS)
{
- int nbNodesS=_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1]-_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
+ ConnType nbNodesS=_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1]-_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
coordsS.resize(SPACEDIM*nbNodesS);
for (ConnType iS=0; iS<nbNodesS; iS++)
for(int idim=0; idim<SPACEDIM; idim++)
coordsS[SPACEDIM*iS+idim]=_coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
}
+
+ /*!
+ * @param icellT id in target mesh in format of MyMeshType.
+ * @param offset is a value in C format that indicates the number of circular permutation.
+ * @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
+ */
+ template<class MyMeshType, class MyMatrix>
+ void PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinatesPermute(ConnType icellT, ConnType offset, std::vector<double>& coordsT)
+ {
+ ConnType nbNodesT=_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1]-_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
+ coordsT.resize(SPACEDIM*nbNodesT);
+ for (ConnType iTTmp=0; iTTmp<nbNodesT; iTTmp++)
+ {
+ ConnType iT=(iTTmp+offset)%nbNodesT;
+ for(int idim=0; idim<SPACEDIM; idim++)
+ coordsT[SPACEDIM*iTTmp+idim]=_coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
+ }
+ }
+
+ /*!
+ * @param icellS id in source mesh in format of MyMeshType.
+ * @param offset is a value in C format that indicates the number of circular permutation.
+ * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
+ */
+ template<class MyMeshType, class MyMatrix>
+ void PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinatesPermute(ConnType icellS, ConnType offset, std::vector<double>& coordsS)
+ {
+ ConnType nbNodesS=_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1]-_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
+ coordsS.resize(SPACEDIM*nbNodesS);
+ for (ConnType iSTmp=0; iSTmp<nbNodesS; iSTmp++)
+ {
+ ConnType iS=(iSTmp+offset)%nbNodesS;
+ for(int idim=0; idim<SPACEDIM; idim++)
+ coordsS[SPACEDIM*iSTmp+idim]=_coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
+ }
+ }
/*!
* @param icellT id in target mesh in format of MyMeshType.
template<class MyMeshType, class MyMatrix>
void PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, std::vector<double>& coordsT, std::vector<double>& coordsS, int& orientation)
{
- /*double epsilon=_precision*_dim_caracteristic;
- coordsT.resize(SPACEDIM*nbNodesT);
- coordsS.resize(SPACEDIM*nbNodesS);
- int nb_dist_NodesT=nbNodesT;
- int nb_dist_NodesS=nbNodesS;
- int i_last = nbNodesT - 1;
- const double * Pi_last=_coordsT +_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+i_last)];
-
- for (int iT=0; iT<nbNodesT; iT++)
- {
- const double * PiT = _coordsT + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)]);
- if(distance2<SPACEDIM>(Pi_last, PiT)>epsilon)
- {
- for (int idim=0; idim<SPACEDIM; idim++)
- coordsT[SPACEDIM*iT+idim]=PiT[idim];
- i_last=iT; Pi_last = PiT;
- }
- else
- nb_dist_NodesT--;
- }
- coordsT.resize(nb_dist_NodesT*SPACEDIM);
- i_last = nbNodesS - 1;
- Pi_last=_coordsS + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+i_last)]);
- for (int iS=0; iS<nbNodesS; iS++)
- {
- const double * PiS=_coordsS+SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)]);
- if(distance2<SPACEDIM>(Pi_last, PiS)>epsilon)
- {
- for (int idim=0; idim<SPACEDIM; idim++)
- coordsS[SPACEDIM*iS+idim]=PiS[idim];
- i_last=iS; Pi_last = PiS;
- }
- else
- nb_dist_NodesS--;
- }
- coordsS.resize(nb_dist_NodesS*SPACEDIM);
- //project cells T and S on the median plane
- // and rotate the median plane
- if(SPACEDIM==3)
- orientation = projectionThis(&coordsT[0], &coordsS[0], nb_dist_NodesT, nb_dist_NodesS);
-
- //DEBUG PRINTS
- if(_print_level >= 3)
- {
- std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl;
- std::cout << std::endl << "icellT= " << icellT << ", nb nodes T= " << nbNodesT << std::endl;
- for(int iT =0; iT< nbNodesT; iT++)
- {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsT[SPACEDIM*iT+idim] << " "; std::cout << std::endl;}
- std::cout << std::endl << "icellS= " << icellS << ", nb nodes S= " << nbNodesS << std::endl;
- for(int iS =0; iS< nbNodesS; iS++)
- {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsS[SPACEDIM*iS+idim]<< " "; std::cout << std::endl; }
- }*/
coordsT.resize(SPACEDIM*nbNodesT);
coordsS.resize(SPACEDIM*nbNodesS);
for (int idim=0; idim<SPACEDIM; idim++)
{for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsS[SPACEDIM*iS+idim]<< " "; std::cout << std::endl; }
}
}
+
+ /*!
+ * Filtering out zero surfaces and badly oriented surfaces
+ * _orientation = -1,0,1,2
+ * -1 : the intersection is taken into account if target and cells have different orientation
+ * 0 : the intersection is always taken into account
+ * 1 : the intersection is taken into account if target and cells have the same orientation
+ * 2 : the absolute value of intersection is always taken into account
+ */
+ template<class MyMeshType, class MyMatrix>
+ double PlanarIntersector<MyMeshType,MyMatrix>::getValueRegardingOption(double val) const
+ {
+ if(_orientation==0)
+ return val;
+ if(_orientation==2)
+ return fabs(val);
+ if (( val > 0.0 && _orientation==1) || ( val < 0.0 && _orientation==-1 ))
+ return _orientation*val;
+ return 0.;
+ }
template<class MyMeshType, class MyMatrix>
- int PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB)
+ int PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(double *Coords_A, double *Coords_B, ConnType nb_NodesA, ConnType nb_NodesB)
{
- return projection(Coords_A,Coords_B,nb_NodesA,nb_NodesB,_dim_caracteristic*_precision,_median_plane,_do_rotate);
+ return Projection(Coords_A,Coords_B,nb_NodesA,nb_NodesB,_dim_caracteristic*_precision,_max_distance_3Dsurf_intersect,_min_dot_btw_3Dsurf_intersect,_median_plane,_do_rotate);
}
template<class MyMeshType, class MyMatrix>
- int PlanarIntersector<MyMeshType,MyMatrix>::projection(double *Coords_A, double *Coords_B,
- int nb_NodesA, int nb_NodesB, double epsilon, double median_plane, bool do_rotate)
+ int PlanarIntersector<MyMeshType,MyMatrix>::Projection(double *Coords_A, double *Coords_B,
+ ConnType nb_NodesA, ConnType nb_NodesB, double epsilon, double md3DSurf, double minDot3DSurf, double median_plane, bool do_rotate)
{
double normal_A[3]={0,0,0};
double normal_B[3]={0,0,0};
bool same_orientation;
//Find the normal to cells A and B
- int i_A1=1;
- while(i_A1<nb_NodesA && distance2<SPACEDIM>(Coords_A,&Coords_A[SPACEDIM*i_A1])< epsilon) i_A1++;
- int i_A2=i_A1+1;
+ ConnType i_A1(1);
+ while(i_A1<nb_NodesA && distance2<SPACEDIM>(Coords_A,&Coords_A[SPACEDIM*i_A1])< epsilon)
+ i_A1++;
+ ConnType i_A2(i_A1+1);
crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
- double normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
+ double normA(sqrt(dotprod<SPACEDIM>(normal_A,normal_A)));
while(i_A2<nb_NodesA && normA < epsilon)
{
crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
}
- int i_B1=1;
- while(i_B1<nb_NodesB && distance2<SPACEDIM>(Coords_B,Coords_B+SPACEDIM*i_B1)< epsilon) i_B1++;
- int i_B2=i_B1+1;
+ ConnType i_B1(1);
+ while(i_B1<nb_NodesB && distance2<SPACEDIM>(Coords_B,Coords_B+SPACEDIM*i_B1)< epsilon)
+ i_B1++;
+ ConnType i_B2(i_B1+1);
crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
- double normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
+ double normB(sqrt(dotprod<SPACEDIM>(normal_B,normal_B)));
while(i_B2<nb_NodesB && normB < epsilon)
{
crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
}
+ //fabien option
+ if(md3DSurf>0.)
+ {
+ double coords_GA[3];
+ for(int i=0;i<3;i++)
+ {
+ coords_GA[i]=0.;
+ for (int j=0;j<nb_NodesA;j++)
+ coords_GA[i]+=Coords_A[3*j+i];
+ coords_GA[i]/=(double)nb_NodesA;
+ }
+ double G1[3],G2[3],G3[3];
+ for(int i=0;i<3;i++)
+ {
+ G1[i]=Coords_B[i]-coords_GA[i];
+ G2[i]=Coords_B[i+3]-coords_GA[i];
+ G3[i]=Coords_B[i+6]-coords_GA[i];
+ }
+ double prodvect[3];
+ prodvect[0]=G1[1]*G2[2]-G1[2]*G2[1];
+ prodvect[1]=G1[2]*G2[0]-G1[0]*G2[2];
+ prodvect[2]=G1[0]*G2[1]-G1[1]*G2[0];
+ double prodscal=prodvect[0]*G3[0]+prodvect[1]*G3[1]+prodvect[2]*G3[2];
+ if(fabs(prodscal)>md3DSurf)
+ return 0;
+ }
if(i_A2<nb_NodesA && i_B2<nb_NodesB)
{
//Build the normal of the median plane
- same_orientation = dotprod<SPACEDIM>(normal_A,normal_B)>=0;
+
+ double dotProd(dotprod<SPACEDIM>(normal_A,normal_B)/(normA*normB));
+
+ if(fabs(dotProd)<minDot3DSurf)
+ return 0;
+
+ same_orientation=(dotProd>=0);
if(!same_orientation)
- for(int idim =0; idim< SPACEDIM; idim++) normal_A[idim] *=-1;
+ for(int idim =0; idim< SPACEDIM; idim++)
+ normal_A[idim] *=-1;
- double normB= sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
+ double normBB(sqrt(dotprod<SPACEDIM>(normal_B,normal_B)));
for(int idim =0; idim< SPACEDIM; idim++)
- linear_comb[idim] = median_plane*normal_A[idim]/normA + (1-median_plane)*normal_B[idim]/normB;
+ linear_comb[idim] = median_plane*normal_A[idim]/normA + (1-median_plane)*normal_B[idim]/normBB;
double norm= sqrt(dotprod<SPACEDIM>(linear_comb,linear_comb));
//Necessarily: norm>epsilon, no need to check
- for(int idim =0; idim< SPACEDIM; idim++) linear_comb[idim]/=norm;
+ for(int idim =0; idim< SPACEDIM; idim++)
+ linear_comb[idim]/=norm;
//Project the nodes of A and B on the median plane
- for(int i_A=0; i_A<nb_NodesA; i_A++)
+ for(ConnType i_A=0; i_A<nb_NodesA; i_A++)
{
proj = dotprod<SPACEDIM>(&Coords_A[SPACEDIM*i_A],linear_comb);
for(int idim =0; idim< SPACEDIM; idim++)
Coords_A[SPACEDIM*i_A+idim] -= proj*linear_comb[idim];
}
- for(int i_B=0; i_B<nb_NodesB; i_B++)
+ for(ConnType i_B=0; i_B<nb_NodesB; i_B++)
{
proj = dotprod<SPACEDIM>(Coords_B+SPACEDIM*i_B,linear_comb);
for(int idim =0; idim< SPACEDIM; idim++)
Coords_B[SPACEDIM*i_B+idim] -= proj*linear_comb[idim];
}
- //Buid the matrix sending A into the Oxy plane and apply it to A and B
+ //Build the matrix sending A into the Oxy plane and apply it to A and B
if(do_rotate)
{
TranslationRotationMatrix rotation;
//rotate3DTriangle(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2], rotation);
- rotate3DTriangle(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2, rotation);
- for (int i=0; i<nb_NodesA; i++) rotation.transform_vector(Coords_A+SPACEDIM*i);
- for (int i=0; i<nb_NodesB; i++) rotation.transform_vector(Coords_B+SPACEDIM*i);
+ Rotate3DTriangle(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2, rotation);
+ for (int i=0; i<nb_NodesA; i++)
+ rotation.transform_vector(Coords_A+SPACEDIM*i);
+ for (int i=0; i<nb_NodesB; i++)
+ rotation.transform_vector(Coords_B+SPACEDIM*i);
}
- if(same_orientation)
- return 1;
- else return -1;
+ return same_orientation?1:-1;
}
else
{
- std::cout << " Maille dégénérée " << "epsilon = " << epsilon << std::endl;
+ std::cout << " Degenerated cell " << "epsilon = " << epsilon << std::endl;
std::cout << " i_A1= " << i_A1 << " i_A2= " << i_A2 << std::endl;
std::cout << " distance2<SPACEDIM>(Coords_A,&Coords_A[i_A1])= " << distance2<SPACEDIM>(Coords_A,&Coords_A[i_A1]) << std::endl;
std::cout << "abs(normal_A) = " << fabs(normal_A[0]) << " ; " <<fabs( normal_A[1]) << " ; " << fabs(normal_A[2]) << std::endl;
std::cout << " i_B1= " << i_B1 << " i_B2= " << i_B2 << std::endl;
std::cout << " distance2<SPACEDIM>(&Coords_B[0],&Coords_B[i_B1])= " << distance2<SPACEDIM>(Coords_B,Coords_B+i_B1) << std::endl;
std::cout << "normal_B = " << normal_B[0] << " ; " << normal_B[1] << " ; " << normal_B[2] << std::endl;
-
return 1;
}
}
template<class MyMeshType, class MyMatrix>
- void PlanarIntersector<MyMeshType,MyMatrix>::rotate3DTriangle(double* PP1, double*PP2, double*PP3,
+ void PlanarIntersector<MyMeshType,MyMatrix>::Rotate3DTriangle(double* PP1, double*PP2, double*PP3,
TranslationRotationMatrix& rotation_matrix)
{
//initializes
