Salome HOME
Indices are stored as mcIdType type instead of int to support switch to 64bits indexing
[tools/medcoupling.git] / src / INTERP_KERNEL / VolSurfUser.txx
index 31988e3e0789833975299db2fd75be962261693d..6e15126030f05eab24de16f8e477892ae6790309 100644 (file)
@@ -1,9 +1,9 @@
-// Copyright (C) 2007-2012  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2019  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.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
@@ -16,6 +16,7 @@
 //
 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
+// Author : Anthony Geay (CEA/DEN)
 #ifndef __VOLSURFUSER_TXX__
 #define __VOLSURFUSER_TXX__
 
 namespace INTERP_KERNEL
 {
   template<class ConnType, NumberingPolicy numPol, int SPACEDIM>
-  double computeVolSurfOfCell(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords)
+  double computeVolSurfOfCell(NormalizedCellType type, const ConnType *connec, mcIdType lgth, const double *coords)
   {
     switch(type)
       {
       case INTERP_KERNEL::NORM_SEG2 :
-      case INTERP_KERNEL::NORM_SEG3 :
       case INTERP_KERNEL::NORM_SEG4 :
         {
-          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
-          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
           return INTERP_KERNEL::calculateLgthForSeg2(coords+(SPACEDIM*N1),coords+(SPACEDIM*N2),SPACEDIM);
         }
+      case INTERP_KERNEL::NORM_SEG3 :
+        {
+          ConnType beginNode = OTT<ConnType,numPol>::coo2C(connec[0]);
+          ConnType endNode = OTT<ConnType,numPol>::coo2C(connec[1]);
+          ConnType middleNode = OTT<ConnType,numPol>::coo2C(connec[2]);
+          return INTERP_KERNEL::calculateLgthForSeg3(coords+(SPACEDIM*beginNode),coords+(SPACEDIM*endNode),coords+(SPACEDIM*middleNode),SPACEDIM);
+        }
       case INTERP_KERNEL::NORM_TRI3 :
-      case INTERP_KERNEL::NORM_TRI6 :
-      case INTERP_KERNEL::NORM_TRI7 :
         {
-          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
-          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
-          int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
+          ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          ConnType N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
               
           return INTERP_KERNEL::calculateAreaForTria(coords+(SPACEDIM*N1),
                                                      coords+(SPACEDIM*N2),
@@ -55,14 +60,25 @@ namespace INTERP_KERNEL
         }
         break;
             
+      case INTERP_KERNEL::NORM_TRI6 :
+      case INTERP_KERNEL::NORM_TRI7 :
+        {
+          const double *pts[6];
+          pts[0] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]);
+          pts[1] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]);
+          pts[2] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]);
+          pts[3] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[3]);
+          pts[4] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[4]);
+          pts[5] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[5]);
+          return INTERP_KERNEL::calculateAreaForQPolyg(pts,6,SPACEDIM);
+        }
+        break;
       case INTERP_KERNEL::NORM_QUAD4 :
-      case INTERP_KERNEL::NORM_QUAD8 :
-      case INTERP_KERNEL::NORM_QUAD9 :
         {
-          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
-          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
-          int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
-          int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
+          ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          ConnType N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
+          ConnType N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
               
           return INTERP_KERNEL::calculateAreaForQuad(coords+SPACEDIM*N1,
                                                      coords+SPACEDIM*N2,
@@ -71,7 +87,21 @@ namespace INTERP_KERNEL
                                                      SPACEDIM);
         }
         break;
-            
+      case INTERP_KERNEL::NORM_QUAD8 :
+      case INTERP_KERNEL::NORM_QUAD9 :  
+        {
+          const double *pts[8];
+          pts[0] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]);
+          pts[1] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]);
+          pts[2] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]);
+          pts[3] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[3]);
+          pts[4] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[4]);
+          pts[5] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[5]);
+          pts[6] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[6]);
+          pts[7] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[7]);
+          return INTERP_KERNEL::calculateAreaForQPolyg(pts,8,SPACEDIM);
+        }
+        break;
       case INTERP_KERNEL::NORM_POLYGON :
         {          
           const double **pts=new const double *[lgth];
@@ -82,13 +112,23 @@ namespace INTERP_KERNEL
           return val;
         }
         break;
