1 // Copyright (C) 2007-2020 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 "MEDCouplingAMRAttribute.hxx"
23 #include "MEDCouplingFieldDouble.hxx"
24 #include "MEDCoupling1GTUMesh.hxx"
25 #include "MEDCouplingIMesh.hxx"
26 #include "MEDCouplingUMesh.hxx"
32 using namespace MEDCoupling;
36 mcIdType MEDCouplingCartesianAMRPatchGen::getNumberOfCellsRecursiveWithOverlap() const
38 return _mesh->getNumberOfCellsRecursiveWithOverlap();
41 mcIdType MEDCouplingCartesianAMRPatchGen::getNumberOfCellsRecursiveWithoutOverlap() const
43 return _mesh->getNumberOfCellsRecursiveWithoutOverlap();
46 mcIdType MEDCouplingCartesianAMRPatchGen::getMaxNumberOfLevelsRelativeToThis() const
48 return _mesh->getMaxNumberOfLevelsRelativeToThis();
51 MEDCouplingCartesianAMRPatchGen::MEDCouplingCartesianAMRPatchGen(MEDCouplingCartesianAMRMeshGen *mesh)
54 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen constructor : input mesh is NULL !");
55 _mesh=mesh; _mesh->incrRef();
58 MEDCouplingCartesianAMRPatchGen::MEDCouplingCartesianAMRPatchGen(const MEDCouplingCartesianAMRPatchGen& other, MEDCouplingCartesianAMRMeshGen *father):RefCountObject(other),_mesh(other._mesh)
60 const MEDCouplingCartesianAMRMeshGen *mesh(other._mesh);
62 _mesh=mesh->deepCopy(father);
65 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatchGen::getMeshSafe() const
67 const MEDCouplingCartesianAMRMeshGen *mesh(_mesh);
69 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen::getMeshSafe const : the mesh is NULL !");
73 MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatchGen::getMeshSafe()
75 MEDCouplingCartesianAMRMeshGen *mesh(_mesh);
77 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen::getMeshSafe : the mesh is NULL !");
81 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRPatchGen::getDirectChildrenWithNull() const
83 std::vector<const BigMemoryObject *> ret;
84 ret.push_back((const MEDCouplingCartesianAMRMeshGen *)_mesh);
89 * \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.
90 * \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,
91 * a the end cell (\b excluded) of the range for the second element of the pair.
93 MEDCouplingCartesianAMRPatch::MEDCouplingCartesianAMRPatch(MEDCouplingCartesianAMRMeshGen *mesh, const std::vector< std::pair<mcIdType,mcIdType> >& bottomLeftTopRight):MEDCouplingCartesianAMRPatchGen(mesh),_bl_tr(bottomLeftTopRight)
95 std::size_t dim(bottomLeftTopRight.size()),dimExp(_mesh->getSpaceDimension());
97 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch constructor : space dimension of father and input bottomLeft/topRight size mismatches !");
100 MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRPatch::deepCopy(MEDCouplingCartesianAMRMeshGen *father) const
102 return new MEDCouplingCartesianAMRPatch(*this,father);
105 void MEDCouplingCartesianAMRPatch::addPatch(const std::vector< std::pair<mcIdType,mcIdType> >& bottomLeftTopRight, const std::vector<mcIdType>& factors)
107 return getMeshSafe()->addPatch(bottomLeftTopRight,factors);
110 mcIdType MEDCouplingCartesianAMRPatch::getNumberOfOverlapedCellsForFather() const
112 return MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(_bl_tr);
116 * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
117 * the must be >= 0. \b WARNING this method only works if \a this and \a other share the same father (no check of this will be done !).
118 * Call isInMyNeighborhoodExt to deal with 2 patches not sharing directly the same father.
120 * \param [in] other - The other patch
121 * \param [in] ghostLev - The size of the neighborhood zone.
123 * \throw if \a this or \a other are invalid (end before start).
124 * \throw if \a ghostLev is \b not >= 0 .
125 * \throw if \a this and \a other have not the same space dimension.
127 * \sa isInMyNeighborhoodExt
129 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhood(const MEDCouplingCartesianAMRPatch *other, mcIdType ghostLev) const
132 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the size of the neighborhood must be >= 0 !");
134 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the input patch is NULL !");
135 const std::vector< std::pair<mcIdType,mcIdType> >& thisp(getBLTRRange());
136 const std::vector< std::pair<mcIdType,mcIdType> >& otherp(other->getBLTRRange());
137 return IsInMyNeighborhood(ghostLev==0?0:1,thisp,otherp);//make hypothesis that nb this->_mesh->getFather->getFactors() is >= ghostLev
141 * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
142 * the must be >= 0. This method works even if \a this and \a other does not share the same father. But the level between their common
143 * ancestor must be the same. If they don't have the same ancestor an exception will be thrown.
145 * \param [in] other - The other patch
146 * \param [in] ghostLev - The size of the neighborhood zone.
148 * \throw if \a this or \a other are invalid (end before start).
149 * \throw if \a ghostLev is \b not >= 0 .
150 * \throw if \a this and \a other have not the same space dimension.
151 * \throw if there is not common ancestor of \a this and \a other.
153 * \sa isInMyNeighborhood, isInMyNeighborhoodDiffLev
155 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt(const MEDCouplingCartesianAMRPatch *other, mcIdType ghostLev) const
158 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt : the size of the neighborhood must be >= 0 !");
160 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt : the input patch is NULL !");
162 const MEDCouplingCartesianAMRMeshGen *com(FindCommonAncestor(this,other,lev));//check that factors are OK
164 return isInMyNeighborhood(other,ghostLev);
165 std::vector<mcIdType> offset(ComputeOffsetFromTwoToOne(com,lev,this,other));
166 const std::vector< std::pair<mcIdType,mcIdType> >& thisp(getBLTRRange());
167 std::vector< std::pair<mcIdType,mcIdType> > otherp(other->getBLTRRange());
168 otherp=MEDCouplingStructuredMesh::TranslateCompactFrmt(otherp,offset);
169 return IsInMyNeighborhood(ghostLev,thisp,otherp);
173 * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
174 * the must be >= 0. This method works even if \a this and \a other does not share the same father.
175 * \a this is expected to be more refined than \a other. That is to say lev of \a this is greater than level of \a other.
177 * \param [in] other - The other patch
178 * \param [in] ghostLev - The size of the neighborhood zone.
180 * \throw if \a this or \a other are invalid (end before start).
181 * \throw if \a ghostLev is \b not >= 0 .
182 * \throw if \a this and \a other have not the same space dimension.
183 * \throw if there is not common ancestor of \a this and \a other.
185 * \sa isInMyNeighborhoodExt
187 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhoodDiffLev(const MEDCouplingCartesianAMRPatch *other, mcIdType ghostLev) const
189 std::vector< std::pair<mcIdType,mcIdType> > thispp,otherpp;
190 std::vector<mcIdType> factors;
191 ComputeZonesOfTwoRelativeToOneDiffLev(ghostLev,this,other,thispp,otherpp,factors);
192 return IsInMyNeighborhood(ghostLev>0?1:0,thispp,otherpp);//1 not ghostLev ! It is not a bug ( I hope :) ) ! Because as \a this is a refinement of \a other ghostLev is supposed to be <= factors
195 std::vector< std::pair<mcIdType,mcIdType> > MEDCouplingCartesianAMRPatch::getBLTRRangeRelativeToGF() const
197 std::vector< std::pair<mcIdType,mcIdType> > ret(_bl_tr);
198 const MEDCouplingCartesianAMRMeshGen *mesh(getMesh());
200 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::getBLTRRangeRelativeToGF : not valid !");
201 const MEDCouplingCartesianAMRMeshGen *fath(mesh->getFather());
203 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::getBLTRRangeRelativeToGF : not valid 2 !");
204 std::vector<mcIdType> factors(fath->getFactors());
205 std::size_t sz(ret.size());
206 for(std::size_t ii=0;ii<sz;ii++)
208 ret[ii].first*=factors[ii];
209 ret[ii].second*=factors[ii];
211 const MEDCouplingCartesianAMRMeshGen *oldFather(fath);
212 fath=oldFather->getFather();
215 mcIdType pos(fath->getPatchIdFromChildMesh(oldFather));
216 const MEDCouplingCartesianAMRPatch *p(fath->getPatch(pos));
217 const std::vector< std::pair<mcIdType,mcIdType> >& tmp(p->getBLTRRange());
218 const std::vector<mcIdType>& factors2(fath->getFactors());
219 std::transform(factors.begin(),factors.end(),factors2.begin(),factors.begin(),std::multiplies<mcIdType>());
220 for(std::size_t ii=0;ii<sz;ii++)
222 mcIdType delta(ret[ii].second-ret[ii].first);
223 ret[ii].first+=tmp[ii].first*factors[ii];
224 ret[ii].second=ret[ii].first+delta;
227 fath=oldFather->getFather();
232 std::vector<mcIdType> MEDCouplingCartesianAMRPatch::computeCellGridSt() const
234 const MEDCouplingCartesianAMRMeshGen *m(getMesh());
236 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::computeCellGridSt : no mesh held by this !");
237 const MEDCouplingCartesianAMRMeshGen *father(m->getFather());
239 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::computeCellGridSt : no father help by underlying mesh !");
240 const std::vector< std::pair<mcIdType,mcIdType> >& bltr(getBLTRRange());
241 const std::vector<mcIdType>& factors(father->getFactors());
242 std::vector<mcIdType> ret(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(bltr));
243 std::transform(ret.begin(),ret.end(),factors.begin(),ret.begin(),std::multiplies<mcIdType>());
247 bool MEDCouplingCartesianAMRPatch::IsInMyNeighborhood(mcIdType ghostLev, const std::vector< std::pair<mcIdType,mcIdType> >& p1, const std::vector< std::pair<mcIdType,mcIdType> >& p2)
249 std::size_t thispsize(p1.size());
250 if(thispsize!=p2.size())
251 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the dimensions must be the same !");
252 for(std::size_t i=0;i<thispsize;i++)
254 const std::pair<mcIdType,mcIdType>& thispp(p1[i]);
255 const std::pair<mcIdType,mcIdType>& otherpp(p2[i]);
256 if(thispp.second<thispp.first)
257 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : this patch is invalid !");
258 if(otherpp.second<otherpp.first)
259 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : this patch is invalid !");
260 if(otherpp.first==thispp.second+ghostLev-1)
262 if(otherpp.second+ghostLev-1==thispp.first)
264 mcIdType start(std::max(thispp.first,otherpp.first)),end(std::min(thispp.second,otherpp.second));
272 * \sa FindNeighborsOfSubPatchesOf
274 std::vector< std::vector< std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *> > > MEDCouplingCartesianAMRPatch::FindNeighborsOfSubPatchesOfSameLev(mcIdType ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2)
277 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindNeighborsOfSubPatchesOfSameLev : the input pointers must be not NULL !");
278 std::vector< std::vector< std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *> > > ret;
279 std::vector< const MEDCouplingCartesianAMRPatch *> p1Work(p1->getMesh()->getPatches()),p2Work(p2->getMesh()->getPatches());
280 while(!p1Work.empty())
282 std::vector< std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *> > retTmp;
283 std::vector<const MEDCouplingCartesianAMRPatch *> p1Work2,p2Work2;
284 for(std::vector<const MEDCouplingCartesianAMRPatch *>::const_iterator it1=p1Work.begin();it1!=p1Work.end();it1++)
286 for(std::vector<const MEDCouplingCartesianAMRPatch *>::const_iterator it2=p2Work.begin();it2!=p2Work.end();it2++)
288 if((*it1)->isInMyNeighborhoodExt(*it2,ghostLev>0?1:0))//1 not ghostLev ! It is not a bug ( I hope :) ) ! Because as \a this is a refinement of \a other ghostLev is supposed to be <= factors
289 retTmp.push_back(std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *>(*it1,*it2));
291 std::vector<const MEDCouplingCartesianAMRPatch *> tmp1((*it1)->getMesh()->getPatches());
292 p1Work2.insert(p1Work2.end(),tmp1.begin(),tmp1.end());
294 for(std::vector<const MEDCouplingCartesianAMRPatch *>::const_iterator it2=p2Work.begin();it2!=p2Work.end();it2++)
296 std::vector<const MEDCouplingCartesianAMRPatch *> tmp2((*it2)->getMesh()->getPatches());
297 p2Work2.insert(p2Work2.end(),tmp2.begin(),tmp2.end());
299 ret.push_back(retTmp);
307 * This method returns all pair of patches (pa,pb) so that pb is in the neighborhood of pa (size of neighborhood is \a ghostLev).
308 * pa is a refinement (a child) of \b p1 and pb is equal to \a p2. So the returned pair do not have the same level as it is the case for
309 * FindNeighborsOfSubPatchesOfSameLev.
311 * \sa FindNeighborsOfSubPatchesOfSameLev
313 void MEDCouplingCartesianAMRPatch::FindNeighborsOfSubPatchesOf(mcIdType ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, std::vector< std::pair<const MEDCouplingCartesianAMRPatch *, const MEDCouplingCartesianAMRPatch *> >& ret)
316 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindNeighborsOfSubPatchesOf : the input pointers must be not NULL !");
317 std::vector< const MEDCouplingCartesianAMRPatch *> p1Work(p1->getMesh()->getPatches());
318 while(!p1Work.empty())
320 std::vector<const MEDCouplingCartesianAMRPatch *> p1Work2;
321 for(std::vector<const MEDCouplingCartesianAMRPatch *>::const_iterator it0=p1Work.begin();it0!=p1Work.end();it0++)
323 if((*it0)->isInMyNeighborhoodDiffLev(p2,ghostLev))
324 ret.push_back(std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *>(*it0,p2));
325 std::vector<const MEDCouplingCartesianAMRPatch *> tmp2((*it0)->getMesh()->getPatches());
326 p1Work2.insert(p1Work2.end(),tmp2.begin(),tmp2.end());
333 * \a p1 and \a p2 are expected to be neighbors (inside the \a ghostLev zone). This method updates \a dataOnP1 only in the ghost part using a part of \a dataOnP2.
335 * \saUpdateNeighborsOfOneWithTwoExt
337 void MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwo(mcIdType ghostLev, const std::vector<mcIdType>& factors, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, DataArrayDouble *dataOnP1, const DataArrayDouble *dataOnP2)
339 const std::vector< std::pair<mcIdType,mcIdType> >& p1BLTR(p1->getBLTRRange());
340 const std::vector< std::pair<mcIdType,mcIdType> >& p2BLTR(p2->getBLTRRange());
341 UpdateNeighborsOfOneWithTwoInternal(ghostLev,factors,p1BLTR,p2BLTR,dataOnP1,dataOnP2);
345 * Idem than UpdateNeighborsOfOneWithTwo, except that here \a p1 and \a p2 are not sharing the same direct father.
347 * \sa UpdateNeighborsOfOneWithTwo
349 void MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwoExt(mcIdType ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, DataArrayDouble *dataOnP1, const DataArrayDouble *dataOnP2)
351 const std::vector< std::pair<mcIdType,mcIdType> >& p1BLTR(p1->getBLTRRange());//p1BLTR=[(10,12),(5,8)]
352 std::vector< std::pair<mcIdType,mcIdType> > p2BLTR(p2->getBLTRRange());//p2BLTR=[(0,1),(0,5)]
354 const MEDCouplingCartesianAMRMeshGen *ca(FindCommonAncestor(p1,p2,lev));
355 std::vector<mcIdType> offset(ComputeOffsetFromTwoToOne(ca,lev,p1,p2));//[12,4]
356 p2BLTR=MEDCouplingStructuredMesh::TranslateCompactFrmt(p2BLTR,offset);//p2BLTR=[(12,13),(4,9)]
357 UpdateNeighborsOfOneWithTwoInternal(ghostLev,p1->getMesh()->getFather()->getFactors(),p1BLTR,p2BLTR,dataOnP1,dataOnP2);
361 * \a p1 is expected to be more refined than \a p2. \a p1 and \a p2 have to share a common ancestor. Compared to UpdateNeighborsOfOneWithTwoExt here \a p1 and \a p2 are \b not at the same level !
363 void MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwoMixedLev(mcIdType ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, DataArrayDouble *dataOnP1, const DataArrayDouble *dataOnP2, bool isConservative)
365 std::vector< std::pair<mcIdType,mcIdType> > p1pp,p2pp;
366 std::vector<mcIdType> factors;
367 ComputeZonesOfTwoRelativeToOneDiffLev(ghostLev,p1,p2,p1pp,p2pp,factors);
369 std::vector<mcIdType> dimsP2NotRefined(p2->computeCellGridSt());
370 std::vector<mcIdType> dimsP2Refined(dimsP2NotRefined);
371 std::transform(dimsP2NotRefined.begin(),dimsP2NotRefined.end(),factors.begin(),dimsP2Refined.begin(),std::multiplies<mcIdType>());
372 std::vector< std::pair<mcIdType,mcIdType> > p2RefinedAbs(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(dimsP2NotRefined));
373 std::vector<mcIdType> dimsP2RefinedGhost(dimsP2Refined.size());
374 std::transform(dimsP2Refined.begin(),dimsP2Refined.end(),dimsP2RefinedGhost.begin(),std::bind2nd(std::plus<mcIdType>(),2*ghostLev));
375 MCAuto<DataArrayDouble> fineP2(DataArrayDouble::New()); fineP2->alloc(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(dimsP2RefinedGhost),dataOnP2->getNumberOfComponents());
376 MEDCouplingIMesh::SpreadCoarseToFineGhost(dataOnP2,dimsP2NotRefined,fineP2,p2RefinedAbs,factors,ghostLev);
379 mcIdType fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(factors));
380 std::transform(fineP2->begin(),fineP2->end(),fineP2->getPointer(),std::bind2nd(std::multiplies<double>(),1./((double)fact)));
383 UpdateNeighborsOfOneWithTwoInternal(ghostLev,p1->getMesh()->getFather()->getFactors(),p1pp,p2pp,dataOnP1,fineP2);
387 * \a p1 is expected to be more refined than \a p2. \a p1 and \a p2 have to share a common ancestor. Compared to UpdateNeighborsOfOneWithTwoExt here \a p1 and \a p2 are \b not at the same level !
388 * This method has 3 outputs. 2 two first are the resp the position of \a p1 and \a p2 relative to \a p1. And \a factToApplyOn2 is the coeff of refinement to be applied on \a p2 to be virtually
389 * on the same level as \a p1.
391 void MEDCouplingCartesianAMRPatch::ComputeZonesOfTwoRelativeToOneDiffLev(mcIdType ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, std::vector< std::pair<mcIdType,mcIdType> >& p1Zone, std::vector< std::pair<mcIdType,mcIdType> >& p2Zone, std::vector<mcIdType>& factToApplyOn2)
393 std::vector<const MEDCouplingCartesianAMRMeshGen *> ancestorsOfThis;
394 const MEDCouplingCartesianAMRMeshGen *work(p1->getMesh()),*work2(0);
395 ancestorsOfThis.push_back(work);
398 work=work->getFather();
400 ancestorsOfThis.push_back(work);
405 std::size_t levThis(0),levOther(0);
406 while(work && !found)
409 work=work->getFather();
413 std::vector<const MEDCouplingCartesianAMRMeshGen *>::iterator it(std::find(ancestorsOfThis.begin(),ancestorsOfThis.end(),work));
414 if(it!=ancestorsOfThis.end())
416 levThis=std::distance(ancestorsOfThis.begin(),it);
422 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeZonesOfTwoRelativeToOneDiffLev : no common ancestor found !");
423 if(levThis<=levOther)
424 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeZonesOfTwoRelativeToOneDiffLev : this method is not called correctly !");
426 const MEDCouplingCartesianAMRMeshGen *comAncestor(ancestorsOfThis[levThis]);
427 mcIdType idThis(comAncestor->getPatchIdFromChildMesh(ancestorsOfThis[levThis-1])),idOther(comAncestor->getPatchIdFromChildMesh(work2));
428 const MEDCouplingCartesianAMRPatch *thisp(comAncestor->getPatch(idThis)),*otherp(comAncestor->getPatch(idOther));
429 std::vector<mcIdType> offset(ComputeOffsetFromTwoToOne(comAncestor,ToIdType(levOther),thisp,otherp));
430 p1Zone=thisp->getBLTRRange(); p2Zone=MEDCouplingStructuredMesh::TranslateCompactFrmt(otherp->getBLTRRange(),offset);
431 factToApplyOn2.resize(p1Zone.size()); std::fill(factToApplyOn2.begin(),factToApplyOn2.end(),1);
433 std::size_t nbOfTurn(levThis-levOther);
434 for(std::size_t i=0;i<nbOfTurn;i++)
436 std::vector< std::pair<mcIdType,mcIdType> > tmp0;
437 MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(p1Zone,p2Zone,tmp0,false);
439 const MEDCouplingCartesianAMRMeshGen *curAncestor(ancestorsOfThis[levThis-i]);
440 ApplyFactorsOnCompactFrmt(p2Zone,curAncestor->getFactors());
441 curAncestor=ancestorsOfThis[levThis-1-i];
442 const std::vector<mcIdType>& factors(curAncestor->getFactors());
443 std::transform(factToApplyOn2.begin(),factToApplyOn2.end(),factors.begin(),factToApplyOn2.begin(),std::multiplies<mcIdType>());
444 mcIdType tmpId(curAncestor->getPatchIdFromChildMesh(ancestorsOfThis[levThis-2-i]));
445 p1Zone=curAncestor->getPatch(tmpId)->getBLTRRange();
449 std::size_t MEDCouplingCartesianAMRPatch::getHeapMemorySizeWithoutChildren() const
451 std::size_t ret(sizeof(MEDCouplingCartesianAMRPatch));
452 ret+=_bl_tr.capacity()*sizeof(std::pair<mcIdType,mcIdType>);
456 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatch::FindCommonAncestor(const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, mcIdType& lev)
458 const MEDCouplingCartesianAMRMeshGen *f1(p1->_mesh),*f2(p2->_mesh);
460 while(f1!=f2 || f1==0 || f2==0)
462 f1=f1->getFather(); f2=f2->getFather();
463 if(f1->getFactors()!=f2->getFactors())
464 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindCommonAncestor : factors differ !");
468 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindCommonAncestor : no common ancestor between p1 and p2 !");
472 std::vector<mcIdType> MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne(const MEDCouplingCartesianAMRMeshGen *comAncestor, mcIdType lev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2)
475 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne : this method is useful only for lev > 0 !");
476 mcIdType zeLev(lev-1);
477 mcIdType dim(p1->getMesh()->getSpaceDimension());
478 if(p2->getMesh()->getSpaceDimension()!=dim)
479 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne : dimension must be the same !");
480 std::vector< mcIdType > ret(dim,0);
481 for(mcIdType i=0;i<zeLev;i++)
483 const MEDCouplingCartesianAMRMeshGen *f1(p1->_mesh),*f2(p2->_mesh);
484 const MEDCouplingCartesianAMRPatch *p1h(0),*p2h(0);
485 for(mcIdType j=0;j<lev-i;j++)
487 const MEDCouplingCartesianAMRMeshGen *f1tmp(f1->getFather()),*f2tmp(f2->getFather());
488 mcIdType pid1(f1tmp->getPatchIdFromChildMesh(f1)),pid2(f2tmp->getPatchIdFromChildMesh(f2));
489 p1h=f1tmp->getPatch(pid1); p2h=f2tmp->getPatch(pid2);
492 std::vector< std::pair<mcIdType,mcIdType> > p2c(p2h->getBLTRRange());
493 for(mcIdType k=0;k<dim;k++)
495 p2c[k].first+=ret[k];
496 p2c[k].second+=ret[k];
498 for(mcIdType k=0;k<dim;k++)
500 ret[k]=p2c[k].first-p1h->getBLTRRange()[k].first;
501 ret[k]*=f1->getFactors()[k];
507 void MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwoInternal(mcIdType ghostLev, const std::vector<mcIdType>& factors, const std::vector< std::pair<mcIdType,mcIdType> >&p1 ,const std::vector< std::pair<mcIdType,mcIdType> >&p2, DataArrayDouble *dataOnP1, const DataArrayDouble *dataOnP2)
508 {//p1=[(1,4),(2,4)] p2=[(4,5),(3,4)]
509 mcIdType dim(ToIdType(factors.size()));
510 std::vector<mcIdType> dimsCoarse(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(p1));//[3,2]
511 std::transform(dimsCoarse.begin(),dimsCoarse.end(),factors.begin(),dimsCoarse.begin(),std::multiplies<mcIdType>());//[12,8]
512 std::transform(dimsCoarse.begin(),dimsCoarse.end(),dimsCoarse.begin(),std::bind2nd(std::plus<mcIdType>(),2*ghostLev));//[14,10]
513 std::vector< std::pair<mcIdType,mcIdType> > rangeCoarse(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(dimsCoarse));//[(0,14),(0,10)]
514 std::vector<mcIdType> fakeFactors(dim,1);
516 std::vector< std::pair<mcIdType,mcIdType> > tmp0,tmp1,tmp2;
517 MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(p1,p2,tmp0,false);//tmp0=[(3,4),(1,2)]
518 ApplyFactorsOnCompactFrmt(tmp0,factors);//tmp0=[(12,16),(4,8)]
519 MEDCouplingStructuredMesh::ApplyGhostOnCompactFrmt(tmp0,ghostLev);//tmp0=[(13,17),(5,9)]
520 std::vector< std::pair<mcIdType,mcIdType> > interstRange(MEDCouplingStructuredMesh::IntersectRanges(tmp0,rangeCoarse));//interstRange=[(13,14),(5,9)]
521 MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(p2,p1,tmp1,false);//tmp1=[(-3,0),(-1,1)]
522 ApplyFactorsOnCompactFrmt(tmp1,factors);//tmp1=[(-12,-4),(-4,0)]
523 MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(tmp1,interstRange,tmp2,false);//tmp2=[(1,2),(1,5)]
525 std::vector< std::pair<mcIdType,mcIdType> > dimsFine(p2);
526 ApplyFactorsOnCompactFrmt(dimsFine,factors);
527 ApplyAllGhostOnCompactFrmt(dimsFine,ghostLev);
529 MCAuto<DataArrayDouble> ghostVals(MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(dimsFine),dataOnP2,tmp2));
530 MEDCouplingIMesh::CondenseFineToCoarse(dimsCoarse,ghostVals,interstRange,fakeFactors,dataOnP1);
533 MEDCouplingCartesianAMRPatch::MEDCouplingCartesianAMRPatch(const MEDCouplingCartesianAMRPatch& other, MEDCouplingCartesianAMRMeshGen *father):MEDCouplingCartesianAMRPatchGen(other,father),_bl_tr(other._bl_tr)
538 * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in refined reference.
539 * \param [in] factors - the factors per axis.
541 void MEDCouplingCartesianAMRPatch::ApplyFactorsOnCompactFrmt(std::vector< std::pair<mcIdType,mcIdType> >& partBeforeFact, const std::vector<mcIdType>& factors)
543 std::size_t sz(factors.size());
544 if(sz!=partBeforeFact.size())
545 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ApplyFactorsOnCompactFrmt : size of input vectors must be the same !");
546 for(std::size_t i=0;i<sz;i++)
548 partBeforeFact[i].first*=factors[i];
549 partBeforeFact[i].second*=factors[i];
554 * This method is different than ApplyGhostOnCompactFrmt. The \a partBeforeFact parameter is enlarger contrary to ApplyGhostOnCompactFrmt.
556 * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in ghost reference.
557 * \param [in] ghostSize - the ghost size of zone for all axis.
559 void MEDCouplingCartesianAMRPatch::ApplyAllGhostOnCompactFrmt(std::vector< std::pair<mcIdType,mcIdType> >& partBeforeFact, mcIdType ghostSize)
562 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ApplyAllGhostOnCompactFrmt : ghost size must be >= 0 !");
563 std::size_t sz(partBeforeFact.size());
564 for(std::size_t i=0;i<sz;i++)
566 partBeforeFact[i].first-=ghostSize;
567 partBeforeFact[i].second+=ghostSize;
571 MEDCouplingCartesianAMRPatchGF::MEDCouplingCartesianAMRPatchGF(MEDCouplingCartesianAMRMesh *mesh):MEDCouplingCartesianAMRPatchGen(mesh)
575 MEDCouplingCartesianAMRPatchGF *MEDCouplingCartesianAMRPatchGF::deepCopy(MEDCouplingCartesianAMRMeshGen *father) const
577 return new MEDCouplingCartesianAMRPatchGF(*this,father);
580 std::size_t MEDCouplingCartesianAMRPatchGF::getHeapMemorySizeWithoutChildren() const
582 return sizeof(MEDCouplingCartesianAMRPatchGF);
585 MEDCouplingCartesianAMRPatchGF::MEDCouplingCartesianAMRPatchGF(const MEDCouplingCartesianAMRPatchGF& other, MEDCouplingCartesianAMRMeshGen *father):MEDCouplingCartesianAMRPatchGen(other,father)
591 int MEDCouplingCartesianAMRMeshGen::getSpaceDimension() const
593 return _mesh->getSpaceDimension();
596 void MEDCouplingCartesianAMRMeshGen::setFactors(const std::vector<mcIdType>& newFactors)
598 if(getSpaceDimension()!=ToIdType(newFactors.size()))
599 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::setFactors : size of input factors is not equal to the space dimension !");
605 if(_factors==newFactors)
607 if(!_patches.empty())
608 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::setFactors : modification of factors is not allowed when presence of patches !");
613 mcIdType MEDCouplingCartesianAMRMeshGen::getMaxNumberOfLevelsRelativeToThis() const
616 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
617 ret=std::max(ret,(*it)->getMaxNumberOfLevelsRelativeToThis()+1);
622 * This method returns the number of cells of \a this with the help of the MEDCouplingIMesh instance representing \a this.
623 * The patches in \a this are ignored here.
624 * \sa getNumberOfCellsAtCurrentLevelGhost, getNumberOfCellsRecursiveWithOverlap
626 mcIdType MEDCouplingCartesianAMRMeshGen::getNumberOfCellsAtCurrentLevel() const
628 return _mesh->getNumberOfCells();
632 * This method returns the number of cells of \a this with the help of the MEDCouplingIMesh instance representing \a this enlarged by \a ghostLev size
633 * to take into account of the ghost cells for future computation.
634 * The patches in \a this are ignored here.
636 * \sa getNumberOfCellsAtCurrentLevel
638 mcIdType MEDCouplingCartesianAMRMeshGen::getNumberOfCellsAtCurrentLevelGhost(mcIdType ghostLev) const
640 MCAuto<MEDCouplingIMesh> tmp(_mesh->buildWithGhost(ghostLev));
641 return tmp->getNumberOfCells();
645 * This method returns the number of cells including the current level but \b also \b including recursively all cells of other levels
646 * starting from this. The set of cells which size is returned here are generally overlapping each other.
648 mcIdType MEDCouplingCartesianAMRMeshGen::getNumberOfCellsRecursiveWithOverlap() const
650 mcIdType ret=_mesh->getNumberOfCells();
651 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
653 ret+=(*it)->getNumberOfCellsRecursiveWithOverlap();
659 * This method returns the max number of cells covering all the space without overlapping.
660 * It returns the number of cells of the mesh with the highest resolution.
661 * The returned value is equal to the number of cells of mesh returned by buildUnstructured.
663 * \sa buildUnstructured
665 mcIdType MEDCouplingCartesianAMRMeshGen::getNumberOfCellsRecursiveWithoutOverlap() const
667 mcIdType ret=_mesh->getNumberOfCells();
668 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
670 ret-=(*it)->getNumberOfOverlapedCellsForFather();
671 ret+=(*it)->getNumberOfCellsRecursiveWithoutOverlap();
677 * This method returns a vector of size equal to getAbsoluteLevelRelativeTo. It allows to find position an absolute position of \a this
678 * relative to \a ref (that is typically the god father).
680 * \sa getPatchAtPosition
682 std::vector<mcIdType> MEDCouplingCartesianAMRMeshGen::getPositionRelativeTo(const MEDCouplingCartesianAMRMeshGen *ref) const
685 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPositionRelativeTo : input pointer is NULL !");
686 std::vector<mcIdType> ret;
687 getPositionRelativeToInternal(ref,ret);
688 std::reverse(ret.begin(),ret.end());
693 * \sa getPositionRelativeTo, getMeshAtPosition
695 const MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRMeshGen::getPatchAtPosition(const std::vector<mcIdType>& pos) const
697 std::size_t sz(pos.size());
699 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPatchAtPosition : empty input -> no patch by definition !");
700 mcIdType patchId(pos[0]);
701 const MEDCouplingCartesianAMRPatch *elt(getPatch(patchId));
704 if(!elt || !elt->getMesh())
705 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPatchAtPosition : NULL element found during walk !");
706 std::vector<mcIdType> pos2(pos.begin()+1,pos.end());
707 return elt->getMesh()->getPatchAtPosition(pos2);
710 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMeshGen::getMeshAtPosition(const std::vector<mcIdType>& pos) const
712 std::size_t sz(pos.size());
715 mcIdType patchId(pos[0]);
716 const MEDCouplingCartesianAMRPatch *elt(getPatch(patchId));
720 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getMeshAtPosition : NULL patch !");
721 return elt->getMesh();
723 if(!elt || !elt->getMesh())
724 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPatchAtPosition : NULL element found during walk !");
725 std::vector<mcIdType> pos2(pos.begin()+1,pos.end());
726 return elt->getMesh()->getMeshAtPosition(pos2);
730 * This method returns grids relative to god father to specified level \a absoluteLev.
732 * \return std::vector<MEDCouplingCartesianAMRPatchGen *> - objects in vector are to be managed (decrRef) by the caller.
734 std::vector<MEDCouplingCartesianAMRPatchGen *> MEDCouplingCartesianAMRMeshGen::retrieveGridsAt(mcIdType absoluteLev) const
737 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::retrieveGridsAt : absolute level must be >=0 !");
738 return getGodFather()->retrieveGridsAt(absoluteLev);
742 * \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,
743 * a the end cell (\b excluded) of the range for the second element of the pair.
744 * \param [in] factors The factor of refinement per axis (different from 0).
746 void MEDCouplingCartesianAMRMeshGen::addPatch(const std::vector< std::pair<mcIdType,mcIdType> >& bottomLeftTopRight, const std::vector<mcIdType>& factors)
748 checkFactorsAndIfNotSetAssign(factors);
749 MCAuto<MEDCouplingIMesh> mesh(static_cast<MEDCouplingIMesh *>(_mesh->buildStructuredSubPart(bottomLeftTopRight)));
750 mesh->refineWithFactor(factors);
751 MCAuto<MEDCouplingCartesianAMRMeshSub> zeMesh(new MEDCouplingCartesianAMRMeshSub(this,mesh));
752 MCAuto<MEDCouplingCartesianAMRPatch> elt(new MEDCouplingCartesianAMRPatch(zeMesh,bottomLeftTopRight));
753 _patches.push_back(elt);
759 class InternalPatch : public RefCountObjectOnly
762 InternalPatch():_nb_of_true(0) { }
763 mcIdType getDimension() const { return ToIdType(_part.size()); }
764 double getEfficiency() const { return (double)_nb_of_true/(double)_crit.size(); }
765 mcIdType getNumberOfCells() const { return ToIdType(_crit.size()); }
766 void setNumberOfTrue(mcIdType nboft) { _nb_of_true=nboft; }
767 std::vector<bool>& getCriterion() { return _crit; }
768 const std::vector<bool>& getConstCriterion() const { return _crit; }
769 void setPart(const std::vector< std::pair<mcIdType,mcIdType> >& part) { _part=part; }
770 std::vector< std::pair<mcIdType,mcIdType> >& getPart() { return _part; }
771 const std::vector< std::pair<mcIdType,mcIdType> >& getConstPart() const { return _part; }
772 bool presenceOfTrue() const { return _nb_of_true>0; }
773 std::vector<mcIdType> computeCGS() const { return MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(_part); }
774 std::vector< std::vector<mcIdType> > computeSignature() const { return MEDCouplingStructuredMesh::ComputeSignaturePerAxisOf(computeCGS(),getConstCriterion()); }
775 double getEfficiencyPerAxis(mcIdType axisId) const { return (double)_nb_of_true/((double)(_part[axisId].second-_part[axisId].first)); }
776 void zipToFitOnCriterion(mcIdType minPatchLgth);
777 void updateNumberOfTrue() const;
778 MCAuto<InternalPatch> extractPart(const std::vector< std::pair<mcIdType,mcIdType> >&partInGlobal) const;
779 MCAuto<InternalPatch> deepCopy() const;
783 mutable mcIdType _nb_of_true;
784 std::vector<bool> _crit;
786 std::vector< std::pair<mcIdType,mcIdType> > _part;
789 void InternalPatch::zipToFitOnCriterion(mcIdType minPatchLgth)
791 std::vector<mcIdType> cgs(computeCGS());
792 std::vector<bool> newCrit;
793 std::vector< std::pair<mcIdType,mcIdType> > newPart,newPart2;
794 mcIdType newNbOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(minPatchLgth,cgs,_crit,newCrit,newPart));
795 MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(_part,newPart,newPart2);
796 if(newNbOfTrue!=_nb_of_true)
797 throw INTERP_KERNEL::Exception("InternalPatch::zipToFitOnCrit : internal error !");
798 _crit=newCrit; _part=newPart2;
801 void InternalPatch::updateNumberOfTrue() const
803 _nb_of_true=ToIdType(std::count(_crit.begin(),_crit.end(),true));
806 MCAuto<InternalPatch> InternalPatch::extractPart(const std::vector< std::pair<mcIdType,mcIdType> >&partInGlobal) const
808 MCAuto<InternalPatch> ret(new InternalPatch);
809 std::vector<mcIdType> cgs(computeCGS());
810 std::vector< std::pair<mcIdType,mcIdType> > newPart;
811 MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(_part,partInGlobal,newPart);
812 MEDCouplingStructuredMesh::ExtractFieldOfBoolFrom(cgs,_crit,newPart,ret->getCriterion());
813 ret->setPart(partInGlobal);
814 ret->updateNumberOfTrue();
818 MCAuto<InternalPatch> InternalPatch::deepCopy() const
820 MCAuto<InternalPatch> ret(new InternalPatch);
825 void DissectBigPatch(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, mcIdType axisId, mcIdType largestLength, mcIdType& cutPlace)
827 mcIdType minimumPatchLength(bso.getMinimumPatchLength());
828 std::vector<double> ratio(largestLength-minimumPatchLength,std::numeric_limits<double>::max());
829 mcIdType index_min = -1;
830 double minSemiEfficiencyRatio(std::numeric_limits<double>::max());
831 double efficiencyPerAxis[2];
833 for(mcIdType i=minimumPatchLength-1;i<largestLength-minimumPatchLength;i++)
835 for(mcIdType h=0;h<2;h++)
837 std::vector< std::pair<mcIdType,mcIdType> > rectH(patchToBeSplit->getConstPart());
839 rectH[axisId].second=patchToBeSplit->getConstPart()[axisId].first+i;
841 rectH[axisId].first=patchToBeSplit->getConstPart()[axisId].first+i;
842 MCAuto<InternalPatch> p(patchToBeSplit->deepCopy());
843 p->zipToFitOnCriterion(bso.getMinimumPatchLength());
844 efficiencyPerAxis[h]=p->getEfficiencyPerAxis(axisId);
846 ratio[i]=std::max(efficiencyPerAxis[0], efficiencyPerAxis[1]) / std::min(efficiencyPerAxis[0], efficiencyPerAxis[1]);
847 if(ratio[i]<minSemiEfficiencyRatio)
849 minSemiEfficiencyRatio = ratio[i];
855 throw INTERP_KERNEL::Exception("DissectBigPatch : just call to Arthur !");
857 cutPlace=index_min+patchToBeSplit->getConstPart()[axisId].first;
860 bool FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, mcIdType axisId, mcIdType& cutPlace)
863 mcIdType minimumPatchLength(bso.getMinimumPatchLength());
864 const mcIdType dim(patchToBeSplit->getDimension());
865 std::vector< std::vector<mcIdType> > signatures(patchToBeSplit->computeSignature());
866 for(mcIdType id=0;id<dim;id++)
868 const std::vector<mcIdType>& signature(signatures[id]);
869 std::vector<mcIdType> hole;
870 std::vector<double> distance;
871 mcIdType len(ToIdType(signature.size()));
872 for(mcIdType i=minimumPatchLength-1;i<len-minimumPatchLength;i++)
877 mcIdType closestHoleToMiddle(hole[0]);
878 mcIdType oldDistanceToMiddle(std::abs(hole[0]-len/2));
879 mcIdType newDistanceToMiddle(oldDistanceToMiddle);
880 for(std::size_t i=0;i<hole.size();i++)
882 newDistanceToMiddle=std::abs(hole[i]-len/2);
883 if(newDistanceToMiddle < oldDistanceToMiddle)
885 oldDistanceToMiddle = newDistanceToMiddle;
886 closestHoleToMiddle = hole[i];
889 cutPlace=closestHoleToMiddle+patchToBeSplit->getConstPart()[axisId].first;
896 bool FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, mcIdType& cutPlace, int& axisId)
898 bool cutFound(false); cutPlace=-1;// do not set axisId before to be sure that cutFound was set to true
899 const std::vector< std::pair<mcIdType,mcIdType> >& part(patchToBeSplit->getConstPart());
900 mcIdType sign,minimumPatchLength(bso.getMinimumPatchLength());
901 const mcIdType dim(patchToBeSplit->getDimension());
903 std::vector<mcIdType> zeroCrossDims(dim,-1);
904 std::vector<mcIdType> zeroCrossVals(dim,-1);
905 std::vector< std::vector<mcIdType> > signatures(patchToBeSplit->computeSignature());
906 for (mcIdType id=0;id<dim;id++)
908 const std::vector<mcIdType>& signature(signatures[id]);
910 std::vector<mcIdType> derivate_second_order,gradient_absolute,zero_cross,edge,max_cross_list ;
911 std::vector<double> distance ;
913 for(std::size_t i=1;i<signature.size()-1;i++)
914 derivate_second_order.push_back(signature[i-1]-2*signature[i]+signature[i+1]) ;
916 // Gradient absolute value
917 for(std::size_t i=1;i<derivate_second_order.size();i++)
918 gradient_absolute.push_back(std::abs(derivate_second_order[i]-derivate_second_order[i-1])) ;
919 if(derivate_second_order.empty())
921 for(std::size_t i=1;i<derivate_second_order.size()-1;i++)
923 if (derivate_second_order[i]*derivate_second_order[i+1] < 0 )
925 if (derivate_second_order[i]*derivate_second_order[i+1] > 0 )
927 if (derivate_second_order[i]*derivate_second_order[i+1] == 0 )
929 if ( sign==0 || sign==-1 )
930 if ( i >= (std::size_t)minimumPatchLength-2 && i <= signature.size()-minimumPatchLength-2 )
932 zero_cross.push_back(ToIdType(i)) ;
933 edge.push_back(gradient_absolute[i]) ;
936 if ( zero_cross.size() > 0 )
938 mcIdType max_cross=*max_element(edge.begin(),edge.end()) ;
939 for (std::size_t i=0;i<edge.size();i++)
940 if (edge[i]==max_cross)
941 max_cross_list.push_back(zero_cross[i]+1) ;
943 double center(static_cast<double>(signature.size())/2.0);
944 for (std::size_t i=0;i<max_cross_list.size();i++)
945 distance.push_back(fabs(FromIdType<double>(max_cross_list[i])+1-center));
947 double distance_min=*min_element(distance.begin(),distance.end()) ;
948 mcIdType pos_distance_min=ToIdType(find(distance.begin(),distance.end(),distance_min)-distance.begin());
949 mcIdType best_place = max_cross_list[pos_distance_min] + part[id].first ;
952 zeroCrossDims[id] = best_place ;
953 zeroCrossVals[id] = max_cross ;
956 derivate_second_order.clear() ;
957 gradient_absolute.clear() ;
960 max_cross_list.clear() ;
964 if ( zeroCrossDims[0]!=-1 || zeroCrossDims[1]!=-1 )
966 mcIdType max_cross_dims = *max_element(zeroCrossVals.begin(),zeroCrossVals.end()) ;
968 if (zeroCrossVals[0]==max_cross_dims && zeroCrossVals[1]==max_cross_dims )
970 mcIdType nl_left(part[0].second-part[0].first);
971 mcIdType nc_left(part[1].second-part[1].first);
972 if ( nl_left >= nc_left )
978 max_cross_dims=ToIdType(std::find(zeroCrossVals.begin(),zeroCrossVals.end(),max_cross_dims)-zeroCrossVals.begin());
980 cutPlace=zeroCrossDims[max_cross_dims];
981 axisId=FromIdType<int>(max_cross_dims);
986 bool TryAction4(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, mcIdType axisId, mcIdType rangeOfAxisId, mcIdType& cutPlace)
988 if(patchToBeSplit->getEfficiency()<=bso.getEfficiencyGoal())
990 if(rangeOfAxisId>=2*bso.getMinimumPatchLength())
992 cutPlace=rangeOfAxisId/2+patchToBeSplit->getConstPart()[axisId].first-1;
999 if(patchToBeSplit->getNumberOfCells()>bso.getMaximumNbOfCellsInPatch() || rangeOfAxisId>bso.getMaximumPatchLength())
1001 DissectBigPatch(bso,patchToBeSplit,axisId,rangeOfAxisId,cutPlace);
1009 MCAuto<InternalPatch> DealWithNoCut(const InternalPatch *patch)
1011 MCAuto<InternalPatch> ret(const_cast<InternalPatch *>(patch));
1016 void DealWithCut(double minPatchLgth, const InternalPatch *patchToBeSplit, int axisId, mcIdType cutPlace, std::vector<MCAuto<InternalPatch> >& listOfPatches)
1018 MCAuto<InternalPatch> leftPart,rightPart;
1019 std::vector< std::pair<mcIdType,mcIdType> > rect(patchToBeSplit->getConstPart());
1020 std::vector< std::pair<mcIdType,mcIdType> > leftRect(rect),rightRect(rect);
1021 leftRect[axisId].second=cutPlace+1;
1022 rightRect[axisId].first=cutPlace+1;
1023 leftPart=patchToBeSplit->extractPart(leftRect);
1024 rightPart=patchToBeSplit->extractPart(rightRect);
1025 leftPart->zipToFitOnCriterion(ToIdType(minPatchLgth)); rightPart->zipToFitOnCriterion(ToIdType(minPatchLgth));
1026 listOfPatches.push_back(leftPart);
1027 listOfPatches.push_back(rightPart);
1032 void MEDCouplingCartesianAMRMeshGen::removeAllPatches()
1038 void MEDCouplingCartesianAMRMeshGen::removePatch(mcIdType patchId)
1040 checkPatchId(patchId);
1041 mcIdType sz(ToIdType(_patches.size())),j(0);
1042 std::vector< MCAuto<MEDCouplingCartesianAMRPatch> > patches(sz-1);
1043 for(mcIdType i=0;i<sz;i++)
1045 patches[j++]=_patches[i];
1046 (const_cast<MEDCouplingCartesianAMRMeshGen *>(_patches[patchId]->getMesh()))->detachFromFather();
1051 mcIdType MEDCouplingCartesianAMRMeshGen::getNumberOfPatches() const
1053 return ToIdType(_patches.size());
1057 * This method is a generic algorithm to create patches in \a this (by destroying the patches if any).
1058 * This method uses \a criterion array as a field on cells on this level.
1059 * This method only create patches at level 0 relative to \a this.
1061 * This generic algorithm can be degenerated into three child ones, depending on the arguments given; in particular depending
1062 * on whether they are equal to 0 or not.
1063 * 1/ If minimumPatchLength = maximumPatchLength = maximumPatchVolume = 0, then we have the Berger-Rigoutsos algorithm.
1064 * This algorithm was developed in 1991 and seems appropriate for sequential programming.
1065 * 2/ If maximumPatchLength = 0, then we have the Livne algorithm.
1066 * This algorithm was developed in 2004 and is an improvement of the Berger-Rigoutsos algorithm.
1067 * 3/ If maximumPatchVolume = 0, the we have the lmin-lmax algorithm.
1068 * This algorithm was developed by Arthur TALPAERT in 2014 and is an improvement of the Livne algorithm. It is especially
1069 * appropriate for parallel computing, where one patch would be given to one CPU. See Arthur TALPAERT's 2014 CANUM poster
1070 * for more information.
1073 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const std::vector<bool>& criterion, const std::vector<mcIdType>& factors)
1075 mcIdType nbCells(getNumberOfCellsAtCurrentLevel());
1076 if(nbCells!=ToIdType(criterion.size()))
1077 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion : the number of tuples of criterion array must be equal to the number of cells at the current level !");
1079 std::vector<mcIdType> cgs(_mesh->getCellGridStructure());
1080 std::vector< MCAuto<InternalPatch> > listOfPatches,listOfPatchesOK;
1082 MCAuto<InternalPatch> p(new InternalPatch);
1083 p->setNumberOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(bso.getMinimumPatchLength(),cgs,criterion,p->getCriterion(),p->getPart()));
1084 if(p->presenceOfTrue())
1085 listOfPatches.push_back(p);
1086 while(!listOfPatches.empty())
1088 std::vector< MCAuto<InternalPatch> > listOfPatchesTmp;
1089 for(std::vector< MCAuto<InternalPatch> >::iterator it=listOfPatches.begin();it!=listOfPatches.end();it++)
1093 mcIdType largestLength,cutPlace;
1094 MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt((*it)->getConstPart(),axisId,largestLength);
1095 if((*it)->getEfficiency()>=bso.getEfficiencyThreshold() && ((*it)->getNumberOfCells()>bso.getMaximumNbOfCellsInPatch() || largestLength>bso.getMaximumPatchLength()))
1097 DissectBigPatch(bso,*it,axisId,largestLength,cutPlace);
1098 DealWithCut(bso.getMinimumPatchLength(),*it,axisId,cutPlace,listOfPatchesTmp);
1101 if(FindHole(bso,*it,axisId,cutPlace))//axisId overwritten here if FindHole equal to true !
1102 { DealWithCut(bso.getMinimumPatchLength(),*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 2
1103 if(FindInflection(bso,*it,cutPlace,axisId))//axisId overwritten here if cutFound equal to true !
1104 { DealWithCut(bso.getMinimumPatchLength(),*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 3
1105 if(TryAction4(bso,*it,axisId,largestLength,cutPlace))
1106 { DealWithCut(bso.getMinimumPatchLength(),*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 4
1108 listOfPatchesOK.push_back(DealWithNoCut(*it));
1110 listOfPatches=listOfPatchesTmp;
1112 for(std::vector< MCAuto<InternalPatch> >::const_iterator it=listOfPatchesOK.begin();it!=listOfPatchesOK.end();it++)
1113 addPatch((*it)->getConstPart(),factors);
1118 * 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.
1119 * This method only create patches at level 0 relative to \a this.
1121 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const DataArrayByte *criterion, const std::vector<mcIdType>& factors)
1123 if(!criterion || !criterion->isAllocated())
1124 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion : the criterion DataArrayByte instance must be allocated and not NULL !");
1125 std::vector<bool> crit(criterion->toVectorOfBool());//check that criterion has one component.
1126 createPatchesFromCriterion(bso,crit,factors);
1130 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const DataArrayDouble *criterion, const std::vector<mcIdType>& factors, double eps)
1133 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion : null criterion pointer !");
1134 std::vector<bool> inp(criterion->toVectorOfBool(eps));
1135 createPatchesFromCriterion(bso,inp,factors);
1138 mcIdType MEDCouplingCartesianAMRMeshGen::getPatchIdFromChildMesh(const MEDCouplingCartesianAMRMeshGen *mesh) const
1141 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ret++)
1143 if((*it)->getMesh()==mesh)
1146 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPatchIdFromChildMesh : no such a mesh in my direct progeny !");
1149 std::vector< const MEDCouplingCartesianAMRPatch *> MEDCouplingCartesianAMRMeshGen::getPatches() const
1151 std::size_t sz(_patches.size());
1152 std::vector< const MEDCouplingCartesianAMRPatch *> ret(sz);
1153 for(std::size_t i=0;i<sz;i++)
1158 const MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRMeshGen::getPatch(mcIdType patchId) const
1160 checkPatchId(patchId);
1161 return _patches[patchId];
1165 * This method states if patch2 (with id \a patchId2) is in the neighborhood of patch1 (with id \a patchId1).
1166 * The neighborhood size is defined by \a ghostLev in the reference of \a this ( \b not in the reference of patches !).
1168 bool MEDCouplingCartesianAMRMeshGen::isPatchInNeighborhoodOf(mcIdType patchId1, mcIdType patchId2, mcIdType ghostLev) const
1170 const MEDCouplingCartesianAMRPatch *p1(getPatch(patchId1)),*p2(getPatch(patchId2));
1171 return p1->isInMyNeighborhood(p2,ghostLev);
1175 * 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.
1176 * 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
1177 * defined by the patch with id \a patchId.
1179 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1180 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1181 * \return DataArrayDouble * - The array of the cell field on the requested patch
1183 * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
1184 * \throw if \a cellFieldOnThis is NULL or not allocated
1185 * \sa fillCellFieldOnPatch, MEDCouplingIMesh::SpreadCoarseToFine
1187 DataArrayDouble *MEDCouplingCartesianAMRMeshGen::createCellFieldOnPatch(mcIdType patchId, const DataArrayDouble *cellFieldOnThis) const
1189 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
1190 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
1191 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1192 const MEDCouplingIMesh *fine(patch->getMesh()->getImageMesh());
1193 MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(fine->getNumberOfCells(),cellFieldOnThis->getNumberOfComponents());
1194 ret->copyStringInfoFrom(*cellFieldOnThis);
1195 MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),ret,patch->getBLTRRange(),getFactors());
1200 * This method is equivalent to MEDCouplingCartesianAMRMesh::createCellFieldOnPatch except that here instead of creating a new instance
1201 * it fills the value into the \a cellFieldOnPatch data.
1203 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1204 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1205 * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
1207 * \sa createCellFieldOnPatch, fillCellFieldComingFromPatch
1209 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatch(mcIdType patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, bool isConservative) const
1211 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
1212 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
1213 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1214 MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors());
1217 mcIdType fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(getFactors()));
1218 std::transform(cellFieldOnPatch->begin(),cellFieldOnPatch->end(),cellFieldOnPatch->getPointer(),std::bind2nd(std::multiplies<double>(),1./((double)fact)));
1223 * This method is the generalization of fillCellFieldOnPatch method. This method only projects coarse to fine without considering the
1224 * potential neighbor patches covered by the ghost cells of patch with id \a patchId.
1226 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1227 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1228 * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
1229 * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
1231 * \sa fillCellFieldOnPatch, fillCellFieldOnPatchGhostAdv
1233 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchGhost(mcIdType patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, mcIdType ghostLev, bool isConservative) const
1235 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
1236 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatchGhost : the input cell field array is NULL or not allocated !");
1237 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1238 MEDCouplingIMesh::SpreadCoarseToFineGhost(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),ghostLev);
1241 mcIdType fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(getFactors()));
1242 std::transform(cellFieldOnPatch->begin(),cellFieldOnPatch->end(),cellFieldOnPatch->getPointer(),std::bind2nd(std::multiplies<double>(),1./((double)fact)));
1247 * This method is equivalent to fillCellFieldOnPatchGhost except that here \b ONLY \b the \b ghost \b zone will be updated
1248 * in \a cellFieldOnPatch.
1250 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1251 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1252 * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled \b only \b in \b the \b ghost \b zone.
1253 * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
1255 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyOnGhostZone(mcIdType patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, mcIdType ghostLev) const
1257 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
1258 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchOnlyOnGhostZone : the input cell field array is NULL or not allocated !");
1259 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1260 MEDCouplingIMesh::SpreadCoarseToFineGhostZone(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),ghostLev);
1264 * This method is a refinement of fillCellFieldOnPatchGhost. fillCellFieldOnPatchGhost is first called.
1265 * Then for all other patches than those pointed by \a patchId that overlap the ghost zone of the patch impact the ghost zone adequately.
1267 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1268 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1269 * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
1270 * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
1271 * \param [in] arrsOnPatches - \b WARNING arrsOnPatches[patchId] is \b NOT \b const. All others are const.
1273 * \sa fillCellFieldOnPatchOnlyGhostAdv
1275 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchGhostAdv(mcIdType patchId, const DataArrayDouble *cellFieldOnThis, mcIdType ghostLev, const std::vector<const DataArrayDouble *>& arrsOnPatches, bool isConservative) const
1277 mcIdType nbp(getNumberOfPatches());
1278 if(nbp!=ToIdType(arrsOnPatches.size()))
1280 std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchGhostAdv : there are " << nbp << " patches in this and " << arrsOnPatches.size() << " arrays in the last parameter !";
1281 throw INTERP_KERNEL::Exception(oss.str().c_str());
1283 DataArrayDouble *theFieldToFill(const_cast<DataArrayDouble *>(arrsOnPatches[patchId]));
1284 // first, do as usual
1285 fillCellFieldOnPatchGhost(patchId,cellFieldOnThis,theFieldToFill,ghostLev,isConservative);
1286 fillCellFieldOnPatchOnlyGhostAdv(patchId,ghostLev,arrsOnPatches);
1290 * This method updates the patch with id \a patchId considering the only the all the patches in \a this to fill ghost zone.
1291 * So \b warning, the DataArrayDouble instance \a arrsOnPatches[patchId] is non const.
1293 * \sa getPatchIdsInTheNeighborhoodOf
1295 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyGhostAdv(mcIdType patchId, mcIdType ghostLev, const std::vector<const DataArrayDouble *>& arrsOnPatches) const
1297 mcIdType nbp(getNumberOfPatches());
1298 if(nbp!=ToIdType(arrsOnPatches.size()))
1300 std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchOnlyGhostAdv : there are " << nbp << " patches in this and " << arrsOnPatches.size() << " arrays in the last parameter !";
1301 throw INTERP_KERNEL::Exception(oss.str().c_str());
1303 const MEDCouplingCartesianAMRPatch *refP(getPatch(patchId));
1304 DataArrayDouble *theFieldToFill(const_cast<DataArrayDouble *>(arrsOnPatches[patchId]));
1305 std::vector<mcIdType> ids(getPatchIdsInTheNeighborhoodOf(patchId,ghostLev));
1306 for(std::vector<mcIdType>::const_iterator it=ids.begin();it!=ids.end();it++)
1308 const MEDCouplingCartesianAMRPatch *otherP(getPatch(*it));
1309 MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwo(ghostLev,_factors,refP,otherP,theFieldToFill,arrsOnPatches[*it]);
1313 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyOnGhostZoneWith(mcIdType ghostLev, const MEDCouplingCartesianAMRPatch *patchToBeModified, const MEDCouplingCartesianAMRPatch *neighborPatch, DataArrayDouble *cellFieldOnPatch, const DataArrayDouble *cellFieldNeighbor) const
1315 MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwo(ghostLev,_factors,patchToBeModified,neighborPatch,cellFieldOnPatch,cellFieldNeighbor);
1319 * This method updates \a cellFieldOnThis part of values coming from the cell field \a cellFieldOnPatch lying on patch having id \a patchId.
1321 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1322 * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
1323 * \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.
1324 * \param [in] isConservative - true if the field needs to be conserved. false if maximum principle has to be applied.
1326 * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
1327 * \throw if \a cellFieldOnPatch is NULL or not allocated
1328 * \sa createCellFieldOnPatch, MEDCouplingIMesh::CondenseFineToCoarse,fillCellFieldComingFromPatchGhost
1330 void MEDCouplingCartesianAMRMeshGen::fillCellFieldComingFromPatch(mcIdType patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis, bool isConservative) const
1332 if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
1333 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch : the input cell field array is NULL or not allocated !");
1334 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1335 MEDCouplingIMesh::CondenseFineToCoarse(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis);
1338 mcIdType fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(getFactors()));
1339 MEDCouplingStructuredMesh::MultiplyPartOf(_mesh->getCellGridStructure(),patch->getBLTRRange(),1./((double)fact),cellFieldOnThis);
1344 * This method is the extension of MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch managing the ghost cells. If this
1345 * method is called with \a ghostLev equal to 0 it behaves exactly as MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch.
1347 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1348 * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
1349 * \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.
1350 * \param [in] ghostLev The size of ghost zone (must be >= 0 !)
1351 * \param [in] isConservative - true if the field needs to be conserved. false if maximum principle has to be applied.
1353 * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
1354 * \throw if \a cellFieldOnPatch is NULL or not allocated
1355 * \sa fillCellFieldComingFromPatch
1357 void MEDCouplingCartesianAMRMeshGen::fillCellFieldComingFromPatchGhost(mcIdType patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis, mcIdType ghostLev, bool isConservative) const
1359 if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
1360 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatchGhost : the input cell field array is NULL or not allocated !");
1361 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1362 MEDCouplingIMesh::CondenseFineToCoarseGhost(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis,ghostLev);
1365 mcIdType fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(getFactors()));
1366 MEDCouplingStructuredMesh::MultiplyPartOfByGhost(_mesh->getCellGridStructure(),patch->getBLTRRange(),ghostLev,1./((double)fact),cellFieldOnThis);
1371 * This method finds all patches (located by their ids) that are in the neighborhood of patch with id \a patchId. The neighborhood size is
1372 * defined by ghostLev.
1374 * \param [in] patchId - the id of the considered patch.
1375 * \param [in] ghostLev - the size of the neighborhood.
1376 * \return DataArrayIdType * - 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.
1378 DataArrayIdType *MEDCouplingCartesianAMRMeshGen::findPatchesInTheNeighborhoodOf(mcIdType patchId, mcIdType ghostLev) const
1380 mcIdType nbp(getNumberOfPatches());
1381 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
1382 for(mcIdType i=0;i<nbp;i++)
1385 if(isPatchInNeighborhoodOf(i,patchId,ghostLev))
1386 ret->pushBackSilent(i);
1391 MEDCouplingUMesh *MEDCouplingCartesianAMRMeshGen::buildUnstructured() const
1393 MCAuto<MEDCouplingUMesh> part(_mesh->buildUnstructured());
1394 std::vector<bool> bs(_mesh->getNumberOfCells(),false);
1395 std::vector<mcIdType> cgs(_mesh->getCellGridStructure());
1396 std::vector< MCAuto<MEDCouplingUMesh> > msSafe(_patches.size()+1);
1398 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
1400 MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
1401 msSafe[ii+1]=(*it)->getMesh()->buildUnstructured();
1403 MCAuto<DataArrayIdType> eltsOff(DataArrayIdType::BuildListOfSwitchedOff(bs));
1404 msSafe[0]=static_cast<MEDCouplingUMesh *>(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
1405 std::vector< const MEDCouplingUMesh * > ms(msSafe.size());
1406 for(std::size_t i=0;i<msSafe.size();i++)
1408 return MEDCouplingUMesh::MergeUMeshes(ms);
1412 * This method returns a mesh containing as cells that there is patches at the current level.
1413 * The patches are seen like 'boxes' that is too say the refinement will not appear here.
1415 * \return MEDCoupling1SGTUMesh * - A new object to be managed by the caller containing as cells as there are patches in \a this.
1417 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMeshGen::buildMeshFromPatchEnvelop() const
1419 std::vector<const MEDCoupling1SGTUMesh *> cells;
1420 std::vector< MCAuto<MEDCoupling1SGTUMesh> > cellsSafe;
1421 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1423 const MEDCouplingCartesianAMRPatch *patch(*it);
1426 MCAuto<MEDCouplingIMesh> cell(patch->getMesh()->getImageMesh()->asSingleCell());
1427 MCAuto<MEDCoupling1SGTUMesh> cell1SGT(cell->build1SGTUnstructured());
1428 cellsSafe.push_back(cell1SGT); cells.push_back(cell1SGT);
1431 return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(cells);
1434 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMeshGen::buildMeshOfDirectChildrenOnly() const
1436 std::vector<const MEDCoupling1SGTUMesh *> patches;
1437 std::vector< MCAuto<MEDCoupling1SGTUMesh> > patchesSafe;
1438 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1440 const MEDCouplingCartesianAMRPatch *patch(*it);
1443 MCAuto<MEDCoupling1SGTUMesh> patchMesh(patch->getMesh()->getImageMesh()->build1SGTUnstructured());
1444 patchesSafe.push_back(patchMesh); patches.push_back(patchMesh);
1447 return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(patches);
1451 * This method works same as buildUnstructured except that arrays are given in input to build a field on cell in output.
1452 * \return MEDCouplingFieldDouble * - a newly created instance the caller has reponsability to deal with.
1453 * \sa buildUnstructured
1455 MEDCouplingFieldDouble *MEDCouplingCartesianAMRMeshGen::buildCellFieldOnRecurseWithoutOverlapWithoutGhost(mcIdType ghostSz, const std::vector<const DataArrayDouble *>& recurseArrs) const
1457 if(recurseArrs.empty())
1458 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::buildCellFieldOnRecurseWithoutOverlapWithoutGhost : array is empty ! Should never happen !");
1460 std::vector<bool> bs(_mesh->getNumberOfCells(),false);
1461 std::vector<mcIdType> cgs(_mesh->getCellGridStructure());
1462 std::vector< MCAuto<MEDCouplingFieldDouble> > msSafe(_patches.size()+1);
1464 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
1466 MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
1467 std::vector<const DataArrayDouble *> tmpArrs(extractSubTreeFromGlobalFlatten((*it)->getMesh(),recurseArrs));
1468 msSafe[ii+1]=(*it)->getMesh()->buildCellFieldOnRecurseWithoutOverlapWithoutGhost(ghostSz,tmpArrs);
1470 MCAuto<DataArrayIdType> eltsOff(DataArrayIdType::BuildListOfSwitchedOff(bs));
1472 MCAuto<MEDCouplingFieldDouble> ret(MEDCouplingFieldDouble::New(ON_CELLS));
1473 MCAuto<DataArrayDouble> arr2(extractGhostFrom(ghostSz,recurseArrs[0]));
1474 arr2=arr2->selectByTupleIdSafe(eltsOff->begin(),eltsOff->end());
1475 ret->setArray(arr2);
1476 ret->setName(arr2->getName());
1477 MCAuto<MEDCouplingUMesh> part(_mesh->buildUnstructured());
1478 MCAuto<MEDCouplingMesh> mesh(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
1482 std::vector< const MEDCouplingFieldDouble * > ms(msSafe.size());
1483 for(std::size_t i=0;i<msSafe.size();i++)
1486 return MEDCouplingFieldDouble::MergeFields(ms);
1490 * This method extracts from \arr arr the part inside \a arr by cutting the \a ghostSz external part.
1491 * \arr is expected to be an array having a number of tuples equal to \c getImageMesh()->buildWithGhost(ghostSz).
1493 DataArrayDouble *MEDCouplingCartesianAMRMeshGen::extractGhostFrom(mcIdType ghostSz, const DataArrayDouble *arr) const
1495 std::vector<mcIdType> st(_mesh->getCellGridStructure());
1496 std::vector< std::pair<mcIdType,mcIdType> > p(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(st));
1497 std::transform(st.begin(),st.end(),st.begin(),std::bind2nd(std::plus<mcIdType>(),2*ghostSz));
1498 MEDCouplingStructuredMesh::ApplyGhostOnCompactFrmt(p,ghostSz);
1499 MCAuto<DataArrayDouble> ret(MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(st,arr,p));
1504 * This method returns all the patches in \a this not equal to \a patchId that are in neighborhood of patch with id \a patchId.
1506 * \sa fillCellFieldOnPatchOnlyGhostAdv
1508 std::vector<mcIdType> MEDCouplingCartesianAMRMeshGen::getPatchIdsInTheNeighborhoodOf(mcIdType patchId, mcIdType ghostLev) const
1510 std::vector<mcIdType> ret;
1511 mcIdType nbp(getNumberOfPatches());
1513 for(mcIdType i=0;i<nbp;i++)
1516 if(isPatchInNeighborhoodOf(i,patchId,ghostLev))
1523 * This method returns a dump python of \a this. It is useful for users of createPatchesFromCriterion method for debugging.
1525 * \sa dumpPatchesOf, createPatchesFromCriterion, createPatchesFromCriterionML
1527 std::string MEDCouplingCartesianAMRMeshGen::buildPythonDumpOfThis() const
1529 std::ostringstream oss;
1530 oss << "amr=MEDCouplingCartesianAMRMesh(\""<< getImageMesh()->getName() << "\"," << getSpaceDimension() << ",[";
1531 std::vector<mcIdType> ngs(getImageMesh()->getNodeGridStructure());
1532 std::vector<double> orig(getImageMesh()->getOrigin()),dxyz(getImageMesh()->getDXYZ());
1533 std::copy(ngs.begin(),ngs.end(),std::ostream_iterator<mcIdType>(oss,","));
1535 std::copy(orig.begin(),orig.end(),std::ostream_iterator<double>(oss,","));
1537 std::copy(dxyz.begin(),dxyz.end(),std::ostream_iterator<double>(oss,","));
1539 dumpPatchesOf("amr",oss);
1543 MEDCouplingCartesianAMRMeshGen::MEDCouplingCartesianAMRMeshGen(const MEDCouplingCartesianAMRMeshGen& other):RefCountObject(other),_mesh(other._mesh),_patches(other._patches),_factors(other._factors)
1545 const MEDCouplingIMesh *mesh(other._mesh);
1547 _mesh=static_cast<MEDCouplingIMesh *>(mesh->deepCopy());
1548 std::size_t sz(other._patches.size());
1549 for(std::size_t i=0;i<sz;i++)
1551 const MEDCouplingCartesianAMRPatch *patch(other._patches[i]);
1553 _patches[i]=patch->deepCopy(this);
1557 MEDCouplingCartesianAMRMeshGen::MEDCouplingCartesianAMRMeshGen(const std::string& meshName, int spaceDim, const mcIdType *nodeStrctStart, const mcIdType *nodeStrctStop,
1558 const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop)
1560 _mesh=MEDCouplingIMesh::New(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
1563 MEDCouplingCartesianAMRMeshGen::MEDCouplingCartesianAMRMeshGen(MEDCouplingIMesh *mesh)
1566 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen(MEDCouplingIMesh *mesh) constructor : The input mesh is null !");
1567 mesh->checkConsistencyLight();
1568 _mesh=mesh; _mesh->incrRef();
1571 void MEDCouplingCartesianAMRMeshGen::checkPatchId(mcIdType patchId) const
1573 mcIdType sz(getNumberOfPatches());
1574 if(patchId<0 || patchId>=sz)
1576 std::ostringstream oss; oss << "MEDCouplingCartesianAMRMeshGen::checkPatchId : invalid patchId (" << patchId << ") ! Must be in [0," << sz << ") !";
1577 throw INTERP_KERNEL::Exception(oss.str().c_str());
1581 void MEDCouplingCartesianAMRMeshGen::checkFactorsAndIfNotSetAssign(const std::vector<mcIdType>& factors)
1583 if(getSpaceDimension()!=ToIdType(factors.size()))
1584 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::checkFactorsAndIfNotSetAssign : invalid size of factors ! size must be equal to the spaceDimension !");
1585 if(_factors.empty())
1591 if(_factors!=factors)
1592 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::checkFactorsAndIfNotSetAssign : the factors ");
1596 void MEDCouplingCartesianAMRMeshGen::retrieveGridsAtInternal(mcIdType lev, std::vector< MCAuto<MEDCouplingCartesianAMRPatchGen> >& grids) const
1600 const MEDCouplingCartesianAMRMesh *thisc(dynamic_cast<const MEDCouplingCartesianAMRMesh *>(this));//tony
1601 MCAuto<MEDCouplingCartesianAMRPatchGF> elt(new MEDCouplingCartesianAMRPatchGF(const_cast<MEDCouplingCartesianAMRMesh *>(thisc)));
1602 grids.push_back(DynamicCastSafe<MEDCouplingCartesianAMRPatchGF,MEDCouplingCartesianAMRPatchGen>(elt));
1606 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1608 const MEDCouplingCartesianAMRPatch *pt(*it);
1611 MCAuto<MEDCouplingCartesianAMRPatch> tmp1(*it);
1612 grids.push_back(DynamicCastSafe<MEDCouplingCartesianAMRPatch,MEDCouplingCartesianAMRPatchGen>(tmp1));
1618 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1620 const MEDCouplingCartesianAMRPatch *pt(*it);
1622 pt->getMesh()->retrieveGridsAtInternal(lev-1,grids);
1627 mcIdType MEDCouplingCartesianAMRMeshGen::GetGhostLevelInFineRef(mcIdType ghostLev, const std::vector<mcIdType>& factors)
1630 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::GetGhostLevelInFineRef : the ghost size must be >=0 !");
1632 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::GetGhostLevelInFineRef : no factors defined !");
1633 mcIdType ghostLevInPatchRef;
1635 ghostLevInPatchRef=0;
1638 ghostLevInPatchRef=(ghostLev-1)/factors[0]+1;
1639 for(std::size_t i=0;i<factors.size();i++)
1640 ghostLevInPatchRef=std::max(ghostLevInPatchRef,(ghostLev-1)/factors[i]+1);
1642 return ghostLevInPatchRef;
1646 * This method returns a sub set of \a all. The subset is defined by the \a head in the tree defined by \a this.
1647 * Elements in \a all are expected to be sorted from god father to most refined structure.
1649 std::vector<const DataArrayDouble *> MEDCouplingCartesianAMRMeshGen::extractSubTreeFromGlobalFlatten(const MEDCouplingCartesianAMRMeshGen *head, const std::vector<const DataArrayDouble *>& all) const
1651 mcIdType maxLev(getMaxNumberOfLevelsRelativeToThis());
1652 std::vector<const DataArrayDouble *> ret;
1653 std::vector<const MEDCouplingCartesianAMRMeshGen *> meshes(1,this);
1655 for(mcIdType i=0;i<maxLev;i++)
1657 std::vector<const MEDCouplingCartesianAMRMeshGen *> meshesTmp;
1658 for(std::vector<const MEDCouplingCartesianAMRMeshGen *>::const_iterator it=meshes.begin();it!=meshes.end();it++)
1660 if((*it)==head || head->isObjectInTheProgeny(*it))
1661 ret.push_back(all[kk]);
1663 std::vector< const MEDCouplingCartesianAMRPatch *> ps((*it)->getPatches());
1664 for(std::vector< const MEDCouplingCartesianAMRPatch *>::const_iterator it0=ps.begin();it0!=ps.end();it0++)
1666 const MEDCouplingCartesianAMRMeshGen *mesh((*it0)->getMesh());
1667 meshesTmp.push_back(mesh);
1673 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::extractSubTreeFromGlobalFlatten : the size of input vector is not compatible with number of leaves in this !");
1677 void MEDCouplingCartesianAMRMeshGen::dumpPatchesOf(const std::string& varName, std::ostream& oss) const
1680 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1682 const MEDCouplingCartesianAMRPatch *patch(*it);
1685 std::ostringstream oss2; oss2 << varName << ".addPatch([";
1686 const std::vector< std::pair<mcIdType,mcIdType> >& bltr(patch->getBLTRRange());
1687 std::size_t sz(bltr.size());
1688 for(std::size_t i=0;i<sz;i++)
1690 oss2 << "(" << bltr[i].first << "," << bltr[i].second << ")";
1695 std::copy(_factors.begin(),_factors.end(),std::ostream_iterator<mcIdType>(oss2,","));
1698 std::ostringstream oss3; oss3 << varName << "[" << j++ << "]";
1699 patch->getMesh()->dumpPatchesOf(oss3.str(),oss);
1704 std::size_t MEDCouplingCartesianAMRMeshGen::getHeapMemorySizeWithoutChildren() const
1706 return sizeof(MEDCouplingCartesianAMRMeshGen);
1709 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMeshGen::getDirectChildrenWithNull() const
1711 std::vector<const BigMemoryObject *> ret;
1712 ret.push_back((const MEDCouplingIMesh *)_mesh);
1713 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1714 ret.push_back((const MEDCouplingCartesianAMRPatch*)*it);
1718 void MEDCouplingCartesianAMRMeshGen::updateTime() const
1720 if((const MEDCouplingIMesh *)_mesh)
1721 updateTimeWith(*_mesh);
1722 for(std::vector< MCAuto<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1724 const MEDCouplingCartesianAMRPatch *elt(*it);
1727 const MEDCouplingCartesianAMRMeshGen *mesh(elt->getMesh());
1729 updateTimeWith(*mesh);
1733 MEDCouplingCartesianAMRMeshSub::MEDCouplingCartesianAMRMeshSub(MEDCouplingCartesianAMRMeshGen *father, MEDCouplingIMesh *mesh):MEDCouplingCartesianAMRMeshGen(mesh),_father(father)
1736 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshSub(MEDCouplingCartesianAMRMeshGen *father, MEDCouplingIMesh *mesh) constructor : empty father !");
1739 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMeshSub::getFather() const
1744 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMeshSub::getGodFather() const
1747 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshSub::getGodFather : Impossible to find a god father because there is a hole in chain !");
1748 return _father->getGodFather();
1752 * This method returns the level of \a this. 0 for god father. 1 for children of god father ...
1754 mcIdType MEDCouplingCartesianAMRMeshSub::getAbsoluteLevel() const
1757 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshSub::getAbsoluteLevel : Impossible to find a god father because there is a hole in chain !");
1758 return _father->getAbsoluteLevel()+1;
1761 void MEDCouplingCartesianAMRMeshSub::detachFromFather()
1767 std::vector< std::pair<mcIdType,mcIdType> > MEDCouplingCartesianAMRMeshSub::positionRelativeToGodFather(std::vector<mcIdType>& st) const
1769 st=_father->getFactors();
1770 std::size_t dim(st.size());
1771 std::vector<mcIdType> prev(st);
1772 mcIdType id(_father->getPatchIdFromChildMesh(this));
1773 const MEDCouplingCartesianAMRPatch *p(_father->getPatch(id));
1774 std::vector< std::pair<mcIdType,mcIdType> > ret(p->getBLTRRange());
1775 std::vector<mcIdType> delta(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(ret)),start(dim);
1776 std::transform(delta.begin(),delta.end(),prev.begin(),delta.begin(),std::multiplies<mcIdType>());
1777 for(std::size_t i=0;i<dim;i++)
1778 start[i]=ret[i].first;
1779 std::transform(start.begin(),start.end(),prev.begin(),start.begin(),std::multiplies<mcIdType>());
1780 const MEDCouplingCartesianAMRMeshGen *it(_father);
1781 while(!dynamic_cast<const MEDCouplingCartesianAMRMesh *>(it))
1783 const MEDCouplingCartesianAMRMeshSub *itc(static_cast<const MEDCouplingCartesianAMRMeshSub *>(it));
1784 mcIdType id2(itc->_father->getPatchIdFromChildMesh(itc));
1785 const MEDCouplingCartesianAMRPatch *p2(itc->_father->getPatch(id2));
1786 const std::vector< std::pair<mcIdType,mcIdType> >& start2(p2->getBLTRRange());
1787 std::vector<mcIdType> tmp(dim);
1788 for(std::size_t i=0;i<dim;i++)
1789 tmp[i]=start2[i].first;
1791 prev=itc->_father->getFactors();
1792 std::transform(st.begin(),st.end(),prev.begin(),st.begin(),std::multiplies<mcIdType>());
1793 std::transform(st.begin(),st.end(),tmp.begin(),tmp.begin(),std::multiplies<mcIdType>());
1794 std::transform(start.begin(),start.end(),tmp.begin(),start.begin(),std::plus<mcIdType>());
1797 for(std::size_t i=0;i<dim;i++)
1799 ret[i].first=start[i];
1800 ret[i].second=start[i]+delta[i];
1805 mcIdType MEDCouplingCartesianAMRMeshSub::getAbsoluteLevelRelativeTo(const MEDCouplingCartesianAMRMeshGen *ref) const
1810 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshSub::getAbsoluteLevelRelativeTo : ref is not in the progeny of this !");
1812 return _father->getAbsoluteLevelRelativeTo(ref)+1;
1815 MEDCouplingCartesianAMRMeshSub::MEDCouplingCartesianAMRMeshSub(const MEDCouplingCartesianAMRMeshSub& other, MEDCouplingCartesianAMRMeshGen *father):MEDCouplingCartesianAMRMeshGen(other),_father(father)
1819 MEDCouplingCartesianAMRMeshSub *MEDCouplingCartesianAMRMeshSub::deepCopy(MEDCouplingCartesianAMRMeshGen *fath) const
1821 return new MEDCouplingCartesianAMRMeshSub(*this,fath);
1825 * \sa getPositionRelativeTo
1827 void MEDCouplingCartesianAMRMeshSub::getPositionRelativeToInternal(const MEDCouplingCartesianAMRMeshGen *ref, std::vector<mcIdType>& ret) const
1832 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshSub::getPositionRelativeToInternal : ref is not in the progeny of this !");
1833 mcIdType myId(_father->getPatchIdFromChildMesh(this));
1834 ret.push_back(myId);
1835 _father->getPositionRelativeToInternal(ref,ret);
1838 MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::New(const std::string& meshName, int spaceDim, const mcIdType *nodeStrctStart, const mcIdType *nodeStrctStop,
1839 const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop)
1841 return new MEDCouplingCartesianAMRMesh(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
1844 MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::New(MEDCouplingIMesh *mesh)
1846 return new MEDCouplingCartesianAMRMesh(mesh);
1849 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMesh::getFather() const
1851 //I'm god father ! No father !
1855 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMesh::getGodFather() const
1860 mcIdType MEDCouplingCartesianAMRMesh::getAbsoluteLevel() const
1865 void MEDCouplingCartesianAMRMesh::detachFromFather()
1866 {//not a bug - do nothing
1869 std::vector< std::pair<mcIdType,mcIdType> > MEDCouplingCartesianAMRMesh::positionRelativeToGodFather(std::vector<mcIdType>& st) const
1871 st=_mesh->getCellGridStructure();
1872 return MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(st);
1875 mcIdType MEDCouplingCartesianAMRMesh::getAbsoluteLevelRelativeTo(const MEDCouplingCartesianAMRMeshGen *ref) const
1879 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::getAbsoluteLevelRelativeTo : ref is not in the progeny of this !");
1882 std::vector<MEDCouplingCartesianAMRPatchGen *> MEDCouplingCartesianAMRMesh::retrieveGridsAt(mcIdType absoluteLev) const
1884 std::vector< MCAuto<MEDCouplingCartesianAMRPatchGen> > rets;
1885 retrieveGridsAtInternal(absoluteLev,rets);
1886 std::vector< MEDCouplingCartesianAMRPatchGen * > ret(rets.size());
1887 for(std::size_t i=0;i<rets.size();i++)
1889 ret[i]=rets[i].retn();
1894 MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::deepCopy(MEDCouplingCartesianAMRMeshGen *father) const
1897 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::deepCopy : specifying a not null father for a God Father object !");
1898 return new MEDCouplingCartesianAMRMesh(*this);
1902 * This method creates a multi level patches split at once.
1903 * This method calls as times as size of \a bso createPatchesFromCriterion. Size of \a bso and size of \a factors must be the same !
1904 * \b WARNING, after the call the number of levels in \a this is equal to bso.size() + 1 !
1907 * \param [in] criterion
1908 * \param [in] factors
1909 * \param [in] eps - See DataArrayDouble::toVectorOfBool for more information about the semantic of eps.
1911 * \sa createPatchesFromCriterion
1913 void MEDCouplingCartesianAMRMesh::createPatchesFromCriterionML(const std::vector<const INTERP_KERNEL::BoxSplittingOptions *>& bso, const DataArrayDouble *criterion, const std::vector< std::vector<mcIdType> >& factors, double eps)
1915 std::size_t nbOfLevs(bso.size());
1916 if(nbOfLevs!=factors.size())
1917 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterionML : size of vectors must be the same !");
1921 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterionML : pointers in 1st arg must be not NULL !");
1922 createPatchesFromCriterion(*bso[0],criterion,factors[0],eps);
1923 for(std::size_t i=1;i<nbOfLevs;i++)
1926 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterionML : presence of a NULL BoxSplittingOptions in input vector !");
1928 std::vector<MEDCouplingCartesianAMRPatchGen *> elts(retrieveGridsAt(ToIdType((i))));
1929 std::size_t sz(elts.size());
1930 std::vector< MCAuto<MEDCouplingCartesianAMRPatchGen> > elts2(sz);
1931 std::vector< MCAuto<DataArrayDouble> > elts3(sz);
1932 for(std::size_t ii=0;ii<sz;ii++)
1935 static const char TMP_STR[]="TMP";
1936 std::vector< std::pair<std::string,int> > fieldNames(1); fieldNames[0].first=TMP_STR; fieldNames[0].second=1;
1937 MCAuto<MEDCouplingAMRAttribute> att(MEDCouplingAMRAttribute::New(this,fieldNames,0));
1939 DataArrayDouble *tmpDa(const_cast<DataArrayDouble *>(att->getFieldOn(this,TMP_STR)));
1940 tmpDa->deepCopyFrom(*criterion);
1941 att->synchronizeCoarseToFine();
1942 for(std::size_t ii=0;ii<sz;ii++)
1944 const DataArrayDouble *critOnLeaf(att->getFieldOn(const_cast<MEDCouplingCartesianAMRMeshGen *>(elts[ii]->getMesh()),TMP_STR));
1945 elts3[ii]=const_cast<DataArrayDouble *>(critOnLeaf); elts3[ii]->incrRef();
1948 for(std::size_t ii=0;ii<sz;ii++)
1949 const_cast<MEDCouplingCartesianAMRMeshGen *>(elts[ii]->getMesh())->createPatchesFromCriterion(*bso[i],elts3[ii],factors[i],eps);
1953 void MEDCouplingCartesianAMRMesh::getPositionRelativeToInternal(const MEDCouplingCartesianAMRMeshGen *ref, std::vector<mcIdType>& ret) const
1958 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(const MEDCouplingCartesianAMRMesh& other):MEDCouplingCartesianAMRMeshGen(other)
1962 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(const std::string& meshName, int spaceDim, const mcIdType *nodeStrctStart, const mcIdType *nodeStrctStop,
1963 const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):MEDCouplingCartesianAMRMeshGen(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop)
1967 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(MEDCouplingIMesh *mesh):MEDCouplingCartesianAMRMeshGen(mesh)
1971 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMesh::getDirectChildrenWithNull() const
1973 std::vector<const BigMemoryObject *> ret(MEDCouplingCartesianAMRMeshGen::getDirectChildrenWithNull());
1977 MEDCouplingCartesianAMRMesh::~MEDCouplingCartesianAMRMesh()