Salome HOME
PR: debug generate flat elements
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index 3fc0ca5a0e5460f0851af9b28ba69d753cfa08c8..49f9c1b3a8776abf98730188553103d7e1ca84fb 100644 (file)
@@ -95,6 +95,7 @@
 #include <numeric>
 #include <limits>
 #include <algorithm>
+#include <sstream>
 
 #define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
 
@@ -364,45 +365,55 @@ int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
   if ( aMesh->ShapeToMesh().IsNull() )
     return 0;
 
-  if ( theElem->GetType() == SMDSAbs_Node )
-    {
-      int aShapeID = theElem->getshapeId();
-      if (aShapeID <= 0)
-        return 0;
-      else
-        return aShapeID;
-    }
-
-  TopoDS_Shape aShape; // the shape a node is on
-  SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
-  while ( nodeIt->more() ) {
-    const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
-    int aShapeID = node->getshapeId();
-    if (aShapeID > 0) {
-      SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID );
-      if ( sm ) {
-        if ( sm->Contains( theElem ))
-          return aShapeID;
-        if ( aShape.IsNull() )
-          aShape = aMesh->IndexToShape( aShapeID );
-      }
-      else {
-        //MESSAGE ( "::FindShape() No SubShape for aShapeID " << aShapeID );
+  int aShapeID = theElem->getshapeId();
+  if ( aShapeID < 1 )
+    return 0;
+
+  if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ))
+    if ( sm->Contains( theElem ))
+      return aShapeID;
+
+  if ( theElem->GetType() == SMDSAbs_Node ) {
+    MESSAGE( ":( Error: invalid myShapeId of node " << theElem->GetID() );
+  }
+  else {
+    MESSAGE( ":( Error: invalid myShapeId of element " << theElem->GetID() );
+  }
+
+  TopoDS_Shape aShape; // the shape a node of theElem is on
+  if ( theElem->GetType() != SMDSAbs_Node )
+  {
+    SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
+    while ( nodeIt->more() ) {
+      const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+      if ((aShapeID = node->getshapeId()) > 0) {
+        if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ) ) {
+          if ( sm->Contains( theElem ))
+            return aShapeID;
+          if ( aShape.IsNull() )
+            aShape = aMesh->IndexToShape( aShapeID );
+        }
       }
     }
   }
 
   // None of nodes is on a proper shape,
   // find the shape among ancestors of aShape on which a node is
