//
// 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
//
// 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
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
- PlanarIntersector<MyMeshType,MyMatrix>::PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, 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):
- _dim_caracteristic(dimCaracteristic),_max_distance_3Dsurf_intersect(md3DSurf),_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();
_do_rotate(doRotate),_orientation(orientation),_print_level(printLevel)
{
_connectT=meshT.getConnectivityPtr();
bbox[2*SPACEDIM*ibox+2*idim+1] = -std::numeric_limits<double>::max();
}
//updating the bounding box with each node of the element
bbox[2*SPACEDIM*ibox+2*idim+1] = -std::numeric_limits<double>::max();
}
//updating the bounding box with each node of the element
{
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++)
{
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++)
*/
template<class MyMeshType, class MyMatrix>
void PlanarIntersector<MyMeshType,MyMatrix>::getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes)
*/
template<class MyMeshType, class MyMatrix>
void PlanarIntersector<MyMeshType,MyMatrix>::getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes)
template<class MyMeshType, class MyMatrix>
void PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT)
{
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++)
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)
{
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.resize(SPACEDIM*nbNodesS);
for (ConnType iS=0; iS<nbNodesS; iS++)
for(int idim=0; idim<SPACEDIM; idim++)
* @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
*/
template<class MyMeshType, class MyMatrix>
* @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, int offset, std::vector<double>& coordsT)
+ void PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinatesPermute(ConnType icellT, ConnType offset, 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)];
* @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
*/
template<class MyMeshType, class MyMatrix>
* @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, int offset, std::vector<double>& coordsS)
+ void PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinatesPermute(ConnType icellS, ConnType offset, 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)];
- 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,_max_distance_3Dsurf_intersect,_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);
- int PlanarIntersector<MyMeshType,MyMatrix>::projection(double *Coords_A, double *Coords_B,
- int nb_NodesA, int nb_NodesB, double epsilon, double md3DSurf, 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};
{
double normal_A[3]={0,0,0};
double normal_B[3]={0,0,0};
crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],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);
while(i_A2<nb_NodesA && normA < epsilon)
{
crossprod<SPACEDIM>(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,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);
while(i_B2<nb_NodesB && normB < epsilon)
{
crossprod<SPACEDIM>(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2,normal_B);
- 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;
+ 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;
- 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);
for(int idim =0; idim< SPACEDIM; idim++)
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] = 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
{
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];
}
{
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];
}
{
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];
}
{
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];
}
if(do_rotate)
{
TranslationRotationMatrix rotation;
//rotate3DTriangle(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2], rotation);
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);
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;
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;
- void PlanarIntersector<MyMeshType,MyMatrix>::rotate3DTriangle(double* PP1, double*PP2, double*PP3,
+ void PlanarIntersector<MyMeshType,MyMatrix>::Rotate3DTriangle(double* PP1, double*PP2, double*PP3,