Salome HOME
Debug on condensation with ghostSz>1 + more synchronization methods.
[modules/med.git] / src / MEDCoupling / MEDCouplingCartesianAMRMesh.cxx
1 // Copyright (C) 2007-2014  CEA/DEN, EDF R&D
2 //
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.
7 //
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.
12 //
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
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 // Author : Anthony Geay
20
21 #include "MEDCouplingCartesianAMRMesh.hxx"
22 #include "MEDCouplingAMRAttribute.hxx"
23 #include "MEDCouplingFieldDouble.hxx"
24 #include "MEDCoupling1GTUMesh.hxx"
25 #include "MEDCouplingIMesh.hxx"
26 #include "MEDCouplingUMesh.hxx"
27
28 #include <limits>
29 #include <sstream>
30 #include <numeric>
31
32 using namespace ParaMEDMEM;
33
34 /// @cond INTERNAL
35
36 int MEDCouplingCartesianAMRPatchGen::getNumberOfCellsRecursiveWithOverlap() const
37 {
38   return _mesh->getNumberOfCellsRecursiveWithOverlap();
39 }
40
41 int MEDCouplingCartesianAMRPatchGen::getNumberOfCellsRecursiveWithoutOverlap() const
42 {
43   return _mesh->getNumberOfCellsRecursiveWithoutOverlap();
44 }
45
46 int MEDCouplingCartesianAMRPatchGen::getMaxNumberOfLevelsRelativeToThis() const
47 {
48   return _mesh->getMaxNumberOfLevelsRelativeToThis();
49 }
50
51 MEDCouplingCartesianAMRPatchGen::MEDCouplingCartesianAMRPatchGen(MEDCouplingCartesianAMRMeshGen *mesh)
52 {
53   if(!mesh)
54     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen constructor : input mesh is NULL !");
55   _mesh=mesh; _mesh->incrRef();
56 }
57
58 MEDCouplingCartesianAMRPatchGen::MEDCouplingCartesianAMRPatchGen(const MEDCouplingCartesianAMRPatchGen& other):RefCountObject(other),_mesh(other._mesh)
59 {
60   const MEDCouplingCartesianAMRMeshGen *mesh(other._mesh);
61   if(mesh)
62     _mesh=mesh->deepCpy();
63 }
64
65 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatchGen::getMeshSafe() const
66 {
67   const MEDCouplingCartesianAMRMeshGen *mesh(_mesh);
68   if(!mesh)
69     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen::getMeshSafe const : the mesh is NULL !");
70   return mesh;
71 }
72
73 MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatchGen::getMeshSafe()
74 {
75   MEDCouplingCartesianAMRMeshGen *mesh(_mesh);
76     if(!mesh)
77       throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen::getMeshSafe : the mesh is NULL !");
78     return mesh;
79 }
80
81 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRPatchGen::getDirectChildren() const
82 {
83   std::vector<const BigMemoryObject *> ret;
84   if((const MEDCouplingCartesianAMRMeshGen *)_mesh)
85     ret.push_back((const MEDCouplingCartesianAMRMeshGen *)_mesh);
86   return ret;
87 }
88
89 /*!
90  * \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.
91  * \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,
92  *                                a the end cell (\b excluded) of the range for the second element of the pair.
93  */
94 MEDCouplingCartesianAMRPatch::MEDCouplingCartesianAMRPatch(MEDCouplingCartesianAMRMeshGen *mesh, const std::vector< std::pair<int,int> >& bottomLeftTopRight):MEDCouplingCartesianAMRPatchGen(mesh),_bl_tr(bottomLeftTopRight)
95 {
96   int dim((int)bottomLeftTopRight.size()),dimExp(_mesh->getSpaceDimension());
97   if(dim!=dimExp)
98     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch constructor : space dimension of father and input bottomLeft/topRight size mismatches !");
99 }
100
101 MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRPatch::deepCpy() const
102 {
103   return new MEDCouplingCartesianAMRPatch(*this);
104 }
105
106 void MEDCouplingCartesianAMRPatch::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
107 {
108   return getMeshSafe()->addPatch(bottomLeftTopRight,factors);
109 }
110
111 int MEDCouplingCartesianAMRPatch::getNumberOfOverlapedCellsForFather() const
112 {
113   return MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(_bl_tr);
114 }
115
116 /*!
117  * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
118  * 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 !).
119  * Call isInMyNeighborhoodExt to deal with 2 patches not sharing directly the same father.
120  *
121  * \param [in] other - The other patch
122  * \param [in] ghostLev - The size of the neighborhood zone.
123  *
124  * \throw if \a this or \a other are invalid (end before start).
125  * \throw if \a ghostLev is \b not >= 0 .
126  * \throw if \a this and \a other have not the same space dimension.
127  *
128  * \sa isInMyNeighborhoodExt
129  */
130 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhood(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const
131 {
132   if(ghostLev<0)
133     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the size of the neighborhood must be >= 0 !");
134   if(!other)
135     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the input patch is NULL !");
136   const std::vector< std::pair<int,int> >& thisp(getBLTRRange());
137   const std::vector< std::pair<int,int> >& otherp(other->getBLTRRange());
138   return IsInMyNeighborhood(ghostLev,thisp,otherp);
139 }
140
141 /*!
142  * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
143  * 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
144  * ancestor must be the same. If they don't have the same ancestor an exception will be thrown.
145  *
146  * \param [in] other - The other patch
147  * \param [in] ghostLev - The size of the neighborhood zone.
148  *
149  * \throw if \a this or \a other are invalid (end before start).
150  * \throw if \a ghostLev is \b not >= 0 .
151  * \throw if \a this and \a other have not the same space dimension.
152  * \throw if there is not common ancestor of \a this and \a other.
153  *
154  * \sa isInMyNeighborhood, isInMyNeighborhoodDiffLev
155  */
156 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const
157 {
158   if(ghostLev<0)
159     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt : the size of the neighborhood must be >= 0 !");
160   if(!other)
161     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt : the input patch is NULL !");
162   int lev;
163   const MEDCouplingCartesianAMRMeshGen *com(FindCommonAncestor(this,other,lev));//check that factors are OK
164   if(lev==0)
165     return isInMyNeighborhood(other,ghostLev);
166   std::vector<int> offset(ComputeOffsetFromTwoToOne(com,lev,this,other));
167   const std::vector< std::pair<int,int> >& thisp(getBLTRRange());
168   std::vector< std::pair<int,int> > otherp(other->getBLTRRange());
169   otherp=MEDCouplingStructuredMesh::TranslateCompactFrmt(otherp,offset);
170   return IsInMyNeighborhood(ghostLev,thisp,otherp);
171 }
172
173 /*!
174  * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
175  * the must be >= 0. This method works even if \a this and \a other does not share the same father.
176  * \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  *
178  * \param [in] other - The other patch
179  * \param [in] ghostLev - The size of the neighborhood zone.
180  *
181  * \throw if \a this or \a other are invalid (end before start).
182  * \throw if \a ghostLev is \b not >= 0 .
183  * \throw if \a this and \a other have not the same space dimension.
184  * \throw if there is not common ancestor of \a this and \a other.
185  *
186  * \sa isInMyNeighborhoodExt
187  */
188 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhoodDiffLev(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const
189 {
190   std::vector< std::pair<int,int> > thispp,otherpp;
191   std::vector<int> factors;
192   ComputeZonesOfTwoRelativeToOneDiffLev(ghostLev,this,other,thispp,otherpp,factors);
193   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
194 }
195
196 std::vector<int> MEDCouplingCartesianAMRPatch::computeCellGridSt() const
197 {
198   const MEDCouplingCartesianAMRMeshGen *m(getMesh());
199   if(!m)
200     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::computeCellGridSt : no mesh held by this !");
201   const MEDCouplingCartesianAMRMeshGen *father(m->getFather());
202   if(!father)
203     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::computeCellGridSt : no father help by underlying mesh !");
204   const std::vector< std::pair<int,int> >& bltr(getBLTRRange());
205   const std::vector<int>& factors(father->getFactors());
206   std::vector<int> ret(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(bltr));
207   std::transform(ret.begin(),ret.end(),factors.begin(),ret.begin(),std::multiplies<int>());
208   return ret;
209 }
210
211 bool MEDCouplingCartesianAMRPatch::IsInMyNeighborhood(int ghostLev, const std::vector< std::pair<int,int> >& p1, const std::vector< std::pair<int,int> >& p2)
212 {
213   std::size_t thispsize(p1.size());
214   if(thispsize!=p2.size())
215     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the dimensions must be the same !");
216   for(std::size_t i=0;i<thispsize;i++)
217     {
218       const std::pair<int,int>& thispp(p1[i]);
219       const std::pair<int,int>& otherpp(p2[i]);
220       if(thispp.second<thispp.first)
221         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : this patch is invalid !");
222       if(otherpp.second<otherpp.first)
223         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : this patch is invalid !");
224       if(otherpp.first==thispp.second+ghostLev-1)
225         continue;
226       if(otherpp.second+ghostLev-1==thispp.first)
227         continue;
228       int start(std::max(thispp.first,otherpp.first)),end(std::min(thispp.second,otherpp.second));
229       if(end<start)
230         return false;
231     }
232   return true;
233 }
234
235 /*!
236  * \sa FindNeighborsOfSubPatchesOf
237  */
238 std::vector< std::vector< std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *> > > MEDCouplingCartesianAMRPatch::FindNeighborsOfSubPatchesOfSameLev(int ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2)
239 {
240   if(!p1 || !p2)
241     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindNeighborsOfSubPatchesOfSameLev : the input pointers must be not NULL !");
242   std::vector< std::vector< std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *> > > ret;
243   std::vector< const MEDCouplingCartesianAMRPatch *> p1Work(p1->getMesh()->getPatches()),p2Work(p2->getMesh()->getPatches());
244   while(!p1Work.empty())
245     {
246       std::vector< std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *> > retTmp;
247       std::vector<const MEDCouplingCartesianAMRPatch *> p1Work2,p2Work2;
248       for(std::vector<const MEDCouplingCartesianAMRPatch *>::const_iterator it1=p1Work.begin();it1!=p1Work.end();it1++)
249         {
250           for(std::vector<const MEDCouplingCartesianAMRPatch *>::const_iterator it2=p2Work.begin();it2!=p2Work.end();it2++)
251             {
252               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
253                 retTmp.push_back(std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *>(*it1,*it2));
254             }
255           std::vector<const MEDCouplingCartesianAMRPatch *> tmp1((*it1)->getMesh()->getPatches());
256           p1Work2.insert(p1Work2.end(),tmp1.begin(),tmp1.end());
257         }
258       for(std::vector<const MEDCouplingCartesianAMRPatch *>::const_iterator it2=p2Work.begin();it2!=p2Work.end();it2++)
259         {
260           std::vector<const MEDCouplingCartesianAMRPatch *> tmp2((*it2)->getMesh()->getPatches());
261           p2Work2.insert(p2Work2.end(),tmp2.begin(),tmp2.end());
262         }
263       ret.push_back(retTmp);
264       p1Work=p1Work2;
265       p2Work=p2Work2;
266     }
267   return ret;
268 }
269
270 /*!
271  * This method returns all pair of patches (pa,pb) so that pb is in the neighborhood of pa (size of neighborhood is \a ghostLev).
272  * 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
273  * FindNeighborsOfSubPatchesOfSameLev.
274  *
275  * \sa FindNeighborsOfSubPatchesOfSameLev
276  */
277 void MEDCouplingCartesianAMRPatch::FindNeighborsOfSubPatchesOf(int ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, std::vector< std::pair<const MEDCouplingCartesianAMRPatch *, const MEDCouplingCartesianAMRPatch *> >& ret)
278 {
279   if(!p1 || !p2)
280     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindNeighborsOfSubPatchesOf : the input pointers must be not NULL !");
281   std::vector< const MEDCouplingCartesianAMRPatch *> p1Work(p1->getMesh()->getPatches());
282   while(!p1Work.empty())
283     {
284       std::vector<const MEDCouplingCartesianAMRPatch *> p1Work2;
285       for(std::vector<const MEDCouplingCartesianAMRPatch *>::const_iterator it0=p1Work.begin();it0!=p1Work.end();it0++)
286         {
287           if((*it0)->isInMyNeighborhoodDiffLev(p2,ghostLev))
288             ret.push_back(std::pair<const MEDCouplingCartesianAMRPatch *,const MEDCouplingCartesianAMRPatch *>(*it0,p2));
289           std::vector<const MEDCouplingCartesianAMRPatch *> tmp2((*it0)->getMesh()->getPatches());
290           p1Work2.insert(p1Work2.end(),tmp2.begin(),tmp2.end());
291         }
292       p1Work=p1Work2;
293     }
294 }
295
296 /*!
297  * \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.
298  *
299  * \saUpdateNeighborsOfOneWithTwoExt
300  */
301 void MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwo(int ghostLev, const std::vector<int>& factors, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, DataArrayDouble *dataOnP1, const DataArrayDouble *dataOnP2)
302 {
303   const std::vector< std::pair<int,int> >& p1BLTR(p1->getBLTRRange());
304   const std::vector< std::pair<int,int> >& p2BLTR(p2->getBLTRRange());
305   UpdateNeighborsOfOneWithTwoInternal(ghostLev,factors,p1BLTR,p2BLTR,dataOnP1,dataOnP2);
306 }
307
308 /*!
309  * Idem than UpdateNeighborsOfOneWithTwo, except that here \a p1 and \a p2 are not sharing the same direct father.
310  *
311  * \sa UpdateNeighborsOfOneWithTwo
312  */
313 void MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwoExt(int ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, DataArrayDouble *dataOnP1, const DataArrayDouble *dataOnP2)
314 {
315   const std::vector< std::pair<int,int> >& p1BLTR(p1->getBLTRRange());//p1BLTR=[(10,12),(5,8)]
316   std::vector< std::pair<int,int> > p2BLTR(p2->getBLTRRange());//p2BLTR=[(0,1),(0,5)]
317   int lev(0);
318   const MEDCouplingCartesianAMRMeshGen *ca(FindCommonAncestor(p1,p2,lev));
319   std::vector<int> offset(ComputeOffsetFromTwoToOne(ca,lev,p1,p2));//[12,4]
320   p2BLTR=MEDCouplingStructuredMesh::TranslateCompactFrmt(p2BLTR,offset);//p2BLTR=[(12,13),(4,9)]
321   UpdateNeighborsOfOneWithTwoInternal(ghostLev,p1->getMesh()->getFather()->getFactors(),p1BLTR,p2BLTR,dataOnP1,dataOnP2);
322 }
323
324 /*!
325  * \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 !
326  */
327 void MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwoMixedLev(int ghostLev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, DataArrayDouble *dataOnP1, const DataArrayDouble *dataOnP2)
328 {
329   std::vector< std::pair<int,int> > p1pp,p2pp;
330   std::vector<int> factors;
331   ComputeZonesOfTwoRelativeToOneDiffLev(ghostLev,p1,p2,p1pp,p2pp,factors);
332   //
333   std::vector<int> dimsP2NotRefined(p2->computeCellGridSt());
334   std::vector<int> dimsP2Refined(dimsP2NotRefined);
335   std::transform(dimsP2NotRefined.begin(),dimsP2NotRefined.end(),factors.begin(),dimsP2Refined.begin(),std::multiplies<int>());
336   std::vector< std::pair<int,int> > p2RefinedAbs(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(dimsP2NotRefined));
337   std::vector<int> dimsP2RefinedGhost(dimsP2Refined.size());
338   std::transform(dimsP2Refined.begin(),dimsP2Refined.end(),dimsP2RefinedGhost.begin(),std::bind2nd(std::plus<int>(),2*ghostLev));
339   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> fineP2(DataArrayDouble::New()); fineP2->alloc(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(dimsP2RefinedGhost),dataOnP2->getNumberOfComponents());
340   MEDCouplingIMesh::SpreadCoarseToFineGhost(dataOnP2,dimsP2NotRefined,fineP2,p2RefinedAbs,factors,ghostLev);
341   //
342   UpdateNeighborsOfOneWithTwoInternal(ghostLev,p1->getMesh()->getFather()->getFactors(),p1pp,p2pp,dataOnP1,fineP2);
343 }
344
345 /*!
346  * \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 !
347  * 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
348  * on the same level as \a p1.
349  */
350 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)
351 {
352   std::vector<const MEDCouplingCartesianAMRMeshGen *> ancestorsOfThis;
353   const MEDCouplingCartesianAMRMeshGen *work(p1->getMesh()),*work2(0);
354   ancestorsOfThis.push_back(work);
355   while(work)
356     {
357       work=work->getFather();
358       if(work)
359         ancestorsOfThis.push_back(work);
360     }
361   //
362   work=p2->getMesh();
363   bool found(false);
364   std::size_t levThis(0),levOther(0);
365   while(work && !found)
366     {
367       work2=work;
368       work=work->getFather();
369       if(work)
370         {
371           levOther++;
372           std::vector<const MEDCouplingCartesianAMRMeshGen *>::iterator it(std::find(ancestorsOfThis.begin(),ancestorsOfThis.end(),work));
373           if(it!=ancestorsOfThis.end())
374             {
375               levThis=std::distance(ancestorsOfThis.begin(),it);
376               found=true;
377             }
378         }
379     }
380   if(!found)
381     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeZonesOfTwoRelativeToOneDiffLev : no common ancestor found !");
382   if(levThis<=levOther)
383     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeZonesOfTwoRelativeToOneDiffLev : this method is not called correctly !");
384   //
385   const MEDCouplingCartesianAMRMeshGen *comAncestor(ancestorsOfThis[levThis]);
386   int idThis(comAncestor->getPatchIdFromChildMesh(ancestorsOfThis[levThis-1])),idOther(comAncestor->getPatchIdFromChildMesh(work2));
387   const MEDCouplingCartesianAMRPatch *thisp(comAncestor->getPatch(idThis)),*otherp(comAncestor->getPatch(idOther));
388   std::vector<int> offset(ComputeOffsetFromTwoToOne(comAncestor,levOther,thisp,otherp));
389   p1Zone=thisp->getBLTRRange(); p2Zone=MEDCouplingStructuredMesh::TranslateCompactFrmt(otherp->getBLTRRange(),offset);
390   factToApplyOn2.resize(p1Zone.size()); std::fill(factToApplyOn2.begin(),factToApplyOn2.end(),1);
391   //
392   std::size_t nbOfTurn(levThis-levOther);
393   for(std::size_t i=0;i<nbOfTurn;i++)
394     {
395       std::vector< std::pair<int,int> > tmp0;
396       MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(p1Zone,p2Zone,tmp0,false);
397       p2Zone=tmp0;
398       const MEDCouplingCartesianAMRMeshGen *curAncestor(ancestorsOfThis[levThis-i]);
399       ApplyFactorsOnCompactFrmt(p2Zone,curAncestor->getFactors());
400       curAncestor=ancestorsOfThis[levThis-1-i];
401       const std::vector<int>& factors(curAncestor->getFactors());
402       std::transform(factToApplyOn2.begin(),factToApplyOn2.end(),factors.begin(),factToApplyOn2.begin(),std::multiplies<int>());
403       int tmpId(curAncestor->getPatchIdFromChildMesh(ancestorsOfThis[levThis-2-i]));
404       p1Zone=curAncestor->getPatch(tmpId)->getBLTRRange();
405     }
406 }
407
408 std::size_t MEDCouplingCartesianAMRPatch::getHeapMemorySizeWithoutChildren() const
409 {
410   std::size_t ret(sizeof(MEDCouplingCartesianAMRPatch));
411   ret+=_bl_tr.capacity()*sizeof(std::pair<int,int>);
412   return ret;
413 }
414
415 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatch::FindCommonAncestor(const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, int& lev)
416 {
417   const MEDCouplingCartesianAMRMeshGen *f1(p1->_mesh),*f2(p2->_mesh);
418   lev=0;
419   while(f1!=f2 || f1==0 || f2==0)
420     {
421       f1=f1->getFather(); f2=f2->getFather();
422       if(f1->getFactors()!=f2->getFactors())
423         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindCommonAncestor : factors differ !");
424       lev++;
425     }
426   if(f1!=f2)
427     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindCommonAncestor : no common ancestor between p1 and p2 !");
428   return f1;
429 }
430
431 std::vector<int> MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne(const MEDCouplingCartesianAMRMeshGen *comAncestor, int lev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2)
432 {
433   if(lev<=0)
434     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne : this method is useful only for lev > 0 !");
435   int zeLev(lev-1);
436   int dim(p1->getMesh()->getSpaceDimension());
437   if(p2->getMesh()->getSpaceDimension()!=dim)
438     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne : dimension must be the same !");
439   std::vector< int > ret(dim,0);
440   for(int i=0;i<zeLev;i++)
441     {
442       const MEDCouplingCartesianAMRMeshGen *f1(p1->_mesh),*f2(p2->_mesh);
443       const MEDCouplingCartesianAMRPatch *p1h(0),*p2h(0);
444       for(int j=0;j<lev-i;j++)
445         {
446           const MEDCouplingCartesianAMRMeshGen *f1tmp(f1->getFather()),*f2tmp(f2->getFather());
447           int pid1(f1tmp->getPatchIdFromChildMesh(f1)),pid2(f2tmp->getPatchIdFromChildMesh(f2));
448           p1h=f1tmp->getPatch(pid1); p2h=f2tmp->getPatch(pid2);
449           f1=f1tmp; f2=f2tmp;
450         }
451       std::vector< std::pair<int,int> > p2c(p2h->getBLTRRange());
452       for(int k=0;k<dim;k++)
453         {
454           p2c[k].first+=ret[k];
455           p2c[k].second+=ret[k];
456         }
457       for(int k=0;k<dim;k++)
458         {
459           ret[k]=p2c[k].first-p1h->getBLTRRange()[k].first;
460           ret[k]*=f1->getFactors()[k];
461         }
462     }
463   return ret;
464 }
465
466 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)
467 {//p1=[(1,4),(2,4)] p2=[(4,5),(3,4)]
468   int dim((int)factors.size());
469   std::vector<int> dimsCoarse(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(p1));//[3,2]
470   std::transform(dimsCoarse.begin(),dimsCoarse.end(),factors.begin(),dimsCoarse.begin(),std::multiplies<int>());//[12,8]
471   std::transform(dimsCoarse.begin(),dimsCoarse.end(),dimsCoarse.begin(),std::bind2nd(std::plus<int>(),2*ghostLev));//[14,10]
472   std::vector< std::pair<int,int> > rangeCoarse(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(dimsCoarse));//[(0,14),(0,10)]
473   std::vector<int> fakeFactors(dim,1);
474   //
475   std::vector< std::pair<int,int> > tmp0,tmp1,tmp2;
476   MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(p1,p2,tmp0,false);//tmp0=[(3,4),(1,2)]
477   ApplyFactorsOnCompactFrmt(tmp0,factors);//tmp0=[(12,16),(4,8)]
478   MEDCouplingStructuredMesh::ApplyGhostOnCompactFrmt(tmp0,ghostLev);//tmp0=[(13,17),(5,9)]
479   std::vector< std::pair<int,int> > interstRange(MEDCouplingStructuredMesh::IntersectRanges(tmp0,rangeCoarse));//interstRange=[(13,14),(5,9)]
480   MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(p2,p1,tmp1,false);//tmp1=[(-3,0),(-1,1)]
481   ApplyFactorsOnCompactFrmt(tmp1,factors);//tmp1=[(-12,-4),(-4,0)]
482   MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(tmp1,interstRange,tmp2,false);//tmp2=[(1,2),(1,5)]
483   //
484   std::vector< std::pair<int,int> > dimsFine(p2);
485   ApplyFactorsOnCompactFrmt(dimsFine,factors);
486   ApplyAllGhostOnCompactFrmt(dimsFine,ghostLev);
487   //
488   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ghostVals(MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(dimsFine),dataOnP2,tmp2));
489   MEDCouplingIMesh::CondenseFineToCoarse(dimsCoarse,ghostVals,interstRange,fakeFactors,dataOnP1);
490 }
491
492 MEDCouplingCartesianAMRPatch::MEDCouplingCartesianAMRPatch(const MEDCouplingCartesianAMRPatch& other):MEDCouplingCartesianAMRPatchGen(other),_bl_tr(other._bl_tr)
493 {
494 }
495
496 /*!
497  * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in refined reference.
498  * \param [in] factors - the factors per axis.
499  */
500 void MEDCouplingCartesianAMRPatch::ApplyFactorsOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, const std::vector<int>& factors)
501 {
502   std::size_t sz(factors.size());
503   if(sz!=partBeforeFact.size())
504     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ApplyFactorsOnCompactFrmt : size of input vectors must be the same !");
505   for(std::size_t i=0;i<sz;i++)
506     {
507       partBeforeFact[i].first*=factors[i];
508       partBeforeFact[i].second*=factors[i];
509     }
510 }
511
512 /*!
513  * This method is different than ApplyGhostOnCompactFrmt. The \a partBeforeFact parameter is enlarger contrary to ApplyGhostOnCompactFrmt.
514  *
515  * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in ghost reference.
516  * \param [in] ghostSize - the ghost size of zone for all axis.
517  */
518 void MEDCouplingCartesianAMRPatch::ApplyAllGhostOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, int ghostSize)
519 {
520   if(ghostSize<0)
521     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ApplyAllGhostOnCompactFrmt : ghost size must be >= 0 !");
522   std::size_t sz(partBeforeFact.size());
523   for(std::size_t i=0;i<sz;i++)
524     {
525       partBeforeFact[i].first-=ghostSize;
526       partBeforeFact[i].second+=ghostSize;
527     }
528 }
529
530 MEDCouplingCartesianAMRPatchGF::MEDCouplingCartesianAMRPatchGF(MEDCouplingCartesianAMRMesh *mesh):MEDCouplingCartesianAMRPatchGen(mesh)
531 {
532 }
533
534 MEDCouplingCartesianAMRPatchGF *MEDCouplingCartesianAMRPatchGF::deepCpy() const
535 {
536   return new MEDCouplingCartesianAMRPatchGF(*this);
537 }
538
539 std::size_t MEDCouplingCartesianAMRPatchGF::getHeapMemorySizeWithoutChildren() const
540 {
541   return sizeof(MEDCouplingCartesianAMRPatchGF);
542 }
543
544 MEDCouplingCartesianAMRPatchGF::MEDCouplingCartesianAMRPatchGF(const MEDCouplingCartesianAMRPatchGF& other):MEDCouplingCartesianAMRPatchGen(other)
545 {
546 }
547
548 /// @endcond
549
550 int MEDCouplingCartesianAMRMeshGen::getSpaceDimension() const
551 {
552   return _mesh->getSpaceDimension();
553 }
554
555 void MEDCouplingCartesianAMRMeshGen::setFactors(const std::vector<int>& newFactors)
556 {
557   if(getSpaceDimension()!=(int)newFactors.size())
558     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::setFactors : size of input factors is not equal to the space dimension !");
559   if(_factors.empty())
560     {
561       _factors=newFactors;
562       return ;
563     }
564   if(_factors==newFactors)
565     return ;
566   if(!_patches.empty())
567     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::setFactors : modification of factors is not allowed when presence of patches !");
568   _factors=newFactors;
569   declareAsNew();
570 }
571
572 int MEDCouplingCartesianAMRMeshGen::getMaxNumberOfLevelsRelativeToThis() const
573 {
574   int ret(1);
575   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
576     ret=std::max(ret,(*it)->getMaxNumberOfLevelsRelativeToThis()+1);
577   return ret;
578 }
579
580 /*!
581  * This method returns the number of cells of \a this with the help of the MEDCouplingIMesh instance representing \a this.
582  * The patches in \a this are ignored here.
583  * \sa getNumberOfCellsAtCurrentLevelGhost, getNumberOfCellsRecursiveWithOverlap
584  */
585 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsAtCurrentLevel() const
586 {
587   return _mesh->getNumberOfCells();
588 }
589
590 /*!
591  * 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
592  * to take into account of the ghost cells for future computation.
593  * The patches in \a this are ignored here.
594  *
595  * \sa getNumberOfCellsAtCurrentLevel
596  */
597 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsAtCurrentLevelGhost(int ghostLev) const
598 {
599   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> tmp(_mesh->buildWithGhost(ghostLev));
600   return tmp->getNumberOfCells();
601 }
602
603 /*!
604  * This method returns the number of cells including the current level but \b also \b including recursively all cells of other levels
605  * starting from this. The set of cells which size is returned here are generally overlapping each other.
606  */
607 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsRecursiveWithOverlap() const
608 {
609   int ret(_mesh->getNumberOfCells());
610   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
611     {
612       ret+=(*it)->getNumberOfCellsRecursiveWithOverlap();
613     }
614   return ret;
615 }
616
617 /*!
618  * This method returns the max number of cells covering all the space without overlapping.
619  * It returns the number of cells of the mesh with the highest resolution.
620  * The returned value is equal to the number of cells of mesh returned by buildUnstructured.
621  *
622  * \sa buildUnstructured
623  */
624 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsRecursiveWithoutOverlap() const
625 {
626   int ret(_mesh->getNumberOfCells());
627   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
628     {
629       ret-=(*it)->getNumberOfOverlapedCellsForFather();
630       ret+=(*it)->getNumberOfCellsRecursiveWithoutOverlap();
631     }
632   return ret;
633 }
634
635 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMeshGen::getFather() const
636 {
637   return _father;
638 }
639
640 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMeshGen::getGodFather() const
641 {
642   if(_father==0)
643     return this;
644   else
645     return _father->getGodFather();
646 }
647
648 /*!
649  * This method returns the level of \a this. 0 for god father. 1 for children of god father ...
650  */
651 int MEDCouplingCartesianAMRMeshGen::getAbsoluteLevel() const
652 {
653   if(_father==0)
654     return 0;
655   else
656     return _father->getAbsoluteLevel()+1;
657 }
658
659 int MEDCouplingCartesianAMRMeshGen::getAbsoluteLevelRelativeTo(const MEDCouplingCartesianAMRMeshGen *ref) const
660 {
661   if(this==ref)
662     return 0;
663   if(_father==0)
664     {
665       if(ref==0)
666         return 0;
667       else
668         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getAbsoluteLevelRelativeTo : ref is not in the progeny of this !");
669     }
670   else
671     return _father->getAbsoluteLevelRelativeTo(ref)+1;
672 }
673
674 /*!
675  * This method returns a vector of size equal to getAbsoluteLevelRelativeTo. It allows to find position an absolute position of \a this
676  * relative to \a ref (that is typically the god father).
677  *
678  * \sa getPatchAtPosition
679  */
680 std::vector<int> MEDCouplingCartesianAMRMeshGen::getPositionRelativeTo(const MEDCouplingCartesianAMRMeshGen *ref) const
681 {
682   if(!ref)
683     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPositionRelativeTo : input pointer is NULL !");
684   std::vector<int> ret;
685   getPositionRelativeToInternal(ref,ret);
686   std::reverse(ret.begin(),ret.end());
687   return ret;
688 }
689
690 /*!
691  * \sa getPositionRelativeTo, getMeshAtPosition
692  */
693 const MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRMeshGen::getPatchAtPosition(const std::vector<int>& pos) const
694 {
695   std::size_t sz(pos.size());
696   if(sz==0)
697     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPatchAtPosition : empty input -> no patch by definition !");
698   int patchId(pos[0]);
699   const MEDCouplingCartesianAMRPatch *elt(getPatch(patchId));
700   if(sz==1)
701     return elt;
702   if(!elt || !elt->getMesh())
703     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPatchAtPosition : NULL element found during walk !");
704   std::vector<int> pos2(pos.begin()+1,pos.end());
705   return elt->getMesh()->getPatchAtPosition(pos2);
706 }
707
708 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMeshGen::getMeshAtPosition(const std::vector<int>& pos) const
709 {
710   std::size_t sz(pos.size());
711   if(sz==0)
712     return this;
713   int patchId(pos[0]);
714   const MEDCouplingCartesianAMRPatch *elt(getPatch(patchId));
715   if(sz==1)
716     {
717       if(!elt)
718         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getMeshAtPosition : NULL patch !");
719       return elt->getMesh();
720     }
721   if(!elt || !elt->getMesh())
722     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPatchAtPosition : NULL element found during walk !");
723   std::vector<int> pos2(pos.begin()+1,pos.end());
724   return elt->getMesh()->getMeshAtPosition(pos2);
725 }
726
727 /*!
728  * This method returns grids relative to god father to specified level \a absoluteLev.
729  *
730  * \return std::vector<MEDCouplingCartesianAMRPatchGen *> - objects in vector are to be managed (decrRef) by the caller.
731  */
732 std::vector<MEDCouplingCartesianAMRPatchGen *> MEDCouplingCartesianAMRMeshGen::retrieveGridsAt(int absoluteLev) const
733 {
734   if(absoluteLev<0)
735     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::retrieveGridsAt : absolute level must be >=0 !");
736   if(_father)
737     return getGodFather()->retrieveGridsAt(absoluteLev);
738   //
739   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGen> > rets;
740   retrieveGridsAtInternal(absoluteLev,rets);
741   std::vector< MEDCouplingCartesianAMRPatchGen * > ret(rets.size());
742   for(std::size_t i=0;i<rets.size();i++)
743     {
744       ret[i]=rets[i].retn();
745     }
746   return ret;
747 }
748
749 void MEDCouplingCartesianAMRMeshGen::detachFromFather()
750 {
751   _father=0;
752   declareAsNew();
753 }
754
755 /*!
756  * \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,
757  *                                a the end cell (\b excluded) of the range for the second element of the pair.
758  * \param [in] factors The factor of refinement per axis (different from 0).
759  */
760 void MEDCouplingCartesianAMRMeshGen::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
761 {
762   checkFactorsAndIfNotSetAssign(factors);
763   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> mesh(static_cast<MEDCouplingIMesh *>(_mesh->buildStructuredSubPart(bottomLeftTopRight)));
764   mesh->refineWithFactor(factors);
765   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRMeshSub> zeMesh(new MEDCouplingCartesianAMRMeshSub(this,mesh));
766   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> elt(new MEDCouplingCartesianAMRPatch(zeMesh,bottomLeftTopRight));
767   _patches.push_back(elt);
768   declareAsNew();
769 }
770
771 /// @cond INTERNAL
772
773 class InternalPatch : public RefCountObjectOnly
774 {
775 public:
776   InternalPatch():_nb_of_true(0) { }
777   int getDimension() const { return (int)_part.size(); }
778   double getEfficiency() const { return (double)_nb_of_true/(double)_crit.size(); }
779   int getNumberOfCells() const { return (int)_crit.size(); }
780   void setNumberOfTrue(int nboft) { _nb_of_true=nboft; }
781   std::vector<bool>& getCriterion() { return _crit; }
782   const std::vector<bool>& getConstCriterion() const { return _crit; }
783   void setPart(const std::vector< std::pair<int,int> >& part) { _part=part; }
784   std::vector< std::pair<int,int> >& getPart() { return _part; }
785   const std::vector< std::pair<int,int> >& getConstPart() const { return _part; }
786   bool presenceOfTrue() const { return _nb_of_true>0; }
787   std::vector<int> computeCGS() const { return MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(_part); }
788   std::vector< std::vector<int> > computeSignature() const { return MEDCouplingStructuredMesh::ComputeSignaturePerAxisOf(computeCGS(),getConstCriterion()); }
789   double getEfficiencyPerAxis(int axisId) const { return (double)_nb_of_true/((double)(_part[axisId].second-_part[axisId].first)); }
790   void zipToFitOnCriterion();
791   void updateNumberOfTrue() const;
792   MEDCouplingAutoRefCountObjectPtr<InternalPatch> extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const;
793   MEDCouplingAutoRefCountObjectPtr<InternalPatch> deepCpy() const;
794 protected:
795   ~InternalPatch() { }
796 private:
797   mutable int _nb_of_true;
798   std::vector<bool> _crit;
799   //! _part is global
800   std::vector< std::pair<int,int> > _part;
801 };
802
803 void InternalPatch::zipToFitOnCriterion()
804 {
805   std::vector<int> cgs(computeCGS());
806   std::vector<bool> newCrit;
807   std::vector< std::pair<int,int> > newPart,newPart2;
808   int newNbOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,_crit,newCrit,newPart));
809   MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(_part,newPart,newPart2);
810   if(newNbOfTrue!=_nb_of_true)
811     throw INTERP_KERNEL::Exception("InternalPatch::zipToFitOnCrit : internal error !");
812   _crit=newCrit; _part=newPart2;
813 }
814
815 void InternalPatch::updateNumberOfTrue() const
816 {
817   _nb_of_true=(int)std::count(_crit.begin(),_crit.end(),true);
818 }
819
820 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const
821 {
822   MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
823   std::vector<int> cgs(computeCGS());
824   std::vector< std::pair<int,int> > newPart;
825   MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(_part,partInGlobal,newPart);
826   MEDCouplingStructuredMesh::ExtractFieldOfBoolFrom(cgs,_crit,newPart,ret->getCriterion());
827   ret->setPart(partInGlobal);
828   ret->updateNumberOfTrue();
829   return ret;
830 }
831
832 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::deepCpy() const
833 {
834   MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
835   (*ret)=*this;
836   return ret;
837 }
838
839 void DissectBigPatch(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
840 {
841   cutFound=false; cutPlace=-1;
842   std::vector<double> ratio(rangeOfAxisId-1);
843   for(int id=0;id<rangeOfAxisId-1;id++)
844     {
845       double efficiency[2];
846       for(int h=0;h<2;h++)
847         {
848           std::vector< std::pair<int,int> > rectH(patchToBeSplit->getConstPart());
849           if(h==0)
850             rectH[axisId].second=patchToBeSplit->getConstPart()[axisId].first+id;
851           else
852             rectH[axisId].first=patchToBeSplit->getConstPart()[axisId].first+id;
853           MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(patchToBeSplit->deepCpy());
854           p->zipToFitOnCriterion();
855           //anouar rectH ?
856           efficiency[h]=p->getEfficiencyPerAxis(axisId);
857         }
858       ratio[id]=std::max(efficiency[0],efficiency[1])/std::min(efficiency[0],efficiency[1]);
859     }
860   int minCellDirection(bso.getMinCellDirection()),indexMin(-1);
861   int dimRatioInner(rangeOfAxisId-1-2*(minCellDirection-1));
862   std::vector<double> ratio_inner(dimRatioInner);
863   double minRatio(1.e10);
864   for(int i=0; i<dimRatioInner; i++)
865     {
866       if(ratio[minCellDirection-1+i]<minRatio)
867         {
868           minRatio=ratio[minCellDirection-1+i];
869           indexMin=i+minCellDirection;
870         }
871     }
872   cutFound=true; cutPlace=indexMin+patchToBeSplit->getConstPart()[axisId].first-1;
873 }
874
875 void FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int& axisId, bool& cutFound, int& cutPlace)
876 {
877   cutPlace=-1; cutFound=false;
878   int minCellDirection(bso.getMinCellDirection());
879   const int dim(patchToBeSplit->getDimension());
880   std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
881   for(int id=0;id<dim;id++)
882     {
883       const std::vector<int>& signature(signatures[id]);
884       std::vector<int> hole;
885       std::vector<double> distance;
886       int len((int)signature.size());
887       for(int i=0;i<len;i++)
888         if(signature[i]==0)
889           if(len>= 2*minCellDirection && i >= minCellDirection-1 && i <= len-minCellDirection-1)
890             hole.push_back(i);
891       if(!hole.empty())
892         {
893           double center(((double)len/2.));
894           for(std::size_t i=0;i<hole.size();i++)
895             distance.push_back(fabs(hole[i]+1.-center));
896
897           std::size_t posDistanceMin(std::distance(distance.begin(),std::min_element(distance.begin(),distance.end())));
898           cutFound=true;
899           axisId=id;
900           cutPlace=hole[posDistanceMin]+patchToBeSplit->getConstPart()[axisId].first+1;
901           return ;
902         }
903     }
904 }
905
906 void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, bool& cutFound, int& cutPlace, int& axisId)
907 {
908   cutFound=false; cutPlace=-1;// do not set axisId before to be sure that cutFound was set to true
909
910   const std::vector< std::pair<int,int> >& part(patchToBeSplit->getConstPart());
911   int sign,minCellDirection(bso.getMinCellDirection());
912   const int dim(patchToBeSplit->getDimension());
913
914   std::vector<int> zeroCrossDims(dim,-1);
915   std::vector<int> zeroCrossVals(dim,-1);
916   std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
917   for (int id=0;id<dim;id++)
918     {
919       const std::vector<int>& signature(signatures[id]);
920
921       std::vector<int> derivate_second_order,gradient_absolute,signe_change,zero_cross,edge,max_cross_list ;
922       std::vector<double> distance ;
923
924       for (unsigned int i=1;i<signature.size()-1;i++)
925         derivate_second_order.push_back(signature[i-1]-2*signature[i]+signature[i+1]) ;
926
927       // Gradient absolute value
928       for ( unsigned int i=1;i<derivate_second_order.size();i++)
929         gradient_absolute.push_back(fabs(derivate_second_order[i]-derivate_second_order[i-1])) ;
930       if(derivate_second_order.empty())
931         continue;
932       for (unsigned int i=0;i<derivate_second_order.size()-1;i++)
933         {
934           if (derivate_second_order[i]*derivate_second_order[i+1] < 0 )
935             sign = -1 ;
936           if (derivate_second_order[i]*derivate_second_order[i+1] > 0 )
937             sign = 1 ;
938           if (derivate_second_order[i]*derivate_second_order[i+1] == 0 )
939             sign = 0 ;
940           if ( sign==0 || sign==-1 )
941             if ( i >= (unsigned int)minCellDirection-2 && i <= signature.size()-minCellDirection-2 )
942               {
943                 zero_cross.push_back(i) ;
944                 edge.push_back(gradient_absolute[i]) ;
945               }
946           signe_change.push_back(sign) ;
947         }
948       if ( zero_cross.size() > 0 )
949         {
950           int max_cross=*max_element(edge.begin(),edge.end()) ;
951           for (unsigned int i=0;i<edge.size();i++)
952             if (edge[i]==max_cross)
953               max_cross_list.push_back(zero_cross[i]+1) ;
954
955           double center((signature.size()/2.0));
956           for (unsigned int i=0;i<max_cross_list.size();i++)
957             distance.push_back(fabs(max_cross_list[i]+1-center));
958
959           float distance_min=*min_element(distance.begin(),distance.end()) ;
960           int pos_distance_min=find(distance.begin(),distance.end(),distance_min)-distance.begin();
961           int best_place = max_cross_list[pos_distance_min] + part[id].first ;
962           if ( max_cross >=0 )
963             {
964               zeroCrossDims[id] = best_place ;
965               zeroCrossVals[id] = max_cross ;
966             }
967         }
968       derivate_second_order.clear() ;
969       gradient_absolute.clear() ;
970       signe_change.clear() ;
971       zero_cross.clear() ;
972       edge.clear() ;
973       max_cross_list.clear() ;
974       distance.clear() ;
975     }
976
977   if ( zeroCrossDims[0]!=-1 || zeroCrossDims[1]!=-1  )
978     {
979       int max_cross_dims = *max_element(zeroCrossVals.begin(),zeroCrossVals.end()) ;
980
981       if (zeroCrossVals[0]==max_cross_dims &&  zeroCrossVals[1]==max_cross_dims )
982         {
983           int nl_left(part[0].second-part[0].first);
984           int nc_left(part[1].second-part[1].first);
985           if ( nl_left >=  nc_left )
986             max_cross_dims = 0 ;
987           else
988             max_cross_dims = 1 ;
989         }
990       else
991         max_cross_dims=std::find(zeroCrossVals.begin(),zeroCrossVals.end(),max_cross_dims)-zeroCrossVals.begin();
992       cutFound=true;
993       cutPlace=zeroCrossDims[max_cross_dims];
994       axisId=max_cross_dims ;
995     }
996 }
997
998 void TryAction4(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
999 {
1000   cutFound=false;
1001   if(patchToBeSplit->getEfficiency()<=bso.getEffeciencySnd())
1002     {
1003       if(rangeOfAxisId>=2*bso.getMinCellDirection())
1004         {
1005           cutFound=true;
1006           cutPlace=rangeOfAxisId/2+patchToBeSplit->getConstPart()[axisId].first-1;
1007         }
1008     }
1009   else
1010     {
1011       if(patchToBeSplit->getNumberOfCells()>bso.getMaxCells())
1012         {
1013           DissectBigPatch(bso,patchToBeSplit,axisId,rangeOfAxisId,cutFound,cutPlace);
1014         }
1015     }
1016 }
1017
1018 MEDCouplingAutoRefCountObjectPtr<InternalPatch> DealWithNoCut(const InternalPatch *patch)
1019 {
1020   MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(const_cast<InternalPatch *>(patch));
1021   ret->incrRef();
1022   return ret;
1023 }
1024
1025 void DealWithCut(const InternalPatch *patchToBeSplit, int axisId, int cutPlace, std::vector<MEDCouplingAutoRefCountObjectPtr<InternalPatch> >& listOfPatches)
1026 {
1027   MEDCouplingAutoRefCountObjectPtr<InternalPatch> leftPart,rightPart;
1028   std::vector< std::pair<int,int> > rect(patchToBeSplit->getConstPart());
1029   std::vector< std::pair<int,int> > leftRect(rect),rightRect(rect);
1030   leftRect[axisId].second=cutPlace+1;
1031   rightRect[axisId].first=cutPlace+1;
1032   leftPart=patchToBeSplit->extractPart(leftRect);
1033   rightPart=patchToBeSplit->extractPart(rightRect);
1034   leftPart->zipToFitOnCriterion(); rightPart->zipToFitOnCriterion();
1035   listOfPatches.push_back(leftPart);
1036   listOfPatches.push_back(rightPart);
1037 }
1038
1039 /// @endcond
1040
1041 /*!
1042  * 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.
1043  * This method only create patches at level 0 relative to \a this.
1044  */
1045 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const std::vector<bool>& criterion, const std::vector<int>& factors)
1046 {
1047   int nbCells(getNumberOfCellsAtCurrentLevel());
1048   if(nbCells!=(int)criterion.size())
1049     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 !");
1050   _patches.clear();
1051   std::vector<int> cgs(_mesh->getCellGridStructure());
1052   std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatches,listOfPatchesOK;
1053   //
1054   MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(new InternalPatch);
1055   p->setNumberOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,criterion,p->getCriterion(),p->getPart()));
1056   if(p->presenceOfTrue())
1057     listOfPatches.push_back(p);
1058   while(!listOfPatches.empty())
1059     {
1060       std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatchesTmp;
1061       for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::iterator it=listOfPatches.begin();it!=listOfPatches.end();it++)
1062         {
1063           //
1064           int axisId,rangeOfAxisId,cutPlace;
1065           bool cutFound;
1066           MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt((*it)->getConstPart(),axisId,rangeOfAxisId);
1067           if((*it)->getEfficiency()>=bso.getEffeciency() && (*it)->getNumberOfCells()<bso.getMaxCells())
1068             { listOfPatchesOK.push_back(DealWithNoCut(*it)); continue; }//action 1
1069           FindHole(bso,*it,axisId,cutFound,cutPlace);//axisId overwritten here if FindHole equal to true !
1070           if(cutFound)
1071             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 2
1072           FindInflection(bso,*it,cutFound,cutPlace,axisId);//axisId overwritten here if cutFound equal to true !
1073           if(cutFound)
1074             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 3
1075           TryAction4(bso,*it,axisId,rangeOfAxisId,cutFound,cutPlace);
1076           if(cutFound)
1077             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 4
1078           listOfPatchesOK.push_back(DealWithNoCut(*it));
1079         }
1080       listOfPatches=listOfPatchesTmp;
1081     }
1082   for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::const_iterator it=listOfPatchesOK.begin();it!=listOfPatchesOK.end();it++)
1083     addPatch((*it)->getConstPart(),factors);
1084   declareAsNew();
1085 }
1086
1087 /*!
1088  * 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.
1089  * This method only create patches at level 0 relative to \a this.
1090  */
1091 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const DataArrayByte *criterion, const std::vector<int>& factors)
1092 {
1093   if(!criterion || !criterion->isAllocated())
1094     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterion : the criterion DataArrayByte instance must be allocated and not NULL !");
1095   std::vector<bool> crit(criterion->toVectorOfBool());//check that criterion has one component.
1096   createPatchesFromCriterion(bso,crit,factors);
1097   declareAsNew();
1098 }
1099
1100 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const DataArrayDouble *criterion, const std::vector<int>& factors, double eps)
1101 {
1102   if(!criterion)
1103     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion : null criterion pointer !");
1104   std::vector<bool> inp(criterion->toVectorOfBool(eps));
1105   createPatchesFromCriterion(bso,inp,factors);
1106 }
1107
1108 /*!
1109  * This method creates a multi level patches split at once.
1110  * This method calls as times as size of \a bso createPatchesFromCriterion. Size of \a bso and size of \a factors must be the same !
1111  * \b WARNING, after the call the number of levels in \a this is equal to bso.size() + 1 !
1112  *
1113  * \param [in] bso
1114  * \param [in] criterion
1115  * \param [in] factors
1116  * \param [in] eps - See DataArrayDouble::toVectorOfBool for more information about the semantic of eps.
1117  *
1118  * \sa createPatchesFromCriterion
1119  */
1120 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterionML(const std::vector<const INTERP_KERNEL::BoxSplittingOptions *>& bso, const DataArrayDouble *criterion, const std::vector< std::vector<int> >& factors, double eps)
1121 {
1122   std::size_t nbOfLevs(bso.size());
1123   if(nbOfLevs!=factors.size())
1124     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterionML : size of vectors must be the same !");
1125   if(nbOfLevs==0)
1126     return ;
1127   if(!bso[0])
1128     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterionML : pointers in 1st arg must be not NULL !");
1129   createPatchesFromCriterion(*bso[0],criterion,factors[0],eps);
1130   for(std::size_t i=1;i<nbOfLevs;i++)
1131     {
1132       if(!bso[i])
1133         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterionML : presence of a NULL BoxSplittingOptions in input vector !");
1134       //
1135       std::vector<MEDCouplingCartesianAMRPatchGen *> elts(retrieveGridsAt((int)(i)));
1136       std::size_t sz(elts.size());
1137       std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGen> > elts2(sz);
1138       std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> > elts3(sz);
1139       for(std::size_t ii=0;ii<sz;ii++)
1140         elts2[ii]=elts[ii];
1141       //
1142       static const char TMP_STR[]="TMP";
1143       std::vector< std::pair<std::string,int> > fieldNames(1); fieldNames[0].first=TMP_STR; fieldNames[0].second=1;
1144       MEDCouplingAutoRefCountObjectPtr<MEDCouplingAMRAttribute> att(MEDCouplingAMRAttribute::New(this,fieldNames,0));
1145       att->alloc();
1146       DataArrayDouble *tmpDa(const_cast<DataArrayDouble *>(att->getFieldOn(this,TMP_STR)));
1147       tmpDa->cpyFrom(*criterion);
1148       att->synchronizeCoarseToFine();
1149       for(std::size_t ii=0;ii<sz;ii++)
1150         {
1151           const DataArrayDouble *critOnLeaf(att->getFieldOn(const_cast<MEDCouplingCartesianAMRMeshGen *>(elts[ii]->getMesh()),TMP_STR));
1152           elts3[ii]=const_cast<DataArrayDouble *>(critOnLeaf); elts3[ii]->incrRef();
1153         }
1154       att=0;
1155       for(std::size_t ii=0;ii<sz;ii++)
1156         const_cast<MEDCouplingCartesianAMRMeshGen *>(elts[ii]->getMesh())->createPatchesFromCriterion(*bso[i],elts3[ii],factors[i],eps);
1157     }
1158 }
1159
1160 void MEDCouplingCartesianAMRMeshGen::removeAllPatches()
1161 {
1162   _patches.clear();
1163   declareAsNew();
1164 }
1165
1166 void MEDCouplingCartesianAMRMeshGen::removePatch(int patchId)
1167 {
1168   checkPatchId(patchId);
1169   int sz((int)_patches.size()),j(0);
1170   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> > patches(sz-1);
1171   for(int i=0;i<sz;i++)
1172     if(i!=patchId)
1173       patches[j++]=_patches[i];
1174   (const_cast<MEDCouplingCartesianAMRMeshGen *>(_patches[patchId]->getMesh()))->detachFromFather();
1175   _patches=patches;
1176   declareAsNew();
1177 }
1178
1179 int MEDCouplingCartesianAMRMeshGen::getNumberOfPatches() const
1180 {
1181   return (int)_patches.size();
1182 }
1183
1184 int MEDCouplingCartesianAMRMeshGen::getPatchIdFromChildMesh(const MEDCouplingCartesianAMRMeshGen *mesh) const
1185 {
1186   int ret(0);
1187   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ret++)
1188     {
1189       if((*it)->getMesh()==mesh)
1190         return ret;
1191     }
1192   throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPatchIdFromChildMesh : no such a mesh in my direct progeny !");
1193 }
1194
1195 std::vector< const MEDCouplingCartesianAMRPatch *> MEDCouplingCartesianAMRMeshGen::getPatches() const
1196 {
1197   std::size_t sz(_patches.size());
1198   std::vector< const MEDCouplingCartesianAMRPatch *> ret(sz);
1199   for(std::size_t i=0;i<sz;i++)
1200     ret[i]=_patches[i];
1201   return ret;
1202 }
1203
1204 const MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRMeshGen::getPatch(int patchId) const
1205 {
1206   checkPatchId(patchId);
1207   return _patches[patchId];
1208 }
1209
1210 /*!
1211  * This method states if patch2 (with id \a patchId2) is in the neighborhood of patch1 (with id \a patchId1).
1212  * The neighborhood size is defined by \a ghostLev in the reference of \a this ( \b not in the reference of patches !).
1213  */
1214 bool MEDCouplingCartesianAMRMeshGen::isPatchInNeighborhoodOf(int patchId1, int patchId2, int ghostLev) const
1215 {
1216   const MEDCouplingCartesianAMRPatch *p1(getPatch(patchId1)),*p2(getPatch(patchId2));
1217   return p1->isInMyNeighborhood(p2,GetGhostLevelInFineRef(ghostLev,_factors));
1218 }
1219
1220 /*!
1221  * 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.
1222  * 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
1223  * defined by the patch with id \a patchId.
1224  *
1225  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1226  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1227  * \return DataArrayDouble * - The array of the cell field on the requested patch
1228  *
1229  * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
1230  * \throw if \a cellFieldOnThis is NULL or not allocated
1231  * \sa fillCellFieldOnPatch, MEDCouplingIMesh::SpreadCoarseToFine
1232  */
1233 DataArrayDouble *MEDCouplingCartesianAMRMeshGen::createCellFieldOnPatch(int patchId, const DataArrayDouble *cellFieldOnThis) const
1234 {
1235   if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
1236     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
1237   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1238   const MEDCouplingIMesh *fine(patch->getMesh()->getImageMesh());
1239   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(fine->getNumberOfCells(),cellFieldOnThis->getNumberOfComponents());
1240   ret->copyStringInfoFrom(*cellFieldOnThis);
1241   MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),ret,patch->getBLTRRange(),getFactors());
1242   return ret.retn();
1243 }
1244
1245 /*!
1246  * This method is equivalent to MEDCouplingCartesianAMRMesh::createCellFieldOnPatch except that here instead of creating a new instance
1247  * it fills the value into the \a cellFieldOnPatch data.
1248  *
1249  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1250  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1251  * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
1252  *
1253  * \sa createCellFieldOnPatch, fillCellFieldComingFromPatch
1254  */
1255 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatch(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, bool isConservative) const
1256 {
1257   if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
1258     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
1259   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1260   MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors());
1261   if(isConservative)
1262     {
1263       int fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(getFactors()));
1264       std::transform(cellFieldOnPatch->begin(),cellFieldOnPatch->end(),cellFieldOnPatch->getPointer(),std::bind2nd(std::multiplies<double>(),1./((double)fact)));
1265     }
1266 }
1267
1268 /*!
1269  * This method is the generalization of fillCellFieldOnPatch method. This method only projects coarse to fine without considering the
1270  * potential neighbor patches covered by the ghost cells of patch with id \a patchId.
1271  *
1272  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1273  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1274  * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
1275  * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
1276  *
1277  * \sa fillCellFieldOnPatch, fillCellFieldOnPatchGhostAdv
1278  */
1279 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchGhost(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, int ghostLev, bool isConservative) const
1280 {
1281   if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
1282     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatchGhost : the input cell field array is NULL or not allocated !");
1283   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1284   MEDCouplingIMesh::SpreadCoarseToFineGhost(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),ghostLev);
1285   if(isConservative)
1286     {
1287       int fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(getFactors()));
1288       std::transform(cellFieldOnPatch->begin(),cellFieldOnPatch->end(),cellFieldOnPatch->getPointer(),std::bind2nd(std::multiplies<double>(),1./((double)fact)));
1289     }
1290 }
1291
1292 /*!
1293  * This method is equivalent to  fillCellFieldOnPatchGhost except that here \b ONLY \b the \b ghost \b zone will be updated
1294  * in \a cellFieldOnPatch.
1295  *
1296  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1297  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1298  * \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.
1299  * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
1300  */
1301 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyOnGhostZone(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, int ghostLev) const
1302 {
1303   if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
1304     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchOnlyOnGhostZone : the input cell field array is NULL or not allocated !");
1305   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1306   MEDCouplingIMesh::SpreadCoarseToFineGhostZone(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),ghostLev);
1307 }
1308
1309 /*!
1310  * This method is a refinement of fillCellFieldOnPatchGhost. fillCellFieldOnPatchGhost is first called.
1311  * Then for all other patches than those pointed by \a patchId that overlap the ghost zone of the patch impact the ghost zone adequately.
1312  *
1313  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1314  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
1315  * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
1316  * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
1317  * \param [in] arrsOnPatches - \b WARNING arrsOnPatches[patchId] is \b NOT \b const. All others are const.
1318  *
1319  * \sa fillCellFieldOnPatchOnlyGhostAdv
1320  */
1321 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchGhostAdv(int patchId, const DataArrayDouble *cellFieldOnThis, int ghostLev, const std::vector<const DataArrayDouble *>& arrsOnPatches, bool isConservative) const
1322 {
1323   int nbp(getNumberOfPatches());
1324   if(nbp!=(int)arrsOnPatches.size())
1325     {
1326       std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchGhostAdv : there are " << nbp << " patches in this and " << arrsOnPatches.size() << " arrays in the last parameter !";
1327       throw INTERP_KERNEL::Exception(oss.str().c_str());
1328     }
1329   DataArrayDouble *theFieldToFill(const_cast<DataArrayDouble *>(arrsOnPatches[patchId]));
1330   // first, do as usual
1331   fillCellFieldOnPatchGhost(patchId,cellFieldOnThis,theFieldToFill,ghostLev,isConservative);
1332   fillCellFieldOnPatchOnlyGhostAdv(patchId,ghostLev,arrsOnPatches);
1333 }
1334
1335 /*!
1336  * This method updates the patch with id \a patchId considering the only the all the patches in \a this to fill ghost zone.
1337  * So \b warning, the DataArrayDouble instance \a arrsOnPatches[patchId] is non const.
1338  *
1339  * \sa getPatchIdsInTheNeighborhoodOf
1340  */
1341 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyGhostAdv(int patchId, int ghostLev, const std::vector<const DataArrayDouble *>& arrsOnPatches) const
1342 {
1343   int nbp(getNumberOfPatches());
1344   if(nbp!=(int)arrsOnPatches.size())
1345     {
1346       std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchOnlyGhostAdv : there are " << nbp << " patches in this and " << arrsOnPatches.size() << " arrays in the last parameter !";
1347       throw INTERP_KERNEL::Exception(oss.str().c_str());
1348     }
1349   const MEDCouplingCartesianAMRPatch *refP(getPatch(patchId));
1350   DataArrayDouble *theFieldToFill(const_cast<DataArrayDouble *>(arrsOnPatches[patchId]));
1351   std::vector<int> ids(getPatchIdsInTheNeighborhoodOf(patchId,ghostLev));
1352   for(std::vector<int>::const_iterator it=ids.begin();it!=ids.end();it++)
1353     {
1354       const MEDCouplingCartesianAMRPatch *otherP(getPatch(*it));
1355       MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwo(ghostLev,_factors,refP,otherP,theFieldToFill,arrsOnPatches[*it]);
1356     }
1357 }
1358
1359 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyOnGhostZoneWith(int ghostLev, const MEDCouplingCartesianAMRPatch *patchToBeModified, const MEDCouplingCartesianAMRPatch *neighborPatch, DataArrayDouble *cellFieldOnPatch, const DataArrayDouble *cellFieldNeighbor) const
1360 {
1361   MEDCouplingCartesianAMRPatch::UpdateNeighborsOfOneWithTwo(ghostLev,_factors,patchToBeModified,neighborPatch,cellFieldOnPatch,cellFieldNeighbor);
1362 }
1363
1364 /*!
1365  * This method updates \a cellFieldOnThis part of values coming from the cell field \a cellFieldOnPatch lying on patch having id \a patchId.
1366  *
1367  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1368  * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
1369  * \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.
1370  * \param [in] isConservative - true if the field needs to be conserved. false if maximum principle has to be applied.
1371  *
1372  * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
1373  * \throw if \a cellFieldOnPatch is NULL or not allocated
1374  * \sa createCellFieldOnPatch, MEDCouplingIMesh::CondenseFineToCoarse,fillCellFieldComingFromPatchGhost
1375  */
1376 void MEDCouplingCartesianAMRMeshGen::fillCellFieldComingFromPatch(int patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis, bool isConservative) const
1377 {
1378   if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
1379       throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch : the input cell field array is NULL or not allocated !");
1380   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1381   MEDCouplingIMesh::CondenseFineToCoarse(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis);
1382   if(!isConservative)
1383     {
1384       int fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(getFactors()));
1385       MEDCouplingStructuredMesh::MultiplyPartOf(_mesh->getCellGridStructure(),patch->getBLTRRange(),1./((double)fact),cellFieldOnThis);
1386     }
1387 }
1388
1389 /*!
1390  * This method is the extension of MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch managing the ghost cells. If this
1391  * method is called with \a ghostLev equal to 0 it behaves exactly as MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch.
1392  *
1393  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
1394  * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
1395  * \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.
1396  * \param [in] ghostLev The size of ghost zone (must be >= 0 !)
1397  * \param [in] isConservative - true if the field needs to be conserved. false if maximum principle has to be applied.
1398  *
1399  * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
1400  * \throw if \a cellFieldOnPatch is NULL or not allocated
1401  * \sa fillCellFieldComingFromPatch
1402  */
1403 void MEDCouplingCartesianAMRMeshGen::fillCellFieldComingFromPatchGhost(int patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis, int ghostLev, bool isConservative) const
1404 {
1405   if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
1406     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatchGhost : the input cell field array is NULL or not allocated !");
1407   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1408   MEDCouplingIMesh::CondenseFineToCoarseGhost(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis,ghostLev);
1409   if(!isConservative)
1410     {
1411       int fact(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(getFactors()));
1412       MEDCouplingStructuredMesh::MultiplyPartOfByGhost(_mesh->getCellGridStructure(),patch->getBLTRRange(),ghostLev,1./((double)fact),cellFieldOnThis);
1413     }
1414 }
1415
1416 /*!
1417  * This method finds all patches (located by their ids) that are in the neighborhood of patch with id \a patchId. The neighborhood size is
1418  * defined by ghostLev.
1419  *
1420  * \param [in] patchId - the id of the considered patch.
1421  * \param [in] ghostLev - the size of the neighborhood.
1422  * \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.
1423  */
1424 DataArrayInt *MEDCouplingCartesianAMRMeshGen::findPatchesInTheNeighborhoodOf(int patchId, int ghostLev) const
1425 {
1426   int nbp(getNumberOfPatches());
1427   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
1428   for(int i=0;i<nbp;i++)
1429     {
1430       if(i!=patchId)
1431         if(isPatchInNeighborhoodOf(i,patchId,ghostLev))
1432           ret->pushBackSilent(i);
1433     }
1434   return ret.retn();
1435 }
1436
1437 MEDCouplingUMesh *MEDCouplingCartesianAMRMeshGen::buildUnstructured() const
1438 {
1439   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> part(_mesh->buildUnstructured());
1440   std::vector<bool> bs(_mesh->getNumberOfCells(),false);
1441   std::vector<int> cgs(_mesh->getCellGridStructure());
1442   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > msSafe(_patches.size()+1);
1443   std::size_t ii(0);
1444   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
1445     {
1446       MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
1447       msSafe[ii+1]=(*it)->getMesh()->buildUnstructured();
1448     }
1449   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsOff(DataArrayInt::BuildListOfSwitchedOff(bs));
1450   msSafe[0]=static_cast<MEDCouplingUMesh *>(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
1451   std::vector< const MEDCouplingUMesh * > ms(msSafe.size());
1452   for(std::size_t i=0;i<msSafe.size();i++)
1453     ms[i]=msSafe[i];
1454   return MEDCouplingUMesh::MergeUMeshes(ms);
1455 }
1456
1457 /*!
1458  * This method returns a mesh containing as cells that there is patches at the current level.
1459  * The patches are seen like 'boxes' that is too say the refinement will not appear here.
1460  *
1461  * \return MEDCoupling1SGTUMesh * - A new object to be managed by the caller containing as cells as there are patches in \a this.
1462  */
1463 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMeshGen::buildMeshFromPatchEnvelop() const
1464 {
1465   std::vector<const MEDCoupling1SGTUMesh *> cells;
1466   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> > cellsSafe;
1467   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1468     {
1469       const MEDCouplingCartesianAMRPatch *patch(*it);
1470       if(patch)
1471         {
1472           MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> cell(patch->getMesh()->getImageMesh()->asSingleCell());
1473           MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> cell1SGT(cell->build1SGTUnstructured());
1474           cellsSafe.push_back(cell1SGT); cells.push_back(cell1SGT);
1475         }
1476     }
1477   return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(cells);
1478 }
1479
1480 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMeshGen::buildMeshOfDirectChildrenOnly() const
1481 {
1482   std::vector<const MEDCoupling1SGTUMesh *> patches;
1483   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> > patchesSafe;
1484   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1485       {
1486         const MEDCouplingCartesianAMRPatch *patch(*it);
1487         if(patch)
1488           {
1489             MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> patchMesh(patch->getMesh()->getImageMesh()->build1SGTUnstructured());
1490             patchesSafe.push_back(patchMesh); patches.push_back(patchMesh);
1491           }
1492       }
1493     return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(patches);
1494 }
1495
1496 /*!
1497  * This method works same as buildUnstructured except that arrays are given in input to build a field on cell in output.
1498  * \return MEDCouplingFieldDouble * - a newly created instance the caller has reponsability to deal with.
1499  * \sa buildUnstructured
1500  */
1501 MEDCouplingFieldDouble *MEDCouplingCartesianAMRMeshGen::buildCellFieldOnRecurseWithoutOverlapWithoutGhost(int ghostSz, const std::vector<const DataArrayDouble *>& recurseArrs) const
1502 {
1503   if(recurseArrs.empty())
1504     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::buildCellFieldOnRecurseWithoutOverlapWithoutGhost : array is empty ! Should never happen !");
1505   //
1506   std::vector<bool> bs(_mesh->getNumberOfCells(),false);
1507   std::vector<int> cgs(_mesh->getCellGridStructure());
1508   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> > msSafe(_patches.size()+1);
1509   std::size_t ii(0);
1510   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
1511     {
1512       MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
1513       std::vector<const DataArrayDouble *> tmpArrs(extractSubTreeFromGlobalFlatten((*it)->getMesh(),recurseArrs));
1514       msSafe[ii+1]=(*it)->getMesh()->buildCellFieldOnRecurseWithoutOverlapWithoutGhost(ghostSz,tmpArrs);
1515     }
1516   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsOff(DataArrayInt::BuildListOfSwitchedOff(bs));
1517   //
1518   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret(MEDCouplingFieldDouble::New(ON_CELLS));
1519   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr2(extractGhostFrom(ghostSz,recurseArrs[0]));
1520   arr2=arr2->selectByTupleIdSafe(eltsOff->begin(),eltsOff->end());
1521   ret->setArray(arr2);
1522   ret->setName(arr2->getName());
1523   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> part(_mesh->buildUnstructured());
1524   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> mesh(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
1525   ret->setMesh(mesh);
1526   msSafe[0]=ret;
1527   //
1528   std::vector< const MEDCouplingFieldDouble * > ms(msSafe.size());
1529   for(std::size_t i=0;i<msSafe.size();i++)
1530     ms[i]=msSafe[i];
1531   //
1532   return MEDCouplingFieldDouble::MergeFields(ms);
1533 }
1534
1535 /*!
1536  * This method extracts from \arr arr the part inside \a arr by cutting the \a ghostSz external part.
1537  * \arr is expected to be an array having a number of tuples equal to \c getImageMesh()->buildWithGhost(ghostSz).
1538  */
1539 DataArrayDouble *MEDCouplingCartesianAMRMeshGen::extractGhostFrom(int ghostSz, const DataArrayDouble *arr) const
1540 {
1541   std::vector<int> st(_mesh->getCellGridStructure());
1542   std::vector< std::pair<int,int> > p(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(st));
1543   std::transform(st.begin(),st.end(),st.begin(),std::bind2nd(std::plus<int>(),2*ghostSz));
1544   MEDCouplingStructuredMesh::ApplyGhostOnCompactFrmt(p,ghostSz);
1545   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(st,arr,p));
1546   return ret.retn();
1547 }
1548
1549 /*!
1550  * This method returns all the patches in \a this not equal to \a patchId that are in neighborhood of patch with id \a patchId.
1551  *
1552  * \sa fillCellFieldOnPatchOnlyGhostAdv
1553  */
1554 std::vector<int> MEDCouplingCartesianAMRMeshGen::getPatchIdsInTheNeighborhoodOf(int patchId, int ghostLev) const
1555 {
1556   std::vector<int> ret;
1557   int nbp(getNumberOfPatches());
1558   //
1559   for(int i=0;i<nbp;i++)
1560     {
1561       if(i!=patchId)
1562         if(isPatchInNeighborhoodOf(i,patchId,ghostLev))
1563           ret.push_back(i);
1564     }
1565   return ret;
1566 }
1567
1568 /*!
1569  * This method returns a dump python of \a this. It is useful for users of createPatchesFromCriterion method for debugging.
1570  *
1571  * \sa dumpPatchesOf, createPatchesFromCriterion, createPatchesFromCriterionML
1572  */
1573 std::string MEDCouplingCartesianAMRMeshGen::buildPythonDumpOfThis() const
1574 {
1575   std::ostringstream oss;
1576   oss << "amr=MEDCouplingCartesianAMRMesh(\""<< getImageMesh()->getName() << "\"," << getSpaceDimension() << ",[";
1577   std::vector<int> ngs(getImageMesh()->getNodeGridStructure());
1578   std::vector<double> orig(getImageMesh()->getOrigin()),dxyz(getImageMesh()->getDXYZ());
1579   std::copy(ngs.begin(),ngs.end(),std::ostream_iterator<int>(oss,","));
1580   oss <<  "],[";
1581   std::copy(orig.begin(),orig.end(),std::ostream_iterator<double>(oss,","));
1582   oss << "],[";
1583   std::copy(dxyz.begin(),dxyz.end(),std::ostream_iterator<double>(oss,","));
1584   oss << "])\n";
1585   dumpPatchesOf("amr",oss);
1586   return oss.str();
1587 }
1588
1589 MEDCouplingCartesianAMRMeshGen::MEDCouplingCartesianAMRMeshGen(const MEDCouplingCartesianAMRMeshGen& other, MEDCouplingCartesianAMRMeshGen *father):_father(father),_mesh(other._mesh),_patches(other._patches),_factors(other._factors)
1590 {
1591   const MEDCouplingIMesh *mesh(other._mesh);
1592   if(mesh)
1593     _mesh=static_cast<MEDCouplingIMesh *>(mesh->deepCpy());
1594   std::size_t sz(other._patches.size());
1595   for(std::size_t i=0;i<sz;i++)
1596     {
1597       const MEDCouplingCartesianAMRPatch *patch(other._patches[i]);
1598       if(patch)
1599         _patches[i]=patch->deepCpy();
1600     }
1601 }
1602
1603 MEDCouplingCartesianAMRMeshGen::MEDCouplingCartesianAMRMeshGen(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
1604                                                                const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):_father(0)
1605 {
1606   _mesh=MEDCouplingIMesh::New(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
1607 }
1608
1609 MEDCouplingCartesianAMRMeshGen::MEDCouplingCartesianAMRMeshGen(MEDCouplingCartesianAMRMeshGen *father, MEDCouplingIMesh *mesh):_father(father)
1610 {
1611   if(!_father)
1612     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen(MEDCouplingIMesh *mesh) constructor : empty father !");
1613   if(!mesh)
1614     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen(MEDCouplingIMesh *mesh) constructor : The input mesh is null !");
1615   mesh->checkCoherency();
1616   _mesh=mesh; _mesh->incrRef();
1617 }
1618
1619 void MEDCouplingCartesianAMRMeshGen::checkPatchId(int patchId) const
1620 {
1621   int sz(getNumberOfPatches());
1622   if(patchId<0 || patchId>=sz)
1623     {
1624       std::ostringstream oss; oss << "MEDCouplingCartesianAMRMeshGen::checkPatchId : invalid patchId (" << patchId << ") ! Must be in [0," << sz << ") !";
1625       throw INTERP_KERNEL::Exception(oss.str().c_str());
1626     }
1627 }
1628
1629 void MEDCouplingCartesianAMRMeshGen::checkFactorsAndIfNotSetAssign(const std::vector<int>& factors)
1630 {
1631   if(getSpaceDimension()!=(int)factors.size())
1632     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::checkFactorsAndIfNotSetAssign : invalid size of factors ! size must be equal to the spaceDimension !");
1633   if(_factors.empty())
1634     {
1635       _factors=factors;
1636     }
1637   else
1638     {
1639       if(_factors!=factors)
1640         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::checkFactorsAndIfNotSetAssign : the factors ");
1641     }
1642 }
1643
1644 void MEDCouplingCartesianAMRMeshGen::retrieveGridsAtInternal(int lev, std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGen> >& grids) const
1645 {
1646   if(lev==0)
1647     {
1648       const MEDCouplingCartesianAMRMesh *thisc(dynamic_cast<const MEDCouplingCartesianAMRMesh *>(this));//tony
1649       MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGF> elt(new MEDCouplingCartesianAMRPatchGF(const_cast<MEDCouplingCartesianAMRMesh *>(thisc)));
1650       grids.push_back(DynamicCastSafe<MEDCouplingCartesianAMRPatchGF,MEDCouplingCartesianAMRPatchGen>(elt));
1651     }
1652   else if(lev==1)
1653     {
1654       for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1655         {
1656           const MEDCouplingCartesianAMRPatch *pt(*it);
1657           if(pt)
1658             {
1659               MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> tmp1(*it);
1660               grids.push_back(DynamicCastSafe<MEDCouplingCartesianAMRPatch,MEDCouplingCartesianAMRPatchGen>(tmp1));
1661             }
1662         }
1663     }
1664   else
1665     {
1666       for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1667         {
1668           const MEDCouplingCartesianAMRPatch *pt(*it);
1669           if(pt)
1670             pt->getMesh()->retrieveGridsAtInternal(lev-1,grids);
1671         }
1672     }
1673 }
1674
1675 int MEDCouplingCartesianAMRMeshGen::GetGhostLevelInFineRef(int ghostLev, const std::vector<int>& factors)
1676 {
1677   if(ghostLev<0)
1678     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::GetGhostLevelInFineRef : the ghost size must be >=0 !");
1679   if(factors.empty())
1680     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::GetGhostLevelInFineRef : no factors defined !");
1681   int ghostLevInPatchRef;
1682   if(ghostLev==0)
1683     ghostLevInPatchRef=0;
1684   else
1685     {
1686       ghostLevInPatchRef=(ghostLev-1)/factors[0]+1;
1687       for(std::size_t i=0;i<factors.size();i++)
1688         ghostLevInPatchRef=std::max(ghostLevInPatchRef,(ghostLev-1)/factors[i]+1);
1689     }
1690   return ghostLevInPatchRef;
1691 }
1692
1693 /*!
1694  * This method returns a sub set of \a all. The subset is defined by the \a head in the tree defined by \a this.
1695  * Elements in \a all are expected to be sorted from god father to most refined structure.
1696  */
1697 std::vector<const DataArrayDouble *> MEDCouplingCartesianAMRMeshGen::extractSubTreeFromGlobalFlatten(const MEDCouplingCartesianAMRMeshGen *head, const std::vector<const DataArrayDouble *>& all) const
1698 {
1699   int maxLev(getMaxNumberOfLevelsRelativeToThis());
1700   std::vector<const DataArrayDouble *> ret;
1701   std::vector<const MEDCouplingCartesianAMRMeshGen *> meshes(1,this);
1702   std::size_t kk(0);
1703   for(int i=0;i<maxLev;i++)
1704     {
1705       std::vector<const MEDCouplingCartesianAMRMeshGen *> meshesTmp;
1706       for(std::vector<const MEDCouplingCartesianAMRMeshGen *>::const_iterator it=meshes.begin();it!=meshes.end();it++)
1707         {
1708           if((*it)==head || head->isObjectInTheProgeny(*it))
1709             ret.push_back(all[kk]);
1710           kk++;
1711           std::vector< const MEDCouplingCartesianAMRPatch *> ps((*it)->getPatches());
1712           for(std::vector< const MEDCouplingCartesianAMRPatch *>::const_iterator it0=ps.begin();it0!=ps.end();it0++)
1713             {
1714               const MEDCouplingCartesianAMRMeshGen *mesh((*it0)->getMesh());
1715               meshesTmp.push_back(mesh);
1716             }
1717         }
1718       meshes=meshesTmp;
1719     }
1720   if(kk!=all.size())
1721     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::extractSubTreeFromGlobalFlatten : the size of input vector is not compatible with number of leaves in this !");
1722   return ret;
1723 }
1724
1725 void MEDCouplingCartesianAMRMeshGen::dumpPatchesOf(const std::string& varName, std::ostream& oss) const
1726 {
1727   std::size_t j(0);
1728   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1729     {
1730       const MEDCouplingCartesianAMRPatch *patch(*it);
1731       if(patch)
1732         {
1733           std::ostringstream oss2; oss2 << varName << ".addPatch([";
1734           const std::vector< std::pair<int,int> >& bltr(patch->getBLTRRange());
1735           std::size_t sz(bltr.size());
1736           for(std::size_t i=0;i<sz;i++)
1737             {
1738               oss2 << "(" << bltr[i].first << "," << bltr[i].second << ")";
1739               if(i!=sz-1)
1740                 oss2 << ",";
1741             }
1742           oss2 << "],[";
1743           std::copy(_factors.begin(),_factors.end(),std::ostream_iterator<int>(oss2,","));
1744           oss2 << "])\n";
1745           oss << oss2.str();
1746           std::ostringstream oss3; oss3 << varName << "[" << j++ << "]";
1747           patch->getMesh()->dumpPatchesOf(oss3.str(),oss);
1748         }
1749     }
1750 }
1751
1752 /*!
1753  * \sa getPositionRelativeTo
1754  */
1755 void MEDCouplingCartesianAMRMeshGen::getPositionRelativeToInternal(const MEDCouplingCartesianAMRMeshGen *ref, std::vector<int>& ret) const
1756 {
1757   if(this==ref)
1758     return ;
1759   if(!_father)
1760     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPositionRelativeToInternal : ref is not in the progeny of this !");
1761   int myId(_father->getPatchIdFromChildMesh(this));
1762   ret.push_back(myId);
1763   _father->getPositionRelativeToInternal(ref,ret);
1764 }
1765
1766 std::size_t MEDCouplingCartesianAMRMeshGen::getHeapMemorySizeWithoutChildren() const
1767 {
1768   return sizeof(MEDCouplingCartesianAMRMeshGen);
1769 }
1770
1771 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMeshGen::getDirectChildren() const
1772 {
1773   std::vector<const BigMemoryObject *> ret;
1774   if((const MEDCouplingIMesh *)_mesh)
1775     ret.push_back((const MEDCouplingIMesh *)_mesh);
1776   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1777     {
1778       if((const MEDCouplingCartesianAMRPatch*)*it)
1779         ret.push_back((const MEDCouplingCartesianAMRPatch*)*it);
1780     }
1781   return ret;
1782 }
1783
1784 void MEDCouplingCartesianAMRMeshGen::updateTime() const
1785 {
1786   if((const MEDCouplingIMesh *)_mesh)
1787     updateTimeWith(*_mesh);
1788   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1789     {
1790       const MEDCouplingCartesianAMRPatch *elt(*it);
1791       if(!elt)
1792         continue;
1793       const MEDCouplingCartesianAMRMeshGen *mesh(elt->getMesh());
1794       if(mesh)
1795         updateTimeWith(*mesh);
1796     }
1797 }
1798
1799 MEDCouplingCartesianAMRMeshSub::MEDCouplingCartesianAMRMeshSub(MEDCouplingCartesianAMRMeshGen *father, MEDCouplingIMesh *mesh):MEDCouplingCartesianAMRMeshGen(father,mesh)
1800 {
1801 }
1802
1803 MEDCouplingCartesianAMRMeshSub::MEDCouplingCartesianAMRMeshSub(const MEDCouplingCartesianAMRMeshSub& other, MEDCouplingCartesianAMRMeshGen *father):MEDCouplingCartesianAMRMeshGen(other,father)
1804 {
1805 }
1806
1807 MEDCouplingCartesianAMRMeshSub *MEDCouplingCartesianAMRMeshSub::deepCpy() const
1808 {
1809   return new MEDCouplingCartesianAMRMeshSub(*this,_father);
1810 }
1811
1812 MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::New(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
1813                                                               const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop)
1814 {
1815   return new MEDCouplingCartesianAMRMesh(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
1816 }
1817
1818 MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::deepCpy() const
1819 {
1820   return new MEDCouplingCartesianAMRMesh(*this);
1821 }
1822
1823 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(const MEDCouplingCartesianAMRMesh& other):MEDCouplingCartesianAMRMeshGen(other,0)
1824 {
1825 }
1826
1827 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
1828                                                          const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):MEDCouplingCartesianAMRMeshGen(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop)
1829 {
1830 }
1831
1832 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMesh::getDirectChildren() const
1833 {
1834   std::vector<const BigMemoryObject *> ret(MEDCouplingCartesianAMRMeshGen::getDirectChildren());
1835   return ret;
1836 }
1837
1838 MEDCouplingCartesianAMRMesh::~MEDCouplingCartesianAMRMesh()
1839 {
1840 }