From: eap <eap@opencascade.com>
Date: Wed, 14 Jan 2009 10:34:53 +0000 (+0000)
Subject: MEDMEM Industrialization 2008
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=5316e403657524ad7ed3fd28fb6e94ffd6731cda;p=tools%2Fmedcoupling.git

MEDMEM Industrialization 2008
  Form interpolation matrix in intersectors
---

diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx
index 43ab7e4cb..435783751 100644
--- a/src/ParaMEDMEM/InterpolationMatrix.cxx
+++ b/src/ParaMEDMEM/InterpolationMatrix.cxx
@@ -46,12 +46,18 @@ namespace ParaMEDMEM
     INTERP_KERNEL::InterpolationOptions(interp_options)
   {
     int nbelems = _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+    if ( getMethod() == "P1" )
+      // a row per node
+      nbelems = _source_support.getNumberOfNodes();
+    if ( getMethod() == "P1d" )
+      // a row per each node of cell
+      nbelems *= ( _source_support.getTypes(MED_EN::MED_CELL)[0] % 100 );
+
     _row_offsets.resize(nbelems+1);
     for (int i=0; i<nbelems+1; i++)
       _row_offsets[i]=0;
 
     _coeffs.resize(nbelems);
-    _coef_values.reserve( target_group.size() );
   }
 
   InterpolationMatrix::~InterpolationMatrix()
@@ -74,16 +80,17 @@ namespace ParaMEDMEM
   void InterpolationMatrix::addContribution (MEDMEM::MESH& distant_support,
                                              int iproc_distant, int* distant_elems)
   {
+    // sleep( _source_group.myRank() );// useful for debug outputting
     if (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");
+      throw MEDMEM::MEDEXCEPTION("local and distant meshes do not have the same dimensions");
 
-    bool isP1Method = ( getMethod() == "P1" || getMethod() == "P1d" );
+    string method = "P0" + getMethod();
 
-    //creating the interpolator structure
-    //vector<map<int,double> > surfaces;
-    INTERP_KERNEL::IntersectionResult surfaces( _source_support.getSpaceDimension() );
-    surfaces.setNeedBaryCentres( isP1Method );
+    //creating the interpolator structure. Use multimap since at "P1" projection, there are
+    // <int,double> pairs with equal int.
+    typedef multimap< int,double> TRow ;
+    vector< TRow > surfaces;
 
     //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);
@@ -93,7 +100,7 @@ namespace ParaMEDMEM
       MEDNormalizedUnstructuredMesh<3,2> target_wrapper (&distant_support);
       MEDNormalizedUnstructuredMesh<3,2> source_wrapper (&_source_support);
       INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
-      interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces);
+      interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,method);
     }
     else if (distant_support.getMeshDimension()==2 && distant_support.getSpaceDimension()==2)
     {
@@ -105,7 +112,7 @@ namespace ParaMEDMEM
            _source_support.getNumberOfPolygons(MED_CELL) > 0 )
         // it must be DualMESH
         interpolator.setIntersectionType( INTERP_KERNEL::Geometric2D );
-      interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces);
+      interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,method);
     }
     else if (distant_support.getMeshDimension()==3 && distant_support.getSpaceDimension()==3)
     {
@@ -113,16 +120,16 @@ namespace ParaMEDMEM
       MEDNormalizedUnstructuredMesh<3,3> source_wrapper (&_source_support);
 
       INTERP_KERNEL::Interpolation3D interpolator (*this);
-      interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces);
+      interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,method);
     }
     else
     {
       throw MEDMEM::MEDEXCEPTION("no interpolator exists for these mesh and space dimensions ");
     }
 
-    if (surfaces.size() != source_size)
+    if (surfaces.size() != _coeffs.size())
     {
-      cout<<"surfaces.size()="<<surfaces.size()<<" source_size="<<source_size<<endl;  
+      cout<<"surfaces.size()="<<surfaces.size()<<" _coeffs.size()="<<_coeffs.size()<<endl;
       throw MEDEXCEPTION("uncoherent number of rows in interpolation matrix");
     }
 
