From eb54dd48c240bd9dc972d01dfa8b7a71c4b439cb Mon Sep 17 00:00:00 2001 From: geay Date: Wed, 18 Jun 2014 17:59:59 +0200 Subject: [PATCH] new algorithm for patch creation with additional options. --- src/INTERP_KERNEL/BoxSplittingOptions.cxx | 22 ++- src/INTERP_KERNEL/BoxSplittingOptions.hxx | 36 ++-- .../MEDCouplingCartesianAMRMesh.cxx | 175 ++++++++++-------- src/MEDCoupling/MEDCouplingStructuredMesh.cxx | 85 ++++++++- src/MEDCoupling/MEDCouplingStructuredMesh.hxx | 3 +- src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 15 +- src/MEDCoupling_Swig/MEDCouplingCommon.i | 19 +- 7 files changed, 233 insertions(+), 122 deletions(-) diff --git a/src/INTERP_KERNEL/BoxSplittingOptions.cxx b/src/INTERP_KERNEL/BoxSplittingOptions.cxx index d932a482d..91f2721f4 100644 --- a/src/INTERP_KERNEL/BoxSplittingOptions.cxx +++ b/src/INTERP_KERNEL/BoxSplittingOptions.cxx @@ -22,24 +22,26 @@ #include -const double INTERP_KERNEL::BoxSplittingOptions::DFT_EFFECIENCY=0.5; +const double INTERP_KERNEL::BoxSplittingOptions::DFT_EFFICIENCY_GOAL=0.5; -const double INTERP_KERNEL::BoxSplittingOptions::DFT_EFFECIENCY_SND=0.5; +const double INTERP_KERNEL::BoxSplittingOptions::DFT_EFFICIENCY_THRESHOLD=0.7; void INTERP_KERNEL::BoxSplittingOptions::init() { - _effeciency=DFT_EFFECIENCY; - _effeciency_snd=DFT_EFFECIENCY_SND; - _min_cell_direction=DFT_MIN_CELL_DIRECTION; - _max_cells=DFT_MAX_CELLS; + _efficiency_goal=DFT_EFFICIENCY_GOAL; + _efficiency_threshold=DFT_EFFICIENCY_THRESHOLD; + _min_patch_length=DFT_MIN_PATCH_LENGTH; + _max_patch_length=DFT_MAX_PATCH_LENGTH; + _max_nb_cells_in_patch=DFT_MAX_PATCH_MEASURE; } std::string INTERP_KERNEL::BoxSplittingOptions::printOptions() const { std::ostringstream oss; - oss << "Efficiency threshold: " << 100*_effeciency << "%"<< std::endl; - oss << "Efficiency Snd threshold: " << 100*_effeciency_snd << "%"<< std::endl; - oss << "Min. box side length: " << _min_cell_direction << std::endl; - oss << "Max. box cells : " << _max_cells << std::endl; + oss << "Efficiency goal: " << 100*_efficiency_goal << "%" << std::endl; + oss << "Efficiency threshold: " << 100*_efficiency_threshold << "%" << std::endl; + oss << "Min. patch side length: " << _min_patch_length << std::endl; + oss << "Max. patch side length: " << _max_patch_length << std::endl; + oss << "Max. patch measure: " << _max_nb_cells_in_patch << std::endl; return oss.str(); } diff --git a/src/INTERP_KERNEL/BoxSplittingOptions.hxx b/src/INTERP_KERNEL/BoxSplittingOptions.hxx index f8e47caed..5637c73f8 100644 --- a/src/INTERP_KERNEL/BoxSplittingOptions.hxx +++ b/src/INTERP_KERNEL/BoxSplittingOptions.hxx @@ -35,28 +35,32 @@ namespace INTERP_KERNEL class INTERPKERNEL_EXPORT BoxSplittingOptions { private: - double _effeciency; - double _effeciency_snd; - int _min_cell_direction; - int _max_cells; + double _efficiency_goal; + double _efficiency_threshold; + int _min_patch_length; + int _max_patch_length; + int _max_nb_cells_in_patch; public: BoxSplittingOptions() { init(); } void init(); - double getEffeciency() const { return _effeciency; } - void setEffeciency(double effeciency) { _effeciency=effeciency; } - double getEffeciencySnd() const { return _effeciency_snd; } - void setEffeciencySnd(double effeciencySnd) { _effeciency_snd=effeciencySnd; } - int getMinCellDirection() const { return _min_cell_direction; } - void setMinCellDirection(int minCellDirection) { _min_cell_direction=minCellDirection; } - int getMaxCells() const { return _max_cells; } - void setMaxCells(int maxCells) { _max_cells=maxCells; } + double getEfficiencyGoal() const { return _efficiency_goal; } + void setEfficiencyGoal(double efficiency) { _efficiency_goal=efficiency; } + double getEfficiencyThreshold() const { return _efficiency_threshold; } + void setEfficiencyThreshold(double efficiencyThreshold) { _efficiency_threshold=efficiencyThreshold; } + int getMinimumPatchLength() const { return _min_patch_length; } + void setMinimumPatchLength(int minPatchLength) { _min_patch_length=minPatchLength; } + int getMaximumPatchLength() const { return _max_patch_length; } + void setMaximumPatchLength(int maxNbCellPatch) { _max_patch_length=maxNbCellPatch; } + int getMaximumNbOfCellsInPatch() const { return _max_nb_cells_in_patch; } + void setMaximumNbOfCellsInPatch(int maxNbCellsInPatch) { _max_nb_cells_in_patch=maxNbCellsInPatch; } void copyOptions(const BoxSplittingOptions & other) { *this=other; } std::string printOptions() const; private: - static const int DFT_MIN_CELL_DIRECTION=3; - static const int DFT_MAX_CELLS=1000; - static const double DFT_EFFECIENCY; - static const double DFT_EFFECIENCY_SND; + static const int DFT_MIN_PATCH_LENGTH=1; + static const int DFT_MAX_PATCH_LENGTH=2147483647; + static const int DFT_MAX_PATCH_MEASURE=2147483647; + static const double DFT_EFFICIENCY_GOAL; + static const double DFT_EFFICIENCY_THRESHOLD; public: static const char PRECISION_STR[]; }; diff --git a/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx b/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx index c8c117530..5fe253057 100644 --- a/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx +++ b/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx @@ -737,7 +737,7 @@ public: std::vector computeCGS() const { return MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(_part); } std::vector< std::vector > computeSignature() const { return MEDCouplingStructuredMesh::ComputeSignaturePerAxisOf(computeCGS(),getConstCriterion()); } double getEfficiencyPerAxis(int axisId) const { return (double)_nb_of_true/((double)(_part[axisId].second-_part[axisId].first)); } - void zipToFitOnCriterion(); + void zipToFitOnCriterion(int minPatchLgth); void updateNumberOfTrue() const; MEDCouplingAutoRefCountObjectPtr extractPart(const std::vector< std::pair >&partInGlobal) const; MEDCouplingAutoRefCountObjectPtr deepCpy() const; @@ -750,12 +750,12 @@ private: std::vector< std::pair > _part; }; -void InternalPatch::zipToFitOnCriterion() +void InternalPatch::zipToFitOnCriterion(int minPatchLgth) { std::vector cgs(computeCGS()); std::vector newCrit; std::vector< std::pair > newPart,newPart2; - int newNbOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,_crit,newCrit,newPart)); + int newNbOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(minPatchLgth,cgs,_crit,newCrit,newPart)); MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(_part,newPart,newPart2); if(newNbOfTrue!=_nb_of_true) throw INTERP_KERNEL::Exception("InternalPatch::zipToFitOnCrit : internal error !"); @@ -786,46 +786,47 @@ MEDCouplingAutoRefCountObjectPtr InternalPatch::deepCpy() const return ret; } -void DissectBigPatch(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace) +void DissectBigPatch(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int largestLength, int& cutPlace) { - cutFound=false; cutPlace=-1; - std::vector ratio(rangeOfAxisId-1); - for(int id=0;id ratio(largestLength-minimumPatchLength,std::numeric_limits::max()); + const int dim(patchToBeSplit->getDimension()); + int index_min = -1; + double minSemiEfficiencyRatio(std::numeric_limits::max()); + double efficiencyPerAxis[2]; + + for(int i=minimumPatchLength-1;i > rectH(patchToBeSplit->getConstPart()); if(h==0) - rectH[axisId].second=patchToBeSplit->getConstPart()[axisId].first+id; + rectH[axisId].second=patchToBeSplit->getConstPart()[axisId].first+i; else - rectH[axisId].first=patchToBeSplit->getConstPart()[axisId].first+id; + rectH[axisId].first=patchToBeSplit->getConstPart()[axisId].first+i; MEDCouplingAutoRefCountObjectPtr p(patchToBeSplit->deepCpy()); - p->zipToFitOnCriterion(); - //anouar rectH ? - efficiency[h]=p->getEfficiencyPerAxis(axisId); + p->zipToFitOnCriterion(bso.getMinimumPatchLength()); + efficiencyPerAxis[h]=p->getEfficiencyPerAxis(axisId); } - ratio[id]=std::max(efficiency[0],efficiency[1])/std::min(efficiency[0],efficiency[1]); - } - int minCellDirection(bso.getMinCellDirection()),indexMin(-1); - int dimRatioInner(rangeOfAxisId-1-2*(minCellDirection-1)); - std::vector ratio_inner(dimRatioInner); - double minRatio(1.e10); - for(int i=0; igetConstPart()[axisId].first-1; + + if(index_min==-1) + throw INTERP_KERNEL::Exception("DissectBigPatch : just call to Arthur !"); + + cutPlace=index_min+patchToBeSplit->getConstPart()[axisId].first; } -void FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int& axisId, bool& cutFound, int& cutPlace) +bool FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int& cutPlace) { - cutPlace=-1; cutFound=false; - int minCellDirection(bso.getMinCellDirection()); + cutPlace=-1; + int minimumPatchLength(bso.getMinimumPatchLength()); const int dim(patchToBeSplit->getDimension()); std::vector< std::vector > signatures(patchToBeSplit->computeSignature()); for(int id=0;id hole; std::vector distance; int len((int)signature.size()); - for(int i=0;i= 2*minCellDirection && i >= minCellDirection-1 && i <= len-minCellDirection-1) - hole.push_back(i); + hole.push_back(i); if(!hole.empty()) { - double center(((double)len/2.)); + int closestHoleToMiddle(hole[0]); + int oldDistanceToMiddle(std::abs(hole[0]-len/2)); + int newDistanceToMiddle(oldDistanceToMiddle); for(std::size_t i=0;igetConstPart()[axisId].first+1; - return ; + { + newDistanceToMiddle=std::abs(hole[i]-len/2); + if(newDistanceToMiddle < oldDistanceToMiddle) + { + oldDistanceToMiddle = newDistanceToMiddle; + closestHoleToMiddle = hole[i]; + } + } + cutPlace=closestHoleToMiddle+patchToBeSplit->getConstPart()[axisId].first; + return true; } } + return false; } -void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, bool& cutFound, int& cutPlace, int& axisId) +bool FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int& cutPlace, int& axisId) { - cutFound=false; cutPlace=-1;// do not set axisId before to be sure that cutFound was set to true - + bool cutFound(false); cutPlace=-1;// do not set axisId before to be sure that cutFound was set to true const std::vector< std::pair >& part(patchToBeSplit->getConstPart()); - int sign,minCellDirection(bso.getMinCellDirection()); + int sign,minimumPatchLength(bso.getMinimumPatchLength()); const int dim(patchToBeSplit->getDimension()); std::vector zeroCrossDims(dim,-1); @@ -868,18 +873,18 @@ void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const Interna { const std::vector& signature(signatures[id]); - std::vector derivate_second_order,gradient_absolute,signe_change,zero_cross,edge,max_cross_list ; + std::vector derivate_second_order,gradient_absolute,zero_cross,edge,max_cross_list ; std::vector distance ; - for (unsigned int i=1;i= (unsigned int)minCellDirection-2 && i <= signature.size()-minCellDirection-2 ) + if ( i >= (std::size_t)minimumPatchLength-2 && i <= signature.size()-minimumPatchLength-2 ) { zero_cross.push_back(i) ; edge.push_back(gradient_absolute[i]) ; } - signe_change.push_back(sign) ; } if ( zero_cross.size() > 0 ) { @@ -906,7 +910,7 @@ void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const Interna for (unsigned int i=0;i=0 ) @@ -917,7 +921,6 @@ void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const Interna } derivate_second_order.clear() ; gradient_absolute.clear() ; - signe_change.clear() ; zero_cross.clear() ; edge.clear() ; max_cross_list.clear() ; @@ -943,26 +946,30 @@ void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const Interna cutPlace=zeroCrossDims[max_cross_dims]; axisId=max_cross_dims ; } + return cutFound; } -void TryAction4(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace) +bool TryAction4(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, int& cutPlace) { - cutFound=false; - if(patchToBeSplit->getEfficiency()<=bso.getEffeciencySnd()) + if(patchToBeSplit->getEfficiency()<=bso.getEfficiencyGoal()) { - if(rangeOfAxisId>=2*bso.getMinCellDirection()) + if(rangeOfAxisId>=2*bso.getMinimumPatchLength()) { - cutFound=true; cutPlace=rangeOfAxisId/2+patchToBeSplit->getConstPart()[axisId].first-1; } + else + return false; } else { - if(patchToBeSplit->getNumberOfCells()>bso.getMaxCells()) + if(patchToBeSplit->getNumberOfCells()>bso.getMaximumNbOfCellsInPatch() || rangeOfAxisId>bso.getMaximumPatchLength()) { - DissectBigPatch(bso,patchToBeSplit,axisId,rangeOfAxisId,cutFound,cutPlace); + DissectBigPatch(bso,patchToBeSplit,axisId,rangeOfAxisId,cutPlace); } + else + return false; } + return true; } MEDCouplingAutoRefCountObjectPtr DealWithNoCut(const InternalPatch *patch) @@ -972,7 +979,7 @@ MEDCouplingAutoRefCountObjectPtr DealWithNoCut(const InternalPatc return ret; } -void DealWithCut(const InternalPatch *patchToBeSplit, int axisId, int cutPlace, std::vector >& listOfPatches) +void DealWithCut(double minPatchLgth, const InternalPatch *patchToBeSplit, int axisId, int cutPlace, std::vector >& listOfPatches) { MEDCouplingAutoRefCountObjectPtr leftPart,rightPart; std::vector< std::pair > rect(patchToBeSplit->getConstPart()); @@ -981,7 +988,7 @@ void DealWithCut(const InternalPatch *patchToBeSplit, int axisId, int cutPlace, rightRect[axisId].first=cutPlace+1; leftPart=patchToBeSplit->extractPart(leftRect); rightPart=patchToBeSplit->extractPart(rightRect); - leftPart->zipToFitOnCriterion(); rightPart->zipToFitOnCriterion(); + leftPart->zipToFitOnCriterion(minPatchLgth); rightPart->zipToFitOnCriterion(minPatchLgth); listOfPatches.push_back(leftPart); listOfPatches.push_back(rightPart); } @@ -1013,8 +1020,21 @@ int MEDCouplingCartesianAMRMeshGen::getNumberOfPatches() const } /*! - * This method creates patches in \a this (by destroying the patches if any). This method uses \a criterion array as a field on cells on this level. + * This method is a generic algorithm to create patches in \a this (by destroying the patches if any). + * This method uses \a criterion array as a field on cells on this level. * This method only create patches at level 0 relative to \a this. + * + * This generic algorithm can be degenerated into three child ones, depending on the arguments given; in particular depending + * on whether they are equal to 0 or not. + * 1/ If minimumPatchLength = maximumPatchLength = maximumPatchVolume = 0, then we have the Berger-Rigoutsos algorithm. + * This algorithm was developed in 1991 and seems appropriate for sequential programming. + * 2/ If maximumPatchLength = 0, then we have the Livne algorithm. + * This algorithm was developed in 2004 and is an improvement of the Berger-Rigoutsos algorithm. + * 3/ If maximumPatchVolume = 0, the we have the lmin-lmax algorithm. + * This algorithm was developed by Arthur TALPAERT in 2014 and is an improvement of the Livne algorithm. It is especially + * appropriate for parallel computing, where one patch would be given to one CPU. See Arthur TALPAERT's 2014 CANUM poster + * for more information. + * */ void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const std::vector& criterion, const std::vector& factors) { @@ -1026,7 +1046,7 @@ void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KER std::vector< MEDCouplingAutoRefCountObjectPtr > listOfPatches,listOfPatchesOK; // MEDCouplingAutoRefCountObjectPtr p(new InternalPatch); - p->setNumberOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,criterion,p->getCriterion(),p->getPart())); + p->setNumberOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(bso.getMinimumPatchLength(),cgs,criterion,p->getCriterion(),p->getPart())); if(p->presenceOfTrue()) listOfPatches.push_back(p); while(!listOfPatches.empty()) @@ -1035,21 +1055,22 @@ void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KER for(std::vector< MEDCouplingAutoRefCountObjectPtr >::iterator it=listOfPatches.begin();it!=listOfPatches.end();it++) { // - int axisId,rangeOfAxisId,cutPlace; - bool cutFound; - MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt((*it)->getConstPart(),axisId,rangeOfAxisId); - if((*it)->getEfficiency()>=bso.getEffeciency() && (*it)->getNumberOfCells()getConstPart(),axisId,largestLength); + if((*it)->getEfficiency()>=bso.getEfficiencyThreshold() && ((*it)->getNumberOfCells()>bso.getMaximumNbOfCellsInPatch() || largestLength>bso.getMaximumPatchLength())) + { + DissectBigPatch(bso,*it,axisId,largestLength,cutPlace); + DealWithCut(bso.getMinimumPatchLength(),*it,axisId,cutPlace,listOfPatchesTmp); + continue; + }//action 1 + if(FindHole(bso,*it,axisId,cutPlace))//axisId overwritten here if FindHole equal to true ! + { DealWithCut(bso.getMinimumPatchLength(),*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 2 + if(FindInflection(bso,*it,cutPlace,axisId))//axisId overwritten here if cutFound equal to true ! + { DealWithCut(bso.getMinimumPatchLength(),*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 3 + if(TryAction4(bso,*it,axisId,largestLength,cutPlace)) + { DealWithCut(bso.getMinimumPatchLength(),*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 4 + else + listOfPatchesOK.push_back(DealWithNoCut(*it)); } listOfPatches=listOfPatchesTmp; } diff --git a/src/MEDCoupling/MEDCouplingStructuredMesh.cxx b/src/MEDCoupling/MEDCouplingStructuredMesh.cxx index 58ec89f01..1e2186855 100644 --- a/src/MEDCoupling/MEDCouplingStructuredMesh.cxx +++ b/src/MEDCoupling/MEDCouplingStructuredMesh.cxx @@ -803,14 +803,17 @@ void MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt(const * \a partCompactFormat that contains all the True in \a crit. The returned vector of boolean is the field reduced to that part. * So the number of True is equal in \a st and in returned vector of boolean. * + * \param [in] minPatchLgth - minimum length that the patch may have for all directions. * \param [in] st - The structure per axis of the structured mesh considered. * \param [in] crit - The field of boolean (for performance reasons) lying on the mesh defined by \a st. * \param [out] partCompactFormat - The minimal part of \a st containing all the true of \a crit. * \param [out] reducedCrit - The reduction of \a criterion on \a partCompactFormat. * \return - The number of True in \a st (that is equal to those in \a reducedCrit) */ -int MEDCouplingStructuredMesh::FindMinimalPartOf(const std::vector& st, const std::vector& crit, std::vector& reducedCrit, std::vector< std::pair >& partCompactFormat) +int MEDCouplingStructuredMesh::FindMinimalPartOf(int minPatchLgth, const std::vector& st, const std::vector& crit, std::vector& reducedCrit, std::vector< std::pair >& partCompactFormat) { + if(minPatchLgth<0) + throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindMinimalPartOf : the input minPatchLgth has to be >=0 !"); if((int)crit.size()!=DeduceNumberOfGivenStructure(st)) throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindMinimalPartOf : size of vector of boolean is invalid regarding the declared structure !"); int ret(-1); @@ -834,6 +837,29 @@ int MEDCouplingStructuredMesh::FindMinimalPartOf(const std::vector& st, con default: throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindMinimalPartOf : only dimension 1, 2 and 3 are supported actually !"); } + std::vector dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(partCompactFormat)); + int i(0); + for(std::vector< std::pair >::iterator it=partCompactFormat.begin();it!=partCompactFormat.end();it++,i++) + { + if(st[i]st[i]) + { + (*it).first-=(*it).second-st[i]; + (*it).second=st[i]; + } + } + } ExtractFieldOfBoolFrom(st,crit,partCompactFormat,reducedCrit); return ret; } @@ -1036,7 +1062,16 @@ int MEDCouplingStructuredMesh::FindMinimalPartOf1D(const std::vector& st, c } } if(ret==0) - return ret; + { + std::size_t sz(st.size()); + partCompactFormat.resize(sz); + for(std::size_t i=0;i& st, c } } if(ret==0) - return ret; + { + std::size_t sz(st.size()); + partCompactFormat.resize(sz); + for(std::size_t i=0;i& st, c } } if(ret==0) - return ret; + { + std::size_t sz(st.size()); + partCompactFormat.resize(sz); + for(std::size_t i=0;i MEDCouplingStructuredMesh::getCellGridStructure() const return ret; } +/*! + * This method returns the squareness of \a this (quadrature). \a this is expected to be with a mesh dimension equal to 2 or 3. + */ +double MEDCouplingStructuredMesh::computeSquareness() const +{ + std::vector cgs(getCellGridStructure()); + if(cgs.empty()) + throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::computeSquareness : empty mesh !"); + std::size_t dim(cgs.size()); + if(dim==1) + throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::computeSquareness : A segment cannot be square !"); + if(dim<4) + { + int minAx(cgs[0]),maxAx(cgs[0]); + for(std::size_t i=1;i getNodeGridStructure() const = 0; MEDCOUPLING_EXPORT std::vector getCellGridStructure() const; + MEDCOUPLING_EXPORT double computeSquareness() const; MEDCOUPLING_EXPORT virtual MEDCouplingStructuredMesh *buildStructuredSubPart(const std::vector< std::pair >& cellPart) const = 0; MEDCOUPLING_EXPORT static std::vector GetSplitVectFromStruct(const std::vector& strct); MEDCOUPLING_EXPORT static bool IsPartStructured(const int *startIds, const int *stopIds, const std::vector& st, std::vector< std::pair >& partCompactFormat); @@ -93,7 +94,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT static int DeduceNumberOfGivenRangeInCompactFrmt(const std::vector< std::pair >& partCompactFormat); MEDCOUPLING_EXPORT static int DeduceNumberOfGivenStructure(const std::vector& st); MEDCOUPLING_EXPORT static void FindTheWidestAxisOfGivenRangeInCompactFrmt(const std::vector< std::pair >& partCompactFormat, int& axisId, int& sizeOfRange); - MEDCOUPLING_EXPORT static int FindMinimalPartOf(const std::vector& st, const std::vector& crit, std::vector& reducedCrit, std::vector< std::pair >& partCompactFormat); + MEDCOUPLING_EXPORT static int FindMinimalPartOf(int minPatchLgth, const std::vector& st, const std::vector& crit, std::vector& reducedCrit, std::vector< std::pair >& partCompactFormat); MEDCOUPLING_EXPORT static std::vector< std::vector > ComputeSignaturePerAxisOf(const std::vector& st, const std::vector& crit); private: static int GetNumberOfCellsOfSubLevelMesh(const std::vector& cgs, int mdim); diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index 1bf0b4efc..31c3d2ac6 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -15124,16 +15124,16 @@ class MEDCouplingBasicsTest(unittest.TestCase): # f.write("test.vti") amr=MEDCouplingCartesianAMRMesh("mesh",2,[51,51],[0.,0.],[0.04,0.04]) arr2=DataArrayByte(im.getNumberOfCells()) ; arr2[:]=0 ; arr2[ids]=1 - bso=BoxSplittingOptions() ; bso.setEffeciency(0.8) ; bso.setEffeciencySnd(0.8) ; bso.setMaxCells(1000) ; bso.setMinCellDirection(3) + bso=BoxSplittingOptions() ; bso.setEfficiencyGoal(0.5); bso.setEfficiencyThreshold(0.8) ; bso.setMaximumNbOfCellsInPatch(3000) ; bso.setMinimumPatchLength(6) ; bso.setMaximumPatchLength(11) amr.createPatchesFromCriterion(bso,arr2,[2,2]) - self.assertEqual(18,amr.getNumberOfPatches()) - exp0=[[(8,14),(19,38)],[(19,31),(8,17)],[(19,31),(33,42)],[(10,14),(12,16)],[(9,14),(16,19)],[(14,19),(9,19)],[(14,17),(19,22)],[(14,19),(31,41)],[(36,42),(19,38)],[(14,15),(22,28)],[(14,17),(28,31)],[(31,36),(9,19)],[(33,36),(19,22)],[(31,36),(31,41)],[(36,40),(12,16)],[(36,41),(16,19)],[(35,36),(22,28)],[(33,36),(28,31)]] + amr.buildMeshFromPatchEnvelop().writeVTK("toto.vtu") + m=amr.getImageMesh() ; m=m.buildUnstructured() ; m.changeSpaceDimension(3,1.); m.writeVTK("grid.vtu") + self.assertEqual(12,amr.getNumberOfPatches()) + exp0=[[(9,19),(9,19)],[(9,19),(31,41)],[(31,41),(9,19)],[(8,17),(19,25)],[(8,17),(25,31)],[(19,25),(8,17)],[(25,31),(8,17)],[(19,25),(33,42)],[(25,31),(33,42)],[(31,41),(31,41)],[(33,42),(19,25)],[(33,42),(25,31)]] for i,bltr in enumerate(exp0): self.assertEqual(amr[i].getBLTRRange(),bltr) pass - m=amr.buildMeshFromPatchEnvelop() - self.assertTrue(m.getNodalConnectivity().isEqual(DataArrayInt([1,0,2,3,5,4,6,7,9,8,10,11,13,12,14,15,17,16,18,19,21,20,22,23,25,24,26,27,29,28,30,31,33,32,34,35,37,36,38,39,41,40,42,43,45,44,46,47,49,48,50,51,53,52,54,55,57,56,58,59,61,60,62,63,65,64,66,67,69,68,70,71]))) - self.assertTrue(m.getCoords().isEqualWithoutConsideringStr(DataArrayDouble([0.32,0.76,0.56,0.76,0.32,1.52,0.56,1.52,0.76,0.32,1.24,0.32,0.76,0.68,1.24,0.68,0.76,1.32,1.24,1.32,0.76,1.68,1.24,1.68,0.4,0.48,0.56,0.48,0.4,0.64,0.56,0.64,0.36,0.64,0.56,0.64,0.36,0.76,0.56,0.76,0.56,0.36,0.76,0.36,0.56,0.76,0.76,0.76,0.56,0.76,0.68,0.76,0.56,0.88,0.68,0.88,0.56,1.24,0.76,1.24,0.56,1.64,0.76,1.64,1.44,0.76,1.68,0.76,1.44,1.52,1.68,1.52,0.56,0.88,0.6,0.88,0.56,1.12,0.6,1.12,0.56,1.12,0.68,1.12,0.56,1.24,0.68,1.24,1.24,0.36,1.44,0.36,1.24,0.76,1.44,0.76,1.32,0.76,1.44,0.76,1.32,0.88,1.44,0.88,1.24,1.24,1.44,1.24,1.24,1.64,1.44,1.64,1.44,0.48,1.6,0.48,1.44,0.64,1.6,0.64,1.44,0.64,1.64,0.64,1.44,0.76,1.64,0.76,1.4,0.88,1.44,0.88,1.4,1.12,1.44,1.12,1.32,1.12,1.44,1.12,1.32,1.24,1.44,1.24],72,2),1e-12)) + self.assertAlmostEqual(0.666666666667,amr[3].getMesh().getImageMesh().computeSquareness(),12) # self.assertEqual(MEDCouplingStructuredMesh.ChangeReferenceToGlobalOfCompactFrmt([(8,32),(4,17)],[(0,24),(2,12)]),[(8,32),(6,16)]) self.assertEqual(MEDCouplingStructuredMesh.ChangeReferenceFromGlobalOfCompactFrmt([(8,32),(4,17)],[(8,32),(6,16)]),[(0,24),(2,12)]) @@ -15254,6 +15254,9 @@ class MEDCouplingBasicsTest(unittest.TestCase): amr[0].addPatch([(10,12),(5,8)],[2,2]) amr[1].addPatch([(0,1),(0,5)],[2,2]) amr[2].addPatch([(3,4),(0,3)],[2,2]) + m=amr.buildMeshFromPatchEnvelop() + self.assertTrue(m.getNodalConnectivity().isEqual(DataArrayInt([1,0,2,3,5,4,6,7,9,8,10,11]))) + self.assertTrue(m.getCoords().isEqualWithoutConsideringStr(DataArrayDouble([1.,2.,4.,2.,1.,4.,4.,4.,4.,3.,5.,3.,4.,5.,5.,5.,0.,4.,1.,4.,0.,6.,1.,6.],12,2),1e-12)) self.assertEqual(3,amr.getMaxNumberOfLevelsRelativeToThis()) att=MEDCouplingAMRAttribute(amr,[("Field",["X"])],ghostSz) att.alloc() diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 189db27d9..df2e0fe53 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -427,14 +427,16 @@ namespace INTERP_KERNEL public: BoxSplittingOptions(); void init() throw(INTERP_KERNEL::Exception); - double getEffeciency() const throw(INTERP_KERNEL::Exception); - void setEffeciency(double effeciency) throw(INTERP_KERNEL::Exception); - double getEffeciencySnd() const throw(INTERP_KERNEL::Exception); - void setEffeciencySnd(double effeciencySnd) throw(INTERP_KERNEL::Exception); - int getMinCellDirection() const throw(INTERP_KERNEL::Exception); - void setMinCellDirection(int minCellDirection) throw(INTERP_KERNEL::Exception); - int getMaxCells() const throw(INTERP_KERNEL::Exception); - void setMaxCells(int maxCells) throw(INTERP_KERNEL::Exception); + double getEfficiencyGoal() const throw(INTERP_KERNEL::Exception); + void setEfficiencyGoal(double efficiency) throw(INTERP_KERNEL::Exception); + double getEfficiencyThreshold() const throw(INTERP_KERNEL::Exception); + void setEfficiencyThreshold(double efficiencyThreshold) throw(INTERP_KERNEL::Exception); + int getMinimumPatchLength() const throw(INTERP_KERNEL::Exception); + void setMinimumPatchLength(int minPatchLength) throw(INTERP_KERNEL::Exception); + int getMaximumPatchLength() const throw(INTERP_KERNEL::Exception); + void setMaximumPatchLength(int maxPatchLength) throw(INTERP_KERNEL::Exception); + int getMaximumNbOfCellsInPatch() const throw(INTERP_KERNEL::Exception); + void setMaximumNbOfCellsInPatch(int maxNbCellsInPatch) throw(INTERP_KERNEL::Exception); void copyOptions(const BoxSplittingOptions & other) throw(INTERP_KERNEL::Exception); std::string printOptions() const throw(INTERP_KERNEL::Exception); %extend @@ -2865,6 +2867,7 @@ namespace ParaMEDMEM int getNodeIdFromPos(int i, int j, int k) const throw(INTERP_KERNEL::Exception); int getNumberOfCellsOfSubLevelMesh() const throw(INTERP_KERNEL::Exception); int getSpaceDimensionOnNodeStruct() const throw(INTERP_KERNEL::Exception); + double computeSquareness() const throw(INTERP_KERNEL::Exception); virtual std::vector getNodeGridStructure() const throw(INTERP_KERNEL::Exception); std::vector getCellGridStructure() const throw(INTERP_KERNEL::Exception); MEDCoupling1SGTUMesh *build1SGTUnstructured() const throw(INTERP_KERNEL::Exception); -- 2.39.2