From 9f36bbc95a1a1726f938bd0a09bfa734c36ad80e Mon Sep 17 00:00:00 2001
From: prascle <prascle>
Date: Sat, 19 Mar 2011 06:55:48 +0000
Subject: [PATCH] PR: double nodes and flat elements for ASTER calculations in
 progress

---
 src/SMDS/SMDS_Downward.hxx         |   2 +
 src/SMDS/SMDS_UnstructuredGrid.cxx |  51 +++++++++
 src/SMDS/SMDS_UnstructuredGrid.hxx |   1 +
 src/SMESH/SMESH_MeshEditor.cxx     | 167 ++++++++++++++++++++++++++++-
 src/SMESH/SMESH_MeshEditor.hxx     |   2 +
 5 files changed, 220 insertions(+), 3 deletions(-)

diff --git a/src/SMDS/SMDS_Downward.hxx b/src/SMDS/SMDS_Downward.hxx
index 173f463f4..2226c55dc 100644
--- a/src/SMDS/SMDS_Downward.hxx
+++ b/src/SMDS/SMDS_Downward.hxx
@@ -61,6 +61,7 @@ public:
   virtual const int* getUpCells(int cellId) = 0;
   virtual const unsigned char* getUpTypes(int cellId) = 0;
   virtual void getNodeIds(int cellId, std::set<int>& nodeSet) = 0;
+  virtual int getNodes(int cellId, int* nodevec) {return 0; }
   int getVtkCellId(int cellId)
   {
     return _vtkCellIds[cellId];
@@ -99,6 +100,7 @@ public:
   virtual const int* getUpCells(int cellId);
   virtual const unsigned char* getUpTypes(int cellId);
   virtual void getNodeIds(int cellId, std::set<int>& nodeSet);
+  virtual int getNodes(int cellId, int* nodevec) { return getNodeSet(cellId, nodevec); }
 protected:
   SMDS_Down1D(SMDS_UnstructuredGrid *grid, int nbDownCells);
   ~SMDS_Down1D();
diff --git a/src/SMDS/SMDS_UnstructuredGrid.cxx b/src/SMDS/SMDS_UnstructuredGrid.cxx
index 0e84f5e33..413b6e644 100644
--- a/src/SMDS/SMDS_UnstructuredGrid.cxx
+++ b/src/SMDS/SMDS_UnstructuredGrid.cxx
@@ -793,6 +793,57 @@ int SMDS_UnstructuredGrid::GetParentVolumes(int* volVtkIds, int vtkId)
   return nbvol;
 }
 
