Added orientation parameter for filtering

diff --git a/src/INTERP_KERNEL/PlanarIntersector.hxx b/src/INTERP_KERNEL/PlanarIntersector.hxx
index 0aa2aa8bfa86c6a58c6a4d20204045dd9020ba47..eb2d371b1af3b51e5f7c654f16bd083f9faafd26 100644 (file)
--- a/src/INTERP_KERNEL/PlanarIntersector.hxx
+++ b/src/INTERP_KERNEL/PlanarIntersector.hxx
@@ -35,7 +35,7 @@ namespace INTERP_KERNEL
     void adjustBoundingBoxes(std::vector<double>& bbox, double Surf3DAdjustmentEps);
     inline void getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes);
   protected :
-    static void projection(std::vector<double>& Coords_A, std::vector<double>& Coords_B, 
+    static int projection(std::vector<double>& Coords_A, std::vector<double>& Coords_B, 
                            int nb_NodesA, int nb_NodesB, double epsilon, double median_plane, bool do_rotate);
     static void rotate3DTriangle( double* PP1, double*PP2, double*PP3,
                                   TranslationRotationMatrix& rotation_matrix);

diff --git a/src/INTERP_KERNEL/PlanarIntersector.txx b/src/INTERP_KERNEL/PlanarIntersector.txx
index f37facc5046167de2a8f192678c8e68a266e6346..280326b9acc18dfc891ba0b68d8310b5e24c00ca 100644 (file)
--- a/src/INTERP_KERNEL/PlanarIntersector.txx
+++ b/src/INTERP_KERNEL/PlanarIntersector.txx
@@ -119,69 +119,83 @@ namespace INTERP_KERNEL
   }
   
   template<class MyMeshType>
-  void PlanarIntersector<MyMeshType>::projection(std::vector< double>& Coords_A, std::vector< double>& Coords_B, 
+  int PlanarIntersector<MyMeshType>::projection(std::vector< double>& Coords_A, std::vector< double>& Coords_B, 
                                                  int nb_NodesA, int nb_NodesB, double epsilon, double median_plane, bool do_rotate)
   {
     double normal_A[3]={0,0,0};
     double normal_B[3]={0,0,0};
     double linear_comb[3];
     double proj;
+               bool same_orientation;
+
     //Find the normal to cells A and B
     int i_A1=1;
     while(i_A1<nb_NodesA && distance2<SPACEDIM>(&Coords_A[0],&Coords_A[SPACEDIM*i_A1])< epsilon) i_A1++;
     int i_A2=i_A1+1;
     crossprod<SPACEDIM>(&Coords_A[0], &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
-    while(i_A2<nb_NodesA && dotprod<SPACEDIM>(normal_A,normal_A)<epsilon)
+               double normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
+    while(i_A2<nb_NodesA && normA < epsilon)
       {
         crossprod<SPACEDIM>(&Coords_A[0], &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2],normal_A);
         i_A2++;
+        normA = sqrt(dotprod<SPACEDIM>(normal_A,normal_A));
+
       }
     int i_B1=1;
     while(i_B1<nb_NodesB && distance2<SPACEDIM>(&Coords_B[0],&Coords_B[SPACEDIM*i_B1])< epsilon) i_B1++;
     int i_B2=i_B1+1;
     crossprod<SPACEDIM>(&Coords_B[0], &Coords_B[SPACEDIM*i_B1], &Coords_B[SPACEDIM*i_B2],normal_B);
-    while(i_B2<nb_NodesB && dotprod<SPACEDIM>(normal_B,normal_B)< epsilon)
+               double normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
+    while(i_B2<nb_NodesB && normB < epsilon)
       {
         crossprod<SPACEDIM>(&Coords_B[0], &Coords_B[SPACEDIM*i_B1], &Coords_B[SPACEDIM*i_B2],normal_B);
         i_B2++;
+                               normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
       }
+
     if(i_A2<nb_NodesA && i_B2<nb_NodesB)
       {
-        //Build the normal of the median plane
-        if(dotprod<SPACEDIM>(normal_A,normal_B)<0)
+                               //Build the normal of the median plane
+                               same_orientation = dotprod<SPACEDIM>(normal_A,normal_B)>=0;
+        
+                               if(!same_orientation)
           for(int idim =0; idim< SPACEDIM; idim++) normal_A[idim] *=-1;
-        for(int idim =0; idim< SPACEDIM; idim++)
-          linear_comb[idim] = median_plane*normal_A[idim] + (1-median_plane)*normal_B[idim];
+                               
+        double normB= 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;
         double norm= sqrt(dotprod<SPACEDIM>(linear_comb,linear_comb));
 
-        if(norm>epsilon)
-          {
-            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++)
-              {
-                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++)
-              {
-                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  
-            if(do_rotate)
-              {
-                TranslationRotationMatrix rotation;
-                //rotate3DTriangle(&Coords_A[0], &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2], rotation);
-                rotate3DTriangle(&Coords_B[0], &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]);
-              }
-          }
+                               //Necessarily: norm>epsilon, no need to check
+                               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++)
+                                       {
+                                               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++)
+                                       {
+                                               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  
+                               if(do_rotate)
+                                       {
+                                               TranslationRotationMatrix rotation;
+                                               //rotate3DTriangle(&Coords_A[0], &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2], rotation);
+                                               rotate3DTriangle(&Coords_B[0], &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;
       }
     else
       {
@@ -192,6 +206,8 @@ namespace INTERP_KERNEL
         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[0],&Coords_B[i_B1]) << std::endl;
         std::cout << "normal_B = " << normal_B[0] << " ; " << normal_B[1] << " ; " << normal_B[2] << std::endl;
+
+                               return 1;
       }
   }