From: eap <eap@opencascade.com>
Date: Wed, 17 Dec 2008 15:56:37 +0000 (+0000)
Subject: MEDMEM Industrialization 2008
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=aaf09124917f5e54a80b9a19c8a9345b5233c9d0;p=tools%2Fmedcoupling.git

MEDMEM Industrialization 2008
    calculate barycentre of intersection
---

diff --git a/src/INTERP_KERNEL/IntersectorTetra.hxx b/src/INTERP_KERNEL/IntersectorTetra.hxx
index a419c43e6..5b4dc72c7 100644
--- a/src/INTERP_KERNEL/IntersectorTetra.hxx
+++ b/src/INTERP_KERNEL/IntersectorTetra.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 __INTERSECTORTETRA_HXX__
 #define __INTERSECTORTETRA_HXX__
 
@@ -175,8 +193,8 @@ namespace INTERP_KERNEL
 
   /** 
    * \brief Class calculating the volume of intersection between a tetrahedral target element and
-   * source elements with triangular or quadratilateral faces.
-   *
+   * source elements with triangular or quadratilateral faces. Optionally it computes barycentre
+   * of intersection
    */
   template<class MyMeshType>
   class INTERPKERNEL_EXPORT IntersectorTetra : public TargetIntersector<typename MyMeshType::MyConnType>
@@ -191,7 +209,7 @@ namespace INTERP_KERNEL
     
     ~IntersectorTetra();
 
-    double intersectSourceCell(typename MyMeshType::MyConnType srcCell);
+    double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0);
 
   private:
     
@@ -224,7 +242,7 @@ namespace INTERP_KERNEL
 
     /// reference to the source mesh
     const MyMeshType& _srcMesh;
-                
+
   };
 
   /**
diff --git a/src/INTERP_KERNEL/IntersectorTetra.txx b/src/INTERP_KERNEL/IntersectorTetra.txx
index b501cec45..4ea8daf51 100644
--- a/src/INTERP_KERNEL/IntersectorTetra.txx
+++ b/src/INTERP_KERNEL/IntersectorTetra.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 __INTERSECTORTETRA_TXX__
 #define __INTERSECTORTETRA_TXX__
 
@@ -7,8 +25,9 @@
 #include "TransformedTriangle.hxx"
 #include "MeshUtils.hxx"
 #include "VectorUtils.hxx"
-#include "CellModel.hxx"
+#include "UniversalCellModel.hxx"
 #include "Log.hxx"
+#include "UnitTetraIntersectionBary.hxx"
 
 #include <cmath>
 #include <cassert>
@@ -83,11 +102,12 @@ namespace INTERP_KERNEL
    * with triangulated faces to avoid having to recalculate these.
    *
    * @param element      global number of the source element (1 <= srcCell < # source cells)
+   * @param element      bary centre of cells intersection
    */
   template<class MyMeshType>
-  double IntersectorTetra<MyMeshType>::intersectSourceCell(typename MyMeshType::MyConnType element)
+  double IntersectorTetra<MyMeshType>::intersectSourceCell(typename MyMeshType::MyConnType element,
+                                                           double* baryCentre)
   {
-        
     //{ could be done on outside?
     // check if we have planar tetra element
     if(_t->determinant() == 0.0)
@@ -98,8 +118,9 @@ namespace INTERP_KERNEL
       }
 
     // get type of cell
-    NormalizedCellType normCellType=_srcMesh.getTypeOfElement(element);
-    const CellModel& cellModelCell=CellModel::getCellModel(normCellType);
+//     NormalizedCellType normCellType=_srcMesh.getTypeOfElement(element);
+//     const CellModel& cellModelCell=CellModel::getCellModel(normCellType);
+    UniversalCellModel<MyMeshType> cellModelCell( _srcMesh, element );
     unsigned nbOfNodes4Type=cellModelCell.getNumberOfNodes();
     // halfspace filtering
     bool isOutside[8] = {true, true, true, true, true, true, true, true};
@@ -134,6 +155,9 @@ namespace INTERP_KERNEL
     
     if(!isTargetOutside)
       {
+        /// calculator of intersection barycentre
+        UnitTetraIntersectionBary baryCalculator;
+        
         for(unsigned ii = 0 ; ii < cellModelCell.getNumberOfSons() ; ++ii)
           {
             NormalizedCellType faceType = cellModelCell.getSonType(ii);
@@ -154,6 +178,8 @@ namespace INTERP_KERNEL
                       TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
                       calculateVolume(tri, key);
                       totalVolume += _volumes[key];
+                      if ( baryCentre )
+                        baryCalculator.addSide( tri );
                     } else {    
                     // count negative as face has reversed orientation
                     totalVolume -= _volumes[key];
@@ -212,6 +238,11 @@ namespace INTERP_KERNEL
               }
             delete [] faceNodes;
           }
+
+        if ( baryCentre ) {
+          baryCalculator.getBary( baryCentre );
+          _t->reverseApply( baryCentre, baryCentre );
+        }
       }
     delete [] cellNodes;
     // reset if it is very small to keep the matrix sparse
diff --git a/src/INTERP_KERNEL/TriangulationIntersector.hxx b/src/INTERP_KERNEL/TriangulationIntersector.hxx
index 04a5ed65c..56704a045 100644
--- a/src/INTERP_KERNEL/TriangulationIntersector.hxx
+++ b/src/INTERP_KERNEL/TriangulationIntersector.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 _TRIANGULATION_INTERSECTOR_HXX_
 #define _TRIANGULATION_INTERSECTOR_HXX_
 
