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);
75 std::size_t MEDCouplingCartesianAMRPatch::getHeapMemorySizeWithoutChildren() const
77 std::size_t ret(sizeof(MEDCouplingCartesianAMRPatch));
78 ret+=_bl_tr.capacity()*sizeof(std::pair<int,int>);
82 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRPatch::getDirectChildren() const
84 std::vector<const BigMemoryObject *> ret;
85 if((const MEDCouplingCartesianAMRMesh *)_mesh)
86 ret.push_back((const MEDCouplingCartesianAMRMesh *)_mesh);
93 MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::New(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
94 const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop)
96 return new MEDCouplingCartesianAMRMesh(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
99 int MEDCouplingCartesianAMRMesh::getSpaceDimension() const
101 return _mesh->getSpaceDimension();
104 void MEDCouplingCartesianAMRMesh::setFactors(const std::vector<int>& newFactors)
106 if(getSpaceDimension()!=(int)newFactors.size())
107 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::setFactors : size of input factors is not equal to the space dimension !");
113 if(_factors==newFactors)
115 if(!_patches.empty())
116 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::setFactors : modification of factors is not allowed when presence of patches !");
120 int MEDCouplingCartesianAMRMesh::getMaxNumberOfLevelsRelativeToThis() const
123 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
124 ret=std::max(ret,(*it)->getMaxNumberOfLevelsRelativeToThis()+1);
128 int MEDCouplingCartesianAMRMesh::getNumberOfCellsAtCurrentLevel() const
130 return _mesh->getNumberOfCells();
133 int MEDCouplingCartesianAMRMesh::getNumberOfCellsRecursiveWithOverlap() const
135 int ret(_mesh->getNumberOfCells());
136 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
138 ret+=(*it)->getNumberOfCellsRecursiveWithOverlap();
143 int MEDCouplingCartesianAMRMesh::getNumberOfCellsRecursiveWithoutOverlap() const
145 int ret(_mesh->getNumberOfCells());
146 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
148 ret-=(*it)->getNumberOfOverlapedCellsForFather();
149 ret+=(*it)->getNumberOfCellsRecursiveWithoutOverlap();
154 const MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::getFather() const
159 const MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::getGodFather() const
164 return _father->getGodFather();
167 void MEDCouplingCartesianAMRMesh::detachFromFather()
173 * \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,
174 * a the end cell (\b excluded) of the range for the second element of the pair.
175 * \param [in] factors The factor of refinement per axis (different from 0).
177 void MEDCouplingCartesianAMRMesh::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
179 checkFactorsAndIfNotSetAssign(factors);
180 MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> mesh(static_cast<MEDCouplingIMesh *>(_mesh->buildStructuredSubPart(bottomLeftTopRight)));
181 mesh->refineWithFactor(factors);
182 MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRMesh> zeMesh(new MEDCouplingCartesianAMRMesh(this,mesh));
183 MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> elt(new MEDCouplingCartesianAMRPatch(zeMesh,bottomLeftTopRight));
184 _patches.push_back(elt);
189 class InternalPatch : public RefCountObjectOnly
192 InternalPatch():_nb_of_true(0) { }
193 int getDimension() const { return (int)_part.size(); }
194 double getEfficiency() const { return (double)_nb_of_true/(double)_crit.size(); }
195 int getNumberOfCells() const { return (int)_crit.size(); }
196 void setNumberOfTrue(int nboft) { _nb_of_true=nboft; }
197 std::vector<bool>& getCriterion() { return _crit; }
198 const std::vector<bool>& getConstCriterion() const { return _crit; }
199 void setPart(const std::vector< std::pair<int,int> >& part) { _part=part; }
200 std::vector< std::pair<int,int> >& getPart() { return _part; }
201 const std::vector< std::pair<int,int> >& getConstPart() const { return _part; }
202 bool presenceOfTrue() const { return _nb_of_true>0; }
203 std::vector<int> computeCGS() const { return MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(_part); }
204 std::vector< std::vector<int> > computeSignature() const { return MEDCouplingStructuredMesh::ComputeSignaturePerAxisOf(computeCGS(),getConstCriterion()); }
205 double getEfficiencyPerAxis(int axisId) const { return (double)_nb_of_true/((double)(_part[axisId].second-_part[axisId].first)); }
206 void zipToFitOnCriterion();
207 void updateNumberOfTrue() const;
208 MEDCouplingAutoRefCountObjectPtr<InternalPatch> extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const;
209 MEDCouplingAutoRefCountObjectPtr<InternalPatch> deepCpy() const;
213 mutable int _nb_of_true;
214 std::vector<bool> _crit;
216 std::vector< std::pair<int,int> > _part;
219 void InternalPatch::zipToFitOnCriterion()
221 std::vector<int> cgs(computeCGS());
222 std::vector<bool> newCrit;
223 std::vector< std::pair<int,int> > newPart,newPart2;
224 int newNbOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,_crit,newCrit,newPart));
225 MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(_part,newPart,newPart2);
226 if(newNbOfTrue!=_nb_of_true)
227 throw INTERP_KERNEL::Exception("InternalPatch::zipToFitOnCrit : internal error !");
228 _crit=newCrit; _part=newPart2;
231 void InternalPatch::updateNumberOfTrue() const
233 _nb_of_true=(int)std::count(_crit.begin(),_crit.end(),true);
236 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const
238 MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
239 std::vector<int> cgs(computeCGS());
240 std::vector< std::pair<int,int> > newPart;
241 MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(_part,partInGlobal,newPart);
242 MEDCouplingStructuredMesh::ExtractFieldOfBoolFrom(cgs,_crit,newPart,ret->getCriterion());
243 ret->setPart(partInGlobal);
244 ret->updateNumberOfTrue();
248 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::deepCpy() const
250 MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
255 void DissectBigPatch(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
257 cutFound=false; cutPlace=-1;
258 std::vector<double> ratio(rangeOfAxisId-1);
259 for(int id=0;id<rangeOfAxisId-1;id++)
261 double efficiency[2];
264 std::vector< std::pair<int,int> > rectH(patchToBeSplit->getConstPart());
266 rectH[axisId].second=patchToBeSplit->getConstPart()[axisId].first+id;
268 rectH[axisId].first=patchToBeSplit->getConstPart()[axisId].first+id;
269 MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(patchToBeSplit->deepCpy());
270 p->zipToFitOnCriterion();
272 efficiency[h]=p->getEfficiencyPerAxis(axisId);
274 ratio[id]=std::max(efficiency[0],efficiency[1])/std::min(efficiency[0],efficiency[1]);
276 int minCellDirection(bso.getMinCellDirection()),indexMin(-1);
277 int dimRatioInner(rangeOfAxisId-1-2*(minCellDirection-1));
278 std::vector<double> ratio_inner(dimRatioInner);
279 double minRatio(1.e10);
280 for(int i=0; i<dimRatioInner; i++)
282 if(ratio[minCellDirection-1+i]<minRatio)
284 minRatio=ratio[minCellDirection-1+i];
285 indexMin=i+minCellDirection;
288 cutFound=true; cutPlace=indexMin+patchToBeSplit->getConstPart()[axisId].first-1;
291 void FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int& axisId, bool& cutFound, int& cutPlace)
293 cutPlace=-1; cutFound=false;
294 int minCellDirection(bso.getMinCellDirection());
295 const int dim(patchToBeSplit->getDimension());
296 std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
297 for(int id=0;id<dim;id++)
299 const std::vector<int>& signature(signatures[id]);
300 std::vector<int> hole;
301 std::vector<double> distance;
302 int len((int)signature.size());
303 for(int i=0;i<len;i++)
305 if(len>= 2*minCellDirection && i >= minCellDirection-1 && i <= len-minCellDirection-1)
309 double center(((double)len/2.));
310 for(std::size_t i=0;i<hole.size();i++)
311 distance.push_back(fabs(hole[i]+1.-center));
313 std::size_t posDistanceMin(std::distance(distance.begin(),std::min_element(distance.begin(),distance.end())));
316 cutPlace=hole[posDistanceMin]+patchToBeSplit->getConstPart()[axisId].first+1;
322 void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, bool& cutFound, int& cutPlace, int& axisId)
324 cutFound=false; cutPlace=-1;// do not set axisId before to be sure that cutFound was set to true
326 const std::vector< std::pair<int,int> >& part(patchToBeSplit->getConstPart());
327 int sign,minCellDirection(bso.getMinCellDirection());
328 const int dim(patchToBeSplit->getDimension());
330 std::vector<int> zeroCrossDims(dim,-1);
331 std::vector<int> zeroCrossVals(dim,-1);
332 std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
333 for (int id=0;id<dim;id++)
335 const std::vector<int>& signature(signatures[id]);
337 std::vector<int> derivate_second_order,gradient_absolute,signe_change,zero_cross,edge,max_cross_list ;
338 std::vector<double> distance ;
340 for (unsigned int i=1;i<signature.size()-1;i++)
341 derivate_second_order.push_back(signature[i-1]-2*signature[i]+signature[i+1]) ;
343 // Gradient absolute value
344 for ( unsigned int i=1;i<derivate_second_order.size();i++)
345 gradient_absolute.push_back(fabs(derivate_second_order[i]-derivate_second_order[i-1])) ;
346 if(derivate_second_order.empty())
348 for (unsigned int i=0;i<derivate_second_order.size()-1;i++)
350 if (derivate_second_order[i]*derivate_second_order[i+1] < 0 )
352 if (derivate_second_order[i]*derivate_second_order[i+1] > 0 )
354 if (derivate_second_order[i]*derivate_second_order[i+1] == 0 )
356 if ( sign==0 || sign==-1 )
357 if ( i >= (unsigned int)minCellDirection-2 && i <= signature.size()-minCellDirection-2 )
359 zero_cross.push_back(i) ;
360 edge.push_back(gradient_absolute[i]) ;
362 signe_change.push_back(sign) ;
364 if ( zero_cross.size() > 0 )
366 int max_cross=*max_element(edge.begin(),edge.end()) ;
367 for (unsigned int i=0;i<edge.size();i++)
368 if (edge[i]==max_cross)
369 max_cross_list.push_back(zero_cross[i]+1) ;
371 double center((signature.size()/2.0));
372 for (unsigned int i=0;i<max_cross_list.size();i++)
373 distance.push_back(fabs(max_cross_list[i]+1-center));
375 float distance_min=*min_element(distance.begin(),distance.end()) ;
376 int pos_distance_min=find(distance.begin(),distance.end(),distance_min)-distance.begin();
377 int best_place = max_cross_list[pos_distance_min] + part[id].first ;
380 zeroCrossDims[id] = best_place ;
381 zeroCrossVals[id] = max_cross ;
384 derivate_second_order.clear() ;
385 gradient_absolute.clear() ;
386 signe_change.clear() ;
389 max_cross_list.clear() ;
393 if ( zeroCrossDims[0]!=-1 || zeroCrossDims[1]!=-1 )
395 int max_cross_dims = *max_element(zeroCrossVals.begin(),zeroCrossVals.end()) ;
397 if (zeroCrossVals[0]==max_cross_dims && zeroCrossVals[1]==max_cross_dims )
399 int nl_left(part[0].second-part[0].first);
400 int nc_left(part[1].second-part[1].first);
401 if ( nl_left >= nc_left )
407 max_cross_dims=std::find(zeroCrossVals.begin(),zeroCrossVals.end(),max_cross_dims)-zeroCrossVals.begin();
409 cutPlace=zeroCrossDims[max_cross_dims];
410 axisId=max_cross_dims ;
414 void TryAction4(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
417 if(patchToBeSplit->getEfficiency()<=bso.getEffeciencySnd())
419 if(rangeOfAxisId>=2*bso.getMinCellDirection())
422 cutPlace=rangeOfAxisId/2+patchToBeSplit->getConstPart()[axisId].first-1;
427 if(patchToBeSplit->getNumberOfCells()>bso.getMaxCells())
429 DissectBigPatch(bso,patchToBeSplit,axisId,rangeOfAxisId,cutFound,cutPlace);
434 MEDCouplingAutoRefCountObjectPtr<InternalPatch> DealWithNoCut(const InternalPatch *patch)
436 MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(const_cast<InternalPatch *>(patch));
441 void DealWithCut(const InternalPatch *patchToBeSplit, int axisId, int cutPlace, std::vector<MEDCouplingAutoRefCountObjectPtr<InternalPatch> >& listOfPatches)
443 MEDCouplingAutoRefCountObjectPtr<InternalPatch> leftPart,rightPart;
444 std::vector< std::pair<int,int> > rect(patchToBeSplit->getConstPart());
445 std::vector< std::pair<int,int> > leftRect(rect),rightRect(rect);
446 leftRect[axisId].second=cutPlace+1;
447 rightRect[axisId].first=cutPlace+1;
448 leftPart=patchToBeSplit->extractPart(leftRect);
449 rightPart=patchToBeSplit->extractPart(rightRect);
450 leftPart->zipToFitOnCriterion(); rightPart->zipToFitOnCriterion();
451 listOfPatches.push_back(leftPart);
452 listOfPatches.push_back(rightPart);
458 * 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.
460 void MEDCouplingCartesianAMRMesh::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const std::vector<bool>& criterion, const std::vector<int>& factors)
462 int nbCells(getNumberOfCellsAtCurrentLevel());
463 if(nbCells!=(int)criterion.size())
464 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 !");
466 std::vector<int> cgs(_mesh->getCellGridStructure());
467 std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatches,listOfPatchesOK;
469 MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(new InternalPatch);
470 p->setNumberOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,criterion,p->getCriterion(),p->getPart()));
471 if(p->presenceOfTrue())
472 listOfPatches.push_back(p);
473 while(!listOfPatches.empty())
475 std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatchesTmp;
476 for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::iterator it=listOfPatches.begin();it!=listOfPatches.end();it++)
479 int axisId,rangeOfAxisId,cutPlace;
481 MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt((*it)->getConstPart(),axisId,rangeOfAxisId);
482 if((*it)->getEfficiency()>=bso.getEffeciency() && (*it)->getNumberOfCells()<bso.getMaxCells())
483 { listOfPatchesOK.push_back(DealWithNoCut(*it)); continue; }//action 1
484 FindHole(bso,*it,axisId,cutFound,cutPlace);//axisId overwritten here if FindHole equal to true !
486 { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 2
487 FindInflection(bso,*it,cutFound,cutPlace,axisId);//axisId overwritten here if cutFound equal to true !
489 { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 3
490 TryAction4(bso,*it,axisId,rangeOfAxisId,cutFound,cutPlace);
492 { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 4
493 listOfPatchesOK.push_back(DealWithNoCut(*it));
495 listOfPatches=listOfPatchesTmp;
497 for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::const_iterator it=listOfPatchesOK.begin();it!=listOfPatchesOK.end();it++)
498 addPatch((*it)->getConstPart(),factors);
502 * 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.
504 void MEDCouplingCartesianAMRMesh::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const DataArrayByte *criterion, const std::vector<int>& factors)
506 if(!criterion || !criterion->isAllocated())
507 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterion : the criterion DataArrayByte instance must be allocated and not NULL !");
508 std::vector<bool> crit(criterion->toVectorOfBool());//check that criterion has one component.
509 createPatchesFromCriterion(bso,crit,factors);
512 void MEDCouplingCartesianAMRMesh::removeAllPatches()
518 void MEDCouplingCartesianAMRMesh::removePatch(int patchId)
520 checkPatchId(patchId);
521 int sz((int)_patches.size()),j(0);
522 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> > patches(sz-1);
523 for(int i=0;i<sz;i++)
525 patches[j++]=_patches[i];
526 (const_cast<MEDCouplingCartesianAMRMesh *>(_patches[patchId]->getMesh()))->detachFromFather();
531 int MEDCouplingCartesianAMRMesh::getNumberOfPatches() const
533 return (int)_patches.size();
536 const MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRMesh::getPatch(int patchId) const
538 checkPatchId(patchId);
539 return _patches[patchId];
543 * 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.
544 * 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
545 * defined by the patch with id \a patchId.
547 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
548 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
549 * \return DataArrayDouble * - The array of the cell field on the requested patch
551 * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
552 * \throw if \a cellFieldOnThis is NULL or not allocated
553 * \sa fillCellFieldOnPatch, MEDCouplingIMesh::SpreadCoarseToFine
555 DataArrayDouble *MEDCouplingCartesianAMRMesh::createCellFieldOnPatch(int patchId, const DataArrayDouble *cellFieldOnThis) const
557 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
558 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
559 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
560 const MEDCouplingIMesh *fine(patch->getMesh()->getImageMesh());
561 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(fine->getNumberOfCells(),cellFieldOnThis->getNumberOfComponents());
562 ret->copyStringInfoFrom(*cellFieldOnThis);
563 MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),ret,patch->getBLTRRange(),getFactors());
568 * This method is equivalent to MEDCouplingCartesianAMRMesh::createCellFieldOnPatch except that here instead of
570 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
571 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
572 * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
574 * \sa createCellFieldOnPatch, fillCellFieldComingFromPatch
576 void MEDCouplingCartesianAMRMesh::fillCellFieldOnPatch(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch) const
578 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
579 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
580 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
581 MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors());
585 * This method updates \a cellFieldOnThis part of values coming from the cell field \a cellFieldOnPatch lying on patch having id \a patchId.
587 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
588 * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
589 * \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.
591 * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
592 * \throw if \a cellFieldOnPatch is NULL or not allocated
593 * \sa createCellFieldOnPatch, MEDCouplingIMesh::CondenseFineToCoarse
595 void MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch(int patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis) const
597 if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
598 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch : the input cell field array is NULL or not allocated !");
599 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
600 MEDCouplingIMesh::CondenseFineToCoarse(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis);
603 MEDCouplingUMesh *MEDCouplingCartesianAMRMesh::buildUnstructured() const
605 MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> part(_mesh->buildUnstructured());
606 std::vector<bool> bs(_mesh->getNumberOfCells(),false);
607 std::vector<int> cgs(_mesh->getCellGridStructure());
608 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > msSafe(_patches.size()+1);
610 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
612 MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
613 msSafe[ii+1]=(*it)->getMesh()->buildUnstructured();
615 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsOff(DataArrayInt::BuildListOfSwitchedOff(bs));
616 msSafe[0]=static_cast<MEDCouplingUMesh *>(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
617 std::vector< const MEDCouplingUMesh * > ms(msSafe.size());
618 for(std::size_t i=0;i<msSafe.size();i++)
620 return MEDCouplingUMesh::MergeUMeshes(ms);
624 * This method returns a mesh containing as cells that there is patches at the current level.
625 * The patches are seen like 'boxes' that is too say the refinement will not appear here.
627 * \return MEDCoupling1SGTUMesh * - A new object to be managed by the caller containing as cells as there are patches in \a this.
629 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMesh::buildMeshFromPatchEnvelop() const
631 std::vector<const MEDCoupling1SGTUMesh *> cells;
632 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> > cellsSafe;
633 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
635 const MEDCouplingCartesianAMRPatch *patch(*it);
638 MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> cell(patch->getMesh()->getImageMesh()->asSingleCell());
639 MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> cell1SGT(cell->build1SGTUnstructured());
640 cellsSafe.push_back(cell1SGT); cells.push_back(cell1SGT);
643 return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(cells);
646 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
647 const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):_father(0)
649 _mesh=MEDCouplingIMesh::New(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
652 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(MEDCouplingCartesianAMRMesh *father, MEDCouplingIMesh *mesh):_father(father)
655 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh(MEDCouplingIMesh *mesh) constructor : empty father !");
657 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh(MEDCouplingIMesh *mesh) constructor : The input mesh is null !");
658 mesh->checkCoherency();
659 _mesh=mesh; _mesh->incrRef();
662 void MEDCouplingCartesianAMRMesh::checkPatchId(int patchId) const
664 int sz(getNumberOfPatches());
665 if(patchId<0 || patchId>=sz)
667 std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::checkPatchId : invalid patchId (" << patchId << ") ! Must be in [0," << sz << ") !";
668 throw INTERP_KERNEL::Exception(oss.str().c_str());
672 void MEDCouplingCartesianAMRMesh::checkFactorsAndIfNotSetAssign(const std::vector<int>& factors)
674 if(getSpaceDimension()!=(int)factors.size())
675 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::checkFactorsAndIfNotSetAssign : invalid size of factors ! size must be equal to the spaceDimension !");
682 if(_factors!=factors)
683 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::checkFactorsAndIfNotSetAssign : the factors ");
687 std::size_t MEDCouplingCartesianAMRMesh::getHeapMemorySizeWithoutChildren() const
689 return sizeof(MEDCouplingCartesianAMRMesh);
692 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMesh::getDirectChildren() const
694 std::vector<const BigMemoryObject *> ret;
695 if((const MEDCouplingIMesh *)_mesh)
696 ret.push_back((const MEDCouplingIMesh *)_mesh);
697 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
699 if((const MEDCouplingCartesianAMRPatch*)*it)
700 ret.push_back((const MEDCouplingCartesianAMRPatch*)*it);
705 void MEDCouplingCartesianAMRMesh::updateTime() const
707 if((const MEDCouplingIMesh *)_mesh)
708 updateTimeWith(*_mesh);
709 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
711 const MEDCouplingCartesianAMRPatch *elt(*it);
714 const MEDCouplingCartesianAMRMesh *mesh(elt->getMesh());
716 updateTimeWith(*mesh);