Salome HOME
Indices are stored as mcIdType type instead of int to support switch to 64bits indexing
[tools/medcoupling.git] / src / INTERP_KERNEL / SplitterTetra.txx
index 4bc794ae6d45e58c16be404d3f68759c9c2e597f..ce0136d99360f1b9c708b70eadd3b703e77ac738 100644 (file)
@@ -1,9 +1,9 @@
-// Copyright (C) 2007-2013  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
@@ -99,9 +99,12 @@ namespace INTERP_KERNEL
    * \param [in] tetraCorners  array 4*3 doubles containing corners of input tetrahedron (P0X,P0Y,P0Y,P1X,P1Y,P1Z,P2X,P2Y,P2Z,P3X,P3Y,P3Z).
    */
   template<class MyMeshType>
-  SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12]): _t(0),_src_mesh(srcMesh)
+  SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12], const ConnType *conn): _t(0),_src_mesh(srcMesh)
   {
-    _conn[0]=0; _conn[1]=1; _conn[2]=2; _conn[3]=3;
+    if(!conn)
+      { _conn[0]=0; _conn[1]=1; _conn[2]=2; _conn[3]=3; }
+    else
+      { _conn[0]=conn[0]; _conn[1]=conn[1]; _conn[2]=conn[2]; _conn[3]=conn[3]; }
     _coords[0]=tetraCorners[0]; _coords[1]=tetraCorners[1]; _coords[2]=tetraCorners[2]; _coords[3]=tetraCorners[3]; _coords[4]=tetraCorners[4]; _coords[5]=tetraCorners[5];
     _coords[6]=tetraCorners[6]; _coords[7]=tetraCorners[7]; _coords[8]=tetraCorners[8]; _coords[9]=tetraCorners[9]; _coords[10]=tetraCorners[10]; _coords[11]=tetraCorners[11];
     // create the affine transform
@@ -118,7 +121,7 @@ namespace INTERP_KERNEL
   SplitterTetra<MyMeshType>::~SplitterTetra()
   {
     delete _t;
-    for(HashMap< int, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
+    for(typename HashMap< ConnType, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
       delete[] iter->second;
   }
 
@@ -188,17 +191,17 @@ namespace INTERP_KERNEL
     // get type of cell
     NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(element));
     const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
-    unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(element)) : cellModelCell.getNumberOfNodes();
+    ConnType nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(element)) : cellModelCell.getNumberOfNodes();
     // halfspace filtering
     bool isOutside[8] = {true, true, true, true, true, true, true, true};
     bool isTargetOutside = false;
 
     // calculate the coordinates of the nodes
-    int *cellNodes=new int[nbOfNodes4Type];
-    for(int i = 0;i<(int)nbOfNodes4Type;++i)
+    ConnType *cellNodes=new ConnType[nbOfNodes4Type];
+    for(ConnType i = 0;i<nbOfNodes4Type;++i)
       {
         // we could store mapping local -> global numbers too, but not sure it is worth it
-        const int globalNodeNum = getGlobalNumberOfNode(i, OTT<ConnType,numPol>::indFC(element), _src_mesh);
+        const ConnType globalNodeNum = getGlobalNumberOfNode(i, OTT<ConnType,numPol>::indFC(element), _src_mesh);
         cellNodes[i]=globalNodeNum;
         if(_nodes.find(globalNodeNum) == _nodes.end()) 
           {
@@ -207,8 +210,7 @@ namespace INTERP_KERNEL
             //std::cout << std::endl << "*** " << globalNodeNum << std::endl;
             calculateNode(globalNodeNum);
           }
-
-        checkIsOutside(_nodes[globalNodeNum], isOutside);       
+        CheckIsOutside(_nodes[globalNodeNum], isOutside);       
       }
 
     // halfspace filtering check
@@ -230,19 +232,19 @@ namespace INTERP_KERNEL
 
         // get nb of sons of a cell
         const ConnType* rawCellConn = _src_mesh.getConnectivityPtr() + OTT<ConnType,numPol>::conn2C( _src_mesh.getConnectivityIndexPtr()[ element ]);
