Salome HOME
PR: big problem of performance in AffectedElemGroupsInRegion
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index 6d3e74d4e62167dc7d24cf358ebb756b1a6bfe30..c6605815ae5311e4e47a368d72882e0799c9edae 100644 (file)
 using namespace std;
 using namespace SMESH::Controls;
 
-typedef map<const SMDS_MeshElement*, list<const SMDS_MeshNode*> >    TElemOfNodeListMap;
-typedef map<const SMDS_MeshElement*, list<const SMDS_MeshElement*> > TElemOfElemListMap;
-
 typedef SMDS_SetIterator< SMDS_pElement, TIDSortedElemSet::const_iterator> TSetIterator;
 
 //=======================================================================
@@ -6028,8 +6025,8 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
     const SMDSAbs_ElementType type = groupDS->GetType();
     SMESHDS_Group* newGroup    = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(), type );
     SMESHDS_Group* newTopGroup = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(), type );
-    groupsByType[ groupDS->GetType() ].push_back( make_tuple( groupDS, newGroup, newTopGroup ));
-    orderedOldNewGroups.push_back( & groupsByType[ groupDS->GetType() ].back() );
+    groupsByType[ type ].push_back( make_tuple( groupDS, newGroup, newTopGroup ));
+    orderedOldNewGroups.push_back( & groupsByType[ type ].back() );
   }
 
   // Loop on nodes and elements to add them in new groups
@@ -6065,7 +6062,26 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
           if ( resElem != sourceElem )
             resultElems.push_back( resElem );
 
-      // add resultElems to groups made by ones the sourceElem belongs to
+      // there must be a top element
+      const SMDS_MeshElement* topElem = 0;
+      if ( isNodes )
+      {
+        topElem = resultElems.back();
+        resultElems.pop_back();
+      }
+      else
+      {
+        list< const SMDS_MeshElement* >::reverse_iterator resElemIt = resultElems.rbegin();
+        for ( ; resElemIt != resultElems.rend() ; ++resElemIt )
+          if ( (*resElemIt)->GetType() == sourceElem->GetType() )
+          {
+            topElem = *resElemIt;
+            resultElems.erase( --(resElemIt.base()) ); // erase *resElemIt
+            break;
+          }
+      }
+
+      // add resultElems to groups originted from ones the sourceElem belongs to
       list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end();
       for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew )
       {
@@ -6074,19 +6090,15 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
         {
           // fill in a new group
           SMDS_MeshGroup & newGroup = gOldNew->get<1>()->SMDSGroup();
-          list< const SMDS_MeshElement* > rejectedElems; // elements of other type
           list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt;
           for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt )
-            if ( !newGroup.Add( *resElemIt ))
-              rejectedElems.push_back( *resElemIt );
+            newGroup.Add( *resElemIt );
 
-          // fill "top" group
-          if ( !rejectedElems.empty() )
+          // fill "top" group
+          if ( topElem )
           {
             SMDS_MeshGroup & newTopGroup = gOldNew->get<2>()->SMDSGroup();
-            resLast = rejectedElems.end();
-            for ( resElemIt = rejectedElems.begin(); resElemIt != resLast; ++resElemIt )
-              !newTopGroup.Add( *resElemIt );
+            newTopGroup.Add( topElem );
           }
         }
       }
