From 9143c19cb83a7a32cceca489743c494f46802dd9 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Mon, 22 Aug 2022 17:54:55 +0200 Subject: [PATCH] WIP --- ...EDCouplingFieldDiscretizationOnNodesFE.cxx | 80 ++++++++++++++----- src/MEDCoupling/MEDCouplingMemArray.txx | 2 +- src/MEDCoupling/MEDCouplingRemapper.cxx | 21 ++++- src/MEDCoupling/MEDCouplingRemapper.hxx | 1 + 4 files changed, 84 insertions(+), 20 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx index 4494bec34..db791fadb 100755 --- a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx @@ -81,25 +81,31 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationOnNodesFE::getMeasureField throw INTERP_KERNEL::Exception("getMeasureField on MEDCouplingFieldDiscretizationOnNodesFE : not implemented yet !"); } +bool IsInside3D(const MEDCouplingGaussLocalization& gl, const std::vector& ptsInCell, const double locInReal[3], double locInRef[3]) +{ + locInRef[0] = 0; locInRef[1] = 0; locInRef[2] = 0; + return true; +} + void MEDCouplingFieldDiscretizationOnNodesFE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const +{ + MCAuto res2( this->getValueOnMulti(arr,mesh,loc,1) ); + std::copy(res2->begin(),res2->end(),res); +} + +/*void GetValueOnMulti() +{ + +}*/ + +DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfTargetPoints) const { const MEDCouplingUMesh *umesh( dynamic_cast(mesh) ); if( !umesh ) THROW_IK_EXCEPTION("getValueOn : not implemented yet for type != MEDCouplingUMesh !"); if(umesh->getSpaceDimension() != 3 || umesh->getMeshDimension() != 3) THROW_IK_EXCEPTION("getValueOn : implemented only for meshes with spacedim == 3 and meshdim == 3 !"); - MEDCouplingNormalizedUnstructuredMesh<3,3> mesh_wrapper(umesh); - BBTreeStandAlone<3,mcIdType> tree( INTERP_KERNEL::BuildBBTree( mesh_wrapper ) ); - std::vector elems; - tree.getElementsAroundPoint(loc,elems); - - /*MCAuto res2(this->getValueOnMulti(arr,mesh,loc,1)); - std::copy(res2->begin(),res2->end(),res);*/ - throw INTERP_KERNEL::Exception("getValueOn on MEDCouplingFieldDiscretizationOnNodesFE : not implemented yet !"); -} - -DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfTargetPoints) const -{ + umesh->checkConsistency(); if(!arr || !arr->isAllocated()) throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti : input array is null or not allocated !"); mcIdType nbOfRows=getNumberOfMeshPlaces(mesh); @@ -107,12 +113,50 @@ DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const { THROW_IK_EXCEPTION( "MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti : input array does not have correct number of tuples ! Excepted " << nbOfRows << " having " << arr->getNumberOfTuples() << " !") } - throw INTERP_KERNEL::Exception("getValueOnMulti on MEDCouplingFieldDiscretizationOnNodesFE : not implemented yet !"); - /*mcIdType nbCols(-1); - std::size_t nbCompo=arr->getNumberOfComponents(); - MCAuto m(computeEvaluationMatrixOnGivenPts(mesh,loc,nbOfTargetPoints,nbCols)); + umesh->checkConsistency(); + std::size_t nbCompo( arr->getNumberOfComponents() ); MCAuto ret(DataArrayDouble::New()); ret->alloc(nbOfTargetPoints,nbCompo); - INTERP_KERNEL::matrixProduct(m->begin(),nbOfTargetPoints,nbCols,arr->begin(),nbOfRows,ToIdType(nbCompo),ret->getPointer()); - return ret.retn();*/ + 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 ) + { + std::vector elems; + tree.getElementsAroundPoint(loc+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(loc,loc+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) + { + { + res[iComp] += resVector->getIJ(0,iPt) * arr->getIJ(conn[iPt],iComp); + } + } + } + found = true; + } + } + if(!found) + THROW_IK_EXCEPTION("getValueOnMulti on MEDCouplingFieldDiscretizationOnNodesFE : fail to locate point #" << iPt << " X=" << loc[0] << " Y=" << loc[1] << " Z=" << loc[2] << " !"); + } + return ret.retn(); } diff --git a/src/MEDCoupling/MEDCouplingMemArray.txx b/src/MEDCoupling/MEDCouplingMemArray.txx index 9cb18e6fc..16a07683f 100755 --- a/src/MEDCoupling/MEDCouplingMemArray.txx +++ b/src/MEDCoupling/MEDCouplingMemArray.txx @@ -5866,7 +5866,7 @@ struct NotInRange throw INTERP_KERNEL::Exception("DataArrayInt::deltaShiftIndex : only single component allowed !"); std::size_t nbOfElements=this->getNumberOfTuples(); if(nbOfElements<2) - throw INTERP_KERNEL::Exception("DataArrayInt::deltaShiftIndex : 1 tuple at least must be present in 'this' !"); + throw INTERP_KERNEL::Exception("DataArrayInt::deltaShiftIndex : 2 tuples at least must be present in 'this' !"); const T *ptr=this->getConstPointer(); DataArrayType *ret=DataArrayType::New(); ret->alloc(nbOfElements-1,1); diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index 8f620ffe6..c85f44ccc 100755 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -181,6 +181,8 @@ int MEDCouplingRemapper::prepareNotInterpKernelOnly() { case 0: return prepareNotInterpKernelOnlyGaussGauss(); + case 1: + return prepareNotInterpKernelOnlyFEFE(); default: { std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !"; @@ -957,6 +959,21 @@ int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss() return 1; } +int MEDCouplingRemapper::prepareNotInterpKernelOnlyFEFE() +{ + if(getIntersectionType()!=INTERP_KERNEL::PointLocator) + throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareNotInterpKernelOnlyFEFE : The intersection type is not supported ! Only PointLocator is supported for FE->FE interpolation ! Please invoke setIntersectionType(PointLocator) on the MEDCouplingRemapper instance !"); + MCAuto trgLoc=_target_ft->getLocalizationOfDiscr(); + mcIdType trgSpaceDim=ToIdType(trgLoc->getNumberOfComponents()); + if(trgSpaceDim!=3) + 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(); + return 1; +} + /*! * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods. * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method. @@ -965,9 +982,11 @@ int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel { if(method=="GAUSSGAUSS") return 0; + if(method=="FEFE") + return 1; std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : "; oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method."; - oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !"; + oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS FEFE !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); } diff --git a/src/MEDCoupling/MEDCouplingRemapper.hxx b/src/MEDCoupling/MEDCouplingRemapper.hxx index 6966f6575..2acdafea1 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.hxx +++ b/src/MEDCoupling/MEDCouplingRemapper.hxx @@ -89,6 +89,7 @@ namespace MEDCoupling // int prepareNotInterpKernelOnly(); int prepareNotInterpKernelOnlyGaussGauss(); + int prepareNotInterpKernelOnlyFEFE(); // static int CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method); // -- 2.39.2