+/*! get the volumes containing a face or an edge of the downward structure
+ * The edge or face does not necessary belong to the vtkUnstructuredGrid
+ * @param volVtkIds vector of parent volume ids to fill (reserve enough space!)
+ * @param downId id in the downward structure
+ * @param downType type of cell
+ */
+int SMDS_UnstructuredGrid::GetParentVolumes(int* volVtkIds, int downId, unsigned char downType)
+{
+  int vtkType = downType;
+  int dim = SMDS_Downward::getCellDimension(vtkType);
+  int nbFaces = 0;
+  int faces[1000];
+  unsigned char cellTypes[1000];
+  int downCellId[1000];
+  if (dim == 1)
+    {
+      nbFaces = _downArray[vtkType]->getNumberOfUpCells(downId);
+      const int *upCells = _downArray[vtkType]->getUpCells(downId);
+      const unsigned char* upTypes = _downArray[vtkType]->getUpTypes(downId);
+      for (int i=0; i< nbFaces; i++)
+        {
+          faces[i] = _downArray[upTypes[i]]->getVtkCellId(upCells[i]);
+          cellTypes[i] = upTypes[i];
+          downCellId[i] = upCells[i];
+        }
+    }
+  else if (dim == 2)
+    {
+      nbFaces = 1;
+      cellTypes[0] = vtkType;
+      downCellId[0] = downId;
+    }
+
+  int nbvol =0;
+  for (int i=0; i<nbFaces; i++)
+    {
+      int vtkTypeFace = cellTypes[i];
+      int downId = downCellId[i];
+      int nv = _downArray[vtkTypeFace]->getNumberOfUpCells(downId);
+      const int *upCells = _downArray[vtkTypeFace]->getUpCells(downId);
+      const unsigned char* upTypes = _downArray[vtkTypeFace]->getUpTypes(downId);
+       for (int j=0; j<nv; j++)
+        {
+          int vtkVolId = _downArray[upTypes[j]]->getVtkCellId(upCells[j]);
+          if (vtkVolId >= 0)
+            volVtkIds[nbvol++] = vtkVolId;
+        }
+    }
+  return nbvol;
+}
+
 /*! get the node id's of a cell.
  * The cell is defined by it's downward connectivity id and type.
  * @param nodeSet set of of vtk node id's to fill.
diff --git a/src/SMDS/SMDS_UnstructuredGrid.hxx b/src/SMDS/SMDS_UnstructuredGrid.hxx
index 8f670af14..f969cab10 100644
--- a/src/SMDS/SMDS_UnstructuredGrid.hxx
+++ b/src/SMDS/SMDS_UnstructuredGrid.hxx
@@ -63,6 +63,7 @@ public:
   void BuildDownwardConnectivity(bool withEdges);
   int GetNeighbors(int* neighborsVtkIds, int* downIds, unsigned char* downTypes, int vtkId);
   int GetParentVolumes(int* volVtkIds, int vtkId);
+  int GetParentVolumes(int* volVtkIds, int downId, unsigned char downType);
   void GetNodeIds(std::set<int>& nodeSet, int downId, unsigned char downType);
   void ModifyCellNodes(int vtkVolId, std::map<int, int> localClonedNodeIds);
   int getOrderedNodesOfFace(int vtkVolId, std::vector<vtkIdType>& orderedNodes);
diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx
index 067e750a4..3fc0ca5a0 100644
--- a/src/SMESH/SMESH_MeshEditor.cxx
+++ b/src/SMESH/SMESH_MeshEditor.cxx
@@ -94,6 +94,7 @@
 #include <set>
 #include <numeric>
 #include <limits>
+#include <algorithm>
 
 #define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
 
@@ -10818,6 +10819,28 @@ bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems,
   return DoubleNodes( theElems, theNodesNot, anAffected );
 }
 
+/*!
+ *  \brief compute an oriented angle between two planes defined by four points.
+ *  The vector (p0,p1) defines the intersection of the 2 planes (p0,p1,g1) and (p0,p1,g2)
+ *  @param p0 base of the rotation axe
+ *  @param p1 extremity of the rotation axe
+ *  @param g1 belongs to the first plane
+ *  @param g2 belongs to the second plane
+ */
+double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const gp_Pnt& g1, const gp_Pnt& g2)
+{
+//  MESSAGE("    p0: " << p0.X() << " " << p0.Y() << " " << p0.Z());
+//  MESSAGE("    p1: " << p1.X() << " " << p1.Y() << " " << p1.Z());
+//  MESSAGE("    g1: " << g1.X() << " " << g1.Y() << " " << g1.Z());
+//  MESSAGE("    g2: " << g2.X() << " " << g2.Y() << " " << g2.Z());
+  gp_Vec vref(p0, p1);
+  gp_Vec v1(p0, g1);
+  gp_Vec v2(p0, g2);
+  gp_Vec n1 = vref.Crossed(v1);
+  gp_Vec n2 = vref.Crossed(v2);
+  return n2.AngleWithRef(n1, vref);
+}
+
 /*!
  * \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.
@@ -10854,6 +10877,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   cellDomains.clear();
   nodeDomains.clear();
   std::map<int,int> emptyMap;
+  std::set<int> emptySet;
   emptyMap.clear();
 
   for (int idom = 0; idom < theElems.size(); idom++)
@@ -10895,7 +10919,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
         }
     }
 
-  MESSAGE("Number of shared faces " << faceDomains.size());
+  //MESSAGE("Number of shared faces " << faceDomains.size());
   std::map<DownIdType, std::map<int, int>, DownIdCompare>::iterator itface;
 
   // --- explore the shared faces domain by domain,
@@ -10946,6 +10970,13 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   //     for each shared face, get the nodes
   //     for each node, for each domain of the face, create a clone of the node
 
+  // --- edges at the intersection of 3 or 4 domains, with the order of domains to build
+  //     junction elements of type prism or hexa. the key is the pair of nodesId (lower first)
+  //     the value is the ordered domain ids. (more than 4 domains not taken into account)
+
+  std::map<std::vector<int>, std::vector<int> > edgesMultiDomains; // nodes of edge --> ordered domains
+  std::map<int, std::vector<int> > mutipleNodes; // nodes muti domains with domain order
+
   for (int idomain = 0; idomain < theElems.size(); idomain++)
     {
       itface = faceDomains.begin();
@@ -10959,6 +10990,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
           std::set<int> oldNodes;
           oldNodes.clear();
           grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+          bool isMultipleDetected = false;
           std::set<int>::iterator itn = oldNodes.begin();
           for (; itn != oldNodes.end(); ++itn)
             {
@@ -10973,13 +11005,112 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                 {
                   int idom = itdom->first;
                   //MESSAGE("         domain " << idom);
-                  if (!nodeDomains[oldId].count(idom))
+                  if (!nodeDomains[oldId].count(idom)) // --- node to clone
                     {
+                      if (nodeDomains[oldId].size() >= 2) // a multiple node
+                        {
+                          vector<int> orderedDoms;
+                          //MESSAGE("multiple node " << oldId);
+                          isMultipleDetected =true;
+                          if (mutipleNodes.count(oldId))
+                            orderedDoms = mutipleNodes[oldId];
+                          else
+                            {
+                              map<int,int>::iterator it = nodeDomains[oldId].begin();
+                              for (; it != nodeDomains[oldId].end(); ++it)
+                                orderedDoms.push_back(it->first);
+                            }
+                          orderedDoms.push_back(idom); // TODO order ==> push_front or back
+                          //stringstream txt;
+                          //for (int i=0; i<orderedDoms.size(); i++)
+                          //  txt << orderedDoms[i] << " ";
+                          //MESSAGE("orderedDoms " << txt.str());
+                          mutipleNodes[oldId] = orderedDoms;
+                        }
                       double *coords = grid->GetPoint(oldId);
                       SMDS_MeshNode *newNode = meshDS->AddNode(coords[0], coords[1], coords[2]);
                       int newId = newNode->getVtkId();
                       nodeDomains[oldId][idom] = newId; // cloned node for other domains
-                      //MESSAGE("         newNode " << newId);
+                      //MESSAGE("   newNode " << newId << " oldNode " << oldId << " size=" <<nodeDomains[oldId].size());
+                    }
+                }
+            }
+          if (isMultipleDetected) // check if an edge of the face is shared between 3 or more domains
+            {
+              //MESSAGE("multiple Nodes detected on a shared face");
+              int downId = itface->first.cellId;
+              unsigned char cellType = itface->first.cellType;
+              int nbEdges = grid->getDownArray(cellType)->getNumberOfDownCells(downId);
+              const int *downEdgeIds = grid->getDownArray(cellType)->getDownCells(downId);
+              const unsigned char* edgeType = grid->getDownArray(cellType)->getDownTypes(downId);
+              for (int ie =0; ie < nbEdges; ie++)
+                {
+                  int nodes[3];
+                  int nbNodes = grid->getDownArray(edgeType[ie])->getNodes(downEdgeIds[ie], nodes);
+                  if (mutipleNodes.count(nodes[0]) && mutipleNodes.count(nodes[nbNodes-1]))
+                    {
+                      vector<int> vn0 = mutipleNodes[nodes[0]];
+                      vector<int> vn1 = mutipleNodes[nodes[nbNodes - 1]];
+                      sort( vn0.begin(), vn0.end() );
+                      sort( vn1.begin(), vn1.end() );
+                      if (vn0 == vn1)
+                        {
+                          //MESSAGE(" detect edgesMultiDomains " << nodes[0] << " " << nodes[nbNodes - 1]);
+                          double *coords = grid->GetPoint(nodes[0]);
+                          gp_Pnt p0(coords[0], coords[1], coords[2]);
+                          coords = grid->GetPoint(nodes[nbNodes - 1]);
+                          gp_Pnt p1(coords[0], coords[1], coords[2]);
+                          gp_Pnt gref;
+                          int vtkVolIds[1000];  // an edge can belong to a lot of volumes
+                          map<int, SMDS_VtkVolume*> domvol; // domain --> a volume with the edge
+                          map<int, double> angleDom; // oriented angles between planes defined by edge and volume centers
+                          int nbvol = grid->GetParentVolumes(vtkVolIds, downEdgeIds[ie], edgeType[ie]);
+                          for (int id=0; id < vn0.size(); id++)
+                            {
+                              int idom = vn0[id];
+                              for (int ivol=0; ivol<nbvol; ivol++)
+                                {
+                                  int smdsId = meshDS->fromVtkToSmds(vtkVolIds[ivol]);
+                                  SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId);
+                                  if (theElems[idom].count(elem))
+                                    {
+                                      SMDS_VtkVolume* svol = dynamic_cast<SMDS_VtkVolume*>(elem);
+                                      domvol[idom] = svol;
+                                      //MESSAGE("  domain " << idom << " volume " << elem->GetID());
+                                      double values[3];
+                                      vtkIdType npts = 0;
+                                      vtkIdType* pts = 0;
+                                      grid->GetCellPoints(vtkVolIds[ivol], npts, pts);
+                                      SMDS_VtkVolume::gravityCenter(grid, pts, npts, values);
+                                      if (id ==0)
+                                        {
+                                          gref.SetXYZ(gp_XYZ(values[0], values[1], values[2]));
+                                          angleDom[idom] = 0;
+                                        }
+                                      else
+                                        {
+                                          gp_Pnt g(values[0], values[1], values[2]);
+                                          angleDom[idom] = OrientedAngle(p0, p1, gref, g); // -pi<angle<+pi
+                                          //MESSAGE("  angle=" << angleDom[idom]);
+                                        }
+                                      break;
+                                    }
+                                }
+                            }
+                          map<double, int> sortedDom; // sort domains by angle
+                          for (map<int, double>::iterator ia = angleDom.begin(); ia != angleDom.end(); ++ia)
+                            sortedDom[ia->second] = ia->first;
+                          vector<int> vnodes;
+                          vector<int> vdom;
+                          for (map<double, int>::iterator ib = sortedDom.begin(); ib != sortedDom.end(); ++ib)
+                            {
+                              vdom.push_back(ib->second);
+                              //MESSAGE("  ordered domain " << ib->second << "  angle " << ib->first);
+                            }
+                          for (int ino = 0; ino < nbNodes; ino++)
+                            vnodes.push_back(nodes[ino]);
+                          edgesMultiDomains[vnodes] = vdom; // nodes vector --> ordered domains
+                        }
                     }
                 }
             }
@@ -11016,6 +11147,36 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
         }
     }
 
+  // --- create volumes on multiple domain intersection if requested
+  //     iterate on edgesMultiDomains
+
+  if (createJointElems)
+    {
+      std::map<std::vector<int>, std::vector<int> >::iterator ite = edgesMultiDomains.begin();
+      for (; ite != edgesMultiDomains.end(); ++ite)
+        {
+          vector<int> nodes = ite->first;
+          vector<int> orderDom = ite->second;
+          vector<int> orderedNodes;
+          if (nodes.size() == 2)
+            {
+              //MESSAGE(" use edgesMultiDomains " << nodes[0] << " " << nodes[1]);
+              for (int ino=0; ino < nodes.size(); ino++)
+                if (orderDom.size() == 3)
+                  for (int idom = 0; idom <orderDom.size(); idom++)
+                    orderedNodes.push_back( nodeDomains[nodes[ino]][orderDom[idom]] );
+                else
+                  for (int idom = orderDom.size()-1; idom >=0; idom--)
+                    orderedNodes.push_back( nodeDomains[nodes[ino]][orderDom[idom]] );
+              this->GetMeshDS()->AddVolumeFromVtkIds(orderedNodes);
+            }
+          else
+            {
+              // TODO quadratic nodes
+            }
+        }
+    }
+
   // --- list the explicit faces and edges of the mesh that need to be modified,
   //     i.e. faces and edges built with one or more duplicated nodes.
   //     associate these faces or edges to their corresponding domain.
diff --git a/src/SMESH/SMESH_MeshEditor.hxx b/src/SMESH/SMESH_MeshEditor.hxx
index 337ce138d..4b52b4f10 100644
--- a/src/SMESH/SMESH_MeshEditor.hxx
+++ b/src/SMESH/SMESH_MeshEditor.hxx
@@ -565,6 +565,8 @@ public:
                             const TIDSortedElemSet& theNodesNot,
                             const TopoDS_Shape&     theShape );
   
+  double OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const gp_Pnt& g1, const gp_Pnt& g2);
+
   bool DoubleNodesOnGroupBoundaries( const std::vector<TIDSortedElemSet>& theElems,
                                      bool createJointElems);
 
-- 
2.39.2