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