+      case INTERP_KERNEL::NORM_QPOLYG :
+        {
+          const double **pts=new const double *[lgth];
+          for(int inod=0;inod<lgth;inod++)
+            pts[inod] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[inod]);
+          double val=INTERP_KERNEL::calculateAreaForQPolyg(pts,lgth,SPACEDIM);
+          delete [] pts;
+          return val;
+        }
+        break;
       case INTERP_KERNEL::NORM_TETRA4 :
       case INTERP_KERNEL::NORM_TETRA10 :
         {
-          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
-          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
-          int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
-          int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
+          ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          ConnType N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
+          ConnType N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
               
           return INTERP_KERNEL::calculateVolumeForTetra(coords+SPACEDIM*N1,
                                                         coords+SPACEDIM*N2,
@@ -100,11 +140,11 @@ namespace INTERP_KERNEL
       case INTERP_KERNEL::NORM_PYRA5 :
       case INTERP_KERNEL::NORM_PYRA13 :
         {
-          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
-          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
-          int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
-          int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
-          int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
+          ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          ConnType N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
+          ConnType N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
+          ConnType N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
               
           return INTERP_KERNEL::calculateVolumeForPyra(coords+SPACEDIM*N1,
                                                        coords+SPACEDIM*N2,
@@ -116,13 +156,14 @@ namespace INTERP_KERNEL
             
       case INTERP_KERNEL::NORM_PENTA6 :
       case INTERP_KERNEL::NORM_PENTA15 :
+      case INTERP_KERNEL::NORM_PENTA18 :
         {
-          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
-          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
-          int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
-          int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
-          int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
-          int N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
+          ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          ConnType N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
+          ConnType N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
+          ConnType N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
+          ConnType N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
               
           return INTERP_KERNEL::calculateVolumeForPenta(coords+SPACEDIM*N1,
                                                         coords+SPACEDIM*N2,
@@ -137,14 +178,14 @@ namespace INTERP_KERNEL
       case INTERP_KERNEL::NORM_HEXA20 :
       case INTERP_KERNEL::NORM_HEXA27 :
         {
-          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
-          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
-          int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
-          int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
-          int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
-          int N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
-          int N7 = OTT<ConnType,numPol>::coo2C(connec[6]);
-          int N8 = OTT<ConnType,numPol>::coo2C(connec[7]);
+          ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          ConnType N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
+          ConnType N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
+          ConnType N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
+          ConnType N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
+          ConnType N7 = OTT<ConnType,numPol>::coo2C(connec[6]);
+          ConnType N8 = OTT<ConnType,numPol>::coo2C(connec[7]);
               
           return INTERP_KERNEL::calculateVolumeForHexa(coords+SPACEDIM*N1,
                                                        coords+SPACEDIM*N2,
@@ -158,7 +199,7 @@ namespace INTERP_KERNEL
         break;
       case INTERP_KERNEL::NORM_HEXGP12:
         {
-          const int connecHexa12[43]={
+          const ConnType connecHexa12[43]={
             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[4]),OTT<ConnType,numPol>::coo2C(connec[5]),-1,
             OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[10]),OTT<ConnType,numPol>::coo2C(connec[9]),OTT<ConnType,numPol>::coo2C(connec[8]),OTT<ConnType,numPol>::coo2C(connec[7]),-1,
             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[7]),OTT<ConnType,numPol>::coo2C(connec[1]),-1,
@@ -180,7 +221,7 @@ namespace INTERP_KERNEL
   }
 
   template<class ConnType, NumberingPolicy numPolConn>
-  double computeVolSurfOfCell2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim)
+  double computeVolSurfOfCell2(NormalizedCellType type, const ConnType *connec, mcIdType lgth, const double *coords, int spaceDim)
   {
     if(spaceDim==3)
       return computeVolSurfOfCell<ConnType,numPolConn,3>(type,connec,lgth,coords);
@@ -193,12 +234,11 @@ namespace INTERP_KERNEL
 
 
   template<class ConnType, NumberingPolicy numPol,int SPACEDIM>
-  void computeBarycenter(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, double *res)
+  void computeBarycenter(NormalizedCellType type, const ConnType *connec, mcIdType lgth, const double *coords, double *res)
   {
     switch(type)
       {
       case NORM_SEG2:
-      case NORM_SEG3:
       case NORM_SEG4:
         {
           std::copy(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]),
@@ -207,8 +247,30 @@ namespace INTERP_KERNEL
           std::transform(res,res+SPACEDIM,res,std::bind2nd(std::multiplies<double>(),0.5));
           break;
         }
+      case NORM_SEG3:
+        {
+          if(SPACEDIM==2)
+            {
+              Edge *ed=Edge::BuildEdgeFrom3Points(coords+2*OTT<ConnType,numPol>::coo2C(connec[0]),coords+2*OTT<ConnType,numPol>::coo2C(connec[2]),coords+2*OTT<ConnType,numPol>::coo2C(connec[1]));
+              ed->getBarycenter(res);
+              ed->decrRef();
+            }
+          else if(SPACEDIM==1)
+            {
+              *res=(coords[OTT<ConnType,numPol>::coo2C(connec[0])]+coords[OTT<ConnType,numPol>::coo2C(connec[1])])/2.;
+            }
+          else if(SPACEDIM==3)
+            {
+              std::copy(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]),
+                        coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]+1),res);
+              std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]),res,std::plus<double>());
+              std::transform(res,res+SPACEDIM,res,std::bind2nd(std::multiplies<double>(),0.5));
+            }
+          else
+            throw INTERP_KERNEL::Exception("computeBarycenter for SEG3 only SPACEDIM 1,2 or 3 supported !");
+          break;
+        }
       case NORM_TRI3:
-      case NORM_TRI6:
       case NORM_TRI7:
         {
           std::copy(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]),
@@ -218,6 +280,25 @@ namespace INTERP_KERNEL
           std::transform(res,res+SPACEDIM,res,std::bind2nd(std::multiplies<double>(),1./3.));
           break;
         }