@@ -6116,7 +6128,9 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
 
         // make a name
         const bool isTop = ( nbNewGroups == 2 &&
-                             newGroupDS->GetType() == oldGroupDS->GetType() );
+                             newGroupDS->GetType() == oldGroupDS->GetType() &&
+                             is2nd );
+
         string name = oldGroupDS->GetStoreName();
         if ( !targetMesh ) {
           string suffix = ( isTop ? "top": postfix.c_str() );
@@ -11103,12 +11117,14 @@ namespace {
 
 //================================================================================
 /*!
-  \brief Identify the elements that will be affected by node duplication (actual duplication is not performed.
+  \brief Identify the elements that will be affected by node duplication (actual duplication is not performed).
   This method is the first step of DoubleNodeElemGroupsInRegion.
   \param theElems - list of groups of elements (edges or faces) to be replicated
   \param theNodesNot - list of groups of nodes not to replicated
   \param theShape - shape to detect affected elements (element which geometric center
-         located on or inside shape).
+         located on or inside shape). If the shape is null, detection is done on faces orientations
+         (select elements with a gravity center on the side given by faces normals).
+         This mode (null shape) is faster, but works only when theElems are faces, with coherents orientations.
          The replicated nodes should be associated to affected elements.
   \return groups of affected elements
   \sa DoubleNodeElemGroupsInRegion()
@@ -11121,44 +11137,145 @@ bool SMESH_MeshEditor::AffectedElemGroupsInRegion( const TIDSortedElemSet& theEl
                                                    TIDSortedElemSet&       theAffectedElems)
 {
   if ( theShape.IsNull() )
-    return false;
-
-  const double aTol = Precision::Confusion();
-  auto_ptr< BRepClass3d_SolidClassifier> bsc3d;
-  auto_ptr<_FaceClassifier>              aFaceClassifier;
-  if ( theShape.ShapeType() == TopAbs_SOLID )
   {
-    bsc3d.reset( new BRepClass3d_SolidClassifier(theShape));;
-    bsc3d->PerformInfinitePoint(aTol);
-  }
-  else if (theShape.ShapeType() == TopAbs_FACE )
-  {
-    aFaceClassifier.reset( new _FaceClassifier(TopoDS::Face(theShape)));
-  }
+    std::set<const SMDS_MeshNode*> alreadyCheckedNodes;
+    std::set<const SMDS_MeshElement*> alreadyCheckedElems;
+    std::set<const SMDS_MeshElement*> edgesToCheck;
+    alreadyCheckedNodes.clear();
+    alreadyCheckedElems.clear();
+    edgesToCheck.clear();
 
-  // iterates on indicated elements and get elements by back references from their nodes
-  TIDSortedElemSet::const_iterator elemItr = theElems.begin();
-  for ( ;  elemItr != theElems.end(); ++elemItr )
+    // --- iterates on elements to be replicated and get elements by back references from their nodes
+
+    TIDSortedElemSet::const_iterator elemItr = theElems.begin();
+    int ielem = 1;
+    for ( ;  elemItr != theElems.end(); ++elemItr )
+    {
+      SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr;
+      if (!anElem || (anElem->GetType() != SMDSAbs_Face))
+        continue;
+      gp_XYZ normal;
+      SMESH_Algo::FaceNormal( anElem, normal, /*normalized=*/true );
+      MESSAGE("element " << ielem++ <<  " normal " << normal.X() << " " << normal.Y() << " " << normal.Z());
+      std::set<const SMDS_MeshNode*> nodesElem;
+      nodesElem.clear();
+      SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
+      while ( nodeItr->more() )
+      {
+        const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
+        nodesElem.insert(aNode);
+      }
+      std::set<const SMDS_MeshNode*>::iterator nodit = nodesElem.begin();
+      for (; nodit != nodesElem.end(); nodit++)
+      {
+        MESSAGE("  noeud ");
+        const SMDS_MeshNode* aNode = *nodit;
+        if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() )
+          continue;
+        if (alreadyCheckedNodes.find(aNode) != alreadyCheckedNodes.end())
+          continue;
+        alreadyCheckedNodes.insert(aNode);
+        SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator();
+        while ( backElemItr->more() )
+        {
+          MESSAGE("    backelem ");
+          const SMDS_MeshElement* curElem = backElemItr->next();
+          if (alreadyCheckedElems.find(curElem) != alreadyCheckedElems.end())
+            continue;
+          if (theElems.find(curElem) != theElems.end())
+            continue;
+          alreadyCheckedElems.insert(curElem);
+          double x=0, y=0, z=0;
+          int nb = 0;
+          SMDS_ElemIteratorPtr nodeItr2 = curElem->nodesIterator();
+          while ( nodeItr2->more() )
+          {
+            const SMDS_MeshNode* anotherNode = cast2Node(nodeItr2->next());
+            x += anotherNode->X();
+            y += anotherNode->Y();
+            z += anotherNode->Z();
+            nb++;
+          }
+          gp_XYZ p;
+          p.SetCoord( x/nb -aNode->X(),
+                      y/nb -aNode->Y(),
+                      z/nb -aNode->Z() );
+          MESSAGE("      check " << p.X() << " " << p.Y() << " " << p.Z());
+          if (normal*p > 0)
+          {
+            MESSAGE("    --- inserted")
+            theAffectedElems.insert( curElem );
+          }
+          else if (curElem->GetType() == SMDSAbs_Edge)
+            edgesToCheck.insert(curElem);
+        }
+      }
+    }
+    // --- add also edges lying on the set of faces (all nodes in alreadyCheckedNodes)
+    std::set<const SMDS_MeshElement*>::iterator eit = edgesToCheck.begin();
+    for( ; eit != edgesToCheck.end(); eit++)
+    {
+      bool onside = true;
+      const SMDS_MeshElement* anEdge = *eit;
+      SMDS_ElemIteratorPtr nodeItr = anEdge->nodesIterator();
+      while ( nodeItr->more() )
+      {
+        const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
+        if (alreadyCheckedNodes.find(aNode) == alreadyCheckedNodes.end())
+        {
+          onside = false;
+          break;
+        }
+      }
+      if (onside)
+      {
+        MESSAGE("    --- edge onside inserted")
+        theAffectedElems.insert(anEdge);
+      }
+    }
+  }
+  else
   {
-    SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr;
-    if (!anElem)
-      continue;
+    const double aTol = Precision::Confusion();
+    auto_ptr< BRepClass3d_SolidClassifier> bsc3d;
+    auto_ptr<_FaceClassifier>              aFaceClassifier;
+    if ( theShape.ShapeType() == TopAbs_SOLID )
+    {
+      bsc3d.reset( new BRepClass3d_SolidClassifier(theShape));;
+      bsc3d->PerformInfinitePoint(aTol);
+    }
+    else if (theShape.ShapeType() == TopAbs_FACE )
+    {
+      aFaceClassifier.reset( new _FaceClassifier(TopoDS::Face(theShape)));
+    }
 
-    SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
-    while ( nodeItr->more() )
+    // iterates on indicated elements and get elements by back references from their nodes
+    TIDSortedElemSet::const_iterator elemItr = theElems.begin();
+    int ielem = 1;
+    for ( ;  elemItr != theElems.end(); ++elemItr )
     {
-      const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
-      if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() )
+      MESSAGE("element " << ielem++);
+      SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr;
+      if (!anElem)
         continue;
-      SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator();
-      while ( backElemItr->more() )
+      SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
+      while ( nodeItr->more() )
       {
-        const SMDS_MeshElement* curElem = backElemItr->next();
-        if ( curElem && theElems.find(curElem) == theElems.end() &&
-             ( bsc3d.get() ?
-               isInside( curElem, *bsc3d, aTol ) :
-               isInside( curElem, *aFaceClassifier, aTol )))
-          theAffectedElems.insert( curElem );
+        MESSAGE("  noeud ");
+        const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
+        if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() )
+          continue;
+        SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator();
+        while ( backElemItr->more() )
+        {
+          MESSAGE("    backelem ");
+          const SMDS_MeshElement* curElem = backElemItr->next();
+          if ( curElem && theElems.find(curElem) == theElems.end() &&
+              ( bsc3d.get() ?
+                isInside( curElem, *bsc3d, aTol ) :
+                isInside( curElem, *aFaceClassifier, aTol )))
+            theAffectedElems.insert( curElem );
+        }
       }
     }
   }
