MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::New(TypeOfField type)
{
switch(type)
- {
+ {
case MEDCouplingFieldDiscretizationP0::TYPE:
return new MEDCouplingFieldDiscretizationP0;
case MEDCouplingFieldDiscretizationP1::TYPE:
return new MEDCouplingFieldDiscretizationKriging;
default:
throw INTERP_KERNEL::Exception("Choosen discretization is not implemented yet.");
- }
+ }
}
TypeOfField MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(const std::string& repr)
}
void MEDCouplingFieldDiscretization::setGaussLocalizationOnType(const MEDCouplingMesh *m, INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
- const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
+ const std::vector<double>& gsCoo, const std::vector<double>& wg)
{
throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
}
void MEDCouplingFieldDiscretization::setGaussLocalizationOnCells(const MEDCouplingMesh *m, const int *begin, const int *end, const std::vector<double>& refCoo,
- const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
+ const std::vector<double>& gsCoo, const std::vector<double>& wg)
{
throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
}
if(newNb>=0)//if newNb<0 the node is considered as out.
{
if(std::find_if(ptToFill+newNb*nbOfComp,ptToFill+(newNb+1)*nbOfComp,std::bind2nd(std::not_equal_to<double>(),std::numeric_limits<double>::max()))
- ==ptToFill+(newNb+1)*nbOfComp)
+ ==ptToFill+(newNb+1)*nbOfComp)
std::copy(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp);
else
{
}
void MEDCouplingFieldDiscretizationP0::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArray *>& arrays,
- const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
+ const int *old2NewBg, bool check)
{
if(!mesh)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::renumberArraysForCell : NULL input mesh !");
}
void MEDCouplingFieldDiscretizationP0::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *tupleIdsBg, const int *tupleIdsEnd,
- DataArrayInt *&cellRestriction, DataArrayInt *&trueTupleRestriction) const throw(INTERP_KERNEL::Exception)
+ DataArrayInt *&cellRestriction, DataArrayInt *&trueTupleRestriction) const
{
if(!mesh)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::computeMeshRestrictionFromTupleIds : NULL input mesh !");
* Nothing to do here.
*/
void MEDCouplingFieldDiscretizationOnNodes::renumberArraysForCell(const MEDCouplingMesh *, const std::vector<DataArray *>& arrays,
- const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
+ const int *old2NewBg, bool check)
{
}
}
void MEDCouplingFieldDiscretizationOnNodes::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *tupleIdsBg, const int *tupleIdsEnd,
- DataArrayInt *&cellRestriction, DataArrayInt *&trueTupleRestriction) const throw(INTERP_KERNEL::Exception)
+ DataArrayInt *&cellRestriction, DataArrayInt *&trueTupleRestriction) const
{
if(!mesh)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodes::computeMeshRestrictionFromTupleIds : NULL input mesh !");
/*!
* This method returns a submesh of 'mesh' instance constituting cell ids contained in array defined as an interval [start;end).
-* @param di is an array returned that specifies entity ids (here nodes ids) in mesh 'mesh' of entity in returned submesh.
+ * @param di is an array returned that specifies entity ids (here nodes ids) in mesh 'mesh' of entity in returned submesh.
* Example : The first node id of returned mesh has the (*di)[0] id in 'mesh'
*/
MEDCouplingMesh *MEDCouplingFieldDiscretizationOnNodes::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
}
void MEDCouplingFieldDiscretizationGauss::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArray *>& arrays,
- const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
+ const int *old2NewBg, bool check)
{
if(!mesh)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::renumberArraysForCell : NULL input mesh !");
INTERP_KERNEL::NormalizedCellType typ=cli.getType();
const std::vector<double>& wg=cli.getWeights();
calculator.addGaussInfo(typ,INTERP_KERNEL::CellModel::GetCellModel(typ).getDimension(),
- &cli.getGaussCoords()[0],(int)wg.size(),&cli.getRefCoords()[0],
- INTERP_KERNEL::CellModel::GetCellModel(typ).getNumberOfNodes());
+ &cli.getGaussCoords()[0],(int)wg.size(),&cli.getRefCoords()[0],
+ INTERP_KERNEL::CellModel::GetCellModel(typ).getNumberOfNodes());
//
int nbt=parts2[i]->getNumberOfTuples();
for(const int *w=parts2[i]->getConstPointer();w!=parts2[i]->getConstPointer()+nbt;w++)
}
void MEDCouplingFieldDiscretizationGauss::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *tupleIdsBg, const int *tupleIdsEnd,
- DataArrayInt *&cellRestriction, DataArrayInt *&trueTupleRestriction) const throw(INTERP_KERNEL::Exception)
+ DataArrayInt *&cellRestriction, DataArrayInt *&trueTupleRestriction) const
{
if(!mesh)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeMeshRestrictionFromTupleIds : NULL input mesh !");
}
void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType(const MEDCouplingMesh *mesh, INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
- const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
+ const std::vector<double>& gsCoo, const std::vector<double>& wg)
{
if(!mesh)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType : NULL input mesh !");
}
void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnCells(const MEDCouplingMesh *mesh, const int *begin, const int *end, const std::vector<double>& refCoo,
- const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
+ const std::vector<double>& gsCoo, const std::vector<double>& wg)
{
if(!mesh)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnCells : NULL input mesh !");
}
void MEDCouplingFieldDiscretizationGaussNE::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArray *>& arrays,
- const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
+ const int *old2NewBg, bool check)
{
if(!mesh)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::renumberArraysForCell : NULL input mesh !");
const double *MEDCouplingFieldDiscretizationGaussNE::GetWeightArrayFromGeometricType(INTERP_KERNEL::NormalizedCellType geoType, std::size_t& lgth)
{
switch(geoType)
- {
+ {
case INTERP_KERNEL::NORM_POINT1:
lgth=(int)sizeof(FGP_POINT1)/sizeof(double);
return FGP_POINT1;
return FGP_PYRA13;
default:
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::GetWeightArrayFromGeometricType : only SEG[2,3,4], TRI[3,6,7], QUAD[4,9], TETRA[4,10], PENTA[6,15], HEXA[8,20,27], PYRA[5,13] supported !");
- }
+ }
}
const double *MEDCouplingFieldDiscretizationGaussNE::GetRefCoordsFromGeometricType(INTERP_KERNEL::NormalizedCellType geoType, std::size_t& lgth)
{
switch(geoType)
- {
+ {
case INTERP_KERNEL::NORM_POINT1:
lgth=0;
return 0;
return REF_PYRA13;
default:
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::GetRefCoordsFromGeometricType : only SEG[2,3,4], TRI[3,6,7], QUAD[4,8,9], TETRA[4,10], PENTA[6,15], HEXA[8,20,27], PYRA[5,13] supported !");
- }
+ }
}
const double *MEDCouplingFieldDiscretizationGaussNE::GetLocsFromGeometricType(INTERP_KERNEL::NormalizedCellType geoType, std::size_t& lgth)
{
switch(geoType)
- {
+ {
case INTERP_KERNEL::NORM_POINT1:
{
lgth=0;
}
default:
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::GetLocsFromGeometricType : only SEG[2,3,4], TRI[3,6,7], QUAD[4,8,9], TETRA[4,10], PENTA[6,15], HEXA[8,20,27], PYRA[5,13] supported !");
- }
+ }
}
void MEDCouplingFieldDiscretizationGaussNE::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *tupleIdsBg, const int *tupleIdsEnd,
- DataArrayInt *&cellRestriction, DataArrayInt *&trueTupleRestriction) const throw(INTERP_KERNEL::Exception)
+ DataArrayInt *&cellRestriction, DataArrayInt *&trueTupleRestriction) const
{
if(!mesh)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::computeMeshRestrictionFromTupleIds : NULL input mesh !");
* \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.
{
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();
}
void MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix(int spaceDimension, int nbOfElems, double *matrixPtr) const
{
switch(spaceDimension)
- {
+ {
case 1:
{
- for(int i=0;i<nbOfElems;i++)
- {
- double val=matrixPtr[i];
- matrixPtr[i]=val*val*val;
- }
+ OperateOnDenseMatrixH3(nbOfElems,matrixPtr);
break;
}
case 2:
{
- for(int i=0;i<nbOfElems;i++)
- {
- double val=matrixPtr[i];
- if(val!=0.)
- matrixPtr[i]=val*val*log(val);
- }
+ OperateOnDenseMatrixH2Ln(nbOfElems,matrixPtr);
break;
}
case 3:
}
default:
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix : only dimension 1, 2 and 3 implemented !");
+ }
+}
+
+void MEDCouplingFieldDiscretizationKriging::OperateOnDenseMatrixH3(int nbOfElems, double *matrixPtr)
+{
+ for(int i=0;i<nbOfElems;i++)
+ {
+ double val=matrixPtr[i];
+ matrixPtr[i]=val*val*val;
+ }
+}
+
+void MEDCouplingFieldDiscretizationKriging::OperateOnDenseMatrixH2Ln(int nbOfElems, double *matrixPtr)
+{
+ for(int i=0;i<nbOfElems;i++)
+ {
+ double val=matrixPtr[i];
+ if(val!=0.)
+ matrixPtr[i]=val*val*log(val);
+ }
+}
+
+/*!
+ * Performs a drift to the rectangular input matrix \a matr.
+ * This method generate a dense matrix starting from an input dense matrix \a matr and input array \a arr.
+ * \param [in] matr The rectangular dense matrix (with only one component). The number of rows of \a matr must be equal to the number of tuples of \a arr
+ * \param [in] arr The array of coords to be appended in the input dense matrix \a matr. Typically arr is an array of coordinates.
+ * \param [out] delta the delta of number of columns between returned dense matrix and input dense matrix \a matr. \a delta is equal to number of components of \a arr + 1.
+ * \sa performDrift
+ */
+DataArrayDouble *MEDCouplingFieldDiscretizationKriging::PerformDriftRect(const DataArrayDouble *matr, const DataArrayDouble *arr, int& delta)
+{
+ if(!matr || !matr->isAllocated() || matr->getNumberOfComponents()!=1)
+ throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::PerformDriftRect : invalid input dense matrix ! Must be allocated not NULL and with exactly one component !");
+ if(!arr || !arr->isAllocated())
+ throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::PerformDriftRect : invalid input array of coordiantes ! Must be allocated and not NULL !");
+ int spaceDimension(arr->getNumberOfComponents()),nbOfPts(arr->getNumberOfTuples()),nbOfEltInMatrx(matr->getNumberOfTuples());
+ delta=spaceDimension+1;
+ int nbOfCols(nbOfEltInMatrx/nbOfPts);
+ if(nbOfEltInMatrx%nbOfPts!=0)
+ throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::PerformDriftRect : size of input dense matrix and input arrays mismatch ! NbOfElems in matrix % nb of tuples in array must be equal to 0 !");
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfPts*(nbOfCols+delta));
+ double *retPtr(ret->getPointer());
+ const double *mPtr(matr->begin()),*aPtr(arr->begin());
+ for(int i=0;i<nbOfPts;i++,aPtr+=spaceDimension,mPtr+=nbOfCols)
+ {
+ retPtr=std::copy(mPtr,mPtr+nbOfCols,retPtr);
+ *retPtr++=1.;
+ retPtr=std::copy(aPtr,aPtr+spaceDimension,retPtr);
}
+ 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();
}
/*!
* \param [in] matr input matrix whose drift part will be added
* \param [out] delta the difference between the size of the output matrix and the input matrix \a matr.
* \return a newly allocated matrix bigger than input matrix \a matr.
+ * \sa MEDCouplingFieldDiscretizationKriging::PerformDriftRect
*/
DataArrayDouble *MEDCouplingFieldDiscretizationKriging::performDrift(const DataArrayDouble *matr, const DataArrayDouble *arr, int& delta) const
{
