Salome HOME
Updated copyright comment
[tools/medcoupling.git] / src / ParaMEDMEM / InterpolationMatrix.cxx
index 021e9ac6eb80cc20987004e1fa7960691590e1cb..b8271a9451a0dbb2c4038f5eda05ad6dfa1ea3c7 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2019  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2024  CEA, EDF
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -48,10 +48,11 @@ namespace MEDCoupling
      Creates an empty matrix structure linking two distributed supports.
      The method must be called by all processors belonging to source
      and target groups.
-     \param source_support local support
+     \param source_field local source field
      \param source_group processor group containing the local processors
      \param target_group processor group containing the distant processors
-     \param method interpolation method
+     \param dec_options DEC options (including projection method P0, P1)
+     \param interp_options interpolation options
   */
   InterpolationMatrix::InterpolationMatrix(const MEDCoupling::ParaFIELD *source_field, 
                                            const ProcessorGroup& source_group,
@@ -639,7 +640,6 @@ namespace MEDCoupling
    * This method introduce global ids aspects in computed 'rowsPartialSumD'.
    * As precondition rowsPartialSumD.size()==policyPartial.size()==globalIdsPartial.size(). Foreach i in [0;rowsPartialSumD.size() ) rowsPartialSumD[i].size()==globalIdsPartial[i].size()
    * @param rowsPartialSumD : in parameter, Partial row sum computed for each lazy procs connected with.
-   * @param rowsPartialSumI : in parameter, Corresponding local ids for each lazy procs connected with.
    * @param globalIdsPartial : in parameter, the global numbering, of elements connected with.
    * @param globalIdsLazySideInteraction : out parameter, constituted from all global ids of lazy procs connected with.
    * @para sumCorresponding : out parameter, relative to 'globalIdsLazySideInteraction'
@@ -849,7 +849,10 @@ namespace MEDCoupling
     _mapping.prepareSendRecv();
   }
 
-
+  MCAuto<DataArrayIdType> InterpolationMatrix::retrieveNonFetchedIdsTarget(mcIdType nbTuples) const
+  {
+    return _mapping.retrieveNonFetchedIdsTarget(nbTuples);
+  }
 
   /*!
      \brief performs t=Ws, where t is the target field, s is the source field
@@ -866,12 +869,10 @@ namespace MEDCoupling
   {
     mcIdType nbcomp = ToIdType(field.getArray()->getNumberOfComponents());
     vector<double> target_value(_col_offsets.size()* nbcomp,0.0);
-
     //computing the matrix multiply on source side
     if (_source_group.containsMyRank())
       {
         mcIdType nbrows = ToIdType(_coeffs.size());
-
         // performing W.S
         // W is the intersection matrix
         // S is the source vector
@@ -894,17 +895,42 @@ namespace MEDCoupling
 
     if (_target_group.containsMyRank())
       {
-        mcIdType nbelems = field.getArray()->getNumberOfTuples() ;
-        double* value = const_cast<double*> (field.getArray()->getPointer());
-        for (mcIdType i=0; i<nbelems*nbcomp; i++)
-          {
-            value[i]=0.0;
-          }
+        field.getArray()->fillWithZero();
       }
 
     //on source side : sending  T=VT^(-1).(W.S)
     //on target side :: receiving T and storing it in field
-    _mapping.sendRecv(&target_value[0],field);
+    _mapping.sendRecv(target_value.data(),field);
+
+    if( _target_group.containsMyRank() )
+    {
+      if( this->_presence_dft_value )
+      {
+        const MCAuto<DataArrayIdType> nonFetchedEntities = _mapping.retrieveNonFetchedIdsTarget(field.getArray()->getNumberOfTuples());
+        double *fieldPtr( field.getArray()->getPointerSilent() );
+        for( const mcIdType *eltId = nonFetchedEntities->begin() ; eltId != nonFetchedEntities->end() ; ++eltId)
+          std::fill( fieldPtr + (*eltId)*nbcomp, fieldPtr + ((*eltId)+1)*nbcomp, this->_dft_value );
+      }
+    }
+  }
+
+  MCAuto<DataArrayIdType> InterpolationMatrix::retrieveNonFetchedIdsSource() const
+  {
+    MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
+    if (_source_group.containsMyRank())
+      {
+        mcIdType nbrows = ToIdType( _coeffs.size() );
+        for (mcIdType irow = 0; irow < nbrows; irow++)
+          {
+            if( _row_offsets[irow+1] == _row_offsets[irow] )
+            {
+              ret->pushBackSilent( irow );
+            }
+          }
+      }
+    else
+      THROW_IK_EXCEPTION("InterpolationMatrix::retrieveNonFetchedIdsSource : not supposed to be called target side !");
+    return ret;
   }
   
 
@@ -941,17 +967,25 @@ namespace MEDCoupling
         //T is the target vector
         for (mcIdType irow = 0; irow < nbrows; irow++)
           {
-            for (mcIdType 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 deno = _deno_reverse_multiply[irow][icol-_row_offsets[irow]];
-                for (std::size_t icomp=0; icomp<nbcomp; icomp++)
-                  {
-                    double coeff_row = source_value[colid*nbcomp+icomp];
-                    array[irow*nbcomp+icomp] += value*coeff_row/deno;
-                  }
-              }
+            if( _row_offsets[irow+1] > _row_offsets[irow] )
+            {
+              for (mcIdType 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 deno = _deno_reverse_multiply[irow][icol-_row_offsets[irow]];
+                  for (std::size_t icomp=0; icomp<nbcomp; icomp++)
+                    {
+                      double coeff_row = source_value[colid*nbcomp+icomp];
+                      array[irow*nbcomp+icomp] += value*coeff_row/deno;
+                    }
+                }
+            }
+            else
+            {
+              if( _presence_dft_value )
+                std::fill(array+irow*nbcomp,array+(irow+1)*nbcomp,this->_dft_value);
+            }
           }
       }
   }