From ec6f4d3964bebc3ae825909313bd1294c61ae7c6 Mon Sep 17 00:00:00 2001
From: eap <eap@opencascade.com>
Date: Wed, 17 Dec 2008 13:32:23 +0000
Subject: [PATCH] MEDMEM Industrialization 2008 -    double
 intersectCells(ConnType icell_A, ConnType icell_B, int nb_NodesA, int
 nb_NodesB); +    double intersectCells(ConnType icell_A, ConnType icell_B,
 int nb_NodesA, int nb_NodesB, +                          double*
 baryCentre=0);

---
 src/INTERP_KERNEL/ConvexIntersector.hxx |  21 ++++-
 src/INTERP_KERNEL/ConvexIntersector.txx | 103 ++++++++++++++----------
 2 files changed, 81 insertions(+), 43 deletions(-)

diff --git a/src/INTERP_KERNEL/ConvexIntersector.hxx b/src/INTERP_KERNEL/ConvexIntersector.hxx
index e86166f4d..f8fa14332 100644
--- a/src/INTERP_KERNEL/ConvexIntersector.hxx
+++ b/src/INTERP_KERNEL/ConvexIntersector.hxx
@@ -1,3 +1,21 @@
+//  Copyright (C) 2007-2008  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 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
+//
+//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
 #ifndef __CONVEXINTERSECTOR_HXX__
 #define __CONVEXINTERSECTOR_HXX__
 
@@ -18,7 +36,8 @@ namespace INTERP_KERNEL
     ConvexIntersector(const MyMeshType& mesh_A, const MyMeshType& mesh_B, 
                       double dimCaracteristic, double precision, double medianPlane,
                       bool doRotate, int printLevel);
-    double intersectCells(ConnType icell_A, ConnType icell_B, int nb_NodesA, int nb_NodesB);
+    double intersectCells(ConnType icell_A, ConnType icell_B, int nb_NodesA, int nb_NodesB,
+                          double* baryCentre=0);
   private :
     double _epsilon;
     const ConnType *_connectA;
diff --git a/src/INTERP_KERNEL/ConvexIntersector.txx b/src/INTERP_KERNEL/ConvexIntersector.txx
index 53af77e48..3816ac2a6 100644
--- a/src/INTERP_KERNEL/ConvexIntersector.txx
+++ b/src/INTERP_KERNEL/ConvexIntersector.txx
@@ -1,3 +1,21 @@
+//  Copyright (C) 2007-2008  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 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
+//
+//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
 #ifndef __CONVEXINTERSECTOR_TXX__
 #define __CONVEXINTERSECTOR_TXX__
 
@@ -34,10 +52,11 @@ namespace INTERP_KERNEL
 
   template<class MyMeshType>
   double ConvexIntersector<MyMeshType>::intersectCells(ConnType icell_A, ConnType icell_B, 
-                                                       int nb_NodesA, int nb_NodesB)
+                                                       int nb_NodesA, int nb_NodesB,
+                                                       double* baryCentre)
   {
     double result = 0;
-		int orientation = 1;
+    int orientation = 1;
 
     /*** Obtain the coordinates of A and B ***/
     std::vector<double> Coords_A (SPACEDIM*nb_NodesA);
@@ -48,50 +67,50 @@ namespace INTERP_KERNEL
     const double * Pi_last= _coordsA + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectA[OTT<ConnType,numPol>::conn2C(_connIndexA[OTT<ConnType,numPol>::ind2C(icell_A)]+i_last)]);
 
     for (int i_A=0; i_A<nb_NodesA; i_A++)
+    {
+      const double * Pi_A = _coordsA + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectA[OTT<ConnType,numPol>::conn2C(_connIndexA[OTT<ConnType,numPol>::ind2C(icell_A)]+i_A)]);
+      if(distance2<SPACEDIM>(Pi_last, Pi_A)>_epsilon)
       {
-        const double * Pi_A = _coordsA + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectA[OTT<ConnType,numPol>::conn2C(_connIndexA[OTT<ConnType,numPol>::ind2C(icell_A)]+i_A)]);
-        if(distance2<SPACEDIM>(Pi_last, Pi_A)>_epsilon)
-          {
-            for (int idim=0; idim<SPACEDIM; idim++)
-              Coords_A[SPACEDIM*i_A+idim]=Pi_A[idim];
-            i_last=i_A; Pi_last = Pi_A;
-          }
-        else 
-          nb_dist_NodesA--;
+        for (int idim=0; idim<SPACEDIM; idim++)
+          Coords_A[SPACEDIM*i_A+idim]=Pi_A[idim];
+        i_last=i_A; Pi_last = Pi_A;
       }
+      else 
+        nb_dist_NodesA--;
+    }
     i_last = nb_NodesB - 1;
     Pi_last=_coordsB + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectB[OTT<ConnType,numPol>::conn2C(_connIndexB[OTT<ConnType,numPol>::ind2C(icell_B)]+i_last)]);
     for (int i_B=0; i_B<nb_NodesB; i_B++)