-        const int rawNbCellNodes = _src_mesh.getConnectivityIndexPtr()[ element+1 ] - _src_mesh.getConnectivityIndexPtr()[ element ];
+        const ConnType rawNbCellNodes = _src_mesh.getConnectivityIndexPtr()[ element+1 ] - _src_mesh.getConnectivityIndexPtr()[ element ];
         unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes);
 
         for(unsigned ii = 0 ; ii < nbOfSons; ++ii)
           {
             // get sons connectivity
             NormalizedCellType faceType;
-            int *faceNodes, nbFaceNodes=-1;
+            ConnType *faceNodes, nbFaceNodes=-1;
             if ( cellModelCell.isDynamic() )
               {
-                faceNodes=new int[nbOfNodes4Type];
+                faceNodes=new ConnType[nbOfNodes4Type];
                 nbFaceNodes = cellModelCell.fillSonCellNodalConnectivity2(ii,rawCellConn,rawNbCellNodes,faceNodes,faceType);
-                for ( int i = 0; i < nbFaceNodes; ++i )
+                for ( ConnType i = 0; i < nbFaceNodes; ++i )
                   faceNodes[i] = OTT<ConnType,numPol>::coo2C(faceNodes[i]);
               }
             else
@@ -250,7 +252,8 @@ namespace INTERP_KERNEL
                 faceType = cellModelCell.getSonType(ii);
                 const CellModel& faceModel=CellModel::GetCellModel(faceType);
                 assert(faceModel.getDimension() == 2);
-                faceNodes=new int[faceModel.getNumberOfNodes()];      
+                nbFaceNodes = cellModelCell.getNumberOfNodesConstituentTheSon(ii);
+                faceNodes = new ConnType[nbFaceNodes];
                 cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
               }
             // intersect a son with the unit tetra
@@ -323,8 +326,8 @@ namespace INTERP_KERNEL
 
               case NORM_POLYGON:
                 {
-                  int nbTria = nbFaceNodes - 2; // split polygon into nbTria triangles
-                  for ( int iTri = 0; iTri < nbTria; ++iTri )
+                  ConnType nbTria = nbFaceNodes - 2; // split polygon into nbTria triangles
+                  for ( ConnType iTri = 0; iTri < nbTria; ++iTri )
                     {
                       TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1+iTri], faceNodes[2+iTri]);
                       if(_volumes.find(key) == _volumes.end())
