MEDMEM Industrialization 2008

  new version

diff --git a/src/INTERP_KERNEL/Interpolation3D.txx b/src/INTERP_KERNEL/Interpolation3D.txx
index b0fb83433ce3e235b5dffe94f91d2ee9755d012d..214d1a60959dd4d4271cd2fa85b1721f9b9ecae6 100644 (file)
--- a/src/INTERP_KERNEL/Interpolation3D.txx
+++ b/src/INTERP_KERNEL/Interpolation3D.txx
@@ -80,7 +80,10 @@ namespace INTERP_KERNEL
    *
    */
   template<class MatrixType, class MyMeshType>
-  void Interpolation3D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result)
+  void Interpolation3D::interpolateMeshes(const MyMeshType& srcMesh,
+                                          const MyMeshType& targetMesh,
+                                          MatrixType& result,
+                                          const std::string& method)
   {
     typedef typename MyMeshType::MyConnType ConnType;
     const NumberingPolicy numPol=MyMeshType::My_numPol;
@@ -107,9 +110,22 @@ namespace INTERP_KERNEL
         targetElems[i] = new MeshElement<ConnType>(OTT<ConnType,numPol>::indFC(i), targetMesh);
       }
     
-    // create empty maps for all source elements
-    result.resize(numTargetElems);
+    bool srcCornerCoeff = ( method[1] == '1' ); // "P1P0" or "P1dP0"
+    bool tgtCornerCoeff = ( method[3] == '1' ); // "P0P1" or "P0P1d"
+    bool nodeCoeff      = ( method.size() < 5 ); // no P1d
+    bool cornerCoeff = ( srcCornerCoeff || tgtCornerCoeff );
+    const int nbCorners = 4;  // coeffs for corners are calculated for tetras only
+    double baryCentreBuf[3], baryCoordBuf[nbCorners];
+    double* baryCentre = cornerCoeff ? baryCentreBuf : 0;
+    double* baryCoords = cornerCoeff ? baryCoordBuf : 0;
 
+    // create empty maps for all source elements
+    if ( tgtCornerCoeff && nodeCoeff)
+      result.resize(targetMesh.getNumberOfNodes());
+    else if ( tgtCornerCoeff && !nodeCoeff)
+      result.resize(numTargetElems * nbCorners);
+    else
+      result.resize(numTargetElems);
 
 #ifdef USE_RECURSIVE_BBOX_FILTER
     
@@ -196,20 +212,56 @@ namespace INTERP_KERNEL
                 iter != currNode->getSrcRegion().getEndElements() ; ++iter)
               {
                 const ConnType srcIdx = (*iter)->getIndex();
-                double* baryCenre =
-                  result.getBaryCentre(OTT<ConnType,numPol>::ind2C(targetIdx),srcIdx);
-                const double vol = intersector->intersectSourceCell(srcIdx,baryCenre);
+                const double vol = intersector->intersectSourceCell(srcIdx,baryCentre);
 
                 if(vol != 0.0)
                   {
-                    result[OTT<ConnType,numPol>::ind2C(targetIdx)].insert(make_pair(srcIdx,vol));
-                    LOG3(2, "Result : V (" << srcIdx << "- " << targetIdx << ") = " <<  result[ OTT<ConnType,numPol>::ind2C(srcIdx) ][targetIdx] );
+                    if ( !baryCentre )
+                    {
+                      result[OTT<ConnType,numPol>::ind2C(targetIdx)].insert(make_pair(srcIdx,vol));
+                      LOG3(2, "Result : V (" << srcIdx << "- " << targetIdx << ") = " <<  result[ OTT<ConnType,numPol>::ind2C(srcIdx) ][targetIdx] );
+                    }
+                    else if ( tgtCornerCoeff )
+                    {
+                      // this code intended for tetras only -> 4 nodes
+                      getBarycentricCoordinates<MyMeshType, 4>( baryCentre, targetIdx,
+                                                                targetMesh, baryCoords );
+                      if ( nodeCoeff ) {  // 1 coeff per target node
+                        for ( int n = 0; n < 4; ++n ) {
+                          ConnType node = getGlobalNumberOfNode(n, targetIdx, targetMesh);
+                          result[node].insert(make_pair(srcIdx,vol*baryCoords[n]));
+                        }
+                      }
+                      else { // several coeffs per target node
+                        int tgtIdx0 = nbCorners * OTT<ConnType,numPol>::ind2C( targetIdx );
+                        result[tgtIdx0+0].insert(make_pair(srcIdx,vol*baryCoords[0]));
+                        result[tgtIdx0+1].insert(make_pair(srcIdx,vol*baryCoords[1]));
+                        result[tgtIdx0+2].insert(make_pair(srcIdx,vol*baryCoords[2]));
+                        result[tgtIdx0+3].insert(make_pair(srcIdx,vol*baryCoords[3]));
+                      }
+                    }
+                    else
+                    {
+                      getBarycentricCoordinates<MyMeshType,4>( baryCentre, srcIdx,
+                                                               srcMesh, baryCoords );
+                      if ( nodeCoeff ) {  // 1 coeff per source node
+                        for ( int n = 0; n < 4; ++n ) {
+                          ConnType node = getGlobalNumberOfNode(n, srcIdx, srcMesh);
+                          result[OTT<ConnType,numPol>::ind2C(targetIdx)].insert(make_pair(OTT<ConnType,numPol>::indFC(node),vol*baryCoords[n]));
+                        }
+                      }
+                      else { // several coeffs per source node
+                        int srcIdx0 = OTT<ConnType,numPol>::indFC( nbCorners * OTT<ConnType,numPol>::ind2C(srcIdx) );
+                        for ( int n = 0; n < 4; ++n )
+                          result[OTT<ConnType,numPol>::ind2C(targetIdx)].insert(std::make_pair(srcIdx0+n,vol*baryCoords[0]));
+                      }
+                    }
                   }
               }
            
             delete intersector;
 
