From: ageay Date: Tue, 20 Apr 2010 07:27:36 +0000 (+0000) Subject: Merge from V5_1_main X-Git-Tag: V5_1_4rc1~7 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=665726568223952c269066ccececcdf6474595d7;p=tools%2Fmedcoupling.git Merge from V5_1_main --- diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx index dccf9b346..51e894818 100644 --- a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx +++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx @@ -19,8 +19,11 @@ #include "MEDCouplingExtrudedMesh.hxx" #include "MEDCouplingUMesh.hxx" #include "MEDCouplingMemArray.hxx" +#include "MEDCouplingFieldDouble.hxx" #include "CellModel.hxx" +#include "InterpolationUtils.hxx" + #include #include #include @@ -332,15 +335,32 @@ int MEDCouplingExtrudedMesh::findCorrespCellByNodalConn(const std::vector& * are created ('m1r' and 'm2r') that can be used for interpolation. * @param m1 input mesh with meshDim==1 and spaceDim==3 * @param m2 input mesh with meshDim==1 and spaceDim==3 - * @param m1r output mesh with ref count equal to 1 with meshDim==1 and spaceDim==2 - * @param m2r output mesh with ref count equal to 1 with meshDim==1 and spaceDim==2 + * @param eps tolerance acceptable to determine compatibility + * @param m1r output mesh with ref count equal to 1 with meshDim==1 and spaceDim==1 + * @param m2r output mesh with ref count equal to 1 with meshDim==1 and spaceDim==1 + * @param v is the output normalized vector of the common direction of 'm1' and 'm2' * @throw in case that m1 and m2 are not compatible each other. */ -void MEDCouplingExtrudedMesh::project1DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, - MEDCouplingUMesh *&m1r, MEDCouplingUMesh *&m2r) throw(INTERP_KERNEL::Exception) +void MEDCouplingExtrudedMesh::project1DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, + MEDCouplingUMesh *&m1r, MEDCouplingUMesh *&m2r, double *v) throw(INTERP_KERNEL::Exception) { - m1r=0; - m2r=0; + if(m1->getSpaceDimension()!=3 || m1->getSpaceDimension()!=3) + throw INTERP_KERNEL::Exception("Input meshes are expected to have a spaceDim==3 for projec1D !"); + m1r=m1->clone(true); + m2r=m2->clone(true); + m1r->changeSpaceDimension(1); + m2r->changeSpaceDimension(1); + std::vector c; + std::vector ref,ref2; + m1->getNodeIdsOfCell(0,c); + m1->getCoordinatesOfNode(c[0],ref); + m1->getCoordinatesOfNode(c[1],ref2); + std::transform(ref2.begin(),ref2.end(),ref.begin(),v,std::minus()); + double n=INTERP_KERNEL::norm<3>(v); + std::transform(v,v+3,v,std::bind2nd(std::multiplies(),1/n)); + m1->project1D(&ref[0],v,eps,m1r->getCoords()->getPointer()); + m2->project1D(&ref[0],v,eps,m2r->getCoords()->getPointer()); + } void MEDCouplingExtrudedMesh::rotate(const double *center, const double *vector, double angle) diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx index f60dea159..0325818ed 100644 --- a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx +++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx @@ -56,8 +56,8 @@ namespace ParaMEDMEM int getCellContainingPoint(const double *pos, double eps) const; static int findCorrespCellByNodalConn(const std::vector& nodalConnec, const int *revNodalPtr, const int *revNodalIndxPtr) throw(INTERP_KERNEL::Exception); - static void project1DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, - MEDCouplingUMesh *&m1r, MEDCouplingUMesh *&m2r) throw(INTERP_KERNEL::Exception); + static void project1DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, + MEDCouplingUMesh *&m1r, MEDCouplingUMesh *&m2r, double *v) throw(INTERP_KERNEL::Exception); void rotate(const double *center, const double *vector, double angle); void translate(const double *vector); MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const; diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx index ed4e77462..9dbff0c21 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx @@ -416,6 +416,13 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::addFields(const MEDCouplingField return ret; } +void MEDCouplingFieldDouble::operator+=(const MEDCouplingFieldDouble& other) +{ + if(!areCompatible(&other)) + throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply += on them !"); + _time_discr->addEqual(other._time_discr); +} + MEDCouplingFieldDouble *MEDCouplingFieldDouble::substractFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) { if(!f1->areCompatible(f2)) @@ -426,16 +433,30 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::substractFields(const MEDCouplin return ret; } +void MEDCouplingFieldDouble::operator-=(const MEDCouplingFieldDouble& other) +{ + if(!areCompatible(&other)) + throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply -= on them !"); + _time_discr->substractEqual(other._time_discr); +} + MEDCouplingFieldDouble *MEDCouplingFieldDouble::multiplyFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) { if(!f1->areCompatibleForMul(f2)) - throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to applymultiplyFields on them !"); + throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply multiplyFields on them !"); MEDCouplingTimeDiscretization *td=f1->_time_discr->multiply(f2->_time_discr); MEDCouplingFieldDouble *ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->getTypeOfField()); ret->setMesh(f1->getMesh()); return ret; } +void MEDCouplingFieldDouble::operator*=(const MEDCouplingFieldDouble& other) +{ + if(!areCompatibleForMul(&other)) + throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply *= on them !"); + _time_discr->multiplyEqual(other._time_discr); +} + MEDCouplingFieldDouble *MEDCouplingFieldDouble::divideFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) { if(!f1->areCompatible(f2)) @@ -445,3 +466,10 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::divideFields(const MEDCouplingFi ret->setMesh(f1->getMesh()); return ret; } + +void MEDCouplingFieldDouble::operator/=(const MEDCouplingFieldDouble& other) +{ + if(!areCompatible(&other)) + throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply /= on them !"); + _time_discr->divideEqual(other._time_discr); +} diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.hxx b/src/MEDCoupling/MEDCouplingFieldDouble.hxx index 8e92075ac..03ff33da0 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.hxx @@ -74,12 +74,16 @@ namespace ParaMEDMEM bool mergeNodes(double eps); static MEDCouplingFieldDouble *mergeFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2); MEDCouplingFieldDouble *operator+(const MEDCouplingFieldDouble& other) const { return addFields(this,&other); } + void operator+=(const MEDCouplingFieldDouble& other); static MEDCouplingFieldDouble *addFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2); MEDCouplingFieldDouble *operator-(const MEDCouplingFieldDouble& other) const { return substractFields(this,&other); } + void operator-=(const MEDCouplingFieldDouble& other); static MEDCouplingFieldDouble *substractFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2); MEDCouplingFieldDouble *operator*(const MEDCouplingFieldDouble& other) const { return multiplyFields(this,&other); } + void operator*=(const MEDCouplingFieldDouble& other); static MEDCouplingFieldDouble *multiplyFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2); MEDCouplingFieldDouble *operator/(const MEDCouplingFieldDouble& other) const { return divideFields(this,&other); } + void operator/=(const MEDCouplingFieldDouble& other); static MEDCouplingFieldDouble *divideFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2); private: MEDCouplingFieldDouble(TypeOfField type, TypeOfTimeDiscretization td); diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index 81c34fc2f..2021c1795 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -146,6 +146,18 @@ DataArrayDouble *DataArrayDouble::add(const DataArrayDouble *a1, const DataArray return ret; } +void DataArrayDouble::addEqual(const DataArrayDouble *other) +{ + int nbOfComp=getNumberOfComponents(); + if(nbOfComp!=other->getNumberOfComponents()) + throw INTERP_KERNEL::Exception("Nb of components mismatch for array add !"); + int nbOfTuple=getNumberOfTuples(); + if(nbOfTuple!=other->getNumberOfTuples()) + throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array add !"); + std::transform(getConstPointer(),getConstPointer()+nbOfTuple*nbOfComp,other->getConstPointer(),getPointer(),std::plus()); + declareAsNew(); +} + DataArrayDouble *DataArrayDouble::substract(const DataArrayDouble *a1, const DataArrayDouble *a2) { int nbOfComp=a1->getNumberOfComponents(); @@ -161,6 +173,18 @@ DataArrayDouble *DataArrayDouble::substract(const DataArrayDouble *a1, const Dat return ret; } +void DataArrayDouble::substractEqual(const DataArrayDouble *other) +{ + int nbOfComp=getNumberOfComponents(); + if(nbOfComp!=other->getNumberOfComponents()) + throw INTERP_KERNEL::Exception("Nb of components mismatch for array substract !"); + int nbOfTuple=getNumberOfTuples(); + if(nbOfTuple!=other->getNumberOfTuples()) + throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array substract !"); + std::transform(getConstPointer(),getConstPointer()+nbOfTuple*nbOfComp,other->getConstPointer(),getPointer(),std::minus()); + declareAsNew(); +} + DataArrayDouble *DataArrayDouble::multiply(const DataArrayDouble *a1, const DataArrayDouble *a2) { int nbOfTuple=a1->getNumberOfTuples(); @@ -208,6 +232,36 @@ DataArrayDouble *DataArrayDouble::multiply(const DataArrayDouble *a1, const Data return ret; } +void DataArrayDouble::multiplyEqual(const DataArrayDouble *other) +{ + int nbOfTuple=getNumberOfTuples(); + int nbOfTuple2=other->getNumberOfTuples(); + int nbOfComp=getNumberOfComponents(); + int nbOfComp2=other->getNumberOfComponents(); + if(nbOfTuple!=nbOfTuple2) + throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array multiplyEqual !"); + DataArrayDouble *ret=0; + if(nbOfComp==nbOfComp2) + { + ret=DataArrayDouble::New(); + ret->alloc(nbOfTuple,nbOfComp); + std::transform(getConstPointer(),getConstPointer()+nbOfTuple*nbOfComp,other->getConstPointer(),getPointer(),std::multiplies()); + } + else + { + if(nbOfComp2==1) + { + const double *ptr=other->getConstPointer(); + double *myPtr=getPointer(); + for(int i=0;i(),ptr[i])); + } + else + throw INTERP_KERNEL::Exception("Nb of components mismatch for array multiplyEqual !"); + } + declareAsNew(); +} + DataArrayDouble *DataArrayDouble::divide(const DataArrayDouble *a1, const DataArrayDouble *a2) { int nbOfComp=a1->getNumberOfComponents(); @@ -223,6 +277,18 @@ DataArrayDouble *DataArrayDouble::divide(const DataArrayDouble *a1, const DataAr return ret; } +void DataArrayDouble::divideEqual(const DataArrayDouble *other) +{ + int nbOfComp=getNumberOfComponents(); + if(nbOfComp!=other->getNumberOfComponents()) + throw INTERP_KERNEL::Exception("Nb of components mismatch for array divideEqual !"); + int nbOfTuple=getNumberOfTuples(); + if(nbOfTuple!=other->getNumberOfTuples()) + throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array divideEqual !"); + std::transform(getConstPointer(),getConstPointer()+nbOfTuple*nbOfComp,other->getConstPointer(),getPointer(),std::divides()); + declareAsNew(); +} + DataArrayInt *DataArrayInt::New() { return new DataArrayInt; @@ -301,3 +367,22 @@ DataArrayInt *DataArrayInt::aggregate(const DataArrayInt *a1, const DataArrayInt return ret; } +int *DataArrayInt::checkAndPreparePermutation(const int *start, const int *end) +{ + int sz=std::distance(start,end); + int *ret=new int[sz]; + int *work=new int[sz]; + std::copy(start,end,work); + std::sort(work,work+sz); + if(std::unique(work,work+sz)!=work+sz) + { + delete [] work; + delete [] ret; + throw INTERP_KERNEL::Exception("Some elements are equals in the specified array !"); + } + int *iter2=ret; + for(const int *iter=start;iter!=end;iter++,iter2++) + *iter2=std::distance(work,std::find(work,work+sz,*iter)); + delete [] work; + return ret; +} diff --git a/src/MEDCoupling/MEDCouplingMemArray.hxx b/src/MEDCoupling/MEDCouplingMemArray.hxx index 68975339c..2ef9a2b28 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.hxx +++ b/src/MEDCoupling/MEDCouplingMemArray.hxx @@ -121,9 +121,13 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void checkNoNullValues() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT static DataArrayDouble *aggregate(const DataArrayDouble *a1, const DataArrayDouble *a2); MEDCOUPLING_EXPORT static DataArrayDouble *add(const DataArrayDouble *a1, const DataArrayDouble *a2); + MEDCOUPLING_EXPORT void addEqual(const DataArrayDouble *other); MEDCOUPLING_EXPORT static DataArrayDouble *substract(const DataArrayDouble *a1, const DataArrayDouble *a2); + MEDCOUPLING_EXPORT void substractEqual(const DataArrayDouble *other); MEDCOUPLING_EXPORT static DataArrayDouble *multiply(const DataArrayDouble *a1, const DataArrayDouble *a2); + MEDCOUPLING_EXPORT void multiplyEqual(const DataArrayDouble *other); MEDCOUPLING_EXPORT static DataArrayDouble *divide(const DataArrayDouble *a1, const DataArrayDouble *a2); + MEDCOUPLING_EXPORT void divideEqual(const DataArrayDouble *other); //! nothing to do here because this class does not aggregate any TimeLabel instance. MEDCOUPLING_EXPORT void updateTime() { } private: @@ -153,6 +157,8 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void writeOnPlace(int id, int element0, const int *others, int sizeOfOthers) { _mem.writeOnPlace(id,element0,others,sizeOfOthers); } //! nothing to do here because this class does not aggregate any TimeLabel instance. MEDCOUPLING_EXPORT void updateTime() { } + public: + MEDCOUPLING_EXPORT static int *checkAndPreparePermutation(const int *start, const int *end); private: DataArrayInt() { } private: diff --git a/src/MEDCoupling/MEDCouplingPointSet.cxx b/src/MEDCoupling/MEDCouplingPointSet.cxx index 89d009306..099b668bd 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.cxx +++ b/src/MEDCoupling/MEDCouplingPointSet.cxx @@ -329,6 +329,44 @@ void MEDCouplingPointSet::scale(const double *point, double factor) updateTime(); } +/*! + * This method is only available for already defined coordinates. + * If not an INTERP_KERNEL::Exception is thrown. The 'newSpaceDim' input must be greater or equal to 1. + * This method simply convert this to newSpaceDim space : + * - by putting a 0. for each \f$ i^{th} \f$ components of each coord of nodes so that i>=getSpaceDim(), if 'newSpaceDim'>getSpaceDimsion() + * - by ignoring each \f$ i^{th} \f$ components of each coord of nodes so that i >= 'newSpaceDim', 'newSpaceDim'=1 !"); + int oldSpaceDim=getSpaceDimension(); + if(newSpaceDim==oldSpaceDim) + return ; + DataArrayDouble *newCoords=DataArrayDouble::New(); + newCoords->alloc(getCoords()->getNumberOfTuples(),newSpaceDim); + const double *oldc=getCoords()->getConstPointer(); + double *nc=newCoords->getPointer(); + int nbOfNodes=getNumberOfNodes(); + int dim=std::min(oldSpaceDim,newSpaceDim); + for(int i=0;isetName(getCoords()->getName().c_str()); + for(int i=0;isetInfoOnComponent(i,getCoords()->getInfoOnComponent(i).c_str()); + setCoords(newCoords); + newCoords->decrRef(); +} + /*! * This method try to substitute this->_coords with other._coords if arrays match. * This method potentially modifies 'this' if it succeeds, otherway an exception is thrown. diff --git a/src/MEDCoupling/MEDCouplingPointSet.hxx b/src/MEDCoupling/MEDCouplingPointSet.hxx index 1b1024aa9..1393966c6 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.hxx +++ b/src/MEDCoupling/MEDCouplingPointSet.hxx @@ -54,6 +54,7 @@ namespace ParaMEDMEM void rotate(const double *center, const double *vector, double angle); void translate(const double *vector); void scale(const double *point, double factor); + void changeSpaceDimension(int newSpaceDim) throw(INTERP_KERNEL::Exception); void tryToShareSameCoords(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception); void findNodesOnPlane(const double *pt, const double *vec, double eps, std::vector& nodes) const throw(INTERP_KERNEL::Exception); static DataArrayDouble *mergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2); diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index 4a21b5c21..60e57d070 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -23,6 +23,7 @@ #include "MEDCouplingExtrudedMesh.hxx" #include "MEDCouplingNormalizedUnstructuredMesh.txx" +#include "Interpolation1D.txx" #include "Interpolation2DCurve.txx" #include "Interpolation2D.txx" #include "Interpolation3D.txx" @@ -293,20 +294,23 @@ int MEDCouplingRemapper::prepareEE(const char *method) if(methC!="P0P0") throw INTERP_KERNEL::Exception("Only P0P0 method implemented for Extruded/Extruded meshes !"); INTERP_KERNEL::Interpolation::checkAndSplitInterpolationMethod(method,_src_method,_target_method); - MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh->getMesh2D()); - MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh->getMesh2D()); - INTERP_KERNEL::Interpolation2D interpolation(*this); + MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D()); + MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D()); + INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this); std::vector > matrix2D; - int nbCols2D=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,method); + int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,method); MEDCouplingUMesh *s1D,*t1D; - MEDCouplingExtrudedMesh::project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),s1D,t1D); - MEDCouplingNormalizedUnstructuredMesh<2,1> s1DWrapper(s1D); - MEDCouplingNormalizedUnstructuredMesh<2,1> t1DWrapper(t1D); + double v[3]; + MEDCouplingExtrudedMesh::project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v); + MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D); + MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D); std::vector > matrix1D; - int nbCols1D=interpolation.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,method); - INTERP_KERNEL::Interpolation2DCurve myInterpolator; - // - _matrix.resize(matrix2D.size()*matrix1D.size()); + INTERP_KERNEL::Interpolation1D interpolation1D(*this); + int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,method); + s1D->decrRef(); + t1D->decrRef(); + buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D, + target_mesh->getMesh3DIds()->getConstPointer()); // _deno_multiply.clear(); _deno_multiply.resize(_matrix.size()); @@ -522,3 +526,41 @@ void MEDCouplingRemapper::computeColSumAndRowSum(const std::vector >& m1D, + const std::vector< std::map >& m2D, + const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc, + const int *corrCellIdTrg) +{ + int nbOf2DCellsTrg=m2D.size(); + int nbOf1DCellsTrg=m1D.size(); + int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg; + _matrix.resize(nbOf3DCellsTrg); + int id2R=0; + for(std::vector< std::map >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++) + { + for(std::map::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++) + { + int id1R=0; + for(std::vector< std::map >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++) + { + for(std::map::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++) + { + _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second); + } + } + } + } +} + +void MEDCouplingRemapper::printMatrix(const std::vector >& m) +{ + int id=0; + for(std::vector >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++) + { + std::cout << "Target Cell # " << id << " : "; + for(std::map::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++) + std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), "; + std::cout << std::endl; + } +} diff --git a/src/MEDCoupling/MEDCouplingRemapper.hxx b/src/MEDCoupling/MEDCouplingRemapper.hxx index 461ce1711..d1b0317d0 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.hxx +++ b/src/MEDCoupling/MEDCouplingRemapper.hxx @@ -49,6 +49,8 @@ namespace ParaMEDMEM MEDCOUPLINGREMAPPER_EXPORT bool setOptionInt(const std::string& key, int value); MEDCOUPLINGREMAPPER_EXPORT bool setOptionDouble(const std::string& key, double value); MEDCOUPLINGREMAPPER_EXPORT bool setOptionString(const std::string& key, std::string& value); + public: + static void printMatrix(const std::vector >& m); private: int prepareUU(const char *method); int prepareEE(const char *method); @@ -58,6 +60,10 @@ namespace ParaMEDMEM void computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField); void computeProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer); void computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer); + void buildFinalInterpolationMatrixByConvolution(const std::vector< std::map >& m1D, + const std::vector< std::map >& m2D, + const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc, + const int *corrCellIdTrg); static void reverseMatrix(const std::vector >& matIn, int nbColsMatIn, std::vector >& matOut); static void computeRowSumAndColSum(const std::vector >& matrixDeno, diff --git a/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx b/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx index 43ecbdca0..7bf11d3ef 100644 --- a/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx @@ -358,7 +358,7 @@ MEDCouplingTimeDiscretization *MEDCouplingNoTimeLabel::aggregate(const MEDCoupli { const MEDCouplingNoTimeLabel *otherC=dynamic_cast(other); if(!otherC) - throw INTERP_KERNEL::Exception("aggregation on mismatched time discretization !"); + throw INTERP_KERNEL::Exception("NoTimeLabel::aggregation on mismatched time discretization !"); MEDCouplingNoTimeLabel *ret=new MEDCouplingNoTimeLabel; ret->setTimeTolerance(getTimeTolerance()); DataArrayDouble *arr=DataArrayDouble::aggregate(getArray(),other->getArray()); @@ -371,7 +371,7 @@ MEDCouplingTimeDiscretization *MEDCouplingNoTimeLabel::add(const MEDCouplingTime { const MEDCouplingNoTimeLabel *otherC=dynamic_cast(other); if(!otherC) - throw INTERP_KERNEL::Exception("add on mismatched time discretization !"); + throw INTERP_KERNEL::Exception("NoTimeLabel::add on mismatched time discretization !"); MEDCouplingNoTimeLabel *ret=new MEDCouplingNoTimeLabel; DataArrayDouble *arr=DataArrayDouble::add(getArray(),other->getArray()); ret->setArray(arr,0); @@ -379,11 +379,19 @@ MEDCouplingTimeDiscretization *MEDCouplingNoTimeLabel::add(const MEDCouplingTime return ret; } +void MEDCouplingNoTimeLabel::addEqual(const MEDCouplingTimeDiscretization *other) +{ + const MEDCouplingNoTimeLabel *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("NoTimeLabel::addEqual on mismatched time discretization !"); + getArray()->addEqual(other->getArray()); +} + MEDCouplingTimeDiscretization *MEDCouplingNoTimeLabel::substract(const MEDCouplingTimeDiscretization *other) const { const MEDCouplingNoTimeLabel *otherC=dynamic_cast(other); if(!otherC) - throw INTERP_KERNEL::Exception("substract on mismatched time discretization !"); + throw INTERP_KERNEL::Exception("NoTimeLabel::substract on mismatched time discretization !"); MEDCouplingNoTimeLabel *ret=new MEDCouplingNoTimeLabel; DataArrayDouble *arr=DataArrayDouble::substract(getArray(),other->getArray()); ret->setArray(arr,0); @@ -391,11 +399,19 @@ MEDCouplingTimeDiscretization *MEDCouplingNoTimeLabel::substract(const MEDCoupli return ret; } +void MEDCouplingNoTimeLabel::substractEqual(const MEDCouplingTimeDiscretization *other) +{ + const MEDCouplingNoTimeLabel *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("NoTimeLabel::substractEqual on mismatched time discretization !"); + getArray()->substractEqual(other->getArray()); +} + MEDCouplingTimeDiscretization *MEDCouplingNoTimeLabel::multiply(const MEDCouplingTimeDiscretization *other) const { const MEDCouplingNoTimeLabel *otherC=dynamic_cast(other); if(!otherC) - throw INTERP_KERNEL::Exception("multiply on mismatched time discretization !"); + throw INTERP_KERNEL::Exception("NoTimeLabel::multiply on mismatched time discretization !"); MEDCouplingNoTimeLabel *ret=new MEDCouplingNoTimeLabel; DataArrayDouble *arr=DataArrayDouble::multiply(getArray(),other->getArray()); ret->setArray(arr,0); @@ -403,6 +419,14 @@ MEDCouplingTimeDiscretization *MEDCouplingNoTimeLabel::multiply(const MEDCouplin return ret; } +void MEDCouplingNoTimeLabel::multiplyEqual(const MEDCouplingTimeDiscretization *other) +{ + const MEDCouplingNoTimeLabel *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("NoTimeLabel::multiplyEqual on mismatched time discretization !"); + getArray()->multiplyEqual(other->getArray()); +} + MEDCouplingTimeDiscretization *MEDCouplingNoTimeLabel::divide(const MEDCouplingTimeDiscretization *other) const { const MEDCouplingNoTimeLabel *otherC=dynamic_cast(other); @@ -415,6 +439,14 @@ MEDCouplingTimeDiscretization *MEDCouplingNoTimeLabel::divide(const MEDCouplingT return ret; } +void MEDCouplingNoTimeLabel::divideEqual(const MEDCouplingTimeDiscretization *other) +{ + const MEDCouplingNoTimeLabel *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("NoTimeLabel::divideEqual on mismatched time discretization !"); + getArray()->divideEqual(other->getArray()); +} + MEDCouplingTimeDiscretization *MEDCouplingNoTimeLabel::performCpy(bool deepCpy) const { return new MEDCouplingNoTimeLabel(*this,deepCpy); @@ -543,7 +575,7 @@ MEDCouplingTimeDiscretization *MEDCouplingWithTimeStep::aggregate(const MEDCoupl { const MEDCouplingWithTimeStep *otherC=dynamic_cast(other); if(!otherC) - throw INTERP_KERNEL::Exception("aggregation on mismatched time discretization !"); + throw INTERP_KERNEL::Exception("WithTimeStep::aggregation on mismatched time discretization !"); MEDCouplingWithTimeStep *ret=new MEDCouplingWithTimeStep; ret->setTimeTolerance(getTimeTolerance()); DataArrayDouble *arr=DataArrayDouble::aggregate(getArray(),other->getArray()); @@ -559,7 +591,7 @@ MEDCouplingTimeDiscretization *MEDCouplingWithTimeStep::add(const MEDCouplingTim { const MEDCouplingWithTimeStep *otherC=dynamic_cast(other); if(!otherC) - throw INTERP_KERNEL::Exception("add on mismatched time discretization !"); + throw INTERP_KERNEL::Exception("WithTimeStep::add on mismatched time discretization !"); MEDCouplingWithTimeStep *ret=new MEDCouplingWithTimeStep; DataArrayDouble *arr=DataArrayDouble::add(getArray(),other->getArray()); ret->setArray(arr,0); @@ -570,11 +602,19 @@ MEDCouplingTimeDiscretization *MEDCouplingWithTimeStep::add(const MEDCouplingTim return ret; } +void MEDCouplingWithTimeStep::addEqual(const MEDCouplingTimeDiscretization *other) +{ + const MEDCouplingWithTimeStep *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("WithTimeStep::addEqual on mismatched time discretization !"); + getArray()->addEqual(other->getArray()); +} + MEDCouplingTimeDiscretization *MEDCouplingWithTimeStep::substract(const MEDCouplingTimeDiscretization *other) const { const MEDCouplingWithTimeStep *otherC=dynamic_cast(other); if(!otherC) - throw INTERP_KERNEL::Exception("substract on mismatched time discretization !"); + throw INTERP_KERNEL::Exception("WithTimeStep::substract on mismatched time discretization !"); MEDCouplingWithTimeStep *ret=new MEDCouplingWithTimeStep; DataArrayDouble *arr=DataArrayDouble::substract(getArray(),other->getArray()); ret->setArray(arr,0); @@ -585,11 +625,19 @@ MEDCouplingTimeDiscretization *MEDCouplingWithTimeStep::substract(const MEDCoupl return ret; } +void MEDCouplingWithTimeStep::substractEqual(const MEDCouplingTimeDiscretization *other) +{ + const MEDCouplingWithTimeStep *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("WithTimeStep::substractEqual on mismatched time discretization !"); + getArray()->substractEqual(other->getArray()); +} + MEDCouplingTimeDiscretization *MEDCouplingWithTimeStep::multiply(const MEDCouplingTimeDiscretization *other) const { const MEDCouplingWithTimeStep *otherC=dynamic_cast(other); if(!otherC) - throw INTERP_KERNEL::Exception("multiply on mismatched time discretization !"); + throw INTERP_KERNEL::Exception("WithTimeStep::multiply on mismatched time discretization !"); MEDCouplingWithTimeStep *ret=new MEDCouplingWithTimeStep; DataArrayDouble *arr=DataArrayDouble::multiply(getArray(),other->getArray()); ret->setArray(arr,0); @@ -600,11 +648,19 @@ MEDCouplingTimeDiscretization *MEDCouplingWithTimeStep::multiply(const MEDCoupli return ret; } +void MEDCouplingWithTimeStep::multiplyEqual(const MEDCouplingTimeDiscretization *other) +{ + const MEDCouplingWithTimeStep *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("WithTimeStep::multiplyEqual on mismatched time discretization !"); + getArray()->multiplyEqual(other->getArray()); +} + MEDCouplingTimeDiscretization *MEDCouplingWithTimeStep::divide(const MEDCouplingTimeDiscretization *other) const { const MEDCouplingWithTimeStep *otherC=dynamic_cast(other); if(!otherC) - throw INTERP_KERNEL::Exception("divide on mismatched time discretization !"); + throw INTERP_KERNEL::Exception("WithTimeStep::divide on mismatched time discretization !"); MEDCouplingWithTimeStep *ret=new MEDCouplingWithTimeStep; DataArrayDouble *arr=DataArrayDouble::divide(getArray(),other->getArray()); ret->setArray(arr,0); @@ -615,6 +671,14 @@ MEDCouplingTimeDiscretization *MEDCouplingWithTimeStep::divide(const MEDCoupling return ret; } +void MEDCouplingWithTimeStep::divideEqual(const MEDCouplingTimeDiscretization *other) +{ + const MEDCouplingWithTimeStep *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("WithTimeStep::divideEqual on mismatched time discretization !"); + getArray()->divideEqual(other->getArray()); +} + MEDCouplingTimeDiscretization *MEDCouplingWithTimeStep::performCpy(bool deepCpy) const { return new MEDCouplingWithTimeStep(*this,deepCpy); @@ -815,7 +879,7 @@ MEDCouplingTimeDiscretization *MEDCouplingConstOnTimeInterval::aggregate(const M { const MEDCouplingConstOnTimeInterval *otherC=dynamic_cast(other); if(!otherC) - throw INTERP_KERNEL::Exception("aggregation on mismatched time discretization !"); + throw INTERP_KERNEL::Exception("ConstOnTimeInterval::aggregation on mismatched time discretization !"); MEDCouplingConstOnTimeInterval *ret=new MEDCouplingConstOnTimeInterval; ret->setTimeTolerance(getTimeTolerance()); DataArrayDouble *arr=DataArrayDouble::aggregate(getArray(),other->getArray()); @@ -833,7 +897,7 @@ MEDCouplingTimeDiscretization *MEDCouplingConstOnTimeInterval::add(const MEDCoup { const MEDCouplingConstOnTimeInterval *otherC=dynamic_cast(other); if(!otherC) - throw INTERP_KERNEL::Exception("add on mismatched time discretization !"); + throw INTERP_KERNEL::Exception("ConstOnTimeInterval::add on mismatched time discretization !"); MEDCouplingConstOnTimeInterval *ret=new MEDCouplingConstOnTimeInterval; DataArrayDouble *arr=DataArrayDouble::add(getArray(),other->getArray()); ret->setArray(arr,0); @@ -845,12 +909,20 @@ MEDCouplingTimeDiscretization *MEDCouplingConstOnTimeInterval::add(const MEDCoup ret->setEndTime(tmp3,tmp1,tmp2); return ret; } + +void MEDCouplingConstOnTimeInterval::addEqual(const MEDCouplingTimeDiscretization *other) +{ + const MEDCouplingConstOnTimeInterval *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("ConstOnTimeInterval::addEqual on mismatched time discretization !"); + getArray()->addEqual(other->getArray()); +} MEDCouplingTimeDiscretization *MEDCouplingConstOnTimeInterval::substract(const MEDCouplingTimeDiscretization *other) const { const MEDCouplingConstOnTimeInterval *otherC=dynamic_cast(other); if(!otherC) - throw INTERP_KERNEL::Exception("substract on mismatched time discretization !"); + throw INTERP_KERNEL::Exception("ConstOnTimeInterval::substract on mismatched time discretization !"); MEDCouplingConstOnTimeInterval *ret=new MEDCouplingConstOnTimeInterval; DataArrayDouble *arr=DataArrayDouble::substract(getArray(),other->getArray()); ret->setArray(arr,0); @@ -863,6 +935,14 @@ MEDCouplingTimeDiscretization *MEDCouplingConstOnTimeInterval::substract(const M return ret; } +void MEDCouplingConstOnTimeInterval::substractEqual(const MEDCouplingTimeDiscretization *other) +{ + const MEDCouplingConstOnTimeInterval *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("ConstOnTimeInterval::substractEqual on mismatched time discretization !"); + getArray()->substractEqual(other->getArray()); +} + MEDCouplingTimeDiscretization *MEDCouplingConstOnTimeInterval::multiply(const MEDCouplingTimeDiscretization *other) const { const MEDCouplingConstOnTimeInterval *otherC=dynamic_cast(other); @@ -880,6 +960,14 @@ MEDCouplingTimeDiscretization *MEDCouplingConstOnTimeInterval::multiply(const ME return ret; } +void MEDCouplingConstOnTimeInterval::multiplyEqual(const MEDCouplingTimeDiscretization *other) +{ + const MEDCouplingConstOnTimeInterval *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("ConstOnTimeInterval::multiplyEqual on mismatched time discretization !"); + getArray()->multiplyEqual(other->getArray()); +} + MEDCouplingTimeDiscretization *MEDCouplingConstOnTimeInterval::divide(const MEDCouplingTimeDiscretization *other) const { const MEDCouplingConstOnTimeInterval *otherC=dynamic_cast(other); @@ -897,6 +985,14 @@ MEDCouplingTimeDiscretization *MEDCouplingConstOnTimeInterval::divide(const MEDC return ret; } +void MEDCouplingConstOnTimeInterval::divideEqual(const MEDCouplingTimeDiscretization *other) +{ + const MEDCouplingConstOnTimeInterval *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("ConstOnTimeInterval::divideEqual on mismatched time discretization !"); + getArray()->divideEqual(other->getArray()); +} + MEDCouplingTwoTimeSteps::MEDCouplingTwoTimeSteps():_start_time(0.),_end_time(0.),_start_dt(-1),_end_dt(-1),_start_it(-1),_end_it(-1),_end_array(0) { } diff --git a/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx b/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx index e13742ee0..e794ae46a 100644 --- a/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx @@ -47,9 +47,13 @@ namespace ParaMEDMEM virtual TypeOfTimeDiscretization getEnum() const = 0; virtual MEDCouplingTimeDiscretization *aggregate(const MEDCouplingTimeDiscretization *other) const = 0; virtual MEDCouplingTimeDiscretization *add(const MEDCouplingTimeDiscretization *other) const = 0; + virtual void addEqual(const MEDCouplingTimeDiscretization *other) = 0; virtual MEDCouplingTimeDiscretization *substract(const MEDCouplingTimeDiscretization *other) const = 0; + virtual void substractEqual(const MEDCouplingTimeDiscretization *other) = 0; virtual MEDCouplingTimeDiscretization *multiply(const MEDCouplingTimeDiscretization *other) const = 0; + virtual void multiplyEqual(const MEDCouplingTimeDiscretization *other) = 0; virtual MEDCouplingTimeDiscretization *divide(const MEDCouplingTimeDiscretization *other) const = 0; + virtual void divideEqual(const MEDCouplingTimeDiscretization *other) = 0; virtual void getTinySerializationIntInformation(std::vector& tinyInfo) const; virtual void getTinySerializationDbleInformation(std::vector& tinyInfo) const; virtual void getTinySerializationStrInformation(std::vector& tinyInfo) const; @@ -99,9 +103,13 @@ namespace ParaMEDMEM TypeOfTimeDiscretization getEnum() const { return DISCRETIZATION; } MEDCouplingTimeDiscretization *aggregate(const MEDCouplingTimeDiscretization *other) const; MEDCouplingTimeDiscretization *add(const MEDCouplingTimeDiscretization *other) const; + void addEqual(const MEDCouplingTimeDiscretization *other); MEDCouplingTimeDiscretization *substract(const MEDCouplingTimeDiscretization *other) const; + void substractEqual(const MEDCouplingTimeDiscretization *other); MEDCouplingTimeDiscretization *multiply(const MEDCouplingTimeDiscretization *other) const; + void multiplyEqual(const MEDCouplingTimeDiscretization *other); MEDCouplingTimeDiscretization *divide(const MEDCouplingTimeDiscretization *other) const; + void divideEqual(const MEDCouplingTimeDiscretization *other); bool isEqual(const MEDCouplingTimeDiscretization *other, double prec) const; bool areCompatible(const MEDCouplingTimeDiscretization *other) const; bool areCompatibleForMul(const MEDCouplingTimeDiscretization *other) const; @@ -133,9 +141,13 @@ namespace ParaMEDMEM TypeOfTimeDiscretization getEnum() const { return DISCRETIZATION; } MEDCouplingTimeDiscretization *aggregate(const MEDCouplingTimeDiscretization *other) const; MEDCouplingTimeDiscretization *add(const MEDCouplingTimeDiscretization *other) const; + void addEqual(const MEDCouplingTimeDiscretization *other); MEDCouplingTimeDiscretization *substract(const MEDCouplingTimeDiscretization *other) const; + void substractEqual(const MEDCouplingTimeDiscretization *other); MEDCouplingTimeDiscretization *multiply(const MEDCouplingTimeDiscretization *other) const; + void multiplyEqual(const MEDCouplingTimeDiscretization *other); MEDCouplingTimeDiscretization *divide(const MEDCouplingTimeDiscretization *other) const; + void divideEqual(const MEDCouplingTimeDiscretization *other); bool isEqual(const MEDCouplingTimeDiscretization *other, double prec) const; bool areCompatible(const MEDCouplingTimeDiscretization *other) const; bool areCompatibleForMul(const MEDCouplingTimeDiscretization *other) const; @@ -184,9 +196,13 @@ namespace ParaMEDMEM TypeOfTimeDiscretization getEnum() const { return DISCRETIZATION; } MEDCouplingTimeDiscretization *aggregate(const MEDCouplingTimeDiscretization *other) const; MEDCouplingTimeDiscretization *add(const MEDCouplingTimeDiscretization *other) const; + void addEqual(const MEDCouplingTimeDiscretization *other); MEDCouplingTimeDiscretization *substract(const MEDCouplingTimeDiscretization *other) const; + void substractEqual(const MEDCouplingTimeDiscretization *other); MEDCouplingTimeDiscretization *multiply(const MEDCouplingTimeDiscretization *other) const; + void multiplyEqual(const MEDCouplingTimeDiscretization *other); MEDCouplingTimeDiscretization *divide(const MEDCouplingTimeDiscretization *other) const; + void divideEqual(const MEDCouplingTimeDiscretization *other); void setStartTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) { _start_time=time; _start_dt=dt; _start_it=it; } void setEndTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) { _end_time=time; _end_dt=dt; _end_it=it; } double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_start_dt; it=_start_it; return _start_time; } diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 360da8429..63183865c 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -25,6 +25,7 @@ #include "BBTree.txx" #include +#include #include #include @@ -604,6 +605,57 @@ void MEDCouplingUMesh::renumberNodes(const int *newNodeNumbers, int newNbOfNodes updateTime(); } +/*! + * This method renumbers cells of 'this' using the array specified by [old2NewBg;old2NewEnd) + * If std::distance(old2NewBg,old2NewEnd)!=this->getNumberOfCells() an INTERP_KERNEL::Exception will be thrown. + * + * If 'check' equals true the method will check that any elements in [old2NewBg;old2NewEnd) is unique ; if not + * an INTERP_KERNEL::Exception will be thrown. When 'check' equals true [old2NewBg;old2NewEnd) is not expected to + * be strictly in [0;this->getNumberOfCells()). + * + * If 'check' equals false the method will not check the content of [old2NewBg;old2NewEnd). + * To avoid any throw of SIGSEGV when 'check' equals false, the elements in [old2NewBg;old2NewEnd) should be unique and + * should be contained in[0;this->getNumberOfCells()). + */ +void MEDCouplingUMesh::renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) +{ + int nbCells=getNumberOfCells(); + const int *array=old2NewBg; + if(std::distance(old2NewBg,old2NewEnd)!=nbCells) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::renumberCells expected to take an array of size getNumberOfCells !"); + if(check) + array=DataArrayInt::checkAndPreparePermutation(old2NewBg,old2NewEnd); + // + const int *conn=_nodal_connec->getConstPointer(); + const int *connI=_nodal_connec_index->getConstPointer(); + DataArrayInt *newConn=DataArrayInt::New(); + newConn->alloc(_nodal_connec->getNumberOfTuples(),_nodal_connec->getNumberOfComponents()); + newConn->copyStringInfoFrom(*_nodal_connec); + DataArrayInt *newConnI=DataArrayInt::New(); + newConnI->alloc(_nodal_connec_index->getNumberOfTuples(),_nodal_connec_index->getNumberOfComponents()); + newConnI->copyStringInfoFrom(*_nodal_connec_index); + // + int *newC=newConn->getPointer(); + int *newCI=newConnI->getPointer(); + int loc=0; + newCI[0]=loc; + for(int i=0;idecrRef(); + newConnI->decrRef(); + if(check) + delete [] (int *)array; +} + /*! * Given a boundary box 'bbox' returns elements 'elems' contained in this 'bbox'. * Warning 'elems' is incremented during the call so if elems is not empty before call returned elements will be @@ -1025,6 +1077,80 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::buildOrthogonalField() const return ret; } +/*! + * This methods returns a vector newly created field on cells that represents the director vector of each 1D cell of this. + * This method is only callable on mesh with meshdim == 1 containing only SEG2. + */ +MEDCouplingFieldDouble *MEDCouplingUMesh::buildLinearField() const +{ + if(getMeshDimension()!=1) + throw INTERP_KERNEL::Exception("Expected a umesh with meshDim == 1 for buildLinearField !"); + if(_types.size()!=1 || *(_types.begin())!=INTERP_KERNEL::NORM_SEG2) + throw INTERP_KERNEL::Exception("Expected a umesh with only NORM_SEG2 type of elements for buildLinearField !"); + MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME); + DataArrayDouble *array=DataArrayDouble::New(); + int nbOfCells=getNumberOfCells(); + int spaceDim=getSpaceDimension(); + array->alloc(nbOfCells,spaceDim); + double *pt=array->getPointer(); + const double *coo=getCoords()->getConstPointer(); + std::vector conn; + conn.reserve(2); + for(int i=0;i()); + } + ret->setArray(array); + array->decrRef(); + ret->setMesh(this); + return ret; +} + +/*! + * This method is only callable on mesh with meshdim == 1 containing only SEG2 and spaceDim==3. + * This method projects this on the 3D line defined by (pt,v). This methods first checks that all SEG2 are along v vector. + * @param pt reference point of the line + * @param v normalized director vector of the line + * @param eps max precision before throwing an exception + * @param res output of size this->getNumberOfCells + */ +void MEDCouplingUMesh::project1D(const double *pt, const double *v, double eps, double *res) const +{ + if(getMeshDimension()!=1) + throw INTERP_KERNEL::Exception("Expected a umesh with meshDim == 1 for project1D !"); + if(_types.size()!=1 || *(_types.begin())!=INTERP_KERNEL::NORM_SEG2) + throw INTERP_KERNEL::Exception("Expected a umesh with only NORM_SEG2 type of elements for project1D !"); + if(getSpaceDimension()!=3) + throw INTERP_KERNEL::Exception("Expected a umesh with spaceDim==3 for project1D !"); + MEDCouplingFieldDouble *f=buildLinearField(); + const double *fPtr=f->getArray()->getConstPointer(); + double tmp[3]; + for(int i=0;i(tmp); + n1/=INTERP_KERNEL::norm<3>(tmp1); + if(n1>eps) + { + f->decrRef(); + throw INTERP_KERNEL::Exception("UMesh::Projection 1D failed !"); + } + } + const double *coo=getCoords()->getConstPointer(); + for(int i=0;i()); + std::transform(tmp,tmp+3,v,tmp,std::multiplies()); + res[i]=std::accumulate(tmp,tmp+3,0.); + } + f->decrRef(); +} + /*! * Returns a cell if any that contains the point located on 'pos' with precison eps. * If 'pos' is outside 'this' -1 is returned. If several cells contain this point the cell with the smallest id is returned. diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 8e71a19f9..0f58f803b 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -70,10 +70,13 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void findBoundaryNodes(std::vector& nodes) const; MEDCOUPLING_EXPORT MEDCouplingPointSet *buildBoundaryMesh(bool keepCoords) const; MEDCOUPLING_EXPORT void renumberNodes(const int *newNodeNumbers, int newNbOfNodes); + MEDCOUPLING_EXPORT void renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check); MEDCOUPLING_EXPORT void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems); MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(bool isAbs) const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildOrthogonalField() const; + MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildLinearField() const; + MEDCOUPLING_EXPORT void project1D(const double *pt, const double *v, double eps, double *res) const; MEDCOUPLING_EXPORT int getCellContainingPoint(const double *pos, double eps) const; MEDCOUPLING_EXPORT void getCellsContainingPoint(const double *pos, double eps, std::vector& elts) const; MEDCOUPLING_EXPORT void getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, std::vector& elts, std::vector& eltsIndex) const; diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx index e187bfcef..71175f08a 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx @@ -61,6 +61,7 @@ namespace ParaMEDMEM CPPUNIT_TEST( testApplyFunc2 ); CPPUNIT_TEST( testOperationsOnFields ); CPPUNIT_TEST( testOperationsOnFields2 ); + CPPUNIT_TEST( testOperationsOnFields3 ); CPPUNIT_TEST( testMergeNodesOnField ); CPPUNIT_TEST( testCheckConsecutiveCellTypes ); CPPUNIT_TEST( testBuildOrthogonalField ); @@ -70,6 +71,8 @@ namespace ParaMEDMEM CPPUNIT_TEST( testScale ); CPPUNIT_TEST( testTryToShareSameCoords ); CPPUNIT_TEST( testFindNodeOnPlane ); + CPPUNIT_TEST( testRenumberCells ); + CPPUNIT_TEST( testChangeSpaceDimension ); CPPUNIT_TEST( test2DInterpP0P0_1 ); CPPUNIT_TEST( test2DInterpP0P0PL_1 ); CPPUNIT_TEST( test2DInterpP0P0PL_2 ); @@ -160,6 +163,7 @@ namespace ParaMEDMEM void testApplyFunc2(); void testOperationsOnFields(); void testOperationsOnFields2(); + void testOperationsOnFields3(); void testMergeNodesOnField(); void testCheckConsecutiveCellTypes(); void testBuildOrthogonalField(); @@ -169,6 +173,8 @@ namespace ParaMEDMEM void testScale(); void testTryToShareSameCoords(); void testFindNodeOnPlane(); + void testRenumberCells(); + void testChangeSpaceDimension(); void test2DInterpP0P0_1(); void test2DInterpP0P0PL_1(); void test2DInterpP0P0PL_2(); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx index d24c3d4da..981e6093a 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx @@ -1487,6 +1487,43 @@ void MEDCouplingBasicsTest::testOperationsOnFields2() m->decrRef(); } +void MEDCouplingBasicsTest::testOperationsOnFields3() +{ + MEDCouplingUMesh *m=build3DSurfTargetMesh_1(); + MEDCouplingFieldDouble *f1=m->fillFromAnalytic(ON_NODES,1,"x+y+z"); + MEDCouplingFieldDouble *f2=m->fillFromAnalytic(ON_NODES,1,"a*a+b+c*c"); + (*f1)/=(*f2); + f1->checkCoherency(); + CPPUNIT_ASSERT(f1->getTypeOfField()==ON_NODES); + CPPUNIT_ASSERT(f1->getTimeDiscretization()==NO_TIME); + const double expected1[9]={-2.4999999999999991, 1.2162162162162162, 0.77868852459016391, + 0.7407407407407407, 1.129032258064516, 0.81632653061224492, + 0.86538461538461531, 1.0919540229885056, 0.84302325581395343}; + CPPUNIT_ASSERT_EQUAL(1,f1->getNumberOfComponents()); + CPPUNIT_ASSERT_EQUAL(9,f1->getNumberOfTuples()); + const double *val=f1->getArray()->getConstPointer(); + for(int i=0;i<9;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],val[i],1.e-12); + f1->decrRef(); + f2->decrRef(); + // + f1=m->buildOrthogonalField(); + f2=m->fillFromAnalytic(ON_CELLS,1,"x"); + (*f1)*=(*f2); + const double expected2[15]={-0.035355339059327376,0.,0.035355339059327376, 0.2592724864350674,0.,-0.2592724864350674, 0.37712361663282529,0.,-0.37712361663282529, -0.035355339059327376,0.,0.035355339059327376, 0.31819805153394637,0.,-0.31819805153394637}; + val=f1->getArray()->getConstPointer(); + for(int i=0;i<15;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],val[i],1.e-12); + f1->decrRef(); + // + f1=m->buildOrthogonalField(); + CPPUNIT_ASSERT_THROW((*f2)*=(*f1),INTERP_KERNEL::Exception); + f1->decrRef(); + f2->decrRef(); + // + m->decrRef(); +} + bool func4(const double *pt, double *res) { res[0]=pt[0]+pt[1]+pt[2]; @@ -1800,3 +1837,44 @@ void MEDCouplingBasicsTest::testFindNodeOnPlane() m3dSurf->decrRef(); mesh->decrRef(); } + +void MEDCouplingBasicsTest::testRenumberCells() +{ + MEDCouplingUMesh *m=build3DSurfTargetMesh_1(); + MEDCouplingUMesh *m2=build3DSurfTargetMesh_1(); + CPPUNIT_ASSERT(m->isEqual(m2,0)); + const int arr[5]={12,3,25,2,26}; + m->renumberCells(arr,arr+5,true); + CPPUNIT_ASSERT(!m->isEqual(m2,0)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_QUAD4,m->getTypeOfCell(0)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(1)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_QUAD4,m->getTypeOfCell(2)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(3)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_QUAD4,m->getTypeOfCell(4)); + const int arr2[5]={5,-1,-5,4,8}; + m->renumberCells(arr2,arr2+5,true); + CPPUNIT_ASSERT(m->isEqual(m2,0)); + m->decrRef(); + m2->decrRef(); +} + +void MEDCouplingBasicsTest::testChangeSpaceDimension() +{ + MEDCouplingUMesh *m1=build3DSurfTargetMesh_1(); + MEDCouplingUMesh *m2=build2DTargetMesh_1(); + // + CPPUNIT_ASSERT_EQUAL(3,m1->getSpaceDimension()); + m1->changeSpaceDimension(2); + CPPUNIT_ASSERT_EQUAL(2,m1->getSpaceDimension()); + m1->setName(m2->getName()); + CPPUNIT_ASSERT(m1->isEqual(m2,1e-12)); + m1->changeSpaceDimension(3); + CPPUNIT_ASSERT_EQUAL(3,m1->getSpaceDimension()); + const double expected[27]={-0.3,-0.3,0., 0.2,-0.3,0., 0.7,-0.3,0., -0.3,0.2,0., 0.2,0.2,0., 0.7,0.2,0., -0.3,0.7,0., 0.2,0.7,0., 0.7,0.7,0.}; + const double *val=m1->getCoords()->getConstPointer(); + for(int i=0;i<27;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[i],val[i],1e-14); + // + m1->decrRef(); + m2->decrRef(); +} diff --git a/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx b/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx index 29ef6031b..a13d013f6 100644 --- a/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx +++ b/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx @@ -26,6 +26,7 @@ #include "MEDCouplingBasicsTest.hxx" #include +#include using namespace ParaMEDMEM; @@ -696,6 +697,176 @@ void MEDCouplingRemapperTest::testExtruded() extT->decrRef(); } +void MEDCouplingRemapperTest::testExtruded2() +{ + MEDCouplingUMesh *meshN,*meshTT,*meshTF; + MEDCouplingBasicsTest::build3DExtrudedUMesh_2(meshN,meshTT,meshTF); + std::vector n; + double pt[3]={300.,300.,0.}; + double v[3]={0.,0.,2.}; + meshN->findNodesOnPlane(pt,v,1e-12,n); + MEDCouplingUMesh *meshN2D=(MEDCouplingUMesh *)meshN->buildFacePartOfMySelfNode(&n[0],&n[0]+n.size(),true); + n.clear(); + bool b=false; + DataArrayInt *da=meshTT->mergeNodes(1e-12,b); + CPPUNIT_ASSERT(b); + da->decrRef(); + meshTT->findNodesOnPlane(pt,v,1e-12,n); + MEDCouplingUMesh *meshTT2D=(MEDCouplingUMesh *)meshTT->buildFacePartOfMySelfNode(&n[0],&n[0]+n.size(),true); + n.clear(); + meshTF->findNodesOnPlane(pt,v,1e-12,n); + MEDCouplingUMesh *meshTF2D=(MEDCouplingUMesh *)meshTF->buildFacePartOfMySelfNode(&n[0],&n[0]+n.size(),true); + n.clear(); + // + MEDCouplingExtrudedMesh *meshNE=MEDCouplingExtrudedMesh::New(meshN,meshN2D,0); + MEDCouplingExtrudedMesh *meshTTE=MEDCouplingExtrudedMesh::New(meshTT,meshTT2D,0); + MEDCouplingExtrudedMesh *meshTFE=MEDCouplingExtrudedMesh::New(meshTF,meshTF2D,0); + // + MEDCouplingRemapper remapper; + remapper.setPrecision(1e-12); + remapper.setIntersectionType(INTERP_KERNEL::Geometric2D); + CPPUNIT_ASSERT_EQUAL(1,remapper.prepare(meshNE,meshTTE,"P0P0")); + MEDCouplingFieldDouble *srcField=MEDCouplingFieldDouble::New(ON_CELLS); + srcField->setNature(IntegralGlobConstraint); + srcField->setMesh(meshNE); + DataArrayDouble *array=DataArrayDouble::New(); + array->alloc(meshNE->getNumberOfCells(),1); + srcField->setArray(array); + double vals1[40]={ + 1000.,1000.,1020.,1030.,1040.,1000.,1000.,1070.,1080.,1090.,1000.,1000.,1120.,1130.,1140.,1000.,1000.,1170.,1180.,1190., + 2000.,2000.,2020.,2030.,2040.,2000.,2000.,2070.,2080.,2090.,2000.,2000.,2120.,2130.,2140.,2000.,2000.,2170.,2180.,2190., + }; + CPPUNIT_ASSERT_EQUAL((int)(sizeof(vals1)/sizeof(double)),meshNE->getNumberOfCells()); + std::copy(vals1,vals1+meshNE->getNumberOfCells(),array->getPointer()); + array->decrRef(); + MEDCouplingFieldDouble *trgField=remapper.transferField(srcField,4.220173); + double expected1[200]={ + 800.,800.,800.,800.,800.,800.,800.,800.,800.,800.,1600.,1600.,1600.,1600.,1600.,1600.,1600.,1600.,1600.,1600., + 102.,102.,102.,102.,102.,102.,102.,102.,102.,102.,202.,202.,202.,202.,202.,202.,202.,202.,202.,202., + 103.,103.,103.,103.,103.,103.,103.,103.,103.,103.,203.,203.,203.,203.,203.,203.,203.,203.,203.,203., + 104.,104.,104.,104.,104.,104.,104.,104.,104.,104.,204.,204.,204.,204.,204.,204.,204.,204.,204.,204., + 219.,219.,219.,219.,219.,219.,219.,219.,219.,219.,419.,419.,419.,419.,419.,419.,419.,419.,419.,419., + 221.,221.,221.,221.,221.,221.,221.,221.,221.,221.,421.,421.,421.,421.,421.,421.,421.,421.,421.,421., + 223.,223.,223.,223.,223.,223.,223.,223.,223.,223.,423.,423.,423.,423.,423.,423.,423.,423.,423.,423., + 117.,117.,117.,117.,117.,117.,117.,117.,117.,117.,217.,217.,217.,217.,217.,217.,217.,217.,217.,217., + 118.,118.,118.,118.,118.,118.,118.,118.,118.,118.,218.,218.,218.,218.,218.,218.,218.,218.,218.,218., + 119.,119.,119.,119.,119.,119.,119.,119.,119.,119.,219.,219.,219.,219.,219.,219.,219.,219.,219.,219. + }; + for(int i=0;i<200;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],trgField->getArray()->getConstPointer()[i],1e-3);//1e-3 precision due to non coincidence in 1D mesh + CPPUNIT_ASSERT_DOUBLES_EQUAL(std::accumulate(expected1,expected1+200,0.),std::accumulate(vals1,vals1+40,0.),1e-10); + CPPUNIT_ASSERT_DOUBLES_EQUAL(std::accumulate(expected1,expected1+200,0.),std::accumulate(trgField->getArray()->getConstPointer(),trgField->getArray()->getConstPointer()+200,0.),1e-10); + trgField->decrRef(); + // + CPPUNIT_ASSERT_EQUAL(1,remapper.prepare(meshNE,meshTFE,"P0P0")); + trgField=remapper.transferField(srcField,4.220173); + double expected2[340]={25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, + 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, + 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, + 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, + 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, + 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, + 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, 29.75, 25.5, 51.25, 51.75, 26., 79., 158.75, 160.25, 80.5, 85.25, 171.25, 172.75, 86.75, 29.25, 58.75, 59.25, + 29.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., + 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, + 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, + 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, + 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, + 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, + 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 50.5, 101.25, 101.75, 51., 154., 308.75, 310.25, 155.5, 160.25, 321.25, + 322.75, 161.75, 54.25, 108.75, 109.25, 54.75, 800., 800., 800., 800., 800., 800., 800., 800., 800., 800., 1600., 1600., 1600., 1600., 1600., 1600., 1600., + 1600., 1600., 1600.}; + for(int i=0;i<340;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],trgField->getArray()->getConstPointer()[i],1e-3);//1e-3 precision due to non coincidence in 1D mesh + CPPUNIT_ASSERT_DOUBLES_EQUAL(std::accumulate(expected2,expected2+340,0.),std::accumulate(vals1,vals1+40,0.),1e-10); + CPPUNIT_ASSERT_DOUBLES_EQUAL(std::accumulate(expected2,expected2+340,0.),std::accumulate(trgField->getArray()->getConstPointer(),trgField->getArray()->getConstPointer()+340,0.),1e-10); + trgField->decrRef(); + srcField->decrRef(); + // + double vals2[200]={ + 100., 200., 300., 400., 500., 600., 700., 800., 900., 1000., 1100., 1200., 1300., 1400., 1500., 1600., 1700., 1800., 1900., 2000, + 101., 201., 301., 401., 501., 601., 701., 801., 901., 1001., 1101., 1201., 1301., 1401., 1501., 1601., 1701., 1801., 1901., 2001, + 102., 202., 302., 402., 502., 602., 702., 802., 902., 1002., 1102., 1202., 1302., 1402., 1502., 1602., 1702., 1802., 1902., 2002, + 103., 203., 303., 403., 503., 603., 703., 803., 903., 1003., 1103., 1203., 1303., 1403., 1503., 1603., 1703., 1803., 1903., 2003, + 104., 204., 304., 404., 504., 604., 704., 804., 904., 1004., 1104., 1204., 1304., 1404., 1504., 1604., 1704., 1804., 1904., 2004, + 105., 205., 305., 405., 505., 605., 705., 805., 905., 1005., 1105., 1205., 1305., 1405., 1505., 1605., 1705., 1805., 1905., 2005, + 106., 206., 306., 406., 506., 606., 706., 806., 906., 1006., 1106., 1206., 1306., 1406., 1506., 1606., 1706., 1806., 1906., 2006, + 107., 207., 307., 407., 507., 607., 707., 807., 907., 1007., 1107., 1207., 1307., 1407., 1507., 1607., 1707., 1807., 1907., 2007, + 108., 208., 308., 408., 508., 608., 708., 808., 908., 1008., 1108., 1208., 1308., 1408., 1508., 1608., 1708., 1808., 1908., 2008, + 109., 209., 309., 409., 509., 609., 709., 809., 909., 1009., 1109., 1209., 1309., 1409., 1509., 1609., 1709., 1809., 1909., 2009. + }; + CPPUNIT_ASSERT_EQUAL(1,remapper.prepare(meshNE,meshTTE,"P0P0")); + trgField=MEDCouplingFieldDouble::New(ON_CELLS); + trgField->setNature(ConservativeVolumic); + trgField->setMesh(meshTTE); + array=DataArrayDouble::New(); + array->alloc(meshTTE->getNumberOfCells(),1); + trgField->setArray(array); + std::copy(vals2,vals2+meshTTE->getNumberOfCells(),array->getPointer()); + array->decrRef(); + srcField=remapper.reverseTransferField(trgField,4.220173); + double expected3[40]={ + 550.,550.,551.,552.,553.,550.,550.,554.,555.,556.,550.,550.,554.,555.,556.,550.,550.,557.,558.,559., + 1550.,1550.,1551.,1552.,1553.,1550.,1550.,1554.,1555.,1556.,1550.,1550.,1554.,1555.,1556.,1550.,1550.,1557.,1558.,1559. + }; + for(int i=0;i<40;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected3[i],srcField->getArray()->getConstPointer()[i],1e-3);//1e-3 precision due to non coincidence in 1D mesh + srcField->decrRef(); + trgField->decrRef(); + // + double vals3[340]={ + 100., 101., 102., 103., 104., 105., 106., 107., 108., 109., 110., 111., 112., 113., 114., 115., + 200., 201., 202., 203., 204., 205., 206., 207., 208., 209., 210., 211., 212., 213., 214., 215., + 300., 301., 302., 303., 304., 305., 306., 307., 308., 309., 310., 311., 312., 313., 314., 315., + 400., 401., 402., 403., 404., 405., 406., 407., 408., 409., 410., 411., 412., 413., 414., 415., + 500., 501., 502., 503., 504., 505., 506., 507., 508., 509., 510., 511., 512., 513., 514., 515., + 600., 601., 602., 603., 604., 605., 606., 607., 608., 609., 610., 611., 612., 613., 614., 615., + 700., 701., 702., 703., 704., 705., 706., 707., 708., 709., 710., 711., 712., 713., 714., 715., + 800., 801., 802., 803., 804., 805., 806., 807., 808., 809., 810., 811., 812., 813., 814., 815., + 900., 901., 902., 903., 904., 905., 906., 907., 908., 909., 910., 911., 912., 913., 914., 915., + 1000., 1001., 1002., 1003., 1004., 1005., 1006., 1007., 1008., 1009., 1010., 1011., 1012., 1013., 1014., 1015., + 1100., 1101., 1102., 1103., 1104., 1105., 1106., 1107., 1108., 1109., 1110., 1111., 1112., 1113., 1114., 1115., + 1200., 1201., 1202., 1203., 1204., 1205., 1206., 1207., 1208., 1209., 1210., 1211., 1212., 1213., 1214., 1215., + 1300., 1301., 1302., 1303., 1304., 1305., 1306., 1307., 1308., 1309., 1310., 1311., 1312., 1313., 1314., 1315., + 1400., 1401., 1402., 1403., 1404., 1405., 1406., 1407., 1408., 1409., 1410., 1411., 1412., 1413., 1414., 1415., + 1500., 1501., 1502., 1503., 1504., 1505., 1506., 1507., 1508., 1509., 1510., 1511., 1512., 1513., 1514., 1515., + 1600., 1601., 1602., 1603., 1604., 1605., 1606., 1607., 1608., 1609., 1610., 1611., 1612., 1613., 1614., 1615., + 1700., 1701., 1702., 1703., 1704., 1705., 1706., 1707., 1708., 1709., 1710., 1711., 1712., 1713., 1714., 1715., + 1800., 1801., 1802., 1803., 1804., 1805., 1806., 1807., 1808., 1809., 1810., 1811., 1812., 1813., 1814., 1815., + 1900., 1901., 1902., 1903., 1904., 1905., 1906., 1907., 1908., 1909., 1910., 1911., 1912., 1913., 1914., 1915., + 2000., 2001., 2002., 2003., 2004., 2005., 2006., 2007., 2008., 2009., 2010., 2011., 2012., 2013., 2014., 2015., + 116.,216.,316.,416.,516.,616.,716.,816.,916.,1016.,1116.,1216.,1316.,1416.,1516.,1616.,1716.,1816.,1916.,2016. + }; + CPPUNIT_ASSERT_EQUAL(1,remapper.prepare(meshNE,meshTFE,"P0P0")); + trgField=MEDCouplingFieldDouble::New(ON_CELLS); + trgField->setNature(ConservativeVolumic); + trgField->setMesh(meshTFE); + array=DataArrayDouble::New(); + array->alloc(meshTFE->getNumberOfCells(),1); + trgField->setArray(array); + std::copy(vals3,vals3+meshTFE->getNumberOfCells(),array->getPointer()); + array->decrRef(); + srcField=remapper.reverseTransferField(trgField,4.220173); + double expected4[40]={ + 566.,566.,552.5,553.5,554.5,566.,566.,554.5,555.5,556.5,566.,566.,558.5,559.5,560.5,566.,566.,560.5,561.5,562.5, + 1566.,1566.,1552.5,1553.5,1554.5,1566.,1566.,1554.5,1555.5,1556.5,1566.,1566.,1558.5,1559.5,1560.5,1566.,1566.,1560.5,1561.5,1562.5 + }; + for(int i=0;i<40;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected4[i],srcField->getArray()->getConstPointer()[i],1e-3);//1e-3 precision due to non coincidence in 1D mesh + srcField->decrRef(); + trgField->decrRef(); + // + meshN2D->decrRef(); + meshTT2D->decrRef(); + meshTF2D->decrRef(); + meshNE->decrRef(); + meshTTE->decrRef(); + meshTFE->decrRef(); + meshN->decrRef(); + meshTT->decrRef(); + meshTF->decrRef(); +} + MEDCouplingUMesh *MEDCouplingRemapperTest::build1DTargetMesh_2() { double targetCoords[20]={ diff --git a/src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx b/src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx index 5bcebfb89..29cf7cbcd 100644 --- a/src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingRemapperTest.hxx @@ -37,6 +37,7 @@ namespace ParaMEDMEM CPPUNIT_TEST( testMultiDimCombi ); CPPUNIT_TEST( testNatureOfField ); CPPUNIT_TEST( testExtruded ); + CPPUNIT_TEST( testExtruded2 ); CPPUNIT_TEST_SUITE_END(); public: void test2DInterpP0P0_1(); @@ -45,6 +46,7 @@ namespace ParaMEDMEM void testMultiDimCombi(); void testNatureOfField(); void testExtruded(); + void testExtruded2(); private: static MEDCouplingUMesh *build1DTargetMesh_2(); static MEDCouplingUMesh *build2DTargetMesh_3();