+      case NORM_TRI6:
+        {
+          if(SPACEDIM==2)
+            {
+              double *pts[6];
+              pts[0] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]));
+              pts[1] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]));
+              pts[2] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]));
+              pts[3] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[3]));
+              pts[4] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[4]));
+              pts[5] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[5]));
+              computeQPolygonBarycenter2D(pts,6,2,res);
+            }
+          else if(SPACEDIM==3)
+            computePolygonBarycenter3D<ConnType,numPol>(connec,lgth/2,coords,res);
+          else
+            throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
+          break;
+        }
       case NORM_QUAD4:
       case NORM_POLYGON:
         {
@@ -232,13 +313,41 @@ namespace INTERP_KERNEL
       case NORM_QUAD8:
         {
           if(SPACEDIM==2)
-            computePolygonBarycenter2D<ConnType,numPol>(connec,lgth/2,coords,res);
+            {
+              double *pts[8];
+              pts[0] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]));
+              pts[1] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]));
+              pts[2] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]));
+              pts[3] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[3]));
+              pts[4] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[4]));
+              pts[5] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[5]));
+              pts[6] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[6]));
+              pts[7] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[7]));
+              computeQPolygonBarycenter2D(pts,8,2,res);
+            }
           else if(SPACEDIM==3)
             computePolygonBarycenter3D<ConnType,numPol>(connec,lgth/2,coords,res);
           else
             throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
           break;
         }
+      case INTERP_KERNEL::NORM_QPOLYG :
+        {
+          if(SPACEDIM==2)
+            {
+              double **pts=new double *[lgth];
+              for(int i=0;i<lgth;i++)
+                pts[i]=const_cast<double *>(coords+2*OTT<ConnType,numPol>::coo2C(connec[i]));
+              computeQPolygonBarycenter2D(pts,lgth,2,res);
+              delete [] pts;
+            }
+          else if(SPACEDIM==3)
+            computePolygonBarycenter3D<ConnType,numPol>(connec,lgth/2,coords,res);
+          else
+            throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
+          break;
+        }
+        break;
       case NORM_TETRA4:
         {
           res[0]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])]; 
@@ -267,7 +376,7 @@ namespace INTERP_KERNEL
         }
       case NORM_HEXA8:
         {
-          const int conn[29]={
+          const ConnType conn[29]={
             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[3]),-1,
             OTT<ConnType,numPol>::coo2C(connec[4]),OTT<ConnType,numPol>::coo2C(connec[7]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[5]),-1,
             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[7]),OTT<ConnType,numPol>::coo2C(connec[4]),-1,
@@ -280,7 +389,7 @@ namespace INTERP_KERNEL
         }
       case NORM_PENTA6:
         {
-          const int conn[22]={
+          const ConnType conn[22]={
             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[2]),-1,
             OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[4]),-1,
             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[3]),-1,
@@ -292,7 +401,7 @@ namespace INTERP_KERNEL
         }
       case INTERP_KERNEL::NORM_HEXGP12:
         {
-          const int connecHexa12[43]={
+          const ConnType connecHexa12[43]={
             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[4]),OTT<ConnType,numPol>::coo2C(connec[5]),-1,
             OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[10]),OTT<ConnType,numPol>::coo2C(connec[9]),OTT<ConnType,numPol>::coo2C(connec[8]),OTT<ConnType,numPol>::coo2C(connec[7]),-1,
             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[7]),OTT<ConnType,numPol>::coo2C(connec[1]),-1,
@@ -315,7 +424,7 @@ namespace INTERP_KERNEL
   }
 
   template<class ConnType, NumberingPolicy numPolConn>
-  void computeBarycenter2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim, double *res)
+  void computeBarycenter2(NormalizedCellType type, const ConnType *connec, mcIdType lgth, const double *coords, int spaceDim, double *res)
   {
     if(spaceDim==3)
       return computeBarycenter<ConnType,numPolConn,3>(type,connec,lgth,coords,res);