From 3624bce9198169126ac994fa31dbc25db2337597 Mon Sep 17 00:00:00 2001 From: ageay Date: Wed, 27 Mar 2013 08:54:41 +0000 Subject: [PATCH] Mesh::ComputeIsoBarycenterOfNodesPerCell --- src/MEDCoupling/MEDCouplingCMesh.cxx | 5 ++ src/MEDCoupling/MEDCouplingCMesh.hxx | 1 + .../MEDCouplingCurveLinearMesh.cxx | 5 ++ .../MEDCouplingCurveLinearMesh.hxx | 1 + src/MEDCoupling/MEDCouplingExtrudedMesh.cxx | 8 +- src/MEDCoupling/MEDCouplingExtrudedMesh.hxx | 1 + src/MEDCoupling/MEDCouplingMesh.hxx | 1 + src/MEDCoupling/MEDCouplingUMesh.cxx | 84 ++++++++++++++++++- src/MEDCoupling/MEDCouplingUMesh.hxx | 1 + src/MEDCoupling/MEDCouplingUMeshDesc.cxx | 4 + src/MEDCoupling/MEDCouplingUMeshDesc.hxx | 1 + src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 26 ++++++ src/MEDCoupling_Swig/MEDCouplingCommon.i | 2 + 13 files changed, 135 insertions(+), 5 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingCMesh.cxx b/src/MEDCoupling/MEDCouplingCMesh.cxx index 1ae822817..1e2e0df00 100644 --- a/src/MEDCoupling/MEDCouplingCMesh.cxx +++ b/src/MEDCoupling/MEDCouplingCMesh.cxx @@ -650,6 +650,11 @@ DataArrayDouble *MEDCouplingCMesh::getBarycenterAndOwner() const return ret; } +DataArrayDouble *MEDCouplingCMesh::computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception) +{ + return MEDCouplingCMesh::getBarycenterAndOwner(); +} + void MEDCouplingCMesh::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception) { throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for CMesh !"); diff --git a/src/MEDCoupling/MEDCouplingCMesh.hxx b/src/MEDCoupling/MEDCouplingCMesh.hxx index e1f25b51a..368f9c335 100644 --- a/src/MEDCoupling/MEDCouplingCMesh.hxx +++ b/src/MEDCoupling/MEDCouplingCMesh.hxx @@ -73,6 +73,7 @@ namespace ParaMEDMEM MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const; DataArrayDouble *getCoordinatesAndOwner() const; DataArrayDouble *getBarycenterAndOwner() const; + DataArrayDouble *computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception); void renumberCells(const int *old2NewBg, bool check=true) throw(INTERP_KERNEL::Exception); //some useful methods void getSplitCellValues(int *res) const; diff --git a/src/MEDCoupling/MEDCouplingCurveLinearMesh.cxx b/src/MEDCoupling/MEDCouplingCurveLinearMesh.cxx index be5a38af3..eb98ac357 100644 --- a/src/MEDCoupling/MEDCouplingCurveLinearMesh.cxx +++ b/src/MEDCoupling/MEDCouplingCurveLinearMesh.cxx @@ -706,6 +706,11 @@ DataArrayDouble *MEDCouplingCurveLinearMesh::getBarycenterAndOwner() const } } +DataArrayDouble *MEDCouplingCurveLinearMesh::computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception) +{ + return MEDCouplingCurveLinearMesh::getBarycenterAndOwner(); +} + /*! * \param [in,out] bary Barycenter array feeded with good values. * \sa MEDCouplingCurveLinearMesh::getBarycenterAndOwner diff --git a/src/MEDCoupling/MEDCouplingCurveLinearMesh.hxx b/src/MEDCoupling/MEDCouplingCurveLinearMesh.hxx index cab1c6739..33de0cafc 100644 --- a/src/MEDCoupling/MEDCouplingCurveLinearMesh.hxx +++ b/src/MEDCoupling/MEDCouplingCurveLinearMesh.hxx @@ -74,6 +74,7 @@ namespace ParaMEDMEM MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const; DataArrayDouble *getCoordinatesAndOwner() const; DataArrayDouble *getBarycenterAndOwner() const; + DataArrayDouble *computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception); void renumberCells(const int *old2NewBg, bool check=true) throw(INTERP_KERNEL::Exception); //some useful methods void getSplitCellValues(int *res) const; diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx index e2a184545..a513cbebe 100644 --- a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx +++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx @@ -734,8 +734,12 @@ DataArrayDouble *MEDCouplingExtrudedMesh::getCoordinatesAndOwner() const DataArrayDouble *MEDCouplingExtrudedMesh::getBarycenterAndOwner() const { - //not yet implemented - return 0; + throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::getBarycenterAndOwner : not yet implemented !"); +} + +DataArrayDouble *MEDCouplingExtrudedMesh::computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception) +{ + throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::computeIsoBarycenterOfNodesPerCell: not yet implemented !"); } void MEDCouplingExtrudedMesh::computeExtrusionAlg(const MEDCouplingUMesh *mesh3D) throw(INTERP_KERNEL::Exception) diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx index 951a591a4..572e82722 100644 --- a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx +++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx @@ -93,6 +93,7 @@ namespace ParaMEDMEM MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const; DataArrayDouble *getCoordinatesAndOwner() const; DataArrayDouble *getBarycenterAndOwner() const; + DataArrayDouble *computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception); //Serialization unserialisation void getTinySerializationInformation(std::vector& tinyInfoD, std::vector& tinyInfo, std::vector& littleStrings) const; void resizeForUnserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector& littleStrings) const; diff --git a/src/MEDCoupling/MEDCouplingMesh.hxx b/src/MEDCoupling/MEDCouplingMesh.hxx index bf508c8eb..501123219 100644 --- a/src/MEDCoupling/MEDCouplingMesh.hxx +++ b/src/MEDCoupling/MEDCouplingMesh.hxx @@ -84,6 +84,7 @@ namespace ParaMEDMEM virtual int getMeshDimension() const = 0; virtual DataArrayDouble *getCoordinatesAndOwner() const = 0; virtual DataArrayDouble *getBarycenterAndOwner() const = 0; + virtual DataArrayDouble *computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception) = 0; virtual DataArrayInt *giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception) = 0; virtual DataArrayInt *computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception) = 0; virtual int getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const = 0; diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 4fa0e852d..5345dd05a 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -6384,9 +6384,13 @@ MEDCouplingMesh *MEDCouplingUMesh::mergeMyselfWith(const MEDCouplingMesh *other) } /*! - * Returns an array with this->getNumberOfCells() tuples and this->getSpaceDimension() dimension. - * The false barycenter is computed that is to say barycenter of a cell is computed using average on each - * components of coordinates of the cell. + * This method is ** very badly named ** : This method computes the center of inertia of each cells in \a this. + * So this method is not a right choice for degnerated meshes (not well oriented, cells with measure close to zero). + * + * \return a newly allocated DataArrayDouble instance that the caller has to deal with. The returned + * DataArrayDouble instance will have \c this->getNumberOfCells() tuples and \c this->getSpaceDimension() components. + * + * \sa MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell */ DataArrayDouble *MEDCouplingUMesh::getBarycenterAndOwner() const { @@ -6408,6 +6412,80 @@ DataArrayDouble *MEDCouplingUMesh::getBarycenterAndOwner() const return ret.retn(); } +/*! + * This method computes for each cell in \a this, the location of the iso barycenter of nodes constituting + * the cell. Contrary to badly named MEDCouplingUMesh::getBarycenterAndOwner method that returns the center of inertia of the + * + * \return a newly allocated DataArrayDouble instance that the caller has to deal with. The returned + * DataArrayDouble instance will have \c this->getNumberOfCells() tuples and \c this->getSpaceDimension() components. + * + * \sa MEDCouplingUMesh::getBarycenterAndOwner + * \throw If \a this is not fully defined (coordinates and connectivity) + * \throw If there is presence in nodal connectivity in \a this of node ids not in [0, \c this->getNumberOfNodes() ) + */ +DataArrayDouble *MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception) +{ + checkFullyDefined(); + MEDCouplingAutoRefCountObjectPtr ret=DataArrayDouble::New(); + int spaceDim=getSpaceDimension(); + int nbOfCells=getNumberOfCells(); + int nbOfNodes=getNumberOfNodes(); + ret->alloc(nbOfCells,spaceDim); + double *ptToFill=ret->getPointer(); + const int *nodal=_nodal_connec->getConstPointer(); + const int *nodalI=_nodal_connec_index->getConstPointer(); + const double *coor=_coords->getConstPointer(); + for(int i=0;i=0 && *conn()); + else + { + std::ostringstream oss; oss << "MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell : on cell #" << i << " presence of nodeId #" << *conn << " should be in [0," << nbOfNodes << ") !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } + } + int nbOfNodesInCell=nodalI[i+1]-nodalI[i]-1; + if(nbOfNodesInCell>0) + std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies(),1./(double)nbOfNodesInCell)); + else + { + std::ostringstream oss; oss << "MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell : on cell #" << i << " presence of cell with no nodes !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } + } + else + { + std::set s(nodal+nodalI[i]+1,nodal+nodalI[i+1]); + s.erase(-1); + for(std::set::const_iterator it=s.begin();it!=s.end();it++) + { + if(*it>=0 && *it()); + else + { + std::ostringstream oss; oss << "MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell : on cell polyhedron cell #" << i << " presence of nodeId #" << *it << " should be in [0," << nbOfNodes << ") !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } + } + if(!s.empty()) + std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies(),1./(double)s.size())); + else + { + std::ostringstream oss; oss << "MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell : on polyhedron cell #" << i << " there are no nodes !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } + } + } + return ret.retn(); +} + /*! * This method is similar to MEDCouplingUMesh::getBarycenterAndOwner except that it works on subPart of 'this' without * building explicitely it. The input part is defined by an array [begin,end). All ids contained in this array should be less than this->getNumberOfCells(). diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 7a85e35b8..98586f3c1 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -212,6 +212,7 @@ namespace ParaMEDMEM // MEDCOUPLING_EXPORT MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const; MEDCOUPLING_EXPORT DataArrayDouble *getBarycenterAndOwner() const; + MEDCOUPLING_EXPORT DataArrayDouble *computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayDouble *getPartBarycenterAndOwner(const int *begin, const int *end) const; MEDCOUPLING_EXPORT static MEDCouplingUMesh *Build0DMeshFromCoords(DataArrayDouble *da) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT static MEDCouplingUMesh *MergeUMeshes(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception); diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx index f492da9f1..5f2de9027 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx @@ -529,3 +529,7 @@ std::string MEDCouplingUMeshDesc::getVTKDataSetType() const throw(INTERP_KERNEL: throw INTERP_KERNEL::Exception("MEDCouplingUMeshDesc::getVTKDataSetType : not implemented yet !"); } +DataArrayDouble *MEDCouplingUMeshDesc::computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception) +{ + throw INTERP_KERNEL::Exception("MEDCouplingUMeshDesc::computeIsoBarycenterOfNodesPerCell : not implemented yet !"); +} diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx index 86ef28f73..ba3f674f8 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx @@ -90,6 +90,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT DataArrayInt *zipCoordsTraducer(); MEDCOUPLING_EXPORT MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const; MEDCOUPLING_EXPORT DataArrayDouble *getBarycenterAndOwner() const; + MEDCOUPLING_EXPORT DataArrayDouble *computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception); private: MEDCouplingUMeshDesc(); ~MEDCouplingUMeshDesc(); diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index 098bdc63e..bde0762f1 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -11839,6 +11839,32 @@ class MEDCouplingBasicsTest(unittest.TestCase): d.alloc(1000,3) ; d.fillWithValue(127) self.assertTrue(len(d.__repr__())<500) pass + + def testSwig2MeshComputeIsoBarycenterOfNodesPerCell1(self): + coo=DataArrayDouble([26.17509821414239,5.0374,200.,26.175098214142388,-5.0374,200.,17.450065476094927,20.1496,200.,8.725032738047464,25.187,200.,43.62516369023732,5.0374,200.,34.90013095218986,10.0748,200.,34.900130952189855,-10.0748,200.,43.625163690237315,-5.0374,200.,26.175098214142402,25.187,200.,26.175098214142395,35.2618,200.,17.45006547609493,40.2992,200.,8.725032738047469,35.2618,200.,26.17509821414239,5.0374,200.,26.175098214142388,-5.0374,200.,17.450065476094927,20.1496,200.,8.725032738047464,25.187,200.,43.62516369023732,5.0374,200.,34.90013095218986,10.0748,200.,34.900130952189855,-10.0748,200.,43.625163690237315,-5.0374,200.,26.175098214142402,25.187,200.,26.175098214142395,35.2618,200.,17.45006547609493,40.2992,200.,8.725032738047469,35.2618,200.],24,3) + m=MEDCouplingUMesh.New("toto",3) + m.allocateCells(0) + m.insertNextCell(NORM_POLYHED,[4,5,0,1,6,7,-1,19,18,13,12,17,16,-1,5,4,16,17,-1,0,5,17,12,-1,1,0,12,13,-1,6,1,13,18,-1,7,6,18,19,-1,4,7,19,16]) + m.insertNextCell(NORM_POLYHED,[9,10,11,3,2,8,-1,20,14,15,23,22,21,-1,10,9,21,22,-1,11,10,22,23,-1,3,11,23,15,-1,2,3,15,14,-1,8,2,14,20,-1,9,8,20,21]) + m.setCoords(coo) + m.checkCoherency1() + # + dReference=DataArrayDouble([(34.900130952189848,0.,200),(17.450065476094931,30.2244,200.)]) + self.assertTrue(m.computeIsoBarycenterOfNodesPerCell().isEqual(dReference,1e-12)) + m.getNodalConnectivity().setIJ(87,0,24) + self.assertRaises(InterpKernelException,m.computeIsoBarycenterOfNodesPerCell) + m.getNodalConnectivity().setIJ(87,0,-2) + self.assertRaises(InterpKernelException,m.computeIsoBarycenterOfNodesPerCell) + m.getNodalConnectivity().setIJ(87,0,21)# put again 21 as at the beginning + # + self.assertTrue(m.unPolyze()) + self.assertEqual([NORM_HEXGP12],m.getAllTypes()) + self.assertTrue(m.computeIsoBarycenterOfNodesPerCell().isEqual(dReference,1e-12)) + m.getNodalConnectivity().setIJ(25,0,24) + self.assertRaises(InterpKernelException,m.computeIsoBarycenterOfNodesPerCell) + m.getNodalConnectivity().setIJ(25,0,-1) + self.assertRaises(InterpKernelException,m.computeIsoBarycenterOfNodesPerCell) + pass def setUp(self): pass diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index f080f8266..ea7851815 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -299,6 +299,7 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::MEDCouplingMesh::giveCellsWithType; %newobject ParaMEDMEM::MEDCouplingMesh::getCoordinatesAndOwner; %newobject ParaMEDMEM::MEDCouplingMesh::getBarycenterAndOwner; +%newobject ParaMEDMEM::MEDCouplingMesh::computeIsoBarycenterOfNodesPerCell; %newobject ParaMEDMEM::MEDCouplingMesh::buildOrthogonalField; %newobject ParaMEDMEM::MEDCouplingMesh::getCellIdsFullyIncludedInNodeIds; %newobject ParaMEDMEM::MEDCouplingMesh::mergeMyselfWith; @@ -546,6 +547,7 @@ namespace ParaMEDMEM virtual int getMeshDimension() const throw(INTERP_KERNEL::Exception); virtual DataArrayDouble *getCoordinatesAndOwner() const throw(INTERP_KERNEL::Exception); virtual DataArrayDouble *getBarycenterAndOwner() const throw(INTERP_KERNEL::Exception); + virtual DataArrayDouble *computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception); virtual DataArrayInt *giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception); virtual DataArrayInt *computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception); virtual int getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception); -- 2.39.2