Salome HOME
Management of the new kid of the block VTK_QUADRATIC_POLYGON in VTK
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingUMesh.cxx
index c04843cf7ec1b11e21997818d4425c4ed080cc2a..28490ba3b08f2d0de29e7261c5ad01d30adf628c 100644 (file)
@@ -30,6 +30,7 @@
 #include "BBTreeDst.txx"
 #include "SplitterTetra.hxx"
 #include "DirectedBoundingBox.hxx"
+#include "InterpKernelMatrixTools.hxx"
 #include "InterpKernelMeshQuality.hxx"
 #include "InterpKernelCellSimplify.hxx"
 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
@@ -418,12 +419,39 @@ MEDCouplingUMeshCellByTypeEntry *MEDCouplingUMesh::cellsByType()
  * Returns a set of all cell types available in \a this mesh.
  * \return std::set<INTERP_KERNEL::NormalizedCellType> - the set of cell types.
  * \warning this method does not throw any exception even if \a this is not defined.
+ * \sa MEDCouplingUMesh::getAllGeoTypesSorted
  */
 std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingUMesh::getAllGeoTypes() const
 {
   return _types;
 }
 
+/*!
+ * This method returns the sorted list of geometric types in \a this.
+ * Sorted means in the same order than the cells in \a this. A single entry in return vector means the maximal chunk of consecutive cells in \a this
+ * having the same geometric type. So a same geometric type can appear more than once if the cells are not sorted per geometric type.
+ *
+ * \throw if connectivity in \a this is not correctly defined.
+ *  
+ * \sa MEDCouplingMesh::getAllGeoTypes
+ */
+std::vector<INTERP_KERNEL::NormalizedCellType> MEDCouplingUMesh::getAllGeoTypesSorted() const
+{
+  std::vector<INTERP_KERNEL::NormalizedCellType> ret;
+  checkConnectivityFullyDefined();
+  int nbOfCells(getNumberOfCells());
+  if(nbOfCells==0)
+    return ret;
+  if(getMeshLength()<1)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getAllGeoTypesSorted : the connectivity in this seems invalid !");
+  const int *c(_nodal_connec->begin()),*ci(_nodal_connec_index->begin());
+  ret.push_back((INTERP_KERNEL::NormalizedCellType)c[*ci++]);
+  for(int i=1;i<nbOfCells;i++,ci++)
+    if(ret.back()!=((INTERP_KERNEL::NormalizedCellType)c[*ci]))
+      ret.push_back((INTERP_KERNEL::NormalizedCellType)c[*ci]);
+  return ret;
+}
+
 /*!
  * This method is a method that compares \a this and \a other.
  * This method compares \b all attributes, even names and component names.
@@ -2909,7 +2937,7 @@ int MEDCouplingUMesh::getNumberOfNodesInCell(int cellId) const
  *         describing the cell types. 
  *  \throw If the coordinates array is not set.
  *  \throw If the nodal connectivity of cells is not defined.
- *  \sa getAllTypes()
+ *  \sa getAllGeoTypes()
  */
 std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingUMesh::getTypesOfPart(const int *begin, const int *end) const
 {
@@ -5244,7 +5272,7 @@ void MEDCouplingUMesh::tessellate2DCurve(double eps)
  *          and \a this->getMeshDimension() != 3. 
  *  \throw If \a policy is not one of the four discussed above.
  *  \throw If the nodal connectivity of cells is not defined.
- * \sa MEDCouplingUMesh::tetrahedrize
+ * \sa MEDCouplingUMesh::tetrahedrize, MEDCoupling1SGTUMesh::sortHexa8EachOther
  */
 DataArrayInt *MEDCouplingUMesh::simplexize(int policy)
 {
@@ -5692,13 +5720,25 @@ void MEDCouplingUMesh::orientCorrectly2DCells(const double *vec, bool polyOnly)
       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)conn[connI[i]];
       if(!polyOnly || (type==INTERP_KERNEL::NORM_POLYGON || type==INTERP_KERNEL::NORM_QPOLYG))
         {
-          bool isQuadratic=INTERP_KERNEL::CellModel::GetCellModel(type).isQuadratic();
+          bool isQuadratic(INTERP_KERNEL::CellModel::GetCellModel(type).isQuadratic());
           if(!IsPolygonWellOriented(isQuadratic,vec,conn+connI[i]+1,conn+connI[i+1],coordsPtr))
             {
               isModified=true;
-              std::vector<int> tmp(connI[i+1]-connI[i]-2);
-              std::copy(conn+connI[i]+2,conn+connI[i+1],tmp.rbegin());
-              std::copy(tmp.begin(),tmp.end(),conn+connI[i]+2);
+              if(!isQuadratic)
+                {
+                  std::vector<int> tmp(connI[i+1]-connI[i]-2);
+                  std::copy(conn+connI[i]+2,conn+connI[i+1],tmp.rbegin());
+                  std::copy(tmp.begin(),tmp.end(),conn+connI[i]+2);
+                }
+              else
+                {
+                  int sz(((int)(connI[i+1]-connI[i]-1))/2);
+                  std::vector<int> tmp0(sz-1),tmp1(sz);
+                  std::copy(conn+connI[i]+2,conn+connI[i]+1+sz,tmp0.rbegin());
+                  std::copy(conn+connI[i]+1+sz,conn+connI[i+1],tmp1.rbegin());
+                  std::copy(tmp0.begin(),tmp0.end(),conn+connI[i]+2);
+                  std::copy(tmp1.begin(),tmp1.end(),conn+connI[i]+1+sz);
+                }
             }
         }
     }