-          } 
+          }
         else // recursion 
           {
 
@@ -334,32 +386,67 @@ namespace INTERP_KERNEL
         IntersectorTetra intersector(srcMesh, targetMesh, targetIdx);
 
         for(vector<int>::const_iterator iter = intersectElems.begin() ; iter != intersectElems.end() ; ++iter)
-          {
+        {
 
-            const int srcIdx = *iter + 1;
-            double baryCenre* = result.getBaryCentre(targetIdx-1,srcIdx);
-            const double vol = intersector.intersectSourceCell(srcIdx, baryCenre);
+          const int srcIdx = *iter + 1;
+          const double vol = intersector.intersectSourceCell(srcIdx, baryCenre);
 
-            if(vol != 0.0)
-              {
-                result[targetIdx - 1].insert(make_pair(srcIdx, vol));
+          if(vol != 0.0)
+          {
+            if ( !baryCentre )
+            {
+              result[targetIdx - 1].insert(make_pair(srcIdx, vol));
+            }
+            else if ( tgtCornerCoeff )
+            {
+              // this code intended for tetras only -> 4 nodes
+              getBarycentricCoordinates<MyMeshType,4>( baryCentre, targetIdx, targetMesh, baryCoords );
+              if ( nodeCoeff ) {  // 1 coeff per target node
+                for ( int n = 0; n < 4; ++n ) {
+                  ConnType node = getGlobalNumberOfNode(n, targetIdx, myMesh_P);
+                  result[node].insert(make_pair(srcIdx,vol*baryCoords[n]));
+                }
               }
-
+              else { // several coeffs per target node
+                int tgtIdx0 = nbCorners * OTT<ConnType,numPol>::ind2C( targetIdx );
+                result[tgtIdx0+0].insert(make_pair(srcIdx,vol*baryCoords[0]));
+                result[tgtIdx0+1].insert(make_pair(srcIdx,vol*baryCoords[1]));
+                result[tgtIdx0+2].insert(make_pair(srcIdx,vol*baryCoords[2]));
+                result[tgtIdx0+3].insert(make_pair(srcIdx,vol*baryCoords[3]));
+              }
+            }
+            else
+            {
+              getBarycentricCoordinates<MyMeshType,4>( baryCentre, srcIdx, srcMesh, baryCoords );
+              if ( nodeCoeff ) {  // 1 coeff per source node
+                for ( n = 0; n < 4; ++n ) {
+                  ConnType node = getGlobalNumberOfNode(n, srcIdx, myMesh_S);
+                  result[OTT<ConnType,numPol>::ind2C(targetIdx)].insert(make_pair(OTT<ConnType,numPol>::indFC(node),vol*baryCoords[n]));
+                }
+              }
+              else { // several coeffs per source node
+                int srcIdx0 = OTT<ConnType,numPol>::indFC( nbCorners * OTT<ConnType,numPol>::ind2C(srcIdx) );
+                result[OTT<ConnType,numPol>::ind2C(targetIdx)].insert(make_pair(srcIdx0+0,vol*baryCoords[0]));
+                result[OTT<ConnType,numPol>::ind2C(targetIdx)].insert(make_pair(srcIdx0+1,vol*baryCoords[1]));
+                result[OTT<ConnType,numPol>::ind2C(targetIdx)].insert(make_pair(srcIdx0+2,vol*baryCoords[2]));
+                result[OTT<ConnType,numPol>::ind2C(targetIdx)].insert(make_pair(srcIdx0+3,vol*baryCoords[3]));
+              }
+            }
           }
+        }
       }
-    
 #endif
 
 
     // free allocated memory
     for(unsigned long i = 0 ; i < numSrcElems ; ++i)
-      {
-        delete srcElems[i];
-      }
+    {
+      delete srcElems[i];
+    }
     for(unsigned long i = 0 ; i < numTargetElems ; ++i)
-      {
-        delete targetElems[i];
-      }
+    {
+      delete targetElems[i];
+    }
 
 
   }

diff --git a/src/INTERP_KERNEL/InterpolationPlanar.txx b/src/INTERP_KERNEL/InterpolationPlanar.txx
index e0559be37c8b6316731922a84c81d389296267f4..f9061110f922b23921033a3982e0749bd1f82fb4 100644 (file)
--- a/src/INTERP_KERNEL/InterpolationPlanar.txx
+++ b/src/INTERP_KERNEL/InterpolationPlanar.txx
@@ -32,6 +32,7 @@
 #include "VectorUtils.hxx"
 #include "BBTree.txx"
 #include<time.h>
+#include "MeshUtils.hxx"
 
 namespace INTERP_KERNEL
 {
@@ -104,7 +105,10 @@ namespace INTERP_KERNEL
       */
   template<class RealPlanar>
   template<class MatrixType, class MyMeshType>
-  void InterpolationPlanar<RealPlanar>::interpolateMeshes(const MyMeshType& myMesh_S, const MyMeshType& myMesh_P, MatrixType& result)
+  void InterpolationPlanar<RealPlanar>::interpolateMeshes(const MyMeshType&  myMesh_S,
+                                                          const MyMeshType&  myMesh_P,
+                                                          MatrixType&        result,
+                                                          const std::string& method)
   {
     static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
     typedef typename MyMeshType::MyConnType ConnType;
@@ -184,7 +188,21 @@ namespace INTERP_KERNEL
 
     long end_filtering=clock();
 
-    result.resize(nbMaille_P);//on initialise.
+    bool srcCornerCoeff = ( method[1] == '1' ); // "P1P0" or "P1dP0"
+    bool tgtCornerCoeff = ( method[3] == '1' ); // "P0P1" or "P0P1d"
+    bool nodeCoeff      = ( method.size() < 5 ); // no P1d
+    bool cornerCoeff = ( srcCornerCoeff || tgtCornerCoeff );
+    const int nbCorners = 3;  // coeffs for corners are calculated for triangles only
+    double baryCentreBuf[SPACEDIM], baryCoordBuf[nbCorners];
+    double* baryCentre = cornerCoeff ? baryCentreBuf : 0;
+    double* baryCoords = cornerCoeff ? baryCoordBuf : 0;
+
+    if ( tgtCornerCoeff && nodeCoeff)
+      result.resize(myMesh_P.getNumberOfNodes());
+    else if ( tgtCornerCoeff && !nodeCoeff)
+      result.resize(nbMaille_P * nbCorners);
+    else
+      result.resize(nbMaille_P);//on initialise.
 
     /****************************************************/
     /* Loop on the target cells - core of the algorithm */
@@ -210,10 +228,9 @@ namespace INTERP_KERNEL
             //BBTree structure returns numbers between 0 and n-1
             int i_S=intersecting_elems[ielem]; //MN: Global number of cell ?
             int nb_nodesS=connIndxS[i_S+1]-connIndxS[i_S];
-            double* baryCenre = result.getBaryCentre( i_P, OTT<ConnType,numPol>::indFC(i_S) );
             double surf=intersector->intersectCells(OTT<ConnType,numPol>::indFC(i_P),
                                                     OTT<ConnType,numPol>::indFC(i_S),
-                                                    nb_nodesP,nb_nodesS,baryCenre);
+                                                    nb_nodesP,nb_nodesS,baryCentre);
 
             //filtering out zero surfaces and badly oriented surfaces
             // orientation = -1,0,1
@@ -221,20 +238,66 @@ namespace INTERP_KERNEL
             // 0 : the intersection is always taken into account
             // 1 : the intersection is taken into account if target and cells have the same orientation
             int orientation=InterpolationOptions::getOrientation();
-            if (( surf > 0.0 && orientation >=0 ) || ( surf < 0.0 && orientation <=0 ))
-              result[i_P].insert(std::make_pair(OTT<ConnType,numPol>::indFC(i_S),surf));
+            if (( surf > 0.0 && orientation >=0 ) || ( surf < 0.0 && orientation <=0 )) {
+              if ( !baryCentre )
+              {
+                result[i_P].insert(std::make_pair(OTT<ConnType,numPol>::indFC(i_S),surf));
+              }
+              else if ( tgtCornerCoeff )
+              {
+                // this code intended for triangles only -> 3 nodes
+                getBarycentricCoordinates<MyMeshType,3>( baryCentre, OTT<ConnType,numPol>::indFC(i_P),
+                                                         myMesh_P, baryCoords );
+                if ( nodeCoeff ) { // 1 coeff per target node
+                  cout << "P elem " << OTT<ConnType,numPol>::indFC(i_P)
+                       << ", S elem " << OTT<ConnType,numPol>::indFC(i_S)
+                       << ", surf: " << surf
+                       << endl;
+                  for ( int n = 0; n < 3; ++n ) {
+                    ConnType node = getGlobalNumberOfNode(n, OTT<ConnType,numPol>::indFC(i_P), myMesh_P);
+                    result[node].insert(make_pair(OTT<ConnType,numPol>::indFC(i_S),surf*baryCoords[n]));
+                    cout << "node " << node << " baryCoord[" << n << "]=" << baryCoords[n]
+                         << ", surf*bc = " << surf*baryCoords[n] << endl;
+                  }
+                }
+                else { // several coeffs per target node
+                  int tgtIdx0 = nbCorners * i_P;
+                  result[tgtIdx0+0].insert(make_pair(OTT<ConnType,numPol>::indFC(i_S),surf*baryCoords[0]));
+                  result[tgtIdx0+1].insert(make_pair(OTT<ConnType,numPol>::indFC(i_S),surf*baryCoords[1]));
+                  result[tgtIdx0+2].insert(make_pair(OTT<ConnType,numPol>::indFC(i_S),surf*baryCoords[2]));
+                }
+              }
+              else
+              {
+                getBarycentricCoordinates<MyMeshType,3>( baryCentre, OTT<ConnType,numPol>::indFC(i_S),
+                                                         myMesh_S, baryCoords );
+                if ( nodeCoeff ) { // 1 coeff per source node
+                  for ( int n = 0; n < 3; ++n ) {
+                    ConnType node = getGlobalNumberOfNode(n, OTT<ConnType,numPol>::indFC(i_S), myMesh_S);
+                    result[i_P].insert(make_pair(OTT<ConnType,numPol>::indFC(node),surf*baryCoords[n]));
+                  }
+                }
+                else { // several coeffs per source node
+                  int srcIdx0 = OTT<ConnType,numPol>::indFC(nbCorners * i_S);
+                  result[i_P].insert(make_pair(srcIdx0+0,surf*baryCoords[0]));
+                  result[i_P].insert(make_pair(srcIdx0+1,surf*baryCoords[1]));
+                  result[i_P].insert(make_pair(srcIdx0+2,surf*baryCoords[2]));
+                }
+              }
+            }
             counter++;
           }
         intersecting_elems.clear();
       }
     delete intersector;
 
+
     /***********************************/
     /*        DEBUG prints             */
     /***********************************/
 
     if (InterpolationOptions::getPrintLevel() >=1)
-      {
+    {
         long end_intersection=clock();
 //         if (_printLevel >=2)
 //           {