From: eap <eap@opencascade.com>
Date: Wed, 19 Jan 2011 17:21:40 +0000 (+0000)
Subject: 0021051: [CEA 432] P2 field evaluation
X-Git-Tag: V6_main_FINAL~1107
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=d5ecebfde747fe83def203153313d99b076787d6;p=tools%2Fmedcoupling.git

0021051: [CEA 432] P2 field evaluation

+  template<class MyMeshType>
+  class PointLocatorInSimplex : public PointLocatorAlgos<MyMeshType>
+  {
+    const MyMeshType& _mesh;
+  public:
+    PointLocatorInSimplex(const MyMeshType& mesh)
+      :PointLocatorAlgos<MyMeshType>(mesh),_mesh(mesh)
+    {
+    }
+
+    //================================================================================
+    /*!
+     * \brief Returns nodes composing the simplex the point x is in
+     */
+    //================================================================================
+
+    virtual std::list<int> locates(const double* x, double eps)
---

diff --git a/src/INTERP_KERNEL/PointLocatorAlgos.txx b/src/INTERP_KERNEL/PointLocatorAlgos.txx
index d7dd7432e..a8c3000c1 100644
--- a/src/INTERP_KERNEL/PointLocatorAlgos.txx
+++ b/src/INTERP_KERNEL/PointLocatorAlgos.txx
@@ -24,6 +24,7 @@
 #include "BBTree.txx"
 
 #include <list>
+#include <set>
 #include <limits>
 
 namespace INTERP_KERNEL
@@ -217,6 +218,98 @@ namespace INTERP_KERNEL
       return (min_sign!=-1 || max_sign!=1);     
     }
   };
+
+  template<class MyMeshType>
+  class PointLocatorInSimplex : public PointLocatorAlgos<MyMeshType>
+  {
+    const MyMeshType& _mesh;
+  public:
+    PointLocatorInSimplex(const MyMeshType& mesh)
+      :PointLocatorAlgos<MyMeshType>(mesh),_mesh(mesh)
+    {
+    }
+
+    //================================================================================
+    /*!
+     * \brief Returns nodes composing the simplex the point x is in
+     */
+    //================================================================================
+
+    virtual std::list<int> locates(const double* x, double eps)
+    {
+      typedef typename MyMeshType::MyConnType ConnType;
+      const NumberingPolicy numPol=MyMeshType::My_numPol;
+
+      std::list<int> simplexNodes;
+      std::list<int> candidates = PointLocatorAlgos<MyMeshType>::locates(x,eps);
+      std::list<int>::iterator eIt = candidates.begin();
+      for ( ; eIt != candidates.end(); ++eIt )
+      {
+        const int i = OTT<ConnType,numPol>::ind2C( *eIt );
+        const double* coords= _mesh.getCoordinatesPtr();
+        const ConnType* conn=_mesh.getConnectivityPtr();
+        const ConnType* conn_index= _mesh.getConnectivityIndexPtr();
+        const ConnType* conn_elem=conn+OTT<ConnType,numPol>::ind2C(conn_index[i]);
+        int conn_elem_sz=conn_index[i+1]-conn_index[i];
+        NormalizedCellType type=_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(i));
+        CellModel cell = CellModel::getCellModel(type);
+
+        if ( cell.isQuadratic() )
+          throw Exception("P2 not implemented yet");
+
+        if ( cell.isSimplex())
+        {
+          for ( int n = 0; n < conn_elem_sz; ++n )
+            simplexNodes.push_back( conn_elem[ n ]);
+        }
+        else
+        {
+          NormalizedCellType simlexType = cell.getDimension()==3 ? NORM_TETRA4 : NORM_TRI3;
+          std::vector<int> sonNodes;
+          NormalizedCellType sonType;
+          const unsigned nbSons = cell.getNumberOfSons2( conn_elem, conn_elem_sz );
+          for ( unsigned s = 0; s < nbSons; ++s )
+          {
+            sonNodes.resize( cell.getNumberOfNodesConstituentTheSon2( s, conn_elem, conn_elem_sz ));
+            cell.fillSonCellNodalConnectivity2( s, conn_elem, conn_elem_sz, &sonNodes[0], sonType );
+            std::set<int> sonNodesSet( sonNodes.begin(), sonNodes.end() );
+
+            std::set< std::set< ConnType > > checkedSonSimplex;
+            for ( unsigned sn = 0; sn < sonNodes.size(); ++sn )
+            {
+              std::vector< ConnType > simplexConn( cell.getDimension() + 1 );
+              unsigned n;
+              for ( n = 0; n < cell.getDimension()-1; ++n )
+                simplexConn[n] = sonNodes[ (sn+n) % sonNodes.size() ];
+
+              for ( unsigned n2 = 0; n2 < sonNodes.size()-cell.getDimension()+1; ++n2 )
+              {
+                simplexConn[n] = sonNodes[ (sn+n+n2) % sonNodes.size() ];
+                std::set< ConnType > sonSimplex( simplexConn.begin(), --simplexConn.end());
+                if ( checkedSonSimplex.insert( sonSimplex ).second )
+                {
+                  for ( unsigned cn = 0; cn < conn_elem_sz; ++cn )
+                    if ( !sonNodesSet.count( conn_elem[cn] ))
+                    {
+                      simplexConn.back() = conn_elem[cn];
+                      if ( isElementContainsPoint( x, simlexType, coords,
+                                                   &simplexConn[0], simplexConn.size(), eps ))
+                      {
+                        simplexNodes.insert( simplexNodes.end(),
+                                             simplexConn.begin(), simplexConn.end());
+                        return simplexNodes;
+                      }
+                    }
+                }
+              }
+            }
+          }
+        }
+      }
+      return simplexNodes;
+    }
+
+  };
 }
 
 #endif