-  if ( aShape.IsNull() ) {
-    //MESSAGE ("::FindShape() - NONE node is on shape")
-    return 0;
+  if ( !aShape.IsNull() ) {
+    TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
+    for ( ; ancIt.More(); ancIt.Next() ) {
+      SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
+      if ( sm && sm->Contains( theElem ))
+        return aMesh->ShapeToIndex( ancIt.Value() );
+    }
   }
-  TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
-  for ( ; ancIt.More(); ancIt.Next() ) {
-    SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
-    if ( sm && sm->Contains( theElem ))
-      return aMesh->ShapeToIndex( ancIt.Value() );
+  else
+  {
+    const map<int,SMESHDS_SubMesh*>& id2sm = GetMeshDS()->SubMeshes();
+    map<int,SMESHDS_SubMesh*>::const_iterator id_sm = id2sm.begin();
+    for ( ; id_sm != id2sm.end(); ++id_sm )
+      if ( id_sm->second->Contains( theElem ))
+        return id_sm->first;
   }
 
   //MESSAGE ("::FindShape() - SHAPE NOT FOUND")
@@ -5307,7 +5318,7 @@ void SMESH_MeshEditor::LinearAngleVariation(const int nbSteps,
  *  \param theCopy - if true, create translated copies of theElems
  *  \param theMakeGroups - if true and theCopy, create translated groups
  *  \param theTargetMesh - mesh to copy translated elements into
- *  \retval SMESH_MeshEditor::PGroupIDs - list of ids of created groups
+ *  \return SMESH_MeshEditor::PGroupIDs - list of ids of created groups
  */
 //================================================================================
 
@@ -9184,7 +9195,7 @@ void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode*        theBetweenNode
 //=======================================================================
 /*!
  * \brief Convert elements contained in a submesh to quadratic
- * \retval int - nb of checked elements
+ * \return int - nb of checked elements
  */
 //=======================================================================
 
@@ -9536,7 +9547,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool        theForce3d,
 //=======================================================================
 /*!
  * \brief Convert quadratic elements to linear ones and remove quadratic nodes
- * \retval int - nb of checked elements
+ * \return int - nb of checked elements
  */
 //=======================================================================
 
@@ -10337,7 +10348,7 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
  * \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1
  * \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1
  * \param nReplaceMap - output map of corresponding nodes
- * \retval bool  - is a success or not
+ * \return bool  - is a success or not
  */
 //================================================================================
 
@@ -10847,6 +10858,7 @@ double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const
  * The nodes of the internal faces at the boundaries of the groups are doubled.
  * In option, the internal faces are replaced by flat elements.
  * Triangles are transformed in prisms, and quadrangles in hexahedrons.
+ * The flat elements are stored in groups of volumes.
  * @param theElems - list of groups of volumes, where a group of volume is a set of
  * SMDS_MeshElements sorted by Id.
  * @param createJointElems - if TRUE, create the elements
@@ -10860,7 +10872,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   MESSAGE("----------------------------------------------");
 
   SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
-  meshDS->BuildDownWardConnectivity(false);
+  meshDS->BuildDownWardConnectivity(true);
   CHRONO(50);
   SMDS_UnstructuredGrid *grid = meshDS->getGrid();
 
@@ -10956,6 +10968,12 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                     continue;
                   int vtkType = grid->GetCellType(vtkId);
                   int downId = grid->CellIdToDownId(vtkId);
+                  if (downId < 0)
+                    {
+                      MESSAGE("doubleNodesOnGroupBoundaries: internal algorithm problem");
+                      continue; // not OK at this stage of the algorithm:
+                                //no cells created after BuildDownWardConnectivity
+                    }
                   DownIdType aCell(downId, vtkType);
                   if (celldom.count(vtkId))
                     continue;
@@ -11033,6 +11051,11 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                       nodeDomains[oldId][idom] = newId; // cloned node for other domains
                       //MESSAGE("   newNode " << newId << " oldNode " << oldId << " size=" <<nodeDomains[oldId].size());
                     }
+                  if (nodeDomains[oldId].size() >= 3)
+                    {
+                      //MESSAGE("confirm multiple node " << oldId);
+                      isMultipleDetected =true;
+                    }
                 }
             }
           if (isMultipleDetected) // check if an edge of the face is shared between 3 or more domains
@@ -11125,11 +11148,12 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   //     (domain1 X domain2) = domain1 + MAXINT*domain2
 
   std::map<int, std::map<long,int> > nodeQuadDomains;
+  std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
 
   if (createJointElems)
     {
       itface = faceDomains.begin();
-      for( ; itface != faceDomains.end();++itface )
+      for (; itface != faceDomains.end(); ++itface)
         {
           DownIdType face = itface->first;
           std::set<int> oldNodes;
@@ -11137,13 +11161,27 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
           oldNodes.clear();
           grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
 
-          std::map<int,int> domvol = itface->second;
-          std::map<int,int>::iterator itdom = domvol.begin();
+          std::map<int, int> domvol = itface->second;
+          std::map<int, int>::iterator itdom = domvol.begin();
           int dom1 = itdom->first;
           int vtkVolId = itdom->second;
           itdom++;
           int dom2 = itdom->first;
-          grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains, nodeQuadDomains);
+          SMDS_MeshVolume *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
+                                                             nodeQuadDomains);
+          stringstream grpname;
+          grpname << "j_";
+          if (dom1 < dom2)
+            grpname << dom1 << "_" << dom2;
+          else
+            grpname << dom2 << "_" << dom1;
+          int idg;
+          string namegrp = grpname.str();
+          if (!mapOfJunctionGroups.count(namegrp))
+            mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+          SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+          if (sgrp)
+            sgrp->Add(vol->GetID());
         }
     }
 
