Salome HOME
Make the medcoupling to VTK accessible from the outside. For MEDWriter.
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingUMesh.cxx
index 043cfaaee13239d1b089039849c0329bb7b0b681..68144cd98eb7afe71b7834e3a01a95a45b4c0f5d 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2015  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2016  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
@@ -55,6 +55,7 @@ double MEDCouplingUMesh::EPS_FOR_POLYH_ORIENTATION=1.e-14;
 
 /// @cond INTERNAL
 const INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::MEDMEM_ORDER[N_MEDMEM_ORDER] = { INTERP_KERNEL::NORM_POINT1, INTERP_KERNEL::NORM_SEG2, INTERP_KERNEL::NORM_SEG3, INTERP_KERNEL::NORM_SEG4, INTERP_KERNEL::NORM_POLYL, INTERP_KERNEL::NORM_TRI3, INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_TRI7, INTERP_KERNEL::NORM_QUAD8, INTERP_KERNEL::NORM_QUAD9, INTERP_KERNEL::NORM_POLYGON, INTERP_KERNEL::NORM_QPOLYG, INTERP_KERNEL::NORM_TETRA4, INTERP_KERNEL::NORM_PYRA5, INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXA8, INTERP_KERNEL::NORM_HEXGP12, INTERP_KERNEL::NORM_TETRA10, INTERP_KERNEL::NORM_PYRA13, INTERP_KERNEL::NORM_PENTA15, INTERP_KERNEL::NORM_HEXA20, INTERP_KERNEL::NORM_HEXA27, INTERP_KERNEL::NORM_POLYHED };
+const int MEDCouplingUMesh::MEDCOUPLING2VTKTYPETRADUCER[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};
 /// @endcond
 
 MEDCouplingUMesh *MEDCouplingUMesh::New()
@@ -1277,6 +1278,8 @@ void MEDCouplingUMesh::convertExtrudedPolyhedra()
  * \warning Cells of the result mesh are \b not sorted by geometric type, hence,
  *          to write this mesh to the MED file, its cells must be sorted using
  *          sortCellsInMEDFileFrmt().
+ * \warning Cells (and most notably polyhedrons) must be correctly oriented for this to work
+ *          properly. See orientCorrectlyPolyhedrons() and arePolyhedronsNotCorrectlyOriented().
  * \return \c true if at least one cell has been converted, \c false else. In the
  *         last case the nodal connectivity remains unchanged.
  * \throw If the coordinates array is not set.