@@ -6268,7 +6308,7 @@ std::vector<int> MEDCouplingUMesh::getDistributionOfTypes() const
   const int *connI=_nodal_connec_index->getConstPointer();
   const int *work=connI;
   int nbOfCells=getNumberOfCells();
-  std::size_t n=getAllTypes().size();
+  std::size_t n=getAllGeoTypes().size();
   std::vector<int> ret(3*n,-1); //ret[3*k+2]==-1 because it has no sense here
   std::set<INTERP_KERNEL::NormalizedCellType> types;
   for(std::size_t i=0;work!=connI+nbOfCells;i++)
@@ -6939,7 +6979,7 @@ DataArrayInt *MEDCouplingUMesh::convertCellArrayPerGeoType(const DataArrayInt *d
   const int *conn=_nodal_connec->getConstPointer();
   const int *connI=_nodal_connec_index->getConstPointer();
   int nbOfCells=getNumberOfCells();
-  std::set<INTERP_KERNEL::NormalizedCellType> types=getAllTypes();
+  std::set<INTERP_KERNEL::NormalizedCellType> types(getAllGeoTypes());
   int *tmp=new int[nbOfCells];
   for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator iter=types.begin();iter!=types.end();iter++)
     {
@@ -7176,6 +7216,56 @@ DataArrayDouble *MEDCouplingUMesh::getPartBarycenterAndOwner(const int *begin, c
   return ret;
 }
 
+/*!
+ * Returns a DataArrayDouble instance giving for each cell in \a this the equation of plane given by "a*X+b*Y+c*Z+d=0".
+ * So the returned instance will have 4 components and \c this->getNumberOfCells() tuples.
+ * So this method expects that \a this has a spaceDimension equal to 3 and meshDimension equal to 2.
+ * The computation of the plane equation is done using each time the 3 first nodes of 2D cells.
+ * This method is useful to detect 2D cells in 3D space that are not coplanar.
+ * 
+ * \return DataArrayDouble * - a new instance of DataArrayDouble having 4 components and a number of tuples equal to number of cells in \a this.
+ * \throw If spaceDim!=3 or meshDim!=2.
+ * \throw If connectivity of \a this is invalid.
+ * \throw If connectivity of a cell in \a this points to an invalid node.
+ */
+DataArrayDouble *MEDCouplingUMesh::computePlaneEquationOf3DFaces() const
+{
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
+  int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
+  if(getSpaceDimension()!=3 || getMeshDimension()!=2)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::computePlaneEquationOf3DFaces : This method must be applied on a mesh having meshDimension equal 2 and a spaceDimension equal to 3 !");
+  ret->alloc(nbOfCells,4);
+  double *retPtr(ret->getPointer());
+  const int *nodal(_nodal_connec->begin()),*nodalI(_nodal_connec_index->begin());
+  const double *coor(_coords->begin());
+  for(int i=0;i<nbOfCells;i++,nodalI++,retPtr+=4)
+    {
+      double matrix[16]={0,0,0,1,0,0,0,1,0,0,0,1,1,1,1,0},matrix2[16];
+      if(nodalI[1]-nodalI[0]>=3)
+        {
+          for(int j=0;j<3;j++)
+            {
+              int nodeId(nodal[nodalI[0]+1+j]);
+              if(nodeId>=0 && nodeId<nbOfNodes)
+                std::copy(coor+nodeId*3,coor+(nodeId+1)*3,matrix+4*j);
+              else
+                {
+                  std::ostringstream oss; oss << "MEDCouplingUMesh::computePlaneEquationOf3DFaces : invalid 2D cell #" << i << " ! This cell points to an invalid nodeId : " << nodeId << " !";
+                  throw INTERP_KERNEL::Exception(oss.str().c_str());
+                }
+            }
+        }
+      else
+        {
+          std::ostringstream oss; oss << "MEDCouplingUMesh::computePlaneEquationOf3DFaces : invalid 2D cell #" << i << " ! Must be constitued by more than 3 nodes !";
+          throw INTERP_KERNEL::Exception(oss.str().c_str());
+        }
+      INTERP_KERNEL::inverseMatrix(matrix,4,matrix2);
+      retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15];
+    }
+  return ret.retn();
+}
+
 /*!
  * This method expects as input a DataArrayDouble non nul instance 'da' that should be allocated. If not an exception is thrown.
  * 
@@ -8080,7 +8170,7 @@ void MEDCouplingUMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData
   int nbOfCells=getNumberOfCells();
   if(nbOfCells<=0)
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::writeVTK : the unstructured mesh has no cells !");
-  static const int PARAMEDMEM2VTKTYPETRADUCER[INTERP_KERNEL::NORM_MAXTYPE+1]={1,3,21,5,9,7,22,34,23,28,-1,-1,-1,-1,10,14,13,-1,12,-1,24,-1,16,27,-1,26,-1,29,-1,-1,25,42,-1,4};
+  static const int PARAMEDMEM2VTKTYPETRADUCER[INTERP_KERNEL::NORM_MAXTYPE+1]={1,3,21,5,9,7,22,34,23,28,-1,-1,-1,-1,10,14,13,-1,12,-1,24,-1,16,27,-1,26,-1,29,-1,-1,25,42,36,4};
   ofs << "  <" << getVTKDataSetType() << ">\n";
   ofs << "    <Piece NumberOfPoints=\"" << getNumberOfNodes() << "\" NumberOfCells=\"" << nbOfCells << "\">\n";
   ofs << "      <PointData>\n" << pointData << std::endl;
@@ -9328,10 +9418,13 @@ DataArrayInt *MEDCouplingUMesh::ComputeRangesFromTypeDistribution(const std::vec
  *          an id of old cell producing it. The caller is to delete this array using
  *         decrRef() as it is no more needed.
  * \return MEDCoupling1SGTUMesh * - the mesh containing only INTERP_KERNEL::NORM_TETRA4 cells.
+ *
+ * \warning This method operates on each cells in this independantly ! So it can leads to non conform mesh in returned value ! If you expect to have a conform mesh in output
+ * the policy PLANAR_FACE_6 should be used on a mesh sorted with MEDCoupling1SGTUMesh::sortHexa8EachOther.
  * 
  * \throw If \a this is not a 3D mesh (spaceDim==3 and meshDim==3).
  * \throw If \a this is not fully constituted with linear 3D cells.
- * \sa MEDCouplingUMesh::simplexize
+ * \sa MEDCouplingUMesh::simplexize, MEDCoupling1SGTUMesh::sortHexa8EachOther
  */
 MEDCoupling1SGTUMesh *MEDCouplingUMesh::tetrahedrize(int policy, DataArrayInt *& n2oCells, int& nbOfAdditionalPoints) const
 {