Salome HOME
Add MEDCouplingMesh.computeMeshCenterOfMass
authorAnthony Geay <anthony.geay@edf.fr>
Tue, 30 Mar 2021 14:05:37 +0000 (16:05 +0200)
committerAnthony Geay <anthony.geay@edf.fr>
Tue, 30 Mar 2021 14:05:37 +0000 (16:05 +0200)
src/MEDCoupling/MEDCouplingMesh.cxx
src/MEDCoupling/MEDCouplingMesh.hxx
src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py
src/MEDCoupling_Swig/MEDCouplingCommon.i

index 6de0b97a7523f71c40c3aaac3d71d7d4b165b1ee..58b55734826cda4818546331bf13c9775a52fed5 100755 (executable)
@@ -674,6 +674,22 @@ void MEDCouplingMesh::getCellsContainingPointsLinearPartOnlyOnNonDynType(const d
   this->getCellsContainingPoints(pos,nbOfPoints,eps,elts,eltsIndex);
 }
 
+/*!
+ * Method computing center of mass of whole mesh (\a this)
+ * \return DataArrayDouble * - a new instance of DataArrayDouble with one tuple of n components where n is space dimension
+ */
+MCAuto<DataArrayDouble> MEDCouplingMesh::computeMeshCenterOfMass() const
+{
+  MCAuto<DataArrayDouble> cellCenters( this->computeCellCenterOfMass() );
+  MCAuto<MEDCouplingFieldDouble> vol( this->getMeasureField(true) );
+  MCAuto<DataArrayDouble> volXCenter( DataArrayDouble::Multiply(cellCenters,vol->getArray()) );
+  MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(1, this->getSpaceDimension());
+  volXCenter->accumulate( ret->getPointer() );
+  double volOfMesh(vol->accumulate(0));
+  ret->applyLin(1.0/volOfMesh,0.0);
+  return ret;
+}
+
 /*!
  * Writes \a this mesh into a VTK format file named as specified.
  *  \param [in] fileName - the name of the file to write in. If the extension is OK the fileName will be used directly.
index 4a3f631331a58c727a758378d94247259c3930dd..f75ca7f4ed021de42acdaa313b2f30444a436ada 100644 (file)
@@ -106,6 +106,7 @@ namespace MEDCoupling
     MEDCOUPLING_EXPORT virtual std::string simpleRepr() const = 0;
     MEDCOUPLING_EXPORT virtual std::string advancedRepr() const = 0;
     // tools
+    MEDCOUPLING_EXPORT virtual MCAuto<DataArrayDouble> computeMeshCenterOfMass() const;
     MEDCOUPLING_EXPORT virtual const DataArrayDouble *getDirectAccessOfCoordsArrIfInStructure() const = 0;
     MEDCOUPLING_EXPORT virtual std::vector<mcIdType> getDistributionOfTypes() const = 0;
     MEDCOUPLING_EXPORT virtual DataArrayIdType *checkTypeConsistencyAndContig(const std::vector<mcIdType>& code, const std::vector<const DataArrayIdType *>& idsPerType) const = 0;
index 6b3955a860bc823c87e0f2e7756ae3f645b8a1ce..70dfa7bf15c82118fd9b9a66efb2705cf3103959 100644 (file)
@@ -956,6 +956,16 @@ class MEDCouplingBasicsTest7(unittest.TestCase):
         self.assertEqual(arr3.getInfoOnComponents(),comps)
         self.assertTrue(arr3.isEqual(arr))
 
+    def testComputeMeshCenterOfMass0(self):
+        #2D
+        arr = DataArrayDouble(5) ; arr.iota()
+        m = MEDCouplingCMesh() ; m.setCoords(arr,arr) ; m=m.buildUnstructured()
+        self.assertTrue( m.computeMeshCenterOfMass().isEqual(DataArrayDouble([2,2],1,2),1e-12) )
+        #3D
+        m = MEDCouplingCMesh() ; m.setCoords(arr,arr,arr) ; m=m.buildUnstructured()
+        self.assertTrue( m.computeMeshCenterOfMass().isEqual(DataArrayDouble([2,2,2],1,3),1e-12) )
+
+
     pass
 
 if __name__ == '__main__':
index ab90fc0f9164a6e8f5b0a7724f8754e580bd2493..1639e0751e394a059ab6de26b5abeb3a058562e9 100644 (file)
@@ -309,6 +309,7 @@ typedef long int mcIdType;
 %newobject MEDCoupling::MEDCouplingMesh::buildPartRange;
 %newobject MEDCoupling::MEDCouplingMesh::giveCellsWithType;
 %newobject MEDCoupling::MEDCouplingMesh::getCoordinatesAndOwner;
+%newobject MEDCoupling::MEDCouplingMesh::computeMeshCenterOfMass;
 %newobject MEDCoupling::MEDCouplingMesh::computeCellCenterOfMass;
 %newobject MEDCoupling::MEDCouplingMesh::computeIsoBarycenterOfNodesPerCell;
 %newobject MEDCoupling::MEDCouplingMesh::buildOrthogonalField;
@@ -750,6 +751,12 @@ namespace MEDCoupling
          {
            return self->simpleRepr();
          }
+         
+          DataArrayDouble *computeMeshCenterOfMass() const
+          {
+            MCAuto<DataArrayDouble> ret(self->computeMeshCenterOfMass());
+            return ret.retn();
+          }
 
          PyObject *getTime()
          {