From: geay <anthony.geay@cea.fr>
Date: Mon, 7 Jul 2014 10:11:49 +0000 (+0200)
Subject: VTHB write
X-Git-Tag: V7_5_0a1~2^2~19
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=a0c3d1050d1756ffb7a93f9369993a6a0e7dcd4e;p=tools%2Fmedcoupling.git

VTHB write
---

diff --git a/src/MEDCoupling/MEDCouplingAMRAttribute.cxx b/src/MEDCoupling/MEDCouplingAMRAttribute.cxx
index d78f2a58d..7f8b5298f 100644
--- a/src/MEDCoupling/MEDCouplingAMRAttribute.cxx
+++ b/src/MEDCoupling/MEDCouplingAMRAttribute.cxx
@@ -24,6 +24,7 @@
 #include "MEDCouplingIMesh.hxx"
 
 #include <sstream>
+#include <fstream>
 
 using namespace ParaMEDMEM;
 
@@ -969,7 +970,6 @@ MEDCouplingFieldDouble *MEDCouplingAMRAttribute::buildCellFieldOnWithGhost(MEDCo
  */
 MEDCouplingFieldDouble *MEDCouplingAMRAttribute::buildCellFieldOnWithoutGhost(MEDCouplingCartesianAMRMeshGen *mesh, const std::string& fieldName) const
 {
-  //tony
   const DataArrayDouble *arr(0);
   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingGridCollection> >::const_iterator it=_levs.begin();it!=_levs.end();it++)
     {
@@ -1000,12 +1000,122 @@ MEDCouplingFieldDouble *MEDCouplingAMRAttribute::buildCellFieldOnWithoutGhost(ME
   return ret.retn();
 }
 
-/*!
- * This method is useful just after a remesh after a time step computation to project values in \a this to the new
- * mesh \a targetGF.
- *
- * This method performs a projection from \a this to a target AMR mesh \a targetGF.
- * This method performs the projection by trying to transfer the finest information to \a targetGF.
+
+std::string MEDCouplingAMRAttribute::writeVTHB(const std::string& fileName) const
+{
+  static const char EXT[]=".vthb";
+  std::string baseName,extName,zeFileName;
+  MEDCouplingMesh::SplitExtension(fileName,baseName,extName);
+  if(extName==EXT)
+    zeFileName=fileName;
+  else
+    { zeFileName=baseName; zeFileName+=EXT; }
+  //
+  std::ofstream ofs(fileName.c_str());
+  ofs << "<VTKFile type=\"vtkOverlappingAMR\" version=\"1.1\" byte_order=\"" << MEDCouplingByteOrderStr() << "\">\n";
+  const MEDCouplingCartesianAMRMesh *gf(getMyGodFather());
+  ofs << "  <vtkOverlappingAMR origin=\"";
+  const MEDCouplingIMesh *gfm(gf->getImageMesh());
+  std::vector<double> orig(gfm->getOrigin());
+  std::vector<double> spacing(gfm->getDXYZ());
+  int dim((int)orig.size());
+  std::copy(orig.begin(),orig.end(),std::ostream_iterator<double>(ofs," ")); ofs << "\" grid_description=\"";
+  for(int i=0;i<dim;i++)
+    {
+      char tmp[2]; tmp[0]='X'+i; tmp[1]='\0';
+      ofs << tmp;
+    }
+  ofs << "\">\n";
+  //
+  int maxLev(gf->getMaxNumberOfLevelsRelativeToThis()),kk(0);
+  for(int i=0;i<maxLev;i++)
+    {
+      std::vector<MEDCouplingCartesianAMRPatchGen *> patches(gf->retrieveGridsAt(i));
+      std::size_t sz(patches.size());
+      std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGen> > patchesSafe(sz);
+      for(std::size_t j=0;j<sz;j++)
+        patchesSafe[j]=patches[j];
+      if(sz==0)
+        continue;
+      ofs << "    <Block level=\"" << i << "\" spacing=\"";
+      std::copy(spacing.begin(),spacing.end(),std::ostream_iterator<double>(ofs," "));
+      ofs << "\">\n";
+      if(i!=maxLev-1)
+        {
+          std::vector<int> factors(patches[0]->getMesh()->getFactors());
+          for(int k=0;k<dim;k++)
+            spacing[k]*=1./((double) factors[k]);
+        }
+      std::size_t jj(0);
+      for(std::vector<MEDCouplingCartesianAMRPatchGen *>::const_iterator it=patches.begin();it!=patches.end();it++,jj++,kk++)
+        {
+          ofs << "      <DataSet index=\"" << jj << "\" amr_box=\"";
+          const MEDCouplingCartesianAMRPatch *patchCast(dynamic_cast<const MEDCouplingCartesianAMRPatch *>(*it));
+          const MEDCouplingCartesianAMRMeshGen *mesh((*it)->getMesh());
+          if(patchCast)
+            {
+              const std::vector< std::pair<int,int> >& bltr(patchCast->getBLTRRangeRelativeToGF());
+              for(int pp=0;pp<dim;pp++)
+                ofs << bltr[pp].first << " " << bltr[pp].second-1 << " ";
+            }
+          else
+            {
+              const MEDCouplingIMesh *im((*it)->getMesh()->getImageMesh());
+              std::vector<int> cgs(im->getCellGridStructure());
+              for(int pp=0;pp<dim;pp++)
+                ofs << "0 " << cgs[pp]-1 << " ";
+            }
+          ofs << "\" file=\"";
+          //
+          int tmp(-1);
+          if(_levs[i]->presenceOf((*it)->getMesh(),tmp))
+            {
+              const DataArrayDoubleCollection& ddc(_levs[i]->getFieldsAt(tmp));
+              std::vector<DataArrayDouble *> arrs(ddc.retrieveFields());
+              std::size_t nbFields(arrs.size());
+              std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> > arrsSafe(nbFields),arrs2Safe(nbFields);
+              std::vector< const MEDCouplingFieldDouble *> fields(nbFields);
+              std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> > fieldsSafe(nbFields);
+              for(std::size_t pp=0;pp<nbFields;pp++)
+                arrsSafe[pp]=arrs[pp];
+              for(std::size_t pp=0;pp<nbFields;pp++)
+                {
+                  MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> im(mesh->getImageMesh()->buildWithGhost(_ghost_lev));
+                  std::vector<int> cgs(mesh->getImageMesh()->getCellGridStructure()),cgsWG(im->getCellGridStructure());
+                  arrs2Safe[pp]=DataArrayDouble::New();
+                  arrs2Safe[pp]->alloc(mesh->getImageMesh()->getNumberOfCells(),arrs[pp]->getNumberOfComponents());
+                  std::vector< std::pair<int,int> > cgs2(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(cgs));
+                  MEDCouplingStructuredMesh::ApplyGhostOnCompactFrmt(cgs2,_ghost_lev);
+                  std::vector<int> fakeFactors(mesh->getImageMesh()->getSpaceDimension(),1);
+                  MEDCouplingIMesh::SpreadCoarseToFine(arrs[pp],cgsWG,arrs2Safe[pp],cgs2,fakeFactors);
+                  arrs2Safe[pp]->copyStringInfoFrom(*arrs[pp]);
+                  //
+                  fieldsSafe[pp]=MEDCouplingFieldDouble::New(ON_CELLS); fields[pp]=fieldsSafe[pp];
+                  fieldsSafe[pp]->setMesh(mesh->getImageMesh());
+                  fieldsSafe[pp]->setArray(arrs2Safe[pp]);
+                  fieldsSafe[pp]->setName(arrs[pp]->getName());
+                }
+              std::ostringstream vtiFileName; vtiFileName << baseName << "_" << kk << ".vti";
+              MEDCouplingFieldDouble::WriteVTK(vtiFileName.str(),fields,true);
+              //
+              ofs << vtiFileName.str() << "\">\n";
+              ofs << "      \n      </DataSet>\n";
+            }
+        }
+      ofs << "    </Block>\n";
+    }
+  //
+  ofs << "  </vtkOverlappingAMR>\n";
+  ofs << "</VTKFile>\n";
+  return zeFileName;
+}
+
+  /*!
+   * This method is useful just after a remesh after a time step computation to project values in \a this to the new
+   * mesh \a targetGF.
+   *
+   * This method performs a projection from \a this to a target AMR mesh \a targetGF.
+   * This method performs the projection by trying to transfer the finest information to \a targetGF.
  * \b WARNING this method does not update the ghost zone, if any.
  * The level0 of \a this god father must have the same structure than those of \a targetGF.
  *
diff --git a/src/MEDCoupling/MEDCouplingAMRAttribute.hxx b/src/MEDCoupling/MEDCouplingAMRAttribute.hxx
index 861c0b5ec..37857698a 100644
--- a/src/MEDCoupling/MEDCouplingAMRAttribute.hxx
+++ b/src/MEDCoupling/MEDCouplingAMRAttribute.hxx
@@ -144,6 +144,7 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildCellFieldOnRecurseWithoutOverlapWithoutGhost(MEDCouplingCartesianAMRMeshGen *mesh, const std::string& fieldName) const;
     MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildCellFieldOnWithGhost(MEDCouplingCartesianAMRMeshGen *mesh, const std::string& fieldName) const;
     MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildCellFieldOnWithoutGhost(MEDCouplingCartesianAMRMeshGen *mesh, const std::string& fieldName) const;
+    MEDCOUPLING_EXPORT std::string writeVTHB(const std::string& fileName) const;
     //
     MEDCOUPLING_EXPORT MEDCouplingAMRAttribute *projectTo(MEDCouplingCartesianAMRMesh *targetGF) const;
     //
diff --git a/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx b/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx
index 68b90c750..84a89c516 100644
--- a/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx
@@ -193,6 +193,43 @@ bool MEDCouplingCartesianAMRPatch::isInMyNeighborhoodDiffLev(const MEDCouplingCa
   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
 }
 
+std::vector< std::pair<int,int> > MEDCouplingCartesianAMRPatch::getBLTRRangeRelativeToGF() const
+{
+  std::vector< std::pair<int,int> > ret(_bl_tr);
+  const MEDCouplingCartesianAMRMeshGen *mesh(getMesh());
+  if(!mesh)
+    throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::getBLTRRangeRelativeToGF : not valid !");
+  const MEDCouplingCartesianAMRMeshGen *fath(mesh->getFather());
+  if(!fath)
+    throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::getBLTRRangeRelativeToGF : not valid 2 !");
+  std::vector<int> factors(fath->getFactors());
+  std::size_t sz(ret.size());
+  for(std::size_t ii=0;ii<sz;ii++)
+    {
+      ret[ii].first*=factors[ii];
+      ret[ii].second*=factors[ii];
+    }
+  const MEDCouplingCartesianAMRMeshGen *oldFather(fath);
+  fath=oldFather->getFather();
+  while(fath)
+    {
+      int pos(fath->getPatchIdFromChildMesh(oldFather));
+      const MEDCouplingCartesianAMRPatch *p(fath->getPatch(pos));
+      const std::vector< std::pair<int,int> >& tmp(p->getBLTRRange());
+      const std::vector<int>& factors2(fath->getFactors());
+      std::transform(factors.begin(),factors.end(),factors2.begin(),factors.begin(),std::multiplies<int>());
+      for(std::size_t ii=0;ii<sz;ii++)
+        {
+          int delta(ret[ii].second-ret[ii].first);
+          ret[ii].first+=tmp[ii].first*factors[ii];
+          ret[ii].second=ret[ii].first+delta;
+        }
+      oldFather=fath;
+      fath=oldFather->getFather();
+    }
+  return ret;
+}
+
 std::vector<int> MEDCouplingCartesianAMRPatch::computeCellGridSt() const
 {
   const MEDCouplingCartesianAMRMeshGen *m(getMesh());
diff --git a/src/MEDCoupling/MEDCouplingCartesianAMRMesh.hxx b/src/MEDCoupling/MEDCouplingCartesianAMRMesh.hxx
index a8eecc4e1..bee12ea85 100644
--- a/src/MEDCoupling/MEDCouplingCartesianAMRMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingCartesianAMRMesh.hxx
@@ -82,6 +82,7 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT bool isInMyNeighborhoodDiffLev(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const;
     // basic set/get
     MEDCOUPLING_EXPORT const std::vector< std::pair<int,int> >& getBLTRRange() const { return _bl_tr; }
+    MEDCOUPLING_EXPORT std::vector< std::pair<int,int> > getBLTRRangeRelativeToGF() const;
     MEDCOUPLING_EXPORT std::vector<int> computeCellGridSt() const;
     MEDCOUPLING_EXPORT static bool IsInMyNeighborhood(int ghostLev, const std::vector< std::pair<int,int> >& p1, const std::vector< std::pair<int,int> >& p2);
     //
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
index a81eebcb3..691a760c1 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
@@ -15626,6 +15626,27 @@ class MEDCouplingBasicsTest(unittest.TestCase):
             self.assertTrue(att.getFieldOn(amr[1].getMesh(),"YY").isEqualWithoutConsideringStr(exp1,1e-12))
             pass
         pass
+    
+    def testSwig2AMR14(self):
+        """ non regression linked to VTHB write."""
+        fact=[2,2] ; fact2=[3,3]
+        amr=MEDCouplingCartesianAMRMesh("mesh",2,[5,5],[0.,0.],[1.,1.])
+        amr.addPatch([(1,3),(0,2)],fact)
+        amr.addPatch([(1,3),(3,4)],fact)
+        amr[0].addPatch([(1,3),(1,3)],fact2)
+        amr[1].addPatch([(1,3),(1,2)],fact2)
+        att=MEDCouplingAMRAttribute(amr,[("YY",1)],2)
+        att.alloc()
+        att.getFieldOn(amr,"YY").iota(0.1)
+        att.getFieldOn(amr[0].getMesh(),"YY").iota(0.2)
+        att.getFieldOn(amr[1].getMesh(),"YY").iota(0.3)
+        att.getFieldOn(amr[0][0].getMesh(),"YY").iota(0.4)
+        att.getFieldOn(amr[1][0].getMesh(),"YY").iota(0.5)
+        self.assertEqual(amr[0].getBLTRRangeRelativeToGF(),[(2,6),(0,4)])
+        self.assertEqual(amr[1].getBLTRRangeRelativeToGF(),[(2,6),(6,8)])
+        self.assertEqual(amr[0][0].getBLTRRangeRelativeToGF(),[(9,15),(3,9)])
+        self.assertEqual(amr[1][0].getBLTRRangeRelativeToGF(),[(9,15),(21,24)])
+        pass
 
     def testSwig2Intersect2DMeshWith1DLine1(self):
         """A basic test with no colinearity between m1 and m2."""
diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i
index 4ff7649a6..cdacf6ec4 100644
--- a/src/MEDCoupling_Swig/MEDCouplingCommon.i
+++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i
@@ -4904,6 +4904,12 @@ namespace ParaMEDMEM
         return convertFromVectorPairInt(ret);
       }
 
+      PyObject *getBLTRRangeRelativeToGF() const throw(INTERP_KERNEL::Exception)
+      {
+        std::vector< std::pair<int,int> > ret(self->getBLTRRangeRelativeToGF());
+        return convertFromVectorPairInt(ret);
+      }
+
       void addPatch(PyObject *bottomLeftTopRight, const std::vector<int>& factors) throw(INTERP_KERNEL::Exception)
       {
         std::vector< std::pair<int,int> > inp;
@@ -5222,6 +5228,7 @@ namespace ParaMEDMEM
     MEDCouplingFieldDouble *buildCellFieldOnWithoutGhost(MEDCouplingCartesianAMRMeshGen *mesh, const std::string& fieldName) const throw(INTERP_KERNEL::Exception);
     bool changeGodFather(MEDCouplingCartesianAMRMesh *gf) throw(INTERP_KERNEL::Exception);
     MEDCouplingAMRAttribute *projectTo(MEDCouplingCartesianAMRMesh *targetGF) const throw(INTERP_KERNEL::Exception);
+    std::string writeVTHB(const std::string& fileName) const throw(INTERP_KERNEL::Exception);
     %extend
     {
       static MEDCouplingAMRAttribute *New(MEDCouplingCartesianAMRMesh *gf, PyObject *fieldNames, int ghostLev) throw(INTERP_KERNEL::Exception)