@@ -133,22 +140,6 @@ 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++)
@@ -160,7 +151,7 @@ namespace ParaMEDMEM
     {
       _row_offsets[ielem+1] += surfaces[ielem].size();
       //_source_indices.push_back(make_pair(iproc_distant, ielem));
-      for (map<int,double>::const_iterator iter = surfaces[ielem].begin();
+      for (TRow::const_iterator iter = surfaces[ielem].begin();
            iter != surfaces[ielem].end();
            iter++)
       {
@@ -192,35 +183,12 @@ 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));
-        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.findBaryCentre( 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]);
-
-//           cout << _source_group.myRank() << "| int area: "<< iter->second << endl;
-//           cout << "bary_centre: " << bary_centre[0] << ", " << bary_centre[1] << ", "
-//                << bary_centre[0] << ") " << endl << "bary_coords: ";
-//           for ( int j = 0; j < bary_coords.size(); ++j )
-//             cout << bary_coords[j] << " ";
-//           cout << endl;
-        }
+        _coeffs[ielem].push_back(make_pair(col_id,iter->second));
+
+//         cout << _source_group.myRank() << "-" << iproc_distant
+//              << "| ielem: "<< ielem
+//              << " | tgt elem: " << iter->first << ", col=" << col_id
+//              << ", vol = " << iter->second << endl;
       }
     }
     delete source_triangle_surf;
@@ -235,7 +203,8 @@ namespace ParaMEDMEM
   */
   void InterpolationMatrix::prepare()
   {
-    int nbelems = _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+    //int nbelems = _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+    int nbelems = _coeffs.size();
     for (int ielem=0; ielem < nbelems; ielem++)
     {
       _row_offsets[ielem+1]+=_row_offsets[ielem];
@@ -275,54 +244,55 @@ namespace ParaMEDMEM
             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;
+              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;
             }
           }
       }
       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)
+//         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);
+        nbrows = _coeffs.size();
+        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++)
             {
-              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;
-//                 cout << _source_group.myRank() << ": row " << irow << ", colid: " << colid
-//                      << " -- coeff_row: " << coeff_row
-//                      << ", value: " << value[ival] << endl;
-              }
+              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;
+//               cout << _source_group.myRank() << ": row " << irow << ", colid: " << colid
+//                    << ", value: " << value
+//                    << ", target_value[" << colid << "]=" << target_value[(colid-1)*nbcomp+icomp]
+//                    << endl;
             }
+          }
       }
       else { // ------------------------------------- P1d - different values per node
 
         int nbValuesPerCell  = _source_support.getTypes(MED_EN::MED_CELL)[0] % 100;
 
+        nbrows = _coeffs.size();
         for (int irow=0; irow<nbrows; irow++)
           for (int icomp=0; icomp< nbcomp; icomp++)
+          {
+            int i = irow/nbValuesPerCell, k = irow%nbValuesPerCell;
+            double coeff_row = field.getValueIJK( i+1, icomp+1, k+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;
-              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;
+              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;
 //                 cout << _source_group.myRank() << ": row " << irow << ", colid: " << colid
 //                      << ", ival: " << ival
 //                      << " -- coeff_row: " << coeff_row
-//                      << ", value: " << value[ival] << endl;
-              }
+//                      << ", value: " << value << endl;
             }
+          }
       }
       
       // performing VT^(-1).(W.S)
@@ -366,34 +336,34 @@ namespace ParaMEDMEM
     //treatment of the transpose matrix multiply on the source side
     if (_source_group.containsMyRank())
     {
-      int nbrows = _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+      if ( getMethod() == "P0" ) {
 
-      vector<double> target_value(nbrows*nbcomp,0.0);
+        int nbrows = _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
 
-      //performing WT.T
-      //WT is W transpose
-      //T is the target vector
-      for (int irow = 0; irow < nbrows; irow++)
-        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 icomp=0; icomp<nbcomp; icomp++)
+        vector<double> target_value(nbrows*nbcomp,0.0);
+
+        //performing WT.T
+        //WT is W transpose
+        //T is the target vector
+        for (int irow = 0; irow < nbrows; irow++)
+          for (int icol = _row_offsets[irow]; icol < _row_offsets[irow+1]; icol++)
           {
-            double coeff_row = source_value[(colid-1)*nbcomp+icomp];
-            target_value[irow*nbcomp+icomp] += value[0]*coeff_row;
-          }
-        }
+            int    colid = _coeffs[irow][icol-_row_offsets[irow]].first;
+            double value = _coeffs[irow][icol-_row_offsets[irow]].second;
 
