From: ageay <ageay>
Date: Tue, 27 Jul 2010 05:53:21 +0000 (+0000)
Subject: Addition of utilities for profile RW in MEDLoader.
X-Git-Tag: V5_1_main_FINAL~73
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=refs%2Fheads%2FBR_HexaSplitTest;p=tools%2Fmedcoupling.git

Addition of utilities for profile RW in MEDLoader.
---

diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx
index 430ca9bee..84ade1a95 100644
--- a/src/MEDCoupling/MEDCouplingMemArray.cxx
+++ b/src/MEDCoupling/MEDCouplingMemArray.cxx
@@ -101,6 +101,34 @@ DataArrayInt *DataArrayDouble::convertToIntArr() const
   return ret;
 }
 
+/*!
+ * This methods has a similar behaviour than std::string::substr. This method returns a newly created DataArrayInt that is part of this with same number of components.
+ * The intervall is specified by [tupleIdBg,tupleIdEnd) except if tupleIdEnd ==-1 in this case the [tupleIdBg,this->end()) will be kept.
+ * This method check that interval is valid regarding this, if not an exception will be thrown.
+ */
+DataArrayDouble *DataArrayDouble::substr(int tupleIdBg, int tupleIdEnd) const throw(INTERP_KERNEL::Exception)
+{
+  int nbt=getNumberOfTuples();
+  if(tupleIdBg<0)
+    throw INTERP_KERNEL::Exception("DataArrayInt::substr : The tupleIdBg parameter must be greater than 0 !");
+  if(tupleIdBg>=nbt)
+    throw INTERP_KERNEL::Exception("DataArrayInt::substr : The tupleIdBg parameter is greater or equal than number of tuples !");
+  int trueEnd=tupleIdEnd;
+  if(tupleIdEnd!=-1)
+    {
+      if(tupleIdEnd>nbt)
+        throw INTERP_KERNEL::Exception("DataArrayInt::substr : The tupleIdBg parameter is greater or equal than number of tuples !");
+    }
+  else
+    trueEnd=nbt;
+  int nbComp=getNumberOfComponents();
+  DataArrayDouble *ret=DataArrayDouble::New();
+  ret->alloc(trueEnd-tupleIdBg,nbComp);
+  ret->copyStringInfoFrom(*this);
+  std::copy(getConstPointer()+tupleIdBg*nbComp,getConstPointer()+trueEnd*nbComp,ret->getPointer());
+  return ret;
+}
+
 void DataArrayDouble::setArrayIn(DataArrayDouble *newArray, DataArrayDouble* &arrayToSet)
 {
   if(newArray!=arrayToSet)
@@ -359,6 +387,34 @@ DataArrayDouble *DataArrayInt::convertToDblArr() const
   return ret;
 }
 
