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