+    {
+      const double * Pi_B=_coordsB+SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectB[OTT<ConnType,numPol>::conn2C(_connIndexB[OTT<ConnType,numPol>::ind2C(icell_B)]+i_B)]);
+      if(distance2<SPACEDIM>(Pi_last, Pi_B)>_epsilon)
       {
-        const double * Pi_B=_coordsB+SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectB[OTT<ConnType,numPol>::conn2C(_connIndexB[OTT<ConnType,numPol>::ind2C(icell_B)]+i_B)]);
-        if(distance2<SPACEDIM>(Pi_last, Pi_B)>_epsilon)
-          {
-            for (int idim=0; idim<SPACEDIM; idim++)
-              Coords_B[SPACEDIM*i_B+idim]=Pi_B[idim];
-            i_last=i_B; Pi_last = Pi_B;
-          }
-        else
-          nb_dist_NodesB--;
+        for (int idim=0; idim<SPACEDIM; idim++)
+          Coords_B[SPACEDIM*i_B+idim]=Pi_B[idim];
+        i_last=i_B; Pi_last = Pi_B;
       }
-      
+      else
+        nb_dist_NodesB--;
+    }
+
     /*** project cells A and B on the median plane ***/
     /***  and rotate the median plane ***/
     if(SPACEDIM==3) 
-			orientation = projection(Coords_A, Coords_B, nb_dist_NodesA, nb_dist_NodesB,_epsilon,
+      orientation = projection(Coords_A, Coords_B, nb_dist_NodesA, nb_dist_NodesB,_epsilon,
                                PlanarIntersector<MyMeshType>::_medianPlane,
                                PlanarIntersector<MyMeshType>::_doRotate);
 
     //DEBUG PRINTS
     if(PlanarIntersector<MyMeshType>::_printLevel >= 3) 
-      {
-        std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl;
-        std::cout << std::endl << "icell_A= " << icell_A << ", nb nodes A= " <<  nb_NodesA << std::endl;
-        for(int i_A =0; i_A< nb_NodesA; i_A++)
-          {for (int idim=0; idim<SPACEDIM; idim++) std::cout << Coords_A[SPACEDIM*i_A+idim] << " "; std::cout << std::endl;}
-        std::cout << std::endl << "icell_B= " << icell_B << ", nb nodes B= " <<  nb_NodesB << std::endl;
-        for(int i_B =0; i_B< nb_NodesB; i_B++)
-          {for (int idim=0; idim<SPACEDIM; idim++) std::cout << Coords_B[SPACEDIM*i_B+idim]<< " "; std::cout << std::endl; }
-      }
+    {
+      std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl;
+      std::cout << std::endl << "icell_A= " << icell_A << ", nb nodes A= " <<  nb_NodesA << std::endl;
+      for(int i_A =0; i_A< nb_NodesA; i_A++)
+      {for (int idim=0; idim<SPACEDIM; idim++) std::cout << Coords_A[SPACEDIM*i_A+idim] << " "; std::cout << std::endl;}
+      std::cout << std::endl << "icell_B= " << icell_B << ", nb nodes B= " <<  nb_NodesB << std::endl;
+      for(int i_B =0; i_B< nb_NodesB; i_B++)
+      {for (int idim=0; idim<SPACEDIM; idim++) std::cout << Coords_B[SPACEDIM*i_B+idim]<< " "; std::cout << std::endl; }
+    }
 
     /*** Compute the intersection area ***/
     INTERP_KERNEL::PolygonAlgorithms<SPACEDIM> P(_epsilon, PlanarIntersector<MyMeshType>::_precision);
@@ -100,20 +119,20 @@ namespace INTERP_KERNEL
     double area[SPACEDIM];
     int nb_inter =((int)inter.size())/SPACEDIM;
     for(int i = 1; i<nb_inter-1; i++)
-      {
-        INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
-        result +=0.5*norm<SPACEDIM>(area);
-      }
+    {
+      INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
+      result +=0.5*norm<SPACEDIM>(area);
+    }
 
     //DEBUG prints
     if(PlanarIntersector<MyMeshType>::_printLevel >= 3)
-      {       
-        std::cout << std::endl << "Number of nodes of the intersection = "<<  nb_inter << std::endl;
-        for(int i=0; i<  nb_inter; i++)
-          {for (int idim=0; idim<SPACEDIM; idim++) std::cout << inter[SPACEDIM*i+idim]<< " "; std::cout << std::endl;}
-        std::cout << std::endl <<"Intersection area = " << result << std::endl;
-      }
-    
+    {       
+      std::cout << std::endl << "Number of nodes of the intersection = "<<  nb_inter << std::endl;
+      for(int i=0; i<  nb_inter; i++)
+      {for (int idim=0; idim<SPACEDIM; idim++) std::cout << inter[SPACEDIM*i+idim]<< " "; std::cout << std::endl;}
+      std::cout << std::endl <<"Intersection area = " << result << std::endl;
+    }
+
     return orientation*result;
   }
 }
-- 
2.39.2