@@ -11157,7 +11195,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
         {
           vector<int> nodes = ite->first;
           vector<int> orderDom = ite->second;
-          vector<int> orderedNodes;
+          vector<vtkIdType> orderedNodes;
           if (nodes.size() == 2)
             {
               //MESSAGE(" use edgesMultiDomains " << nodes[0] << " " << nodes[1]);
@@ -11168,10 +11206,22 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                 else
                   for (int idom = orderDom.size()-1; idom >=0; idom--)
                     orderedNodes.push_back( nodeDomains[nodes[ino]][orderDom[idom]] );
-              this->GetMeshDS()->AddVolumeFromVtkIds(orderedNodes);
+              SMDS_MeshVolume* vol = this->GetMeshDS()->AddVolumeFromVtkIds(orderedNodes);
+
+              stringstream grpname;
+              grpname << "mj_";
+              grpname << 0 << "_" << 0;
+              int idg;
+              string namegrp = grpname.str();
+              if (!mapOfJunctionGroups.count(namegrp))
+                mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+              SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+              if (sgrp)
+                sgrp->Add(vol->GetID());
             }
           else
             {
+              //MESSAGE("Quadratic multiple joints not implemented");
               // TODO quadratic nodes
             }
         }
@@ -11200,6 +11250,8 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
               int vtkId = l.cells[i];
               int vtkType = grid->GetCellType(vtkId);
               int downId = grid->CellIdToDownId(vtkId);
+              if (downId < 0)
+                continue; // new cells: not to be modified
               DownIdType aCell(downId, vtkType);
               int volParents[1000];
               int nbvol = grid->GetParentVolumes(volParents, vtkId);
@@ -11258,6 +11310,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
         }
     }
 
+  meshDS->CleanDownWardConnectivity(); // Mesh has been modified, downward connectivity is no more usable, free memory
   grid->BuildLinks();
 
   CHRONOSTOP(50);
@@ -11265,6 +11318,167 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   return true;
 }
 
