From: ageay Date: Mon, 19 Apr 2010 08:55:54 +0000 (+0000) Subject: *** empty log message *** X-Git-Tag: V5_1_main_FINAL~123 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=c6e0d0d939cbc2a1fca0bee9673373dee530ec2c;p=tools%2Fmedcoupling.git *** empty log message *** --- diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx index f0f093140..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 @@ -341,10 +344,23 @@ int MEDCouplingExtrudedMesh::findCorrespCellByNodalConn(const std::vector& void MEDCouplingExtrudedMesh::project1DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, MEDCouplingUMesh *&m1r, MEDCouplingUMesh *&m2r, double *v) throw(INTERP_KERNEL::Exception) { + 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/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index acfb5511c..f7306bc36 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" @@ -295,19 +296,21 @@ int MEDCouplingRemapper::prepareEE(const char *method) 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); + INTERP_KERNEL::Interpolation2D 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; double v[3]; MEDCouplingExtrudedMesh::project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v); - MEDCouplingNormalizedUnstructuredMesh<2,1> s1DWrapper(s1D); - MEDCouplingNormalizedUnstructuredMesh<2,1> t1DWrapper(t1D); + 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()); @@ -523,3 +526,29 @@ 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); + } + } + } + } +} diff --git a/src/MEDCoupling/MEDCouplingRemapper.hxx b/src/MEDCoupling/MEDCouplingRemapper.hxx index 461ce1711..5ecef0758 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.hxx +++ b/src/MEDCoupling/MEDCouplingRemapper.hxx @@ -58,6 +58,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/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index c39d36250..63183865c 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -25,6 +25,7 @@ #include "BBTree.txx" #include +#include #include #include @@ -1101,12 +1102,55 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::buildLinearField() const getNodeIdsOfCell(i,conn); pt=std::transform(coo+conn[1]*spaceDim,coo+(conn[1]+1)*spaceDim,coo+conn[0]*spaceDim,pt,std::minus()); } - array->alloc(nbOfCells,spaceDim); + 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 1195a327d..0f58f803b 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -76,6 +76,7 @@ namespace ParaMEDMEM 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;