Addition of tools for combination of part definition for // manip.

[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingFieldDiscretization.cxx

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index 83720b7abc7bd5bcb896dc87c0974308302f195b..22a30a28a1bd649a28ab19b49cc3a386a23f2d3e 100644 (file)
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -207,7 +207,7 @@ std::size_t MEDCouplingFieldDiscretization::getHeapMemorySizeWithoutChildren() c
   return 0;
 }
 
-std::vector<const BigMemoryObject *> MEDCouplingFieldDiscretization::getDirectChildren() const
+std::vector<const BigMemoryObject *> MEDCouplingFieldDiscretization::getDirectChildrenWithNull() const
 {
   return std::vector<const BigMemoryObject *>();
 }
@@ -1108,11 +1108,10 @@ std::size_t MEDCouplingFieldDiscretizationPerCell::getHeapMemorySizeWithoutChild
   return ret;
 }
 
-std::vector<const BigMemoryObject *> MEDCouplingFieldDiscretizationPerCell::getDirectChildren() const
+std::vector<const BigMemoryObject *> MEDCouplingFieldDiscretizationPerCell::getDirectChildrenWithNull() const
 {
-  std::vector<const BigMemoryObject *> ret(MEDCouplingFieldDiscretization::getDirectChildren());
-  if(_discr_per_cell)
-    ret.push_back(_discr_per_cell);
+  std::vector<const BigMemoryObject *> ret(MEDCouplingFieldDiscretization::getDirectChildrenWithNull());
+  ret.push_back(_discr_per_cell);
   return ret;
 }
 
@@ -2822,24 +2821,36 @@ DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeEvaluationMatrixO
  * \param [out] isDrift return if drift coefficients are present in the returned vector of coefficients. If different from 0 there is presence of drift coefficients.
  * \param [out] matSz the size of returned square matrix
  * \return the new result matrix to be deallocated by the caller.
+ * \sa computeMatrix
  */
 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeInverseMatrix(const MEDCouplingMesh *mesh, int& isDrift, int& matSz) const
 {
-  if(!mesh)
-    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficients : NULL input mesh !");
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
-  int nbOfPts=coords->getNumberOfTuples();
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix=coords->buildEuclidianDistanceDenseMatrix();
-  operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer());
-  // Drift
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixWithDrift=performDrift(matrix,coords,isDrift);
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv=DataArrayDouble::New();
-  matSz=nbOfPts+isDrift;
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixWithDrift(computeMatrix(mesh,isDrift,matSz));
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv(DataArrayDouble::New());
   matrixInv->alloc(matSz*matSz,1);
   INTERP_KERNEL::inverseMatrix(matrixWithDrift->getConstPointer(),matSz,matrixInv->getPointer());
   return matrixInv.retn();
 }
 
+/*!
+ * This method computes the kriging matrix.
+ * \return the new result matrix to be deallocated by the caller.
+ * \sa computeInverseMatrix
+ */
+DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeMatrix(const MEDCouplingMesh *mesh, int& isDrift, int& matSz) const
+{
+  if(!mesh)
+      throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::computeMatrix : NULL input mesh !");
+    MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords(getLocalizationOfDiscValues(mesh));
+    int nbOfPts(coords->getNumberOfTuples());
+    MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix(coords->buildEuclidianDistanceDenseMatrix());
+    operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer());
+    // Drift
+    MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixWithDrift(performDrift(matrix,coords,isDrift));
+    matSz=nbOfPts+isDrift;
+    return matrixWithDrift.retn();
+}
+
 /*!
  * This method computes coefficients to apply to each representing points of \a mesh, that is to say the nodes of \a mesh given a field array \a arr whose
  * number of tuples should be equal to the number of representing points in \a mesh.
@@ -2854,13 +2865,10 @@ DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficie
 {
   int nbRows(-1);
   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv(computeInverseMatrix(mesh,isDrift,nbRows));
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK=DataArrayDouble::New();
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK(DataArrayDouble::New());
   KnewiK->alloc(nbRows*1,1);
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr2=DataArrayDouble::New();
-  arr2->alloc(nbRows*1,1);
-  double *work=std::copy(arr->begin(),arr->end(),arr2->getPointer());
-  std::fill(work,work+isDrift,0.);
-  INTERP_KERNEL::matrixProduct(matrixInv->getConstPointer(),nbRows,nbRows,arr2->getConstPointer(),nbRows,1,KnewiK->getPointer());
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr2(PerformDriftOfVec(arr,isDrift));
+  INTERP_KERNEL::matrixProduct(matrixInv->getConstPointer(),nbRows,nbRows,arr2->getConstPointer(),arr2->getNumberOfTuples(),1,KnewiK->getPointer());
   return KnewiK.retn();
 }
 
@@ -2945,6 +2953,23 @@ DataArrayDouble *MEDCouplingFieldDiscretizationKriging::PerformDriftRect(const D
   return ret.retn();
 }
 
+/*!
+ * \return a newly allocated array having \a isDrift more tuples than \a arr.
+ * \sa computeVectorOfCoefficients
+ */
+DataArrayDouble *MEDCouplingFieldDiscretizationKriging::PerformDriftOfVec(const DataArrayDouble *arr, int isDrift)
+{
+  if(!arr || !arr->isAllocated() || arr->getNumberOfComponents()!=1)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::PerformDriftOfVec : input array must be not NULL allocated and with one component !");
+  if(isDrift<0)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::PerformDriftOfVec : isDrift parameter must be >=0 !");
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr2(DataArrayDouble::New());
+  arr2->alloc((arr->getNumberOfTuples()+isDrift)*1,1);
+  double *work(std::copy(arr->begin(),arr->end(),arr2->getPointer()));
+  std::fill(work,work+isDrift,0.);
+  return arr2.retn();
+}
+
 /*!
  * Starting from a square matrix \a matr, this method returns a newly allocated dense square matrix whose \a matr is included in returned matrix
  * in the top left corner, and in the remaining returned matrix the parameters to take into account about the kriging drift.