+/*!
+ * This methods has a similar behaviour than std::string::substr. This method returns a newly created DataArrayInt that is part of this with same number of components.
+ * The intervall is specified by [tupleIdBg,tupleIdEnd) except if tupleIdEnd ==-1 in this case the [tupleIdBg,this->end()) will be kept.
+ * This method check that interval is valid regarding this, if not an exception will be thrown.
+ */
+DataArrayInt *DataArrayInt::substr(int tupleIdBg, int tupleIdEnd) const throw(INTERP_KERNEL::Exception)
+{
+  int nbt=getNumberOfTuples();
+  if(tupleIdBg<0)
+    throw INTERP_KERNEL::Exception("DataArrayInt::substr : The tupleIdBg parameter must be greater than 0 !");
+  if(tupleIdBg>=nbt)
+    throw INTERP_KERNEL::Exception("DataArrayInt::substr : The tupleIdBg parameter is greater or equal than number of tuples !");
+  int trueEnd=tupleIdEnd;
+  if(tupleIdEnd!=-1)
+    {
+      if(tupleIdEnd>nbt)
+        throw INTERP_KERNEL::Exception("DataArrayInt::substr : The tupleIdBg parameter is greater or equal than number of tuples !");
+    }
+  else
+    trueEnd=nbt;
+  int nbComp=getNumberOfComponents();
+  DataArrayInt *ret=DataArrayInt::New();
+  ret->alloc(trueEnd-tupleIdBg,nbComp);
+  ret->copyStringInfoFrom(*this);
+  std::copy(getConstPointer()+tupleIdBg*nbComp,getConstPointer()+trueEnd*nbComp,ret->getPointer());
+  return ret;
+}
+
 void DataArrayInt::reAlloc(int nbOfTuples)
 {
   _mem.reAlloc(_info_on_compo.size()*nbOfTuples);
diff --git a/src/MEDCoupling/MEDCouplingMemArray.hxx b/src/MEDCoupling/MEDCouplingMemArray.hxx
index 9974aa563..eec90bf5e 100644
--- a/src/MEDCoupling/MEDCouplingMemArray.hxx
+++ b/src/MEDCoupling/MEDCouplingMemArray.hxx
@@ -113,6 +113,7 @@ namespace ParaMEDMEM
     //!alloc or useArray should have been called before.
     MEDCOUPLING_EXPORT void reAlloc(int nbOfTuples);
     MEDCOUPLING_EXPORT DataArrayInt *convertToIntArr() const;
+    MEDCOUPLING_EXPORT DataArrayDouble *substr(int tupleIdBg, int tupleIdEnd=-1) const throw(INTERP_KERNEL::Exception);
     MEDCOUPLING_EXPORT void getTuple(int tupleId, double *res) const { std::copy(_mem.getConstPointerLoc(tupleId*_info_on_compo.size()),_mem.getConstPointerLoc((tupleId+1)*_info_on_compo.size()),res); }
     MEDCOUPLING_EXPORT double getIJ(int tupleId, int compoId) const { return _mem[tupleId*_info_on_compo.size()+compoId]; }
     MEDCOUPLING_EXPORT void setIJ(int tupleId, int compoId, double newVal) { _mem[tupleId*_info_on_compo.size()+compoId]=newVal; }
@@ -150,6 +151,7 @@ namespace ParaMEDMEM
     //!alloc or useArray should have been called before.
     MEDCOUPLING_EXPORT void reAlloc(int nbOfTuples);
     MEDCOUPLING_EXPORT DataArrayDouble *convertToDblArr() const;
+    MEDCOUPLING_EXPORT DataArrayInt *substr(int tupleIdBg, int tupleIdEnd=-1) const throw(INTERP_KERNEL::Exception);
     MEDCOUPLING_EXPORT void getTuple(int tupleId, int *res) const { std::copy(_mem.getConstPointerLoc(tupleId*_info_on_compo.size()),_mem.getConstPointerLoc((tupleId+1)*_info_on_compo.size()),res); }
     MEDCOUPLING_EXPORT int getIJ(int tupleId, int compoId) const { return _mem[tupleId*_info_on_compo.size()+compoId]; }
     MEDCOUPLING_EXPORT void setIJ(int tupleId, int compoId, int newVal) { _mem[tupleId*_info_on_compo.size()+compoId]=newVal; }
diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx
index 342fc33a8..26f538942 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -1837,7 +1837,7 @@ DataArrayInt *MEDCouplingUMesh::rearrange2ConsecutiveCellTypes()
 /*!
  * This methods split this into as mush as untructured meshes that consecutive set of same type cells.
  * So this method has typically a sense if MEDCouplingUMesh::checkConsecutiveCellTypes has a sense.
- * This method makes asumption (no check) that connectivity is correctly set before calling.
+ * This method makes asumption that connectivity is correctly set before calling.
  */
 std::vector<MEDCouplingUMesh *> MEDCouplingUMesh::splitByType() const
 {
@@ -1863,6 +1863,67 @@ std::vector<MEDCouplingUMesh *> MEDCouplingUMesh::splitByType() const
   return ret;
 }
 
+/*!
+ * This method makes the assumption that da->getNumberOfTuples()<this->getNumberOfCells(). This method makes the assumption that ids contained in 'da'
+ * are in [0:getNumberOfCells())
+ */
+DataArrayInt *MEDCouplingUMesh::convertCellArrayPerGeoType(const DataArrayInt *da) const throw(INTERP_KERNEL::Exception)
+{
+  checkFullyDefined();
+  const int *conn=_nodal_connec->getConstPointer();
+  const int *connI=_nodal_connec_index->getConstPointer();
+  int nbOfCells=getNumberOfCells();
+  std::set<INTERP_KERNEL::NormalizedCellType> types=getAllTypes();
+  int *tmp=new int[nbOfCells];
+  for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator iter=types.begin();iter!=types.end();iter++)
+    {
+      int j=0;
+      for(const int *i=connI;i!=connI+nbOfCells;i++)
+        if((INTERP_KERNEL::NormalizedCellType)conn[*i]==(*iter))
+          tmp[std::distance(connI,i)]=j++;
+    }
+  DataArrayInt *ret=DataArrayInt::New();
+  ret->alloc(da->getNumberOfTuples(),da->getNumberOfComponents());
+  ret->copyStringInfoFrom(*da);
+  int *retPtr=ret->getPointer();
+  const int *daPtr=da->getConstPointer();
+  int nbOfElems=da->getNbOfElems();
+  for(int k=0;k<nbOfElems;k++)
+    retPtr[k]=tmp[daPtr[k]];
+  delete [] tmp;
+  return ret;
+}
+
+/*!
+ * This method reduced number of cells of this by keeping cells whose type is different from 'type' and if type=='type'
+ * cells whose ids is in 'idsPerGeoType' array.
+ * This method conserves coords and name of mesh.
+ */
+MEDCouplingUMesh *MEDCouplingUMesh::keepSpecifiedCells(INTERP_KERNEL::NormalizedCellType type, const std::vector<int>& idsPerGeoType)
+{
+  std::vector<int> idsTokeep;
+  int nbOfCells=getNumberOfCells();
+  int j=0;
+  for(int i=0;i<nbOfCells;i++)
+    if(getTypeOfCell(i)!=type)
+      idsTokeep.push_back(i);
+    else
+      {
+        if(std::find(idsPerGeoType.begin(),idsPerGeoType.end(),j)!=idsPerGeoType.end())
+          idsTokeep.push_back(i);
+        j++;
+      }
+  MEDCouplingPointSet *ret=buildPartOfMySelf(&idsTokeep[0],&idsTokeep[0]+idsTokeep.size(),true);
+  MEDCouplingUMesh *ret2=dynamic_cast<MEDCouplingUMesh *>(ret);
+  if(!ret2)
+    {
+      ret->decrRef();
+      return 0;
+    }
+  ret2->setName(getName());
+  return ret2;
+}
+
 /*!
  * Returns a newly created mesh (with ref count ==1) that contains merge of 'this' and 'other'.
  */
diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx
index 3dde3bde2..f3fdf9a77 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.hxx
@@ -96,6 +96,8 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT bool checkConsecutiveCellTypes() const;
     MEDCOUPLING_EXPORT DataArrayInt *rearrange2ConsecutiveCellTypes();
     MEDCOUPLING_EXPORT std::vector<MEDCouplingUMesh *> splitByType() const;