+/*!
+ * \brief Double nodes on some external faces and create flat elements.
+ * Flat elements are mainly used by some types of mechanic calculations.
+ *
+ * Each group of the list must be constituted of faces.
+ * Triangles are transformed in prisms, and quadrangles in hexahedrons.
+ * @param theElems - list of groups of faces, where a group of faces is a set of
+ * SMDS_MeshElements sorted by Id.
+ * @return TRUE if operation has been completed successfully, FALSE otherwise
+ */
+bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector<TIDSortedElemSet>& theElems)
+{
+  MESSAGE("-------------------------------------------------");
+  MESSAGE("SMESH_MeshEditor::CreateFlatElementsOnFacesGroups");
+  MESSAGE("-------------------------------------------------");
+
+  SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
+
+  // --- For each group of faces
+  //     duplicate the nodes, create a flat element based on the face
+  //     replace the nodes of the faces by their clones
+
+  std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> clonedNodes;
+  std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> intermediateNodes;
+  clonedNodes.clear();
+  intermediateNodes.clear();
+  std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
+  mapOfJunctionGroups.clear();
+
+  for (int idom = 0; idom < theElems.size(); idom++)
+    {
+      const TIDSortedElemSet& domain = theElems[idom];
+      TIDSortedElemSet::const_iterator elemItr = domain.begin();
+      for (; elemItr != domain.end(); ++elemItr)
+        {
+          SMDS_MeshElement* anElem = (SMDS_MeshElement*) *elemItr;
+          SMDS_MeshFace* aFace = dynamic_cast<SMDS_MeshFace*> (anElem);
+          if (!aFace)
+            continue;
+          // MESSAGE("aFace=" << aFace->GetID());
+          bool isQuad = aFace->IsQuadratic();
+          vector<const SMDS_MeshNode*> ln0, ln1, ln2, ln3, ln4;
+
+          // --- clone the nodes, create intermediate nodes for non medium nodes of a quad face
+
+          SMDS_ElemIteratorPtr nodeIt = aFace->nodesIterator();
+          while (nodeIt->more())
+            {
+              const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*> (nodeIt->next());
+              bool isMedium = isQuad && (aFace->IsMediumNode(node));
+              if (isMedium)
+                ln2.push_back(node);
+              else
+                ln0.push_back(node);
+
+              const SMDS_MeshNode* clone = 0;
+              if (!clonedNodes.count(node))
+                {
+                  clone = meshDS->AddNode(node->X(), node->Y(), node->Z());
+                  clonedNodes[node] = clone;
+                }
+              else
+                clone = clonedNodes[node];
+
+              if (isMedium)
+                ln3.push_back(clone);
+              else
+                ln1.push_back(clone);
+
+              const SMDS_MeshNode* inter = 0;
+              if (isQuad && (!isMedium))
+                {
+                  if (!intermediateNodes.count(node))
+                    {
+                      inter = meshDS->AddNode(node->X(), node->Y(), node->Z());
+                      intermediateNodes[node] = inter;
+                    }
+                  else
+                    inter = intermediateNodes[node];
+                  ln4.push_back(inter);
+                }
+            }
+
+          // --- extrude the face
+
+          vector<const SMDS_MeshNode*> ln;
+          SMDS_MeshVolume* vol = 0;
+          vtkIdType aType = aFace->GetVtkType();
+          switch (aType)
+          {
+            case VTK_TRIANGLE:
+              vol = meshDS->AddVolume(ln0[2], ln0[1], ln0[0], ln1[2], ln1[1], ln1[0]);
+              // MESSAGE("vol prism " << vol->GetID());
+              ln.push_back(ln1[0]);
+              ln.push_back(ln1[1]);
+              ln.push_back(ln1[2]);
+              break;
+            case VTK_QUAD:
+              vol = meshDS->AddVolume(ln0[3], ln0[2], ln0[1], ln0[0], ln1[3], ln1[2], ln1[1], ln1[0]);
+              // MESSAGE("vol hexa " << vol->GetID());
+              ln.push_back(ln1[0]);
+              ln.push_back(ln1[1]);
+              ln.push_back(ln1[2]);
+              ln.push_back(ln1[3]);
+              break;
+            case VTK_QUADRATIC_TRIANGLE:
+              vol = meshDS->AddVolume(ln1[0], ln1[1], ln1[2], ln0[0], ln0[1], ln0[2], ln3[0], ln3[1], ln3[2],
+                                      ln2[0], ln2[1], ln2[2], ln4[0], ln4[1], ln4[2]);
+              // MESSAGE("vol quad prism " << vol->GetID());
+              ln.push_back(ln1[0]);
+              ln.push_back(ln1[1]);
+              ln.push_back(ln1[2]);
+              ln.push_back(ln3[0]);
+              ln.push_back(ln3[1]);
+              ln.push_back(ln3[2]);
+              break;
+            case VTK_QUADRATIC_QUAD:
+//              vol = meshDS->AddVolume(ln0[0], ln0[1], ln0[2], ln0[3], ln1[0], ln1[1], ln1[2], ln1[3],
+//                                      ln2[0], ln2[1], ln2[2], ln2[3], ln3[0], ln3[1], ln3[2], ln3[3],
+//                                      ln4[0], ln4[1], ln4[2], ln4[3]);
+              vol = meshDS->AddVolume(ln1[0], ln1[1], ln1[2], ln1[3], ln0[0], ln0[1], ln0[2], ln0[3],
+                                      ln3[0], ln3[1], ln3[2], ln3[3], ln2[0], ln2[1], ln2[2], ln2[3],
+                                      ln4[0], ln4[1], ln4[2], ln4[3]);
+              // MESSAGE("vol quad hexa " << vol->GetID());
+              ln.push_back(ln1[0]);
+              ln.push_back(ln1[1]);
+              ln.push_back(ln1[2]);
+              ln.push_back(ln1[3]);
+              ln.push_back(ln3[0]);
+              ln.push_back(ln3[1]);
+              ln.push_back(ln3[2]);
+              ln.push_back(ln3[3]);
+              break;
+            case VTK_POLYGON:
+              break;
+            default:
+              break;
+          }
+
+          if (vol)
+            {
+              stringstream grpname;
+              grpname << "jf_";
+              grpname << idom;
+              int idg;
+              string namegrp = grpname.str();
+              if (!mapOfJunctionGroups.count(namegrp))
+                mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+              SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+              if (sgrp)
+                sgrp->Add(vol->GetID());
+            }
+
+          // --- modify the face
+
+          aFace->ChangeNodes(&ln[0], ln.size());
+        }
+    }
+  return true;
+}
+
 //================================================================================
 /*!
  * \brief Generates skin mesh (containing 2D cells) from 3D mesh