1 // Copyright (C) 2007-2014 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay
21 #include "MEDCouplingCartesianAMRMesh.hxx"
22 #include "MEDCoupling1GTUMesh.hxx"
23 #include "MEDCouplingIMesh.hxx"
24 #include "MEDCouplingUMesh.hxx"
30 using namespace ParaMEDMEM;
35 * \param [in] mesh not null pointer of refined mesh replacing the cell range of \a father defined by the bottom left and top right just after.
36 * \param [in] bottomLeftTopRight a vector equal to the space dimension of \a mesh that specifies for each dimension, the included cell start of the range for the first element of the pair,
37 * a the end cell (\b excluded) of the range for the second element of the pair.
39 MEDCouplingCartesianAMRPatch::MEDCouplingCartesianAMRPatch(MEDCouplingCartesianAMRMesh *mesh, const std::vector< std::pair<int,int> >& bottomLeftTopRight)
42 throw INTERP_KERNEL::Exception("EDCouplingCartesianAMRPatch constructor : input mesh is NULL !");
43 _mesh=mesh; _mesh->incrRef();
44 int dim((int)bottomLeftTopRight.size()),dimExp(_mesh->getSpaceDimension());
46 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch constructor : space dimension of father and input bottomLeft/topRight size mismatches !");
47 _bl_tr=bottomLeftTopRight;
50 int MEDCouplingCartesianAMRPatch::getNumberOfCellsRecursiveWithOverlap() const
52 return _mesh->getNumberOfCellsRecursiveWithOverlap();
55 int MEDCouplingCartesianAMRPatch::getNumberOfCellsRecursiveWithoutOverlap() const
57 return _mesh->getNumberOfCellsRecursiveWithoutOverlap();
60 int MEDCouplingCartesianAMRPatch::getMaxNumberOfLevelsRelativeToThis() const
62 return _mesh->getMaxNumberOfLevelsRelativeToThis();
65 void MEDCouplingCartesianAMRPatch::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
67 return _mesh->addPatch(bottomLeftTopRight,factors);
70 int MEDCouplingCartesianAMRPatch::getNumberOfOverlapedCellsForFather() const
72 return MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(_bl_tr);
76 * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
79 * \param [in] other - The other patch
80 * \param [in] ghostLev - The size of the neighborhood zone.
82 * \throw if \a this or \a other are invalid (end before start).
83 * \throw if \a ghostLev is \b not >= 0 .
84 * \throw if \a this and \a other have not the same space dimension.
86 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhood(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const
89 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the size of the neighborhood must be >= 0 !");
91 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the input patch is NULL !");
92 const std::vector< std::pair<int,int> >& thisp(getBLTRRange());
93 const std::vector< std::pair<int,int> >& otherp(other->getBLTRRange());
94 std::size_t thispsize(thisp.size());
95 if(thispsize!=otherp.size())
96 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the dimensions must be the same !");
97 for(std::size_t i=0;i<thispsize;i++)
99 const std::pair<int,int>& thispp(thisp[i]);
100 const std::pair<int,int>& otherpp(otherp[i]);
101 if(thispp.second<thispp.first)
102 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : this patch is invalid !");
103 if(otherpp.second<otherpp.first)
104 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : this patch is invalid !");
105 if(otherpp.first==thispp.second+ghostLev-1)
107 if(otherpp.second+ghostLev-1==thispp.first)
114 std::size_t MEDCouplingCartesianAMRPatch::getHeapMemorySizeWithoutChildren() const
116 std::size_t ret(sizeof(MEDCouplingCartesianAMRPatch));
117 ret+=_bl_tr.capacity()*sizeof(std::pair<int,int>);
121 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRPatch::getDirectChildren() const
123 std::vector<const BigMemoryObject *> ret;
124 if((const MEDCouplingCartesianAMRMesh *)_mesh)
125 ret.push_back((const MEDCouplingCartesianAMRMesh *)_mesh);
132 MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::New(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
133 const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop)
135 return new MEDCouplingCartesianAMRMesh(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
138 int MEDCouplingCartesianAMRMesh::getSpaceDimension() const
140 return _mesh->getSpaceDimension();
143 void MEDCouplingCartesianAMRMesh::setFactors(const std::vector<int>& newFactors)
145 if(getSpaceDimension()!=(int)newFactors.size())
146 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::setFactors : size of input factors is not equal to the space dimension !");
152 if(_factors==newFactors)
154 if(!_patches.empty())
155 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::setFactors : modification of factors is not allowed when presence of patches !");
159 int MEDCouplingCartesianAMRMesh::getMaxNumberOfLevelsRelativeToThis() const
162 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
163 ret=std::max(ret,(*it)->getMaxNumberOfLevelsRelativeToThis()+1);
167 int MEDCouplingCartesianAMRMesh::getNumberOfCellsAtCurrentLevel() const
169 return _mesh->getNumberOfCells();
172 int MEDCouplingCartesianAMRMesh::getNumberOfCellsRecursiveWithOverlap() const
174 int ret(_mesh->getNumberOfCells());
175 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
177 ret+=(*it)->getNumberOfCellsRecursiveWithOverlap();
182 int MEDCouplingCartesianAMRMesh::getNumberOfCellsRecursiveWithoutOverlap() const
184 int ret(_mesh->getNumberOfCells());
185 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
187 ret-=(*it)->getNumberOfOverlapedCellsForFather();
188 ret+=(*it)->getNumberOfCellsRecursiveWithoutOverlap();
193 const MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::getFather() const
198 const MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::getGodFather() const
203 return _father->getGodFather();
206 void MEDCouplingCartesianAMRMesh::detachFromFather()
212 * \param [in] bottomLeftTopRight a vector equal to the space dimension of \a mesh that specifies for each dimension, the included cell start of the range for the first element of the pair,
213 * a the end cell (\b excluded) of the range for the second element of the pair.
214 * \param [in] factors The factor of refinement per axis (different from 0).
216 void MEDCouplingCartesianAMRMesh::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
218 checkFactorsAndIfNotSetAssign(factors);
219 MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> mesh(static_cast<MEDCouplingIMesh *>(_mesh->buildStructuredSubPart(bottomLeftTopRight)));
220 mesh->refineWithFactor(factors);
221 MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRMesh> zeMesh(new MEDCouplingCartesianAMRMesh(this,mesh));
222 MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> elt(new MEDCouplingCartesianAMRPatch(zeMesh,bottomLeftTopRight));
223 _patches.push_back(elt);
228 class InternalPatch : public RefCountObjectOnly
231 InternalPatch():_nb_of_true(0) { }
232 int getDimension() const { return (int)_part.size(); }
233 double getEfficiency() const { return (double)_nb_of_true/(double)_crit.size(); }
234 int getNumberOfCells() const { return (int)_crit.size(); }
235 void setNumberOfTrue(int nboft) { _nb_of_true=nboft; }
236 std::vector<bool>& getCriterion() { return _crit; }
237 const std::vector<bool>& getConstCriterion() const { return _crit; }
238 void setPart(const std::vector< std::pair<int,int> >& part) { _part=part; }
239 std::vector< std::pair<int,int> >& getPart() { return _part; }
240 const std::vector< std::pair<int,int> >& getConstPart() const { return _part; }
241 bool presenceOfTrue() const { return _nb_of_true>0; }
242 std::vector<int> computeCGS() const { return MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(_part); }
243 std::vector< std::vector<int> > computeSignature() const { return MEDCouplingStructuredMesh::ComputeSignaturePerAxisOf(computeCGS(),getConstCriterion()); }
244 double getEfficiencyPerAxis(int axisId) const { return (double)_nb_of_true/((double)(_part[axisId].second-_part[axisId].first)); }
245 void zipToFitOnCriterion();
246 void updateNumberOfTrue() const;
247 MEDCouplingAutoRefCountObjectPtr<InternalPatch> extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const;
248 MEDCouplingAutoRefCountObjectPtr<InternalPatch> deepCpy() const;
252 mutable int _nb_of_true;
253 std::vector<bool> _crit;
255 std::vector< std::pair<int,int> > _part;
258 void InternalPatch::zipToFitOnCriterion()
260 std::vector<int> cgs(computeCGS());
261 std::vector<bool> newCrit;
262 std::vector< std::pair<int,int> > newPart,newPart2;
263 int newNbOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,_crit,newCrit,newPart));
264 MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(_part,newPart,newPart2);
265 if(newNbOfTrue!=_nb_of_true)
266 throw INTERP_KERNEL::Exception("InternalPatch::zipToFitOnCrit : internal error !");
267 _crit=newCrit; _part=newPart2;
270 void InternalPatch::updateNumberOfTrue() const
272 _nb_of_true=(int)std::count(_crit.begin(),_crit.end(),true);
275 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const
277 MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
278 std::vector<int> cgs(computeCGS());
279 std::vector< std::pair<int,int> > newPart;
280 MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(_part,partInGlobal,newPart);
281 MEDCouplingStructuredMesh::ExtractFieldOfBoolFrom(cgs,_crit,newPart,ret->getCriterion());
282 ret->setPart(partInGlobal);
283 ret->updateNumberOfTrue();
287 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::deepCpy() const
289 MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
294 void DissectBigPatch(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
296 cutFound=false; cutPlace=-1;
297 std::vector<double> ratio(rangeOfAxisId-1);
298 for(int id=0;id<rangeOfAxisId-1;id++)
300 double efficiency[2];
303 std::vector< std::pair<int,int> > rectH(patchToBeSplit->getConstPart());
305 rectH[axisId].second=patchToBeSplit->getConstPart()[axisId].first+id;
307 rectH[axisId].first=patchToBeSplit->getConstPart()[axisId].first+id;
308 MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(patchToBeSplit->deepCpy());
309 p->zipToFitOnCriterion();
311 efficiency[h]=p->getEfficiencyPerAxis(axisId);
313 ratio[id]=std::max(efficiency[0],efficiency[1])/std::min(efficiency[0],efficiency[1]);
315 int minCellDirection(bso.getMinCellDirection()),indexMin(-1);
316 int dimRatioInner(rangeOfAxisId-1-2*(minCellDirection-1));
317 std::vector<double> ratio_inner(dimRatioInner);
318 double minRatio(1.e10);
319 for(int i=0; i<dimRatioInner; i++)
321 if(ratio[minCellDirection-1+i]<minRatio)
323 minRatio=ratio[minCellDirection-1+i];
324 indexMin=i+minCellDirection;
327 cutFound=true; cutPlace=indexMin+patchToBeSplit->getConstPart()[axisId].first-1;
330 void FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int& axisId, bool& cutFound, int& cutPlace)
332 cutPlace=-1; cutFound=false;
333 int minCellDirection(bso.getMinCellDirection());
334 const int dim(patchToBeSplit->getDimension());
335 std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
336 for(int id=0;id<dim;id++)
338 const std::vector<int>& signature(signatures[id]);
339 std::vector<int> hole;
340 std::vector<double> distance;
341 int len((int)signature.size());
342 for(int i=0;i<len;i++)
344 if(len>= 2*minCellDirection && i >= minCellDirection-1 && i <= len-minCellDirection-1)
348 double center(((double)len/2.));
349 for(std::size_t i=0;i<hole.size();i++)
350 distance.push_back(fabs(hole[i]+1.-center));
352 std::size_t posDistanceMin(std::distance(distance.begin(),std::min_element(distance.begin(),distance.end())));
355 cutPlace=hole[posDistanceMin]+patchToBeSplit->getConstPart()[axisId].first+1;
361 void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, bool& cutFound, int& cutPlace, int& axisId)
363 cutFound=false; cutPlace=-1;// do not set axisId before to be sure that cutFound was set to true
365 const std::vector< std::pair<int,int> >& part(patchToBeSplit->getConstPart());
366 int sign,minCellDirection(bso.getMinCellDirection());
367 const int dim(patchToBeSplit->getDimension());
369 std::vector<int> zeroCrossDims(dim,-1);
370 std::vector<int> zeroCrossVals(dim,-1);
371 std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
372 for (int id=0;id<dim;id++)
374 const std::vector<int>& signature(signatures[id]);
376 std::vector<int> derivate_second_order,gradient_absolute,signe_change,zero_cross,edge,max_cross_list ;
377 std::vector<double> distance ;
379 for (unsigned int i=1;i<signature.size()-1;i++)
380 derivate_second_order.push_back(signature[i-1]-2*signature[i]+signature[i+1]) ;
382 // Gradient absolute value
383 for ( unsigned int i=1;i<derivate_second_order.size();i++)
384 gradient_absolute.push_back(fabs(derivate_second_order[i]-derivate_second_order[i-1])) ;
385 if(derivate_second_order.empty())
387 for (unsigned int i=0;i<derivate_second_order.size()-1;i++)
389 if (derivate_second_order[i]*derivate_second_order[i+1] < 0 )
391 if (derivate_second_order[i]*derivate_second_order[i+1] > 0 )
393 if (derivate_second_order[i]*derivate_second_order[i+1] == 0 )
395 if ( sign==0 || sign==-1 )
396 if ( i >= (unsigned int)minCellDirection-2 && i <= signature.size()-minCellDirection-2 )
398 zero_cross.push_back(i) ;
399 edge.push_back(gradient_absolute[i]) ;
401 signe_change.push_back(sign) ;
403 if ( zero_cross.size() > 0 )
405 int max_cross=*max_element(edge.begin(),edge.end()) ;
406 for (unsigned int i=0;i<edge.size();i++)
407 if (edge[i]==max_cross)
408 max_cross_list.push_back(zero_cross[i]+1) ;
410 double center((signature.size()/2.0));
411 for (unsigned int i=0;i<max_cross_list.size();i++)
412 distance.push_back(fabs(max_cross_list[i]+1-center));
414 float distance_min=*min_element(distance.begin(),distance.end()) ;
415 int pos_distance_min=find(distance.begin(),distance.end(),distance_min)-distance.begin();
416 int best_place = max_cross_list[pos_distance_min] + part[id].first ;
419 zeroCrossDims[id] = best_place ;
420 zeroCrossVals[id] = max_cross ;
423 derivate_second_order.clear() ;
424 gradient_absolute.clear() ;
425 signe_change.clear() ;
428 max_cross_list.clear() ;
432 if ( zeroCrossDims[0]!=-1 || zeroCrossDims[1]!=-1 )
434 int max_cross_dims = *max_element(zeroCrossVals.begin(),zeroCrossVals.end()) ;
436 if (zeroCrossVals[0]==max_cross_dims && zeroCrossVals[1]==max_cross_dims )
438 int nl_left(part[0].second-part[0].first);
439 int nc_left(part[1].second-part[1].first);
440 if ( nl_left >= nc_left )
446 max_cross_dims=std::find(zeroCrossVals.begin(),zeroCrossVals.end(),max_cross_dims)-zeroCrossVals.begin();
448 cutPlace=zeroCrossDims[max_cross_dims];
449 axisId=max_cross_dims ;
453 void TryAction4(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
456 if(patchToBeSplit->getEfficiency()<=bso.getEffeciencySnd())
458 if(rangeOfAxisId>=2*bso.getMinCellDirection())
461 cutPlace=rangeOfAxisId/2+patchToBeSplit->getConstPart()[axisId].first-1;
466 if(patchToBeSplit->getNumberOfCells()>bso.getMaxCells())
468 DissectBigPatch(bso,patchToBeSplit,axisId,rangeOfAxisId,cutFound,cutPlace);
473 MEDCouplingAutoRefCountObjectPtr<InternalPatch> DealWithNoCut(const InternalPatch *patch)
475 MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(const_cast<InternalPatch *>(patch));
480 void DealWithCut(const InternalPatch *patchToBeSplit, int axisId, int cutPlace, std::vector<MEDCouplingAutoRefCountObjectPtr<InternalPatch> >& listOfPatches)
482 MEDCouplingAutoRefCountObjectPtr<InternalPatch> leftPart,rightPart;
483 std::vector< std::pair<int,int> > rect(patchToBeSplit->getConstPart());
484 std::vector< std::pair<int,int> > leftRect(rect),rightRect(rect);
485 leftRect[axisId].second=cutPlace+1;
486 rightRect[axisId].first=cutPlace+1;
487 leftPart=patchToBeSplit->extractPart(leftRect);
488 rightPart=patchToBeSplit->extractPart(rightRect);
489 leftPart->zipToFitOnCriterion(); rightPart->zipToFitOnCriterion();
490 listOfPatches.push_back(leftPart);
491 listOfPatches.push_back(rightPart);
497 * 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.
499 void MEDCouplingCartesianAMRMesh::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const std::vector<bool>& criterion, const std::vector<int>& factors)
501 int nbCells(getNumberOfCellsAtCurrentLevel());
502 if(nbCells!=(int)criterion.size())
503 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterion : the number of tuples of criterion array must be equal to the number of cells at the current level !");
505 std::vector<int> cgs(_mesh->getCellGridStructure());
506 std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatches,listOfPatchesOK;
508 MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(new InternalPatch);
509 p->setNumberOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,criterion,p->getCriterion(),p->getPart()));
510 if(p->presenceOfTrue())
511 listOfPatches.push_back(p);
512 while(!listOfPatches.empty())
514 std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatchesTmp;
515 for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::iterator it=listOfPatches.begin();it!=listOfPatches.end();it++)
518 int axisId,rangeOfAxisId,cutPlace;
520 MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt((*it)->getConstPart(),axisId,rangeOfAxisId);
521 if((*it)->getEfficiency()>=bso.getEffeciency() && (*it)->getNumberOfCells()<bso.getMaxCells())
522 { listOfPatchesOK.push_back(DealWithNoCut(*it)); continue; }//action 1
523 FindHole(bso,*it,axisId,cutFound,cutPlace);//axisId overwritten here if FindHole equal to true !
525 { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 2
526 FindInflection(bso,*it,cutFound,cutPlace,axisId);//axisId overwritten here if cutFound equal to true !
528 { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 3
529 TryAction4(bso,*it,axisId,rangeOfAxisId,cutFound,cutPlace);
531 { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 4
532 listOfPatchesOK.push_back(DealWithNoCut(*it));
534 listOfPatches=listOfPatchesTmp;
536 for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::const_iterator it=listOfPatchesOK.begin();it!=listOfPatchesOK.end();it++)
537 addPatch((*it)->getConstPart(),factors);
541 * 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.
543 void MEDCouplingCartesianAMRMesh::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const DataArrayByte *criterion, const std::vector<int>& factors)
545 if(!criterion || !criterion->isAllocated())
546 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterion : the criterion DataArrayByte instance must be allocated and not NULL !");
547 std::vector<bool> crit(criterion->toVectorOfBool());//check that criterion has one component.
548 createPatchesFromCriterion(bso,crit,factors);
551 void MEDCouplingCartesianAMRMesh::removeAllPatches()
557 void MEDCouplingCartesianAMRMesh::removePatch(int patchId)
559 checkPatchId(patchId);
560 int sz((int)_patches.size()),j(0);
561 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> > patches(sz-1);
562 for(int i=0;i<sz;i++)
564 patches[j++]=_patches[i];
565 (const_cast<MEDCouplingCartesianAMRMesh *>(_patches[patchId]->getMesh()))->detachFromFather();
570 int MEDCouplingCartesianAMRMesh::getNumberOfPatches() const
572 return (int)_patches.size();
575 const MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRMesh::getPatch(int patchId) const
577 checkPatchId(patchId);
578 return _patches[patchId];
582 * This method states if patch2 (with id \a patchId2) is in the neighborhood of patch1 (with id \a patchId1).
583 * The neighborhood size is defined by \a ghostLev in the reference of \a this ( \b not in the reference of patches !).
585 bool MEDCouplingCartesianAMRMesh::isPatchInNeighborhoodOf(int patchId1, int patchId2, int ghostLev) const
588 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::isPatchInNeighborhoodOf : the ghost size must be >=0 !");
589 const MEDCouplingCartesianAMRPatch *p1(getPatch(patchId1)),*p2(getPatch(patchId2));
591 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::isPatchInNeighborhoodOf : no factors defined !");
592 int ghostLevInPatchRef;
594 ghostLevInPatchRef=0;
597 ghostLevInPatchRef=(ghostLev-1)/_factors[0]+1;
598 for(std::size_t i=0;i<_factors.size();i++)
599 ghostLevInPatchRef=std::max(ghostLevInPatchRef,_factors[i]/ghostLev+1);
601 return p1->isInMyNeighborhood(p2,ghostLevInPatchRef);
605 * This method creates a new cell field array on given \a patchId patch in \a this starting from a coarse cell field on \a this \a cellFieldOnThis.
606 * This method can be seen as a fast projection from the cell field \a cellFieldOnThis on \c this->getImageMesh() to a refined part of \a this
607 * defined by the patch with id \a patchId.
609 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
610 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
611 * \return DataArrayDouble * - The array of the cell field on the requested patch
613 * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
614 * \throw if \a cellFieldOnThis is NULL or not allocated
615 * \sa fillCellFieldOnPatch, MEDCouplingIMesh::SpreadCoarseToFine
617 DataArrayDouble *MEDCouplingCartesianAMRMesh::createCellFieldOnPatch(int patchId, const DataArrayDouble *cellFieldOnThis) const
619 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
620 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
621 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
622 const MEDCouplingIMesh *fine(patch->getMesh()->getImageMesh());
623 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(fine->getNumberOfCells(),cellFieldOnThis->getNumberOfComponents());
624 ret->copyStringInfoFrom(*cellFieldOnThis);
625 MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),ret,patch->getBLTRRange(),getFactors());
630 * This method is equivalent to MEDCouplingCartesianAMRMesh::createCellFieldOnPatch except that here instead of creating a new instance
631 * it fills the value into the \a cellFieldOnPatch data.
633 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
634 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
635 * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
637 * \sa createCellFieldOnPatch, fillCellFieldComingFromPatch
639 void MEDCouplingCartesianAMRMesh::fillCellFieldOnPatch(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch) const
641 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
642 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
643 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
644 MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors());
648 * This method is the generalization of fillCellFieldOnPatch method. This method only projects coarse to fine without considering the
649 * potential neighbor patches covered by the ghost cells of patch with id \a patchId.
651 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
652 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
653 * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
654 * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
656 * \sa fillCellFieldOnPatch, fillCellFieldOnPatchGhostAdv
658 void MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchGhost(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, int ghostLev) const
660 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
661 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatchGhost : the input cell field array is NULL or not allocated !");
662 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
663 MEDCouplingIMesh::SpreadCoarseToFineGhost(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),ghostLev);
667 * This method is a refinement of fillCellFieldOnPatchGhost. fillCellFieldOnPatchGhost is first called.
668 * Then for all other patches than those pointed by \a patchId that overlap the ghost zone of the patch impact the ghost zone adequately.
670 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
671 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
672 * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
673 * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
674 * \param [in] arrsOnPatches - \b WARNING arrsOnPatches[patchId] is \b NOT \b const. All others are const.
676 void MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchGhostAdv(int patchId, const DataArrayDouble *cellFieldOnThis, int ghostLev, const std::vector<const DataArrayDouble *>& arrsOnPatches) const
678 int nbp(getNumberOfPatches());
679 if(nbp!=(int)arrsOnPatches.size())
681 std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchGhostAdv : there are " << nbp << " patches in this and " << arrsOnPatches.size() << " arrays in the last parameter !";
682 throw INTERP_KERNEL::Exception(oss.str().c_str());
684 // first, do as usual
685 fillCellFieldOnPatchGhost(patchId,cellFieldOnThis,const_cast<DataArrayDouble *>(arrsOnPatches[patchId]),ghostLev);
687 const MEDCouplingCartesianAMRPatch *refP(getPatch(patchId));
688 const std::vector< std::pair<int,int> >& refBLTR(refP->getBLTRRange());
689 for(int i=0;i<nbp;i++)
692 if(isPatchInNeighborhoodOf(i,patchId,ghostLev))
694 const MEDCouplingCartesianAMRPatch *otherP(getPatch(i));
695 const std::vector< std::pair<int,int> >& otherBLTR(otherP->getBLTRRange());
696 std::vector< std::pair<int,int> > tmp0;
697 MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(refBLTR,otherBLTR,tmp0);
698 ApplyFactorsOnCompactFrmt(tmp0,_factors);
699 ApplyGhostOnCompactFrmt(tmp0,ghostLev);
700 //std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
706 * This method updates \a cellFieldOnThis part of values coming from the cell field \a cellFieldOnPatch lying on patch having id \a patchId.
708 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
709 * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
710 * \param [in,out] cellFieldOnThis The array of the cell field on \a this to be updated only on the part concerning the patch with id \a patchId.
712 * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
713 * \throw if \a cellFieldOnPatch is NULL or not allocated
714 * \sa createCellFieldOnPatch, MEDCouplingIMesh::CondenseFineToCoarse,fillCellFieldComingFromPatchGhost
716 void MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch(int patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis) const
718 if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
719 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch : the input cell field array is NULL or not allocated !");
720 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
721 MEDCouplingIMesh::CondenseFineToCoarse(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis);
725 * This method is the extension of MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch managing the ghost cells. If this
726 * method is called with \a ghostLev equal to 0 it behaves exactly as MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch.
728 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
729 * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
730 * \param [in,out] cellFieldOnThis The array of the cell field on \a this to be updated only on the part concerning the patch with id \a patchId.
731 * \param [in] ghostLev The size of ghost zone (must be >= 0 !)
733 * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
734 * \throw if \a cellFieldOnPatch is NULL or not allocated
735 * \sa fillCellFieldComingFromPatch
737 void MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatchGhost(int patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis, int ghostLev) const
739 if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
740 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatchGhost : the input cell field array is NULL or not allocated !");
741 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
742 MEDCouplingIMesh::CondenseFineToCoarseGhost(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis,ghostLev);
746 * This method finds all patches (located by their ids) that are in the neighborhood of patch with id \a patchId. The neighborhood size is
747 * defined by ghostLev.
749 * \param [in] patchId - the id of the considered patch.
750 * \param [in] ghostLev - the size of the neighborhood.
751 * \return DataArrayInt * - the newly allocated array containing the list of patches in the neighborhood of the considered patch. This array is to be deallocated by the caller.
753 DataArrayInt *MEDCouplingCartesianAMRMesh::findPatchesInTheNeighborhoodOf(int patchId, int ghostLev) const
755 int nbp(getNumberOfPatches());
756 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
757 for(int i=0;i<nbp;i++)
760 if(isPatchInNeighborhoodOf(i,patchId,ghostLev))
761 ret->pushBackSilent(i);
766 MEDCouplingUMesh *MEDCouplingCartesianAMRMesh::buildUnstructured() const
768 MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> part(_mesh->buildUnstructured());
769 std::vector<bool> bs(_mesh->getNumberOfCells(),false);
770 std::vector<int> cgs(_mesh->getCellGridStructure());
771 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > msSafe(_patches.size()+1);
773 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
775 MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
776 msSafe[ii+1]=(*it)->getMesh()->buildUnstructured();
778 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsOff(DataArrayInt::BuildListOfSwitchedOff(bs));
779 msSafe[0]=static_cast<MEDCouplingUMesh *>(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
780 std::vector< const MEDCouplingUMesh * > ms(msSafe.size());
781 for(std::size_t i=0;i<msSafe.size();i++)
783 return MEDCouplingUMesh::MergeUMeshes(ms);
787 * This method returns a mesh containing as cells that there is patches at the current level.
788 * The patches are seen like 'boxes' that is too say the refinement will not appear here.
790 * \return MEDCoupling1SGTUMesh * - A new object to be managed by the caller containing as cells as there are patches in \a this.
792 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMesh::buildMeshFromPatchEnvelop() const
794 std::vector<const MEDCoupling1SGTUMesh *> cells;
795 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> > cellsSafe;
796 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
798 const MEDCouplingCartesianAMRPatch *patch(*it);
801 MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> cell(patch->getMesh()->getImageMesh()->asSingleCell());
802 MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> cell1SGT(cell->build1SGTUnstructured());
803 cellsSafe.push_back(cell1SGT); cells.push_back(cell1SGT);
806 return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(cells);
809 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
810 const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):_father(0)
812 _mesh=MEDCouplingIMesh::New(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
815 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(MEDCouplingCartesianAMRMesh *father, MEDCouplingIMesh *mesh):_father(father)
818 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh(MEDCouplingIMesh *mesh) constructor : empty father !");
820 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh(MEDCouplingIMesh *mesh) constructor : The input mesh is null !");
821 mesh->checkCoherency();
822 _mesh=mesh; _mesh->incrRef();
825 void MEDCouplingCartesianAMRMesh::checkPatchId(int patchId) const
827 int sz(getNumberOfPatches());
828 if(patchId<0 || patchId>=sz)
830 std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::checkPatchId : invalid patchId (" << patchId << ") ! Must be in [0," << sz << ") !";
831 throw INTERP_KERNEL::Exception(oss.str().c_str());
835 void MEDCouplingCartesianAMRMesh::checkFactorsAndIfNotSetAssign(const std::vector<int>& factors)
837 if(getSpaceDimension()!=(int)factors.size())
838 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::checkFactorsAndIfNotSetAssign : invalid size of factors ! size must be equal to the spaceDimension !");
845 if(_factors!=factors)
846 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::checkFactorsAndIfNotSetAssign : the factors ");
851 * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in refined reference.
852 * \param [in] factors - the factors per axis.
854 void MEDCouplingCartesianAMRMesh::ApplyFactorsOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, const std::vector<int>& factors)
856 std::size_t sz(factors.size());
857 if(sz!=partBeforeFact.size())
858 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::ApplyFactorsOnCompactFrmt : size of input vectors must be the same !");
859 for(std::size_t i=0;i<sz;i++)
861 partBeforeFact[i].first*=factors[i];
862 partBeforeFact[i].second*=factors[i];
867 * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in ghost reference.
868 * \param [in] ghostSize - the ghost size of zone for all axis.
870 void MEDCouplingCartesianAMRMesh::ApplyGhostOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, int ghostSize)
873 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::ApplyGhostOnCompactFrmt : ghost size must be >= 0 !");
874 std::size_t sz(partBeforeFact.size());
875 for(std::size_t i=0;i<sz;i++)
877 partBeforeFact[i].first+=ghostSize;
878 partBeforeFact[i].second+=ghostSize;
882 std::size_t MEDCouplingCartesianAMRMesh::getHeapMemorySizeWithoutChildren() const
884 return sizeof(MEDCouplingCartesianAMRMesh);
887 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMesh::getDirectChildren() const
889 std::vector<const BigMemoryObject *> ret;
890 if((const MEDCouplingIMesh *)_mesh)
891 ret.push_back((const MEDCouplingIMesh *)_mesh);
892 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
894 if((const MEDCouplingCartesianAMRPatch*)*it)
895 ret.push_back((const MEDCouplingCartesianAMRPatch*)*it);
900 void MEDCouplingCartesianAMRMesh::updateTime() const
902 if((const MEDCouplingIMesh *)_mesh)
903 updateTimeWith(*_mesh);
904 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
906 const MEDCouplingCartesianAMRPatch *elt(*it);
909 const MEDCouplingCartesianAMRMesh *mesh(elt->getMesh());
911 updateTimeWith(*mesh);