Updated copyright comment

[tools/medcoupling.git] / src / ParaMEDMEM / InterpolationMatrix.cxx

diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx
index 021e9ac6eb80cc20987004e1fa7960691590e1cb..b8271a9451a0dbb2c4038f5eda05ad6dfa1ea3c7 100644 (file)
--- a/src/ParaMEDMEM/InterpolationMatrix.cxx
+++ b/src/ParaMEDMEM/InterpolationMatrix.cxx
@@ -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
 //
 // 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.
      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 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,
   */
   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.
    * 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'
    * @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();
   }
 
     _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
 
   /*!
      \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);
   {
     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());
     //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
         // performing W.S
         // W is the intersection matrix
         // S is the source vector


@@ -894,17 +895,42 @@ namespace MEDCoupling
 
     if (_target_group.containsMyRank())
       {
 
     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
       }
 
     //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++)
           {
         //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);
+            }

           }
       }
   }
           }
       }
   }