@@ -16,7 +34,8 @@ namespace INTERP_KERNEL
   public:
     TriangulationIntersector(const MyMeshType& mesh_A, const MyMeshType& mesh_B,
                              double dimCaracteristic, double precision, double medianPlane, 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 :
     const ConnType *_connectA;
     const ConnType *_connectB;
diff --git a/src/INTERP_KERNEL/TriangulationIntersector.txx b/src/INTERP_KERNEL/TriangulationIntersector.txx
index ac4114189..9a2ae9653 100644
--- a/src/INTERP_KERNEL/TriangulationIntersector.txx
+++ b/src/INTERP_KERNEL/TriangulationIntersector.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 __TRIANGULATIONINTERSECTOR_TXX__
 #define __TRIANGULATIONINTERSECTOR_TXX__
 
@@ -31,72 +49,78 @@ namespace INTERP_KERNEL
 
   template<class MyMeshType>
   double TriangulationIntersector<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);
     std::vector<double> Coords_B (SPACEDIM*nb_NodesB);
     for (int idim=0; idim<SPACEDIM; idim++)
-      {
-        for (ConnType i_A=0; i_A<nb_NodesA; i_A++)
-          Coords_A[SPACEDIM*i_A+idim] = _coordsA[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectA[OTT<ConnType,numPol>::conn2C(_connIndexA[OTT<ConnType,numPol>::ind2C(icell_A)]+i_A)])+idim];
-        for (ConnType i_B=0; i_B<nb_NodesB; i_B++)
-          Coords_B[SPACEDIM*i_B+idim] = _coordsB[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectB[OTT<ConnType,numPol>::conn2C(_connIndexB[OTT<ConnType,numPol>::ind2C(icell_B)]+i_B)])+idim];
-      }
-    
+    {
+      for (ConnType i_A=0; i_A<nb_NodesA; i_A++)
+        Coords_A[SPACEDIM*i_A+idim] = _coordsA[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectA[OTT<ConnType,numPol>::conn2C(_connIndexA[OTT<ConnType,numPol>::ind2C(icell_A)]+i_A)])+idim];
+      for (ConnType i_B=0; i_B<nb_NodesB; i_B++)
+        Coords_B[SPACEDIM*i_B+idim] = _coordsB[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectB[OTT<ConnType,numPol>::conn2C(_connIndexB[OTT<ConnType,numPol>::ind2C(icell_B)]+i_B)])+idim];
+    }
+
     //project cells A and B on the median plane and rotate the median plane
     if(SPACEDIM==3)
       orientation = projection(Coords_A, Coords_B, nb_NodesA, nb_NodesB,
-															 PlanarIntersector<MyMeshType>::_dimCaracteristic * PlanarIntersector<MyMeshType>::_precision,
-															 PlanarIntersector<MyMeshType>::_medianPlane, PlanarIntersector<MyMeshType>::_doRotate);
-    
+                               PlanarIntersector<MyMeshType>::_dimCaracteristic * PlanarIntersector<MyMeshType>::_precision,
+                               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
     double area[SPACEDIM];
     for(ConnType i_A = 1; i_A<nb_NodesA-1; i_A++)
+    {
+      for(ConnType i_B = 1; i_B<nb_NodesB-1; i_B++)
       {
-        for(ConnType i_B = 1; i_B<nb_NodesB-1; i_B++)
-          {
-            std::vector<double> inter;
-            INTERP_KERNEL::intersec_de_triangle(&Coords_A[0],&Coords_A[SPACEDIM*i_A],&Coords_A[SPACEDIM*(i_A+1)],
-                                                &Coords_B[0],&Coords_B[SPACEDIM*i_B],&Coords_B[SPACEDIM*(i_B+1)],
-                                                inter, PlanarIntersector<MyMeshType>::_dimCaracteristic,
-                                                PlanarIntersector<MyMeshType>::_precision);
-            ConnType nb_inter=((ConnType)inter.size())/2;
-            if(nb_inter >3) inter=reconstruct_polygon(inter);
-            for(ConnType i = 1; i<nb_inter-1; i++)
-              {
-                INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
-                result +=0.5*fabs(area[0]);
-              }
-            //DEBUG prints
-            if(PlanarIntersector<MyMeshType>::_printLevel >= 3)
-              {
-                std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl;
-                for(ConnType i=0; i< nb_inter; i++)
-                  {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; }
-              }
-          }
+        std::vector<double> inter;
+        INTERP_KERNEL::intersec_de_triangle(&Coords_A[0],&Coords_A[SPACEDIM*i_A],&Coords_A[SPACEDIM*(i_A+1)],
+                                            &Coords_B[0],&Coords_B[SPACEDIM*i_B],&Coords_B[SPACEDIM*(i_B+1)],
+                                            inter, PlanarIntersector<MyMeshType>::_dimCaracteristic,
+                                            PlanarIntersector<MyMeshType>::_precision);
+        ConnType nb_inter=((ConnType)inter.size())/2;
+        if(nb_inter >3) inter=reconstruct_polygon(inter);
+        for(ConnType i = 1; i<nb_inter-1; i++)
+        {
+          INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
+          result +=0.5*fabs(area[0]);
+        }
+        if ( baryCentre ) {
+          std::vector<double> Bary=INTERP_KERNEL::bary_poly(inter);
+          baryCentre[0] = Bary[0];
+          baryCentre[1] = Bary[1];
+        }
+        //DEBUG prints
+        if(PlanarIntersector<MyMeshType>::_printLevel >= 3)
+        {
+          std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl;
+          for(ConnType i=0; i< nb_inter; i++)
+          {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; }
+        }
       }
-    
+    }
+
     //DEBUG PRINTS    
     if(PlanarIntersector<MyMeshType>::_printLevel >= 3) 
       std::cout << std::endl <<"Intersection area = " << result << std::endl;
-    
+
     return orientation*result;
   }
 }