+    MEDCOUPLING_EXPORT DataArrayInt *convertCellArrayPerGeoType(const DataArrayInt *da) const throw(INTERP_KERNEL::Exception);
+    MEDCOUPLING_EXPORT MEDCouplingUMesh *keepSpecifiedCells(INTERP_KERNEL::NormalizedCellType type, const std::vector<int>& idsPerGeoType);
     //
     MEDCOUPLING_EXPORT MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const;
     MEDCOUPLING_EXPORT DataArrayDouble *getBarycenterAndOwner() const;
diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
index 7307e3834..503eba3c0 100644
--- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
+++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
@@ -34,6 +34,7 @@ namespace ParaMEDMEM
     CPPUNIT_TEST_SUITE(MEDCouplingBasicsTest);
     CPPUNIT_TEST( testArray );
     CPPUNIT_TEST( testArray2 );
+    CPPUNIT_TEST( testArray3 );
     CPPUNIT_TEST( testMesh );
     CPPUNIT_TEST( testMeshPointsCloud );
     CPPUNIT_TEST( testMeshM1D );
@@ -144,6 +145,7 @@ namespace ParaMEDMEM
   public:
     void testArray();
     void testArray2();
+    void testArray3();
     void testMesh();
     void testMeshPointsCloud();
     void testMeshM1D();
diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx
index 46e79f9ec..3ce905d5e 100644
--- a/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx
+++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx
@@ -78,6 +78,49 @@ void MEDCouplingBasicsTest::testArray2()
   arr->decrRef();
 }
 
+void MEDCouplingBasicsTest::testArray3()
+{
+  DataArrayInt *arr1=DataArrayInt::New();
+  arr1->alloc(7,2);
+  int *tmp=arr1->getPointer();
+  const int arr1Ref[14]={0,10,1,11,2,12,3,13,4,14,5,15,6,16};
+  std::copy(arr1Ref,arr1Ref+14,tmp);
+  CPPUNIT_ASSERT_EQUAL(7,arr1->getNumberOfTuples());
+  CPPUNIT_ASSERT_EQUAL(2,arr1->getNumberOfComponents());
+  CPPUNIT_ASSERT(std::equal(arr1Ref,arr1Ref+14,arr1->getConstPointer()));
+  DataArrayInt *arr2=arr1->substr(3);
+  CPPUNIT_ASSERT_EQUAL(4,arr2->getNumberOfTuples());
+  CPPUNIT_ASSERT_EQUAL(2,arr2->getNumberOfComponents());
+  CPPUNIT_ASSERT(std::equal(arr1Ref+6,arr1Ref+14,arr2->getConstPointer()));
+  arr2->decrRef();
+  DataArrayInt *arr3=arr1->substr(2,5);
+  CPPUNIT_ASSERT_EQUAL(3,arr3->getNumberOfTuples());
+  CPPUNIT_ASSERT_EQUAL(2,arr3->getNumberOfComponents());
+  CPPUNIT_ASSERT(std::equal(arr1Ref+4,arr1Ref+10,arr3->getConstPointer()));
+  arr1->decrRef();
+  arr3->decrRef();
+  //
+  DataArrayDouble *arr4=DataArrayDouble::New();
+  arr4->alloc(7,2);
+  double *tmp2=arr4->getPointer();
+  const int arr4Ref[14]={0.8,10.8,1.9,11.9,2.1,12.1,3.2,13.2,4.3,14.3,5.4,15.4,6.5,16.5};
+  std::copy(arr4Ref,arr4Ref+14,tmp2);
+  CPPUNIT_ASSERT_EQUAL(7,arr4->getNumberOfTuples());
+  CPPUNIT_ASSERT_EQUAL(2,arr4->getNumberOfComponents());
+  CPPUNIT_ASSERT(std::equal(arr4Ref,arr4Ref+14,arr4->getConstPointer()));
+  DataArrayDouble *arr5=arr4->substr(3);
+  CPPUNIT_ASSERT_EQUAL(4,arr5->getNumberOfTuples());
+  CPPUNIT_ASSERT_EQUAL(2,arr5->getNumberOfComponents());
+  CPPUNIT_ASSERT(std::equal(arr4Ref+6,arr4Ref+14,arr5->getConstPointer()));
+  arr5->decrRef();
+  DataArrayDouble *arr6=arr4->substr(2,5);
+  CPPUNIT_ASSERT_EQUAL(3,arr6->getNumberOfTuples());
+  CPPUNIT_ASSERT_EQUAL(2,arr6->getNumberOfComponents());
+  CPPUNIT_ASSERT(std::equal(arr4Ref+4,arr4Ref+10,arr6->getConstPointer()));
+  arr4->decrRef();
+  arr6->decrRef();
+}
+
 void MEDCouplingBasicsTest::testMesh()
 {
   const int nbOfCells=6;
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
index f6254d4fd..8fef761ba 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
@@ -35,6 +35,48 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         arr3=arr2.convertToDblArr();
         self.failUnless(arr.isEqual(arr3,1e-14))
         pass
+
+    def testArray3(self):
+        arr1=DataArrayInt.New();
+        arr1Ref=[0,10,1,11,2,12,3,13,4,14,5,15,6,16]
+        arr1.setValues(arr1Ref,7,2);
+        self.failUnlessEqual(7,arr1.getNumberOfTuples());
+        self.failUnlessEqual(2,arr1.getNumberOfComponents());
+        self.failUnlessEqual(arr1Ref,arr1.getValues());
+        arr2=arr1.substr(3);
+        self.failUnlessEqual(4,arr2.getNumberOfTuples());
+        self.failUnlessEqual(2,arr2.getNumberOfComponents());
+        self.failUnlessEqual(arr1Ref[6:],arr2.getValues());
+        arr3=arr1.substr(2,5);
+        self.failUnlessEqual(3,arr3.getNumberOfTuples());
+        self.failUnlessEqual(2,arr3.getNumberOfComponents());
+        self.failUnlessEqual(arr1Ref[4:10],arr3.getValues());
+        #
+        arr4=DataArrayDouble.New();
+        arr4Ref=[0.8,10.8,1.9,11.9,2.1,12.1,3.2,13.2,4.3,14.3,5.4,15.4,6.5,16.5]
+        arr4.setValues(arr4Ref,7,2);
+        self.failUnlessEqual(7,arr4.getNumberOfTuples());
+        self.failUnlessEqual(2,arr4.getNumberOfComponents());
+        tmp=arr4.getValues()
+        for i in xrange(14):
+            self.failUnless(abs(arr4Ref[i]-tmp[i])<1e-14);
+            pass
+        arr5=arr4.substr(3);
+        self.failUnlessEqual(4,arr5.getNumberOfTuples());
+        self.failUnlessEqual(2,arr5.getNumberOfComponents());
+        tmp=arr5.getValues()
+        for i in xrange(8):
+            self.failUnless(abs(arr4Ref[6+i]-tmp[i])<1e-14);
+            pass
+        arr6=arr4.substr(2,5);
+        self.failUnlessEqual(3,arr6.getNumberOfTuples());
+        self.failUnlessEqual(2,arr6.getNumberOfComponents());
+        tmp=arr6.getValues()
+        for i in xrange(6):
+            self.failUnless(abs(arr4Ref[4+i]-tmp[i])<1e-14);
+            pass
+        pass
+
     def testMesh(self):
         tab4=[1, 2, 8, 7, 2, 3, 9, 8, 3,
               4, 10, 9, 4, 5, 11, 10, 5,
diff --git a/src/MEDCoupling_Swig/MEDCouplingTypemaps.i b/src/MEDCoupling_Swig/MEDCouplingTypemaps.i
index 6583ae968..5bc13876d 100644
--- a/src/MEDCoupling_Swig/MEDCouplingTypemaps.i
+++ b/src/MEDCoupling_Swig/MEDCouplingTypemaps.i
@@ -76,6 +76,34 @@ static int *convertPyToNewIntArr2(PyObject *pyLi, int *size)
     }
 }
 
