#include <sstream>
-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();
}
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[];
};
std::vector<int> computeCGS() const { return MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(_part); }
std::vector< std::vector<int> > 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<InternalPatch> extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const;
MEDCouplingAutoRefCountObjectPtr<InternalPatch> deepCpy() const;
std::vector< std::pair<int,int> > _part;
};
-void InternalPatch::zipToFitOnCriterion()
+void InternalPatch::zipToFitOnCriterion(int minPatchLgth)
{
std::vector<int> cgs(computeCGS());
std::vector<bool> newCrit;
std::vector< std::pair<int,int> > 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 !");
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<double> ratio(rangeOfAxisId-1);
- for(int id=0;id<rangeOfAxisId-1;id++)
+ int minimumPatchLength(bso.getMinimumPatchLength());
+ std::vector<double> ratio(largestLength-minimumPatchLength,std::numeric_limits<double>::max());
+ const int dim(patchToBeSplit->getDimension());
+ int index_min = -1;
+ double minSemiEfficiencyRatio(std::numeric_limits<double>::max());
+ double efficiencyPerAxis[2];
+
+ for(int i=minimumPatchLength-1;i<largestLength-minimumPatchLength;i++)
{
- double efficiency[2];
+ int numberOfFlags_h;
for(int h=0;h<2;h++)
{
std::vector< std::pair<int,int> > 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<InternalPatch> 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<double> ratio_inner(dimRatioInner);
- double minRatio(1.e10);
- for(int i=0; i<dimRatioInner; i++)
- {
- if(ratio[minCellDirection-1+i]<minRatio)
+ ratio[i]=std::max(efficiencyPerAxis[0], efficiencyPerAxis[1]) / std::min(efficiencyPerAxis[0], efficiencyPerAxis[1]);
+ if(ratio[i]<minSemiEfficiencyRatio)
{
- minRatio=ratio[minCellDirection-1+i];
- indexMin=i+minCellDirection;
+ minSemiEfficiencyRatio = ratio[i];
+ index_min = i;
}
}
- cutFound=true; cutPlace=indexMin+patchToBeSplit->getConstPart()[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<int> > signatures(patchToBeSplit->computeSignature());
for(int id=0;id<dim;id++)
std::vector<int> hole;
std::vector<double> distance;
int len((int)signature.size());
- for(int i=0;i<len;i++)
+ for(int i=minimumPatchLength-1;i<len-minimumPatchLength;i++)
if(signature[i]==0)
- if(len>= 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;i<hole.size();i++)
- distance.push_back(fabs(hole[i]+1.-center));
-
- std::size_t posDistanceMin(std::distance(distance.begin(),std::min_element(distance.begin(),distance.end())));
- cutFound=true;
- axisId=id;
- cutPlace=hole[posDistanceMin]+patchToBeSplit->getConstPart()[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<int,int> >& part(patchToBeSplit->getConstPart());
- int sign,minCellDirection(bso.getMinCellDirection());
+ int sign,minimumPatchLength(bso.getMinimumPatchLength());
const int dim(patchToBeSplit->getDimension());
std::vector<int> zeroCrossDims(dim,-1);
{
const std::vector<int>& signature(signatures[id]);
- std::vector<int> derivate_second_order,gradient_absolute,signe_change,zero_cross,edge,max_cross_list ;
+ std::vector<int> derivate_second_order,gradient_absolute,zero_cross,edge,max_cross_list ;
std::vector<double> distance ;
- for (unsigned int i=1;i<signature.size()-1;i++)
+ for(std::size_t i=1;i<signature.size()-1;i++)
derivate_second_order.push_back(signature[i-1]-2*signature[i]+signature[i+1]) ;
// Gradient absolute value
- for ( unsigned int i=1;i<derivate_second_order.size();i++)
+ for(std::size_t i=1;i<derivate_second_order.size();i++)
gradient_absolute.push_back(fabs(derivate_second_order[i]-derivate_second_order[i-1])) ;
if(derivate_second_order.empty())
continue;
- for (unsigned int i=0;i<derivate_second_order.size()-1;i++)
+ for(std::size_t i=1;i<derivate_second_order.size()-1;i++)
{
if (derivate_second_order[i]*derivate_second_order[i+1] < 0 )
sign = -1 ;
if (derivate_second_order[i]*derivate_second_order[i+1] == 0 )
sign = 0 ;
if ( sign==0 || sign==-1 )
- if ( 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 )
{
for (unsigned int i=0;i<max_cross_list.size();i++)
distance.push_back(fabs(max_cross_list[i]+1-center));
- float distance_min=*min_element(distance.begin(),distance.end()) ;
+ double distance_min=*min_element(distance.begin(),distance.end()) ;
int pos_distance_min=find(distance.begin(),distance.end(),distance_min)-distance.begin();
int best_place = max_cross_list[pos_distance_min] + part[id].first ;
if ( max_cross >=0 )
}
derivate_second_order.clear() ;
gradient_absolute.clear() ;
- signe_change.clear() ;
zero_cross.clear() ;
edge.clear() ;
max_cross_list.clear() ;
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<InternalPatch> DealWithNoCut(const InternalPatch *patch)
return ret;
}
-void DealWithCut(const InternalPatch *patchToBeSplit, int axisId, int cutPlace, std::vector<MEDCouplingAutoRefCountObjectPtr<InternalPatch> >& listOfPatches)
+void DealWithCut(double minPatchLgth, const InternalPatch *patchToBeSplit, int axisId, int cutPlace, std::vector<MEDCouplingAutoRefCountObjectPtr<InternalPatch> >& listOfPatches)
{
MEDCouplingAutoRefCountObjectPtr<InternalPatch> leftPart,rightPart;
std::vector< std::pair<int,int> > rect(patchToBeSplit->getConstPart());
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);
}
}
/*!
- * 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<bool>& criterion, const std::vector<int>& factors)
{
std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatches,listOfPatchesOK;
//
MEDCouplingAutoRefCountObjectPtr<InternalPatch> 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())
for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::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()<bso.getMaxCells())
- { listOfPatchesOK.push_back(DealWithNoCut(*it)); continue; }//action 1
- FindHole(bso,*it,axisId,cutFound,cutPlace);//axisId overwritten here if FindHole equal to true !
- if(cutFound)
- { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 2
- FindInflection(bso,*it,cutFound,cutPlace,axisId);//axisId overwritten here if cutFound equal to true !
- if(cutFound)
- { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 3
- TryAction4(bso,*it,axisId,rangeOfAxisId,cutFound,cutPlace);
- if(cutFound)
- { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 4
- listOfPatchesOK.push_back(DealWithNoCut(*it));
+ int axisId,largestLength,cutPlace;
+ MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt((*it)->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;
}
* \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<int>& st, const std::vector<bool>& crit, std::vector<bool>& reducedCrit, std::vector< std::pair<int,int> >& partCompactFormat)
+int MEDCouplingStructuredMesh::FindMinimalPartOf(int minPatchLgth, const std::vector<int>& st, const std::vector<bool>& crit, std::vector<bool>& reducedCrit, std::vector< std::pair<int,int> >& 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);
default:
throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindMinimalPartOf : only dimension 1, 2 and 3 are supported actually !");
}
+ std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(partCompactFormat));
+ int i(0);
+ for(std::vector< std::pair<int,int> >::iterator it=partCompactFormat.begin();it!=partCompactFormat.end();it++,i++)
+ {
+ if(st[i]<minPatchLgth)
+ throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindMinimalPartOf : the input patch is tinier than the min length constraint !");
+ int start((*it).first),stop((*it).second),middle((start+stop)/2);
+ if(stop-start<minPatchLgth)
+ {
+ (*it).first=middle-minPatchLgth/2;
+ (*it).second=middle+minPatchLgth-minPatchLgth/2;
+ if((*it).first<0)
+ {
+ (*it).second+=-(*it).first;
+ (*it).first=0;
+ }
+ if((*it).second>st[i])
+ {
+ (*it).first-=(*it).second-st[i];
+ (*it).second=st[i];
+ }
+ }
+ }
ExtractFieldOfBoolFrom(st,crit,partCompactFormat,reducedCrit);
return ret;
}
}
}
if(ret==0)
- return ret;
+ {
+ std::size_t sz(st.size());
+ partCompactFormat.resize(sz);
+ for(std::size_t i=0;i<sz;i++)
+ {
+ partCompactFormat[i].first=st[i]/2;
+ partCompactFormat[i].second=st[i]/2;
+ }
+ return ret;
+ }
partCompactFormat.resize(1);
partCompactFormat[0].first=nxMin; partCompactFormat[0].second=nxMax+1;
return ret;
}
}
if(ret==0)
- return ret;
+ {
+ std::size_t sz(st.size());
+ partCompactFormat.resize(sz);
+ for(std::size_t i=0;i<sz;i++)
+ {
+ partCompactFormat[i].first=st[i]/2;
+ partCompactFormat[i].second=st[i]/2;
+ }
+ return ret;
+ }
partCompactFormat.resize(2);
partCompactFormat[0].first=nxMin; partCompactFormat[0].second=nxMax+1;
partCompactFormat[1].first=nyMin; partCompactFormat[1].second=nyMax+1;
}
}
if(ret==0)
- return ret;
+ {
+ std::size_t sz(st.size());
+ partCompactFormat.resize(sz);
+ for(std::size_t i=0;i<sz;i++)
+ {
+ partCompactFormat[i].first=st[i]/2;
+ partCompactFormat[i].second=st[i]/2;
+ }
+ return ret;
+ }
partCompactFormat.resize(3);
partCompactFormat[0].first=nxMin; partCompactFormat[0].second=nxMax+1;
partCompactFormat[1].first=nyMin; partCompactFormat[1].second=nyMax+1;
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<int> 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<dim;i++)
+ {
+ minAx=std::min(minAx,cgs[i]);
+ maxAx=std::max(maxAx,cgs[i]);
+ }
+ return (double)minAx/(double)maxAx;
+ }
+ throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::computeSquareness : only dimension 2 and 3 supported !");
+}
+
/*!
* Given a struct \a strct it returns a split vector [1,strct[0],strct[0]*strct[1]...]
* This decomposition allows to quickly find i,j,k given a global id.
MEDCOUPLING_EXPORT virtual void getSplitNodeValues(int *res) const;
MEDCOUPLING_EXPORT virtual std::vector<int> getNodeGridStructure() const = 0;
MEDCOUPLING_EXPORT std::vector<int> getCellGridStructure() const;
+ MEDCOUPLING_EXPORT double computeSquareness() const;
MEDCOUPLING_EXPORT virtual MEDCouplingStructuredMesh *buildStructuredSubPart(const std::vector< std::pair<int,int> >& cellPart) const = 0;
MEDCOUPLING_EXPORT static std::vector<int> GetSplitVectFromStruct(const std::vector<int>& strct);
MEDCOUPLING_EXPORT static bool IsPartStructured(const int *startIds, const int *stopIds, const std::vector<int>& st, std::vector< std::pair<int,int> >& partCompactFormat);
MEDCOUPLING_EXPORT static int DeduceNumberOfGivenRangeInCompactFrmt(const std::vector< std::pair<int,int> >& partCompactFormat);
MEDCOUPLING_EXPORT static int DeduceNumberOfGivenStructure(const std::vector<int>& st);
MEDCOUPLING_EXPORT static void FindTheWidestAxisOfGivenRangeInCompactFrmt(const std::vector< std::pair<int,int> >& partCompactFormat, int& axisId, int& sizeOfRange);
- MEDCOUPLING_EXPORT static int FindMinimalPartOf(const std::vector<int>& st, const std::vector<bool>& crit, std::vector<bool>& reducedCrit, std::vector< std::pair<int,int> >& partCompactFormat);
+ MEDCOUPLING_EXPORT static int FindMinimalPartOf(int minPatchLgth, const std::vector<int>& st, const std::vector<bool>& crit, std::vector<bool>& reducedCrit, std::vector< std::pair<int,int> >& partCompactFormat);
MEDCOUPLING_EXPORT static std::vector< std::vector<int> > ComputeSignaturePerAxisOf(const std::vector<int>& st, const std::vector<bool>& crit);
private:
static int GetNumberOfCellsOfSubLevelMesh(const std::vector<int>& cgs, int mdim);
# 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)])
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()
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
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<int> getNodeGridStructure() const throw(INTERP_KERNEL::Exception);
std::vector<int> getCellGridStructure() const throw(INTERP_KERNEL::Exception);
MEDCoupling1SGTUMesh *build1SGTUnstructured() const throw(INTERP_KERNEL::Exception);