@@ -1898,7 +1901,7 @@ void MEDCouplingUMesh::FindCommonCellsAlg(int compType, int startCellId, const D
 /*!
  * Checks if \a this mesh includes all cells of an \a other mesh, and returns an array
  * giving for each cell of the \a other an id of a cell in \a this mesh. A value larger
- * than \a other->getNumberOfCells() in the returned array means that there is no
+ * than \a this->getNumberOfCells() in the returned array means that there is no
  * corresponding cell in \a this mesh.
  * It is expected that \a this and \a other meshes share the same node coordinates
  * array, if it is not so an exception is thrown. 
@@ -2437,6 +2440,7 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On
                                             DataArrayInt *& cellIdsNeededToBeRenum, DataArrayInt *& cellIdsNotModified) const
 {
   typedef MCAuto<DataArrayInt> DAInt;
+  typedef MCAuto<MEDCouplingUMesh> MCUMesh;
 
   checkFullyDefined();
   otherDimM1OnSameCoords.checkFullyDefined();
@@ -2444,21 +2448,62 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findNodesToDuplicate : meshes do not share the same coords array !");
   if(otherDimM1OnSameCoords.getMeshDimension()!=getMeshDimension()-1)
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findNodesToDuplicate : the mesh given in other parameter must have this->getMeshDimension()-1 !");
-  DataArrayInt *cellIdsRk0=0,*cellIdsRk1=0;
-  findCellIdsLyingOn(otherDimM1OnSameCoords,cellIdsRk0,cellIdsRk1);
-  DAInt cellIdsRk0Auto(cellIdsRk0),cellIdsRk1Auto(cellIdsRk1);
-  DAInt s0=cellIdsRk1->buildComplement(cellIdsRk0->getNumberOfTuples());
-  s0->transformWithIndArr(cellIdsRk0Auto->begin(),cellIdsRk0Auto->end());
-  MCAuto<MEDCouplingUMesh> m0Part=static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(s0->begin(),s0->end(),true));
-  DAInt s1=m0Part->computeFetchedNodeIds();
-  DAInt s2=otherDimM1OnSameCoords.computeFetchedNodeIds();
-  DAInt s3=s2->buildSubstraction(s1);
-  cellIdsRk1->transformWithIndArr(cellIdsRk0Auto->begin(),cellIdsRk0Auto->end());
+
+  // Checking star-shaped M1 group:
+  DAInt dt0=DataArrayInt::New(),dit0=DataArrayInt::New(),rdt0=DataArrayInt::New(),rdit0=DataArrayInt::New();
+  MCUMesh meshM2 = otherDimM1OnSameCoords.buildDescendingConnectivity(dt0, dit0, rdt0, rdit0);
+  DAInt dsi = rdit0->deltaShiftIndex();
+  DAInt idsTmp0 = dsi->findIdsNotInRange(-1, 3);
+  if(idsTmp0->getNumberOfTuples())
+    throw INTERP_KERNEL::Exception("MEDFileUMesh::buildInnerBoundaryAlongM1Group: group is too complex: some points (or edges) have more than two connected segments (or faces)!");
+  dt0=0; dit0=0; rdt0=0; rdit0=0; idsTmp0=0;
+
+  // Get extreme nodes from the group (they won't be duplicated), ie nodes belonging to boundary cells of M1
+  DAInt xtremIdsM2 = dsi->findIdsEqual(1); dsi = 0;
+  MCUMesh meshM2Part = static_cast<MEDCouplingUMesh *>(meshM2->buildPartOfMySelf(xtremIdsM2->begin(), xtremIdsM2->end(),true));
+  DAInt xtrem = meshM2Part->computeFetchedNodeIds();
+  // Remove from the list points on the boundary of the M0 mesh (those need duplication!)
+  dt0=DataArrayInt::New(),dit0=DataArrayInt::New(),rdt0=DataArrayInt::New(),rdit0=DataArrayInt::New();
+  MCUMesh m0desc = buildDescendingConnectivity(dt0, dit0, rdt0, rdit0); dt0=0; dit0=0; rdt0=0;
+  dsi = rdit0->deltaShiftIndex();
+  DAInt boundSegs = dsi->findIdsEqual(1);   // boundary segs/faces of the M0 mesh
+  MCUMesh m0descSkin = static_cast<MEDCouplingUMesh *>(m0desc->buildPartOfMySelf(boundSegs->begin(),boundSegs->end(), true));
+  DAInt fNodes = m0descSkin->computeFetchedNodeIds();
+  // In 3D, some points on the boundary of M0 still need duplication:
+  DAInt notDup = 0;
+  if (getMeshDimension() == 3)
+    {
+      DAInt dnu1=DataArrayInt::New(), dnu2=DataArrayInt::New(), dnu3=DataArrayInt::New(), dnu4=DataArrayInt::New();
+      MCUMesh m0descSkinDesc = m0descSkin->buildDescendingConnectivity(dnu1, dnu2, dnu3, dnu4);
+      dnu1=0;dnu2=0;dnu3=0;dnu4=0;
+      DataArrayInt * corresp=0;
+      meshM2->areCellsIncludedIn(m0descSkinDesc,2,corresp);
+      DAInt validIds = corresp->findIdsInRange(0, meshM2->getNumberOfCells());
+      corresp->decrRef();
+      if (validIds->getNumberOfTuples())
+        {
+          MCUMesh m1IntersecSkin = static_cast<MEDCouplingUMesh *>(m0descSkinDesc->buildPartOfMySelf(validIds->begin(), validIds->end(), true));
+          DAInt notDuplSkin = m1IntersecSkin->findBoundaryNodes();
+          DAInt fNodes1 = fNodes->buildSubstraction(notDuplSkin);
+          notDup = xtrem->buildSubstraction(fNodes1);
+        }
+      else
+        notDup = xtrem->buildSubstraction(fNodes);
+    }
+  else
+    notDup = xtrem->buildSubstraction(fNodes);
+
+  // Now compute cells around group (i.e. cells where we will do the propagation to identify the two sub-sets delimited by the group)
+  DAInt m1Nodes = otherDimM1OnSameCoords.computeFetchedNodeIds();
+  DAInt dupl = m1Nodes->buildSubstraction(notDup);
+  DAInt cellsAroundGroup = getCellIdsLyingOnNodes(dupl->begin(), dupl->end(), false);  // false= take cell in, even if not all nodes are in notDup
+
   //
-  MCAuto<MEDCouplingUMesh> m0Part2=static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(cellIdsRk1->begin(),cellIdsRk1->end(),true));
+  MCUMesh m0Part2=static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(cellsAroundGroup->begin(),cellsAroundGroup->end(),true));
   int nCells2 = m0Part2->getNumberOfCells();
   DAInt desc00=DataArrayInt::New(),descI00=DataArrayInt::New(),revDesc00=DataArrayInt::New(),revDescI00=DataArrayInt::New();
-  MCAuto<MEDCouplingUMesh> m01=m0Part2->buildDescendingConnectivity(desc00,descI00,revDesc00,revDescI00);
+  MCUMesh m01=m0Part2->buildDescendingConnectivity(desc00,descI00,revDesc00,revDescI00);
+
   // Neighbor information of the mesh without considering the crack (serves to count how many connex pieces it is made of)
   DataArrayInt *tmp00=0,*tmp11=0;
   MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(desc00,descI00,revDesc00,revDescI00, tmp00, tmp11);
@@ -2468,8 +2513,6 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On
   DataArrayInt *idsTmp=0;
   bool b=m01->areCellsIncludedIn(&otherDimM1OnSameCoords,2,idsTmp);
   DAInt ids(idsTmp);
-  if(!b)
-    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findNodesToDuplicate : the given mdim-1 mesh in other is not a constituent of this !");
   // In the neighbor information remove the connection between high dimension cells and its low level constituents which are part
   // of the frontier given in parameter (i.e. the cells of low dimension from the group delimiting the crack):
   MEDCouplingUMesh::RemoveIdsFromIndexedArrays(ids->begin(),ids->end(),desc00,descI00);
@@ -2515,12 +2558,12 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findNodesToDuplicate(): internal error - too many iterations.");
 
   DAInt cellsToModifyConn1_torenum=cellsToModifyConn0_torenum->buildComplement(neighI00->getNumberOfTuples()-1);
-  cellsToModifyConn0_torenum->transformWithIndArr(cellIdsRk1->begin(),cellIdsRk1->end());
-  cellsToModifyConn1_torenum->transformWithIndArr(cellIdsRk1->begin(),cellIdsRk1->end());
+  cellsToModifyConn0_torenum->transformWithIndArr(cellsAroundGroup->begin(),cellsAroundGroup->end());
+  cellsToModifyConn1_torenum->transformWithIndArr(cellsAroundGroup->begin(),cellsAroundGroup->end());
   //
   cellIdsNeededToBeRenum=cellsToModifyConn0_torenum.retn();
   cellIdsNotModified=cellsToModifyConn1_torenum.retn();
-  nodeIdsToDuplicate=s3.retn();
+  nodeIdsToDuplicate=dupl.retn();
 }
 
 /*!
@@ -8895,7 +8938,6 @@ 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,36,4};
   ofs << "  <" << getVTKDataSetType() << ">\n";
   ofs << "    <Piece NumberOfPoints=\"" << getNumberOfNodes() << "\" NumberOfCells=\"" << nbOfCells << "\">\n";
   ofs << "      <PointData>\n" << pointData << std::endl;
@@ -8938,7 +8980,7 @@ void MEDCouplingUMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData
           w4=std::copy(c.begin(),c.end(),w4);
         }
     }
-  types->transformWithIndArr(PARAMEDMEM2VTKTYPETRADUCER,PARAMEDMEM2VTKTYPETRADUCER+INTERP_KERNEL::NORM_MAXTYPE+1);
+  types->transformWithIndArr(MEDCOUPLING2VTKTYPETRADUCER,MEDCOUPLING2VTKTYPETRADUCER+INTERP_KERNEL::NORM_MAXTYPE+1);
   types->writeVTK(ofs,8,"UInt8","types",byteData);
   offsets->writeVTK(ofs,8,"Int32","offsets",byteData);
   if(szFaceOffsets!=0)