+static void convertPyToNewIntArr3(PyObject *pyLi, std::vector<int>& arr)
+{
+  if(PyList_Check(pyLi))
+    {
+      int size=PyList_Size(pyLi);
+      arr.resize(size);
+      for(int i=0;i<size;i++)
+        {
+          PyObject *o=PyList_GetItem(pyLi,i);
+          if(PyInt_Check(o))
+            {
+              int val=(int)PyInt_AS_LONG(o);
+              arr[i]=val;
+            }
+          else
+            {
+              PyErr_SetString(PyExc_TypeError,"list must contain integers only");
+              PyErr_Print();
+            }
+        }
+    }
+  else
+    {
+      PyErr_SetString(PyExc_TypeError,"convertPyToNewIntArr : not a list");
+      PyErr_Print();
+    }
+}
+
 static PyObject *convertDblArrToPyList(const double *ptr, int size)
 {
   PyObject *ret=PyList_New(size);
diff --git a/src/MEDCoupling_Swig/libMEDCoupling_Swig.i b/src/MEDCoupling_Swig/libMEDCoupling_Swig.i
index e5a317ce0..2f20bf5e5 100644
--- a/src/MEDCoupling_Swig/libMEDCoupling_Swig.i
+++ b/src/MEDCoupling_Swig/libMEDCoupling_Swig.i
@@ -65,11 +65,13 @@ using namespace INTERP_KERNEL;
 %newobject ParaMEDMEM::DataArrayDouble::performCpy;
 %newobject ParaMEDMEM::DataArrayInt::deepCopy;
 %newobject ParaMEDMEM::DataArrayInt::performCpy;
+%newobject ParaMEDMEM::DataArrayInt::substr;
 %newobject ParaMEDMEM::DataArrayDouble::aggregate;
 %newobject ParaMEDMEM::DataArrayDouble::add;
 %newobject ParaMEDMEM::DataArrayDouble::substract;
 %newobject ParaMEDMEM::DataArrayDouble::multiply;
 %newobject ParaMEDMEM::DataArrayDouble::divide;
+%newobject ParaMEDMEM::DataArrayDouble::substr;
 %newobject ParaMEDMEM::MEDCouplingFieldDouble::clone;
 %newobject ParaMEDMEM::MEDCouplingFieldDouble::buildNewTimeReprFromThis;
 %newobject ParaMEDMEM::MEDCouplingMesh::buildOrthogonalField;
@@ -84,6 +86,7 @@ using namespace INTERP_KERNEL;
 %newobject ParaMEDMEM::MEDCouplingUMesh::mergeUMeshes;
 %newobject ParaMEDMEM::MEDCouplingUMesh::buildNewNumberingFromCommNodesFrmt;
 %newobject ParaMEDMEM::MEDCouplingUMesh::rearrange2ConsecutiveCellTypes;
+%newobject ParaMEDMEM::MEDCouplingUMesh::convertCellArrayPerGeoType;
 %newobject ParaMEDMEM::MEDCouplingExtrudedMesh::New;
 %newobject ParaMEDMEM::MEDCouplingCMesh::New;
 %feature("unref") DataArrayDouble "$this->decrRef();"
@@ -355,6 +358,7 @@ namespace ParaMEDMEM
     //tools
     bool checkConsecutiveCellTypes() const;
     DataArrayInt *rearrange2ConsecutiveCellTypes();
+    DataArrayInt *convertCellArrayPerGeoType(const DataArrayInt *da) const throw(INTERP_KERNEL::Exception);
     DataArrayInt *zipConnectivityTraducer(int compType);
     void getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const;
     MEDCouplingUMesh *buildDescendingConnectivity(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const;
@@ -416,6 +420,14 @@ namespace ParaMEDMEM
         return ret;
       }
 
+      PyObject *keepSpecifiedCells(INTERP_KERNEL::NormalizedCellType type, PyObject *ids)
+      {
+        std::vector<int> idsPerGeoType;
+        convertPyToNewIntArr3(ids,idsPerGeoType);
+        MEDCouplingUMesh *ret=self->keepSpecifiedCells(type,idsPerGeoType);
+        return SWIG_NewPointerObj(SWIG_as_voidptr(ret),SWIGTYPE_p_ParaMEDMEM__MEDCouplingUMesh, SWIG_POINTER_OWN | 0 );
+      }
+
       PyObject *getCellsContainingPoints(PyObject *p, int nbOfPoints, double eps) const
       {
         int sz;
diff --git a/src/MEDLoader/MEDLoader.cxx b/src/MEDLoader/MEDLoader.cxx
index fa0724c10..be1aa373f 100644
--- a/src/MEDLoader/MEDLoader.cxx
+++ b/src/MEDLoader/MEDLoader.cxx
@@ -109,6 +109,10 @@ med_geometrie_element typmai3[32] = { MED_POINT1,//0
                                       MED_POLYEDRE//31
 };
 
+double MEDLoader::_EPS_FOR_NODE_COMP=1.e-12;
+
+int MEDLoader::_COMP_FOR_CELL=0;
+
 using namespace ParaMEDMEM;
 
 namespace MEDLoaderNS
@@ -161,8 +165,11 @@ namespace MEDLoaderNS
                                                   INTERP_KERNEL::NormalizedCellType type, std::vector<int>& conn4MEDFile, std::vector<int>& fam4MEDFile);
   ParaMEDMEM::MEDCouplingFieldDouble *readFieldDoubleLev1(const char *fileName, const char *meshName, int meshDimRelToMax, const char *fieldName, int iteration, int order,
                                                           ParaMEDMEM::TypeOfField typeOfOutField);
+  med_idt appendFieldSimpleAtt(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, med_int& numdt, med_int& numo, med_float& dt);
   void appendFieldDirectly(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f);
-  void prepareCellFieldDoubleForWriting(const ParaMEDMEM::MEDCouplingFieldDouble *f, std::list<MEDLoader::MEDFieldDoublePerCellType>& split);
+  void appendNodeProfileField(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, const int *thisMeshNodeIds);
+  void appendCellProfileField(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, const int *thisMeshCellIds);
+  void prepareCellFieldDoubleForWriting(const ParaMEDMEM::MEDCouplingFieldDouble *f, const int *cellIds, std::list<MEDLoader::MEDFieldDoublePerCellType>& split);
   void writeUMeshDirectly(const char *fileName, ParaMEDMEM::MEDCouplingUMesh *mesh, const DataArrayInt *families, bool forceFromScratch);
   void writeUMeshesDirectly(const char *fileName, const char *meshName, const std::vector<ParaMEDMEM::MEDCouplingUMesh *>& meshes, bool forceFromScratch);
   void writeFieldAndMeshDirectly(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, bool forceFromScratch);
@@ -171,6 +178,22 @@ namespace MEDLoaderNS
 
 const char WHITE_SPACES[]=" \n";
 
+/*!
+ * This method set the epsilon value used for node comparison when trying to buid a profile for a field on node/cell on an already written mesh.
+ */
+void MEDLoader::setEpsilonForNodeComp(double val)
+{
+  _EPS_FOR_NODE_COMP=val;
+}
+
+/*!
+ * This method set the policy comparison when trying to fit the already written mesh on a field. The semantic of the policy is specified in MEDCouplingUMesh::zipConnectivityTraducer.
+ */
+void MEDLoader::setCompPolicyForCell(int val)
+{
+  _COMP_FOR_CELL=val;
+}
+
 /*!
  * @param lgth is the size of fam tab. For classical types conn is size of 'lgth'*number_of_nodes_in_type.
  * @param index is optionnal only for polys. Set it to 0 if it is not the case.
@@ -201,8 +224,10 @@ void MEDLoader::MEDConnOfOneElemType::releaseArray()
   delete [] _global;
 }
 
-MEDLoader::MEDFieldDoublePerCellType::MEDFieldDoublePerCellType(INTERP_KERNEL::NormalizedCellType type, double *values, int ncomp, int ntuple):_ntuple(ntuple),_ncomp(ncomp),_values(values),_type(type)
+MEDLoader::MEDFieldDoublePerCellType::MEDFieldDoublePerCellType(INTERP_KERNEL::NormalizedCellType type, double *values, int ncomp, int ntuple, const int *cellIdPerType):_ntuple(ntuple),_ncomp(ncomp),_values(values),_type(type)
 {
+  if(cellIdPerType)
+    _cell_id_per_type.insert(_cell_id_per_type.end(),cellIdPerType,cellIdPerType+ntuple);
 }
 
 void MEDLoader::MEDFieldDoublePerCellType::releaseArray()
@@ -571,7 +596,14 @@ void MEDLoaderNS::readFieldDoubleDataInMedFile(const char *fileName, const char
                     }
                   MEDchampLire(fid,(char *)meshName,(char *)fieldName,(unsigned char*)valr,MED_FULL_INTERLACE,MED_ALL,locname,
                                pflname,MED_COMPACT,tabEnt[typeOfOutField],tabType[typeOfOutField][j],iteration,order);
-                  field.push_back(MEDLoader::MEDFieldDoublePerCellType(typmai2[j],valr,ncomp,nval));
+                  int *pfl=0;
+                  if(pflname[0]!='\0')
+                    {
+                      pfl=new int[nval];
+                      MEDprofilLire(fid,pfl,pflname);
+                    }
+                  field.push_back(MEDLoader::MEDFieldDoublePerCellType(typmai2[j],valr,ncomp,nval,pfl));
+                  delete [] pfl;
                   delete [] dt_unit;
                   delete [] maa_ass;
                 }
@@ -1208,6 +1240,20 @@ ParaMEDMEM::MEDCouplingFieldDouble *MEDLoaderNS::readFieldDoubleLev1(const char
   ParaMEDMEM::MEDCouplingUMesh *mesh=readUMeshFromFileLev1(fileName,meshName,meshDimRelToMax,familiesToKeep,typesToKeep,meshDim);
   if(typeOfOutField==ON_CELLS)
     MEDLoaderNS::keepSpecifiedMeshDim<MEDLoader::MEDFieldDoublePerCellType>(fieldPerCellType,meshDim);
+  //for profiles
+  for(std::list<MEDLoader::MEDFieldDoublePerCellType>::const_iterator iter=fieldPerCellType.begin();iter!=fieldPerCellType.end();iter++)
+    {
+      const std::vector<int>& cellIds=(*iter).getCellIdPerType();
+      if(!cellIds.empty())
+        {
+          std::vector<int> ci(cellIds.size());
+          std::transform(cellIds.begin(),cellIds.end(),ci.begin(),std::bind2nd(std::plus<int>(),-1));
+          ParaMEDMEM::MEDCouplingUMesh *mesh2=mesh->keepSpecifiedCells((*iter).getType(),ci);
+          mesh->decrRef();
+          mesh=mesh2;
+        }
+    }
+  //
   ParaMEDMEM::MEDCouplingFieldDouble *ret=ParaMEDMEM::MEDCouplingFieldDouble::New(typeOfOutField,ONE_TIME);
   ret->setName(fieldName);
   ret->setTime(time,iteration,order);
@@ -1409,15 +1455,68 @@ void MEDLoaderNS::writeUMeshesDirectly(const char *fileName, const char *meshNam
   m->decrRef();
 }
 
-void MEDLoaderNS::appendFieldDirectly(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f)
+/*!
+ * This method makes the assumption that f->getMesh() nodes are fully included in already written mesh in 'fileName'.
+ * @param thisMeshNodeIds points to a tab of size f->getMesh()->getNumberOfNodes() that says for a node i in f->getMesh() that its id is thisMeshNodeIds[i] is already written mesh.
+ */
+void MEDLoaderNS::appendNodeProfileField(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, const int *thisMeshNodeIds)
+{
+  //not implemented yet.
+  med_int numdt,numo;
+  med_float dt;
+  int nbComp=f->getNumberOfComponents();
+  med_idt fid=appendFieldSimpleAtt(fileName,f,numdt,numo,dt);
+  MEDfermer(fid);
+}
+
+/*!
+ * This method makes the assumption that f->getMesh() cells are fully included in already written mesh in 'fileName'.
+ * @param thisMeshCellIdsPerType points to a tab of size f->getMesh()->getNumberOfCells() that says for a cell i in f->getMesh() that its id is thisMeshCellIds[i] of corresponding type is already written mesh.
+ */
+void MEDLoaderNS::appendCellProfileField(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, const int *thisMeshCellIdsPerType)
+{
+  //not implemented yet.
+  med_int numdt,numo;
+  med_float dt;
+  int nbComp=f->getNumberOfComponents();
+  med_idt fid=appendFieldSimpleAtt(fileName,f,numdt,numo,dt);
+  std::list<MEDLoader::MEDFieldDoublePerCellType> split;
+  prepareCellFieldDoubleForWriting(f,thisMeshCellIdsPerType,split);
+  const double *pt=f->getArray()->getConstPointer();
+  int number=0;
+  for(std::list<MEDLoader::MEDFieldDoublePerCellType>::const_iterator iter=split.begin();iter!=split.end();iter++)
+    {
+      char *nommaa=MEDLoaderBase::buildEmptyString(MED_TAILLE_NOM);
+      strcpy(nommaa,f->getMesh()->getName());
+      char *profileName=MEDLoaderBase::buildEmptyString(MED_TAILLE_NOM);
+      std::ostringstream oss; oss << "Pfl" << f->getName() << "_" << number++;
+      strcpy(profileName,oss.str().c_str());
+      const std::vector<int>& ids=(*iter).getCellIdPerType();
+      int *profile=new int [ids.size()];
+      std::transform(ids.begin(),ids.end(),profile,std::bind2nd(std::plus<int>(),1));
+      MEDprofilEcr(fid,profile,ids.size(),profileName);
+      delete [] profile;
+      MEDchampEcr(fid,nommaa,(char *)f->getName(),(unsigned char*)pt,MED_FULL_INTERLACE,(*iter).getNbOfTuple(),
+                  (char *)MED_NOGAUSS,MED_ALL,profileName,MED_COMPACT,MED_MAILLE,
+                  typmai3[(int)(*iter).getType()],numdt,(char *)"",dt,numo);
+      delete [] profileName;
+      delete [] nommaa;
+      pt+=(*iter).getNbOfTuple()*nbComp;
+    }
+  MEDfermer(fid);
+}
+
+/*!
+ * This method performs the classical job for fields before any values setting.
+ */
+med_idt MEDLoaderNS::appendFieldSimpleAtt(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, med_int& numdt, med_int& numo, med_float& dt)
 {
   med_idt fid=MEDouvrir((char *)fileName,MED_LECTURE_ECRITURE);
   int nbComp=f->getNumberOfComponents();
   char *comp=MEDLoaderBase::buildEmptyString(nbComp*MED_TAILLE_PNOM);
   char *unit=MEDLoaderBase::buildEmptyString(nbComp*MED_TAILLE_PNOM);
   MEDchampCr(fid,(char *)f->getName(),MED_FLOAT64,comp,unit,nbComp);
-  med_int numdt,numo;
-  med_float dt;
+  
   ParaMEDMEM::TypeOfTimeDiscretization td=f->getTimeDiscretization();
   if(td==ParaMEDMEM::NO_TIME)
     {
@@ -1430,13 +1529,24 @@ void MEDLoaderNS::appendFieldDirectly(const char *fileName, ParaMEDMEM::MEDCoupl
       numdt=(med_int)tmp1; numo=(med_int)tmp2;
       dt=(med_float)tmp0;
     }
+  delete [] comp;
+  delete [] unit;
+  return fid;
+}
+
+void MEDLoaderNS::appendFieldDirectly(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f)
+{
+  med_int numdt,numo;
+  med_float dt;
+  int nbComp=f->getNumberOfComponents();
+  med_idt fid=appendFieldSimpleAtt(fileName,f,numdt,numo,dt);
   const double *pt=f->getArray()->getConstPointer();
   switch(f->getTypeOfField())
     {
     case ParaMEDMEM::ON_CELLS:
       {
         std::list<MEDLoader::MEDFieldDoublePerCellType> split;
-        prepareCellFieldDoubleForWriting(f,split);
+        prepareCellFieldDoubleForWriting(f,0,split);
         for(std::list<MEDLoader::MEDFieldDoublePerCellType>::const_iterator iter=split.begin();iter!=split.end();iter++)
           {
             char *nommaa=MEDLoaderBase::buildEmptyString(MED_TAILLE_NOM);
@@ -1462,12 +1572,14 @@ void MEDLoaderNS::appendFieldDirectly(const char *fileName, ParaMEDMEM::MEDCoupl
     default:
       throw INTERP_KERNEL::Exception("Not managed this type of FIELD !");
     }
-  delete [] comp;
-  delete [] unit;
   MEDfermer(fid);
 }
 
-void MEDLoaderNS::prepareCellFieldDoubleForWriting(const ParaMEDMEM::MEDCouplingFieldDouble *f, std::list<MEDLoader::MEDFieldDoublePerCellType>& split)
+/*!
+ * This method splits field 'f' into types to be ready for writing.
+ * @param cellIdsPerType this parameter can be 0 if not in profile mode. If it is != 0 this array is of size f->getMesh()->getNumberOfCells().
+ */
+void MEDLoaderNS::prepareCellFieldDoubleForWriting(const ParaMEDMEM::MEDCouplingFieldDouble *f, const int *cellIdsPerType, std::list<MEDLoader::MEDFieldDoublePerCellType>& split)
 {
   int nbComp=f->getNumberOfComponents();
   const MEDCouplingMesh *mesh=f->getMesh();
@@ -1480,11 +1592,18 @@ void MEDLoaderNS::prepareCellFieldDoubleForWriting(const ParaMEDMEM::MEDCoupling
   const int *conn=meshC->getNodalConnectivity()->getConstPointer();
   int nbOfCells=meshC->getNumberOfCells();
   INTERP_KERNEL::NormalizedCellType curType;
+  const int *wCellIdsPT=cellIdsPerType;
   for(const int *pt=connI;pt!=connI+nbOfCells;)
     {
       curType=(INTERP_KERNEL::NormalizedCellType)conn[*pt];
       const int *pt2=std::find_if(pt+1,connI+nbOfCells,ConnReaderML(conn,(int)curType));
-      split.push_back(MEDLoader::MEDFieldDoublePerCellType(curType,0,nbComp,pt2-pt));
+      if(!cellIdsPerType)
+        split.push_back(MEDLoader::MEDFieldDoublePerCellType(curType,0,nbComp,pt2-pt,0));
+      else
+        {
+          split.push_back(MEDLoader::MEDFieldDoublePerCellType(curType,0,nbComp,pt2-pt,wCellIdsPT));
+          wCellIdsPT+=std::distance(pt,pt2);
+        }
       pt=pt2;
     }
 }
@@ -1521,6 +1640,55 @@ void MEDLoaderNS::writeFieldTryingToFitExistingMesh(const char *fileName, ParaME
       throw INTERP_KERNEL::Exception(oss.str().c_str());
     }
   MEDCouplingUMesh *m=MEDLoader::ReadUMeshFromFile(fileName,f->getMesh()->getName(),f2);
+  MEDCouplingUMesh *m2=MEDCouplingUMesh::mergeUMeshes(m,(MEDCouplingUMesh *)f->getMesh());
+  bool areNodesMerged;
+  DataArrayInt *da=m2->mergeNodes(MEDLoader::_EPS_FOR_NODE_COMP,areNodesMerged);
+  if(!areNodesMerged || m2->getNumberOfNodes()!=m->getNumberOfNodes())
+    {
+      da->decrRef();
+      m2->decrRef();
+      m->decrRef();
+      std::ostringstream oss; oss << "Nodes in already written mesh \"" << f->getMesh()->getName() << "\" in file \"" << fileName << "\" does not fit coordinates of unstructured grid f->getMesh() !";
+      throw INTERP_KERNEL::Exception(oss.str().c_str());
+    }
+  switch(f->getTypeOfField())
+    {
+    case ParaMEDMEM::ON_CELLS:
+      {
+        da->decrRef();
+        DataArrayInt *da2=m2->zipConnectivityTraducer(MEDLoader::_COMP_FOR_CELL);
+        if(m2->getNumberOfCells()!=m->getNumberOfCells())
+          {
+            da2->decrRef();
+            m2->decrRef();
+            m->decrRef();
+            std::ostringstream oss1; oss1 << "Cells in already written mesh \"" << f->getMesh()->getName() << "\" in file \"" << fileName << "\" does not fit connectivity of unstructured grid f->getMesh() !";
+            throw INTERP_KERNEL::Exception(oss1.str().c_str());
+          }
+        da=m2->convertCellArrayPerGeoType(da2);
+        DataArrayInt *da3=da->substr(m2->getNumberOfCells());
+        da2->decrRef();
+        da2=m2->convertCellArrayPerGeoType(da3);
+        da3->decrRef();
+        appendCellProfileField(fileName,f,da2->getConstPointer());
+        da2->decrRef();
+        break;
+      }
+    case ParaMEDMEM::ON_NODES:
+      {
+        appendNodeProfileField(fileName,f,da->getConstPointer()+m->getNumberOfNodes());
+        break;
+      }
+    default:
+      {
+        da->decrRef();
+        m2->decrRef();
+        m->decrRef();
+        throw INTERP_KERNEL::Exception("Not implemented other profile fitting from already written mesh for fields than on NODES and on CELLS.");
+      }
+    }
+  da->decrRef();
+  m2->decrRef();
   m->decrRef();
 }
 
diff --git a/src/MEDLoader/MEDLoader.hxx b/src/MEDLoader/MEDLoader.hxx
index 965fa0584..83dc10ecb 100644
--- a/src/MEDLoader/MEDLoader.hxx
+++ b/src/MEDLoader/MEDLoader.hxx
@@ -37,7 +37,7 @@ namespace ParaMEDMEM
 
 class MEDLOADER_EXPORT MEDLoader
 {
-public:
+ public:
   class MEDConnOfOneElemType
   {
   public:
@@ -63,20 +63,24 @@ public:
   class MEDFieldDoublePerCellType
   {
   public:
-    MEDFieldDoublePerCellType(INTERP_KERNEL::NormalizedCellType type, double *values, int ncomp, int ntuple);
+    MEDFieldDoublePerCellType(INTERP_KERNEL::NormalizedCellType type, double *values, int ncomp, int ntuple, const int *cellIdPerType);
     INTERP_KERNEL::NormalizedCellType getType() const { return _type; }
     int getNbComp() const { return _ncomp; }
     int getNbOfTuple() const { return _ntuple; }
     int getNbOfValues() const { return _ncomp*_ntuple; }
     double *getArray() const { return _values; }
+    const std::vector<int>& getCellIdPerType() const { return _cell_id_per_type; }
     void releaseArray();
   private:
     int _ntuple;
     int _ncomp;
     double *_values;
+    std::vector<int> _cell_id_per_type;
     INTERP_KERNEL::NormalizedCellType _type;
   };
   //
+  static void setEpsilonForNodeComp(double val);
+  static void setCompPolicyForCell(int val);
   static std::vector<std::string> GetMeshNames(const char *fileName);
   static std::vector<std::string> GetMeshGroupsNames(const char *fileName, const char *meshName);
   static std::vector<std::string> GetMeshFamilyNames(const char *fileName, const char *meshName);
@@ -97,8 +101,11 @@ public:
   static void WriteUMeshes(const char *fileName, const char *meshName, const std::vector<ParaMEDMEM::MEDCouplingUMesh *>& meshes, bool writeFromScratch);
   static void WriteField(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, bool writeFromScratch);
   static void WriteFieldUsingAlreadyWrittenMesh(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f);
-private:
+ private:
   MEDLoader();
+ public:
+  static double _EPS_FOR_NODE_COMP;
+  static int _COMP_FOR_CELL;
 };
 
 #endif
diff --git a/src/MEDLoader/Test/MEDLoaderTest.cxx b/src/MEDLoader/Test/MEDLoaderTest.cxx
index 1ce22526d..341e52cee 100644
--- a/src/MEDLoader/Test/MEDLoaderTest.cxx
+++ b/src/MEDLoader/Test/MEDLoaderTest.cxx
@@ -331,6 +331,9 @@ void MEDLoaderTest::testFieldProfilRW1()
 {
   const char fileName[]="file12.med";
   MEDCouplingUMesh *mesh1=build3DMesh_1();
+  bool b;
+  DataArrayInt *da=mesh1->mergeNodes(1e-12,b);
+  da->decrRef();
   MEDLoader::WriteUMesh(fileName,mesh1,true);
   const int part1[5]={1,2,4,13,15};
   MEDCouplingUMesh *mesh2=(MEDCouplingUMesh *)mesh1->buildPartOfMySelf(part1,part1+5,true);
@@ -352,6 +355,12 @@ void MEDLoaderTest::testFieldProfilRW1()
   f1->checkCoherency();
   //
   MEDLoader::WriteField(fileName,f1,false);//<- false important for the test
+  //
+  MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldDoubleCell(fileName,f1->getMesh()->getName(),0,f1->getName(),2,7);
+  f2->checkCoherency();
+  CPPUNIT_ASSERT(f1->isEqual(f2,1e-12,1e-12));
+  //
+  f2->decrRef();
   f1->decrRef();
   mesh1->decrRef();
   mesh2->decrRef();