From b607ffc713080a567fb90595118069ac18181e99 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Tue, 21 Apr 2020 08:45:31 +0200 Subject: [PATCH] ParaUMesh.redistributeCells implementation. --- src/MEDCoupling/MCAuto.hxx | 30 ++++- src/MEDCoupling/MEDCouplingMemArray.hxx | 14 ++- src/MEDCoupling/MEDCouplingMemArray.txx | 54 ++++++--- src/MEDCoupling/MEDCouplingSkyLineArray.cxx | 71 ++++++++++++ src/MEDCoupling/MEDCouplingSkyLineArray.hxx | 27 ++++- src/MEDCoupling/MEDCouplingTraits.hxx | 2 +- src/MEDCoupling/MEDCouplingUMesh.hxx | 30 ++++- .../MEDCouplingBasicsTest7.py | 35 +++++- src/MEDCoupling_Swig/MEDCouplingCommon.i | 11 ++ .../MEDCouplingDataArrayTypemaps.i | 11 ++ src/MEDCoupling_Swig/MEDCouplingMemArray.i | 2 + src/MEDLoader/MEDFileBasis.hxx | 1 + src/ParaMEDMEM/CommInterface.cxx | 74 ++++++------ src/ParaMEDMEM/CommInterface.hxx | 102 ++++++++++++++++- src/ParaMEDMEM/ParaUMesh.cxx | 108 +++++++++++++++--- src/ParaMEDMEM/ParaUMesh.hxx | 14 ++- src/ParaMEDMEM_Swig/ParaMEDMEMCommon.i | 57 ++++++++- 17 files changed, 549 insertions(+), 94 deletions(-) diff --git a/src/MEDCoupling/MCAuto.hxx b/src/MEDCoupling/MCAuto.hxx index 1d83037b9..7e92f0653 100644 --- a/src/MEDCoupling/MCAuto.hxx +++ b/src/MEDCoupling/MCAuto.hxx @@ -32,14 +32,14 @@ namespace MEDCoupling class MCAuto { public: - static MCAuto TakeRef(T *ptr) { MCAuto ret(MCAuto(nullptr)); ret.takeRef(ptr); return ret; } - MCAuto(const MCAuto& other):_ptr(0) { referPtr(other._ptr); } - MCAuto(T *ptr=0):_ptr(ptr) { } + static MCAuto TakeRef(T *ptr) { MCAuto ret; ret.takeRef(ptr); return ret; } + MCAuto(const MCAuto& other):_ptr(nullptr) { referPtr(other._ptr); } + MCAuto(T *ptr=nullptr):_ptr(ptr) { } ~MCAuto() { destroyPtr(); } void checkNotNull() const { if(!_ptr) throw INTERP_KERNEL::Exception("Pointer is nullptr !"); } - bool isNull() const { return _ptr==0; } + bool isNull() const { return _ptr==nullptr; } bool isNotNull() const { return !isNull(); } - void nullify() { destroyPtr(); _ptr=0; } + void nullify() { destroyPtr(); _ptr=nullptr; } bool operator==(const MCAuto& other) const { return _ptr==other._ptr; } bool operator==(const T *other) const { return _ptr==other; } MCAuto &operator=(const MCAuto& other) { if(_ptr!=other._ptr) { destroyPtr(); referPtr(other._ptr); } return *this; } @@ -71,6 +71,26 @@ namespace MEDCoupling return ret; } + template + std::vector> FromVecToVecAuto(const std::vector& inputVec) + { + std::size_t size(inputVec.size()); + std::vector> ret(size); + typename std::vector>::iterator itArrays(ret.begin()); + std::for_each(inputVec.begin(),inputVec.end(),[&itArrays](T *elt) { (*itArrays).takeRef(elt); itArrays++; }); + return ret; + } + + template + std::vector> FromVecConstToVecAuto(const std::vector& inputVec) + { + std::size_t size(inputVec.size()); + std::vector> ret(size); + typename std::vector>::iterator itArrays(ret.begin()); + std::for_each(inputVec.begin(),inputVec.end(),[&itArrays](const T *elt) { (*itArrays).takeRef(const_cast(elt)); itArrays++; }); + return ret; + } + template typename MEDCoupling::MCAuto DynamicCast(typename MEDCoupling::MCAuto& autoSubPtr) throw() { diff --git a/src/MEDCoupling/MEDCouplingMemArray.hxx b/src/MEDCoupling/MEDCouplingMemArray.hxx index f9a847cc4..501cea75e 100755 --- a/src/MEDCoupling/MEDCouplingMemArray.hxx +++ b/src/MEDCoupling/MEDCouplingMemArray.hxx @@ -165,6 +165,7 @@ namespace MEDCoupling void setPartOfValuesBase3(const DataArray *aBase, const mcIdType *bgTuples, const mcIdType *endTuples, mcIdType bgComp, mcIdType endComp, mcIdType stepComp, bool strictCompoCompare=true); virtual void *getVoidStarPointer() = 0; virtual DataArray *deepCopy() const = 0; + virtual DataArray *copySorted(bool asc=true) const = 0; virtual DataArray *buildNewEmptyInstance() const = 0; virtual bool isAllocated() const = 0; virtual void checkAllocated() const = 0; @@ -231,6 +232,7 @@ namespace MEDCoupling typedef T Type; public: static MCAuto< typename Traits::ArrayTypeCh > NewFromStdVector(const typename std::vector& v); + static MCAuto< typename Traits::ArrayTypeCh > NewFromArray(const T *arrBegin, const T *arrEnd); std::vector< MCAuto< typename Traits::ArrayTypeCh > > explodeComponents() const; // std::size_t getHeapMemorySizeWithoutChildren() const; @@ -310,6 +312,7 @@ namespace MEDCoupling MemArray& accessToMemArray() { return _mem; } const MemArray& accessToMemArray() const { return _mem; } protected: + typename Traits::ArrayTypeCh *copySortedImpl(bool asc) const; typename Traits::ArrayType *mySelectByTupleId(const mcIdType *new2OldBg, const mcIdType *new2OldEnd) const; typename Traits::ArrayType *mySelectByTupleId(const DataArrayIdType& di) const; typename Traits::ArrayType *mySelectByTupleIdSafe(const mcIdType *new2OldBg, const mcIdType *new2OldEnd) const; @@ -392,6 +395,7 @@ namespace MEDCoupling static DataArrayFloat *New(); public:// abstract method overload DataArrayFloat *deepCopy() const; + DataArrayFloat *copySorted(bool asc=true) const override { return this->copySortedImpl(asc); } std::string getClassName() const override { return std::string("DataArrayFloat"); } DataArrayFloat *buildNewEmptyInstance() const { return DataArrayFloat::New(); } DataArrayFloat *selectByTupleRanges(const std::vector >& ranges) const { return DataArrayTemplateFP::mySelectByTupleRanges(ranges); } @@ -423,6 +427,7 @@ namespace MEDCoupling static DataArrayDouble *New(); double doubleValue() const; DataArrayDouble *deepCopy() const; + DataArrayDouble *copySorted(bool asc=true) const override { return this->copySortedImpl(asc); } std::string getClassName() const override { return std::string("DataArrayDouble"); } DataArrayDouble *buildNewEmptyInstance() const { return DataArrayDouble::New(); } void checkMonotonic(bool increasing, double eps) const; @@ -544,7 +549,6 @@ namespace MEDCoupling using DataArrayType = typename Traits::ArrayType; public: static DataArrayType *New(); - static MCAuto NewFromArray(const T *arrBegin, const T *arrEnd); T intValue() const; bool isEqual(const DataArrayDiscrete& other) const; bool isEqualIfNotWhy(const DataArrayDiscrete& other, std::string& reason) const; @@ -703,8 +707,9 @@ namespace MEDCoupling { friend class DataArrayDiscrete; public: - DataArrayInt32 *deepCopy() const;//ok - DataArrayInt32 *buildNewEmptyInstance() const { return DataArrayInt32::New(); }//ok + DataArrayInt32 *deepCopy() const; + DataArrayInt32 *copySorted(bool asc=true) const override { return this->copySortedImpl(asc); } + DataArrayInt32 *buildNewEmptyInstance() const { return DataArrayInt32::New(); } public: DataArrayInt32 *selectByTupleId(const mcIdType *new2OldBg, const mcIdType *new2OldEnd) const { return this->mySelectByTupleId(new2OldBg,new2OldEnd); } DataArrayInt32 *selectByTupleId(const DataArrayIdType& di) const { return this->mySelectByTupleId(di); } @@ -725,6 +730,7 @@ namespace MEDCoupling friend class DataArrayDiscrete; public: DataArrayInt64 *deepCopy() const; + DataArrayInt64 *copySorted(bool asc=true) const override { return this->copySortedImpl(asc); } DataArrayInt64 *buildNewEmptyInstance() const { return DataArrayInt64::New(); }//ok public: DataArrayInt64 *selectByTupleId(const mcIdType *new2OldBg, const mcIdType *new2OldEnd) const { return this->mySelectByTupleId(new2OldBg,new2OldEnd); } @@ -814,6 +820,7 @@ namespace MEDCoupling DataArrayChar *buildEmptySpecializedDAChar() const; DataArrayByteIterator *iterator(); DataArrayByte *deepCopy() const; + DataArrayByte *copySorted(bool asc=true) const override { return this->copySortedImpl(asc); } DataArrayByte *performCopyOrIncrRef(bool deepCopy) const; DataArrayByte *buildNewEmptyInstance() const { return DataArrayByte::New(); } char byteValue() const; @@ -843,6 +850,7 @@ namespace MEDCoupling DataArrayChar *buildEmptySpecializedDAChar() const; DataArrayAsciiCharIterator *iterator(); DataArrayAsciiChar *deepCopy() const; + DataArrayAsciiChar *copySorted(bool asc=true) const override { throw INTERP_KERNEL::Exception("DataArrayAsciiChar::copySorted : not implemented for DataArrayByte"); } DataArrayAsciiChar *performCopyOrIncrRef(bool deepCopy) const; DataArrayAsciiChar *buildNewEmptyInstance() const { return DataArrayAsciiChar::New(); } char asciiCharValue() const; diff --git a/src/MEDCoupling/MEDCouplingMemArray.txx b/src/MEDCoupling/MEDCouplingMemArray.txx index 373e9201b..70aefc7a5 100755 --- a/src/MEDCoupling/MEDCouplingMemArray.txx +++ b/src/MEDCoupling/MEDCouplingMemArray.txx @@ -672,6 +672,20 @@ namespace MEDCoupling std::copy(v.begin(),v.end(),pt); return ret; } + + /*! + * Returns a newly created array containing a copy of the input array defined by [ \a arrBegin, \a arrEnd ) + */ + template + MCAuto< typename Traits::ArrayTypeCh > DataArrayTemplate::NewFromArray(const T *arrBegin, const T *arrEnd) + { + using DataArrayT = typename Traits::ArrayTypeCh; + MCAuto< DataArrayT > ret(DataArrayT::New()); + std::size_t nbElts(std::distance(arrBegin,arrEnd)); + ret->alloc(nbElts,1); + std::copy(arrBegin,arrEnd,ret->getPointer()); + return ret; + } template std::vector< MCAuto< typename Traits::ArrayTypeCh > > DataArrayTemplate::explodeComponents() const @@ -1087,10 +1101,12 @@ namespace MEDCoupling } /*! - * Sorts values of the array. + * Sorts values of the array. \b Warning, this method is not const, it alterates \a this content. + * * \param [in] asc - \a true means ascending order, \a false, descending. * \throw If \a this is not allocated. * \throw If \a this->getNumberOfComponents() != 1. + * \sa copySorted */ template void DataArrayTemplate::sort(bool asc) @@ -1105,6 +1121,23 @@ namespace MEDCoupling declareAsNew(); } + /*! + * Sorts values of the array and put the result in a newly allocated returned array. + * This method does not alterate \a this content. + * + * \param [in] asc - \a true means ascending order, \a false, descending. + * \throw If \a this is not allocated. + * \throw If \a this->getNumberOfComponents() != 1. + * \sa sort + */ + template + typename Traits::ArrayTypeCh *DataArrayTemplate::copySortedImpl(bool asc) const + { + MCAuto::ArrayTypeCh> ret(static_cast::ArrayTypeCh *>(this->deepCopy())); + ret->sort(asc); + return ret.retn(); + } + /*! * Returns a copy of \a this array with values permuted as required by \a old2New array. * The values are permuted so that \c new[ \a old2New[ i ]] = \c old[ i ]. @@ -3634,19 +3667,6 @@ struct NotInRange return new typename Traits::ArrayType; } - /*! - * Returns a newly created array containing a copy of the input array defined by [ \a arrBegin, \a arrEnd ) - */ - template - MCAuto< typename Traits::ArrayType > DataArrayDiscrete::NewFromArray(const T *arrBegin, const T *arrEnd) - { - MCAuto< typename Traits::ArrayType > ret(DataArrayDiscrete::New()); - std::size_t nbElts(std::distance(arrBegin,arrEnd)); - ret->alloc(nbElts,1); - std::copy(arrBegin,arrEnd,ret->getPointer()); - return ret; - } - /*! * Checks if values of \a this and another DataArrayInt are equal. For more info see * \ref MEDCouplingArrayBasicsCompare. @@ -5715,8 +5735,8 @@ struct NotInRange * - \a this : [0, 1, 1, 2, 2, 3, 4, 4, 5, 5, 5, 11] * - \a return is : [0, 1, 3, 5, 6, 8, 11, 12] * - * \return a newly allocated array containing the indexed array of - * \throw if \a this is not allocated or if \a this has not exactly one component or if number of tuples is equal to 0. + * \return a newly allocated array containing the indexed array format of groups by same consecutive value. + * \throw if \a this is not allocated or if \a this has not exactly one component. * \sa DataArrayInt::buildUnique, MEDCouplingSkyLineArray::groupPacks */ template @@ -5725,8 +5745,6 @@ struct NotInRange this->checkAllocated(); if(this->getNumberOfComponents()!=1) throw INTERP_KERNEL::Exception("DataArrayInt::indexOfSameConsecutiveValueGroups : only single component allowed !"); - if(this->getNumberOfTuples()==0) - throw INTERP_KERNEL::Exception("DataArrayInt::indexOfSameConsecutiveValueGroups : number of tuples must be > 0 !"); const T *pt(this->begin()); const T *const ptEnd(this->end()) , * const ptBg(this->begin()); const T *oldPt(pt); diff --git a/src/MEDCoupling/MEDCouplingSkyLineArray.cxx b/src/MEDCoupling/MEDCouplingSkyLineArray.cxx index 5fdc66fa8..8ab51c2d6 100755 --- a/src/MEDCoupling/MEDCouplingSkyLineArray.cxx +++ b/src/MEDCoupling/MEDCouplingSkyLineArray.cxx @@ -213,6 +213,19 @@ DataArrayIdType* MEDCouplingSkyLineArray::getValuesArray() const return const_cast(this)->_values; } +MEDCouplingSkyLineArray *MEDCouplingSkyLineArray::deepCopy() const +{ + MCAuto indexCpy(this->_index->deepCopy()); + MCAuto valuesCpy(this->_values->deepCopy()); + MCAuto ret(MEDCouplingSkyLineArray::New(indexCpy,valuesCpy)); + if(_super_index.isNotNull()) + { + MCAuto superIndexCpy(this->_super_index->deepCopy()); + ret->_super_index = superIndexCpy; + } + return ret.retn(); +} + void MEDCouplingSkyLineArray::checkSuperIndex(const std::string& func) const { if (!_super_index->getNbOfElems()) @@ -347,6 +360,64 @@ MEDCouplingSkyLineArray *MEDCouplingSkyLineArray::uniqueNotSortedByPack() const MCAuto ret(MEDCouplingSkyLineArray::New(retIndex,retValues)); return ret.retn(); } +/*! + * Take as input skylinearrays containing the same number of packs ( \a this->getNumberOf ). + * For each packs, this method aggregates corresponding pack in \a sks. + * + * \throw if either a lenght of not nullptr instances in \a sks is zero or if number of packs of not nullptr instances in \a sks is not the same. + * \return a newly allocated skyline array that is the aggregated packs of inputs + */ +MEDCouplingSkyLineArray *MEDCouplingSkyLineArray::AggregatePacks(const std::vector& sks) +{ + std::vectorsksEff; + mcIdType nbOfPacks(std::numeric_limits::max()); + constexpr char MSG[]="MEDCouplingSkyLineArray::AggregatePacks : "; + for(auto sk : sks) + { + if(sk) + { + mcIdType curNbPacks(sk->getNumberOf()); + if(sksEff.empty()) + nbOfPacks = curNbPacks; + if(nbOfPacks != curNbPacks) + { + std::ostringstream oss; oss << MSG << "first not null input ska has " << nbOfPacks << " whereas there is presence of ska with " << curNbPacks << " !"; + throw INTERP_KERNEL::Exception(oss.str()); + } + sksEff.push_back(sk); + } + } + if(sksEff.empty()) + { + std::ostringstream oss; oss << MSG << "input vector contains no not nullptr elements !"; + throw INTERP_KERNEL::Exception(oss.str()); + } + // + MCAuto index(DataArrayIdType::New()); index->alloc(nbOfPacks+1,1); + mcIdType *indexPtr(index->getPointer()); *indexPtr=0; + std::vector indicesIn(SkyLineArrayIndexIterator(0,&sksEff),SkyLineArrayIndexIterator(sksEff.size(),&sksEff)); + for( mcIdType packId = 0 ; packId < nbOfPacks ; ++packId, ++indexPtr ) + { + mcIdType nbOfAggPacks(0); + std::for_each(indicesIn.begin(),indicesIn.end(),[packId,&nbOfAggPacks](const mcIdType *elt) { nbOfAggPacks+=elt[packId+1]-elt[packId]; }); + indexPtr[1] = indexPtr[0] + nbOfAggPacks; + } + mcIdType nbOfTuplesOut(index->back()); + MCAuto values(DataArrayIdType::New()); values->alloc(nbOfTuplesOut,1); + mcIdType *valuesPtr(values->getPointer()); + // let's go to populate values array + std::vector valuesIn(SkyLineArrayValuesIterator(0,&sksEff),SkyLineArrayValuesIterator(sksEff.size(),&sksEff)); + for( mcIdType packId = 0 ; packId < nbOfPacks ; ++packId ) + { + std::size_t pos(0); + std::for_each(valuesIn.begin(),valuesIn.end(),[packId,&indicesIn,&valuesPtr,&pos](const mcIdType *elt) + { valuesPtr=std::copy(elt+indicesIn[pos][packId],elt+indicesIn[pos][packId+1],valuesPtr); ++pos; } + ); + } + // + MCAuto ret(MEDCouplingSkyLineArray::New(index,values)); + return ret.retn(); +} /** * For a 2- or 3-level SkyLine array, return a copy of the absolute pack with given identifier. diff --git a/src/MEDCoupling/MEDCouplingSkyLineArray.hxx b/src/MEDCoupling/MEDCouplingSkyLineArray.hxx index a56d94b9a..b4f0d2509 100644 --- a/src/MEDCoupling/MEDCouplingSkyLineArray.hxx +++ b/src/MEDCoupling/MEDCouplingSkyLineArray.hxx @@ -17,8 +17,7 @@ // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -#ifndef __PARAMEDMEM_MEDCOUPLINGSKYLINEARRAY_HXX__ -#define __PARAMEDMEM_MEDCOUPLINGSKYLINEARRAY_HXX__ +#pragma once #include "MEDCoupling.hxx" #include "MEDCouplingMemArray.hxx" @@ -105,10 +104,13 @@ namespace MEDCoupling DataArrayIdType* getIndexArray() const; DataArrayIdType* getValuesArray() const; + MEDCouplingSkyLineArray *deepCopy() const; + std::string simpleRepr() const; MEDCouplingSkyLineArray *groupPacks(const DataArrayIdType *indexedPacks) const; MEDCouplingSkyLineArray *uniqueNotSortedByPack() const; + static MEDCouplingSkyLineArray *AggregatePacks(const std::vector& sks); void getSimplePackSafe(const mcIdType absolutePackId, std::vector & pack) const; const mcIdType * getSimplePackSafePtr(const mcIdType absolutePackId, mcIdType & packSize) const; @@ -141,5 +143,24 @@ namespace MEDCoupling MCAuto _values; }; + template + class SkyLineArrayGenIterator : public std::iterator< std::input_iterator_tag, const mcIdType *, mcIdType, const mcIdType **, const mcIdType *> + { + std::size_t _num = 0; + std::vector *_data = nullptr; + public: + explicit SkyLineArrayGenIterator(std::size_t num , std::vector *data) : _num(num),_data(data) {} + SkyLineArrayGenIterator& operator++() { ++_num; return *this; } + bool operator==(const SkyLineArrayGenIterator& other) const { return _num == other._num; } + bool operator!=(const SkyLineArrayGenIterator& other) const { return !(*this == other); } + reference operator*() const { T tt; return tt((*_data)[_num]); } + }; + + struct SkyLineArrayIndexPtrFunctor { const mcIdType *operator()(const MEDCouplingSkyLineArray *ska) { return ska->getIndex(); } }; + + using SkyLineArrayIndexIterator = SkyLineArrayGenIterator; + + struct SkyLineArrayValuesPtrFunctor { const mcIdType *operator()(const MEDCouplingSkyLineArray *ska) { return ska->getValues(); } }; + + using SkyLineArrayValuesIterator = SkyLineArrayGenIterator; } -# endif diff --git a/src/MEDCoupling/MEDCouplingTraits.hxx b/src/MEDCoupling/MEDCouplingTraits.hxx index fad3d7cf7..806c58efa 100644 --- a/src/MEDCoupling/MEDCouplingTraits.hxx +++ b/src/MEDCoupling/MEDCouplingTraits.hxx @@ -30,7 +30,7 @@ namespace MEDCoupling template struct MEDCOUPLING_EXPORT Traits { - typedef T EltType; + using EltType = T; }; class DataArrayInt32; diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 778411e29..c09784c41 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -18,8 +18,7 @@ // // Author : Anthony Geay (CEA/DEN) -#ifndef __PARAMEDMEM_MEDCOUPLINGUMESH_HXX__ -#define __PARAMEDMEM_MEDCOUPLINGUMESH_HXX__ +#pragma once #include "MEDCoupling.hxx" #include "MEDCouplingPointSet.hxx" @@ -442,6 +441,31 @@ namespace MEDCoupling mcIdType _conn_lgth; static const int NOTICABLE_FIRST_VAL=-7; }; + + template + class UMeshGenIterator : public std::iterator< std::input_iterator_tag, const TOUT *, mcIdType, const TOUT **, const TOUT *> + { + std::size_t _num = 0; + std::vector *_data = nullptr; + using my_reference = typename std::iterator< std::input_iterator_tag, const TOUT *, mcIdType, const TOUT **, const TOUT *>::reference; + public: + explicit UMeshGenIterator(std::size_t num , std::vector *data) : _num(num),_data(data) {} + UMeshGenIterator& operator++() { ++_num; return *this; } + bool operator==(const UMeshGenIterator& other) const { return _num == other._num; } + bool operator!=(const UMeshGenIterator& other) const { return !(*this == other); } + my_reference operator*() const { T tt; return tt((*_data)[_num]); } + }; + + struct UMeshIndexConnectivityFunctor { const DataArrayIdType *operator()(const MEDCouplingUMesh *um) { return um->getNodalConnectivityIndex(); } }; + + using UMeshConnectivityIndexIterator = UMeshGenIterator; + + struct UMeshConnectivityFunctor { const DataArrayIdType *operator()(const MEDCouplingUMesh *um) { return um->getNodalConnectivity(); } }; + + using UMeshConnectivityIterator = UMeshGenIterator; + + struct UMeshCoordsFunctor { const DataArrayDouble *operator()(const MEDCouplingUMesh *um) { return um->getCoords(); } }; + + using UMeshCoordsIterator = UMeshGenIterator; } -#endif diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py index ed473e56b..702c80303 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py @@ -1,10 +1,10 @@ # -*- coding: utf-8 -*- -# Copyright (C) 2007-2020 CEA/DEN, EDF R&D +# Copyright (C) 2007-2020 CEA/DEN,EDF R&D # # This library is free software; you can redistribute it and/or # modify it under the terms of the GNU Lesser General Public # License as published by the Free Software Foundation; either -# version 2.1 of the License, or (at your option) any later version. +# version 2.1 of the License,or (at your option) any later version. # # This library is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -12,8 +12,8 @@ # Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public -# License along with this library; if not, write to the Free Software -# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# License along with this library; if not,write to the Free Software +# Foundation,Inc.,59 Temple Place,Suite 330,Boston,MA 02111-1307 USA # # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com # @@ -28,7 +28,7 @@ import unittest from math import pi,e,sqrt,cos,sin from datetime import datetime from MEDCouplingDataForTest import MEDCouplingDataForTest -import rlcompleter,readline # this line has to be here, to ensure a usability of MEDCoupling/MEDLoader. B4 removing it please notify to anthony.geay@cea.fr +import rlcompleter,readline # this line has to be here,to ensure a usability of MEDCoupling/MEDLoader. B4 removing it please notify to anthony.geay@cea.fr class MEDCouplingBasicsTest7(unittest.TestCase): @@ -184,7 +184,7 @@ class MEDCouplingBasicsTest7(unittest.TestCase): def testDAIAggregateMulti1(self): a=DataArrayInt64.New() - a.setValues(list(range(4)), 2, 2) + a.setValues(list(range(4)),2, 2) a.setName("aa") b=DataArrayInt64.New() b.setValues(list(range(6)), 3, 2) @@ -785,6 +785,29 @@ class MEDCouplingBasicsTest7(unittest.TestCase): sk2 = sk.uniqueNotSortedByPack() self.assertTrue(sk2.getIndexArray().isEqual(DataArrayInt([0,3,8,13,16,21,29,37,42]))) self.assertTrue(sk2.getValuesArray().isEqual(DataArrayInt([1,4,5,0,2,4,5,6,1,3,5,6,7,2,6,7,0,1,5,8,9,0,1,2,4,6,8,9,10,1,2,3,5,7,9,10,11,2,3,6,10,11]))) + + def testSkyLineAggregatePacks1(self): + arr = DataArrayDouble(3) ; arr.iota() + m = MEDCouplingCMesh() ; m.setCoords(arr,arr) ; m = m.buildUnstructured() + a,b = m.computeEnlargedNeighborsOfNodes() + sk = MEDCouplingSkyLineArray(b,a) + sk1 = sk.deepCopy() + sk1.getValuesArray()[:] *= 2 + sk2 = sk.deepCopy() + sk2.getValuesArray()[:] *= 3 + skOut = MEDCouplingSkyLineArray.AggregatePacks([sk,sk1,sk2]) + self.assertTrue(skOut.getIndexArray().isEqual(DataArrayInt([0,9,24,33,48,72,87,96,111,120]))) + self.assertTrue(skOut.getValuesArray().isEqual(DataArrayInt([1,3,4,2,6,8,3,9,12,0,2,3,4,5,0,4,6,8,10,0,6,9,12,15,1,4,5,2,8,10,3,12,15,0,1,4,6,7,0,2,8,12,14,0,3,12,18,21,0,1,2,3,5,6,7,8,0,2,4,6,10,12,14,16,0,3,6,9,15,18,21,24,1,2,4,7,8,2,4,8,14,16,3,6,12,21,24,3,4,7,6,8,14,9,12,21,3,4,5,6,8,6,8,10,12,16,9,12,15,18,24,4,5,7,8,10,14,12,15,21]))) + + def testDACopySorted1(self): + d = DataArrayInt32([5,1,100,20]) + self.assertTrue(d.copySorted().isEqual(DataArrayInt32([1,5,20,100]))) + d = DataArrayInt64([5,1,100,20]) + self.assertTrue(d.copySorted().isEqual(DataArrayInt64([1,5,20,100]))) + d = DataArrayInt([5,1,100,20]) + self.assertTrue(d.copySorted().isEqual(DataArrayInt([1,5,20,100]))) + d = DataArrayDouble([5,1,100,20]) + self.assertTrue(d.copySorted().isEqual(DataArrayDouble([1,5,20,100]),1e-10)) pass diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index d2143962a..6fc12afea 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -468,6 +468,8 @@ typedef long int mcIdType; %newobject MEDCoupling::MEDCouplingSkyLineArray::getValuesArray; %newobject MEDCoupling::MEDCouplingSkyLineArray::groupPacks; %newobject MEDCoupling::MEDCouplingSkyLineArray::uniqueNotSortedByPack; +%newobject MEDCoupling::MEDCouplingSkyLineArray::AggregatePacks; +%newobject MEDCoupling::MEDCouplingSkyLineArray::deepCopy; %feature("unref") MEDCouplingPointSet "$this->decrRef();" %feature("unref") MEDCouplingMesh "$this->decrRef();" @@ -1301,6 +1303,8 @@ namespace MEDCoupling MEDCouplingSkyLineArray *groupPacks(const DataArrayIdType *indexedPacks) const; MEDCouplingSkyLineArray *uniqueNotSortedByPack() const; + + MEDCouplingSkyLineArray *deepCopy() const; %extend { @@ -1391,6 +1395,13 @@ namespace MEDCoupling convertFromPyObjVectorOfObj(listePacks,SWIGTITraits::TI,"DataArrayIdType",packs); self->replaceSimplePacks(idx, packs); } + + static MEDCouplingSkyLineArray *AggregatePacks(PyObject *sks) + { + std::vector sksCpp; + convertFromPyObjVectorOfObj(sks,SWIGTYPE_p_MEDCoupling__MEDCouplingSkyLineArray,"MEDCouplingSkyLineArray",sksCpp); + return MEDCoupling::MEDCouplingSkyLineArray::AggregatePacks(sksCpp); + } void replacePack(const mcIdType superIdx, const mcIdType idx, PyObject *pack) { diff --git a/src/MEDCoupling_Swig/MEDCouplingDataArrayTypemaps.i b/src/MEDCoupling_Swig/MEDCouplingDataArrayTypemaps.i index f6d8969d4..a3f52dec0 100644 --- a/src/MEDCoupling_Swig/MEDCouplingDataArrayTypemaps.i +++ b/src/MEDCoupling_Swig/MEDCouplingDataArrayTypemaps.i @@ -1389,6 +1389,17 @@ static void convertFromPyObjVectorOfObj(PyObject *pyLi, swig_type_info *ty, cons throw INTERP_KERNEL::Exception("convertFromPyObjVectorOfObj : not a list nor a tuple"); } +//convertFromVectorAutoObjToPyObj(inpv,SWIGTYPE_p_MEDCoupling__MEDCouplingUMesh) +template +static PyObject *convertFromVectorAutoObjToPyObj(std::vector< MCAuto >& inpVector, swig_type_info *ty) +{ + std::size_t sz(inpVector.size()); + PyObject *ret = PyList_New(sz); + for(std::size_t i=0;i cpp int sw=1 * if python list[int] -> cpp vector sw=2 diff --git a/src/MEDCoupling_Swig/MEDCouplingMemArray.i b/src/MEDCoupling_Swig/MEDCouplingMemArray.i index daf6d12a9..5325ef184 100644 --- a/src/MEDCoupling_Swig/MEDCouplingMemArray.i +++ b/src/MEDCoupling_Swig/MEDCouplingMemArray.i @@ -73,6 +73,7 @@ //$$$$$$$$$$$$$$$$$$ %newobject MEDCoupling::DataArray::deepCopy; +%newobject MEDCoupling::DataArray::copySorted; %newobject MEDCoupling::DataArray::buildNewEmptyInstance; %newobject MEDCoupling::DataArray::selectByTupleRanges; %newobject MEDCoupling::DataArray::selectByTupleId; @@ -493,6 +494,7 @@ typedef DataArrayInt64 DataArrayIdType; virtual std::size_t getNbOfElems() const; virtual std::size_t getNbOfElemAllocated() const; virtual DataArray *deepCopy() const; + virtual DataArray *copySorted() const; virtual DataArray *buildNewEmptyInstance() const; virtual DataArray *selectByTupleIdSafeSlice(int bg, int end2, int step) const; virtual void rearrange(int newNbOfCompo); diff --git a/src/MEDLoader/MEDFileBasis.hxx b/src/MEDLoader/MEDFileBasis.hxx index 24e82e90d..4ab7b043c 100644 --- a/src/MEDLoader/MEDFileBasis.hxx +++ b/src/MEDLoader/MEDFileBasis.hxx @@ -79,6 +79,7 @@ namespace MEDCoupling //DataArrayMedInt *buildNewEmptyInstance() const { return new DataArrayMedInt(); }//ko DataArray *buildNewEmptyInstance() const { if ( sizeof(med_int)==sizeof(int)) return DataArrayInt32::New(); return DataArrayInt64::New(); } public: + DataArrayMedInt *copySorted(bool asc=true) const override { MCAuto ret(this->deepCopy()); ret->sort(); return ret.retn(); } DataArray *selectByTupleId(const mcIdType *new2OldBg, const mcIdType *new2OldEnd) const { return this->mySelectByTupleId(new2OldBg,new2OldEnd); } DataArray *selectByTupleId(const DataArrayIdType& di) const { return this->mySelectByTupleId(di); } DataArray *selectByTupleIdSafe(const mcIdType *new2OldBg, const mcIdType *new2OldEnd) const { return this->mySelectByTupleIdSafe(new2OldBg,new2OldEnd); } diff --git a/src/ParaMEDMEM/CommInterface.cxx b/src/ParaMEDMEM/CommInterface.cxx index 9aaebdaec..6fda70ddd 100644 --- a/src/ParaMEDMEM/CommInterface.cxx +++ b/src/ParaMEDMEM/CommInterface.cxx @@ -19,10 +19,14 @@ #include "CommInterface.hxx" -#include - namespace MEDCoupling { + MPI_Datatype ParaTraits::MPIDataType = MPI_DOUBLE; + + MPI_Datatype ParaTraits::MPIDataType = MPI_INT; + + MPI_Datatype ParaTraits::MPIDataType = MPI_LONG; + /*! \anchor CommInterface-det \class CommInterface @@ -62,7 +66,25 @@ namespace MEDCoupling * Generalized AllGather collective communication. * This method send input \a array to all procs. */ - void CommInterface::allGatherArrays(MPI_Comm comm, const DataArrayIdType *array, std::unique_ptr& result, std::unique_ptr& resultIndex) const + void CommInterface::allGatherArrays(MPI_Comm comm, const DataArrayIdType *array, std::vector< MCAuto >& arraysOut) const + { + std::unique_ptr result, resultIndex; + int size(this->allGatherArrays(comm,array,result,resultIndex)); + arraysOut.resize(size); + for(int i = 0 ; i < size ; ++i) + { + arraysOut[i] = DataArrayIdType::New(); + mcIdType nbOfEltPack(resultIndex[i+1]-resultIndex[i]); + arraysOut[i]->alloc(nbOfEltPack,1); + std::copy(result.get()+resultIndex[i],result.get()+resultIndex[i+1],arraysOut[i]->getPointer()); + } + } + + /*! + * Generalized AllGather collective communication. + * This method send input \a array to all procs. + */ + int CommInterface::allGatherArrays(MPI_Comm comm, const DataArrayIdType *array, std::unique_ptr& result, std::unique_ptr& resultIndex) const { int size; this->commSize(comm,&size); @@ -74,7 +96,20 @@ namespace MEDCoupling std::unique_ptr nbOfElemsInt( CommInterface::ToIntArray(nbOfElems,size) ); std::unique_ptr offsetsIn( CommInterface::ComputeOffset(nbOfElemsInt,size) ); this->allGatherV(array->begin(),nbOfCellsRequested,MPI_ID_TYPE,result.get(),nbOfElemsInt.get(),offsetsIn.get(),MPI_ID_TYPE,comm); - resultIndex = std::move(nbOfElems); + resultIndex = ComputeOffsetFull(nbOfElems,size); + return size; + } + + void CommInterface::allToAllArrays(MPI_Comm comm, const std::vector< MCAuto >& arrays, std::vector< MCAuto >& arraysOut) const + { + this->allToAllArraysT(comm,arrays,arraysOut); + } + + void CommInterface::allToAllArrays(MPI_Comm comm, const std::vector< MCAuto >& arrays, MCAuto& arraysOut) const + { + std::unique_ptr notUsed1; + mcIdType notUsed2; + this->allToAllArraysT2(comm,arrays,arraysOut,notUsed1,notUsed2); } /*! @@ -82,34 +117,7 @@ namespace MEDCoupling */ void CommInterface::allToAllArrays(MPI_Comm comm, const std::vector< MCAuto >& arrays, std::vector< MCAuto >& arraysOut) const { - int size; - this->commSize(comm,&size); - if( arrays.size() != ToSizeT(size) ) - throw INTERP_KERNEL::Exception("AllToAllArrays : internal error ! Invalid size of input array."); - - std::vector< const DataArrayIdType *> arraysBis(FromVecAutoToVecOfConst(arrays)); - std::unique_ptr nbOfElems2(new mcIdType[size]),nbOfElems3(new mcIdType[size]); - for(int curRk = 0 ; curRk < size ; ++curRk) - { - nbOfElems3[curRk] = arrays[curRk]->getNumberOfTuples(); - } - this->allToAll(nbOfElems3.get(),1,MPI_ID_TYPE,nbOfElems2.get(),1,MPI_ID_TYPE,comm); - mcIdType nbOfCellIdsSum(std::accumulate(nbOfElems2.get(),nbOfElems2.get()+size,0)); - MCAuto cellIdsFromProcs(DataArrayIdType::New()); - cellIdsFromProcs->alloc(nbOfCellIdsSum,1); - std::unique_ptr nbOfElemsInt( CommInterface::ToIntArray(nbOfElems3,size) ),nbOfElemsOutInt( CommInterface::ToIntArray(nbOfElems2,size) ); - std::unique_ptr offsetsIn( CommInterface::ComputeOffset(nbOfElemsInt,size) ), offsetsOut( CommInterface::ComputeOffset(nbOfElemsOutInt,size) ); - { - MCAuto arraysAcc(DataArrayIdType::Aggregate(arraysBis)); - this->allToAllV(arraysAcc->begin(),nbOfElemsInt.get(),offsetsIn.get(),MPI_ID_TYPE, - cellIdsFromProcs->getPointer(),nbOfElemsOutInt.get(),offsetsOut.get(),MPI_ID_TYPE,comm); - } - std::unique_ptr offsetsOutIdType( CommInterface::ComputeOffset(nbOfElems2,size) ); - // build output arraysOut by spliting cellIdsFromProcs into parts - arraysOut.resize(size); - for(int curRk = 0 ; curRk < size ; ++curRk) - { - arraysOut[curRk] = DataArrayIdType::NewFromArray(cellIdsFromProcs->begin()+offsetsOutIdType[curRk],cellIdsFromProcs->begin()+offsetsOutIdType[curRk]+nbOfElems2[curRk]); - } + this->allToAllArraysT(comm,arrays,arraysOut); } + } diff --git a/src/ParaMEDMEM/CommInterface.hxx b/src/ParaMEDMEM/CommInterface.hxx index d132d09db..e969f9394 100644 --- a/src/ParaMEDMEM/CommInterface.hxx +++ b/src/ParaMEDMEM/CommInterface.hxx @@ -25,9 +25,33 @@ #include #include +#include namespace MEDCoupling { + template + struct ParaTraits + { + using EltType = T; + }; + + template<> + struct ParaTraits + { + static MPI_Datatype MPIDataType; + }; + + template<> + struct ParaTraits + { + static MPI_Datatype MPIDataType; + }; + + template<> + struct ParaTraits + { + static MPI_Datatype MPIDataType; + }; class CommInterface { @@ -94,8 +118,69 @@ namespace MEDCoupling MPI_Op op, int root, MPI_Comm comm) const { return MPI_Reduce(sendbuf, recvbuf, count, datatype, op, root, comm); } int allReduce(void* sendbuf, void* recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm) const { return MPI_Allreduce(sendbuf, recvbuf, count, datatype, op, comm); } public: - void allGatherArrays(MPI_Comm comm, const DataArrayIdType *array, std::unique_ptr& result, std::unique_ptr& resultIndex) const; + void allGatherArrays(MPI_Comm comm, const DataArrayIdType *array, std::vector< MCAuto >& arraysOut) const; + int allGatherArrays(MPI_Comm comm, const DataArrayIdType *array, std::unique_ptr& result, std::unique_ptr& resultIndex) const; void allToAllArrays(MPI_Comm comm, const std::vector< MCAuto >& arrays, std::vector< MCAuto >& arraysOut) const; + void allToAllArrays(MPI_Comm comm, const std::vector< MCAuto >& arrays, std::vector< MCAuto >& arraysOut) const; + void allToAllArrays(MPI_Comm comm, const std::vector< MCAuto >& arrays, MCAuto& arraysOut) const; + template + int allToAllArraysT2(MPI_Comm comm, const std::vector< MCAuto::ArrayType> >& arrays, MCAuto::ArrayType>& arrayOut, std::unique_ptr& nbOfElems2, mcIdType& nbOfComponents) const + { + using DataArrayT = typename Traits::ArrayType; + int size; + this->commSize(comm,&size); + if( arrays.size() != ToSizeT(size) ) + throw INTERP_KERNEL::Exception("AllToAllArrays : internal error ! Invalid size of input array."); + + std::vector< const DataArrayT *> arraysBis(FromVecAutoToVecOfConst(arrays)); + std::unique_ptr nbOfElems3(new mcIdType[size]); + nbOfElems2.reset(new mcIdType[size]); + nbOfComponents = std::numeric_limits::max(); + for(int curRk = 0 ; curRk < size ; ++curRk) + { + mcIdType curNbOfCompo( ToIdType( arrays[curRk]->getNumberOfComponents() ) ); + if(nbOfComponents != std::numeric_limits::max()) + { + if( nbOfComponents != curNbOfCompo ) + throw INTERP_KERNEL::Exception("AllToAllArrays : internal error ! Nb of components is not homogeneous !"); + } + else + { + nbOfComponents = curNbOfCompo; + } + nbOfElems3[curRk] = arrays[curRk]->getNbOfElems(); + } + this->allToAll(nbOfElems3.get(),1,MPI_ID_TYPE,nbOfElems2.get(),1,MPI_ID_TYPE,comm); + mcIdType nbOfCellIdsSum(std::accumulate(nbOfElems2.get(),nbOfElems2.get()+size,0)); + arrayOut = DataArrayT::New(); + arrayOut->alloc(nbOfCellIdsSum/nbOfComponents,nbOfComponents); + std::unique_ptr nbOfElemsInt( CommInterface::ToIntArray(nbOfElems3,size) ),nbOfElemsOutInt( CommInterface::ToIntArray(nbOfElems2,size) ); + std::unique_ptr offsetsIn( CommInterface::ComputeOffset(nbOfElemsInt,size) ), offsetsOut( CommInterface::ComputeOffset(nbOfElemsOutInt,size) ); + { + MCAuto arraysAcc(DataArrayT::Aggregate(arraysBis)); + this->allToAllV(arraysAcc->begin(),nbOfElemsInt.get(),offsetsIn.get(),ParaTraits::MPIDataType, + arrayOut->getPointer(),nbOfElemsOutInt.get(),offsetsOut.get(),ParaTraits::MPIDataType,comm); + } + return size; + } + + template + void allToAllArraysT(MPI_Comm comm, const std::vector< MCAuto::ArrayType> >& arrays, std::vector< MCAuto::ArrayType> >& arraysOut) const + { + using DataArrayT = typename Traits::ArrayType; + MCAuto cellIdsFromProcs; + std::unique_ptr nbOfElems2; + mcIdType nbOfComponents(0); + int size(this->allToAllArraysT2(comm,arrays,cellIdsFromProcs,nbOfElems2,nbOfComponents)); + std::unique_ptr offsetsOutIdType( CommInterface::ComputeOffset(nbOfElems2,size) ); + // build output arraysOut by spliting cellIdsFromProcs into parts + arraysOut.resize(size); + for(int curRk = 0 ; curRk < size ; ++curRk) + { + arraysOut[curRk] = DataArrayT::NewFromArray(cellIdsFromProcs->begin()+offsetsOutIdType[curRk],cellIdsFromProcs->begin()+offsetsOutIdType[curRk]+nbOfElems2[curRk]); + arraysOut[curRk]->rearrange(nbOfComponents); + } + } public: /*! @@ -138,5 +223,20 @@ namespace MEDCoupling } return ret; } + + /*! + * Helper of alltoallv and allgatherv + */ + template + static std::unique_ptr ComputeOffsetFull(const std::unique_ptr& counts, std::size_t sizeOfCounts) + { + std::unique_ptr ret(new T[sizeOfCounts+1]); + ret[0] = static_cast(0); + for(std::size_t i = 1 ; i < sizeOfCounts+1 ; ++i) + { + ret[i] = ret[i-1] + counts[i-1]; + } + return ret; + } }; } diff --git a/src/ParaMEDMEM/ParaUMesh.cxx b/src/ParaMEDMEM/ParaUMesh.cxx index 812a7e4d8..60d44a323 100644 --- a/src/ParaMEDMEM/ParaUMesh.cxx +++ b/src/ParaMEDMEM/ParaUMesh.cxx @@ -37,6 +37,11 @@ using namespace MEDCoupling; +ParaUMesh *ParaUMesh::New(MEDCouplingUMesh *mesh, DataArrayIdType *globalCellIds, DataArrayIdType *globalNodeIds) +{ + return new ParaUMesh(mesh,globalCellIds,globalNodeIds); +} + ParaUMesh::ParaUMesh(MEDCouplingUMesh *mesh, DataArrayIdType *globalCellIds, DataArrayIdType *globalNodeIds) { _mesh.takeRef(mesh); @@ -52,6 +57,16 @@ ParaUMesh::ParaUMesh(MEDCouplingUMesh *mesh, DataArrayIdType *globalCellIds, Dat throw INTERP_KERNEL::Exception("ParaUMesh constructor : mismatch between # cells and len of global # cells."); } +std::size_t ParaUMesh::getHeapMemorySizeWithoutChildren() const +{ + return 0; +} + +std::vector ParaUMesh::getDirectChildrenWithNull() const +{ + return {_mesh,_cell_global,_node_global}; +} + /*! * This method computes the cells part of distributed mesh lying on \a globalNodeIds nodes. * The input \a globalNodeIds are not supposed to reside on the current process. @@ -117,33 +132,25 @@ MCAuto ParaUMesh::getCellIdsLyingOnNodes(const DataArrayIdType } /*! + * Return part of \a this mesh split over COMM_WORLD. Part is defined by global cell ids array \a globaCellIds. */ -MCAuto ParaUMesh::redistributeCells(const DataArrayIdType *globalCellIds) const +ParaUMesh *ParaUMesh::redistributeCells(const DataArrayIdType *globalCellIds) const { MPI_Comm comm(MPI_COMM_WORLD); CommInterface ci; - int size; - ci.commSize(comm,&size); - std::unique_ptr nbOfElems(new mcIdType[size]); - mcIdType nbOfCellsRequested(globalCellIds->getNumberOfTuples()); - ci.allGather(&nbOfCellsRequested,1,MPI_ID_TYPE,nbOfElems.get(),1,MPI_ID_TYPE,comm); - mcIdType nbOfCellIdsSum(std::accumulate(nbOfElems.get(),nbOfElems.get()+size,0)); - std::unique_ptr allGlobalCellIds(new mcIdType[nbOfCellIdsSum]); - std::unique_ptr nbOfElemsInt( CommInterface::ToIntArray(nbOfElems,size) ); - std::unique_ptr offsetsIn( CommInterface::ComputeOffset(nbOfElemsInt,size) ); - ci.allGatherV(globalCellIds->begin(),nbOfCellsRequested,MPI_ID_TYPE,allGlobalCellIds.get(),nbOfElemsInt.get(),offsetsIn.get(),MPI_ID_TYPE,comm); - mcIdType offset(0); + std::unique_ptr allGlobalCellIds,allGlobalCellIdsIndex; + int size(ci.allGatherArrays(comm,globalCellIds,allGlobalCellIds,allGlobalCellIdsIndex)); // Prepare ParaUMesh parts to be sent : compute for each proc the contribution of current rank. std::vector< MCAuto > globalCellIdsToBeSent(size),globalNodeIdsToBeSent(size); std::vector< MCAuto > meshPartsToBeSent(size); for(int curRk = 0 ; curRk < size ; ++curRk) { + mcIdType offset(allGlobalCellIdsIndex[curRk]); MCAuto globalCellIdsOfCurProc(DataArrayIdType::New()); - globalCellIdsOfCurProc->useArray(allGlobalCellIds.get()+offset,false,DeallocType::CPP_DEALLOC,nbOfElems[curRk],1); - offset += nbOfElems[curRk]; + globalCellIdsOfCurProc->useArray(allGlobalCellIds.get()+offset,false,DeallocType::CPP_DEALLOC,allGlobalCellIdsIndex[curRk+1]-offset,1); // the key call is here : compute for rank curRk the cells to be sent MCAuto globalCellIdsCaptured(_cell_global->buildIntersection(globalCellIdsOfCurProc));// OK for the global cellIds - MCAuto localCellIdsCaptured(_node_global->findIdForEach(globalCellIdsCaptured->begin(),globalCellIdsCaptured->end())); + MCAuto localCellIdsCaptured(_cell_global->findIdForEach(globalCellIdsCaptured->begin(),globalCellIdsCaptured->end())); MCAuto meshPart(_mesh->buildPartOfMySelf(localCellIdsCaptured->begin(),localCellIdsCaptured->end(),true)); MCAuto o2n(meshPart->zipCoordsTraducer());// OK for the mesh MCAuto n2o(o2n->invertArrayO2N2N2O(meshPart->getNumberOfNodes())); @@ -152,5 +159,74 @@ MCAuto ParaUMesh::redistributeCells(const DataArrayIdType *globalCell globalCellIdsToBeSent[curRk] = globalCellIdsCaptured; globalNodeIdsToBeSent[curRk] = globalNodeIdsPart; } - // Receive + // Receive + std::vector< MCAuto > globalCellIdsReceived,globalNodeIdsReceived; + ci.allToAllArrays(comm,globalCellIdsToBeSent,globalCellIdsReceived); + ci.allToAllArrays(comm,globalNodeIdsToBeSent,globalNodeIdsReceived); + //now exchange the 3 arrays for the umesh : connectivity, connectivityindex and coordinates + std::vector meshPartsToBeSent2(FromVecAutoToVecOfConst(meshPartsToBeSent)); + //connectivityindex + std::vector< MCAuto > connectivityIndexReceived,connectivityReceived; + { + std::vector connectivityIndexToBeSent(UMeshConnectivityIndexIterator(0,&meshPartsToBeSent2),UMeshConnectivityIndexIterator(meshPartsToBeSent2.size(),&meshPartsToBeSent2)); + ci.allToAllArrays(comm,FromVecConstToVecAuto(connectivityIndexToBeSent),connectivityIndexReceived); + } + //connectivity + { + std::vector connectivityToBeSent(UMeshConnectivityIterator(0,&meshPartsToBeSent2),UMeshConnectivityIterator(meshPartsToBeSent2.size(),&meshPartsToBeSent2)); + ci.allToAllArrays(comm,FromVecConstToVecAuto(connectivityToBeSent),connectivityReceived); + } + //coordinates + MCAuto coords; + { + std::vector coordsToBeSent(UMeshCoordsIterator(0,&meshPartsToBeSent2),UMeshCoordsIterator(meshPartsToBeSent2.size(),&meshPartsToBeSent2)); + ci.allToAllArrays(comm,FromVecConstToVecAuto(coordsToBeSent),coords); + } + /////// Sort it all ! + // firstly deal with nodes. + MCAuto aggregatedNodeIds( DataArrayIdType::Aggregate(FromVecAutoToVecOfConst(globalNodeIdsReceived)) ); + MCAuto aggregatedNodeIdsSorted(aggregatedNodeIds->copySorted()); + MCAuto nodeIdsIntoAggregatedIds(DataArrayIdType::FindPermutationFromFirstToSecondDuplicate(aggregatedNodeIdsSorted,aggregatedNodeIds)); + MCAuto idxOfSameNodeIds(aggregatedNodeIdsSorted->indexOfSameConsecutiveValueGroups()); + MCAuto n2o_nodes(nodeIdsIntoAggregatedIds->selectByTupleIdSafe(idxOfSameNodeIds->begin(),idxOfSameNodeIds->end()-1));//new == new ordering so that global node ids are sorted . old == coarse ordering implied by the aggregation + MCAuto finalGlobalNodeIds(aggregatedNodeIdsSorted->selectByTupleIdSafe(idxOfSameNodeIds->begin(),idxOfSameNodeIds->end()-1)); + MCAuto finalCoords(coords->selectByTupleIdSafe(n2o_nodes->begin(),n2o_nodes->end())); + finalCoords->copyStringInfoFrom(*_mesh->getCoords()); + // secondly renumbering of node ids in connectivityReceived + for(int curRk = 0 ; curRk < size ; ++curRk) + { + auto current(globalNodeIdsReceived[curRk]); + MCAuto aa(finalGlobalNodeIds->findIdForEach(current->begin(),current->end())); + // work on connectivityReceived[curRk] with transformWithIndArr but do not forget type of cells that should be excluded ! + auto connectivityToModify(connectivityReceived[curRk]); + auto connectivityIndex(connectivityIndexReceived[curRk]); + MCAuto types(connectivityToModify->selectByTupleIdSafe(connectivityIndex->begin(),connectivityIndex->end()-1)); + connectivityToModify->setPartOfValuesSimple3(0,connectivityIndex->begin(),connectivityIndex->end()-1,0,1,1); + connectivityToModify->transformWithIndArr(aa->begin(),aa->end()); + connectivityToModify->setPartOfValues3(types,connectivityIndex->begin(),connectivityIndex->end()-1,0,1,1,true); + } + // thirdly renumber cells + MCAuto aggregatedCellIds( DataArrayIdType::Aggregate(FromVecAutoToVecOfConst(globalCellIdsReceived)) ); + MCAuto aggregatedCellIdsSorted(aggregatedCellIds->copySorted()); + MCAuto idsIntoAggregatedIds(DataArrayIdType::FindPermutationFromFirstToSecondDuplicate(aggregatedCellIdsSorted,aggregatedCellIds)); + MCAuto cellIdsOfSameNodeIds(aggregatedCellIdsSorted->indexOfSameConsecutiveValueGroups()); + MCAuto n2o_cells(idsIntoAggregatedIds->selectByTupleIdSafe(cellIdsOfSameNodeIds->begin(),cellIdsOfSameNodeIds->end()-1));//new == new ordering so that global cell ids are sorted . old == coarse ordering implied by the aggregation + // TODO : check coordsReceived==globalCellIds + MCAuto connSorted,indicesSorted; + { + MCAuto conn(DataArrayIdType::Aggregate(FromVecAutoToVecOfConst(connectivityReceived))); + MCAuto connIndex(DataArrayIdType::AggregateIndexes(FromVecAutoToVecOfConst(connectivityIndexReceived))); + { + DataArrayIdType *indicesSortedTmp(nullptr),*valuesSortedTmp(nullptr); + DataArrayIdType::ExtractFromIndexedArrays(n2o_cells->begin(),n2o_cells->end(),conn,connIndex,valuesSortedTmp,indicesSortedTmp); + indicesSorted = indicesSortedTmp; connSorted=valuesSortedTmp; + } + } + // finalize all + MCAuto mesh(MEDCouplingUMesh::New(_mesh->getName(),_mesh->getMeshDimension())); + mesh->setConnectivity(connSorted,indicesSorted,true); + mesh->setCoords(finalCoords); + mesh->setDescription(_mesh->getDescription()); + MCAuto ret(ParaUMesh::New(mesh,aggregatedCellIdsSorted,finalGlobalNodeIds)); + return ret.retn(); } diff --git a/src/ParaMEDMEM/ParaUMesh.hxx b/src/ParaMEDMEM/ParaUMesh.hxx index 3d0bff392..5247b3716 100644 --- a/src/ParaMEDMEM/ParaUMesh.hxx +++ b/src/ParaMEDMEM/ParaUMesh.hxx @@ -34,13 +34,21 @@ namespace MEDCoupling * * This class is very specific to the requirement of parallel code computations. */ - class ParaUMesh + class ParaUMesh : public RefCountObject { public: - ParaUMesh(MEDCouplingUMesh *mesh, DataArrayIdType *globalCellIds, DataArrayIdType *globalNodeIds); + static ParaUMesh *New(MEDCouplingUMesh *mesh, DataArrayIdType *globalCellIds, DataArrayIdType *globalNodeIds); MCAuto getCellIdsLyingOnNodes(const DataArrayIdType *globalNodeIds, bool fullyIn) const; - MCAuto redistributeCells(const DataArrayIdType *globalCellIds) const; + ParaUMesh *redistributeCells(const DataArrayIdType *globalCellIds) const; + MEDCouplingUMesh *getMesh() { return _mesh; } + DataArrayIdType *getGlobalCellIds() { return _cell_global; } + DataArrayIdType *getGlobalNodeIds() { return _node_global; } + protected: virtual ~ParaUMesh() { } + ParaUMesh(MEDCouplingUMesh *mesh, DataArrayIdType *globalCellIds, DataArrayIdType *globalNodeIds); + std::string getClassName() const override { return "ParaUMesh"; } + std::size_t getHeapMemorySizeWithoutChildren() const override; + std::vector getDirectChildrenWithNull() const override; private: MCAuto _mesh; MCAuto _cell_global; diff --git a/src/ParaMEDMEM_Swig/ParaMEDMEMCommon.i b/src/ParaMEDMEM_Swig/ParaMEDMEMCommon.i index 452834e76..b1fcf0103 100644 --- a/src/ParaMEDMEM_Swig/ParaMEDMEMCommon.i +++ b/src/ParaMEDMEM_Swig/ParaMEDMEMCommon.i @@ -58,12 +58,19 @@ using namespace ICoCo; %rename(ICoCoMEDField) ICoCo::MEDField; %include "ICoCoMEDField.hxx" +%newobject MEDCoupling::ParaUMesh::New; +%newobject MEDCoupling::ParaUMesh::getMesh; +%newobject MEDCoupling::ParaUMesh::getGlobalCellIds; +%newobject MEDCoupling::ParaUMesh::getGlobalNodeIds; %newobject MEDCoupling::ParaUMesh::getCellIdsLyingOnNodes; +%newobject MEDCoupling::ParaUMesh::redistributeCells; +%newobject MEDCoupling::ParaSkyLineArray::New; %newobject MEDCoupling::ParaSkyLineArray::equiRedistribute; %newobject MEDCoupling::ParaSkyLineArray::getSkyLineArray; %newobject MEDCoupling::ParaSkyLineArray::getGlobalIdsArray; %feature("unref") ParaSkyLineArray "$this->decrRef();" +%feature("unref") ParaUMesh "$this->decrRef();" %nodefaultctor; @@ -127,14 +134,60 @@ namespace MEDCoupling int reduce(void* sendbuf, void* recvbuf, int count, MPI_Datatype datatype, MPI_Op op, int root, MPI_Comm comm) const; int allReduce(void* sendbuf, void* recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm) const; + %extend + { + PyObject *allGatherArrays(const DataArrayIdType *array) const + { + std::vector< MCAuto > ret; + self->allGatherArrays(MPI_COMM_WORLD,array,ret); + return convertFromVectorAutoObjToPyObj(ret,SWIGTITraits::TI); + } + + PyObject *allToAllArrays(PyObject *arrays) const + { + std::vector< DataArrayIdType * > arraysIn; + std::vector< MCAuto > arrayOut; + convertFromPyObjVectorOfObj(arrays,SWIGTITraits::TI,"DataArrayIdType",arraysIn); + std::vector< MCAuto > arraysIn2(FromVecToVecAuto(arraysIn)); + self->allToAllArrays(MPI_COMM_WORLD,arraysIn2,arrayOut); + return convertFromVectorAutoObjToPyObj(arrayOut,SWIGTITraits::TI); + } + } }; - class ParaUMesh + class ParaUMesh : public RefCountObject { public: - ParaUMesh(MEDCouplingUMesh *mesh, DataArrayIdType *globalCellIds, DataArrayIdType *globalNodeIds); + static ParaUMesh *New(MEDCouplingUMesh *mesh, DataArrayIdType *globalCellIds, DataArrayIdType *globalNodeIds); + ParaUMesh *redistributeCells(const DataArrayIdType *globalCellIds) const; %extend { + ParaUMesh(MEDCouplingUMesh *mesh, DataArrayIdType *globalCellIds, DataArrayIdType *globalNodeIds) + { + return ParaUMesh::New(mesh,globalCellIds,globalNodeIds); + } + + MEDCouplingUMesh *getMesh() + { + MEDCouplingUMesh *ret(self->getMesh()); + if(ret) ret->incrRef(); + return ret; + } + + DataArrayIdType *getGlobalCellIds() + { + DataArrayIdType *ret(self->getGlobalCellIds()); + if(ret) ret->incrRef(); + return ret; + } + + DataArrayIdType *getGlobalNodeIds() + { + DataArrayIdType *ret(self->getGlobalNodeIds()); + if(ret) ret->incrRef(); + return ret; + } + DataArrayIdType *getCellIdsLyingOnNodes(const DataArrayIdType *globalNodeIds, bool fullyIn) const { MCAuto ret(self->getCellIdsLyingOnNodes(globalNodeIds,fullyIn)); -- 2.39.2