From 7cc89f462a7eb82e219de6eb96241bc0ebd2c623 Mon Sep 17 00:00:00 2001 From: geay Date: Tue, 3 Jun 2014 18:26:47 +0200 Subject: [PATCH] stash 4 --- src/MEDCoupling/MEDCouplingAMRAttribute.cxx | 21 +- .../MEDCouplingCartesianAMRMesh.cxx | 201 ++++++++++++++---- .../MEDCouplingCartesianAMRMesh.hxx | 6 + src/MEDCoupling_Swig/MEDCouplingCommon.i | 8 + 4 files changed, 189 insertions(+), 47 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingAMRAttribute.cxx b/src/MEDCoupling/MEDCouplingAMRAttribute.cxx index 1a91be1b9..a3c652165 100644 --- a/src/MEDCoupling/MEDCouplingAMRAttribute.cxx +++ b/src/MEDCoupling/MEDCouplingAMRAttribute.cxx @@ -19,7 +19,9 @@ // Author : Anthony Geay #include "MEDCouplingAMRAttribute.hxx" +#include "MEDCouplingFieldDouble.hxx" #include "MEDCouplingMemArray.hxx" +#include "MEDCouplingIMesh.hxx" #include @@ -527,22 +529,31 @@ MEDCouplingFieldDouble *MEDCouplingAMRAttribute::buildCellFieldOnRecurseWithoutO /*! * This method builds a newly created field on cell just lying on mesh \a mesh without its eventual refinement. + * The output field also displays ghost cells. * * \return MEDCouplingFieldDouble * - a field on cells that the caller has to deal with (deallocate it). * */ MEDCouplingFieldDouble *MEDCouplingAMRAttribute::buildCellFieldOnWithGhost(int ghostLev, MEDCouplingCartesianAMRMeshGen *mesh, const std::string& fieldName) const { + const DataArrayDouble *arr(0); for(std::vector< MEDCouplingAutoRefCountObjectPtr >::const_iterator it=_levs.begin();it!=_levs.end();it++) { int tmp(-1); if((*it)->presenceOf(mesh,tmp)) { const DataArrayDoubleCollection& ddc((*it)->retrieveFieldsAt(tmp)); - const DataArrayDouble *arr(ddc.getFieldWithName(fieldName)); + arr=ddc.getFieldWithName(fieldName); } } - throw INTERP_KERNEL::Exception("MEDCouplingAMRAttribute::buildCellFieldOnWithGhost : the mesh specified is not in the progeny of this !"); + if(!arr) + throw INTERP_KERNEL::Exception("MEDCouplingAMRAttribute::buildCellFieldOnWithGhost : the mesh specified is not in the progeny of this !"); + MEDCouplingAutoRefCountObjectPtr im(mesh->getImageMesh()->buildWithGhost(ghostLev)); + MEDCouplingAutoRefCountObjectPtr ret(MEDCouplingFieldDouble::New(ON_CELLS)); + ret->setMesh(im); + ret->setArray(const_cast(arr)); + ret->setName(arr->getName()); + return ret.retn(); } /*! @@ -575,7 +586,7 @@ void MEDCouplingAMRAttribute::synchronizeCoarseToFine(int ghostLev) // for(std::size_t i=1;isynchronizeFineEachOther(ghostLev); diff --git a/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx b/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx index 719dd0646..9a71e5a92 100644 --- a/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx +++ b/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx @@ -102,7 +102,8 @@ int MEDCouplingCartesianAMRPatch::getNumberOfOverlapedCellsForFather() const /*! * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter - * the must be >= 0. + * the must be >= 0. \b WARNING this method only works if \a this and \a other share the same father (no check of this will be done !). + * Call isInMyNeighborhoodExt to deal with 2 patches not sharing directly the same father. * * \param [in] other - The other patch * \param [in] ghostLev - The size of the neighborhood zone. @@ -110,6 +111,8 @@ int MEDCouplingCartesianAMRPatch::getNumberOfOverlapedCellsForFather() const * \throw if \a this or \a other are invalid (end before start). * \throw if \a ghostLev is \b not >= 0 . * \throw if \a this and \a other have not the same space dimension. + * + * \sa isInMyNeighborhoodExt */ bool MEDCouplingCartesianAMRPatch::isInMyNeighborhood(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const { @@ -119,13 +122,54 @@ bool MEDCouplingCartesianAMRPatch::isInMyNeighborhood(const MEDCouplingCartesian throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the input patch is NULL !"); const std::vector< std::pair >& thisp(getBLTRRange()); const std::vector< std::pair >& otherp(other->getBLTRRange()); - std::size_t thispsize(thisp.size()); - if(thispsize!=otherp.size()) + return IsInMyNeighborhood(ghostLev,thisp,otherp); +} + +/*! + * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter + * the must be >= 0. This method works even if \a this and \a other does not share the same father. + * + * \param [in] other - The other patch + * \param [in] ghostLev - The size of the neighborhood zone. + * + * \throw if \a this or \a other are invalid (end before start). + * \throw if \a ghostLev is \b not >= 0 . + * \throw if \a this and \a other have not the same space dimension. + * \throw if there is not common ancestor of \a this and \a other. + * + * \sa isInMyNeighborhood + */ +bool MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const +{ + if(ghostLev<0) + throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt : the size of the neighborhood must be >= 0 !"); + if(!other) + throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt : the input patch is NULL !"); + int lev; + const MEDCouplingCartesianAMRMeshGen *com(FindCommonAncestor(this,other,lev));//check that factors are OK + if(lev==0) + return isInMyNeighborhood(other,ghostLev); + std::vector offset(ComputeOffsetFromTwoToOne(com,lev,this,other)); + const std::vector< std::pair >& thisp(getBLTRRange()); + std::vector< std::pair > otherp(other->getBLTRRange()); + std::size_t sz(offset.size()); + for(std::size_t i=0;i >& p1, const std::vector< std::pair >& p2) +{ + std::size_t thispsize(p1.size()); + if(thispsize!=p2.size()) throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the dimensions must be the same !"); for(std::size_t i=0;i& thispp(thisp[i]); - const std::pair& otherpp(otherp[i]); + const std::pair& thispp(p1[i]); + const std::pair& otherpp(p2[i]); if(thispp.second_mesh),*f2(p2->_mesh); + lev=0; + while(f1!=f2 || f1==0 || f2==0) + { + f1=f1->getFather(); f2=f2->getFather(); + if(f1->getFactors()!=f2->getFactors()) + throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindCommonAncestor : factors differ !"); + lev++; + } + if(f1!=f2) + throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindCommonAncestor : no common ancestor between p1 and p2 !"); + return f1; +} + +std::vector MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne(const MEDCouplingCartesianAMRMeshGen *comAncestor, int lev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2) +{ + if(lev<=0) + throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne : this method is useful only for lev > 0 !"); + int dim(p1->getMesh()->getSpaceDimension()); + if(p2->getMesh()->getSpaceDimension()!=dim) + throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne : dimension must be the same !"); + std::vector< int > ret(dim,0); + + for(int i=0;i_mesh),*f2(p2->_mesh); + const MEDCouplingCartesianAMRPatch *p1h(0),*p2h(0); + for(int j=0;jgetFather()),*f2tmp(f2->getFather()); + int pid1(f1tmp->getPatchIdFromChildMesh(f1)),pid2(f2tmp->getPatchIdFromChildMesh(f2)); + p1h=f1tmp->getPatch(pid1); p2h=f2tmp->getPatch(pid2); + f1=f1tmp; f2=f2tmp; + } + std::vector< std::pair > p2c(p2h->getBLTRRange()); + for(int k=0;kgetBLTRRange()[k].first; + ret[k]*=f1->getFactors()[k]; + } + } + return ret; +} + MEDCouplingCartesianAMRPatchGF::MEDCouplingCartesianAMRPatchGF(MEDCouplingCartesianAMRMesh *mesh):MEDCouplingCartesianAMRPatchGen(mesh) { } @@ -717,21 +812,8 @@ const MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRMeshGen::getPatch(int */ bool MEDCouplingCartesianAMRMeshGen::isPatchInNeighborhoodOf(int patchId1, int patchId2, int ghostLev) const { - if(ghostLev<0) - throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::isPatchInNeighborhoodOf : the ghost size must be >=0 !"); const MEDCouplingCartesianAMRPatch *p1(getPatch(patchId1)),*p2(getPatch(patchId2)); - if(_factors.empty()) - throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::isPatchInNeighborhoodOf : no factors defined !"); - int ghostLevInPatchRef; - if(ghostLev==0) - ghostLevInPatchRef=0; - else - { - ghostLevInPatchRef=(ghostLev-1)/_factors[0]+1; - for(std::size_t i=0;i<_factors.size();i++) - ghostLevInPatchRef=std::max(ghostLevInPatchRef,(ghostLev-1)/_factors[i]+1); - } - return p1->isInMyNeighborhood(p2,ghostLevInPatchRef); + return p1->isInMyNeighborhood(p2,GetGhostLevelInFineRef(ghostLev,_factors)); } /*! @@ -842,6 +924,8 @@ void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchGhostAdv(int patchId, c /*! * This method updates the patch with id \a patchId considering the only the all the patches in \a this to fill ghost zone. * So \b warning, the DataArrayDouble instance \a arrsOnPatches[patchId] is non const. + * + * \sa getPatchIdsInTheNeighborhoodOf */ void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyGhostAdv(int patchId, int ghostLev, const std::vector& arrsOnPatches) const { @@ -858,31 +942,27 @@ void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyGhostAdv(int patchI std::transform(dimsCoarse.begin(),dimsCoarse.end(),_factors.begin(),dimsCoarse.begin(),std::multiplies());//[12,8] std::transform(dimsCoarse.begin(),dimsCoarse.end(),dimsCoarse.begin(),std::bind2nd(std::plus(),2*ghostLev));//[14,10] std::vector< std::pair > rangeCoarse(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(dimsCoarse));//[(0,14),(0,10)] - std::vector fakeFactors(dim,1); + std::vector fakeFactors(dim,1),ids(getPatchIdsInTheNeighborhoodOf(patchId,ghostLev)); // - for(int i=0;i::const_iterator it=ids.begin();it!=ids.end();it++) { - if(i!=patchId) - if(isPatchInNeighborhoodOf(i,patchId,ghostLev)) - { - const MEDCouplingCartesianAMRPatch *otherP(getPatch(i)); - const std::vector< std::pair >& otherBLTR(otherP->getBLTRRange());//[(4,5),(3,4)] - std::vector< std::pair > tmp0,tmp1,tmp2; - MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(refBLTR,otherBLTR,tmp0,false);//tmp0=[(3,4),(1,2)] - ApplyFactorsOnCompactFrmt(tmp0,_factors);//tmp0=[(12,16),(4,8)] - ApplyGhostOnCompactFrmt(tmp0,ghostLev);//tmp0=[(13,17),(5,9)] - std::vector< std::pair > interstRange(MEDCouplingStructuredMesh::IntersectRanges(tmp0,rangeCoarse));//interstRange=[(13,14),(5,9)] - MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(otherBLTR,refBLTR,tmp1,false);//tmp1=[(-3,0),(-1,1)] - ApplyFactorsOnCompactFrmt(tmp1,_factors);//tmp1=[(-12,-4),(-4,0)] - MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(tmp1,interstRange,tmp2,false);//tmp2=[(1,2),(1,5)] - // - std::vector< std::pair > dimsFine(otherBLTR); - ApplyFactorsOnCompactFrmt(dimsFine,_factors); - ApplyAllGhostOnCompactFrmt(dimsFine,ghostLev); - // - MEDCouplingAutoRefCountObjectPtr ghostVals(MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(dimsFine),arrsOnPatches[i],tmp2)); - MEDCouplingIMesh::CondenseFineToCoarse(dimsCoarse,ghostVals,interstRange,fakeFactors,theFieldToFill); - } + const MEDCouplingCartesianAMRPatch *otherP(getPatch(*it)); + const std::vector< std::pair >& otherBLTR(otherP->getBLTRRange());//[(4,5),(3,4)] + std::vector< std::pair > tmp0,tmp1,tmp2; + MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(refBLTR,otherBLTR,tmp0,false);//tmp0=[(3,4),(1,2)] + ApplyFactorsOnCompactFrmt(tmp0,_factors);//tmp0=[(12,16),(4,8)] + ApplyGhostOnCompactFrmt(tmp0,ghostLev);//tmp0=[(13,17),(5,9)] + std::vector< std::pair > interstRange(MEDCouplingStructuredMesh::IntersectRanges(tmp0,rangeCoarse));//interstRange=[(13,14),(5,9)] + MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(otherBLTR,refBLTR,tmp1,false);//tmp1=[(-3,0),(-1,1)] + ApplyFactorsOnCompactFrmt(tmp1,_factors);//tmp1=[(-12,-4),(-4,0)] + MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(tmp1,interstRange,tmp2,false);//tmp2=[(1,2),(1,5)] + // + std::vector< std::pair > dimsFine(otherBLTR); + ApplyFactorsOnCompactFrmt(dimsFine,_factors); + ApplyAllGhostOnCompactFrmt(dimsFine,ghostLev); + // + MEDCouplingAutoRefCountObjectPtr ghostVals(MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(dimsFine),arrsOnPatches[*it],tmp2)); + MEDCouplingIMesh::CondenseFineToCoarse(dimsCoarse,ghostVals,interstRange,fakeFactors,theFieldToFill); } } @@ -1059,6 +1139,25 @@ DataArrayDouble *MEDCouplingCartesianAMRMeshGen::extractGhostFrom(int ghostSz, c return ret.retn(); } +/*! + * This method returns all the patches in \a this not equal to \a patchId that are in neighborhood of patch with id \a patchId. + * + * \sa fillCellFieldOnPatchOnlyGhostAdv + */ +std::vector MEDCouplingCartesianAMRMeshGen::getPatchIdsInTheNeighborhoodOf(int patchId, int ghostLev) const +{ + std::vector ret; + int nbp(getNumberOfPatches()); + // + for(int i=0;i& factors) +{ + if(ghostLev<0) + throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::GetGhostLevelInFineRef : the ghost size must be >=0 !"); + if(factors.empty()) + throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::GetGhostLevelInFineRef : no factors defined !"); + int ghostLevInPatchRef; + if(ghostLev==0) + ghostLevInPatchRef=0; + else + { + ghostLevInPatchRef=(ghostLev-1)/factors[0]+1; + for(std::size_t i=0;i >& getBLTRRange() const { return _bl_tr; } + MEDCOUPLING_EXPORT static bool IsInMyNeighborhood(int ghostLev, const std::vector< std::pair >& p1, const std::vector< std::pair >& p2); private: std::size_t getHeapMemorySizeWithoutChildren() const; + static const MEDCouplingCartesianAMRMeshGen *FindCommonAncestor(const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, int& lev); + static std::vector ComputeOffsetFromTwoToOne(const MEDCouplingCartesianAMRMeshGen *comAncestor, int lev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2); private: //! bottom left/top right cell range relative to \a _father std::vector< std::pair > _bl_tr; @@ -169,6 +173,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT MEDCoupling1SGTUMesh *buildMeshOfDirectChildrenOnly() const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildCellFieldOnRecurseWithoutOverlapWithoutGhost(int ghostSz, const std::vector& recurseArrs) const; MEDCOUPLING_EXPORT DataArrayDouble *extractGhostFrom(int ghostSz, const DataArrayDouble *arr) const; + MEDCOUPLING_EXPORT std::vector getPatchIdsInTheNeighborhoodOf(int patchId, int ghostLev) const; protected: MEDCouplingCartesianAMRMeshGen(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop, const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop); @@ -179,6 +184,7 @@ namespace ParaMEDMEM static void ApplyFactorsOnCompactFrmt(std::vector< std::pair >& partBeforeFact, const std::vector& factors); static void ApplyGhostOnCompactFrmt(std::vector< std::pair >& partBeforeFact, int ghostSize); static void ApplyAllGhostOnCompactFrmt(std::vector< std::pair >& partBeforeFact, int ghostSize); + static int GetGhostLevelInFineRef(int ghostLev, const std::vector& factors); std::vector extractSubTreeFromGlobalFlatten(const MEDCouplingCartesianAMRMeshGen *head, const std::vector& all) const; protected: MEDCOUPLING_EXPORT std::size_t getHeapMemorySizeWithoutChildren() const; diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 3a9582be5..b2e9d6e6c 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -4885,6 +4885,7 @@ namespace ParaMEDMEM int getPatchIdFromChildMesh(const MEDCouplingCartesianAMRMeshGen *mesh) const throw(INTERP_KERNEL::Exception); MEDCouplingUMesh *buildUnstructured() const throw(INTERP_KERNEL::Exception); DataArrayDouble *extractGhostFrom(int ghostSz, const DataArrayDouble *arr) const throw(INTERP_KERNEL::Exception); + std::vector getPatchIdsInTheNeighborhoodOf(int patchId, int ghostLev) const throw(INTERP_KERNEL::Exception); MEDCoupling1SGTUMesh *buildMeshFromPatchEnvelop() const throw(INTERP_KERNEL::Exception); MEDCoupling1SGTUMesh *buildMeshOfDirectChildrenOnly() const throw(INTERP_KERNEL::Exception); void removeAllPatches() throw(INTERP_KERNEL::Exception); @@ -4932,6 +4933,13 @@ namespace ParaMEDMEM return ret; } + MEDCouplingFieldDouble *buildCellFieldOnRecurseWithoutOverlapWithoutGhost(int ghostSz, PyObject *recurseArrs) const + { + std::vector inp; + convertFromPyObjVectorOfObj(recurseArrs,SWIGTYPE_p_ParaMEDMEM__DataArrayDouble,"DataArrayDouble",inp); + return self->buildCellFieldOnRecurseWithoutOverlapWithoutGhost(ghostSz,inp); + } + MEDCouplingCartesianAMRMeshGen *getFather() const throw(INTERP_KERNEL::Exception) { MEDCouplingCartesianAMRMeshGen *ret(const_cast(self->getFather())); -- 2.39.2