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 "MEDCouplingAMRAttribute.hxx"
23 #include "MEDCouplingFieldDouble.hxx"
24 #include "MEDCoupling1GTUMesh.hxx"
25 #include "MEDCouplingIMesh.hxx"
26 #include "MEDCouplingUMesh.hxx"
32 using namespace ParaMEDMEM;
36 int MEDCouplingCartesianAMRPatchGen::getNumberOfCellsRecursiveWithOverlap() const
38 return _mesh->getNumberOfCellsRecursiveWithOverlap();
41 int MEDCouplingCartesianAMRPatchGen::getNumberOfCellsRecursiveWithoutOverlap() const
43 return _mesh->getNumberOfCellsRecursiveWithoutOverlap();
46 int 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 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatchGen::getMeshSafe() const
60 const MEDCouplingCartesianAMRMeshGen *mesh(_mesh);
62 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen::getMeshSafe const : the mesh is NULL !");
66 MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatchGen::getMeshSafe()
68 MEDCouplingCartesianAMRMeshGen *mesh(_mesh);
70 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen::getMeshSafe : the mesh is NULL !");
74 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRPatchGen::getDirectChildren() const
76 std::vector<const BigMemoryObject *> ret;
77 if((const MEDCouplingCartesianAMRMeshGen *)_mesh)
78 ret.push_back((const MEDCouplingCartesianAMRMeshGen *)_mesh);
83 * \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.
84 * \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,
85 * a the end cell (\b excluded) of the range for the second element of the pair.
87 MEDCouplingCartesianAMRPatch::MEDCouplingCartesianAMRPatch(MEDCouplingCartesianAMRMeshGen *mesh, const std::vector< std::pair<int,int> >& bottomLeftTopRight):MEDCouplingCartesianAMRPatchGen(mesh),_bl_tr(bottomLeftTopRight)
89 int dim((int)bottomLeftTopRight.size()),dimExp(_mesh->getSpaceDimension());
91 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch constructor : space dimension of father and input bottomLeft/topRight size mismatches !");
94 void MEDCouplingCartesianAMRPatch::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
96 return getMeshSafe()->addPatch(bottomLeftTopRight,factors);
99 int MEDCouplingCartesianAMRPatch::getNumberOfOverlapedCellsForFather() const
101 return MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(_bl_tr);
105 * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
106 * 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 !).
107 * Call isInMyNeighborhoodExt to deal with 2 patches not sharing directly the same father.
109 * \param [in] other - The other patch
110 * \param [in] ghostLev - The size of the neighborhood zone.
112 * \throw if \a this or \a other are invalid (end before start).
113 * \throw if \a ghostLev is \b not >= 0 .
114 * \throw if \a this and \a other have not the same space dimension.
116 * \sa isInMyNeighborhoodExt
118 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhood(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const
121 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the size of the neighborhood must be >= 0 !");
123 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the input patch is NULL !");
124 const std::vector< std::pair<int,int> >& thisp(getBLTRRange());
125 const std::vector< std::pair<int,int> >& otherp(other->getBLTRRange());
126 return IsInMyNeighborhood(ghostLev,thisp,otherp);
130 * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
131 * 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
132 * ancestor must be the same. If they don't have the same ancestor an exception will be thrown.
134 * \param [in] other - The other patch
135 * \param [in] ghostLev - The size of the neighborhood zone.
137 * \throw if \a this or \a other are invalid (end before start).
138 * \throw if \a ghostLev is \b not >= 0 .
139 * \throw if \a this and \a other have not the same space dimension.
140 * \throw if there is not common ancestor of \a this and \a other.
142 * \sa isInMyNeighborhood, isInMyNeighborhoodDiffLev
144 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const
147 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt : the size of the neighborhood must be >= 0 !");
149 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt : the input patch is NULL !");
151 const MEDCouplingCartesianAMRMeshGen *com(FindCommonAncestor(this,other,lev));//check that factors are OK
153 return isInMyNeighborhood(other,ghostLev);
154 std::vector<int> offset(ComputeOffsetFromTwoToOne(com,lev,this,other));
155 const std::vector< std::pair<int,int> >& thisp(getBLTRRange());
156 std::vector< std::pair<int,int> > otherp(other->getBLTRRange());
157 otherp=MEDCouplingStructuredMesh::TranslateCompactFrmt(otherp,offset);
158 return IsInMyNeighborhood(ghostLev,thisp,otherp);
162 * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
163 * the must be >= 0. This method works even if \a this and \a other does not share the same father.
164 * \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.
166 * \param [in] other - The other patch
167 * \param [in] ghostLev - The size of the neighborhood zone.
169 * \throw if \a this or \a other are invalid (end before start).
170 * \throw if \a ghostLev is \b not >= 0 .
171 * \throw if \a this and \a other have not the same space dimension.
172 * \throw if there is not common ancestor of \a this and \a other.
174 * \sa isInMyNeighborhoodExt
176 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhoodDiffLev(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const
178 std::vector< std::pair<int,int> > thispp,otherpp;
179 std::vector<int> factors;
180 ComputeZonesOfTwoRelativeToOneDiffLev(ghostLev,this,other,thispp,otherpp,factors);
181 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
184 std::vector<int> MEDCouplingCartesianAMRPatch::computeCellGridSt() const
186 const MEDCouplingCartesianAMRMeshGen *m(getMesh());
188 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::computeCellGridSt : no mesh held by this !");
189 const MEDCouplingCartesianAMRMeshGen *father(m->getFather());
191 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::computeCellGridSt : no father help by underlying mesh !");
192 const std::vector< std::pair<int,int> >& bltr(getBLTRRange());
193 const std::vector<int>& factors(father->getFactors());
194 std::vector<int> ret(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(bltr));
195 std::transform(ret.begin(),ret.end(),factors.begin(),ret.begin(),std::multiplies<int>());
199 bool MEDCouplingCartesianAMRPatch::IsInMyNeighborhood(int ghostLev, const std::vector< std::pair<int,int> >& p1, const std::vector< std::pair<int,int> >& p2)
201 std::size_t thispsize(p1.size());
202 if(thispsize!=p2.size())
203 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the dimensions must be the same !");
204 for(std::size_t i=0;i<thispsize;i++)
206 const std::pair<int,int>& thispp(p1[i]);
207 const std::pair<int,int>& otherpp(p2[i]);
208 if(thispp.second<thispp.first)
209 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : this patch is invalid !");
210 if(otherpp.second<otherpp.first)
211 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : this patch is invalid !");
212 if(otherpp.first==thispp.second+ghostLev-1)
214 if(otherpp.second+ghostLev-1==thispp.first)
216 int start(std::max(thispp.first,otherpp.first)),end(std::min(thispp.second,otherpp.second));
224 * \sa FindNeighborsOfSubPatchesOf
226 std::vector< std::vector< std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *> > > MEDCouplingCartesianAMRPatch::FindNeighborsOfSubPatchesOfSameLev(int ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2)
229 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindNeighborsOfSubPatchesOfSameLev : the input pointers must be not NULL !");
230 std::vector< std::vector< std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *> > > ret;
231 std::vector< const MEDCouplingCartesianAMRPatch *> p1Work(p1->getMesh()->getPatches()),p2Work(p2->getMesh()->getPatches());
232 while(!p1Work.empty())
234 std::vector< std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *> > retTmp;
235 std::vector<const MEDCouplingCartesianAMRPatch *> p1Work2,p2Work2;
236 for(std::vector<const MEDCouplingCartesianAMRPatch *>::const_iterator it1=p1Work.begin();it1!=p1Work.end();it1++)
238 for(std::vector<const MEDCouplingCartesianAMRPatch *>::const_iterator it2=p2Work.begin();it2!=p2Work.end();it2++)
240 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
241 retTmp.push_back(std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *>(*it1,*it2));
243 std::vector<const MEDCouplingCartesianAMRPatch *> tmp1((*it1)->getMesh()->getPatches());
244 p1Work2.insert(p1Work2.end(),tmp1.begin(),tmp1.end());
246 for(std::vector<const MEDCouplingCartesianAMRPatch *>::const_iterator it2=p2Work.begin();it2!=p2Work.end();it2++)
248 std::vector<const MEDCouplingCartesianAMRPatch *> tmp2((*it2)->getMesh()->getPatches());
249 p2Work2.insert(p2Work2.end(),tmp2.begin(),tmp2.end());
251 ret.push_back(retTmp);
259 * This method returns all pair of patches (pa,pb) so that pb is in the neighborhood of pa (size of neighborhood is \a ghostLev).
260 * 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
261 * FindNeighborsOfSubPatchesOfSameLev.
263 * \sa FindNeighborsOfSubPatchesOfSameLev
265 void MEDCouplingCartesianAMRPatch::FindNeighborsOfSubPatchesOf(int ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, std::vector< std::pair<const MEDCouplingCartesianAMRPatch *, const MEDCouplingCartesianAMRPatch *> >& ret)
268 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindNeighborsOfSubPatchesOf : the input pointers must be not NULL !");
269 std::vector< const MEDCouplingCartesianAMRPatch *> p1Work(p1->getMesh()->getPatches());
270 while(!p1Work.empty())
272 std::vector<const MEDCouplingCartesianAMRPatch *> p1Work2;
273 for(std::vector<const MEDCouplingCartesianAMRPatch *>::const_iterator it0=p1Work.begin();it0!=p1Work.end();it0++)
275 if((*it0)->isInMyNeighborhoodDiffLev(p2,ghostLev))
276 ret.push_back(std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *>(*it0,p2));
277 std::vector<const MEDCouplingCartesianAMRPatch *> tmp2((*it0)->getMesh()->getPatches());
278 p1Work2.insert(p1Work2.end(),tmp2.begin(),tmp2.end());
285 * \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.
287 * \saUpdateNeighborsOfOneWithTwoExt
289 void MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwo(int ghostLev, const std::vector<int>& factors, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, DataArrayDouble *dataOnP1, const DataArrayDouble *dataOnP2)
291 const std::vector< std::pair<int,int> >& p1BLTR(p1->getBLTRRange());
292 const std::vector< std::pair<int,int> >& p2BLTR(p2->getBLTRRange());
293 UpdateNeighborsOfOneWithTwoInternal(ghostLev,factors,p1BLTR,p2BLTR,dataOnP1,dataOnP2);
297 * Idem than UpdateNeighborsOfOneWithTwo, except that here \a p1 and \a p2 are not sharing the same direct father.
299 * \sa UpdateNeighborsOfOneWithTwo
301 void MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwoExt(int ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, DataArrayDouble *dataOnP1, const DataArrayDouble *dataOnP2)
303 const std::vector< std::pair<int,int> >& p1BLTR(p1->getBLTRRange());//p1BLTR=[(10,12),(5,8)]
304 std::vector< std::pair<int,int> > p2BLTR(p2->getBLTRRange());//p2BLTR=[(0,1),(0,5)]
306 const MEDCouplingCartesianAMRMeshGen *ca(FindCommonAncestor(p1,p2,lev));
307 std::vector<int> offset(ComputeOffsetFromTwoToOne(ca,lev,p1,p2));//[12,4]
308 p2BLTR=MEDCouplingStructuredMesh::TranslateCompactFrmt(p2BLTR,offset);//p2BLTR=[(12,13),(4,9)]
309 UpdateNeighborsOfOneWithTwoInternal(ghostLev,p1->getMesh()->getFather()->getFactors(),p1BLTR,p2BLTR,dataOnP1,dataOnP2);
313 * \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 !
315 void MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwoMixedLev(int ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, DataArrayDouble *dataOnP1, const DataArrayDouble *dataOnP2)
317 std::vector< std::pair<int,int> > p1pp,p2pp;
318 std::vector<int> factors;
319 ComputeZonesOfTwoRelativeToOneDiffLev(ghostLev,p1,p2,p1pp,p2pp,factors);
321 std::vector<int> dimsP2NotRefined(p2->computeCellGridSt());
322 std::vector<int> dimsP2Refined(dimsP2NotRefined);
323 std::transform(dimsP2NotRefined.begin(),dimsP2NotRefined.end(),factors.begin(),dimsP2Refined.begin(),std::multiplies<int>());
324 std::vector< std::pair<int,int> > p2RefinedAbs(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(dimsP2NotRefined));
325 std::vector<int> dimsP2RefinedGhost(dimsP2Refined.size());
326 std::transform(dimsP2Refined.begin(),dimsP2Refined.end(),dimsP2RefinedGhost.begin(),std::bind2nd(std::plus<int>(),2*ghostLev));
327 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> fineP2(DataArrayDouble::New()); fineP2->alloc(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(dimsP2RefinedGhost),dataOnP2->getNumberOfComponents());
328 MEDCouplingIMesh::SpreadCoarseToFineGhost(dataOnP2,dimsP2NotRefined,fineP2,p2RefinedAbs,factors,ghostLev);
330 UpdateNeighborsOfOneWithTwoInternal(ghostLev,p1->getMesh()->getFather()->getFactors(),p1pp,p2pp,dataOnP1,fineP2);
334 * \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 !
335 * 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 virtualy
336 * on the same level as \a p1.
338 void MEDCouplingCartesianAMRPatch::ComputeZonesOfTwoRelativeToOneDiffLev(int ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, std::vector< std::pair<int,int> >& p1Zone, std::vector< std::pair<int,int> >& p2Zone, std::vector<int>& factToApplyOn2)
340 std::vector<const MEDCouplingCartesianAMRMeshGen *> ancestorsOfThis;
341 const MEDCouplingCartesianAMRMeshGen *work(p1->getMesh()),*work2(0);
342 ancestorsOfThis.push_back(work);
345 work=work->getFather();
347 ancestorsOfThis.push_back(work);
352 std::size_t levThis(0),levOther(0);
353 while(work && !found)
356 work=work->getFather();
360 std::vector<const MEDCouplingCartesianAMRMeshGen *>::iterator it(std::find(ancestorsOfThis.begin(),ancestorsOfThis.end(),work));
361 if(it!=ancestorsOfThis.end())
363 levThis=std::distance(ancestorsOfThis.begin(),it);
369 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeZonesOfTwoRelativeToOneDiffLev : no common ancestor found !");
370 if(levThis<=levOther)
371 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeZonesOfTwoRelativeToOneDiffLev : this method is not called correctly !");
373 const MEDCouplingCartesianAMRMeshGen *comAncestor(ancestorsOfThis[levThis]);
374 int idThis(comAncestor->getPatchIdFromChildMesh(ancestorsOfThis[levThis-1])),idOther(comAncestor->getPatchIdFromChildMesh(work2));
375 const MEDCouplingCartesianAMRPatch *thisp(comAncestor->getPatch(idThis)),*otherp(comAncestor->getPatch(idOther));
376 std::vector<int> offset(ComputeOffsetFromTwoToOne(comAncestor,levOther,thisp,otherp));
377 p1Zone=thisp->getBLTRRange(); p2Zone=MEDCouplingStructuredMesh::TranslateCompactFrmt(otherp->getBLTRRange(),offset);
378 factToApplyOn2.resize(p1Zone.size()); std::fill(factToApplyOn2.begin(),factToApplyOn2.end(),1);
380 std::size_t nbOfTurn(levThis-levOther);
381 for(std::size_t i=0;i<nbOfTurn;i++)
383 std::vector< std::pair<int,int> > tmp0;
384 MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(p1Zone,p2Zone,tmp0,false);
386 const MEDCouplingCartesianAMRMeshGen *curAncestor(ancestorsOfThis[levThis-i]);
387 ApplyFactorsOnCompactFrmt(p2Zone,curAncestor->getFactors());
388 curAncestor=ancestorsOfThis[levThis-1-i];
389 const std::vector<int>& factors(curAncestor->getFactors());
390 std::transform(factToApplyOn2.begin(),factToApplyOn2.end(),factors.begin(),factToApplyOn2.begin(),std::multiplies<int>());
391 int tmpId(curAncestor->getPatchIdFromChildMesh(ancestorsOfThis[levThis-2-i]));
392 p1Zone=curAncestor->getPatch(tmpId)->getBLTRRange();
396 std::size_t MEDCouplingCartesianAMRPatch::getHeapMemorySizeWithoutChildren() const
398 std::size_t ret(sizeof(MEDCouplingCartesianAMRPatch));
399 ret+=_bl_tr.capacity()*sizeof(std::pair<int,int>);
403 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatch::FindCommonAncestor(const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, int& lev)
405 const MEDCouplingCartesianAMRMeshGen *f1(p1->_mesh),*f2(p2->_mesh);
407 while(f1!=f2 || f1==0 || f2==0)
409 f1=f1->getFather(); f2=f2->getFather();
410 if(f1->getFactors()!=f2->getFactors())
411 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindCommonAncestor : factors differ !");
415 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindCommonAncestor : no common ancestor between p1 and p2 !");
419 std::vector<int> MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne(const MEDCouplingCartesianAMRMeshGen *comAncestor, int lev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2)
422 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne : this method is useful only for lev > 0 !");
424 int dim(p1->getMesh()->getSpaceDimension());
425 if(p2->getMesh()->getSpaceDimension()!=dim)
426 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne : dimension must be the same !");
427 std::vector< int > ret(dim,0);
428 for(int i=0;i<zeLev;i++)
430 const MEDCouplingCartesianAMRMeshGen *f1(p1->_mesh),*f2(p2->_mesh);
431 const MEDCouplingCartesianAMRPatch *p1h(0),*p2h(0);
432 for(int j=0;j<lev-i;j++)
434 const MEDCouplingCartesianAMRMeshGen *f1tmp(f1->getFather()),*f2tmp(f2->getFather());
435 int pid1(f1tmp->getPatchIdFromChildMesh(f1)),pid2(f2tmp->getPatchIdFromChildMesh(f2));
436 p1h=f1tmp->getPatch(pid1); p2h=f2tmp->getPatch(pid2);
439 std::vector< std::pair<int,int> > p2c(p2h->getBLTRRange());
440 for(int k=0;k<dim;k++)
442 p2c[k].first+=ret[k];
443 p2c[k].second+=ret[k];
445 for(int k=0;k<dim;k++)
447 ret[k]=p2c[k].first-p1h->getBLTRRange()[k].first;
448 ret[k]*=f1->getFactors()[k];
454 void MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwoInternal(int ghostLev, const std::vector<int>& factors, const std::vector< std::pair<int,int> >&p1 ,const std::vector< std::pair<int,int> >&p2, DataArrayDouble *dataOnP1, const DataArrayDouble *dataOnP2)
455 {//p1=[(1,4),(2,4)] p2=[(4,5),(3,4)]
456 int dim((int)factors.size());
457 std::vector<int> dimsCoarse(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(p1));//[3,2]
458 std::transform(dimsCoarse.begin(),dimsCoarse.end(),factors.begin(),dimsCoarse.begin(),std::multiplies<int>());//[12,8]
459 std::transform(dimsCoarse.begin(),dimsCoarse.end(),dimsCoarse.begin(),std::bind2nd(std::plus<int>(),2*ghostLev));//[14,10]
460 std::vector< std::pair<int,int> > rangeCoarse(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(dimsCoarse));//[(0,14),(0,10)]
461 std::vector<int> fakeFactors(dim,1);
463 std::vector< std::pair<int,int> > tmp0,tmp1,tmp2;
464 MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(p1,p2,tmp0,false);//tmp0=[(3,4),(1,2)]
465 ApplyFactorsOnCompactFrmt(tmp0,factors);//tmp0=[(12,16),(4,8)]
466 MEDCouplingStructuredMesh::ApplyGhostOnCompactFrmt(tmp0,ghostLev);//tmp0=[(13,17),(5,9)]
467 std::vector< std::pair<int,int> > interstRange(MEDCouplingStructuredMesh::IntersectRanges(tmp0,rangeCoarse));//interstRange=[(13,14),(5,9)]
468 MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(p2,p1,tmp1,false);//tmp1=[(-3,0),(-1,1)]
469 ApplyFactorsOnCompactFrmt(tmp1,factors);//tmp1=[(-12,-4),(-4,0)]
470 MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(tmp1,interstRange,tmp2,false);//tmp2=[(1,2),(1,5)]
472 std::vector< std::pair<int,int> > dimsFine(p2);
473 ApplyFactorsOnCompactFrmt(dimsFine,factors);
474 ApplyAllGhostOnCompactFrmt(dimsFine,ghostLev);
476 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ghostVals(MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(dimsFine),dataOnP2,tmp2));
477 MEDCouplingIMesh::CondenseFineToCoarse(dimsCoarse,ghostVals,interstRange,fakeFactors,dataOnP1);
481 * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in refined reference.
482 * \param [in] factors - the factors per axis.
484 void MEDCouplingCartesianAMRPatch::ApplyFactorsOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, const std::vector<int>& factors)
486 std::size_t sz(factors.size());
487 if(sz!=partBeforeFact.size())
488 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ApplyFactorsOnCompactFrmt : size of input vectors must be the same !");
489 for(std::size_t i=0;i<sz;i++)
491 partBeforeFact[i].first*=factors[i];
492 partBeforeFact[i].second*=factors[i];
497 * This method is different than ApplyGhostOnCompactFrmt. The \a partBeforeFact parameter is enlarger contrary to ApplyGhostOnCompactFrmt.
499 * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in ghost reference.
500 * \param [in] ghostSize - the ghost size of zone for all axis.
502 void MEDCouplingCartesianAMRPatch::ApplyAllGhostOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, int ghostSize)
505 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ApplyAllGhostOnCompactFrmt : ghost size must be >= 0 !");
506 std::size_t sz(partBeforeFact.size());
507 for(std::size_t i=0;i<sz;i++)
509 partBeforeFact[i].first-=ghostSize;
510 partBeforeFact[i].second+=ghostSize;
516 int MEDCouplingCartesianAMRMeshGen::getSpaceDimension() const
518 return _mesh->getSpaceDimension();
521 void MEDCouplingCartesianAMRMeshGen::setFactors(const std::vector<int>& newFactors)
523 if(getSpaceDimension()!=(int)newFactors.size())
524 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::setFactors : size of input factors is not equal to the space dimension !");
530 if(_factors==newFactors)
532 if(!_patches.empty())
533 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::setFactors : modification of factors is not allowed when presence of patches !");
538 int MEDCouplingCartesianAMRMeshGen::getMaxNumberOfLevelsRelativeToThis() const
541 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
542 ret=std::max(ret,(*it)->getMaxNumberOfLevelsRelativeToThis()+1);
547 * This method returns the number of cells of \a this with the help of the MEDCouplingIMesh instance representing \a this.
548 * The patches in \a this are ignored here.
549 * \sa getNumberOfCellsAtCurrentLevelGhost, getNumberOfCellsRecursiveWithOverlap
551 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsAtCurrentLevel() const
553 return _mesh->getNumberOfCells();
557 * 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
558 * to take into account of the ghost cells for future computation.
559 * The patches in \a this are ignored here.
561 * \sa getNumberOfCellsAtCurrentLevel
563 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsAtCurrentLevelGhost(int ghostLev) const
565 MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> tmp(_mesh->buildWithGhost(ghostLev));
566 return tmp->getNumberOfCells();
570 * This method returns the number of cells including the current level but \b also \b including recursively all cells of other levels
571 * starting from this. The set of cells which size is returned here are generally overlapping each other.
573 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsRecursiveWithOverlap() const
575 int ret(_mesh->getNumberOfCells());
576 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
578 ret+=(*it)->getNumberOfCellsRecursiveWithOverlap();
584 * This method returns the max number of cells covering all the space without overlapping.
585 * It returns the number of cells of the mesh with the highest resolution.
586 * The returned value is equal to the number of cells of mesh returned by buildUnstructured.
588 * \sa buildUnstructured
590 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsRecursiveWithoutOverlap() const
592 int ret(_mesh->getNumberOfCells());
593 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
595 ret-=(*it)->getNumberOfOverlapedCellsForFather();
596 ret+=(*it)->getNumberOfCellsRecursiveWithoutOverlap();
601 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMeshGen::getFather() const
606 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMeshGen::getGodFather() const
611 return _father->getGodFather();
615 * This method returns the level of \a this. 0 for god father. -1 for children of god father ...
617 int MEDCouplingCartesianAMRMeshGen::getAbsoluteLevel() const
622 return _father->getAbsoluteLevel()-1;
626 * This method returns grids relative to god father to specified level \a absoluteLev.
628 * \return std::vector<MEDCouplingCartesianAMRPatchGen *> - objects in vector are to be managed (decrRef) by the caller.
630 std::vector<MEDCouplingCartesianAMRPatchGen *> MEDCouplingCartesianAMRMeshGen::retrieveGridsAt(int absoluteLev) const
633 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::retrieveGridsAt : absolute level must be >=0 !");
635 return getGodFather()->retrieveGridsAt(absoluteLev);
637 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGen> > rets;
638 retrieveGridsAtInternal(absoluteLev,rets);
639 std::vector< MEDCouplingCartesianAMRPatchGen * > ret(rets.size());
640 for(std::size_t i=0;i<rets.size();i++)
642 ret[i]=rets[i].retn();
647 void MEDCouplingCartesianAMRMeshGen::detachFromFather()
654 * \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,
655 * a the end cell (\b excluded) of the range for the second element of the pair.
656 * \param [in] factors The factor of refinement per axis (different from 0).
658 void MEDCouplingCartesianAMRMeshGen::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
660 checkFactorsAndIfNotSetAssign(factors);
661 MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> mesh(static_cast<MEDCouplingIMesh *>(_mesh->buildStructuredSubPart(bottomLeftTopRight)));
662 mesh->refineWithFactor(factors);
663 MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRMeshSub> zeMesh(new MEDCouplingCartesianAMRMeshSub(this,mesh));
664 MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> elt(new MEDCouplingCartesianAMRPatch(zeMesh,bottomLeftTopRight));
665 _patches.push_back(elt);
671 class InternalPatch : public RefCountObjectOnly
674 InternalPatch():_nb_of_true(0) { }
675 int getDimension() const { return (int)_part.size(); }
676 double getEfficiency() const { return (double)_nb_of_true/(double)_crit.size(); }
677 int getNumberOfCells() const { return (int)_crit.size(); }
678 void setNumberOfTrue(int nboft) { _nb_of_true=nboft; }
679 std::vector<bool>& getCriterion() { return _crit; }
680 const std::vector<bool>& getConstCriterion() const { return _crit; }
681 void setPart(const std::vector< std::pair<int,int> >& part) { _part=part; }
682 std::vector< std::pair<int,int> >& getPart() { return _part; }
683 const std::vector< std::pair<int,int> >& getConstPart() const { return _part; }
684 bool presenceOfTrue() const { return _nb_of_true>0; }
685 std::vector<int> computeCGS() const { return MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(_part); }
686 std::vector< std::vector<int> > computeSignature() const { return MEDCouplingStructuredMesh::ComputeSignaturePerAxisOf(computeCGS(),getConstCriterion()); }
687 double getEfficiencyPerAxis(int axisId) const { return (double)_nb_of_true/((double)(_part[axisId].second-_part[axisId].first)); }
688 void zipToFitOnCriterion();
689 void updateNumberOfTrue() const;
690 MEDCouplingAutoRefCountObjectPtr<InternalPatch> extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const;
691 MEDCouplingAutoRefCountObjectPtr<InternalPatch> deepCpy() const;
695 mutable int _nb_of_true;
696 std::vector<bool> _crit;
698 std::vector< std::pair<int,int> > _part;
701 void InternalPatch::zipToFitOnCriterion()
703 std::vector<int> cgs(computeCGS());
704 std::vector<bool> newCrit;
705 std::vector< std::pair<int,int> > newPart,newPart2;
706 int newNbOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,_crit,newCrit,newPart));
707 MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(_part,newPart,newPart2);
708 if(newNbOfTrue!=_nb_of_true)
709 throw INTERP_KERNEL::Exception("InternalPatch::zipToFitOnCrit : internal error !");
710 _crit=newCrit; _part=newPart2;
713 void InternalPatch::updateNumberOfTrue() const
715 _nb_of_true=(int)std::count(_crit.begin(),_crit.end(),true);
718 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const
720 MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
721 std::vector<int> cgs(computeCGS());
722 std::vector< std::pair<int,int> > newPart;
723 MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(_part,partInGlobal,newPart);
724 MEDCouplingStructuredMesh::ExtractFieldOfBoolFrom(cgs,_crit,newPart,ret->getCriterion());
725 ret->setPart(partInGlobal);
726 ret->updateNumberOfTrue();
730 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::deepCpy() const
732 MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
737 void DissectBigPatch(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
739 cutFound=false; cutPlace=-1;
740 std::vector<double> ratio(rangeOfAxisId-1);
741 for(int id=0;id<rangeOfAxisId-1;id++)
743 double efficiency[2];
746 std::vector< std::pair<int,int> > rectH(patchToBeSplit->getConstPart());
748 rectH[axisId].second=patchToBeSplit->getConstPart()[axisId].first+id;
750 rectH[axisId].first=patchToBeSplit->getConstPart()[axisId].first+id;
751 MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(patchToBeSplit->deepCpy());
752 p->zipToFitOnCriterion();
754 efficiency[h]=p->getEfficiencyPerAxis(axisId);
756 ratio[id]=std::max(efficiency[0],efficiency[1])/std::min(efficiency[0],efficiency[1]);
758 int minCellDirection(bso.getMinCellDirection()),indexMin(-1);
759 int dimRatioInner(rangeOfAxisId-1-2*(minCellDirection-1));
760 std::vector<double> ratio_inner(dimRatioInner);
761 double minRatio(1.e10);
762 for(int i=0; i<dimRatioInner; i++)
764 if(ratio[minCellDirection-1+i]<minRatio)
766 minRatio=ratio[minCellDirection-1+i];
767 indexMin=i+minCellDirection;
770 cutFound=true; cutPlace=indexMin+patchToBeSplit->getConstPart()[axisId].first-1;
773 void FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int& axisId, bool& cutFound, int& cutPlace)
775 cutPlace=-1; cutFound=false;
776 int minCellDirection(bso.getMinCellDirection());
777 const int dim(patchToBeSplit->getDimension());
778 std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
779 for(int id=0;id<dim;id++)
781 const std::vector<int>& signature(signatures[id]);
782 std::vector<int> hole;
783 std::vector<double> distance;
784 int len((int)signature.size());
785 for(int i=0;i<len;i++)
787 if(len>= 2*minCellDirection && i >= minCellDirection-1 && i <= len-minCellDirection-1)
791 double center(((double)len/2.));
792 for(std::size_t i=0;i<hole.size();i++)
793 distance.push_back(fabs(hole[i]+1.-center));
795 std::size_t posDistanceMin(std::distance(distance.begin(),std::min_element(distance.begin(),distance.end())));
798 cutPlace=hole[posDistanceMin]+patchToBeSplit->getConstPart()[axisId].first+1;
804 void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, bool& cutFound, int& cutPlace, int& axisId)
806 cutFound=false; cutPlace=-1;// do not set axisId before to be sure that cutFound was set to true
808 const std::vector< std::pair<int,int> >& part(patchToBeSplit->getConstPart());
809 int sign,minCellDirection(bso.getMinCellDirection());
810 const int dim(patchToBeSplit->getDimension());
812 std::vector<int> zeroCrossDims(dim,-1);
813 std::vector<int> zeroCrossVals(dim,-1);
814 std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
815 for (int id=0;id<dim;id++)
817 const std::vector<int>& signature(signatures[id]);
819 std::vector<int> derivate_second_order,gradient_absolute,signe_change,zero_cross,edge,max_cross_list ;
820 std::vector<double> distance ;
822 for (unsigned int i=1;i<signature.size()-1;i++)
823 derivate_second_order.push_back(signature[i-1]-2*signature[i]+signature[i+1]) ;
825 // Gradient absolute value
826 for ( unsigned int i=1;i<derivate_second_order.size();i++)
827 gradient_absolute.push_back(fabs(derivate_second_order[i]-derivate_second_order[i-1])) ;
828 if(derivate_second_order.empty())
830 for (unsigned int i=0;i<derivate_second_order.size()-1;i++)
832 if (derivate_second_order[i]*derivate_second_order[i+1] < 0 )
834 if (derivate_second_order[i]*derivate_second_order[i+1] > 0 )
836 if (derivate_second_order[i]*derivate_second_order[i+1] == 0 )
838 if ( sign==0 || sign==-1 )
839 if ( i >= (unsigned int)minCellDirection-2 && i <= signature.size()-minCellDirection-2 )
841 zero_cross.push_back(i) ;
842 edge.push_back(gradient_absolute[i]) ;
844 signe_change.push_back(sign) ;
846 if ( zero_cross.size() > 0 )
848 int max_cross=*max_element(edge.begin(),edge.end()) ;
849 for (unsigned int i=0;i<edge.size();i++)
850 if (edge[i]==max_cross)
851 max_cross_list.push_back(zero_cross[i]+1) ;
853 double center((signature.size()/2.0));
854 for (unsigned int i=0;i<max_cross_list.size();i++)
855 distance.push_back(fabs(max_cross_list[i]+1-center));
857 float distance_min=*min_element(distance.begin(),distance.end()) ;
858 int pos_distance_min=find(distance.begin(),distance.end(),distance_min)-distance.begin();
859 int best_place = max_cross_list[pos_distance_min] + part[id].first ;
862 zeroCrossDims[id] = best_place ;
863 zeroCrossVals[id] = max_cross ;
866 derivate_second_order.clear() ;
867 gradient_absolute.clear() ;
868 signe_change.clear() ;
871 max_cross_list.clear() ;
875 if ( zeroCrossDims[0]!=-1 || zeroCrossDims[1]!=-1 )
877 int max_cross_dims = *max_element(zeroCrossVals.begin(),zeroCrossVals.end()) ;
879 if (zeroCrossVals[0]==max_cross_dims && zeroCrossVals[1]==max_cross_dims )
881 int nl_left(part[0].second-part[0].first);
882 int nc_left(part[1].second-part[1].first);
883 if ( nl_left >= nc_left )
889 max_cross_dims=std::find(zeroCrossVals.begin(),zeroCrossVals.end(),max_cross_dims)-zeroCrossVals.begin();
891 cutPlace=zeroCrossDims[max_cross_dims];
892 axisId=max_cross_dims ;
896 void TryAction4(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
899 if(patchToBeSplit->getEfficiency()<=bso.getEffeciencySnd())
901 if(rangeOfAxisId>=2*bso.getMinCellDirection())
904 cutPlace=rangeOfAxisId/2+patchToBeSplit->getConstPart()[axisId].first-1;
909 if(patchToBeSplit->getNumberOfCells()>bso.getMaxCells())
911 DissectBigPatch(bso,patchToBeSplit,axisId,rangeOfAxisId,cutFound,cutPlace);
916 MEDCouplingAutoRefCountObjectPtr<InternalPatch> DealWithNoCut(const InternalPatch *patch)
918 MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(const_cast<InternalPatch *>(patch));
923 void DealWithCut(const InternalPatch *patchToBeSplit, int axisId, int cutPlace, std::vector<MEDCouplingAutoRefCountObjectPtr<InternalPatch> >& listOfPatches)
925 MEDCouplingAutoRefCountObjectPtr<InternalPatch> leftPart,rightPart;
926 std::vector< std::pair<int,int> > rect(patchToBeSplit->getConstPart());
927 std::vector< std::pair<int,int> > leftRect(rect),rightRect(rect);
928 leftRect[axisId].second=cutPlace+1;
929 rightRect[axisId].first=cutPlace+1;
930 leftPart=patchToBeSplit->extractPart(leftRect);
931 rightPart=patchToBeSplit->extractPart(rightRect);
932 leftPart->zipToFitOnCriterion(); rightPart->zipToFitOnCriterion();
933 listOfPatches.push_back(leftPart);
934 listOfPatches.push_back(rightPart);
940 * 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.
941 * This method only create patches at level 0 relative to \a this.
943 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const std::vector<bool>& criterion, const std::vector<int>& factors)
945 int nbCells(getNumberOfCellsAtCurrentLevel());
946 if(nbCells!=(int)criterion.size())
947 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 !");
949 std::vector<int> cgs(_mesh->getCellGridStructure());
950 std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatches,listOfPatchesOK;
952 MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(new InternalPatch);
953 p->setNumberOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,criterion,p->getCriterion(),p->getPart()));
954 if(p->presenceOfTrue())
955 listOfPatches.push_back(p);
956 while(!listOfPatches.empty())
958 std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatchesTmp;
959 for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::iterator it=listOfPatches.begin();it!=listOfPatches.end();it++)
962 int axisId,rangeOfAxisId,cutPlace;
964 MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt((*it)->getConstPart(),axisId,rangeOfAxisId);
965 if((*it)->getEfficiency()>=bso.getEffeciency() && (*it)->getNumberOfCells()<bso.getMaxCells())
966 { listOfPatchesOK.push_back(DealWithNoCut(*it)); continue; }//action 1
967 FindHole(bso,*it,axisId,cutFound,cutPlace);//axisId overwritten here if FindHole equal to true !
969 { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 2
970 FindInflection(bso,*it,cutFound,cutPlace,axisId);//axisId overwritten here if cutFound equal to true !
972 { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 3
973 TryAction4(bso,*it,axisId,rangeOfAxisId,cutFound,cutPlace);
975 { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 4
976 listOfPatchesOK.push_back(DealWithNoCut(*it));
978 listOfPatches=listOfPatchesTmp;
980 for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::const_iterator it=listOfPatchesOK.begin();it!=listOfPatchesOK.end();it++)
981 addPatch((*it)->getConstPart(),factors);
986 * 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.
987 * This method only create patches at level 0 relative to \a this.
989 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const DataArrayByte *criterion, const std::vector<int>& factors)
991 if(!criterion || !criterion->isAllocated())
992 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterion : the criterion DataArrayByte instance must be allocated and not NULL !");
993 std::vector<bool> crit(criterion->toVectorOfBool());//check that criterion has one component.
994 createPatchesFromCriterion(bso,crit,factors);
998 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const DataArrayDouble *criterion, const std::vector<int>& factors, double eps)
1001 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion : null criterion pointer !");
1002 std::vector<bool> inp(criterion->toVectorOfBool(eps));
1003 createPatchesFromCriterion(bso,inp,factors);
1007 * This method creates a multi level patches split at once.
1008 * This method calls as times as size of \a bso createPatchesFromCriterion. Size of \a bso and size of \a factors must be the same !
1011 * \param [in] criterion
1012 * \param [in] factors
1013 * \param [in] eps - See DataArrayDouble::toVectorOfBool for more information about the semantic of eps.
1015 * \sa createPatchesFromCriterion
1017 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterionML(const std::vector<const INTERP_KERNEL::BoxSplittingOptions *>& bso, const DataArrayDouble *criterion, const std::vector< std::vector<int> >& factors, double eps)
1019 std::size_t nbOfLevs(bso.size());
1020 if(nbOfLevs!=factors.size())
1021 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterionML : size of vectors must be the same !");
1025 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterionML : pointers in 1st arg must be not NULL !");
1026 createPatchesFromCriterion(*bso[0],criterion,factors[0],eps);
1027 for(std::size_t i=1;i<nbOfLevs;i++)
1030 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterionML : presence of a NULL BoxSplittingOptions in input vector !");
1032 std::vector<MEDCouplingCartesianAMRPatchGen *> elts(retrieveGridsAt((int)(i)));
1033 std::size_t sz(elts.size());
1034 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGen> > elts2(sz);
1035 std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> > elts3(sz);
1036 for(std::size_t ii=0;ii<sz;ii++)
1039 static const char TMP_STR[]="TMP";
1040 std::vector< std::pair<std::string,int> > fieldNames(1); fieldNames[0].first=TMP_STR; fieldNames[0].second=1;
1041 MEDCouplingAutoRefCountObjectPtr<MEDCouplingAMRAttribute> att(MEDCouplingAMRAttribute::New(this,fieldNames,0));
1043 DataArrayDouble *tmpDa(const_cast<DataArrayDouble *>(att->getFieldOn(this,TMP_STR)));
1044 tmpDa->cpyFrom(*criterion);
1045 att->synchronizeCoarseToFine();
1046 for(std::size_t ii=0;ii<sz;ii++)
1048 const DataArrayDouble *critOnLeaf(att->getFieldOn(const_cast<MEDCouplingCartesianAMRMeshGen *>(elts[ii]->getMesh()),TMP_STR));
1049 elts3[ii]=const_cast<DataArrayDouble *>(critOnLeaf); elts3[ii]->incrRef();
1052 for(std::size_t ii=0;ii<sz;ii++)
1053 const_cast<MEDCouplingCartesianAMRMeshGen *>(elts[ii]->getMesh())->createPatchesFromCriterion(*bso[i],elts3[ii],factors[i],eps);
1057 void MEDCouplingCartesianAMRMeshGen::removeAllPatches()
1063 void MEDCouplingCartesianAMRMeshGen::removePatch(int patchId)
1065 checkPatchId(patchId);
1066 int sz((int)_patches.size()),j(0);
1067 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> > patches(sz-1);
1068 for(int i=0;i<sz;i++)
1070 patches[j++]=_patches[i];
1071 (const_cast<MEDCouplingCartesianAMRMeshGen *>(_patches[patchId]->getMesh()))->detachFromFather();
1076 int MEDCouplingCartesianAMRMeshGen::getNumberOfPatches() const
1078 return (int)_patches.size();
1081 int MEDCouplingCartesianAMRMeshGen::getPatchIdFromChildMesh(const MEDCouplingCartesianAMRMeshGen *mesh) const
1084 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ret++)
1086 if((*it)->getMesh()==mesh)
1089 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPatchIdFromChildMesh : no such a mesh in my direct progeny !");
1092 std::vector< const MEDCouplingCartesianAMRPatch *> MEDCouplingCartesianAMRMeshGen::getPatches() const
1094 std::size_t sz(_patches.size());
1095 std::vector< const MEDCouplingCartesianAMRPatch *> ret(sz);
1096 for(std::size_t i=0;i<sz;i++)
1101 const MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRMeshGen::getPatch(int patchId) const
1103 checkPatchId(patchId);
1104 return _patches[patchId];
1108 * This method states if patch2 (with id \a patchId2) is in the neighborhood of patch1 (with id \a patchId1).
1109 * The neighborhood size is defined by \a ghostLev in the reference of \a this ( \b not in the reference of patches !).
1111 bool MEDCouplingCartesianAMRMeshGen::isPatchInNeighborhoodOf(int patchId1, int patchId2, int ghostLev) const
1113 const MEDCouplingCartesianAMRPatch *p1(getPatch(patchId1)),*p2(getPatch(patchId2));
1114 return p1->isInMyNeighborhood(p2,GetGhostLevelInFineRef(ghostLev,_factors));
1118 * 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.
1119 * 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
1120 * defined by the patch with id \a patchId.
1122 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1123 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1124 * \return DataArrayDouble * - The array of the cell field on the requested patch
1126 * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
1127 * \throw if \a cellFieldOnThis is NULL or not allocated
1128 * \sa fillCellFieldOnPatch, MEDCouplingIMesh::SpreadCoarseToFine
1130 DataArrayDouble *MEDCouplingCartesianAMRMeshGen::createCellFieldOnPatch(int patchId, const DataArrayDouble *cellFieldOnThis) const
1132 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
1133 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
1134 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1135 const MEDCouplingIMesh *fine(patch->getMesh()->getImageMesh());
1136 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(fine->getNumberOfCells(),cellFieldOnThis->getNumberOfComponents());
1137 ret->copyStringInfoFrom(*cellFieldOnThis);
1138 MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),ret,patch->getBLTRRange(),getFactors());
1143 * This method is equivalent to MEDCouplingCartesianAMRMesh::createCellFieldOnPatch except that here instead of creating a new instance
1144 * it fills the value into the \a cellFieldOnPatch data.
1146 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1147 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1148 * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
1150 * \sa createCellFieldOnPatch, fillCellFieldComingFromPatch
1152 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatch(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, bool isConservative) const
1154 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
1155 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
1156 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1157 MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors());
1160 int fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(getFactors()));
1161 std::transform(cellFieldOnPatch->begin(),cellFieldOnPatch->end(),cellFieldOnPatch->getPointer(),std::bind2nd(std::multiplies<double>(),1./((double)fact)));
1166 * This method is the generalization of fillCellFieldOnPatch method. This method only projects coarse to fine without considering the
1167 * potential neighbor patches covered by the ghost cells of patch with id \a patchId.
1169 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1170 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1171 * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
1172 * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
1174 * \sa fillCellFieldOnPatch, fillCellFieldOnPatchGhostAdv
1176 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchGhost(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, int ghostLev, bool isConservative) const
1178 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
1179 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatchGhost : the input cell field array is NULL or not allocated !");
1180 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1181 MEDCouplingIMesh::SpreadCoarseToFineGhost(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),ghostLev);
1184 int fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(getFactors()));
1185 std::transform(cellFieldOnPatch->begin(),cellFieldOnPatch->end(),cellFieldOnPatch->getPointer(),std::bind2nd(std::multiplies<double>(),1./((double)fact)));
1190 * This method is equivalent to fillCellFieldOnPatchGhost except that here \b ONLY \b the \b ghost \b zone will be updated
1191 * in \a cellFieldOnPatch.
1193 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1194 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1195 * \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.
1196 * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
1198 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyOnGhostZone(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, int ghostLev) const
1200 if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
1201 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchOnlyOnGhostZone : the input cell field array is NULL or not allocated !");
1202 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1203 MEDCouplingIMesh::SpreadCoarseToFineGhostZone(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),ghostLev);
1207 * This method is a refinement of fillCellFieldOnPatchGhost. fillCellFieldOnPatchGhost is first called.
1208 * Then for all other patches than those pointed by \a patchId that overlap the ghost zone of the patch impact the ghost zone adequately.
1210 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1211 * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1212 * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
1213 * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
1214 * \param [in] arrsOnPatches - \b WARNING arrsOnPatches[patchId] is \b NOT \b const. All others are const.
1216 * \sa fillCellFieldOnPatchOnlyGhostAdv
1218 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchGhostAdv(int patchId, const DataArrayDouble *cellFieldOnThis, int ghostLev, const std::vector<const DataArrayDouble *>& arrsOnPatches, bool isConservative) const
1220 int nbp(getNumberOfPatches());
1221 if(nbp!=(int)arrsOnPatches.size())
1223 std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchGhostAdv : there are " << nbp << " patches in this and " << arrsOnPatches.size() << " arrays in the last parameter !";
1224 throw INTERP_KERNEL::Exception(oss.str().c_str());
1226 DataArrayDouble *theFieldToFill(const_cast<DataArrayDouble *>(arrsOnPatches[patchId]));
1227 // first, do as usual
1228 fillCellFieldOnPatchGhost(patchId,cellFieldOnThis,theFieldToFill,ghostLev,isConservative);
1229 fillCellFieldOnPatchOnlyGhostAdv(patchId,ghostLev,arrsOnPatches);
1233 * This method updates the patch with id \a patchId considering the only the all the patches in \a this to fill ghost zone.
1234 * So \b warning, the DataArrayDouble instance \a arrsOnPatches[patchId] is non const.
1236 * \sa getPatchIdsInTheNeighborhoodOf
1238 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyGhostAdv(int patchId, int ghostLev, const std::vector<const DataArrayDouble *>& arrsOnPatches) const
1240 int nbp(getNumberOfPatches());
1241 if(nbp!=(int)arrsOnPatches.size())
1243 std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchOnlyGhostAdv : there are " << nbp << " patches in this and " << arrsOnPatches.size() << " arrays in the last parameter !";
1244 throw INTERP_KERNEL::Exception(oss.str().c_str());
1246 const MEDCouplingCartesianAMRPatch *refP(getPatch(patchId));
1247 DataArrayDouble *theFieldToFill(const_cast<DataArrayDouble *>(arrsOnPatches[patchId]));
1248 std::vector<int> ids(getPatchIdsInTheNeighborhoodOf(patchId,ghostLev));
1249 for(std::vector<int>::const_iterator it=ids.begin();it!=ids.end();it++)
1251 const MEDCouplingCartesianAMRPatch *otherP(getPatch(*it));
1252 MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwo(ghostLev,_factors,refP,otherP,theFieldToFill,arrsOnPatches[*it]);
1256 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyOnGhostZoneWith(int ghostLev, const MEDCouplingCartesianAMRPatch *patchToBeModified, const MEDCouplingCartesianAMRPatch *neighborPatch, DataArrayDouble *cellFieldOnPatch, const DataArrayDouble *cellFieldNeighbor) const
1258 MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwo(ghostLev,_factors,patchToBeModified,neighborPatch,cellFieldOnPatch,cellFieldNeighbor);
1262 * This method updates \a cellFieldOnThis part of values coming from the cell field \a cellFieldOnPatch lying on patch having id \a patchId.
1264 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1265 * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
1266 * \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.
1267 * \param [in] isConservative - true if the field needs to be conserved. false if maximum principle has to be applied.
1269 * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
1270 * \throw if \a cellFieldOnPatch is NULL or not allocated
1271 * \sa createCellFieldOnPatch, MEDCouplingIMesh::CondenseFineToCoarse,fillCellFieldComingFromPatchGhost
1273 void MEDCouplingCartesianAMRMeshGen::fillCellFieldComingFromPatch(int patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis, bool isConservative) const
1275 if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
1276 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch : the input cell field array is NULL or not allocated !");
1277 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1278 MEDCouplingIMesh::CondenseFineToCoarse(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis);
1281 int fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(getFactors()));
1282 MEDCouplingStructuredMesh::MultiplyPartOf(_mesh->getCellGridStructure(),patch->getBLTRRange(),1./((double)fact),cellFieldOnThis);
1287 * This method is the extension of MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch managing the ghost cells. If this
1288 * method is called with \a ghostLev equal to 0 it behaves exactly as MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch.
1290 * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1291 * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
1292 * \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.
1293 * \param [in] ghostLev The size of ghost zone (must be >= 0 !)
1294 * \param [in] isConservative - true if the field needs to be conserved. false if maximum principle has to be applied.
1296 * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
1297 * \throw if \a cellFieldOnPatch is NULL or not allocated
1298 * \sa fillCellFieldComingFromPatch
1300 void MEDCouplingCartesianAMRMeshGen::fillCellFieldComingFromPatchGhost(int patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis, int ghostLev, bool isConservative) const
1302 if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
1303 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatchGhost : the input cell field array is NULL or not allocated !");
1304 const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1305 MEDCouplingIMesh::CondenseFineToCoarseGhost(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis,ghostLev);
1308 int fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(getFactors()));
1309 MEDCouplingStructuredMesh::MultiplyPartOfByGhost(_mesh->getCellGridStructure(),patch->getBLTRRange(),ghostLev,1./((double)fact),cellFieldOnThis);
1314 * This method finds all patches (located by their ids) that are in the neighborhood of patch with id \a patchId. The neighborhood size is
1315 * defined by ghostLev.
1317 * \param [in] patchId - the id of the considered patch.
1318 * \param [in] ghostLev - the size of the neighborhood.
1319 * \return DataArrayInt * - the newly allocated array containing the list of patches in the neighborhood of the considered patch. This array is to be deallocated by the caller.
1321 DataArrayInt *MEDCouplingCartesianAMRMeshGen::findPatchesInTheNeighborhoodOf(int patchId, int ghostLev) const
1323 int nbp(getNumberOfPatches());
1324 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
1325 for(int i=0;i<nbp;i++)
1328 if(isPatchInNeighborhoodOf(i,patchId,ghostLev))
1329 ret->pushBackSilent(i);
1334 MEDCouplingUMesh *MEDCouplingCartesianAMRMeshGen::buildUnstructured() const
1336 MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> part(_mesh->buildUnstructured());
1337 std::vector<bool> bs(_mesh->getNumberOfCells(),false);
1338 std::vector<int> cgs(_mesh->getCellGridStructure());
1339 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > msSafe(_patches.size()+1);
1341 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
1343 MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
1344 msSafe[ii+1]=(*it)->getMesh()->buildUnstructured();
1346 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsOff(DataArrayInt::BuildListOfSwitchedOff(bs));
1347 msSafe[0]=static_cast<MEDCouplingUMesh *>(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
1348 std::vector< const MEDCouplingUMesh * > ms(msSafe.size());
1349 for(std::size_t i=0;i<msSafe.size();i++)
1351 return MEDCouplingUMesh::MergeUMeshes(ms);
1355 * This method returns a mesh containing as cells that there is patches at the current level.
1356 * The patches are seen like 'boxes' that is too say the refinement will not appear here.
1358 * \return MEDCoupling1SGTUMesh * - A new object to be managed by the caller containing as cells as there are patches in \a this.
1360 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMeshGen::buildMeshFromPatchEnvelop() const
1362 std::vector<const MEDCoupling1SGTUMesh *> cells;
1363 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> > cellsSafe;
1364 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1366 const MEDCouplingCartesianAMRPatch *patch(*it);
1369 MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> cell(patch->getMesh()->getImageMesh()->asSingleCell());
1370 MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> cell1SGT(cell->build1SGTUnstructured());
1371 cellsSafe.push_back(cell1SGT); cells.push_back(cell1SGT);
1374 return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(cells);
1377 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMeshGen::buildMeshOfDirectChildrenOnly() const
1379 std::vector<const MEDCoupling1SGTUMesh *> patches;
1380 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> > patchesSafe;
1381 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1383 const MEDCouplingCartesianAMRPatch *patch(*it);
1386 MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> patchMesh(patch->getMesh()->getImageMesh()->build1SGTUnstructured());
1387 patchesSafe.push_back(patchMesh); patches.push_back(patchMesh);
1390 return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(patches);
1394 * This method works same as buildUnstructured except that arrays are given in input to build a field on cell in output.
1395 * \return MEDCouplingFieldDouble * - a newly created instance the caller has reponsability to deal with.
1396 * \sa buildUnstructured
1398 MEDCouplingFieldDouble *MEDCouplingCartesianAMRMeshGen::buildCellFieldOnRecurseWithoutOverlapWithoutGhost(int ghostSz, const std::vector<const DataArrayDouble *>& recurseArrs) const
1400 if(recurseArrs.empty())
1401 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::buildCellFieldOnRecurseWithoutOverlapWithoutGhost : array is empty ! Should never happen !");
1403 std::vector<bool> bs(_mesh->getNumberOfCells(),false);
1404 std::vector<int> cgs(_mesh->getCellGridStructure());
1405 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> > msSafe(_patches.size()+1);
1407 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
1409 MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
1410 std::vector<const DataArrayDouble *> tmpArrs(extractSubTreeFromGlobalFlatten((*it)->getMesh(),recurseArrs));
1411 msSafe[ii+1]=(*it)->getMesh()->buildCellFieldOnRecurseWithoutOverlapWithoutGhost(ghostSz,tmpArrs);
1413 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsOff(DataArrayInt::BuildListOfSwitchedOff(bs));
1415 MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret(MEDCouplingFieldDouble::New(ON_CELLS));
1416 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr2(extractGhostFrom(ghostSz,recurseArrs[0]));
1417 arr2=arr2->selectByTupleIdSafe(eltsOff->begin(),eltsOff->end());
1418 ret->setArray(arr2);
1419 ret->setName(arr2->getName());
1420 MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> part(_mesh->buildUnstructured());
1421 MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> mesh(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
1425 std::vector< const MEDCouplingFieldDouble * > ms(msSafe.size());
1426 for(std::size_t i=0;i<msSafe.size();i++)
1429 return MEDCouplingFieldDouble::MergeFields(ms);
1433 * This method extracts from \arr arr the part inside \a arr by cutting the \a ghostSz external part.
1434 * \arr is expected to be an array having a number of tuples equal to \c getImageMesh()->buildWithGhost(ghostSz).
1436 DataArrayDouble *MEDCouplingCartesianAMRMeshGen::extractGhostFrom(int ghostSz, const DataArrayDouble *arr) const
1438 std::vector<int> st(_mesh->getCellGridStructure());
1439 std::vector< std::pair<int,int> > p(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(st));
1440 std::transform(st.begin(),st.end(),st.begin(),std::bind2nd(std::plus<int>(),2*ghostSz));
1441 MEDCouplingStructuredMesh::ApplyGhostOnCompactFrmt(p,ghostSz);
1442 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(st,arr,p));
1447 * This method returns all the patches in \a this not equal to \a patchId that are in neighborhood of patch with id \a patchId.
1449 * \sa fillCellFieldOnPatchOnlyGhostAdv
1451 std::vector<int> MEDCouplingCartesianAMRMeshGen::getPatchIdsInTheNeighborhoodOf(int patchId, int ghostLev) const
1453 std::vector<int> ret;
1454 int nbp(getNumberOfPatches());
1456 for(int i=0;i<nbp;i++)
1459 if(isPatchInNeighborhoodOf(i,patchId,ghostLev))
1466 * This method returns a dump python of \a this. It is useful for users of createPatchesFromCriterion method for debugging.
1468 * \sa dumpPatchesOf, createPatchesFromCriterion, createPatchesFromCriterionML
1470 std::string MEDCouplingCartesianAMRMeshGen::buildPythonDumpOfThis() const
1472 std::ostringstream oss;
1473 oss << "amr=MEDCouplingCartesianAMRMesh(\""<< getImageMesh()->getName() << "\"," << getSpaceDimension() << ",[";
1474 std::vector<int> ngs(getImageMesh()->getNodeGridStructure());
1475 std::vector<double> orig(getImageMesh()->getOrigin()),dxyz(getImageMesh()->getDXYZ());
1476 std::copy(ngs.begin(),ngs.end(),std::ostream_iterator<int>(oss,","));
1478 std::copy(orig.begin(),orig.end(),std::ostream_iterator<double>(oss,","));
1480 std::copy(dxyz.begin(),dxyz.end(),std::ostream_iterator<double>(oss,","));
1482 dumpPatchesOf("amr",oss);
1486 MEDCouplingCartesianAMRMeshGen::MEDCouplingCartesianAMRMeshGen(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
1487 const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):_father(0)
1489 _mesh=MEDCouplingIMesh::New(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
1492 MEDCouplingCartesianAMRMeshGen::MEDCouplingCartesianAMRMeshGen(MEDCouplingCartesianAMRMeshGen *father, MEDCouplingIMesh *mesh):_father(father)
1495 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen(MEDCouplingIMesh *mesh) constructor : empty father !");
1497 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen(MEDCouplingIMesh *mesh) constructor : The input mesh is null !");
1498 mesh->checkCoherency();
1499 _mesh=mesh; _mesh->incrRef();
1502 void MEDCouplingCartesianAMRMeshGen::checkPatchId(int patchId) const
1504 int sz(getNumberOfPatches());
1505 if(patchId<0 || patchId>=sz)
1507 std::ostringstream oss; oss << "MEDCouplingCartesianAMRMeshGen::checkPatchId : invalid patchId (" << patchId << ") ! Must be in [0," << sz << ") !";
1508 throw INTERP_KERNEL::Exception(oss.str().c_str());
1512 void MEDCouplingCartesianAMRMeshGen::checkFactorsAndIfNotSetAssign(const std::vector<int>& factors)
1514 if(getSpaceDimension()!=(int)factors.size())
1515 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::checkFactorsAndIfNotSetAssign : invalid size of factors ! size must be equal to the spaceDimension !");
1516 if(_factors.empty())
1522 if(_factors!=factors)
1523 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::checkFactorsAndIfNotSetAssign : the factors ");
1527 void MEDCouplingCartesianAMRMeshGen::retrieveGridsAtInternal(int lev, std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGen> >& grids) const
1531 const MEDCouplingCartesianAMRMesh *thisc(dynamic_cast<const MEDCouplingCartesianAMRMesh *>(this));//tony
1532 MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGF> elt(new MEDCouplingCartesianAMRPatchGF(const_cast<MEDCouplingCartesianAMRMesh *>(thisc)));
1533 grids.push_back(DynamicCastSafe<MEDCouplingCartesianAMRPatchGF,MEDCouplingCartesianAMRPatchGen>(elt));
1537 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1539 const MEDCouplingCartesianAMRPatch *pt(*it);
1542 MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> tmp1(*it);
1543 grids.push_back(DynamicCastSafe<MEDCouplingCartesianAMRPatch,MEDCouplingCartesianAMRPatchGen>(tmp1));
1549 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1551 const MEDCouplingCartesianAMRPatch *pt(*it);
1553 pt->getMesh()->retrieveGridsAtInternal(lev-1,grids);
1558 int MEDCouplingCartesianAMRMeshGen::GetGhostLevelInFineRef(int ghostLev, const std::vector<int>& factors)
1561 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::GetGhostLevelInFineRef : the ghost size must be >=0 !");
1563 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::GetGhostLevelInFineRef : no factors defined !");
1564 int ghostLevInPatchRef;
1566 ghostLevInPatchRef=0;
1569 ghostLevInPatchRef=(ghostLev-1)/factors[0]+1;
1570 for(std::size_t i=0;i<factors.size();i++)
1571 ghostLevInPatchRef=std::max(ghostLevInPatchRef,(ghostLev-1)/factors[i]+1);
1573 return ghostLevInPatchRef;
1577 * This method returns a sub set of \a all. The subset is defined by the \a head in the tree defined by \a this.
1578 * Elements in \a all are expected to be sorted from god father to most refined structure.
1580 std::vector<const DataArrayDouble *> MEDCouplingCartesianAMRMeshGen::extractSubTreeFromGlobalFlatten(const MEDCouplingCartesianAMRMeshGen *head, const std::vector<const DataArrayDouble *>& all) const
1582 int maxLev(getMaxNumberOfLevelsRelativeToThis());
1583 std::vector<const DataArrayDouble *> ret;
1584 std::vector<const MEDCouplingCartesianAMRMeshGen *> meshes(1,this);
1586 for(int i=0;i<maxLev;i++)
1588 std::vector<const MEDCouplingCartesianAMRMeshGen *> meshesTmp;
1589 for(std::vector<const MEDCouplingCartesianAMRMeshGen *>::const_iterator it=meshes.begin();it!=meshes.end();it++)
1591 if((*it)==head || head->isObjectInTheProgeny(*it))
1592 ret.push_back(all[kk]);
1594 std::vector< const MEDCouplingCartesianAMRPatch *> ps((*it)->getPatches());
1595 for(std::vector< const MEDCouplingCartesianAMRPatch *>::const_iterator it0=ps.begin();it0!=ps.end();it0++)
1597 const MEDCouplingCartesianAMRMeshGen *mesh((*it0)->getMesh());
1598 meshesTmp.push_back(mesh);
1604 throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::extractSubTreeFromGlobalFlatten : the size of input vector is not compatible with number of leaves in this !");
1608 void MEDCouplingCartesianAMRMeshGen::dumpPatchesOf(const std::string& varName, std::ostream& oss) const
1611 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1613 const MEDCouplingCartesianAMRPatch *patch(*it);
1616 std::ostringstream oss2; oss2 << varName << ".addPatch([";
1617 const std::vector< std::pair<int,int> >& bltr(patch->getBLTRRange());
1618 std::size_t sz(bltr.size());
1619 for(std::size_t i=0;i<sz;i++)
1621 oss2 << "(" << bltr[i].first << "," << bltr[i].second << ")";
1626 std::copy(_factors.begin(),_factors.end(),std::ostream_iterator<int>(oss2,","));
1629 std::ostringstream oss3; oss3 << varName << "[" << j++ << "]";
1630 patch->getMesh()->dumpPatchesOf(oss3.str(),oss);
1635 std::size_t MEDCouplingCartesianAMRMeshGen::getHeapMemorySizeWithoutChildren() const
1637 return sizeof(MEDCouplingCartesianAMRMeshGen);
1640 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMeshGen::getDirectChildren() const
1642 std::vector<const BigMemoryObject *> ret;
1643 if((const MEDCouplingIMesh *)_mesh)
1644 ret.push_back((const MEDCouplingIMesh *)_mesh);
1645 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1647 if((const MEDCouplingCartesianAMRPatch*)*it)
1648 ret.push_back((const MEDCouplingCartesianAMRPatch*)*it);
1653 void MEDCouplingCartesianAMRMeshGen::updateTime() const
1655 if((const MEDCouplingIMesh *)_mesh)
1656 updateTimeWith(*_mesh);
1657 for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1659 const MEDCouplingCartesianAMRPatch *elt(*it);
1662 const MEDCouplingCartesianAMRMeshGen *mesh(elt->getMesh());
1664 updateTimeWith(*mesh);
1668 MEDCouplingCartesianAMRMeshSub::MEDCouplingCartesianAMRMeshSub(MEDCouplingCartesianAMRMeshGen *father, MEDCouplingIMesh *mesh):MEDCouplingCartesianAMRMeshGen(father,mesh)
1672 MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::New(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
1673 const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop)
1675 return new MEDCouplingCartesianAMRMesh(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
1678 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
1679 const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):MEDCouplingCartesianAMRMeshGen(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop)
1683 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMesh::getDirectChildren() const
1685 std::vector<const BigMemoryObject *> ret(MEDCouplingCartesianAMRMeshGen::getDirectChildren());
1689 MEDCouplingCartesianAMRMesh::~MEDCouplingCartesianAMRMesh()