@@ -11251,11 +11368,14 @@ double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const
 
 /*!
  * \brief Double nodes on shared faces between groups of volumes and create flat elements on demand.
- * The list of groups must describe a partition of the mesh volumes.
- * 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.
+ *  The list of groups must contain at least two groups. The groups have to be disjoint: no common element into two different groups.
+ * The nodes of the internal faces at the boundaries of the groups are doubled. Optionally, the internal faces are replaced by flat elements.
+ * Triangles are transformed into prisms, and quadrangles into hexahedrons.
+ * The flat elements are stored in groups of volumes. These groups are named according to the position of the group in the list:
+ * the group j_n_p is the group of the flat elements that are built between the group #n and the group #p in the list.
+ * If there is no shared faces between the group #n and the group #p in the list, the group j_n_p is not created.
+ * All the flat elements are gathered into the group named "joints3D" (or "joints2D" in 2D situation).
+ * The flat element of the multiple junctions between the simple junction are stored in a group named "jointsMultiples".
  * @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
@@ -11289,6 +11409,32 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   std::set<int> emptySet;
   emptyMap.clear();
 
+  MESSAGE(".. Number of domains :"<<theElems.size());
+
+  // Check if the domains do not share an element
+  for (int idom = 0; idom < theElems.size()-1; idom++)
+    {
+//       MESSAGE("... Check of domain #" << idom);
+      const TIDSortedElemSet& domain = theElems[idom];
+      TIDSortedElemSet::const_iterator elemItr = domain.begin();
+      for (; elemItr != domain.end(); ++elemItr)
+        {
+          SMDS_MeshElement* anElem = (SMDS_MeshElement*) *elemItr;
+          int idombisdeb = idom + 1 ;
+          for (int idombis = idombisdeb; idombis < theElems.size(); idombis++) // check if the element belongs to a domain further in the list
+          {
+            const TIDSortedElemSet& domainbis = theElems[idombis];
+            if ( domainbis.count(anElem) )
+            {
+              MESSAGE(".... Domain #" << idom);
+              MESSAGE(".... Domain #" << idombis);
+              throw SALOME_Exception("The domains are not disjoint.");
+              return false ;
+            }
+          }
+        }
+    }
+
   for (int idom = 0; idom < theElems.size(); idom++)
     {
 
@@ -11297,7 +11443,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       //     and corresponding volume of this domain, for each shared face.
       //     a volume has a face shared by 2 domains if it has a neighbor which is not in his domain.
 
-      //MESSAGE("Domain " << idom);
+      MESSAGE("... Neighbors of domain #" << idom);
       const TIDSortedElemSet& domain = theElems[idom];
       TIDSortedElemSet::const_iterator elemItr = domain.begin();
       for (; elemItr != domain.end(); ++elemItr)
@@ -11317,15 +11463,25 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
               const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
               if (! domain.count(elem)) // neighbor is in another domain : face is shared
                 {
-                  DownIdType face(downIds[n], downTypes[n]);
-                  if (!faceDomains.count(face))
-                    faceDomains[face] = emptyMap; // create an empty entry for face
-                  if (!faceDomains[face].count(idom))
-                    {
-                      faceDomains[face][idom] = vtkId; // volume associated to face in this domain
-                      celldom[vtkId] = idom;
-                      //MESSAGE("       cell with a border " << vtkId << " domain " << idom);
-                    }
+                  bool ok = false ;
+                  for (int idombis = 0; idombis < theElems.size(); idombis++) // check if the neighbor belongs to another domain of the list
+                  {
+                    // MESSAGE("Domain " << idombis);
+                    const TIDSortedElemSet& domainbis = theElems[idombis];
+                    if ( domainbis.count(elem)) ok = true ; // neighbor is in a correct domain : face is kept
+                  }
+                  if ( ok ) // the characteristics of the face is stored
+                  {
+                    DownIdType face(downIds[n], downTypes[n]);
+                    if (!faceDomains.count(face))
+                      faceDomains[face] = emptyMap; // create an empty entry for face
+                    if (!faceDomains[face].count(idom))
+                      {
+                        faceDomains[face][idom] = vtkId; // volume associated to face in this domain
+                        celldom[vtkId] = idom;
+                        //MESSAGE("       cell with a border " << vtkId << " domain " << idom);
+                      }
+                  }
                 }
             }
         }
@@ -11396,6 +11552,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   std::map<int, std::vector<int> > mutipleNodes; // nodes multi domains with domain order
   std::map<int, std::vector<int> > mutipleNodesToFace; // nodes multi domains with domain order to transform in Face (junction between 3 or more 2D domains)
 
+  MESSAGE(".. Duplication of the nodes");
   for (int idomain = 0; idomain < theElems.size(); idomain++)
     {
       itface = faceDomains.begin();
@@ -11458,6 +11615,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
         }
     }
 
+  MESSAGE(".. Creation of elements");
   for (int idomain = 0; idomain < theElems.size(); idomain++)
     {
       itface = faceDomains.begin();
@@ -11587,6 +11745,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   std::map<int, std::map<long,int> > nodeQuadDomains;
   std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
 
+  MESSAGE(".. Creation of elements: simple junction");
   if (createJointElems)
     {
       int idg;
@@ -11637,6 +11796,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   //     iterate on mutipleNodesToFace
   //     iterate on edgesMultiDomains
 
+  MESSAGE(".. Creation of elements: multiple junction");
   if (createJointElems)
     {
       // --- iterate on mutipleNodesToFace
@@ -11709,6 +11869,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   faceOrEdgeDom.clear();
   feDom.clear();
 
+  MESSAGE(".. Modification of elements");
   for (int idomain = 0; idomain < theElems.size(); idomain++)
     {
       std::map<int, std::map<int, int> >::const_iterator itnod = nodeDomains.begin();