From 1530ae5c2c38df33901fbaa4e1ad6c95a39339ae Mon Sep 17 00:00:00 2001
From: eap <eap@opencascade.com>
Date: Wed, 17 Dec 2008 17:46:51 +0000
Subject: [PATCH] MEDMEM Industrialization 2008    support P0->P1 etc
 projections

---
 src/ParaMEDMEM/InterpolationMatrix.cxx | 131 +++++++++++++++++++++----
 src/ParaMEDMEM/InterpolationMatrix.hxx |   4 +-
 2 files changed, 114 insertions(+), 21 deletions(-)

diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx
index cdc3691f9..9fba6504a 100644
--- a/src/ParaMEDMEM/InterpolationMatrix.cxx
+++ b/src/ParaMEDMEM/InterpolationMatrix.cxx
@@ -9,6 +9,7 @@
 #include "Interpolation3D.txx"
 #include "MEDNormalizedUnstructuredMesh.hxx"
 #include "InterpolationOptions.hxx"
+#include "IntersectionResult.hxx"
 
 /*! \class InterpolationMatrix
 
@@ -50,6 +51,7 @@ namespace ParaMEDMEM
       _row_offsets[i]=0;
 
     _coeffs.resize(nbelems);
+    _coef_values.reserve( target_group.size() );
   }
 
   InterpolationMatrix::~InterpolationMatrix()
@@ -73,11 +75,16 @@ namespace ParaMEDMEM
                                              int iproc_distant, int* distant_elems)
   {
     if (distant_support.getMeshDimension() != _source_support.getMeshDimension() ||
-        distant_support.getMeshDimension() != _source_support.getMeshDimension() )
+        distant_support.getSpaceDimension() != _source_support.getSpaceDimension() )
       throw MEDMEM::MEDEXCEPTION("local and distant meshes do not have the same space and mesh dimensions");
 
+    bool isP1Method = ( getMethod() == "P1" || getMethod() == "P1d" );
+
     //creating the interpolator structure
-    vector<map<int,double> > surfaces;
+    //vector<map<int,double> > surfaces;
+    INTERP_KERNEL::IntersectionResult surfaces( _source_support.getSpaceDimension() );
+    surfaces.setNeedBaryCentres( isP1Method );
+
     //computation of the intersection volumes between source and target elements
     int source_size = _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
 
@@ -94,6 +101,10 @@ namespace ParaMEDMEM
       MEDNormalizedUnstructuredMesh<2,2> source_wrapper(&_source_support);
 
       INTERP_KERNEL::Interpolation2D interpolator (*this);
+      if ( distant_support.getNumberOfPolygons(MED_CELL) > 0 ||
+           _source_support.getNumberOfPolygons(MED_CELL) > 0 )
+        // it must be DualMESH
+        interpolator.setIntersectionType( INTERP_KERNEL::Geometric2D );
       interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces);
     }
     else if (distant_support.getMeshDimension()==3 && distant_support.getSpaceDimension()==3)
@@ -122,6 +133,22 @@ namespace ParaMEDMEM
                                     "all cells", MED_EN::MED_CELL);
     MEDMEM::FIELD<double>* source_triangle_surf = getSupportVolumes(source_support);
 
+    // allocate vector of coefficients
+    int nbNodesPerCell = source_support.getTypes()[0] % 100;
+    int nbValuesPerCell = isP1Method ? nbNodesPerCell : 1;
+    int nbIntersections = surfaces.getNumberOfIntersections();
+    _coef_values.push_back(vector<double>());
+    vector<double>& coef_values = _coef_values.back();
+    coef_values.reserve( nbIntersections * nbValuesPerCell );
+
+    // prepare to getting coordinates of nodes of a source cell
+    // needed for comuting barycentric coordinates of intersection gravity center
+    const int* src_conn = _source_support.getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,
+                                                          MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+    const double * src_coords = _source_support.getCoordinates(MED_EN::MED_FULL_INTERLACE);
+    vector<const double*> cell_coords( nbNodesPerCell );
+    vector<double> bary_coords( nbNodesPerCell );
+
     //storing the source volumes
     _source_volume.resize(source_size);
     for (int i=0; i<source_size; i++)
@@ -165,13 +192,34 @@ namespace ParaMEDMEM
         //col_id is the number of the column
         //iter->second is the value of the coefficient
 
-        _coeffs[ielem].push_back(make_pair(col_id,iter->second));
+        //_coeffs[ielem].push_back(make_pair(col_id,iter->second));
+        if ( nbValuesPerCell == 1 )
+        {
+          coef_values.push_back(iter->second);
+          _coeffs[ielem].push_back(make_pair(col_id, &coef_values.back() ));
+        }
+        else
+        {
+          // barycenter of intersection of cells ielem and iter->first
+          double* bary_centre = surfaces.getBaryCentre( ielem, iter->first );
+          // coordinates of a source cell
+          const int* cell_conn = src_conn + nbNodesPerCell * ielem;
+          for ( int i = 0; i < nbNodesPerCell; ++i, ++cell_conn )
+            cell_coords[i] = src_coords + _source_support.getSpaceDimension() * ( *cell_conn-1 );
+          // barycentric coodrs
+          INTERP_KERNEL::barycentric_coodrs( cell_coords, bary_centre, &bary_coords[0]);
+          // store (surface * bary_coords)
+          coef_values.push_back(iter->second * bary_coords[0]);
+          _coeffs[ielem].push_back(make_pair(col_id, &coef_values.back() ));
+          for ( int j = 1; j < bary_coords.size(); ++j )
+            coef_values.push_back(iter->second * bary_coords[j]);
+        }
       }
     }
     delete source_triangle_surf;
     delete target_triangle_surf;
   }
-	
+
   /*! The call to this method updates the arrays on the target side
     so that they know which amount of data from which processor they 
     should expect. 
@@ -212,18 +260,57 @@ namespace ParaMEDMEM
       // W is the intersection matrix
       // S is the source vector
 
-      for (int irow=0; irow<nbrows; irow++)
-        for (int icomp=0; icomp< nbcomp; icomp++)
-        {
-          double coeff_row = field.getValueIJ(irow+1,icomp+1);
-          for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
+      if ( getMethod() == "P0" ) { // --------------- P0 - one value per cell
+
+        for (int irow=0; irow<nbrows; irow++)
+          for (int icomp=0; icomp< nbcomp; icomp++)
           {
-            int colid= _coeffs[irow][icol-_row_offsets[irow]].first;
-            double value = _coeffs[irow][icol-_row_offsets[irow]].second;
-            target_value[(colid-1)*nbcomp+icomp]+=value*coeff_row;
+            double coeff_row = field.getValueIJ(irow+1,icomp+1);
+            for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
+            {
+              int     colid = _coeffs[irow][icol-_row_offsets[irow]].first;
+              double* value = _coeffs[irow][icol-_row_offsets[irow]].second;
+              target_value[(colid-1)*nbcomp+icomp]+=value[0]*coeff_row;
+            }
           }
-        }
-
+      }
+      else if ( getMethod() == "P1" ) { // ---------- P1 - one value per node
+
+        int nbValuesPerCell  = _source_support.getTypes(MED_EN::MED_CELL)[0] % 100;
+        const int* cell_conn = _source_support.getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,
+                                                               MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+
+        for (int irow=0; irow<nbrows; irow++, cell_conn += nbValuesPerCell)
+          for (int icomp=0; icomp< nbcomp; icomp++)
+            for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
+            {
+              int     colid = _coeffs[irow][icol-_row_offsets[irow]].first;
+              double* value = _coeffs[irow][icol-_row_offsets[irow]].second;
+              for ( int ival=0; ival<nbValuesPerCell; ++ival)
+              {
+                double coeff_row = field.getValueIJ(cell_conn[ival],icomp+1);
+                target_value[(colid-1)*nbcomp+icomp]+=value[ival]*coeff_row;
+              }
+            }
+      }
+      else { // ------------------------------------- P1d - different values per node
+
+        int nbValuesPerCell  = _source_support.getTypes(MED_EN::MED_CELL)[0] % 100;
+
+        for (int irow=0; irow<nbrows; irow++)
+          for (int icomp=0; icomp< nbcomp; icomp++)
+            for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
+            {
+              int     colid = _coeffs[irow][icol-_row_offsets[irow]].first;
+              double* value = _coeffs[irow][icol-_row_offsets[irow]].second;
+              for ( int ival=0; ival<nbValuesPerCell; ++ival)
+              {
+                double coeff_row = field.getValueIJK(irow+1,icomp+1,ival+1);
+                target_value[(colid-1)*nbcomp+icomp]+=value[ival]*coeff_row;
+              }
+            }
+      }
+      
       // performing VT^(-1).(W.S)
       // where VT^(-1) is the inverse of the diagonal matrix containing 
       // the volumes of target cells
@@ -276,12 +363,12 @@ namespace ParaMEDMEM
         for (int icol = _row_offsets[irow]; icol < _row_offsets[irow+1]; icol++)
         {
           int colid = _coeffs[irow][icol-_row_offsets[irow]].first;
-          double value = _coeffs[irow][icol-_row_offsets[irow]].second;
+          double* value = _coeffs[irow][icol-_row_offsets[irow]].second;
           
           for (int icomp=0; icomp<nbcomp; icomp++)
           {
             double coeff_row = source_value[(colid-1)*nbcomp+icomp];
-            target_value[irow*nbcomp+icomp] += value*coeff_row;
+            target_value[irow*nbcomp+icomp] += value[0]*coeff_row;
           }
         }
 
@@ -292,10 +379,14 @@ namespace ParaMEDMEM
         for (int icomp=0; icomp<nbcomp; icomp++)
           target_value[i*nbcomp+icomp] /= _source_volume[i];
 
-      //storing S= VS^(-1).(WT.T)
-      for (int irow=0; irow<nbrows; irow++)
-        for (int icomp=0; icomp<nbcomp; icomp++)
-          field.setValueIJ(irow+1,icomp+1,target_value[irow*nbcomp+icomp]);
+      if ( getMethod() == "P0" ) {
+
+        //storing S= VS^(-1).(WT.T)
+        for (int irow=0; irow<nbrows; irow++)
+          for (int icomp=0; icomp<nbcomp; icomp++)
+            field.setValueIJ(irow+1,icomp+1,target_value[irow*nbcomp+icomp]);
+
+      }
     }
   }
 
diff --git a/src/ParaMEDMEM/InterpolationMatrix.hxx b/src/ParaMEDMEM/InterpolationMatrix.hxx
index 09cd332f5..9faaef41e 100644
--- a/src/ParaMEDMEM/InterpolationMatrix.hxx
+++ b/src/ParaMEDMEM/InterpolationMatrix.hxx
@@ -42,7 +42,9 @@ namespace ParaMEDMEM
     const ProcessorGroup& _target_group;
     vector<double> _target_volume;
     vector<double> _source_volume;
-    vector<vector<pair<int,double> > > _coeffs;
+//     vector<vector<pair<int,double> > > _coeffs;
+    vector<vector<pair<int,double*> > > _coeffs;
+    vector< vector<double> > _coef_values; // values per target mesh
   };
   
 }
-- 
2.39.2