]> SALOME platform Git repositories - modules/med.git/commitdiff
Salome HOME
MEDCouplingFieldDouble::nodeToCellDiscretization
authorageay <ageay>
Tue, 17 Sep 2013 08:43:17 +0000 (08:43 +0000)
committerageay <ageay>
Tue, 17 Sep 2013 08:43:17 +0000 (08:43 +0000)
src/MEDCoupling/MEDCouplingField.cxx
src/MEDCoupling/MEDCouplingFieldDouble.cxx
src/MEDCoupling/MEDCouplingFieldDouble.hxx
src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
src/MEDCoupling_Swig/MEDCouplingCommon.i

index 2931cac0a67037e3e9c2908c5e73a03ab0d42fa1..a4a0534f6dc53b633db384fd09553cf6b45a5c37 100644 (file)
@@ -174,10 +174,13 @@ std::vector<const BigMemoryObject *> MEDCouplingField::getDirectChildren() const
 /*!
  * Returns a type of \ref MEDCouplingSpatialDisc "spatial discretization" of \a this
  * field in terms of enum ParaMEDMEM::TypeOfField. 
- *  \return ParaMEDMEM::TypeOfField - the type of \a this field.
+ * \return ParaMEDMEM::TypeOfField - the type of \a this field.
+ * \throw If the geometric type is empty.
  */
 TypeOfField MEDCouplingField::getTypeOfField() const
 {
+  if(!((const MEDCouplingFieldDiscretization *)_type))
+    throw INTERP_KERNEL::Exception("MEDCouplingField::getTypeOfField : spatial discretization is null !");
   return _type->getEnum();
 }
 
index b18199dd6cef18142c0b096afe2ec57897222e8f..e721d1423edb000fda95a7a592f8c3c516c5cd50 100644 (file)
@@ -217,6 +217,62 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildNewTimeReprFromThis(TypeOfT
   return ret.retn();
 }
 
