From: ndjinga <ndjinga>
Date: Fri, 26 Sep 2008 14:46:05 +0000 (+0000)
Subject: Added orientation parameter for filtering
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=3add054f31f85c27020af3d2f5a135685b9b9890;p=tools%2Fmedcoupling.git

Added orientation parameter for filtering
---

diff --git a/src/INTERP_KERNEL/PlanarIntersector.hxx b/src/INTERP_KERNEL/PlanarIntersector.hxx
index 0aa2aa8bf..eb2d371b1 100644
--- 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 f37facc50..280326b9a 100644
--- 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;
       }
   }