-      //performing VS^(-1).(WT.T)
-      //VS^(-1) is the inverse of the diagonal matrix storing
-      //volumes of the source cells
-      for (int i=0; i<nbrows; i++)
-        for (int icomp=0; icomp<nbcomp; icomp++)
-          target_value[i*nbcomp+icomp] /= _source_volume[i];
+            for (int icomp=0; icomp<nbcomp; icomp++)
+            {
+              double coeff_row = source_value[(colid-1)*nbcomp+icomp];
+              target_value[irow*nbcomp+icomp] += value*coeff_row;
+            }
+          }
 
-      if ( getMethod() == "P0" ) {
+        //performing VS^(-1).(WT.T)
+        //VS^(-1) is the inverse of the diagonal matrix storing
+        //volumes of the source cells
+        for (int i=0; i<nbrows; i++)
+          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++)
diff --git a/src/ParaMEDMEM/InterpolationMatrix.hxx b/src/ParaMEDMEM/InterpolationMatrix.hxx
index 9faaef41e..904f17c7c 100644
--- a/src/ParaMEDMEM/InterpolationMatrix.hxx
+++ b/src/ParaMEDMEM/InterpolationMatrix.hxx
@@ -42,9 +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<double> > _coef_values; // values per target mesh
+    vector<vector<pair<int,double> > > _coeffs;
+    //    vector<vector<pair<int,double*> > > _coeffs;
+    //    vector< vector<double> > _coef_values; // values per target mesh
   };
   
 }
diff --git a/src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx b/src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx
index fd4cf270e..6f3345e28 100644
--- a/src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx
+++ b/src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx
@@ -125,7 +125,7 @@ void ParaMEDMEMTest::testIntersectionDEC_2D_(const string& srcMethod/*="P0"*/,
   if (tmp_dir == "")
     tmp_dir = "/tmp";
   string filename_xml1              = data_dir + "/share/salome/resources/med/square1_split";
-  string filename_xml2              = data_dir + "/share/salome/resources/med/square2_split"; 
+  string filename_xml2              = data_dir + "/share/salome/resources/med/square2_split";
 
   // To remove tmp files from disk
   ParaMEDMEMTest_TmpFilesRemover aRemover;
@@ -386,7 +386,7 @@ void ParaMEDMEMTest::testIntersectionDEC_3D_(const string& srcMethod/*="P0"*/,
   if (source_group->containsMyRank())
   {
     string master = filename_xml1;
-    string meshname = "Box1Moderate";
+    string meshname = "Mesh"; //"Box1Moderate";
 
     CPPUNIT_ASSERT_NO_THROW(mesh = new MESH(MED_DRIVER,master,meshname));
     if ( srcMethod == "P1" )
@@ -431,7 +431,7 @@ void ParaMEDMEMTest::testIntersectionDEC_3D_(const string& srcMethod/*="P0"*/,
   if (target_group->containsMyRank())
   {
     string master= filename_xml2;
-    string meshname = "Box2Moderate";
+    string meshname = "Mesh"; //"Box2Moderate";
 
     CPPUNIT_ASSERT_NO_THROW(mesh = new MESH(MED_DRIVER,master,meshname));
     if ( tgtMethod == "P1" )
@@ -523,6 +523,15 @@ void ParaMEDMEMTest::testIntersectionDEC_3D_(const string& srcMethod/*="P0"*/,
     //       MPI_Bcast(&field_after_int,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
 
     //                  CPPUNIT_ASSERT_DOUBLES_EQUAL(field_before_int, field_after_int, 1e-6);
+    {
+      MEDMEM::FIELD<double>* field = parafield->getField();
+      int nb_local=field->getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
+      //cout << rank << ": Nb values = " << nb_local << " " <<srcMethod<<"->"<<tgtMethod<<endl;
+      for ( int i=1; i <=nb_local; ++i ) {
+        //cout <<rank << ": " << i << " " <<field->getValueIJ(i,1) << endl;
+        CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., field->getValueIJ(i,1), 1e-6);
+      }
+    }
 
   }