From: Anthony Geay Date: Tue, 27 Sep 2016 13:34:09 +0000 (+0200) Subject: Implementation of planar symmetry + aggregation methods for MEDFileData+MEDFileUMesh X-Git-Tag: V8_2_0a1~22 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=cdd09520be1ff9d51b7f67e39fb0866bb71db901;p=tools%2Fmedcoupling.git Implementation of planar symmetry + aggregation methods for MEDFileData+MEDFileUMesh --- diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index f8209fd68..2c3cb6501 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -946,6 +946,24 @@ namespace INTERP_KERNEL } } + /*! + * Find a vector orthogonal to the input vector + */ + inline void orthogonalVect3(const double inpVect[3], double outVect[3]) + { + std::vector sw(3,false); + double inpVect2[3]; + std::transform(inpVect,inpVect+3,inpVect2,std::ptr_fun(fabs)); + std::size_t posMin(std::distance(inpVect2,std::min_element(inpVect2,inpVect2+3))); + sw[posMin]=true; + std::size_t posMax(std::distance(inpVect2,std::max_element(inpVect2,inpVect2+3))); + if(posMax==posMin) + { posMax=(posMin+1)%3; } + sw[posMax]=true; + std::size_t posMid(std::distance(sw.begin(),std::find(sw.begin(),sw.end(),false))); + outVect[posMin]=0.; outVect[posMid]=1.; outVect[posMax]=-inpVect[posMid]/inpVect[posMax]; + } + /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ /* Computes the dot product of a and b */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ diff --git a/src/MEDCoupling/MEDCouplingCurveLinearMesh.cxx b/src/MEDCoupling/MEDCouplingCurveLinearMesh.cxx index d5621e903..d16420e8b 100644 --- a/src/MEDCoupling/MEDCouplingCurveLinearMesh.cxx +++ b/src/MEDCoupling/MEDCouplingCurveLinearMesh.cxx @@ -637,9 +637,9 @@ void MEDCouplingCurveLinearMesh::rotate(const double *center, const double *vect int nbNodes=_coords->getNumberOfTuples(); double *coords=_coords->getPointer(); if(spaceDim==3) - MEDCouplingPointSet::Rotate3DAlg(center,vector,angle,nbNodes,coords); + DataArrayDouble::Rotate3DAlg(center,vector,angle,nbNodes,coords,coords); else if(spaceDim==2) - MEDCouplingPointSet::Rotate2DAlg(center,angle,nbNodes,coords); + DataArrayDouble::Rotate2DAlg(center,angle,nbNodes,coords,coords); else throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::rotate : invalid space dim for rotation must be 2 or 3"); _coords->declareAsNew(); diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index 35853e0a2..194cdccf4 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -19,7 +19,6 @@ // Author : Anthony Geay (CEA/DEN) #include "MEDCouplingMemArray.txx" -#include "MCAuto.hxx" #include "BBTree.txx" #include "GenMathFormulae.hxx" @@ -3217,6 +3216,21 @@ void DataArrayDouble::applyFuncFast64(const std::string& func) declareAsNew(); } +/*! + * \return a new object that is the result of the symmetry along 3D plane defined by its normal vector \a normalVector and a point \a point. + */ +MCAuto DataArrayDouble::symmetry3DPlane(const double point[3], const double normalVector[3]) const +{ + checkAllocated(); + if(getNumberOfComponents()!=3) + throw INTERP_KERNEL::Exception("DataArrayDouble::symmetry3DPlane : this is excepted to have 3 components !"); + int nbTuples(getNumberOfTuples()); + MCAuto ret(DataArrayDouble::New()); + ret->alloc(nbTuples,3); + Symmetry3DPlane(point,normalVector,nbTuples,begin(),ret->getPointer()); + return ret; +} + DataArrayDoubleIterator *DataArrayDouble::iterator() { return new DataArrayDoubleIterator(this); @@ -4356,6 +4370,112 @@ void DataArrayDouble::finishUnserialization(const std::vector& tinyInfoI, c } } +/*! + * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in \a coordsIn + * around an axe ( \a center, \a vect) and with angle \a angle. + */ +void DataArrayDouble::Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, const double *coordsIn, double *coordsOut) +{ + if(!center || !vect) + throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : null vector in input !"); + double sina(sin(angle)); + double cosa(cos(angle)); + double vectorNorm[3]; + double matrix[9]; + double matrixTmp[9]; + double norm(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2])); + if(norm::min()) + throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : magnitude of input vector is too close of 0. !"); + std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies(),1/norm)); + //rotation matrix computation + matrix[0]=cosa; matrix[1]=0.; matrix[2]=0.; matrix[3]=0.; matrix[4]=cosa; matrix[5]=0.; matrix[6]=0.; matrix[7]=0.; matrix[8]=cosa; + matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2]; + matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2]; + matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2]; + std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies(),1-cosa)); + std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus()); + matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1]; + matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0]; + matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.; + std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies(),sina)); + std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus()); + //rotation matrix computed. + double tmp[3]; + for(int i=0; i()); + coordsOut[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0]; + coordsOut[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1]; + coordsOut[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2]; + } +} + +void DataArrayDouble::Symmetry3DPlane(const double point[3], const double normalVector[3], int nbNodes, const double *coordsIn, double *coordsOut) +{ + double matrix[9],matrix2[9],matrix3[9]; + double vect[3],crossVect[3]; + INTERP_KERNEL::orthogonalVect3(normalVector,vect); + crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1]; + crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2]; + crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0]; + double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect)); + matrix[0]=vect[0]/nv; matrix[1]=crossVect[0]/nc; matrix[2]=-normalVector[0]/ni; + matrix[3]=vect[1]/nv; matrix[4]=crossVect[1]/nc; matrix[5]=-normalVector[1]/ni; + matrix[6]=vect[2]/nv; matrix[7]=crossVect[2]/nc; matrix[8]=-normalVector[2]/ni; + matrix2[0]=vect[0]/nv; matrix2[1]=vect[1]/nv; matrix2[2]=vect[2]/nv; + matrix2[3]=crossVect[0]/nc; matrix2[4]=crossVect[1]/nc; matrix2[5]=crossVect[2]/nc; + matrix2[6]=normalVector[0]/ni; matrix2[7]=normalVector[1]/ni; matrix2[8]=normalVector[2]/ni; + for(int i=0;i<3;i++) + for(int j=0;j<3;j++) + { + double val(0.); + for(int k=0;k<3;k++) + val+=matrix[3*i+k]*matrix2[3*k+j]; + matrix3[3*i+j]=val; + } + //rotation matrix computed. + double tmp[3]; + for(int i=0; i()); + coordsOut[i*3]=matrix3[0]*tmp[0]+matrix3[1]*tmp[1]+matrix3[2]*tmp[2]+point[0]; + coordsOut[i*3+1]=matrix3[3]*tmp[0]+matrix3[4]*tmp[1]+matrix3[5]*tmp[2]+point[1]; + coordsOut[i*3+2]=matrix3[6]*tmp[0]+matrix3[7]*tmp[1]+matrix3[8]*tmp[2]+point[2]; + } +} + +void DataArrayDouble::GiveBaseForPlane(const double normalVector[3], double baseOfPlane[9]) +{ + double vect[3],crossVect[3]; + INTERP_KERNEL::orthogonalVect3(normalVector,vect); + crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1]; + crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2]; + crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0]; + double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect)); + baseOfPlane[0]=vect[0]/nv; baseOfPlane[1]=vect[1]/nv; baseOfPlane[2]=vect[2]/nv; + baseOfPlane[3]=crossVect[0]/nc; baseOfPlane[4]=crossVect[1]/nc; baseOfPlane[5]=crossVect[2]/nc; + baseOfPlane[6]=normalVector[0]/ni; baseOfPlane[7]=normalVector[1]/ni; baseOfPlane[8]=normalVector[2]/ni; +} + +/*! + * Low static method that operates 3D rotation of \a nbNodes 3D nodes whose coordinates are arranged in \a coords + * around the center point \a center and with angle \a angle. + */ +void DataArrayDouble::Rotate2DAlg(const double *center, double angle, int nbNodes, const double *coordsIn, double *coordsOut) +{ + double cosa=cos(angle); + double sina=sin(angle); + double matrix[4]; + matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa; + double tmp[2]; + for(int i=0; i()); + coordsOut[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0]; + coordsOut[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1]; + } +} + DataArrayDoubleIterator::DataArrayDoubleIterator(DataArrayDouble *da):_da(da),_tuple_id(0),_nb_comp(0),_nb_tuple(0) { if(_da) diff --git a/src/MEDCoupling/MEDCouplingMemArray.hxx b/src/MEDCoupling/MEDCouplingMemArray.hxx index 77f89b9c7..991673bf9 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.hxx +++ b/src/MEDCoupling/MEDCouplingMemArray.hxx @@ -22,6 +22,7 @@ #define __MEDCOUPLING_MEDCOUPLINGMEMARRAY_HXX__ #include "MEDCoupling.hxx" +#include "MCAuto.hxx" #include "MEDCouplingTimeLabel.hxx" #include "MEDCouplingRefCountObject.hxx" #include "InterpKernelException.hxx" @@ -385,6 +386,7 @@ namespace MEDCoupling MEDCOUPLING_EXPORT DataArrayDouble *applyFuncNamedCompo(int nbOfComp, const std::vector& varsOrder, const std::string& func, bool isSafe=true) const; MEDCOUPLING_EXPORT void applyFuncFast32(const std::string& func); MEDCOUPLING_EXPORT void applyFuncFast64(const std::string& func); + MEDCOUPLING_EXPORT MCAuto symmetry3DPlane(const double point[3], const double normalVector[3]) const; MEDCOUPLING_EXPORT DataArrayInt *findIdsInRange(double vmin, double vmax) const; MEDCOUPLING_EXPORT DataArrayInt *findIdsNotInRange(double vmin, double vmax) const; MEDCOUPLING_EXPORT static DataArrayDouble *Aggregate(const DataArrayDouble *a1, const DataArrayDouble *a2); @@ -409,6 +411,10 @@ namespace MEDCoupling MEDCOUPLING_EXPORT MemArray& accessToMemArray() { return _mem; } MEDCOUPLING_EXPORT const MemArray& accessToMemArray() const { return _mem; } MEDCOUPLING_EXPORT std::vector toVectorOfBool(double eps) const; + MEDCOUPLING_EXPORT static void Rotate2DAlg(const double *center, double angle, int nbNodes, const double *coordsIn, double *coordsOut); + MEDCOUPLING_EXPORT static void Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, const double *coordsIn, double *coordsOut); + MEDCOUPLING_EXPORT static void Symmetry3DPlane(const double point[3], const double normalVector[3], int nbNodes, const double *coordsIn, double *coordsOut); + MEDCOUPLING_EXPORT static void GiveBaseForPlane(const double normalVector[3], double baseOfPlane[9]); public: MEDCOUPLING_EXPORT void getTinySerializationIntInformation(std::vector& tinyInfo) const; MEDCOUPLING_EXPORT void getTinySerializationStrInformation(std::vector& tinyInfo) const; diff --git a/src/MEDCoupling/MEDCouplingPointSet.cxx b/src/MEDCoupling/MEDCouplingPointSet.cxx index 299d21948..b3761c2d9 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.cxx +++ b/src/MEDCoupling/MEDCouplingPointSet.cxx @@ -994,8 +994,8 @@ bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* */ bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bb1, const double* bb2, int dim, double eps) { - double* bbtemp = new double[2*dim]; - double deltamax=0.0; + double* bbtemp(new double[2*dim]); + double deltamax(0.0); for (int i=0; i< dim; i++) { @@ -1011,7 +1011,7 @@ bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBou bbtemp[i*2+1]=bb2[i*2+1]+deltamax*eps; } - bool intersects = !bb1.isDisjointWith( bbtemp ); + bool intersects(!bb1.isDisjointWith(bbtemp)); delete [] bbtemp; return intersects; } @@ -1021,49 +1021,9 @@ bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBou */ void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle) { - double *coords=_coords->getPointer(); - int nbNodes=getNumberOfNodes(); - Rotate3DAlg(center,vect,angle,nbNodes,coords); -} - -/*! - * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords' - * around an axe ('center','vect') and with angle 'angle'. - */ -void MEDCouplingPointSet::Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords) -{ - if(!center || !vect) - throw INTERP_KERNEL::Exception("MEDCouplingPointSet::Rotate3DAlg : null vector in input !"); - double sina=sin(angle); - double cosa=cos(angle); - double vectorNorm[3]; - double matrix[9]; - double matrixTmp[9]; - double norm=sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]); - if(norm::min()) - throw INTERP_KERNEL::Exception("MEDCouplingPointSet::Rotate3DAlg : magnitude of input vector is too close of 0. !"); - std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies(),1/norm)); - //rotation matrix computation - matrix[0]=cosa; matrix[1]=0.; matrix[2]=0.; matrix[3]=0.; matrix[4]=cosa; matrix[5]=0.; matrix[6]=0.; matrix[7]=0.; matrix[8]=cosa; - matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2]; - matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2]; - matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2]; - std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies(),1-cosa)); - std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus()); - matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1]; - matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0]; - matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.; - std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies(),sina)); - std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus()); - //rotation matrix computed. - double tmp[3]; - for(int i=0; i()); - coords[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0]; - coords[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1]; - coords[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2]; - } + double *coords(_coords->getPointer()); + int nbNodes(getNumberOfNodes()); + DataArrayDouble::Rotate3DAlg(center,vect,angle,nbNodes,coords,coords); } /*! @@ -1157,7 +1117,7 @@ MEDCouplingMesh *MEDCouplingPointSet::buildPartRange(int beginCellIds, int endCe */ MEDCouplingMesh *MEDCouplingPointSet::buildPartRangeAndReduceNodes(int beginCellIds, int endCellIds, int stepCellIds, int& beginOut, int& endOut, int& stepOut, DataArrayInt*& arr) const { - MCAuto ret=buildPartOfMySelfSlice(beginCellIds,endCellIds,stepCellIds,true); + MCAuto ret(buildPartOfMySelfSlice(beginCellIds,endCellIds,stepCellIds,true)); arr=ret->zipCoordsTraducer(); return ret.retn(); } @@ -1167,28 +1127,9 @@ MEDCouplingMesh *MEDCouplingPointSet::buildPartRangeAndReduceNodes(int beginCell */ void MEDCouplingPointSet::rotate2D(const double *center, double angle) { - double *coords=_coords->getPointer(); - int nbNodes=getNumberOfNodes(); - Rotate2DAlg(center,angle,nbNodes,coords); -} - -/*! - * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords' - * around the center point 'center' and with angle 'angle'. - */ -void MEDCouplingPointSet::Rotate2DAlg(const double *center, double angle, int nbNodes, double *coords) -{ - double cosa=cos(angle); - double sina=sin(angle); - double matrix[4]; - matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa; - double tmp[2]; - for(int i=0; i()); - coords[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0]; - coords[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1]; - } + double *coords(_coords->getPointer()); + int nbNodes(getNumberOfNodes()); + DataArrayDouble::Rotate2DAlg(center,angle,nbNodes,coords,coords); } /// @cond INTERNAL @@ -1212,8 +1153,8 @@ public: */ void MEDCouplingPointSet::project2DCellOnXY(const int *startConn, const int *endConn, std::vector& res) const { - const double *coords=_coords->getConstPointer(); - int spaceDim=getSpaceDimension(); + const double *coords(_coords->getConstPointer()); + int spaceDim(getSpaceDimension()); for(const int *it=startConn;it!=endConn;it++) res.insert(res.end(),coords+spaceDim*(*it),coords+spaceDim*(*it+1)); if(spaceDim==2) @@ -1239,7 +1180,7 @@ void MEDCouplingPointSet::project2DCellOnXY(const int *startConn, const int *end */ bool MEDCouplingPointSet::isButterfly2DCell(const std::vector& res, bool isQuad, double eps) { - std::size_t nbOfNodes=res.size()/2; + std::size_t nbOfNodes(res.size()/2); std::vector nodes(nbOfNodes); for(std::size_t i=0;i& res, bool pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes); else pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes); - bool ret=pol->isButterflyAbs(); + bool ret(pol->isButterflyAbs()); delete pol; return ret; } @@ -1269,7 +1210,7 @@ bool MEDCouplingPointSet::areCellsFrom2MeshEqual(const MEDCouplingPointSet *othe std::vector c1,c2; getNodeIdsOfCell(cellId,c1); other->getNodeIdsOfCell(cellId,c2); - std::size_t sz=c1.size(); + std::size_t sz(c1.size()); if(sz!=c2.size()) return false; for(std::size_t i=0;i& ms); MEDCOUPLING_EXPORT static MEDCouplingPointSet *BuildInstanceFromMeshType(MEDCouplingMeshType type); - MEDCOUPLING_EXPORT static void Rotate2DAlg(const double *center, double angle, int nbNodes, double *coords); - MEDCOUPLING_EXPORT static void Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords); MEDCOUPLING_EXPORT static DataArrayInt *ComputeNbOfInteractionsWithSrcCells(const MEDCouplingPointSet *srcMesh, const MEDCouplingPointSet *trgMesh, double eps); MEDCOUPLING_EXPORT MEDCouplingMesh *buildPart(const int *start, const int *end) const; MEDCOUPLING_EXPORT MEDCouplingMesh *buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const; diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index ea521f675..a3bbcc180 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -4009,7 +4009,7 @@ DataArrayInt *MEDCouplingUMesh::getCellIdsCrossingPlane(const double *origin, co MCAuto coo=_coords->deepCopy(); double normm2(sqrt(vec2[0]*vec2[0]+vec2[1]*vec2[1]+vec2[2]*vec2[2])); if(normm2/normm>1e-6) - MEDCouplingPointSet::Rotate3DAlg(origin,vec2,angle,coo->getNumberOfTuples(),coo->getPointer()); + DataArrayDouble::Rotate3DAlg(origin,vec2,angle,coo->getNumberOfTuples(),coo->getPointer(),coo->getPointer()); MCAuto mw=clone(false);//false -> shallow copy mw->setCoords(coo); mw->getBoundingBox(bbox); @@ -7963,7 +7963,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshes(const MEDCouplingUMesh *mesh1, * \throw If \a a[ *i* ]->getMeshDimension() < 0. * \throw If the meshes in \a a are of different dimension (getMeshDimension()). */ -MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshes(std::vector& a) +MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshes(const std::vector& a) { std::size_t sz=a.size(); if(sz==0) @@ -7996,7 +7996,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshes(std::vector& a) +MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshesLL(const std::vector& a) { if(a.empty()) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::MergeUMeshes : input array must be NON EMPTY !"); diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 68093553c..9a8fc33c2 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -232,7 +232,7 @@ namespace MEDCoupling MEDCOUPLING_EXPORT int split2DCells(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *midOpt=0, const DataArrayInt *midOptI=0); MEDCOUPLING_EXPORT static MEDCouplingUMesh *Build0DMeshFromCoords(DataArrayDouble *da); MEDCOUPLING_EXPORT static MEDCouplingUMesh *MergeUMeshes(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2); - MEDCOUPLING_EXPORT static MEDCouplingUMesh *MergeUMeshes(std::vector& a); + MEDCOUPLING_EXPORT static MEDCouplingUMesh *MergeUMeshes(const std::vector& a); MEDCOUPLING_EXPORT static MEDCouplingUMesh *MergeUMeshesOnSameCoords(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2); MEDCOUPLING_EXPORT static MEDCouplingUMesh *MergeUMeshesOnSameCoords(const std::vector& meshes); MEDCOUPLING_EXPORT static MEDCouplingUMesh *FuseUMeshesOnSameCoords(const std::vector& meshes, int compType, std::vector& corr); @@ -311,7 +311,7 @@ namespace MEDCoupling void getCellsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints, double eps, MCAuto& elts, MCAuto& eltsIndex) const; /// @cond INTERNAL - static MEDCouplingUMesh *MergeUMeshesLL(std::vector& a); + static MEDCouplingUMesh *MergeUMeshesLL(const std::vector& a); typedef int (*DimM1DescNbrer)(int id, unsigned nb, const INTERP_KERNEL::CellModel& cm, bool compute, const int *conn1, const int *conn2); template MEDCouplingUMesh *buildDescendingConnectivityGen(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx, DimM1DescNbrer nbrer) const; diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx index 1373d886e..0348dafb8 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx @@ -2175,7 +2175,7 @@ void MEDCouplingBasicsTest1::testGetCellsContainingPoint() CPPUNIT_ASSERT(std::equal(t2->begin(),t2->end(),expectedValues2)); //2D with no help of bounding box. double center[2]={0.2,0.2}; - MEDCouplingPointSet::Rotate2DAlg(center,0.78539816339744830962,6,pos); + DataArrayDouble::Rotate2DAlg(center,0.78539816339744830962,6,pos,pos); targetMesh->rotate(center,0,0.78539816339744830962); targetMesh->getCellsContainingPoints(pos,6,1e-12,t1,t2); CPPUNIT_ASSERT_EQUAL(6,(int)t1->getNbOfElems()); diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py index ad084d3d0..57951ce92 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py @@ -4375,6 +4375,41 @@ class MEDCouplingBasicsTest5(unittest.TestCase): self.assertEqual(f1.isEqual(f4,1e-12,0.),expected) pass pass + + def testDADSymmetry1(self): + arr=DataArrayDouble([2,3,4],1,3) + res=arr.symmetry3DPlane([0.,0.,0.],[0.,0.,2.]) + self.assertTrue(res.isEqual(DataArrayDouble([2,3,-4],1,3),1e-14)) + # + res=arr.symmetry3DPlane([-1000,100,-1],[0.,0.,2.]) + self.assertTrue(res.isEqual(DataArrayDouble([2,3,-6],1,3),1e-14)) + # + res=arr.symmetry3DPlane([0,0,0],[1.,0.,0.]) + self.assertTrue(res.isEqual(DataArrayDouble([-2,3,4],1,3),1e-14)) + # + res=arr.symmetry3DPlane([0,0,0],[0.,1.,0.]) + self.assertTrue(res.isEqual(DataArrayDouble([2,-3,4],1,3),1e-14)) + # + res=arr.symmetry3DPlane([0,0,0],[-1.,1.,0.]) + self.assertTrue(res.isEqual(DataArrayDouble([3,2,4],1,3),1e-14)) + # + plane=[5.,4.,-7.] + a=DataArrayDouble(DataArrayDouble.GiveBaseForPlane(plane)) + self.assertAlmostEqual(DataArrayDouble.Dot(a[0],a[1]).magnitude()[0],0.,13) + self.assertAlmostEqual(DataArrayDouble.Dot(a[0],a[2]).magnitude()[0],0.,13) + self.assertAlmostEqual(DataArrayDouble.Dot(a[1],a[2]).magnitude()[0],0.,13) + coo=DataArrayDouble.Aggregate([10*a[0]+10*a[1],-10*a[0]+10*a[1],-10*a[0]-10*a[1],10*a[0]-10*a[1]]) + m=MEDCouplingUMesh("",2) ; m.setCoords(coo) ; m.allocateCells() + m.insertNextCell(NORM_QUAD4,[0,1,2,3]) + d,_=m.distanceToPoint(arr) + res=arr.symmetry3DPlane([0.,0.,0.],plane) # + d2,_=m.distanceToPoint(res) + self.assertAlmostEqual(abs(d-d2),0.,12) + self.assertAlmostEqual(DataArrayDouble.Dot(res-arr,a[0])[0],0.,12) + self.assertAlmostEqual(DataArrayDouble.Dot(res-arr,a[1])[0],0.,12) + self.assertAlmostEqual((res-arr).magnitude()[0]-2*d,0.,12) + self.assertTrue(res.isEqual(DataArrayDouble([2.666666666666667,3.5333333333333333,3.0666666666666666],1,3),1e-12)) + pass pass diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index eebe898c7..cb145e880 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -1621,7 +1621,7 @@ namespace MEDCoupling int sz; INTERP_KERNEL::AutoCPtr c=convertPyToNewDblArr2(center,&sz); INTERP_KERNEL::AutoCPtr coo=convertPyToNewDblArr2(coords,&sz); - MEDCoupling::MEDCouplingPointSet::Rotate2DAlg(c,angle,nbNodes,coo); + MEDCoupling::DataArrayDouble::Rotate2DAlg(c,angle,nbNodes,coo,coo); for(int i=0;i(coo)); + MEDCoupling::DataArrayDouble::Rotate2DAlg(c,angle,nbNodes,coo,const_cast(coo)); } static void Rotate3DAlg(PyObject *center, PyObject *vect, double angle, int nbNodes, PyObject *coords) throw(INTERP_KERNEL::Exception) @@ -1646,7 +1646,7 @@ namespace MEDCoupling INTERP_KERNEL::AutoCPtr c=convertPyToNewDblArr2(center,&sz); INTERP_KERNEL::AutoCPtr coo=convertPyToNewDblArr2(coords,&sz); INTERP_KERNEL::AutoCPtr v=convertPyToNewDblArr2(vect,&sz2); - MEDCoupling::MEDCouplingPointSet::Rotate3DAlg(c,v,angle,nbNodes,coo); + MEDCoupling::DataArrayDouble::Rotate3DAlg(c,v,angle,nbNodes,coo,coo); for(int i=0;i v=convertPyToNewDblArr2(vect,&sz2); - MEDCoupling::MEDCouplingPointSet::Rotate3DAlg(c,v,angle,nbNodes,const_cast(coo)); + MEDCoupling::DataArrayDouble::Rotate3DAlg(c,v,angle,nbNodes,coo,const_cast(coo)); } } }; diff --git a/src/MEDCoupling_Swig/MEDCouplingMemArray.i b/src/MEDCoupling_Swig/MEDCouplingMemArray.i index 28e239007..bb63fbb5b 100644 --- a/src/MEDCoupling_Swig/MEDCouplingMemArray.i +++ b/src/MEDCoupling_Swig/MEDCouplingMemArray.i @@ -151,6 +151,7 @@ %newobject MEDCoupling::DataArrayDouble::Multiply; %newobject MEDCoupling::DataArrayDouble::Divide; %newobject MEDCoupling::DataArrayDouble::Pow; +%newobject MEDCoupling::DataArrayDouble::symmetry3DPlane; %newobject MEDCoupling::DataArrayDouble::subArray; %newobject MEDCoupling::DataArrayDouble::changeNbOfComponents; %newobject MEDCoupling::DataArrayDouble::accumulatePerChunck; @@ -870,6 +871,34 @@ namespace MEDCoupling return convertDblArrToPyListOfTuple(vals,nbOfComp,nbOfTuples); } + DataArrayDouble *symmetry3DPlane(PyObject *point, PyObject *normalVector) throw(INTERP_KERNEL::Exception) + { + const char msg[]="Python wrap of DataArrayDouble::symmetry3DPlane : "; + double val,val2; + DataArrayDouble *a,*a2; + DataArrayDoubleTuple *aa,*aa2; + std::vector bb,bb2; + int sw; + const double *centerPtr(convertObjToPossibleCpp5_Safe(point,sw,val,a,aa,bb,msg,1,3,true)); + const double *vectorPtr(convertObjToPossibleCpp5_Safe(normalVector,sw,val2,a2,aa2,bb2,msg,1,3,true)); + MCAuto ret(self->symmetry3DPlane(centerPtr,vectorPtr)); + return ret.retn(); + } + + static PyObject *GiveBaseForPlane(PyObject *normalVector) throw(INTERP_KERNEL::Exception) + { + const char msg[]="Python wrap of DataArrayDouble::GiveBaseForPlane : "; + double val,val2; + DataArrayDouble *a,*a2; + DataArrayDoubleTuple *aa,*aa2; + std::vector bb,bb2; + int sw; + const double *vectorPtr(convertObjToPossibleCpp5_Safe(normalVector,sw,val,a,aa,bb,msg,1,3,true)); + double res[9]; + DataArrayDouble::GiveBaseForPlane(vectorPtr,res); + return convertDblArrToPyListOfTuple(res,3,3); + } + DataArrayDouble *renumber(PyObject *li) throw(INTERP_KERNEL::Exception) { void *da=0; diff --git a/src/MEDLoader/MEDFileData.cxx b/src/MEDLoader/MEDFileData.cxx index 96e77a72c..72606a8d8 100644 --- a/src/MEDLoader/MEDFileData.cxx +++ b/src/MEDLoader/MEDFileData.cxx @@ -219,6 +219,74 @@ bool MEDFileData::unPolyzeMeshes() return !meshesImpacted.empty(); } +/*! + * Precondition : all instances in \a mfds should have a single mesh with fields on it. If there is an instance with not exactly one mesh an exception will be thrown. + * You can invoke MEDFileFields::partOfThisLyingOnSpecifiedMeshName method to make it work. + */ +MCAuto MEDFileData::Aggregate(const std::vector& mfds) +{ + if(mfds.empty()) + throw INTERP_KERNEL::Exception("MEDFileData::Aggregate : empty vector !"); + std::size_t sz(mfds.size()),i(0); + MCAuto ret(MEDFileData::New()); + std::vector ms(sz); + std::vector< std::vector< std::pair > > dts(sz); + for(std::vector::const_iterator it=mfds.begin();it!=mfds.end();it++,i++) + { + const MEDFileData *elt(*it); + if(!elt) + throw INTERP_KERNEL::Exception("MEDFileData::Aggregate : presence of NULL pointer !"); + const MEDFileMeshes *meshes(elt->getMeshes()); + if(!meshes) + throw INTERP_KERNEL::Exception("MEDFileData::Aggregate : presence of an instance with no meshes attached on it !"); + if(meshes->getNumberOfMeshes()!=1) + throw INTERP_KERNEL::Exception("MEDFileData::Aggregate : all instances in input vector must lie on exactly one mesh ! To have it you can invoke partOfThisLyingOnSpecifiedMeshName method."); + const MEDFileMesh *mesh(meshes->getMeshAtPos(0)); + if(!mesh) + throw INTERP_KERNEL::Exception("MEDFileData::Aggregate : presence of null mesh in a MEDFileData instance among input vector !"); + const MEDFileUMesh *umesh(dynamic_cast(mesh)); + if(!umesh) + throw INTERP_KERNEL::Exception("MEDFileData::Aggregate : works only for unstructured meshes !"); + ms[i]=umesh; + dts[i]=umesh->getAllDistributionOfTypes(); + } + MCAuto agg_m(MEDFileUMesh::Aggregate(ms)); + MCAuto mss(MEDFileMeshes::New()); mss->pushMesh(agg_m); + ret->setMeshes(mss); + // fields + std::vector fieldNames(mfds[0]->getFields()->getFieldsNames()); + std::set fieldNamess(fieldNames.begin(),fieldNames.end()); + if(fieldNames.size()!=fieldNamess.size()) + throw INTERP_KERNEL::Exception("MEDFileData::Aggregate : field names must be different each other !"); + std::vector< std::vector > vectOfFields(fieldNames.size()); + std::vector< std::vector< MCAuto< MEDFileAnyTypeFieldMultiTS > > > vectOfFields2(fieldNames.size()); + MCAuto fss(MEDFileFields::New()); + for(std::vector::const_iterator it=mfds.begin();it!=mfds.end();it++) + { + std::vector fieldNames0((*it)->getFields()->getFieldsNames()); + std::set fieldNamess0(fieldNames0.begin(),fieldNames0.end()); + if(fieldNamess!=fieldNamess0) + throw INTERP_KERNEL::Exception("MEDFileData::Aggregate : field names must be the same for all input instances !"); + i=0; + for(std::vector::const_iterator it1=fieldNames.begin();it1!=fieldNames.end();it1++,i++) + { + MCAuto fmts((*it)->getFields()->getFieldWithName(*it1)); + if(fmts.isNull()) + throw INTERP_KERNEL::Exception("MEDFileData::Aggregate : internal error 1 !"); + vectOfFields2[i].push_back(fmts); vectOfFields[i].push_back(fmts); + } + } + i=0; + for(std::vector::const_iterator it1=fieldNames.begin();it1!=fieldNames.end();it1++,i++) + { + MCAuto fmts(MEDFileAnyTypeFieldMultiTS::Aggregate(vectOfFields[i],dts)); + fmts->setMeshName(agg_m->getName()); + fss->pushField(fmts); + } + ret->setFields(fss); + return ret; +} + MEDFileData::MEDFileData() { } diff --git a/src/MEDLoader/MEDFileData.hxx b/src/MEDLoader/MEDFileData.hxx index d09cfd33e..bdb21eb97 100644 --- a/src/MEDLoader/MEDFileData.hxx +++ b/src/MEDLoader/MEDFileData.hxx @@ -53,6 +53,7 @@ namespace MEDCoupling MEDLOADER_EXPORT bool changeMeshNames(const std::vector< std::pair >& modifTab); MEDLOADER_EXPORT bool changeMeshName(const std::string& oldMeshName, const std::string& newMeshName); MEDLOADER_EXPORT bool unPolyzeMeshes(); + MEDLOADER_EXPORT static MCAuto Aggregate(const std::vector& mfds); // MEDLOADER_EXPORT void write(const std::string& fileName, int mode) const; private: diff --git a/src/MEDLoader/MEDFileField.cxx b/src/MEDLoader/MEDFileField.cxx index 612cb82d9..81292d840 100644 --- a/src/MEDLoader/MEDFileField.cxx +++ b/src/MEDLoader/MEDFileField.cxx @@ -18,9 +18,10 @@ // // Author : Anthony Geay (CEA/DEN) -#include "MEDFileField.hxx" +#include "MEDFileField.txx" #include "MEDFileMesh.hxx" #include "MEDLoaderBase.hxx" +#include "MEDLoaderTraits.hxx" #include "MEDFileUtilities.hxx" #include "MEDFileSafeCaller.txx" #include "MEDFileFieldOverView.hxx" @@ -42,6 +43,9 @@ extern med_geometry_type typmai3[34]; using namespace MEDCoupling; +template class MEDFileField1TSTemplateWithoutSDA; +template class MEDFileField1TSTemplateWithoutSDA; + const char MEDFileField1TSWithoutSDA::TYPE_STR[]="FLOAT64"; const char MEDFileIntField1TSWithoutSDA::TYPE_STR[]="INT32"; @@ -1579,10 +1583,6 @@ bool MEDFileFieldPerMeshPerType::keepOnlyGaussDiscretization(std::size_t idOfDis return true; } -MEDFileFieldPerMeshPerType::MEDFileFieldPerMeshPerType(MEDFileFieldPerMesh *fath, INTERP_KERNEL::NormalizedCellType geoType):_father(fath),_geo_type(geoType) -{ -} - MEDFileFieldPerMeshPerType::MEDFileFieldPerMeshPerType(med_idt fid, MEDFileFieldPerMesh *fath, TypeOfField type, INTERP_KERNEL::NormalizedCellType geoType, const MEDFileFieldNameScope& nasc, const PartDefinition *pd):_father(fath),_geo_type(geoType) { INTERP_KERNEL::AutoPtr pflName=MEDLoaderBase::buildEmptyString(MED_NAME_SIZE); @@ -2341,6 +2341,90 @@ const MEDFileFieldPerMeshPerTypePerDisc *MEDFileFieldPerMesh::getLeafGivenTypeAn throw INTERP_KERNEL::Exception(oss.str()); } +/*! + * \param [in,out] start - Integer that gives the current position in the final aggregated array + * \param [in] pms - list of elements to aggregate. integer gives the mesh id + * \param [in] dts - (Distribution of types) = level 1 : meshes to aggregate. Level 2 : all geo type. Level 3 pair specifying geo type and number of elem in geotype. + * \param [out] extractInfo - Gives information about the where the data comes from. It is a vector of triplet. First element in the triplet the mesh pos. The 2nd one the start pos. The 3rd the end pos. + */ +MCAuto MEDFileFieldPerMeshPerTypePerDisc::Aggregate(int &start, const std::vector< std::pair >& pms, const std::vector< std::vector< std::pair > >& dts, TypeOfField tof, MEDFileFieldPerMeshPerType *father, std::vector > >& extractInfo) +{ + MCAuto ret(new MEDFileFieldPerMeshPerTypePerDisc(father,tof)); + if(pms.empty()) + throw INTERP_KERNEL::Exception("MEDFileFieldPerMeshPerTypePerDisc::Aggregate : empty input vector !"); + for(std::vector >::const_iterator it=pms.begin();it!=pms.end();it++) + { + if(!(*it).second) + throw INTERP_KERNEL::Exception("MEDFileFieldPerMeshPerTypePerDisc::Aggregate : presence of null pointer !"); + if(!(*it).second->getProfile().empty()) + throw INTERP_KERNEL::Exception("MEDFileFieldPerMeshPerTypePerDisc::Aggregate : not implemented yet for profiles !"); + if(!(*it).second->getLocalization().empty()) + throw INTERP_KERNEL::Exception("MEDFileFieldPerMeshPerTypePerDisc::Aggregate : not implemented yet for gauss pts !"); + } + INTERP_KERNEL::NormalizedCellType gt(pms[0].second->getGeoType()); + std::size_t i(0); + std::vector< std::pair > filteredDTS; + for(std::vector< std::vector< std::pair > >::const_iterator it=dts.begin();it!=dts.end();it++,i++) + for(std::vector< std::pair >::const_iterator it2=(*it).begin();it2!=(*it).end();it2++) + if((*it2).first==gt) + filteredDTS.push_back(std::pair(i,(*it2).second)); + if(pms.size()!=filteredDTS.size()) + throw INTERP_KERNEL::Exception("MEDFileFieldPerMeshPerTypePerDisc::Aggregate : not implemented yet for generated profiles !"); + std::vector >::const_iterator it1(pms.begin()); + std::vector< std::pair >::const_iterator it2(filteredDTS.begin()); + int zeStart(start),nval(0); + for(;it1!=pms.end();it1++,it2++) + { + if((*it1).first!=(*it2).first) + throw INTERP_KERNEL::Exception("MEDFileFieldPerMeshPerTypePerDisc::Aggregate : not implemented yet for generated profiles 2 !"); + int s1((*it1).second->getStart()),e1((*it1).second->getEnd()); + extractInfo.push_back(std::pair >((*it1).first,std::pair(s1,e1))); + start+=e1-s1; + nval+=((*it1).second)->getNumberOfVals(); + } + ret->_start=zeStart; ret->_end=start; ret->_nval=nval; + return ret; +} + +MCAuto MEDFileFieldPerMeshPerType::Aggregate(int &start, const std::vector >& pms, const std::vector< std::vector< std::pair > >& dts, INTERP_KERNEL::NormalizedCellType gt, MEDFileFieldPerMesh *father, std::vector > >& extractInfo) +{ + MCAuto ret(new MEDFileFieldPerMeshPerType(father,gt)); + std::map > > m; + for(std::vector >::const_iterator it=pms.begin();it!=pms.end();it++) + { + for(std::vector< MCAuto >::const_iterator it2=(*it).second->_field_pm_pt_pd.begin();it2!=(*it).second->_field_pm_pt_pd.end();it2++) + m[(*it2)->getType()].push_back(std::pair((*it).first,*it2)); + } + for(std::map > >::const_iterator it=m.begin();it!=m.end();it++) + { + MCAuto agg(MEDFileFieldPerMeshPerTypePerDisc::Aggregate(start,(*it).second,dts,(*it).first,ret,extractInfo)); + ret->_field_pm_pt_pd.push_back(agg); + } + return ret; +} + +MCAuto MEDFileFieldPerMesh::Aggregate(int &start, const std::vector& pms, const std::vector< std::vector< std::pair > >& dts, MEDFileAnyTypeField1TSWithoutSDA *father, std::vector > >& extractInfo) +{ + MCAuto ret(new MEDFileFieldPerMesh(father,pms[0]->getMeshName(),pms[0]->getMeshIteration(),pms[0]->getMeshOrder())); + std::map > > m; + std::size_t i(0); + for(std::vector::const_iterator it=pms.begin();it!=pms.end();it++,i++) + { + const std::vector< MCAuto< MEDFileFieldPerMeshPerType > >& v((*it)->_field_pm_pt); + for(std::vector< MCAuto< MEDFileFieldPerMeshPerType > >::const_iterator it2=v.begin();it2!=v.end();it2++) + { + INTERP_KERNEL::NormalizedCellType gt((*it2)->getGeoType()); + m[gt].push_back(std::pair(i,*it2)); + } + } + for(std::map > >::const_iterator it=m.begin();it!=m.end();it++) + { + MCAuto agg(MEDFileFieldPerMeshPerType::Aggregate(start,(*it).second,dts,(*it).first,ret,extractInfo)); + ret->_field_pm_pt.push_back(agg); + } + return ret; +} + int MEDFileFieldPerMesh::addNewEntryIfNecessary(INTERP_KERNEL::NormalizedCellType type) { int i=0; @@ -5114,22 +5198,6 @@ std::vector< std::vector > MEDFileField1TSWithoutSDA::getFiel return ret; } -/*! - * Returns a pointer to the underground DataArrayDouble instance. So the - * caller should not decrRef() it. This method allows for a direct access to the field - * values. This method is quite unusable if there is more than a nodal field or a cell - * field on single geometric cell type. - * \return DataArrayDouble * - the pointer to the field values array. - */ -DataArrayDouble *MEDFileField1TSWithoutSDA::getUndergroundDataArrayDouble() const -{ - const DataArrayDouble *ret=_arr; - if(ret) - return const_cast(ret); - else - return 0; -} - const char *MEDFileField1TSWithoutSDA::getTypeStr() const { return TYPE_STR; @@ -5149,18 +5217,6 @@ MEDFileIntField1TSWithoutSDA *MEDFileField1TSWithoutSDA::convertToInt() const return ret.retn(); } -/*! - * Returns a pointer to the underground DataArrayDouble instance. So the - * caller should not decrRef() it. This method allows for a direct access to the field - * values. This method is quite unusable if there is more than a nodal field or a cell - * field on single geometric cell type. - * \return DataArrayDouble * - the pointer to the field values array. - */ -DataArray *MEDFileField1TSWithoutSDA::getUndergroundDataArray() const -{ - return getUndergroundDataArrayDouble(); -} - /*! * Returns a pointer to the underground DataArrayDouble instance and a * sequence describing parameters of a support of each part of \a this field. The @@ -5184,7 +5240,7 @@ DataArrayDouble *MEDFileField1TSWithoutSDA::getUndergroundDataArrayDoubleExt(std if(_field_per_mesh[0]==0) throw INTERP_KERNEL::Exception("MEDFileField1TSWithoutSDA::getUndergroundDataArrayExt : no field specified !"); _field_per_mesh[0]->getUndergroundDataArrayExt(entries); - return getUndergroundDataArrayDouble(); + return getUndergroundDataArrayTemplate(); } /*! @@ -5208,13 +5264,13 @@ DataArray *MEDFileField1TSWithoutSDA::getUndergroundDataArrayExt(std::vector< st return getUndergroundDataArrayDoubleExt(entries); } -MEDFileField1TSWithoutSDA::MEDFileField1TSWithoutSDA(const std::string& fieldName, int csit, int iteration, int order, const std::vector& infos):MEDFileAnyTypeField1TSWithoutSDA(fieldName,csit,iteration,order) +MEDFileField1TSWithoutSDA::MEDFileField1TSWithoutSDA(const std::string& fieldName, int csit, int iteration, int order, const std::vector& infos):MEDFileField1TSTemplateWithoutSDA(fieldName,csit,iteration,order) { - DataArrayDouble *arr(getOrCreateAndGetArrayDouble()); + DataArrayDouble *arr(getOrCreateAndGetArrayTemplate()); arr->setInfoAndChangeNbOfCompo(infos); } -MEDFileField1TSWithoutSDA::MEDFileField1TSWithoutSDA():MEDFileAnyTypeField1TSWithoutSDA() +MEDFileField1TSWithoutSDA::MEDFileField1TSWithoutSDA():MEDFileField1TSTemplateWithoutSDA() { } @@ -5233,57 +5289,6 @@ MEDFileAnyTypeField1TSWithoutSDA *MEDFileField1TSWithoutSDA::deepCopy() const return ret.retn(); } -void MEDFileField1TSWithoutSDA::setArray(DataArray *arr) -{ - if(!arr) - { - _nb_of_tuples_to_be_allocated=-1; - _arr=0; - return ; - } - DataArrayDouble *arrC=dynamic_cast(arr); - if(!arrC) - throw INTERP_KERNEL::Exception("MEDFileField1TSWithoutSDA::setArray : the input not null array is not of type DataArrayDouble !"); - else - _nb_of_tuples_to_be_allocated=-3; - arrC->incrRef(); - _arr=arrC; -} - -DataArray *MEDFileField1TSWithoutSDA::createNewEmptyDataArrayInstance() const -{ - return DataArrayDouble::New(); -} - -DataArrayDouble *MEDFileField1TSWithoutSDA::getOrCreateAndGetArrayDouble() -{ - DataArrayDouble *ret=_arr; - if(ret) - return ret; - _arr=DataArrayDouble::New(); - return _arr; -} - -DataArray *MEDFileField1TSWithoutSDA::getOrCreateAndGetArray() -{ - return getOrCreateAndGetArrayDouble(); -} - -const DataArrayDouble *MEDFileField1TSWithoutSDA::getOrCreateAndGetArrayDouble() const -{ - const DataArrayDouble *ret=_arr; - if(ret) - return ret; - DataArrayDouble *ret2=DataArrayDouble::New(); - const_cast(this)->_arr=DataArrayDouble::New(); - return ret2; -} - -const DataArray *MEDFileField1TSWithoutSDA::getOrCreateAndGetArray() const -{ - return getOrCreateAndGetArrayDouble(); -} - //= MEDFileIntField1TSWithoutSDA MEDFileIntField1TSWithoutSDA *MEDFileIntField1TSWithoutSDA::New(const std::string& fieldName, int csit, int iteration, int order, const std::vector& infos) @@ -5291,14 +5296,14 @@ MEDFileIntField1TSWithoutSDA *MEDFileIntField1TSWithoutSDA::New(const std::strin return new MEDFileIntField1TSWithoutSDA(fieldName,csit,iteration,order,infos); } -MEDFileIntField1TSWithoutSDA::MEDFileIntField1TSWithoutSDA():MEDFileAnyTypeField1TSWithoutSDA() +MEDFileIntField1TSWithoutSDA::MEDFileIntField1TSWithoutSDA():MEDFileField1TSTemplateWithoutSDA() { } MEDFileIntField1TSWithoutSDA::MEDFileIntField1TSWithoutSDA(const std::string& fieldName, int csit, int iteration, int order, - const std::vector& infos):MEDFileAnyTypeField1TSWithoutSDA(fieldName,csit,iteration,order) + const std::vector& infos):MEDFileField1TSTemplateWithoutSDA(fieldName,csit,iteration,order) { - DataArrayInt *arr(getOrCreateAndGetArrayInt()); + DataArrayInt *arr(getOrCreateAndGetArrayTemplate()); arr->setInfoAndChangeNbOfCompo(infos); } @@ -5321,34 +5326,6 @@ MEDFileField1TSWithoutSDA *MEDFileIntField1TSWithoutSDA::convertToDouble() const return ret.retn(); } -/*! - * Returns a pointer to the underground DataArrayInt instance. So the - * caller should not decrRef() it. This method allows for a direct access to the field - * values. This method is quite unusable if there is more than a nodal field or a cell - * field on single geometric cell type. - * \return DataArrayInt * - the pointer to the field values array. - */ -DataArray *MEDFileIntField1TSWithoutSDA::getUndergroundDataArray() const -{ - return getUndergroundDataArrayInt(); -} - -/*! - * Returns a pointer to the underground DataArrayInt instance. So the - * caller should not decrRef() it. This method allows for a direct access to the field - * values. This method is quite unusable if there is more than a nodal field or a cell - * field on single geometric cell type. - * \return DataArrayInt * - the pointer to the field values array. - */ -DataArrayInt *MEDFileIntField1TSWithoutSDA::getUndergroundDataArrayInt() const -{ - const DataArrayInt *ret=_arr; - if(ret) - return const_cast(ret); - else - return 0; -} - /*! * Returns a pointer to the underground DataArrayInt instance and a * sequence describing parameters of a support of each part of \a this field. The @@ -5393,7 +5370,7 @@ DataArrayInt *MEDFileIntField1TSWithoutSDA::getUndergroundDataArrayIntExt(std::v if(_field_per_mesh[0]==0) throw INTERP_KERNEL::Exception("MEDFileField1TSWithoutSDA::getUndergroundDataArrayExt : no field specified !"); _field_per_mesh[0]->getUndergroundDataArrayExt(entries); - return getUndergroundDataArrayInt(); + return getUndergroundDataArrayTemplate(); } MEDFileAnyTypeField1TSWithoutSDA *MEDFileIntField1TSWithoutSDA::shallowCpy() const @@ -5411,57 +5388,6 @@ MEDFileAnyTypeField1TSWithoutSDA *MEDFileIntField1TSWithoutSDA::deepCopy() const return ret.retn(); } -void MEDFileIntField1TSWithoutSDA::setArray(DataArray *arr) -{ - if(!arr) - { - _nb_of_tuples_to_be_allocated=-1; - _arr=0; - return ; - } - DataArrayInt *arrC=dynamic_cast(arr); - if(!arrC) - throw INTERP_KERNEL::Exception("MEDFileIntField1TSWithoutSDA::setArray : the input not null array is not of type DataArrayInt !"); - else - _nb_of_tuples_to_be_allocated=-3; - arrC->incrRef(); - _arr=arrC; -} - -DataArray *MEDFileIntField1TSWithoutSDA::createNewEmptyDataArrayInstance() const -{ - return DataArrayInt::New(); -} - -DataArrayInt *MEDFileIntField1TSWithoutSDA::getOrCreateAndGetArrayInt() -{ - DataArrayInt *ret=_arr; - if(ret) - return ret; - _arr=DataArrayInt::New(); - return _arr; -} - -DataArray *MEDFileIntField1TSWithoutSDA::getOrCreateAndGetArray() -{ - return getOrCreateAndGetArrayInt(); -} - -const DataArrayInt *MEDFileIntField1TSWithoutSDA::getOrCreateAndGetArrayInt() const -{ - const DataArrayInt *ret=_arr; - if(ret) - return ret; - DataArrayInt *ret2=DataArrayInt::New(); - const_cast(this)->_arr=DataArrayInt::New(); - return ret2; -} - -const DataArray *MEDFileIntField1TSWithoutSDA::getOrCreateAndGetArray() const -{ - return getOrCreateAndGetArrayInt(); -} - MEDFileAnyTypeField1TS::MEDFileAnyTypeField1TS() { } @@ -6328,7 +6254,7 @@ MEDFileField1TS *MEDFileField1TS::New(const std::string& fileName, const std::st */ MEDFileField1TS *MEDFileField1TS::New(const MEDFileField1TSWithoutSDA& other, bool shallowCopyOfContent) { - MCAuto ret=new MEDFileField1TS(other,shallowCopyOfContent); + MCAuto ret(new MEDFileField1TS(other,shallowCopyOfContent)); ret->contentNotNull(); return ret.retn(); } @@ -6340,7 +6266,7 @@ MEDFileField1TS *MEDFileField1TS::New(const MEDFileField1TSWithoutSDA& other, bo */ MEDFileField1TS *MEDFileField1TS::New() { - MCAuto ret=new MEDFileField1TS; + MCAuto ret(new MEDFileField1TS); ret->contentNotNull(); return ret.retn(); } @@ -6761,7 +6687,7 @@ MEDFileAnyTypeField1TS *MEDFileField1TS::shallowCpy() const DataArrayDouble *MEDFileField1TS::getUndergroundDataArray() const { - return contentNotNull()->getUndergroundDataArrayDouble(); + return contentNotNull()->getUndergroundDataArrayTemplate(); } DataArrayDouble *MEDFileField1TS::getUndergroundDataArrayExt(std::vector< std::pair,std::pair > >& entries) const @@ -7168,7 +7094,7 @@ MEDFileIntField1TSWithoutSDA *MEDFileIntField1TS::contentNotNull() DataArrayInt *MEDFileIntField1TS::getUndergroundDataArray() const { - return contentNotNull()->getUndergroundDataArrayInt(); + return contentNotNull()->getUndergroundDataArrayTemplate(); } //= MEDFileAnyTypeFieldMultiTSWithoutSDA @@ -7489,6 +7415,8 @@ void MEDFileAnyTypeFieldMultiTSWithoutSDA::pushBackTimeStep(MCAutogetInfo()); } checkThatComponentsMatch(tse2->getInfo()); + if(getDtUnit().empty() && !tse->getDtUnit().empty()) + setDtUnit(tse->getDtUnit()); _time_steps.push_back(tse); } @@ -7540,6 +7468,7 @@ void MEDFileAnyTypeFieldMultiTSWithoutSDA::loadStructureOrStructureAndBigArraysR _time_steps[i]->loadStructureAndBigArraysRecursively(fid,*this,ms,entities); else _time_steps[i]->loadOnlyStructureOfDataRecursively(fid,*this,ms,entities); + synchronizeNameScope(); } } @@ -9094,6 +9023,97 @@ MEDFileAnyTypeFieldMultiTS *MEDFileAnyTypeFieldMultiTS::extractPart(const std::m return fmtsOut.retn(); } +template +MCAuto AggregateHelperF1TS(const std::vector< typename MLFieldTraits::F1TSType const *>& f1tss, const std::vector< std::vector< std::pair > >& dts) +{ + MCAuto< typename MLFieldTraits::F1TSType > ret(MLFieldTraits::F1TSType::New()); + if(f1tss.empty()) + throw INTERP_KERNEL::Exception("AggregateHelperF1TS : empty vector !"); + std::size_t sz(f1tss.size()),i(0); + std::vector< typename MLFieldTraits::F1TSWSDAType const *> f1tsw(sz); + for(typename std::vector< typename MLFieldTraits::F1TSType const *>::const_iterator it=f1tss.begin();it!=f1tss.end();it++,i++) + { + typename MLFieldTraits::F1TSType const *elt(*it); + if(!elt) + throw INTERP_KERNEL::Exception("AggregateHelperF1TS : presence of a null pointer !"); + f1tsw[i]=dynamic_cast::F1TSWSDAType const *>(elt->contentNotNullBase()); + } + typename MLFieldTraits::F1TSWSDAType *retc(dynamic_cast::F1TSWSDAType *>(ret->contentNotNullBase())); + if(!retc) + throw INTERP_KERNEL::Exception("AggregateHelperF1TS : internal error 1 !"); + retc->aggregate(f1tsw,dts); + ret->setDtUnit(f1tss[0]->getDtUnit()); + return DynamicCast::F1TSType , MEDFileAnyTypeField1TS>(ret); +} + +template +MCAuto< MEDFileAnyTypeFieldMultiTS > AggregateHelperFMTS(const std::vector< typename MLFieldTraits::FMTSType const *>& fmtss, const std::vector< std::vector< std::pair > >& dts) +{ + MCAuto< typename MLFieldTraits::FMTSType > ret(MLFieldTraits::FMTSType::New()); + if(fmtss.empty()) + throw INTERP_KERNEL::Exception("AggregateHelperFMTS : empty vector !"); + std::size_t sz(fmtss.size()); + for(typename std::vector< typename MLFieldTraits::FMTSType const *>::const_iterator it=fmtss.begin();it!=fmtss.end();it++) + { + typename MLFieldTraits::FMTSType const *elt(*it); + if(!elt) + throw INTERP_KERNEL::Exception("AggregateHelperFMTS : presence of null pointer !"); + } + int nbTS(fmtss[0]->getNumberOfTS()); + for(typename std::vector< typename MLFieldTraits::FMTSType const *>::const_iterator it=fmtss.begin();it!=fmtss.end();it++) + if((*it)->getNumberOfTS()!=nbTS) + throw INTERP_KERNEL::Exception("AggregateHelperFMTS : all fields must have the same number of TS !"); + for(int iterTS=0;iterTS::F1TSType const *> f1tss(sz); + std::vector< MCAuto::F1TSType> > f1tss2(sz); + for(typename std::vector< typename MLFieldTraits::FMTSType const *>::const_iterator it=fmtss.begin();it!=fmtss.end();it++,i++) + { f1tss2[i]=(*it)->getTimeStepAtPos(iterTS); f1tss[i]=f1tss2[i]; } + MCAuto f1ts(AggregateHelperF1TS(f1tss,dts)); + ret->pushBackTimeStep(f1ts); + ret->setDtUnit(f1ts->getDtUnit()); + } + return DynamicCast::FMTSType , MEDFileAnyTypeFieldMultiTS>(ret); +} + +/*! + * \a dts and \a ftmss are expected to have same size. + */ +MCAuto MEDFileAnyTypeFieldMultiTS::Aggregate(const std::vector& fmtss, const std::vector< std::vector< std::pair > >& dts) +{ + if(fmtss.empty()) + throw INTERP_KERNEL::Exception("MEDFileAnyTypeFieldMultiTS::Aggregate : input vector is empty !"); + std::size_t sz(fmtss.size()); + std::vector fmtss1; + std::vector fmtss2; + for(std::vector::const_iterator it=fmtss.begin();it!=fmtss.end();it++) + { + if(!(*it)) + throw INTERP_KERNEL::Exception("MEDFileAnyTypeFieldMultiTS::Aggregate : presence of null instance in input vector !"); + const MEDFileFieldMultiTS *elt1(dynamic_cast(*it)); + if(elt1) + { + fmtss1.push_back(elt1); + continue; + } + const MEDFileIntFieldMultiTS *elt2(dynamic_cast(*it)); + if(elt2) + { + fmtss2.push_back(elt2); + continue; + } + throw INTERP_KERNEL::Exception("MEDFileAnyTypeFieldMultiTS::Aggregate : not recognized type !"); + } + if(fmtss1.size()!=sz && fmtss2.size()!=sz) + throw INTERP_KERNEL::Exception("MEDFileAnyTypeFieldMultiTS::Aggregate : type of data is not homogeneous !"); + if(fmtss1.size()==sz) + return AggregateHelperFMTS(fmtss1,dts); + if(fmtss2.size()!=sz) + return AggregateHelperFMTS(fmtss2,dts); + throw INTERP_KERNEL::Exception("MEDFileAnyTypeFieldMultiTS::Aggregate : not implemented yet !"); +} + MEDFileAnyTypeFieldMultiTSIterator *MEDFileAnyTypeFieldMultiTS::iterator() { return new MEDFileAnyTypeFieldMultiTSIterator(this); @@ -9899,7 +9919,7 @@ DataArrayInt *MEDFileIntFieldMultiTS::getFieldWithProfile(TypeOfField type, int * delete this field using decrRef() as it is no more needed. * \throw If \a pos is not a valid time step id. */ -MEDFileAnyTypeField1TS *MEDFileIntFieldMultiTS::getTimeStepAtPos(int pos) const +MEDFileIntField1TS *MEDFileIntFieldMultiTS::getTimeStepAtPos(int pos) const { const MEDFileAnyTypeField1TSWithoutSDA *item=contentNotNullBase()->getTimeStepAtPos2(pos); if(!item) diff --git a/src/MEDLoader/MEDFileField.hxx b/src/MEDLoader/MEDFileField.hxx index a748d0b09..d26ae7a69 100644 --- a/src/MEDLoader/MEDFileField.hxx +++ b/src/MEDLoader/MEDFileField.hxx @@ -27,6 +27,8 @@ #include "MEDFileUtilities.hxx" #include "MCAuto.hxx" +#include "MEDLoaderTraits.hxx" +#include "MEDCouplingTraits.hxx" #include "MEDCouplingRefCountObject.hxx" #include "MEDCouplingFieldInt.hxx" #include "MEDCouplingMemArray.hxx" @@ -127,6 +129,7 @@ namespace MEDCoupling int getNumberOfTuples() const; int getStart() const { return _start; } int getEnd() const { return _end; } + int getNumberOfVals() const { return _nval; } DataArray *getOrCreateAndGetArray(); const DataArray *getOrCreateAndGetArray() const; const std::vector& getInfo() const; @@ -152,6 +155,8 @@ namespace MEDCoupling static MEDFileFieldPerMeshPerTypePerDisc *NewObjectOnSameDiscThanPool(TypeOfField typeF, INTERP_KERNEL::NormalizedCellType geoType, DataArrayInt *idsOfMeshElt, bool isPfl, int nbi, int offset, std::list< const MEDFileFieldPerMeshPerTypePerDisc *>& entriesOnSameDisc, MEDFileFieldGlobsReal& glob, bool ¬InExisting); + static MCAuto Aggregate(int &start, const std::vector >& pms, const std::vector< std::vector< std::pair > >& dts, TypeOfField tof, MEDFileFieldPerMeshPerType *father, std::vector > >& extractInfo); + MEDFileFieldPerMeshPerTypePerDisc(MEDFileFieldPerMeshPerType *fath, TypeOfField type):_type(type),_father(fath),_start(-1),_end(-1),_nval(-1),_loc_id(-5),_profile_it(-1) { } private: MEDFileFieldPerMeshPerTypePerDisc(MEDFileFieldPerMeshPerType *fath, TypeOfField type, int profileIt, const PartDefinition *pd); MEDFileFieldPerMeshPerTypePerDisc(MEDFileFieldPerMeshPerType *fath, TypeOfField type, int profileIt, const std::string& dummy); @@ -216,19 +221,21 @@ namespace MEDCoupling void changeLocsRefsNamesGen(const std::vector< std::pair, std::string > >& mapOfModif); MEDFileFieldPerMeshPerTypePerDisc *getLeafGivenLocId(int locId); const MEDFileFieldPerMeshPerTypePerDisc *getLeafGivenLocId(int locId) const; + int getNumberOfLoc() const { return _field_pm_pt_pd.size(); } void getFieldAtLevel(int meshDim, TypeOfField type, const MEDFileFieldGlobsReal *glob, std::vector< std::pair >& dads, std::vector& pfls, std::vector& locs, std::vector& geoTypes) const; void fillValues(int& startEntryId, std::vector< std::pair,std::pair > >& entries) const; void setLeaves(const std::vector< MCAuto< MEDFileFieldPerMeshPerTypePerDisc > >& leaves); bool keepOnlySpatialDiscretization(TypeOfField tof, int &globalNum, std::vector< std::pair >& its); bool keepOnlyGaussDiscretization(std::size_t idOfDisc, int &globalNum, std::vector< std::pair >& its); static med_entity_type ConvertIntoMEDFileType(TypeOfField ikType, INTERP_KERNEL::NormalizedCellType ikGeoType, med_geometry_type& medfGeoType); + static MCAuto Aggregate(int &start, const std::vector< std::pair >& pms, const std::vector< std::vector< std::pair > >& dts, INTERP_KERNEL::NormalizedCellType gt, MEDFileFieldPerMesh *father, std::vector > >& extractInfo); + MEDFileFieldPerMeshPerType(MEDFileFieldPerMesh *father, INTERP_KERNEL::NormalizedCellType gt):_father(father),_geo_type(gt) { } private: std::vector addNewEntryIfNecessary(const MEDCouplingFieldDouble *field, int offset, int nbOfCells); std::vector addNewEntryIfNecessaryGauss(const MEDCouplingFieldDouble *field, int offset, int nbOfCells); std::vector addNewEntryIfNecessary(const MEDCouplingFieldDouble *field, const DataArrayInt *subCells); std::vector addNewEntryIfNecessaryGauss(const MEDCouplingFieldDouble *field, const DataArrayInt *subCells); MEDFileFieldPerMeshPerType(med_idt fid, MEDFileFieldPerMesh *fath, TypeOfField type, INTERP_KERNEL::NormalizedCellType geoType, const MEDFileFieldNameScope& nasc, const PartDefinition *pd); - MEDFileFieldPerMeshPerType(MEDFileFieldPerMesh *fath, INTERP_KERNEL::NormalizedCellType geoType); private: MEDFileFieldPerMesh *_father; std::vector< MCAuto > _field_pm_pt_pd; @@ -284,6 +291,7 @@ namespace MEDCoupling void getUndergroundDataArrayExt(std::vector< std::pair,std::pair > >& entries) const; MEDFileFieldPerMeshPerTypePerDisc *getLeafGivenTypeAndLocId(INTERP_KERNEL::NormalizedCellType typ, int locId); const MEDFileFieldPerMeshPerTypePerDisc *getLeafGivenTypeAndLocId(INTERP_KERNEL::NormalizedCellType typ, int locId) const; + static MCAuto Aggregate(int &start, const std::vector& pms, const std::vector< std::vector< std::pair > >& dts, MEDFileAnyTypeField1TSWithoutSDA *father, std::vector > >& extractInfo); private: int addNewEntryIfNecessary(INTERP_KERNEL::NormalizedCellType type); MEDCouplingFieldDouble *finishField(TypeOfField type, const MEDFileFieldGlobsReal *glob, @@ -303,6 +311,7 @@ namespace MEDCoupling static int ComputeNbOfElems(const MEDFileFieldGlobsReal *glob, TypeOfField type, const std::vector& geoTypes, const std::vector< std::pair >& dads, const std::vector& locs); MEDFileFieldPerMesh(med_idt fid, MEDFileAnyTypeField1TSWithoutSDA *fath, int meshCsit, int meshIteration, int meshOrder, const MEDFileFieldNameScope& nasc, const MEDFileMesh *mm, const std::vector< std::pair > *entities); MEDFileFieldPerMesh(MEDFileAnyTypeField1TSWithoutSDA *fath, const MEDCouplingMesh *mesh); + MEDFileFieldPerMesh(MEDFileAnyTypeField1TSWithoutSDA *fath, const std::string& meshName, int meshIt, int meshOrd):_father(fath),_mesh_name(meshName),_mesh_iteration(meshIt),_mesh_order(meshOrd) { } private: std::string _mesh_name; int _mesh_iteration; @@ -562,16 +571,34 @@ namespace MEDCoupling class MEDFileIntField1TSWithoutSDA; + template + class MEDFileField1TSTemplateWithoutSDA : public MEDFileAnyTypeField1TSWithoutSDA + { + protected: + MEDFileField1TSTemplateWithoutSDA(const std::string& fieldName, int csit, int iteration, int order):MEDFileAnyTypeField1TSWithoutSDA(fieldName,csit,iteration,order) { } + MEDFileField1TSTemplateWithoutSDA():MEDFileAnyTypeField1TSWithoutSDA() { } + public: + MEDLOADER_EXPORT void setArray(DataArray *arr); + MEDLOADER_EXPORT DataArray *createNewEmptyDataArrayInstance() const; + MEDLOADER_EXPORT typename Traits::ArrayType *getOrCreateAndGetArrayTemplate(); + MEDLOADER_EXPORT typename Traits::ArrayType const *getOrCreateAndGetArrayTemplate() const; + MEDLOADER_EXPORT typename Traits::ArrayType *getUndergroundDataArrayTemplate() const; + MEDLOADER_EXPORT DataArray *getOrCreateAndGetArray(); + MEDLOADER_EXPORT const DataArray *getOrCreateAndGetArray() const; + MEDLOADER_EXPORT DataArray *getUndergroundDataArray() const; + MEDLOADER_EXPORT void aggregate(const typename std::vector< typename MLFieldTraits::F1TSWSDAType const * >& f1tss, const std::vector< std::vector< std::pair > >& dts); + protected: + MCAuto< typename Traits::ArrayType > _arr; + }; + /*! * SDA is for Shared Data Arrays such as profiles. */ - class MEDFileField1TSWithoutSDA : public MEDFileAnyTypeField1TSWithoutSDA + class MEDFileField1TSWithoutSDA : public MEDFileField1TSTemplateWithoutSDA { public: MEDLOADER_EXPORT const char *getTypeStr() const; - MEDLOADER_EXPORT DataArray *getUndergroundDataArray() const; MEDLOADER_EXPORT DataArray *getUndergroundDataArrayExt(std::vector< std::pair,std::pair > >& entries) const; - MEDLOADER_EXPORT DataArrayDouble *getUndergroundDataArrayDouble() const; MEDLOADER_EXPORT DataArrayDouble *getUndergroundDataArrayDoubleExt(std::vector< std::pair,std::pair > >& entries) const; MEDLOADER_EXPORT std::vector< std::vector > getFieldSplitedByType2(const std::string& mname, std::vector& types, std::vector< std::vector >& typesF, std::vector< std::vector >& pfls, std::vector< std::vector >& locs) const; MEDLOADER_EXPORT static void CheckMeshDimRel(int meshDimRelToMax); @@ -582,15 +609,7 @@ namespace MEDCoupling MEDLOADER_EXPORT MEDFileField1TSWithoutSDA(const std::string& fieldName, int csit, int iteration, int order, const std::vector& infos); MEDLOADER_EXPORT MEDFileAnyTypeField1TSWithoutSDA *shallowCpy() const; MEDLOADER_EXPORT MEDFileAnyTypeField1TSWithoutSDA *deepCopy() const; - MEDLOADER_EXPORT void setArray(DataArray *arr); - MEDLOADER_EXPORT DataArray *createNewEmptyDataArrayInstance() const; - MEDLOADER_EXPORT DataArray *getOrCreateAndGetArray(); - MEDLOADER_EXPORT const DataArray *getOrCreateAndGetArray() const; - MEDLOADER_EXPORT DataArrayDouble *getOrCreateAndGetArrayDouble(); - MEDLOADER_EXPORT const DataArrayDouble *getOrCreateAndGetArrayDouble() const; MEDLOADER_EXPORT MEDFileIntField1TSWithoutSDA *convertToInt() const; - protected: - MCAuto< DataArrayDouble > _arr; public: static const char TYPE_STR[]; }; @@ -598,7 +617,7 @@ namespace MEDCoupling /*! * SDA is for Shared Data Arrays such as profiles. */ - class MEDFileIntField1TSWithoutSDA : public MEDFileAnyTypeField1TSWithoutSDA + class MEDFileIntField1TSWithoutSDA : public MEDFileField1TSTemplateWithoutSDA { public: MEDLOADER_EXPORT MEDFileIntField1TSWithoutSDA(); @@ -606,21 +625,11 @@ namespace MEDCoupling MEDLOADER_EXPORT MEDFileAnyTypeField1TSWithoutSDA *deepCopy() const; MEDLOADER_EXPORT MEDFileAnyTypeField1TSWithoutSDA *shallowCpy() const; MEDLOADER_EXPORT const char *getTypeStr() const; - MEDLOADER_EXPORT DataArray *getUndergroundDataArray() const; MEDLOADER_EXPORT DataArray *getUndergroundDataArrayExt(std::vector< std::pair,std::pair > >& entries) const; - MEDLOADER_EXPORT void setArray(DataArray *arr); - MEDLOADER_EXPORT DataArray *createNewEmptyDataArrayInstance() const; - MEDLOADER_EXPORT DataArray *getOrCreateAndGetArray(); - MEDLOADER_EXPORT const DataArray *getOrCreateAndGetArray() const; - MEDLOADER_EXPORT DataArrayInt *getOrCreateAndGetArrayInt(); - MEDLOADER_EXPORT const DataArrayInt *getOrCreateAndGetArrayInt() const; - MEDLOADER_EXPORT DataArrayInt *getUndergroundDataArrayInt() const; MEDLOADER_EXPORT DataArrayInt *getUndergroundDataArrayIntExt(std::vector< std::pair,std::pair > >& entries) const; MEDLOADER_EXPORT MEDFileField1TSWithoutSDA *convertToDouble() const; protected: MEDFileIntField1TSWithoutSDA(const std::string& fieldName, int csit, int iteration, int order, const std::vector& infos); - protected: - MCAuto< DataArrayInt > _arr; public: MEDLOADER_EXPORT static const char TYPE_STR[]; }; @@ -1000,6 +1009,7 @@ namespace MEDCoupling public: MEDLOADER_EXPORT virtual MEDFileAnyTypeFieldMultiTS *buildNewEmpty() const = 0; MEDLOADER_EXPORT MEDFileAnyTypeFieldMultiTS *extractPart(const std::map >& extractDef, MEDFileMesh *mm) const; + MEDLOADER_EXPORT static MCAuto Aggregate(const std::vector& fmtss, const std::vector< std::vector< std::pair > >& dts); public: MEDLOADER_EXPORT std::vector getPflsReallyUsed() const; MEDLOADER_EXPORT std::vector getLocsReallyUsed() const; @@ -1076,7 +1086,7 @@ namespace MEDCoupling MEDLOADER_EXPORT static MEDFileIntFieldMultiTS *LoadSpecificEntities(const std::string& fileName, const std::string& fieldName, const std::vector< std::pair >& entities, bool loadAll=true); MEDLOADER_EXPORT MEDFileAnyTypeFieldMultiTS *shallowCpy() const; MEDLOADER_EXPORT void checkCoherencyOfType(const MEDFileAnyTypeField1TS *f1ts) const; - MEDLOADER_EXPORT MEDFileAnyTypeField1TS *getTimeStepAtPos(int pos) const; + MEDLOADER_EXPORT MEDFileIntField1TS *getTimeStepAtPos(int pos) const; MEDLOADER_EXPORT MEDFileFieldMultiTS *convertToDouble(bool isDeepCpyGlobs=true) const; // MEDLOADER_EXPORT MEDCouplingFieldInt *field(int iteration, int order, const MEDFileMesh *mesh) const; diff --git a/src/MEDLoader/MEDFileField.txx b/src/MEDLoader/MEDFileField.txx new file mode 100644 index 000000000..ae2b3a93a --- /dev/null +++ b/src/MEDLoader/MEDFileField.txx @@ -0,0 +1,172 @@ +// Copyright (C) 2007-2016 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// 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 +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay (EDF R&D) + +#ifndef __MEDFILEFIELD_TXX__ +#define __MEDFILEFIELD_TXX__ + +#include "MEDFileField.hxx" +#include "MEDCouplingTraits.hxx" + +namespace MEDCoupling +{ + template + void MEDFileField1TSTemplateWithoutSDA::setArray(DataArray *arr) + { + if(!arr) + { + _nb_of_tuples_to_be_allocated=-1; + _arr=0; + return ; + } + typename Traits::ArrayType *arrC=dynamic_cast::ArrayType *>(arr); + if(!arrC) + throw INTERP_KERNEL::Exception("MEDFileField1TSTemplateWithoutSDA::setArray : the input not null array is not of type DataArrayDouble !"); + else + _nb_of_tuples_to_be_allocated=-3; + arrC->incrRef(); + _arr=arrC; + } + + /*! + * Returns a pointer to the underground DataArrayTemplate instance. So the + * caller should not decrRef() it. This method allows for a direct access to the field + * values. This method is quite unusable if there is more than a nodal field or a cell + * field on single geometric cell type. + * \return DataArrayTemplate * - the pointer to the field values array. + */ + template + DataArray *MEDFileField1TSTemplateWithoutSDA::getOrCreateAndGetArray() + { + return getOrCreateAndGetArrayTemplate(); + } + + template + const DataArray *MEDFileField1TSTemplateWithoutSDA::getOrCreateAndGetArray() const + { + return getOrCreateAndGetArrayTemplate(); + } + + template + DataArray *MEDFileField1TSTemplateWithoutSDA::createNewEmptyDataArrayInstance() const + { + return Traits::ArrayType::New(); + } + + template + typename Traits::ArrayType const *MEDFileField1TSTemplateWithoutSDA::getOrCreateAndGetArrayTemplate() const + { + typename Traits::ArrayType const *ret(_arr); + if(ret) + return ret; + (const_cast< MEDFileField1TSTemplateWithoutSDA *>(this))->_arr=Traits::ArrayType::New(); + return _arr; + } + + template + typename Traits::ArrayType *MEDFileField1TSTemplateWithoutSDA::getOrCreateAndGetArrayTemplate() + { + typename Traits::ArrayType *ret(_arr); + if(ret) + return ret; + _arr=Traits::ArrayType::New(); + return _arr; + } + + /*! + * Returns a pointer to the underground DataArrayTemplate instance. So the + * caller should not decrRef() it. This method allows for a direct access to the field + * values. This method is quite unusable if there is more than a nodal field or a cell + * field on single geometric cell type. + * \return DataArrayTemplate * - the pointer to the field values array. + */ + template + typename Traits::ArrayType *MEDFileField1TSTemplateWithoutSDA::getUndergroundDataArrayTemplate() const + { + typename Traits::ArrayType const *ret(_arr); + if(ret) + return const_cast::ArrayType *>(ret); + else + return 0; + } + + /*! + * Returns a pointer to the underground DataArrayDouble instance. So the + * caller should not decrRef() it. This method allows for a direct access to the field + * values. This method is quite unusable if there is more than a nodal field or a cell + * field on single geometric cell type. + * \return DataArrayDouble * - the pointer to the field values array. + */ + template + DataArray *MEDFileField1TSTemplateWithoutSDA::getUndergroundDataArray() const + { + return getUndergroundDataArrayTemplate(); + } + + template + void MEDFileField1TSTemplateWithoutSDA::aggregate(const std::vector::F1TSWSDAType const *>& f1tss, const std::vector< std::vector< std::pair > >& dts) + { + if(f1tss.empty()) + throw INTERP_KERNEL::Exception("MEDFileField1TSTemplateWithoutSDA::aggregate : empty vector !"); + std::size_t sz(f1tss.size()),ii(0); + std::vector pms; + std::vector das(sz); + for(typename std::vector::F1TSWSDAType const *>::const_iterator it=f1tss.begin();it!=f1tss.end();it++,ii++) + { + if(!*it) + throw INTERP_KERNEL::Exception("MEDFileField1TSTemplateWithoutSDA::aggregate : presence of null pointer in input vector !"); + if((*it)->_field_per_mesh.empty()) + throw INTERP_KERNEL::Exception("MEDFileField1TSTemplateWithoutSDA::aggregate : no info !"); + const typename Traits::ArrayType *arr((*it)->getUndergroundDataArrayTemplate()); + if(!arr) + throw INTERP_KERNEL::Exception("MEDFileField1TSTemplateWithoutSDA::aggregate : presence of null array !"); + das[ii]=arr; + pms.push_back((*it)->_field_per_mesh[0]); + } + typename MLFieldTraits::F1TSWSDAType const *refPt(f1tss[0]); + setName(refPt->getName()); + + const DataArray *arr(refPt->getUndergroundDataArray()); + int nbCompo(arr->getNumberOfComponents()); + for(typename std::vector::F1TSWSDAType const *>::const_iterator it=f1tss.begin();it!=f1tss.end();it++) + { + const typename Traits::ArrayType *myArr((*it)->getUndergroundDataArrayTemplate()); + if(myArr->getNumberOfComponents()!=nbCompo) + throw INTERP_KERNEL::Exception("MEDFileField1TSTemplateWithoutSDA::aggregate : arrays must have same number of components !"); + } + std::vector > > extractInfo; + int start(0); + MCAuto fpm(MEDFileFieldPerMesh::Aggregate(start,pms,dts,this,extractInfo)); + _field_per_mesh.push_back(fpm); + int iteration,order; + double tv(f1tss[0]->getTime(iteration,order)); + _iteration=iteration; _order=order; _dt=tv; + _arr=Traits::ArrayType::New(); + _arr->alloc(start,nbCompo); _arr->copyStringInfoFrom(*arr); + start=0; + for(std::vector > >::const_iterator it=extractInfo.begin();it!=extractInfo.end();it++) + { + const DataArray *zeArr(das[(*it).first]); + _arr->setContigPartOfSelectedValuesSlice(start,zeArr,(*it).second.first,(*it).second.second,1); + start+=(*it).second.second-(*it).second.first; + } + } +} + +#endif diff --git a/src/MEDLoader/MEDFileMesh.cxx b/src/MEDLoader/MEDFileMesh.cxx index d187a455e..0332e42d6 100644 --- a/src/MEDLoader/MEDFileMesh.cxx +++ b/src/MEDLoader/MEDFileMesh.cxx @@ -2090,6 +2090,9 @@ std::vector MEDFileMesh::getAllGeoTypes() con return ret; } +/*! + * \sa getAllDistributionOfTypes + */ std::vector MEDFileMesh::getDistributionOfTypes(int meshDimRelToMax) const { MCAuto mLev(getMeshAtLevel(meshDimRelToMax)); @@ -2335,7 +2338,7 @@ std::vector MEDFileUMesh::getDirectChildrenWithNull() c return ret; } -MEDFileMesh *MEDFileUMesh::shallowCpy() const +MEDFileUMesh *MEDFileUMesh::shallowCpy() const { MCAuto ret(new MEDFileUMesh(*this)); return ret.retn(); @@ -2346,7 +2349,7 @@ MEDFileMesh *MEDFileUMesh::createNewEmpty() const return new MEDFileUMesh; } -MEDFileMesh *MEDFileUMesh::deepCopy() const +MEDFileUMesh *MEDFileUMesh::deepCopy() const { MCAuto ret(new MEDFileUMesh(*this)); ret->deepCpyEquivalences(*this); @@ -3588,6 +3591,30 @@ MEDCoupling1GTUMesh *MEDFileUMesh::getDirectUndergroundSingleGeoTypeMesh(INTERP_ return sp->getDirectUndergroundSingleGeoTypeMesh(gt); } +/*! + * This method returns for each geo types in \a this number of cells with this geo type. + * This method returns info as a vector of pair. The first element of pair is geo type and the second the number of cells associated. + * This method also returns the number of nodes of \a this (key associated is NORM_ERROR) + * + * \sa getDistributionOfTypes + */ +std::vector< std::pair > MEDFileUMesh::getAllDistributionOfTypes() const +{ + std::vector< std::pair > ret; + std::vector nel(getNonEmptyLevels()); + for(std::vector::reverse_iterator it=nel.rbegin();it!=nel.rend();it++) + { + std::vector gt(getGeoTypesAtLevel(*it)); + for(std::vector::const_iterator it1=gt.begin();it1!=gt.end();it1++) + { + int nbCells(getNumberOfCellsWithType(*it1)); + ret.push_back(std::pair(*it1,nbCells)); + } + } + ret.push_back(std::pair(INTERP_KERNEL::NORM_ERROR,getNumberOfNodes())); + return ret; +} + /*! * Given a relative level \a meshDimRelToMax it returns the sorted vector of geometric types present in \a this. * \throw if the reqsuested \a meshDimRelToMax does not exist. @@ -3703,12 +3730,42 @@ void MEDFileUMesh::setCoords(DataArrayDouble *coords) if(coords==(DataArrayDouble *)_coords) return ; coords->checkAllocated(); - int nbOfTuples=coords->getNumberOfTuples(); + int nbOfTuples(coords->getNumberOfTuples()); _coords=coords; coords->incrRef(); _fam_coords=DataArrayInt::New(); _fam_coords->alloc(nbOfTuples,1); _fam_coords->fillWithZero(); + _num_coords=0; _rev_num_coords=0; _name_coords=0; + for(std::vector< MCAuto >::iterator it=_ms.begin();it!=_ms.end();it++) + if((MEDFileUMeshSplitL1 *)(*it)) + (*it)->setCoords(coords); +} + +/*! + * Change coords without changing anything concerning families and numbering on nodes. + */ +void MEDFileUMesh::setCoordsForced(DataArrayDouble *coords) +{ + if(!coords) + throw INTERP_KERNEL::Exception("MEDFileUMesh::setCoordsForced : null pointer in input !"); + if(coords==(DataArrayDouble *)_coords) + return ; + coords->checkAllocated(); + int nbOfTuples(coords->getNumberOfTuples()); + if(_coords.isNull()) + { + _coords=coords; + coords->incrRef(); + } + else + { + int oldNbTuples(_coords->getNumberOfTuples()); + if(oldNbTuples!=nbOfTuples) + throw INTERP_KERNEL::Exception("MEDFileUMesh::setCoordsForced : number of tuples is not the same -> invoke setCoords instead !"); + _coords=coords; + coords->incrRef(); + } for(std::vector< MCAuto >::iterator it=_ms.begin();it!=_ms.end();it++) if((MEDFileUMeshSplitL1 *)(*it)) (*it)->setCoords(coords); @@ -4453,6 +4510,120 @@ MEDFileUMesh *MEDFileUMesh::quadraticToLinear(double eps) const return ret.retn(); } +/*! + * Computes the symmetry of \a this. + * \return a new object. + */ +MCAuto MEDFileUMesh::symmetry3DPlane(const double point[3], const double normalVector[3]) const +{ + MCAuto ret(deepCopy()); + DataArrayDouble *myCoo(getCoords()); + if(myCoo) + { + MCAuto newCoo(myCoo->symmetry3DPlane(point,normalVector)); + ret->setCoordsForced(newCoo); + } + return ret; +} + +MCAuto MEDFileUMesh::Aggregate(const std::vector& meshes) +{ + if(meshes.empty()) + throw INTERP_KERNEL::Exception("MEDFileUMesh::Aggregate : empty input vector !"); + std::size_t sz(meshes.size()),i(0); + std::vector coos(sz); + std::vector fam_coos(sz),num_coos(sz); + for(std::vector::const_iterator it=meshes.begin();it!=meshes.end();it++,i++) + { + if(!(*it)) + throw INTERP_KERNEL::Exception("MEDFileUMesh::Aggregate : presence of NULL pointer in input vector !"); + coos[i]=(*it)->getCoords(); + fam_coos[i]=(*it)->getFamilyFieldAtLevel(1); + num_coos[i]=(*it)->getNumberFieldAtLevel(1); + } + const MEDFileUMesh *ref(meshes[0]); + int spaceDim(ref->getSpaceDimension()),meshDim(ref->getMeshDimension()); + std::vector levs(ref->getNonEmptyLevels()); + std::map > m_fam,m_renum; + std::map > > m_mesh2; + std::map > m_mesh; + std::map map1; + std::map > map2; + for(std::vector::const_iterator it=meshes.begin();it!=meshes.end();it++,i++) + { + if((*it)->getSpaceDimension()!=spaceDim) + throw INTERP_KERNEL::Exception("MEDFileUMesh::Aggregate : space dimension must be homogeneous !"); + if((*it)->getMeshDimension()!=meshDim) + throw INTERP_KERNEL::Exception("MEDFileUMesh::Aggregate : mesh dimension must be homogeneous !"); + if((*it)->getNonEmptyLevels()!=levs) + throw INTERP_KERNEL::Exception("MEDFileUMesh::Aggregate : levels must be the same for elements in input vector !"); + for(std::vector::const_iterator it2=levs.begin();it2!=levs.end();it2++) + { + MCAuto locMesh((*it)->getMeshAtLevel(*it2)); + m_mesh[*it2].push_back(locMesh); m_mesh2[*it2].push_back(locMesh); + m_fam[*it2].push_back((*it)->getFamilyFieldAtLevel(*it2)); + m_renum[*it2].push_back((*it)->getNumberFieldAtLevel(*it2)); + } + const std::map& locMap1((*it)->getFamilyInfo()); + for(std::map::const_iterator it3=locMap1.begin();it3!=locMap1.end();it3++) + map1[(*it3).first]=(*it3).second; + const std::map >& locMap2((*it)->getGroupInfo()); + for(std::map >::const_iterator it4=locMap2.begin();it4!=locMap2.end();it4++) + map2[(*it4).first]=(*it4).second; + } + // Easy part : nodes + MCAuto ret(MEDFileUMesh::New()); + MCAuto coo(DataArrayDouble::Aggregate(coos)); + ret->setCoords(coo); + if(std::find(fam_coos.begin(),fam_coos.end(),(const DataArrayInt *)0)==fam_coos.end()) + { + MCAuto fam_coo(DataArrayInt::Aggregate(fam_coos)); + ret->setFamilyFieldArr(1,fam_coo); + } + if(std::find(num_coos.begin(),num_coos.end(),(const DataArrayInt *)0)==num_coos.end()) + { + MCAuto num_coo(DataArrayInt::Aggregate(num_coos)); + ret->setRenumFieldArr(1,num_coo); + } + // cells + for(std::vector::const_iterator it=levs.begin();it!=levs.end();it++) + { + std::map >::const_iterator it2(m_mesh.find(*it)); + if(it2==m_mesh.end()) + throw INTERP_KERNEL::Exception("MEDFileUMesh::Aggregate : internal error 1 !"); + MCAuto mesh(MEDCouplingUMesh::MergeUMeshes((*it2).second)); + mesh->setCoords(coo); mesh->setName(ref->getName()); + MCAuto renum(mesh->sortCellsInMEDFileFrmt()); + ret->setMeshAtLevel(*it,mesh); + std::map >::const_iterator it3(m_fam.find(*it)),it4(m_renum.find(*it)); + if(it3!=m_fam.end()) + { + const std::vector& fams((*it3).second); + if(std::find(fams.begin(),fams.end(),(const DataArrayInt *)0)==fams.end()) + { + MCAuto famm(DataArrayInt::Aggregate(fams)); + famm->renumberInPlace(renum->begin()); + ret->setFamilyFieldArr(*it,famm); + } + } + if(it4!=m_renum.end()) + { + const std::vector& renums((*it4).second); + if(std::find(renums.begin(),renums.end(),(const DataArrayInt *)0)==renums.end()) + { + MCAuto renumm(DataArrayInt::Aggregate(renums)); + renumm->renumberInPlace(renum->begin()); + ret->setRenumFieldArr(*it,renumm); + } + } + } + // + ret->setFamilyInfo(map1); + ret->setGroupInfo(map2); + ret->setName(ref->getName()); + return ret; +} + void MEDFileUMesh::serialize(std::vector& tinyDouble, std::vector& tinyInt, std::vector& tinyStr, std::vector< MCAuto >& bigArraysI, MCAuto& bigArrayD) { clearNonDiscrAttributes(); @@ -6305,7 +6476,7 @@ std::string MEDFileCMesh::advancedRepr() const return simpleRepr(); } -MEDFileMesh *MEDFileCMesh::shallowCpy() const +MEDFileCMesh *MEDFileCMesh::shallowCpy() const { MCAuto ret(new MEDFileCMesh(*this)); return ret.retn(); @@ -6316,7 +6487,7 @@ MEDFileMesh *MEDFileCMesh::createNewEmpty() const return new MEDFileCMesh; } -MEDFileMesh *MEDFileCMesh::deepCopy() const +MEDFileCMesh *MEDFileCMesh::deepCopy() const { MCAuto ret(new MEDFileCMesh(*this)); ret->deepCpyEquivalences(*this); @@ -6547,7 +6718,7 @@ std::vector MEDFileCurveLinearMesh::getDirectChildrenWi return ret; } -MEDFileMesh *MEDFileCurveLinearMesh::shallowCpy() const +MEDFileCurveLinearMesh *MEDFileCurveLinearMesh::shallowCpy() const { MCAuto ret(new MEDFileCurveLinearMesh(*this)); return ret.retn(); @@ -6558,7 +6729,7 @@ MEDFileMesh *MEDFileCurveLinearMesh::createNewEmpty() const return new MEDFileCurveLinearMesh; } -MEDFileMesh *MEDFileCurveLinearMesh::deepCopy() const +MEDFileCurveLinearMesh *MEDFileCurveLinearMesh::deepCopy() const { MCAuto ret(new MEDFileCurveLinearMesh(*this)); ret->deepCpyEquivalences(*this); diff --git a/src/MEDLoader/MEDFileMesh.hxx b/src/MEDLoader/MEDFileMesh.hxx index 9f0231313..1e0f48bb2 100644 --- a/src/MEDLoader/MEDFileMesh.hxx +++ b/src/MEDLoader/MEDFileMesh.hxx @@ -252,8 +252,8 @@ namespace MEDCoupling MEDLOADER_EXPORT std::size_t getHeapMemorySizeWithoutChildren() const; MEDLOADER_EXPORT std::vector getDirectChildrenWithNull() const; MEDLOADER_EXPORT MEDFileMesh *createNewEmpty() const; - MEDLOADER_EXPORT MEDFileMesh *deepCopy() const; - MEDLOADER_EXPORT MEDFileMesh *shallowCpy() const; + MEDLOADER_EXPORT MEDFileUMesh *deepCopy() const; + MEDLOADER_EXPORT MEDFileUMesh *shallowCpy() const; MEDLOADER_EXPORT bool isEqual(const MEDFileMesh *other, double eps, std::string& what) const; MEDLOADER_EXPORT void checkConsistency() const; MEDLOADER_EXPORT void checkSMESHConsistency() const; @@ -299,6 +299,7 @@ namespace MEDCoupling MEDLOADER_EXPORT DataArrayInt *getFamiliesArr(int meshDimRelToMaxExt, const std::vector& fams, bool renum=false) const; MEDLOADER_EXPORT MEDCouplingUMesh *getMeshAtLevel(int meshDimRelToMax, bool renum=false) const; MEDLOADER_EXPORT std::vector getDistributionOfTypes(int meshDimRelToMax) const; + MEDLOADER_EXPORT std::vector< std::pair > getAllDistributionOfTypes() const; MEDLOADER_EXPORT MEDCouplingUMesh *getLevel0Mesh(bool renum=false) const; MEDLOADER_EXPORT MEDCouplingUMesh *getLevelM1Mesh(bool renum=false) const; MEDLOADER_EXPORT MEDCouplingUMesh *getLevelM2Mesh(bool renum=false) const; @@ -312,6 +313,7 @@ namespace MEDCoupling // MEDLOADER_EXPORT void setFamilyNameAttachedOnId(int id, const std::string& newFamName); MEDLOADER_EXPORT void setCoords(DataArrayDouble *coords); + MEDLOADER_EXPORT void setCoordsForced(DataArrayDouble *coords); MEDLOADER_EXPORT void eraseGroupsAtLevel(int meshDimRelToMaxExt); MEDLOADER_EXPORT void setFamilyFieldArr(int meshDimRelToMaxExt, DataArrayInt *famArr); MEDLOADER_EXPORT void setRenumFieldArr(int meshDimRelToMaxExt, DataArrayInt *renumArr); @@ -334,6 +336,8 @@ namespace MEDCoupling MEDLOADER_EXPORT MEDFileUMesh *buildExtrudedMesh(const MEDCouplingUMesh *m1D, int policy) const; MEDLOADER_EXPORT MEDFileUMesh *linearToQuadratic(int conversionType=0, double eps=1e-12) const; MEDLOADER_EXPORT MEDFileUMesh *quadraticToLinear(double eps=1e-12) const; + MEDLOADER_EXPORT MCAuto symmetry3DPlane(const double point[3], const double normalVector[3]) const; + MEDLOADER_EXPORT static MCAuto Aggregate(const std::vector& meshes); // serialization MEDLOADER_EXPORT void serialize(std::vector& tinyDouble, std::vector& tinyInt, std::vector& tinyStr, std::vector< MCAuto >& bigArraysI, MCAuto& bigArrayD); @@ -447,8 +451,8 @@ namespace MEDCoupling MEDLOADER_EXPORT std::size_t getHeapMemorySizeWithoutChildren() const; MEDLOADER_EXPORT std::vector getDirectChildrenWithNull() const; MEDLOADER_EXPORT MEDFileMesh *createNewEmpty() const; - MEDLOADER_EXPORT MEDFileMesh *deepCopy() const; - MEDLOADER_EXPORT MEDFileMesh *shallowCpy() const; + MEDLOADER_EXPORT MEDFileCMesh *deepCopy() const; + MEDLOADER_EXPORT MEDFileCMesh *shallowCpy() const; MEDLOADER_EXPORT bool isEqual(const MEDFileMesh *other, double eps, std::string& what) const; MEDLOADER_EXPORT int getMeshDimension() const; MEDLOADER_EXPORT int getSpaceDimension() const; @@ -480,8 +484,8 @@ namespace MEDCoupling MEDLOADER_EXPORT std::size_t getHeapMemorySizeWithoutChildren() const; MEDLOADER_EXPORT std::vector getDirectChildrenWithNull() const; MEDLOADER_EXPORT MEDFileMesh *createNewEmpty() const; - MEDLOADER_EXPORT MEDFileMesh *deepCopy() const; - MEDLOADER_EXPORT MEDFileMesh *shallowCpy() const; + MEDLOADER_EXPORT MEDFileCurveLinearMesh *deepCopy() const; + MEDLOADER_EXPORT MEDFileCurveLinearMesh *shallowCpy() const; MEDLOADER_EXPORT bool isEqual(const MEDFileMesh *other, double eps, std::string& what) const; MEDLOADER_EXPORT int getMeshDimension() const; MEDLOADER_EXPORT std::string simpleRepr() const; diff --git a/src/MEDLoader/MEDLoaderTraits.hxx b/src/MEDLoader/MEDLoaderTraits.hxx new file mode 100644 index 000000000..25cdfb3ee --- /dev/null +++ b/src/MEDLoader/MEDLoaderTraits.hxx @@ -0,0 +1,58 @@ +// Copyright (C) 2016 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// 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 +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay (EDF R&D) + +#ifndef __MEDLOADERTRAITS_HXX__ +#define __MEDLOADERTRAITS_HXX__ + +#include "MEDLoaderDefines.hxx" + +namespace MEDCoupling +{ + template + struct MEDLOADER_EXPORT MLFieldTraits + { + typedef T EltType; + }; + + class MEDFileFieldMultiTS; + class MEDFileField1TS; + class MEDFileIntFieldMultiTS; + class MEDFileIntField1TS; + class MEDFileField1TSWithoutSDA; + class MEDFileIntField1TSWithoutSDA; + + template<> + struct MEDLOADER_EXPORT MLFieldTraits + { + typedef MEDFileFieldMultiTS FMTSType; + typedef MEDFileField1TS F1TSType; + typedef MEDFileField1TSWithoutSDA F1TSWSDAType; + }; + + template<> + struct MEDLOADER_EXPORT MLFieldTraits + { + typedef MEDFileIntFieldMultiTS FMTSType; + typedef MEDFileIntField1TS F1TSType; + typedef MEDFileIntField1TSWithoutSDA F1TSWSDAType; + }; +} + +#endif diff --git a/src/MEDLoader/Swig/MEDLoaderCommon.i b/src/MEDLoader/Swig/MEDLoaderCommon.i index 818147457..3a348ebde 100644 --- a/src/MEDLoader/Swig/MEDLoaderCommon.i +++ b/src/MEDLoader/Swig/MEDLoaderCommon.i @@ -125,6 +125,8 @@ using namespace MEDCoupling; %newobject MEDCoupling::MEDFileUMesh::buildExtrudedMesh; %newobject MEDCoupling::MEDFileUMesh::linearToQuadratic; %newobject MEDCoupling::MEDFileUMesh::quadraticToLinear; +%newobject MEDCoupling::MEDFileUMesh::symmetry3DPlane; +%newobject MEDCoupling::MEDFileUMesh::Aggregate; %newobject MEDCoupling::MEDFileCMesh::New; %newobject MEDCoupling::MEDFileCurveLinearMesh::New; %newobject MEDCoupling::MEDFileMeshMultiTS::New; @@ -205,6 +207,7 @@ using namespace MEDCoupling; %newobject MEDCoupling::MEDFileData::getMeshes; %newobject MEDCoupling::MEDFileData::getFields; %newobject MEDCoupling::MEDFileData::getParams; +%newobject MEDCoupling::MEDFileData::Aggregate; %newobject MEDCoupling::MEDFileParameterDouble1TS::New; %newobject MEDCoupling::MEDFileParameterDouble1TS::deepCopy; @@ -1241,6 +1244,7 @@ namespace MEDCoupling // void setFamilyNameAttachedOnId(int id, const std::string& newFamName) throw(INTERP_KERNEL::Exception); void setCoords(DataArrayDouble *coords) throw(INTERP_KERNEL::Exception); + void setCoordsForced(DataArrayDouble *coords) throw(INTERP_KERNEL::Exception); void eraseGroupsAtLevel(int meshDimRelToMaxExt) throw(INTERP_KERNEL::Exception); void removeMeshAtLevel(int meshDimRelToMax) throw(INTERP_KERNEL::Exception); void setMeshAtLevel(int meshDimRelToMax, MEDCoupling1GTUMesh *m) throw(INTERP_KERNEL::Exception); @@ -1394,6 +1398,34 @@ namespace MEDCoupling self->removeMeshAtLevel(meshDimRelToMax); } + MEDFileUMesh *symmetry3DPlane(PyObject *point, PyObject *normalVector) const throw(INTERP_KERNEL::Exception) + { + const char msg[]="Python wrap of MEDFileUMesh::symmetry3DPlane : "; + double val,val2; + DataArrayDouble *a,*a2; + DataArrayDoubleTuple *aa,*aa2; + std::vector bb,bb2; + int sw; + const double *centerPtr(convertObjToPossibleCpp5_Safe(point,sw,val,a,aa,bb,msg,1,3,true)); + const double *vectorPtr(convertObjToPossibleCpp5_Safe(normalVector,sw,val2,a2,aa2,bb2,msg,1,3,true)); + MCAuto ret(self->symmetry3DPlane(centerPtr,vectorPtr)); + return ret.retn(); + } + + static MEDFileUMesh *Aggregate(PyObject *meshes) throw(INTERP_KERNEL::Exception) + { + std::vector meshesCpp; + convertFromPyObjVectorOfObj(meshes,SWIGTYPE_p_MEDCoupling__MEDFileUMesh,"MEDFileUMesh",meshesCpp); + MCAuto ret(MEDFileUMesh::Aggregate(meshesCpp)); + return ret.retn(); + } + + PyObject *getAllDistributionOfTypes() const throw(INTERP_KERNEL::Exception) + { + std::vector< std::pair > ret(self->getAllDistributionOfTypes()); + return convertVecPairIntToPy(ret); + } + DataArrayInt *deduceNodeSubPartFromCellSubPart(PyObject *extractDef) const throw(INTERP_KERNEL::Exception) { std::map > extractDefCpp; @@ -2304,17 +2336,8 @@ namespace MEDCoupling PyObject *getIterations() const throw(INTERP_KERNEL::Exception) { - std::vector< std::pair > res=self->getIterations(); - PyObject *ret=PyList_New(res.size()); - int rk=0; - for(std::vector< std::pair >::const_iterator iter=res.begin();iter!=res.end();iter++,rk++) - { - PyObject *elt=PyTuple_New(2); - PyTuple_SetItem(elt,0,SWIG_From_int((*iter).first)); - PyTuple_SetItem(elt,1,SWIG_From_int((*iter).second)); - PyList_SetItem(ret,rk,elt); - } - return ret; + std::vector< std::pair > res(self->getIterations()); + return convertVecPairIntToPy(res); } PyObject *getTimeSteps() const throw(INTERP_KERNEL::Exception) @@ -3460,6 +3483,14 @@ namespace MEDCoupling std::vector< std::pair > modifTab=convertVecPairStStFromPy(li); return self->changeMeshNames(modifTab); } + + static MEDFileData *Aggregate(PyObject *mfds) throw(INTERP_KERNEL::Exception) + { + std::vector mfdsCpp; + convertFromPyObjVectorOfObj(mfds,SWIGTYPE_p_MEDCoupling__MEDFileData,"MEDFileData",mfdsCpp); + MCAuto ret(MEDFileData::Aggregate(mfdsCpp)); + return ret.retn(); + } } }; diff --git a/src/MEDLoader/Swig/MEDLoaderTest3.py b/src/MEDLoader/Swig/MEDLoaderTest3.py index 9ec3797c1..7dc76b2a2 100644 --- a/src/MEDLoader/Swig/MEDLoaderTest3.py +++ b/src/MEDLoader/Swig/MEDLoaderTest3.py @@ -5527,6 +5527,222 @@ class MEDLoaderTest3(unittest.TestCase): self.assertEqual(fPart.getTime(),list(tt)) pass pass + + def test42(self): + """ Testing of MEDFileData::Aggregate and MEDFileUMesh::Aggregate and MEDFileUMesh::getAllDistributionOfType """ + fname1="Pyfile106_1.med" + fname2="Pyfile106_2.med" + fname3="Pyfile106_3.med" + meshName="mesh" + mm1=MEDFileUMesh() + da1=DataArrayDouble([1,2,10,3,4,11,5,6,12,7,8,13],4,3) ; da1.setInfoOnComponents(["aa [m]","bbb [kg]","cccc [MW]"]) + mm1.setCoords(da1) + mm1_0=MEDCouplingUMesh(meshName,3) ; mm1_0.allocateCells() + mm1_0.setCoords(da1) + mm1_0.insertNextCell(NORM_TETRA4,[0,1,2,3]) + mm1_0.insertNextCell(NORM_TETRA4,[4,5,6,7]) + mm1_0.insertNextCell(NORM_PENTA6,[8,9,10,11,12,13]) + mm1_0.insertNextCell(NORM_PENTA6,[14,15,16,17,18,19]) + mm1_0.insertNextCell(NORM_PENTA6,[20,21,22,23,24,25]) + mm1[0]=mm1_0 + mm1.setFamilyFieldArr(0,DataArrayInt([1,2,3,4,5])) + mm1.setRenumFieldArr(0,DataArrayInt([11,12,13,14,15])) + # + mm1_1=MEDCouplingUMesh(meshName,2) ; mm1_1.allocateCells() + mm1_1.setCoords(da1) + mm1_1.insertNextCell(NORM_TRI3,[0,1,2]) + mm1_1.insertNextCell(NORM_TRI3,[3,4,5]) + mm1_1.insertNextCell(NORM_QUAD4,[6,7,8,9]) + mm1_1.insertNextCell(NORM_QUAD4,[10,11,12,13]) + mm1_1.insertNextCell(NORM_QUAD4,[14,15,16,17]) + mm1_1.insertNextCell(NORM_QUAD4,[18,19,20,21]) + mm1[-1]=mm1_1 + mm1.setFamilyFieldArr(-1,DataArrayInt([6,7,8,9,10,11])) + mm1.setRenumFieldArr(-1,DataArrayInt([16,17,18,19,20,21])) + for i in range(1,10): + mm1.setFamilyId("F%d"%i,i) + mm1.setFamilyId("FAMILLE_ZERO",0) + mm1.setFamilyId("H1",100) + mm1.setFamiliesOnGroup("myGRP",["F2","F6"]) + mm1.setFamiliesOnGroup("myGRP1",["F2","F6"]) + mm1.setFamilyFieldArr(1,DataArrayInt([12,13,14,15])) + mm1.setRenumFieldArr(1,DataArrayInt([22,23,24,25])) + ############## + mm2=MEDFileUMesh() + da1=DataArrayDouble([9,10,30,11,12,31,13,14,32,15,16,33,17,18,34],5,3) ; da1.setInfoOnComponents(["aa [m]","bbb [kg]","cccc [MW]"]) + mm2.setCoords(da1) + mm2_0=MEDCouplingUMesh(meshName,3) ; mm2_0.allocateCells() + mm2_0.setCoords(da1) + mm2_0.insertNextCell(NORM_TETRA4,[100,101,102,103]) + mm2_0.insertNextCell(NORM_TETRA4,[104,105,106,107]) + mm2_0.insertNextCell(NORM_TETRA4,[108,109,110,111]) + mm2_0.insertNextCell(NORM_PENTA6,[112,113,114,115,116,117]) + mm2[0]=mm2_0 + mm2.setFamilyFieldArr(0,DataArrayInt([40,41,42,43])) + mm2.setRenumFieldArr(0,DataArrayInt([50,51,52,53])) + # + mm2_1=MEDCouplingUMesh(meshName,2) ; mm2_1.allocateCells() + mm2_1.setCoords(da1) + mm2_1.insertNextCell(NORM_TRI3,[100,101,102]) + mm2_1.insertNextCell(NORM_TRI3,[103,104,105]) + mm2_1.insertNextCell(NORM_TRI3,[106,107,108]) + mm2_1.insertNextCell(NORM_QUAD4,[109,110,111,112]) + mm2_1.insertNextCell(NORM_QUAD4,[113,114,115,116]) + mm2_1.insertNextCell(NORM_QUAD4,[117,118,119,120]) + mm2_1.insertNextCell(NORM_QUAD4,[121,122,123,124]) + mm2_1.insertNextCell(NORM_QUAD4,[125,126,127,128]) + mm2[-1]=mm2_1 + mm2.setFamilyFieldArr(-1,DataArrayInt([200,201,202,203,204,205,206,207])) + mm2.setRenumFieldArr(-1,DataArrayInt([300,301,302,303,304,305,306,307])) + for i in range(1,12): + mm2.setFamilyId("G%d"%i,i+30) + mm2.setFamilyId("H1",100) + mm2.setFamilyId("FAMILLE_ZERO",0) + mm2.setFamiliesOnGroup("myGRP",["G2","G6"]) + mm2.setFamiliesOnGroup("myGRP2",["G4","G7"]) + mm2.setFamilyFieldArr(1,DataArrayInt([112,113,114,115,116])) + mm2.setRenumFieldArr(1,DataArrayInt([122,123,124,125,126])) + # + mm=MEDFileUMesh.Aggregate([mm1,mm2]) + ####### + def CheckMesh(tester,mm): + cooExp=DataArrayDouble([(1,2,10),(3,4,11),(5,6,12),(7,8,13),(9,10,30),(11,12,31),(13,14,32),(15,16,33),(17,18,34)]) ; cooExp.setInfoOnComponents(["aa [m]","bbb [kg]","cccc [MW]"]) + tester.assertTrue(mm.getCoords().isEqual(cooExp,1e-12)) + tester.assertTrue(mm[0].getNodalConnectivity().isEqual(DataArrayInt([14,0,1,2,3,14,4,5,6,7,14,104,105,106,107,14,108,109,110,111,14,112,113,114,115,16,8,9,10,11,12,13,16,14,15,16,17,18,19,16,20,21,22,23,24,25,16,116,117,118,119,120,121]))) + tester.assertTrue(mm[0].getNodalConnectivityIndex().isEqual(DataArrayInt([0,5,10,15,20,25,32,39,46,53]))) + tester.assertTrue(mm[-1].getNodalConnectivity().isEqual(DataArrayInt([3,0,1,2,3,3,4,5,3,104,105,106,3,107,108,109,3,110,111,112,4,6,7,8,9,4,10,11,12,13,4,14,15,16,17,4,18,19,20,21,4,113,114,115,116,4,117,118,119,120,4,121,122,123,124,4,125,126,127,128,4,129,130,131,132]))) + tester.assertTrue(mm[-1].getNodalConnectivityIndex().isEqual(DataArrayInt([0,4,8,12,16,20,25,30,35,40,45,50,55,60,65]))) + tester.assertTrue(mm.getFamilyFieldAtLevel(0).isEqual(DataArrayInt([1,2,40,41,42,3,4,5,43]))) + tester.assertTrue(mm.getNumberFieldAtLevel(0).isEqual(DataArrayInt([11,12,50,51,52,13,14,15,53]))) + tester.assertTrue(mm.getFamilyFieldAtLevel(-1).isEqual(DataArrayInt([6,7,200,201,202,8,9,10,11,203,204,205,206,207]))) + tester.assertTrue(mm.getNumberFieldAtLevel(-1).isEqual(DataArrayInt([16,17,300,301,302,18,19,20,21,303,304,305,306,307]))) + refFamIds=[("FAMILLE_ZERO",0),('F1',1),('F2',2),('F3',3),('F4',4),('F5',5),('F6',6),('F7',7),('F8',8),('F9',9),('G1',31),('G10',40),('G11',41),('G2',32),('G3',33),('G4',34),('G5',35),('G6',36),('G7',37),('G8',38),('G9',39),("H1",100)] + tester.assertEqual(set(mm.getFamiliesNames()),set([elt[0] for elt in refFamIds])) + tester.assertEqual(set([mm.getFamilyId(elt) for elt in mm.getFamiliesNames()]),set([elt[1] for elt in refFamIds])) + tester.assertEqual(mm.getGroupsNames(),('myGRP','myGRP1','myGRP2')) + tester.assertEqual(mm.getAllDistributionOfTypes(),[(NORM_TRI3,5),(NORM_QUAD4,9),(NORM_TETRA4,5),(NORM_PENTA6,4),(NORM_ERROR,9)]) + pass + CheckMesh(self,mm) + ## + fieldName="zeField" + t1=(2.3,3,5) + t2=(5.6,7,12) + infoc=["dd [W]","eee [kA]"] + ## + fmts1=MEDFileFieldMultiTS() + f1ts1=MEDFileField1TS() + f1_1=MEDCouplingFieldDouble(ON_CELLS) ; f1_1.setMesh(mm1[0]) ; f1_1.setName(fieldName) + arr1=DataArrayDouble([(10,110),(11,111),(12,112),(13,113),(14,114)]) + arr1.setInfoOnComponents(infoc) + f1_1.setArray(arr1) ; f1_1.setTime(*t1) ; f1_1.setTimeUnit("ms") + f1_1.checkConsistencyLight() + f1ts1.setFieldNoProfileSBT(f1_1) + # + f1_2=MEDCouplingFieldDouble(ON_CELLS) ; f1_2.setMesh(mm1[-1]) ; f1_2.setName(fieldName) + arr2=DataArrayDouble([(15,115),(16,116),(17,117),(18,118),(19,119),(20,120)]) + arr2.setInfoOnComponents(infoc) + f1_2.setArray(arr2) ; f1_2.setTime(*t1) ; f1_2.setTimeUnit("ms") + f1_2.checkConsistencyLight() + f1ts1.setFieldNoProfileSBT(f1_2) + f1_3=MEDCouplingFieldDouble(ON_NODES) ; f1_3.setMesh(mm1[0]) ; f1_3.setName(fieldName) + arr3=DataArrayDouble([(21,121),(22,122),(23,123),(24,124)]) + arr3.setInfoOnComponents(infoc) + f1_3.setArray(arr3) ; f1_3.setTime(*t1) ; f1_3.setTimeUnit("ms") + f1_3.checkConsistencyLight() + f1ts1.setFieldNoProfileSBT(f1_3) + fmts1.pushBackTimeStep(f1ts1) + # + f1ts2=f1ts1.deepCopy() + f1ts2.setTime(t2[1],t2[2],t2[0]) + f1ts2.getUndergroundDataArray()[:]+=2000 + fmts1.pushBackTimeStep(f1ts2) + ### fmts2 + fmts2=MEDFileFieldMultiTS() + f1ts3=MEDFileField1TS() + f2_1=MEDCouplingFieldDouble(ON_CELLS) ; f2_1.setMesh(mm2[0]) ; f2_1.setName(fieldName) + arr4=DataArrayDouble([(50,150),(51,151),(52,152),(53,153)]) + arr4.setInfoOnComponents(infoc) + f2_1.setArray(arr4) ; f2_1.setTime(*t1) ; f2_1.setTimeUnit("ms") + f2_1.checkConsistencyLight() + f1ts3.setFieldNoProfileSBT(f2_1) + f2_2=MEDCouplingFieldDouble(ON_CELLS) ; f2_2.setMesh(mm2[-1]) ; f2_2.setName(fieldName) + arr5=DataArrayDouble([(54,154),(55,155),(56,156),(57,157),(158,158),(59,159),(60,160),(61,161)]) + arr5.setInfoOnComponents(infoc) + f2_2.setArray(arr5) ; f2_2.setTime(*t1) ; f2_2.setTimeUnit("ms") + f2_2.checkConsistencyLight() + f1ts3.setFieldNoProfileSBT(f2_2) + f2_3=MEDCouplingFieldDouble(ON_NODES) ; f2_3.setMesh(mm2[0]) ; f2_3.setName(fieldName) + arr6=DataArrayDouble([(62,162),(63,163),(64,164),(65,165),(66,166)]) + arr6.setInfoOnComponents(infoc) + f2_3.setArray(arr6) ; f2_3.setTime(*t1) ; f2_3.setTimeUnit("ms") + f2_3.checkConsistencyLight() + f1ts3.setFieldNoProfileSBT(f2_3) + fmts2.pushBackTimeStep(f1ts3) + # + f1ts4=f1ts3.deepCopy() + f1ts4.setTime(t2[1],t2[2],t2[0]) + f1ts4.getUndergroundDataArray()[:]+=2000 + fmts2.pushBackTimeStep(f1ts4) + # + mfd1=MEDFileData() + mfd1.setMeshes(MEDFileMeshes()) + mfd1.getMeshes().pushMesh(mm1) + mfd1.setFields(MEDFileFields()) + mfd1.getFields().pushField(fmts1) + # + mfd2=MEDFileData() + mfd2.setMeshes(MEDFileMeshes()) + mfd2.getMeshes().pushMesh(mm2) + mfd2.setFields(MEDFileFields()) + mfd2.getFields().pushField(fmts2) + # ze Call ! + mfd=MEDFileData.Aggregate([mfd1,mfd2]) + def CheckMFD(tester,mfd): + tester.assertEqual(len(mfd.getMeshes()),1) + tester.assertEqual(len(mfd.getFields()),1) + CheckMesh(self,mfd.getMeshes()[0]) + tester.assertEqual(len(mfd.getFields()[0]),2) + zeF1=mfd.getFields()[0][0] + zeF1_1=zeF1.getFieldOnMeshAtLevel(ON_CELLS,0,mfd.getMeshes()[0]) + ref=MEDCouplingFieldDouble.MergeFields([f1_1,f2_1]) + o2n=ref.getMesh().deepCopy().sortCellsInMEDFileFrmt() + ref.renumberCells(o2n) + tester.assertTrue(ref.isEqual(zeF1_1,1e-12,1e-12)) + zeF1_2=zeF1.getFieldOnMeshAtLevel(ON_CELLS,-1,mfd.getMeshes()[0]) + ref=MEDCouplingFieldDouble.MergeFields([f1_2,f2_2]) + o2n=ref.getMesh().deepCopy().sortCellsInMEDFileFrmt() + ref.renumberCells(o2n) + tester.assertTrue(ref.isEqual(zeF1_2,1e-12,1e-12)) + zeF1_3=zeF1.getFieldOnMeshAtLevel(ON_NODES,0,mfd.getMeshes()[0]) + ref=MEDCouplingFieldDouble.MergeFields([f1_3,f2_3]) + o2n=ref.getMesh().deepCopy().sortCellsInMEDFileFrmt() + ref.renumberCells(o2n) + tester.assertTrue(ref.isEqual(zeF1_3,1e-12,1e-12)) + # + zeF2=mfd.getFields()[0][1] + zeF2_1=zeF2.getFieldOnMeshAtLevel(ON_CELLS,0,mfd.getMeshes()[0]) + ref=MEDCouplingFieldDouble.MergeFields([f1_1,f2_1]) + o2n=ref.getMesh().deepCopy().sortCellsInMEDFileFrmt() + ref.renumberCells(o2n) + ref.setTime(*t2) ; ref.getArray()[:]+=2000 + tester.assertTrue(ref.isEqual(zeF2_1,1e-12,1e-12)) + zeF2_2=zeF2.getFieldOnMeshAtLevel(ON_CELLS,-1,mfd.getMeshes()[0]) + ref=MEDCouplingFieldDouble.MergeFields([f1_2,f2_2]) + o2n=ref.getMesh().deepCopy().sortCellsInMEDFileFrmt() + ref.renumberCells(o2n) + ref.setTime(*t2) ; ref.getArray()[:]+=2000 + tester.assertTrue(ref.isEqual(zeF2_2,1e-12,1e-12)) + zeF2_3=zeF2.getFieldOnMeshAtLevel(ON_NODES,0,mfd.getMeshes()[0]) + ref=MEDCouplingFieldDouble.MergeFields([f1_3,f2_3]) + o2n=ref.getMesh().deepCopy().sortCellsInMEDFileFrmt() + ref.renumberCells(o2n) + ref.setTime(*t2) ; ref.getArray()[:]+=2000 + tester.assertTrue(ref.isEqual(zeF2_3,1e-12,1e-12)) + CheckMFD(self,mfd) + mfd1.write(fname1,2) ; mfd2.write(fname2,2) + mfd=MEDFileData.Aggregate([MEDFileData(fname1),MEDFileData(fname2)]) + CheckMFD(self,mfd) + pass pass diff --git a/src/MEDLoader/Swig/MEDLoaderTypemaps.i b/src/MEDLoader/Swig/MEDLoaderTypemaps.i index 6c0fb4ae7..9df1cb801 100644 --- a/src/MEDLoader/Swig/MEDLoaderTypemaps.i +++ b/src/MEDLoader/Swig/MEDLoaderTypemaps.i @@ -197,6 +197,20 @@ static PyObject *convertFieldDoubleVecToPy(const std::vector >& vec) +{ + PyObject *ret(PyList_New(vec.size())); + int rk=0; + for(std::vector< std::pair >::const_iterator iter=vec.begin();iter!=vec.end();iter++,rk++) + { + PyObject *elt=PyTuple_New(2); + PyTuple_SetItem(elt,0,SWIG_From_int((*iter).first)); + PyTuple_SetItem(elt,1,SWIG_From_int((*iter).second)); + PyList_SetItem(ret,rk,elt); + } + return ret; +} + PyObject *convertVecPairVecStToPy(const std::vector< std::pair, std::string > >& vec) { int sz=(int)vec.size();