+/*!
+ * This method converts a field on nodes (\a this) to a cell field (returned field). The convertion is a \b non \b conservative remapping !
+ * This method is useful only for users that need a fast convertion from node to cell spatial discretization. The algorithm applied is simply to attach
+ * to each cell the average of values on nodes constituting this cell.
+ *
+ * \return MEDCouplingFieldDouble* - a new instance of MEDCouplingFieldDouble. The
+ *         caller is to delete this field using decrRef() as it is no more needed. The returned field will share the same mesh object object than those in \a this.
+ * \throw If \a this spatial discretization is empty or not ON_NODES.
+ * \throw If \a this is not coherent (see MEDCouplingFieldDouble::checkCoherency).
+ * 
+ * \warning This method is a \b non \b conservative method of remapping from node spatial discretization to cell spatial discretization.
+ * If a conservative method of interpolation is required ParaMEDMEM::MEDCouplingRemapper class should be used instead with "P1P0" method.
+ */
+MEDCouplingFieldDouble *MEDCouplingFieldDouble::nodeToCellDiscretization() const
+{
+  checkCoherency();
+  TypeOfField tf(getTypeOfField());
+  if(tf!=ON_NODES)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::nodeToCellDiscretization : this field is expected to be on ON_NODES !");
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret(clone(false));
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDiscretizationP0> nsp(new MEDCouplingFieldDiscretizationP0);
+  ret->setDiscretization(nsp);
+  const MEDCouplingMesh *m(getMesh());//m is non empty thanks to checkCoherency call
+  int nbCells(m->getNumberOfCells()),nbNodes(m->getNumberOfNodes());
+  std::vector<DataArrayDouble *> arrs(getArrays());
+  std::size_t sz(arrs.size());
+  std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> > outArrsSafe(sz); std::vector<DataArrayDouble *> outArrs(sz);
+  for(std::size_t j=0;j<sz;j++)
+    {
+      int nbCompo(arrs[j]->getNumberOfComponents());
+      outArrsSafe[j]=DataArrayDouble::New(); outArrsSafe[j]->alloc(nbCells,nbCompo);
+      outArrsSafe[j]->copyStringInfoFrom(*arrs[j]);
+      outArrs[j]=outArrsSafe[j];
+      double *pt(outArrsSafe[j]->getPointer());
+      const double *srcPt(arrs[j]->begin());
+      for(int i=0;i<nbCells;i++,pt+=nbCompo)
+        {
+          std::vector<int> nodeIds;
+          m->getNodeIdsOfCell(i,nodeIds);
+          std::fill(pt,pt+nbCompo,0.);
+          std::size_t nbNodesInCell(nodeIds.size());
+          for(std::size_t k=0;k<nbNodesInCell;k++)
+            std::transform(srcPt+nodeIds[k]*nbCompo,srcPt+(nodeIds[k]+1)*nbCompo,pt,pt,std::plus<double>());
+          if(nbNodesInCell!=0)
+            std::transform(pt,pt+nbCompo,pt,std::bind2nd(std::multiplies<double>(),1./((double)nbNodesInCell)));
+          else
+            {
+              std::ostringstream oss; oss << "MEDCouplingFieldDouble::nodeToCellDiscretization : Cell id #" << i << " has been detected to have no nodes !";
+              throw INTERP_KERNEL::Exception(oss.str().c_str());
+            }
+        }
+    }
+  ret->setArrays(outArrs);
+  return ret.retn();
+}
+
 /*!
  * Copies tiny info (component names, name and description) from an \a other field to
  * \a this one.
index d0696a8042c6c69f167e0ac9d1834823cc7289a3..b64819edca9f96b709938472c48a4e44d1b92f42 100644 (file)
@@ -63,6 +63,7 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT MEDCouplingFieldDouble *clone(bool recDeepCpy) const;
     MEDCOUPLING_EXPORT MEDCouplingFieldDouble *cloneWithMesh(bool recDeepCpy) const;
     MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildNewTimeReprFromThis(TypeOfTimeDiscretization td, bool deepCopy) const;
+    MEDCOUPLING_EXPORT MEDCouplingFieldDouble *nodeToCellDiscretization() const;
     MEDCOUPLING_EXPORT TypeOfTimeDiscretization getTimeDiscretization() const;
     MEDCOUPLING_EXPORT void checkCoherency() const;
     MEDCOUPLING_EXPORT void setNature(NatureOfField nat);
index 9b2a1527b7ed840408fca4290b05d8bf0575104c..b560dbafbae48c50c42e7723e27b7cd74f12a247 100644 (file)
@@ -13830,6 +13830,27 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertEqual(f4.getMesh(),None)
         pass
 
+    # test a simple node to cell convertion of a field
+    def testSwig2NodeToCellDiscretization1(self):
+        f=MEDCouplingFieldDouble(ON_NODES) ; f.setTime(1.1,2,3)
+        a1=DataArrayDouble(4) ; a1.iota()
+        a2=DataArrayDouble(3) ; a2.iota()
+        m=MEDCouplingCMesh() ; m.setCoords(a1,a2)
+        f.setMesh(m)
+        arr=DataArrayDouble([21.,121.,20.,120.,19.,119.,18.,118.,17.,117.,16.,116.,15.,115.,14.,114.,13.,113.,12.,112.,11.,111.,10.,110.],12,2) ; arr.setInfoOnComponents(["aa [km]","bbb [kJ]"])
+        f.setArray(arr) ; f.setName("toto")
+        #
+        f2=f.nodeToCellDiscretization()
+        self.assertEqual(ON_CELLS,f2.getTypeOfField())
+        self.assertEqual("toto",f2.getName())
+        self.assertEqual([1.1,2,3],f2.getTime())
+        self.assertEqual(["aa [km]","bbb [kJ]"],f2.getArray().getInfoOnComponents())
+        self.assertEqual(6,f2.getArray().getNumberOfTuples())
+        self.assertEqual(f.getMesh().getHiddenCppPointer(),f2.getMesh().getHiddenCppPointer())
+        exp=DataArrayDouble([18.5,118.5,17.5,117.5,16.5,116.5,14.5,114.5,13.5,113.5,12.5,112.5],6,2) ; exp.setInfoOnComponents(["aa [km]","bbb [kJ]"])
+        self.assertTrue(f2.getArray().isEqual(exp,1e-13))
+        pass
+
     def setUp(self):
         pass
     pass
index a511180cd2de7b0ab68a1816177efc38f92b959b..5f62bec585d63fa041f10eb6d4748cb67fe49416 100644 (file)
@@ -207,6 +207,7 @@ using namespace INTERP_KERNEL;
 %newobject ParaMEDMEM::MEDCouplingFieldDouble::cloneWithMesh;
 %newobject ParaMEDMEM::MEDCouplingFieldDouble::deepCpy;
 %newobject ParaMEDMEM::MEDCouplingFieldDouble::buildNewTimeReprFromThis;
+%newobject ParaMEDMEM::MEDCouplingFieldDouble::nodeToCellDiscretization;
 %newobject ParaMEDMEM::MEDCouplingFieldDouble::getValueOnMulti;
 %newobject ParaMEDMEM::MEDCouplingFieldTemplate::New;
 %newobject ParaMEDMEM::DataArray::deepCpy;
@@ -3391,6 +3392,7 @@ namespace ParaMEDMEM
     MEDCouplingFieldDouble *cloneWithMesh(bool recDeepCpy) const;
     MEDCouplingFieldDouble *deepCpy() const;
     MEDCouplingFieldDouble *buildNewTimeReprFromThis(TypeOfTimeDiscretization td, bool deepCpy) const throw(INTERP_KERNEL::Exception);
+    MEDCouplingFieldDouble *nodeToCellDiscretization() const throw(INTERP_KERNEL::Exception);
     TypeOfTimeDiscretization getTimeDiscretization() const throw(INTERP_KERNEL::Exception);
     double getIJ(int tupleId, int compoId) const throw(INTERP_KERNEL::Exception);
     double getIJK(int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception);