Salome HOME
Update copyrights
[tools/medcoupling.git] / src / INTERP_KERNEL / SplitterTetra.txx
index 9312b5fa6eb8d31dfc81259d5e712d3eab7825d2..83b2628f9840640d53cfd9c89ea1502bc60332b6 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
@@ -84,6 +84,32 @@ namespace INTERP_KERNEL
     // create the affine transform
     _t=new TetraAffineTransform(_coords);
   }
+
+  /**
+   * SplitterTetra class computes for a list of cell ids of a given mesh \a srcMesh (badly named) the intersection with a 
+   * single TETRA4 cell given by \a tetraCorners (of length 4) and \a nodesId (of length 4 too). \a nodedIds is given only to establish
+   * if a partial computation of a triangle has already been performed (to increase performance).
+   *
+   * The \a srcMesh can contain polyhedron cells.
+   * 
+   * 
+   * Constructor creating object from the four corners of the tetrahedron.
+   *
+   * \param [in] srcMesh       mesh containing the source elements
+   * \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], const int *conn): _t(0),_src_mesh(srcMesh)
+  {
+    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
+    _t=new TetraAffineTransform(_coords);
+  }
   
   /**
    * Destructor
@@ -184,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
@@ -227,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 int[nbFaceNodes];
                 cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
               }
             // intersect a son with the unit tetra
@@ -368,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};
@@ -557,7 +582,6 @@ namespace INTERP_KERNEL
    * @param polyNodesNbr number of the nodes of the polygon source face
    * @param polyNodes numbers of the nodes of the polygon source face
    * @param polyCoords coordinates of the nodes of the polygon source face
-   * @param polyCoords coordinates of the nodes of the polygon source face
    * @param dimCaracteristic characteristic size of the meshes containing the triangles
    * @param precision precision for double float data used for comparison
    * @param listOfTetraFacesTreated list of tetra faces treated
@@ -573,8 +597,6 @@ namespace INTERP_KERNEL
                                                         std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
                                                         std::set<TriangleFaceKey>& listOfTetraFacesColinear)
   {
-    typedef typename MyMeshType::MyConnType ConnType;
-
     double totalSurface = 0.0;
 
     // check if we have planar tetra element
@@ -600,8 +622,8 @@ namespace INTERP_KERNEL
             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
@@ -826,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
@@ -899,6 +921,46 @@ namespace INTERP_KERNEL
       }
     _nodes.clear();
   }
+  
+  /*!
+   * \param [in] targetCell in C mode.
+   * \param [out] tetra is the output result tetra containers.
+   */
+  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]);
+    INTERP_KERNEL::NormalizedCellType gt(_target_mesh.getTypeOfElement(targetCell));
+    std::vector<int> 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];
+    int tmp2[4];
+    for(std::size_t i=0;i<nbTetras;i++)
+      {
+        for(int j=0;j<4;j++)
+          {
+            int cellId(tetrasNodalConn[4*i+j]);
+            tmp2[j]=cellId;
+            if(cellId>=0)
+              {
+                tmp[j*3+0]=coords[3*cellId+0];
+                tmp[j*3+1]=coords[3*cellId+1];
+                tmp[j*3+2]=coords[3*cellId+2];
+              }
+            else
+              {
+                tmp[j*3+0]=addCoords[3*(-cellId-1)+0];
+                tmp[j*3+1]=addCoords[3*(-cellId-1)+1];
+                tmp[j*3+2]=addCoords[3*(-cellId-1)+2];
+              }
+          }
+        tetra[i]=new SplitterTetra<MyMeshTypeS>(_src_mesh,tmp,tmp2);
+      }
+  }
 
   /*!
    * @param targetCell in C mode.
@@ -1046,7 +1108,7 @@ 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];
@@ -1085,25 +1147,10 @@ namespace INTERP_KERNEL
    */
   template<class MyMeshTypeT, class MyMeshTypeS>
   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
-  {
-    // Define 8 hexahedral subzones as in Grandy, p449
-    // the values correspond to the nodes that correspond to nodes 1,2,3,4,5,6,7,8 in the subcell
-    // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
-    static const int subZones[64] = 
-      {
-        0,8,21,12,9,20,26,22,
-        8,1,13,21,20,10,23,26,
-        12,21,16,3,22,26,25,17,
-        21,13,2,16,26,23,18,25,
-        9,20,26,22,4,11,24,14,
-        20,10,23,26,11,5,15,24,
-        22,26,25,17,14,24,19,7,
-        26,23,18,25,24,15,6,19
-      };
-    
+  { 
     for(int i = 0; i < 8; ++i)
       {
-        sixSplit(&subZones[8*i],tetra);
+        sixSplit(GENERAL_48_SUBZONES+8*i,tetra);
       }
   }
   
@@ -1163,6 +1210,8 @@ namespace INTERP_KERNEL
     while ( allNodeIndices.size() < nbOfCellNodes )
       allNodeIndices.push_back( allNodeIndices.size() );
     std::vector<int> classicFaceNodes(4);
+    if(cellModelCell.isQuadratic())
+      throw INTERP_KERNEL::Exception("SplitterTetra2::splitConvex : quadratic 3D cells are not implemented yet !");
     int* faceNodes = cellModelCell.isDynamic() ? &allNodeIndices[0] : &classicFaceNodes[0];
 
     // nodes of tetrahedron
@@ -1238,33 +1287,6 @@ namespace INTERP_KERNEL
 
           case GENERAL_48:
             {
-              // Each sub-node is the barycenter of two other nodes.
-              // For the edges, these lie on the original mesh.
-              // For the faces, these are the edge sub-nodes.
-              // For the cell these are two face sub-nodes.
-              static const int GENERAL_48_SUB_NODES[38] = 
-                {
-                  0,1,   // sub-node 9 (edge)
-                  0,4,   // sub-node 10 (edge)
-                  1,5,   // sub-node 11 (edge)
-                  4,5,   // sub-node 12 (edge)
-                  0,3,   // sub-node 13 (edge)
-                  1,2,   // sub-node 14 (edge)
-                  4,7,   // sub-node 15 (edge)
-                  5,6,   // sub-node 16 (edge)
-                  2,3,   // sub-node 17 (edge)
-                  3,7,   // sub-node 18 (edge)
-                  2,6,   // sub-node 19 (edge)
-                  6,7,   // sub-node 20 (edge)
-                  8,11,  // sub-node 21 (face)
-                  12,13, // sub-node 22 (face)
-                  9,17,  // sub-node 23 (face)
-                  10,18, // sub-node 24 (face)
-                  14,15, // sub-node 25 (face)
-                  16,19, // sub-node 26 (face)
-                  20,25  // sub-node 27 (cell)
-                };
-
               for(int i = 0; i < 19; ++i)
                 {
                   double* barycenter = new double[3];