From 47a73b15cbb7c005984d36dc9bf9e0d2bd348a42 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Tue, 23 Aug 2022 17:05:35 +0200 Subject: [PATCH] Remapper takes advantage of implementation of MEDCouplingFieldDiscretizationOnNodesFE::getValueOn --- ...EDCouplingFieldDiscretizationOnNodesFE.cxx | 88 ++++++++++--------- ...EDCouplingFieldDiscretizationOnNodesFE.hxx | 5 ++ src/MEDCoupling/MEDCouplingRemapper.cxx | 76 ++++++++-------- src/MEDCoupling/MEDCouplingRemapper.hxx | 1 + 4 files changed, 94 insertions(+), 76 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx index f769818e9..4e0138911 100755 --- a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx @@ -146,16 +146,46 @@ bool IsInside3D(const MEDCouplingGaussLocalization& gl, const std::vector res2( this->getValueOnMulti(arr,mesh,loc,1) ); + MCAuto res2( this->getValueOnMulti(arr,mesh,ptsCoo,1) ); std::copy(res2->begin(),res2->end(),res); } -/*void GetValueOnMulti() +void MEDCouplingFieldDiscretizationOnNodesFE::GetRefCoordOfListOf3DPtsIn3D(const MEDCouplingUMesh *umesh, const double *ptsCoo, mcIdType nbOfPts, + std::function&)> customFunc) { - -}*/ + const double *coordsOfMesh( umesh->getCoords()->begin() ); + MEDCouplingNormalizedUnstructuredMesh<3,3> mesh_wrapper(umesh); + BBTreeStandAlone<3,mcIdType> tree( INTERP_KERNEL::BuildBBTree( mesh_wrapper ) ); + for(mcIdType iPt = 0 ; iPt < nbOfPts ; ++iPt) + { + std::vector elems; + tree.getElementsAroundPoint(ptsCoo+3*iPt,elems); + bool found(false); + for(auto cellId = elems.cbegin() ; cellId != elems.cend() && !found ; ++cellId) + { + INTERP_KERNEL::NormalizedCellType gt( umesh->getTypeOfCell(*cellId) ); + std::vector conn; + umesh->getNodeIdsOfCell(*cellId,conn); + MCAuto refCoo( MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(gt) ); + std::vector refCooCpp(refCoo->begin(),refCoo->end()); + std::vector gsCoo(ptsCoo + iPt*3,ptsCoo + (iPt+1)*3); + MEDCouplingGaussLocalization gl(gt,refCooCpp,{0,0,0},{1.}); + std::vector ptsInCell; ptsInCell.reserve(conn.size()*gl.getDimension()); + std::for_each( conn.cbegin(), conn.cend(), [coordsOfMesh,&ptsInCell](mcIdType c) { ptsInCell.insert(ptsInCell.end(),coordsOfMesh+c*3,coordsOfMesh+(c+1)*3); } ); + std::vector locInRef(3); + if( IsInside3D(gl,ptsInCell,gsCoo.data(),locInRef.data()) ) + { + gl.setGaussCoords(locInRef); + customFunc(gl,conn); + found = true; + } + } + if(!found) + THROW_IK_EXCEPTION("getValueOnMulti on MEDCouplingFieldDiscretizationOnNodesFE : fail to locate point #" << iPt << " X=" << ptsCoo[0] << " Y=" << ptsCoo[1] << " Z=" << ptsCoo[2] << " !"); + } +} DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfTargetPoints) const { @@ -176,46 +206,24 @@ DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const std::size_t nbCompo( arr->getNumberOfComponents() ); MCAuto ret(DataArrayDouble::New()); ret->alloc(nbOfTargetPoints,nbCompo); - const double *coordsOfMesh( umesh->getCoords()->begin() ); - MEDCouplingNormalizedUnstructuredMesh<3,3> mesh_wrapper(umesh); - BBTreeStandAlone<3,mcIdType> tree( INTERP_KERNEL::BuildBBTree( mesh_wrapper ) ); double *res( ret->getPointer() ); - for(mcIdType iPt = 0 ; iPt < nbOfTargetPoints ; ++iPt, res+=nbCompo ) + + auto arrayFeeder = [arr,&res,nbCompo](const MEDCouplingGaussLocalization& gl, const std::vector& conn) { - std::vector elems; - tree.getElementsAroundPoint(loc+3*iPt,elems); - bool found(false); - for(auto cellId = elems.cbegin() ; cellId != elems.cend() && !found ; ++cellId) + MCAuto resVector( gl.getShapeFunctionValues() ); { - INTERP_KERNEL::NormalizedCellType gt( umesh->getTypeOfCell(*cellId) ); - std::vector conn; - umesh->getNodeIdsOfCell(*cellId,conn); - MCAuto refCoo( MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(gt) ); - std::vector refCooCpp(refCoo->begin(),refCoo->end()); - std::vector gsCoo(loc + iPt*3,loc + (iPt+1)*3); - MEDCouplingGaussLocalization gl(gt,refCooCpp,{0,0,0},{1.}); - std::vector ptsInCell; ptsInCell.reserve(conn.size()*gl.getDimension()); - std::for_each( conn.cbegin(), conn.cend(), [coordsOfMesh,&ptsInCell](mcIdType c) { ptsInCell.insert(ptsInCell.end(),coordsOfMesh+c*3,coordsOfMesh+(c+1)*3); } ); - std::vector locInRef(3); - if( IsInside3D(gl,ptsInCell,gsCoo.data(),locInRef.data()) ) - { - gl.setGaussCoords(locInRef); - MCAuto resVector( gl.getShapeFunctionValues() ); + std::for_each(res,res+nbCompo,[](double& v) { v = 0.0; }); + for(std::size_t iComp = 0 ; iComp < nbCompo ; ++iComp) + for(int iPt = 0 ; iPt < gl.getNumberOfPtsInRefCell(); ++iPt) { - std::for_each(res,res+nbCompo,[](double& v) { v = 0.0; }); - for(std::size_t iComp = 0 ; iComp < nbCompo ; ++iComp) - for(int iPt = 0 ; iPt < gl.getNumberOfPtsInRefCell(); ++iPt) - { - { - res[iComp] += resVector->getIJ(0,iPt) * arr->getIJ(conn[iPt],iComp); - } - } + { + res[iComp] += resVector->getIJ(0,iPt) * arr->getIJ(conn[iPt],iComp); + } } - found = true; - } + res += nbCompo; } - if(!found) - THROW_IK_EXCEPTION("getValueOnMulti on MEDCouplingFieldDiscretizationOnNodesFE : fail to locate point #" << iPt << " X=" << loc[0] << " Y=" << loc[1] << " Z=" << loc[2] << " !"); - } + }; + + GetRefCoordOfListOf3DPtsIn3D(umesh,loc,nbOfTargetPoints,arrayFeeder); return ret.retn(); } diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.hxx index 04559381b..3aa28072b 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.hxx @@ -22,6 +22,8 @@ #include "MEDCouplingFieldDiscretization.hxx" +#include + namespace MEDCoupling { /*! @@ -42,6 +44,9 @@ namespace MEDCoupling MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const override; MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const override; MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const override; + public: + MEDCOUPLING_EXPORT static void GetRefCoordOfListOf3DPtsIn3D(const MEDCouplingUMesh *umesh, const double *ptsCoo, mcIdType nbOfPts, + std::function&)> customFunc); public: static const char REPR[]; static constexpr TypeOfField TYPE = ON_NODES_FE; diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index c85f44ccc..231f29bba 100755 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -27,6 +27,7 @@ #include "MEDCouplingCMesh.hxx" #include "MEDCouplingNormalizedUnstructuredMesh.txx" #include "MEDCouplingNormalizedCartesianMesh.txx" +#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx" #include "Interpolation1D.txx" #include "Interpolation2DCurve.hxx" @@ -123,10 +124,7 @@ void MEDCouplingRemapper::setCrudeMatrixEx(const MEDCouplingFieldTemplate *src, } } _matrix=m; - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(srcNbElem); + synchronizeSizeOfSideMatricesAfterMatrixComputation(srcNbElem); } int MEDCouplingRemapper::prepareInterpKernelOnly() @@ -691,11 +689,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() } else throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension"); - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(nbCols); - declareAsNew(); + synchronizeSizeOfSideMatricesAfterMatrixComputation(nbCols); return 1; } @@ -729,11 +723,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyEE() buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D, target_mesh->getMesh3DIds()->getConstPointer()); // - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(nbCols2D*nbCols1D); - declareAsNew(); + synchronizeSizeOfSideMatricesAfterMatrixComputation(nbCols2D*nbCols1D); return 1; } @@ -785,11 +775,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUC() ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix); nullifiedTinyCoeffInCrudeMatrixAbs(0.); // - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(src_mesh->getNumberOfCells()); - declareAsNew(); + synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells()); return 1; } @@ -839,11 +825,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyCU() } nullifiedTinyCoeffInCrudeMatrixAbs(0.); // - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(src_mesh->getNumberOfCells()); - declareAsNew(); + synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells()); return 1; } @@ -892,11 +874,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyCC() } nullifiedTinyCoeffInCrudeMatrixAbs(0.); // - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(src_mesh->getNumberOfCells()); - declareAsNew(); + synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells()); return 1; } @@ -951,11 +929,7 @@ int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss() for(const mcIdType *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++) _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.; } - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples()); - declareAsNew(); + synchronizeSizeOfSideMatricesAfterMatrixComputation( srcLoc->getNumberOfTuples() ); return 1; } @@ -969,8 +943,29 @@ int MEDCouplingRemapper::prepareNotInterpKernelOnlyFEFE() THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only spacedim 3 supported for target !") if(_src_ft->getMesh()->getSpaceDimension() != 3) THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only spacedim 3 supported for source !") - - declareAsNew(); + const MEDCouplingUMesh *srcUMesh( dynamic_cast(_src_ft->getMesh()) ); + const MEDCouplingPointSet *trgMesh( dynamic_cast(_target_ft->getMesh()) ); + if( !srcUMesh ) + THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only 3D UMesh supported as source !"); + if( !trgMesh ) + THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only 3D PointSet mesh supported as target !"); + + _matrix.clear(); + _matrix.resize(trgMesh->getNumberOfNodes()); + mcIdType rowId(0); + + auto matrixFeeder = [this,&rowId](const MEDCouplingGaussLocalization& gl, const std::vector& conn) + { + auto& row = this->_matrix[rowId++]; + MCAuto resVector( gl.getShapeFunctionValues() ); + for(int iPt = 0 ; iPt < gl.getNumberOfPtsInRefCell(); ++iPt) + { + row[ conn[iPt] ] = resVector->getIJ(0,iPt); + } + }; + + MEDCouplingFieldDiscretizationOnNodesFE::GetRefCoordOfListOf3DPtsIn3D(srcUMesh,trgMesh->getCoords()->begin(),trgMesh->getNumberOfNodes(),matrixFeeder); + synchronizeSizeOfSideMatricesAfterMatrixComputation( srcUMesh->getNumberOfNodes() ); return 1; } @@ -1048,6 +1043,15 @@ void MEDCouplingRemapper::checkPrepare() const throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !"); } +void MEDCouplingRemapper::synchronizeSizeOfSideMatricesAfterMatrixComputation(mcIdType nbOfColsInMatrix) +{ + _deno_multiply.clear(); + _deno_multiply.resize(_matrix.size()); + _deno_reverse_multiply.clear(); + _deno_reverse_multiply.resize(nbOfColsInMatrix); + declareAsNew(); +} + /*! * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft. * This method returns 3 information (2 in output parameters and 1 in return). diff --git a/src/MEDCoupling/MEDCouplingRemapper.hxx b/src/MEDCoupling/MEDCouplingRemapper.hxx index 2acdafea1..23bf893b3 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.hxx +++ b/src/MEDCoupling/MEDCouplingRemapper.hxx @@ -96,6 +96,7 @@ namespace MEDCoupling bool isInterpKernelOnlyOrNotOnly() const; void updateTime() const; void checkPrepare() const; + void synchronizeSizeOfSideMatricesAfterMatrixComputation(mcIdType nbOfColsInMatrix); std::string checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const; void releaseData(bool matrixSuppression); void restartUsing(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target); -- 2.39.2