@@ -391,7 +394,6 @@ namespace INTERP_KERNEL
   {
     typedef typename MyMeshType::MyConnType ConnType;
     typedef double Vect2[2];
-    typedef double Vect3[3];
     typedef double Triangle2[3][2];
 
     const double *const tri0[3] = {p1, p2, p3};
@@ -587,16 +589,14 @@ namespace INTERP_KERNEL
    */
   template<class MyMeshType>
   double SplitterTetra<MyMeshType>::intersectSourceFace(const NormalizedCellType polyType,
-                                                        const int polyNodesNbr,
-                                                        const int *const polyNodes,
+                                                        const ConnType polyNodesNbr,
+                                                        const ConnType *const polyNodes,
                                                         const double *const *const polyCoords,
                                                         const double dimCaracteristic,
                                                         const double precision,
                                                         std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
                                                         std::set<TriangleFaceKey>& listOfTetraFacesColinear)
   {
-    typedef typename MyMeshType::MyConnType ConnType;
-
     double totalSurface = 0.0;
 
     // check if we have planar tetra element
@@ -614,16 +614,16 @@ namespace INTERP_KERNEL
     bool isTargetOutside = false;
 
     // calculate the coordinates of the nodes
-    for(int i = 0;i<(int)polyNodesNbr;++i)
+    for(ConnType i = 0;i<polyNodesNbr;++i)
       {
-        const int globalNodeNum = polyNodes[i];
+        const ConnType globalNodeNum = polyNodes[i];
         if(_nodes.find(globalNodeNum) == _nodes.end())
           {
             calculateNode2(globalNodeNum, polyCoords[i]);
           }
 
-        checkIsStrictlyOutside(_nodes[globalNodeNum], isStrictlyOutside, precision);
-        checkIsOutside(_nodes[globalNodeNum], isOutside, precision);
+        CheckIsStrictlyOutside(_nodes[globalNodeNum], isStrictlyOutside, precision);
+        CheckIsOutside(_nodes[globalNodeNum], isOutside, precision);
       }
 
     // halfspace filtering check
@@ -684,8 +684,8 @@ namespace INTERP_KERNEL
                     double planeConstant = dot(planeNormal, coordsTetraTriNode1);
                     if (IsFacesCoplanar(planeNormal, planeConstant, polyCoords, precision))
                       {
-                        int nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
-                        for (int iTri = 0; iTri < nbrPolyTri; ++iTri)
+                        ConnType nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
+                        for (ConnType iTri = 0; iTri < nbrPolyTri; ++iTri)
                           {
                             double volume = CalculateIntersectionSurfaceOfCoplanarTriangles(planeNormal,
                                                                                             planeConstant,
@@ -782,8 +782,8 @@ namespace INTERP_KERNEL
 
                 case NORM_POLYGON:
                   {
-                    int nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
-                    for (int iTri = 0; iTri < nbrPolyTri; ++iTri)
+                    ConnType nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
+                    for (ConnType iTri = 0; iTri < nbrPolyTri; ++iTri)
                       {
                         TriangleFaceKey key = TriangleFaceKey(polyNodes[0], polyNodes[1 + iTri], polyNodes[2 + iTri]);
                         if (_volumes.find(key) == _volumes.end())
@@ -848,7 +848,7 @@ namespace INTERP_KERNEL
     for(int i = 0;i<(int)nbOfNodes4Type;++i)
     {
       _t->apply(nodes[i], tetraCorners[i]);
-      checkIsOutside(nodes[i], isOutside);
+      CheckIsOutside(nodes[i], isOutside);
     }
 
     // halfspace filtering check
@@ -866,7 +866,7 @@ namespace INTERP_KERNEL
     if(!isTargetOutside)
     {
       const CellModel& cellModelCell=CellModel::GetCellModel(NORM_TETRA4);
-      int cellNodes[4] = { 0, 1, 2, 3 }, faceNodes[3];
+      ConnType cellNodes[4] = { 0, 1, 2, 3 }, faceNodes[3];
 
       for(unsigned ii = 0 ; ii < 4 ; ++ii)
       {
@@ -909,9 +909,9 @@ namespace INTERP_KERNEL
   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::releaseArrays()
   {
     // free potential sub-mesh nodes that have been allocated
-    typename MyMeshTypeT::MyConnType nbOfNodesT = _node_ids.size();// Issue 0020634.
-    if((int)_nodes.size()>=/*8*/nbOfNodesT)
+    if(_nodes.size()>=/*8*/_node_ids.size())
       {
+        typename MyMeshTypeT::MyConnType nbOfNodesT = static_cast<typename MyMeshTypeT::MyConnType>(_node_ids.size());
         std::vector<const double*>::iterator iter = _nodes.begin() + /*8*/nbOfNodesT;
         while(iter != _nodes.end())
           {
@@ -929,20 +929,23 @@ namespace INTERP_KERNEL
   template<class MyMeshTypeT, class MyMeshTypeS>
   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitTargetCell2(typename MyMeshTypeT::MyConnType targetCell, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
   {
-    const int *refConn(_target_mesh.getConnectivityPtr());
-    const int *cellConn(refConn+_target_mesh.getConnectivityIndexPtr()[targetCell]);
+    typedef typename MyMeshTypeT::MyConnType TConnType;
+    const TConnType *refConn(_target_mesh.getConnectivityPtr());
+    const TConnType *cellConn(refConn+_target_mesh.getConnectivityIndexPtr()[targetCell]);
     INTERP_KERNEL::NormalizedCellType gt(_target_mesh.getTypeOfElement(targetCell));
-    std::vector<int> tetrasNodalConn;
+    std::vector<TConnType> tetrasNodalConn;
     std::vector<double> addCoords;
     const double *coords(_target_mesh.getCoordinatesPtr());
     SplitIntoTetras(_splitting_pol,gt,cellConn,refConn+_target_mesh.getConnectivityIndexPtr()[targetCell+1],coords,tetrasNodalConn,addCoords);
     std::size_t nbTetras(tetrasNodalConn.size()/4); tetra.resize(nbTetras);
     double tmp[12];
+    typename MyMeshTypeS::MyConnType tmp2[4];
     for(std::size_t i=0;i<nbTetras;i++)
       {
         for(int j=0;j<4;j++)
           {
-            int cellId(tetrasNodalConn[4*i+j]);
+            TConnType cellId(tetrasNodalConn[4*i+j]);
+            tmp2[j]=cellId;
             if(cellId>=0)
               {
                 tmp[j*3+0]=coords[3*cellId+0];
@@ -956,7 +959,7 @@ namespace INTERP_KERNEL
                 tmp[j*3+2]=addCoords[3*(-cellId-1)+2];
               }
           }
-        tetra[i]=new SplitterTetra<MyMeshTypeS>(_src_mesh,tmp);
+        tetra[i]=new SplitterTetra<MyMeshTypeS>(_src_mesh,tmp,tmp2);
       }
   }
 
@@ -978,7 +981,7 @@ namespace INTERP_KERNEL
         _node_ids.resize(8);
         tetra.reserve(1);
         const double *nodes[4];
-        int conn[4];
+        ConnType conn[4];
         for(int node = 0; node < 4 ; ++node)
           {
             nodes[node]=getCoordsOfNode2(node, OTT<ConnType,numPol>::indFC(targetCell),_target_mesh,conn[node]);
@@ -1059,7 +1062,7 @@ namespace INTERP_KERNEL
     for(int i = 0; i < 5; ++i)
       {
         const double* nodes[4];
-        int conn[4];
+        typename MyMeshTypeS::MyConnType conn[4];
         for(int j = 0; j < 4; ++j)
           {
             conn[j] = subZone[ SPLIT_NODES_5[4*i+j] ];
@@ -1083,7 +1086,7 @@ namespace INTERP_KERNEL
     for(int i = 0; i < 6; ++i)
       {
         const double* nodes[4];
-        int conn[4];
+        typename MyMeshTypeS::MyConnType conn[4];
         for(int j = 0; j < 4; ++j)
           {
             conn[j] = subZone[SPLIT_NODES_6[4*i+j]];
@@ -1106,11 +1109,11 @@ namespace INTERP_KERNEL
   {
     // The two nodes of the original mesh cell used in each tetrahedron.
     // The tetrahedra all have nodes (cellCenter, faceCenter, edgeNode1, edgeNode2)
-    // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
+    // For the correspondence of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
     
     // nodes to use for tetrahedron
     const double* nodes[4];
-    int conn[4];
+    typename MyMeshTypeS::MyConnType conn[4];
     // get the cell center
     conn[0] = 14;
     nodes[0] = getCoordsOfSubNode(conn[0]);
@@ -1170,7 +1173,7 @@ namespace INTERP_KERNEL
     
     // create tetrahedra
     const double* nodes[4];
-    int conn[4];
+    typename MyMeshTypeS::MyConnType conn[4];
     for(int i = 0; i < 2; ++i)
       {
         for(int j = 0; j < 4; ++j)
@@ -1196,22 +1199,24 @@ namespace INTERP_KERNEL
     // get type of cell and nb of cell nodes
     NormalizedCellType normCellType=_target_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(targetCell));
     const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
-    unsigned nbOfCellNodes=cellModelCell.isDynamic() ? _target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell)) : cellModelCell.getNumberOfNodes();
+    ConnType nbOfCellNodes=cellModelCell.isDynamic() ? _target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell)) : cellModelCell.getNumberOfNodes();
 
     // get nb of cell sons (faces)
     const ConnType* rawCellConn = _target_mesh.getConnectivityPtr() + OTT<ConnType,numPol>::conn2C( _target_mesh.getConnectivityIndexPtr()[ targetCell ]);
-    const int rawNbCellNodes = _target_mesh.getConnectivityIndexPtr()[ targetCell+1 ] - _target_mesh.getConnectivityIndexPtr()[ targetCell ];
+    const ConnType rawNbCellNodes = _target_mesh.getConnectivityIndexPtr()[ targetCell+1 ] - _target_mesh.getConnectivityIndexPtr()[ targetCell ];
     unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes);
 
     // indices of nodes of a son
-    static std::vector<int> allNodeIndices; // == 0,1,2,...,nbOfCellNodes-1
-    while ( allNodeIndices.size() < nbOfCellNodes )
-      allNodeIndices.push_back( allNodeIndices.size() );
-    std::vector<int> classicFaceNodes(4);
-    int* faceNodes = cellModelCell.isDynamic() ? &allNodeIndices[0] : &classicFaceNodes[0];
+    static std::vector<ConnType> allNodeIndices; // == 0,1,2,...,nbOfCellNodes-1
+    while ( allNodeIndices.size() < (std::size_t)nbOfCellNodes )
+      allNodeIndices.push_back( static_cast<ConnType>(allNodeIndices.size()) );
+    std::vector<ConnType> classicFaceNodes(4);
+    if(cellModelCell.isQuadratic())
+      throw INTERP_KERNEL::Exception("SplitterTetra2::splitConvex : quadratic 3D cells are not implemented yet !");
+    ConnType* faceNodes = cellModelCell.isDynamic() ? &allNodeIndices[0] : &classicFaceNodes[0];
 
     // nodes of tetrahedron
-    int conn[4];
+    typename MyMeshTypeS::MyConnType conn[4];
     const double* nodes[4];
     nodes[3] = getCoordsOfSubNode2( nbOfCellNodes,conn[3]); // barycenter
 
@@ -1256,7 +1261,7 @@ namespace INTERP_KERNEL
   {
     // retrieve real mesh nodes
     
-    typename MyMeshTypeT::MyConnType nbOfNodesT = _node_ids.size();// Issue 0020634. _node_ids.resize(8);
+    typename MyMeshTypeT::MyConnType nbOfNodesT = static_cast<typename MyMeshTypeT::MyConnType>(_node_ids.size());// Issue 0020634. _node_ids.resize(8);
     for(int node = 0; node < nbOfNodesT ; ++node)
       {
         // calculate only normal nodes