From 46ddce3caf8076fa140934bb3ce88c5d8040efb9 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Tue, 2 Jan 2018 16:32:59 +0100 Subject: [PATCH] Implemenatation of MEDFileFields::linearToQuadratic + Remapper is now dealing with 1D->0D (3D space) P1P1 PointLocator --- .../Bases/InterpKernelAssert.hxx | 34 +++ src/INTERP_KERNEL/BoundingBox.cxx | 10 + src/INTERP_KERNEL/BoundingBox.hxx | 2 + src/INTERP_KERNEL/CMakeLists.txt | 1 + src/INTERP_KERNEL/Interpolation1D0D.cxx | 57 +++++ src/INTERP_KERNEL/Interpolation1D0D.hxx | 45 ++++ src/INTERP_KERNEL/Interpolation1D0D.txx | 123 +++++++++++ src/INTERP_KERNEL/InterpolationUtils.hxx | 34 ++- src/INTERP_KERNEL/MeshElement.txx | 2 - src/MEDCoupling/MEDCouplingMemArray.cxx | 5 + .../MEDCouplingNormalizedUnstructuredMesh.txx | 2 + src/MEDCoupling/MEDCouplingRefCountObject.cxx | 9 + src/MEDCoupling/MEDCouplingRefCountObject.hxx | 3 +- src/MEDCoupling/MEDCouplingRemapper.cxx | 10 + .../Test/MEDCouplingBasicsTestInterp.cxx | 2 +- .../MEDCouplingRefCountObject.i | 3 +- .../MEDCouplingRemapperTest.py | 45 ++++ src/MEDLoader/MEDFileField.cxx | 206 ++++++++++++++++++ src/MEDLoader/MEDFileField.hxx | 1 + src/MEDLoader/MEDFileField.txx | 22 ++ src/MEDLoader/MEDFileField1TS.hxx | 3 + src/MEDLoader/MEDFileFieldGlobs.cxx | 6 +- src/MEDLoader/MEDFileFieldInternal.cxx | 19 ++ src/MEDLoader/MEDFileFieldInternal.hxx | 2 + src/MEDLoader/MEDFileFieldVisitor.hxx | 1 + src/MEDLoader/Swig/MEDLoaderCommon.i | 10 + src/MEDLoader/Swig/MEDLoaderTest3.py | 42 ++++ .../TestPyWrapGathered_medcoupling.py | 21 ++ src/PyWrapping/medcoupling_pycode | 11 +- 29 files changed, 716 insertions(+), 15 deletions(-) create mode 100644 src/INTERP_KERNEL/Bases/InterpKernelAssert.hxx create mode 100644 src/INTERP_KERNEL/Interpolation1D0D.cxx create mode 100644 src/INTERP_KERNEL/Interpolation1D0D.hxx create mode 100644 src/INTERP_KERNEL/Interpolation1D0D.txx diff --git a/src/INTERP_KERNEL/Bases/InterpKernelAssert.hxx b/src/INTERP_KERNEL/Bases/InterpKernelAssert.hxx new file mode 100644 index 000000000..51daf4284 --- /dev/null +++ b/src/INTERP_KERNEL/Bases/InterpKernelAssert.hxx @@ -0,0 +1,34 @@ +// Copyright (C) 2018 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 __INTERPKERNELASSERT_HXX__ +#define __INTERPKERNELASSERT_HXX__ + +#include "InterpKernelException.hxx" + +#include + +#define IKAssert(a) { bool verdict(a); \ + if(!verdict) { std::ostringstream osszz; osszz << "Assertion \"" << #a << "\" failed into " << __FILE__ << " at line " << __LINE__ << " !"; throw INTERP_KERNEL::Exception(osszz.str()); } } + +#define IKAssertMsg(a,b) { bool verdict(a); \ + if(!verdict) { std::ostringstream osszz; osszz << "Assertion \"" << #a << "\" failed into " << __FILE__ << " at line " << __LINE__ << " with message \"" << b << "\" !"; throw INTERP_KERNEL::Exception(osszz.str()); } } + +#endif diff --git a/src/INTERP_KERNEL/BoundingBox.cxx b/src/INTERP_KERNEL/BoundingBox.cxx index f8c640b25..b20ca4100 100644 --- a/src/INTERP_KERNEL/BoundingBox.cxx +++ b/src/INTERP_KERNEL/BoundingBox.cxx @@ -162,4 +162,14 @@ namespace INTERP_KERNEL return valid; } + void BoundingBox::toCompactData(double data[6]) const + { + data[0]=_coords[XMIN]; + data[1]=_coords[XMAX]; + data[2]=_coords[YMIN]; + data[3]=_coords[YMAX]; + data[4]=_coords[ZMIN]; + data[5]=_coords[ZMAX]; + } + } diff --git a/src/INTERP_KERNEL/BoundingBox.hxx b/src/INTERP_KERNEL/BoundingBox.hxx index f64da9a2c..450e4bcf9 100644 --- a/src/INTERP_KERNEL/BoundingBox.hxx +++ b/src/INTERP_KERNEL/BoundingBox.hxx @@ -53,6 +53,8 @@ namespace INTERP_KERNEL inline void dumpCoords() const; + void toCompactData(double data[6]) const; + private: bool isValid() const; diff --git a/src/INTERP_KERNEL/CMakeLists.txt b/src/INTERP_KERNEL/CMakeLists.txt index 7fe1bd116..159c4edef 100644 --- a/src/INTERP_KERNEL/CMakeLists.txt +++ b/src/INTERP_KERNEL/CMakeLists.txt @@ -37,6 +37,7 @@ SET(interpkernel_SOURCES Interpolation3D.cxx Interpolation2D3D.cxx Interpolation3D1D.cxx + Interpolation1D0D.cxx MeshElement.cxx InterpKernelMeshQuality.cxx InterpKernelCellSimplify.cxx diff --git a/src/INTERP_KERNEL/Interpolation1D0D.cxx b/src/INTERP_KERNEL/Interpolation1D0D.cxx new file mode 100644 index 000000000..1e7862c37 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation1D0D.cxx @@ -0,0 +1,57 @@ +// Copyright (C) 2018 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) + +#include "Interpolation1D0D.hxx" +#include "Interpolation1D0D.txx" + +namespace INTERP_KERNEL +{ + /** + * \class Interpolation1D0D + * \brief Class used to calculate the interpolation between a 1D mesh and 0D mesh (in 3D space) + * Can be seen as a specialization of Interpolation3D, and allows notably the adjustment of bounind boxes. + * + */ + + Interpolation1D0D::Interpolation1D0D() + {} + + Interpolation1D0D::Interpolation1D0D(const InterpolationOptions& io):Interpolation(io) + {} + + /** + * Inspired from PlanarIntersector::adjustBoundingBoxes + */ + void Interpolation1D0D::adjustBoundingBoxes(std::vector& bbox) + { + const int SPACE_DIM = 3; + const double adj(getPrecision());// here precision is used instead of getBoundingBoxAdjustment and getBoundingBoxAdjustmentAbs because in the context only precision is relevant + + long size = bbox.size()/(2*SPACE_DIM); + for (int i=0; i + +namespace INTERP_KERNEL +{ + class INTERPKERNEL_EXPORT Interpolation1D0D : public Interpolation + { + public: + Interpolation1D0D(); + Interpolation1D0D(const InterpolationOptions& io); + template + int interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method); + private: + void adjustBoundingBoxes(std::vector& bbox); + }; +} + +#endif diff --git a/src/INTERP_KERNEL/Interpolation1D0D.txx b/src/INTERP_KERNEL/Interpolation1D0D.txx new file mode 100644 index 000000000..7c62140fe --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation1D0D.txx @@ -0,0 +1,123 @@ +// Copyright (C) 2018 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 __INTERPOLATION1D0D_TXX__ +#define __INTERPOLATION1D0D_TXX__ + +#include "Interpolation1D0D.hxx" +#include "Interpolation.txx" +#include "MeshElement.txx" +#include "PointLocator3DIntersectorP0P0.txx" +#include "PointLocator3DIntersectorP0P1.txx" +#include "PointLocator3DIntersectorP1P0.txx" +#include "PointLocator3DIntersectorP1P1.txx" +#include "Log.hxx" + +#include "BBTree.txx" + +#include "InterpKernelAssert.hxx" + +namespace INTERP_KERNEL +{ + /** + * Very similar to Interpolation3D::interpolateMeshes, except for the bounding boxes that can be + * adjusted in a similar fashion as in InterpolationPlanar::performAdjustmentOfBB() + **/ + template + int Interpolation1D0D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method) + { +#if __cplusplus >= 201103L + constexpr int SPACEDIM=MyMeshType::MY_SPACEDIM; + using ConnType=typename MyMeshType::MyConnType; + IKAssert(SPACEDIM==3); + + if(InterpolationOptions::getIntersectionType() != PointLocator) + INTERP_KERNEL::Exception("Invalid 1D/0D intersection type specified : must be PointLocator."); + + std::string methC ( InterpolationOptions::filterInterpolationMethod(method) ); + if(methC!="P1P1") + throw Exception("Invalid method chosen must be in \"P1P1\"."); + + const double epsilon(getPrecision()); + // create MeshElement objects corresponding to each element of the two meshes + const unsigned long numSrcElems(srcMesh.getNumberOfElements()), numTargetElems(targetMesh.getNumberOfElements()); + + LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements "); + + std::vector*> srcElems(numSrcElems); + + std::map*, int> indices; + + for(unsigned long i = 0 ; i < numSrcElems ; ++i) + srcElems[i] = new MeshElement(i, srcMesh); + + // create empty maps for all source elements + result.resize(targetMesh.getNumberOfNodes()); + + // create BBTree structure + // - get bounding boxes + std::vector bboxes(2*SPACEDIM*numSrcElems); + int* srcElemIdx = new int[numSrcElems]; + for(unsigned long i = 0; i < numSrcElems ; ++i) + { + // get source bboxes in right order + srcElems[i]->getBoundingBox()->toCompactData(bboxes.data()+6*i); + srcElemIdx[i] = srcElems[i]->getIndex(); + } + + adjustBoundingBoxes(bboxes); + const double *bboxPtr(nullptr); + if(numSrcElems>0) + bboxPtr=&bboxes[0]; + BBTree tree(bboxPtr, srcElemIdx, 0, numSrcElems); + const ConnType *trgConnPtr(targetMesh.getConnectivityPtr()),*trgConnIPtr(targetMesh.getConnectivityIndexPtr()); + const ConnType *srcConnPtr(srcMesh.getConnectivityPtr()),*srcConnIPtr(srcMesh.getConnectivityIndexPtr()); + const double *trgCooPtr(targetMesh.getCoordinatesPtr()),*srcCooPtr(srcMesh.getCoordinatesPtr()); + for(unsigned long i = 0; i < numTargetElems; ++i) + { + IKAssert(trgConnIPtr[i+1]==i+1 && trgConnIPtr[i]==i); + std::vector srcSegCondidates; + const double *trgCellPosition(trgCooPtr+SPACEDIM*trgConnPtr[i]); + typename MatrixType::value_type& resRow(result[trgConnPtr[i]]); + tree.getElementsAroundPoint(trgCellPosition, srcSegCondidates); + for(auto srcSeg: srcSegCondidates) + { + IKAssertMsg(srcConnIPtr[srcSeg+1]==2*(srcSeg+1) && srcConnIPtr[srcSeg]==2*srcSeg,"Only implemented for linear 1D source"); + double bc0(0.),bc1(0.); + ConnType srcNode0(srcConnPtr[2*srcSeg]),srcNode1(srcConnPtr[2*srcSeg+1]); + if(IsPointOn3DSeg(srcCooPtr+SPACEDIM*srcNode0,srcCooPtr+SPACEDIM*srcNode1,trgCellPosition,epsilon,bc0,bc1)) + { + resRow.insert(std::make_pair(srcNode0,bc0)); + resRow.insert(std::make_pair(srcNode1,bc1)); + continue; + } + } + } + delete [] srcElemIdx; + for(unsigned long i = 0 ; i < numSrcElems ; ++i) + delete srcElems[i]; + return srcMesh.getNumberOfNodes(); + } +#else + throw INTERP_KERNEL::Exception("Go buying a C++11 compiler :)"); +#endif +} + +#endif diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index e222e1b07..fad2faea7 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -68,7 +68,7 @@ namespace INTERP_KERNEL /* calcul la surface d'un triangle */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ - inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3) + inline double Surf_Tri(const double *P_1,const double *P_2,const double *P_3) { double A=(P_3[1]-P_1[1])*(P_2[0]-P_1[0])-(P_2[1]-P_1[1])*(P_3[0]-P_1[0]); double Surface = 0.5*fabs(A); @@ -83,9 +83,9 @@ namespace INTERP_KERNEL //fonction qui calcul le determinant des vecteurs: P3P1 et P3P2 //(cf doc CGAL). - inline double mon_determinant(const double* P_1, - const double* P_2, - const double* P_3) + inline double mon_determinant(const double *P_1, + const double *P_2, + const double *P_3) { double mon_det=(P_1[0]-P_3[0])*(P_2[1]-P_3[1])-(P_2[0]-P_3[0])*(P_1[1]-P_3[1]); return mon_det; @@ -98,8 +98,7 @@ namespace INTERP_KERNEL { double X=P_1[0]-P_2[0]; double Y=P_1[1]-P_2[1]; - double norme=sqrt(X*X+Y*Y); - return norme; + return sqrt(X*X+Y*Y); } /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ @@ -536,6 +535,29 @@ namespace INTERP_KERNEL } } + /*! + * Compute barycentric coords of point \a point relative to segment S [segStart,segStop] in 3D space. + * If point \a point is further from S than eps false is returned. + * If point \a point projection on S is outside S false is also returned. + * If point \a point is closer from S than eps and its projection inside S true is returned and \a bc0 and \a bc1 output parameter set. + */ + inline bool IsPointOn3DSeg(const double segStart[3], const double segStop[3], const double point[3], double eps, double& bc0, double& bc1) + { + double AB[3]={segStop[0]-segStart[0],segStop[1]-segStart[1],segStop[2]-segStart[2]},AP[3]={point[0]-segStart[0],point[1]-segStart[1],point[2]-segStart[2]}; + double l_AB(sqrt(AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2])); + double AP_dot_AB((AP[0]*AB[0]+AP[1]*AB[1]+AP[2]*AB[2])/(l_AB*l_AB)); + double projOfPOnAB[3]={segStart[0]+AP_dot_AB*AB[0],segStart[1]+AP_dot_AB*AB[1],segStart[2]+AP_dot_AB*AB[2]}; + double V_dist_P_AB[3]={point[0]-projOfPOnAB[0],point[1]-projOfPOnAB[1],point[2]-projOfPOnAB[2]}; + double dist_P_AB(sqrt(V_dist_P_AB[0]*V_dist_P_AB[0]+V_dist_P_AB[1]*V_dist_P_AB[1]+V_dist_P_AB[2]*V_dist_P_AB[2])); + if(dist_P_AB>=eps) + return false;//to far from segment [segStart,segStop] + if(AP_dot_AB<-eps || AP_dot_AB>1.+eps) + return false; + AP_dot_AB=std::max(AP_dot_AB,0.); AP_dot_AB=std::min(AP_dot_AB,1.); + bc0=1.-AP_dot_AB; bc1=AP_dot_AB; + return true; + } + /*! * Calculate pseudo barycentric coordinates of a point p with respect to the quadrangle vertices. * This method makes the assumption that: diff --git a/src/INTERP_KERNEL/MeshElement.txx b/src/INTERP_KERNEL/MeshElement.txx index e2796abad..552ab2625 100644 --- a/src/INTERP_KERNEL/MeshElement.txx +++ b/src/INTERP_KERNEL/MeshElement.txx @@ -61,8 +61,6 @@ namespace INTERP_KERNEL delete _box; } - - ///////////////////////////////////////////////////////////////////// /// ElementBBoxOrder ///////////// ///////////////////////////////////////////////////////////////////// diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index 1c20ae77a..9f04fa9f6 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -4257,6 +4257,11 @@ DataArrayInt *DataArrayInt::checkAndPreparePermutation() const * inversely. * In case of success (no throw) : \c ids1->renumber(ret)->isEqual(ids2) where \a ret is the return of this method. * + * \b Example: + * - \a ids1 : [3,1,103,4,6,10,-7,205] + * - \a ids2 : [-7,1,205,10,6,3,103,4] + * - \a return is : [5,1,6,7,4,3,0,2] because ids2[5]==ids1[0], ids2[1]==ids1[1], ids2[6]==ids1[2]... + * * \return DataArrayInt * - a new instance of DataArrayInt. The caller is to delete this * array using decrRef() as it is no more needed. * \throw If either ids1 or ids2 is null not allocated or not with one components. diff --git a/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.txx b/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.txx index ce207679e..111b6186c 100644 --- a/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.txx +++ b/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.txx @@ -22,6 +22,7 @@ #define __MEDCOUPLINGNORMALIZEDUNSTRUCTUREDMESH_TXX__ #include "MEDCouplingNormalizedUnstructuredMesh.hxx" +#include "InterpKernelAssert.hxx" #include "MEDCouplingUMesh.hxx" #include "MEDCoupling1GTUMesh.hxx" @@ -124,6 +125,7 @@ MEDCouplingNormalizedUnstructuredMesh::~MEDCouplingNormalizedU template void MEDCouplingNormalizedUnstructuredMesh::prepare() { + IKAssert(_mesh->getSpaceDimension()==SPACEDIM); const MEDCoupling::MEDCouplingUMesh *m1(dynamic_cast(_mesh)); if(m1) { diff --git a/src/MEDCoupling/MEDCouplingRefCountObject.cxx b/src/MEDCoupling/MEDCouplingRefCountObject.cxx index 725c1cf19..b7e08b45a 100644 --- a/src/MEDCoupling/MEDCouplingRefCountObject.cxx +++ b/src/MEDCoupling/MEDCouplingRefCountObject.cxx @@ -75,6 +75,15 @@ const char *MEDCoupling::MEDCouplingByteOrderStr() return BIGENDIAN_STR; } +bool MEDCoupling::IsCXX11Compiled() +{ +#if __cplusplus >= 201103L + return true; +#else + return false; +#endif +} + //= std::size_t BigMemoryObject::getHeapMemorySize() const diff --git a/src/MEDCoupling/MEDCouplingRefCountObject.hxx b/src/MEDCoupling/MEDCouplingRefCountObject.hxx index 7fdf91f53..b0fb54ce1 100644 --- a/src/MEDCoupling/MEDCouplingRefCountObject.hxx +++ b/src/MEDCoupling/MEDCouplingRefCountObject.hxx @@ -64,7 +64,8 @@ namespace MEDCoupling MEDCOUPLING_EXPORT int MEDCouplingSizeOfVoidStar(); MEDCOUPLING_EXPORT bool MEDCouplingByteOrder(); MEDCOUPLING_EXPORT const char *MEDCouplingByteOrderStr(); - + MEDCOUPLING_EXPORT bool IsCXX11Compiled(); + class BigMemoryObject { public: diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index b6407f88b..808a25e03 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -36,6 +36,7 @@ #include "Interpolation2D1D.txx" #include "Interpolation2D3D.txx" #include "Interpolation3D1D.txx" +#include "Interpolation1D0D.txx" #include "InterpolationCU.txx" #include "InterpolationCC.txx" @@ -440,6 +441,15 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() INTERP_KERNEL::Interpolation3D1D interpolation(*this); nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method); } + else if(srcMeshDim==1 && trgMeshDim==0 && srcSpaceDim==3) + { + if(getIntersectionType()!=INTERP_KERNEL::PointLocator) + throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 3D space ! Select PointLocator as intersection type !"); + MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh); + MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh); + INTERP_KERNEL::Interpolation1D0D interpolation(*this); + nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method); + } else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3) { if(getIntersectionType()!=INTERP_KERNEL::PointLocator) diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx index 71d30669a..683a8d195 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx @@ -1465,7 +1465,7 @@ void MEDCouplingBasicsTestInterp::test3DInterpP0P0Empty() sourceMesh->allocateCells(0); sourceMesh->finishInsertingCells(); DataArrayDouble *myCoords=DataArrayDouble::New(); - myCoords->alloc(0,0); + myCoords->alloc(0,2); sourceMesh->setCoords(myCoords); myCoords->decrRef(); MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New(); diff --git a/src/MEDCoupling_Swig/MEDCouplingRefCountObject.i b/src/MEDCoupling_Swig/MEDCouplingRefCountObject.i index fa29e6032..031f0c380 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRefCountObject.i +++ b/src/MEDCoupling_Swig/MEDCouplingRefCountObject.i @@ -62,7 +62,8 @@ namespace MEDCoupling int MEDCouplingSizeOfVoidStar(); bool MEDCouplingByteOrder(); const char *MEDCouplingByteOrderStr(); - + bool IsCXX11Compiled(); + class BigMemoryObject { public: diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py index b00dc71ef..a0a0b01a3 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py @@ -1155,6 +1155,51 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertEqual(1, len(rmp.getCrudeMatrix()[0])) pass + @unittest.skipUnless(MEDCouplingHasNumPyBindings() and MEDCouplingHasSciPyBindings() and IsCXX11Compiled(),"requires numpy AND scipy AND C++11") + def testP1P1PL3DSpaceFrom1DTo0D(self): + from scipy.sparse import csr_matrix + from numpy import array + + def generateTrg(eps): + trgArr=DataArrayDouble([(0.5,0.5,0.5),(0.2,0.2,0.2),(0.9,0.9,0.9),(0.7+eps*sqrt(3),0.7-eps*sqrt(3),0.7)]) + trg=MEDCouplingUMesh("trg",0) ; trg.setCoords(trgArr) + trg.allocateCells() + RenumTrg=[2,3,0,1] + for rt in RenumTrg: + trg.insertNextCell(NORM_POINT1,[rt]) + return trg + + srcArr=DataArrayDouble([(0.,0.,1.),(0.,0.,0.),(1.,1.,1.)]) + src=MEDCouplingUMesh("src",1) ; src.setCoords(srcArr) + src.allocateCells() + src.insertNextCell(NORM_SEG2,[1,2]) + # + trg=generateTrg(1e-7)# trg point 3 of trg cell 1 is NOT closer enough to source edge #1 -> not intercepted + # + rem=MEDCouplingRemapper() + rem.setIntersectionType(PointLocator) + self.assertEqual(rem.prepare(src,trg,"P1P1"),1) + mat=rem.getCrudeCSRMatrix() + row=array([2,2, 0,0, 1,1]) # here no ref to point 3 ! + col=array([1,2, 1,2, 1,2]) + data=array([0.1,0.9, 0.5,0.5, 0.8,0.2]) + mExp=csr_matrix((data,(row,col)),shape=(4,3)) + delta=abs(mExp-mat) + self.assertAlmostEqual(delta.sum(),0.,14) + # + trg=generateTrg(1e-14) # trg point 3 of trg cell 1 is closer enough to source edge #1 -> intercepted + rem=MEDCouplingRemapper() + rem.setIntersectionType(PointLocator) + self.assertEqual(rem.prepare(src,trg,"P1P1"),1) + mat=rem.getCrudeCSRMatrix() + row=array([2,2, 3,3, 0,0, 1,1]) # here ref to target point 3 + col=array([1,2, 1,2, 1,2, 1,2]) + data=array([0.1,0.9, 0.3,0.7, 0.5,0.5, 0.8,0.2]) + mExp2=csr_matrix((data,(row,col)),shape=(4,3)) + delta2=abs(mExp2-mat) + self.assertAlmostEqual(delta2.sum(),0.,14) + pass + def checkMatrix(self,mat1,mat2,nbCols,eps): self.assertEqual(len(mat1),len(mat2)) for i in range(len(mat1)): diff --git a/src/MEDLoader/MEDFileField.cxx b/src/MEDLoader/MEDFileField.cxx index f4c906f8b..c310994f7 100644 --- a/src/MEDLoader/MEDFileField.cxx +++ b/src/MEDLoader/MEDFileField.cxx @@ -607,6 +607,212 @@ void MEDFileFields::accept(MEDFileFieldVisitor& visitor) const } } +class MEDFileFieldLin2QuadVisitor : public MEDFileFieldVisitor +{ +public: + MEDFileFieldLin2QuadVisitor(const MEDFileUMesh *lin, const MEDFileUMesh *quad, const MEDFileFieldGlobsReal *linGlobs, MEDFileFields* outFs):_lin(lin),_quad(quad),_lin_globs(linGlobs),_out_fs(outFs),_gt(INTERP_KERNEL::NORM_ERROR),_1ts_update_requested(false) { } + void newFieldEntry(const MEDFileAnyTypeFieldMultiTSWithoutSDA *field) { if(field->getMeshName()!=_lin->getName()) return; _cur_fmts=MEDFileFieldMultiTS::New(); } + void endFieldEntry(const MEDFileAnyTypeFieldMultiTSWithoutSDA *field) { if(_cur_fmts.isNotNull()) { if(_cur_fmts->getNumberOfTS()>0) _out_fs->pushField(_cur_fmts); } } + // + void newTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts); + void endTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts); + // + void newMeshEntry(const MEDFileFieldPerMesh *fpm); + void endMeshEntry(const MEDFileFieldPerMesh *fpm) { } + // + void newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt); + void endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt) { } + // + void newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd); +private: + void updateData(MEDFileFieldPerMeshPerTypePerDisc *pmtd); +private: + const MEDFileUMesh *_lin; + const MEDFileUMesh *_quad; + const MEDFileFieldGlobsReal *_lin_globs; + MEDFileFields *_out_fs; + MCAuto _cur_fmts; + MCAuto _cur_f1ts; + INTERP_KERNEL::NormalizedCellType _gt; + // Info on 1TS modification + bool _1ts_update_requested; + // Cache of matrix to compute faster the values on newly created points + std::string _pfl; + MCAuto _matrix; + MCAuto _new_pts_ids; +}; + +void MEDFileFieldLin2QuadVisitor::newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd) +{ + if(_cur_f1ts.isNull()) + return; + if(pmptpd->getType()!=ON_NODES) + throw INTERP_KERNEL::Exception("Not managed yet for ON_CELLS ON_GAUSS_NE and ON_GAUSS_PT"); + _1ts_update_requested=true; + MEDFileAnyTypeField1TSWithoutSDA *ct(_cur_f1ts->contentNotNullBase()); + int locId(pmptpd->getFather()->locIdOfLeaf(pmptpd)); + MEDFileFieldPerMeshPerTypePerDisc *pmtdToModify(ct->getLeafGivenMeshAndTypeAndLocId(_lin->getName(),_gt,locId)); + std::string pflName(pmptpd->getProfile()); + if(pflName==_pfl && _matrix.isNotNull()) + { + updateData(pmtdToModify); + return ; + } + _pfl=pflName; _matrix.nullify(); _new_pts_ids.nullify(); + MCAuto pfl; + if(pflName.empty()) + pfl=DataArrayInt::Range(0,pmptpd->getNumberOfVals(),1); + else + pfl=_lin_globs->getProfile(pflName)->deepCopy(); + { + MCAuto mesh3D(_lin->getMeshAtLevel(0)),mesh3DQuadratic(_quad->getMeshAtLevel(0)); + MCAuto cellIds(mesh3D->getCellIdsLyingOnNodes(pfl->begin(),pfl->end(),true)); + MCAuto mesh3DQuadraticRestricted(mesh3DQuadratic->buildPartOfMySelf(cellIds->begin(),cellIds->end(),true)); + MCAuto mesh3DQuadraticRestrictedNodeIds(mesh3DQuadraticRestricted->computeFetchedNodeIds()); + MCAuto orphansNodes; + { + MCAuto tmp1(mesh3D->buildPartOfMySelf(cellIds->begin(),cellIds->end(),true)); + MCAuto tmp2(tmp1->computeFetchedNodeIds()); + orphansNodes=pfl->buildSubstraction(tmp2); + } + mesh3DQuadraticRestrictedNodeIds->checkMonotonic(true); + _new_pts_ids=mesh3DQuadraticRestrictedNodeIds->buildSubstraction(pfl); + MCAuto allSeg3; + { + MCAuto a,b,c,d; + MCAuto seg3Tmp(mesh3DQuadraticRestricted->explodeIntoEdges(a,b,c,d)); + allSeg3=MEDCoupling1SGTUMesh::New(seg3Tmp); + } + if(allSeg3->getCellModelEnum()!=INTERP_KERNEL::NORM_SEG3) + throw INTERP_KERNEL::Exception("MEDFileFieldLin2QuadVisitor::newPerMeshPerTypePerDisc : invalid situation where SEG3 expected !"); + MCAuto midPts,cellSeg3Ids; + { + DataArrayInt *nodeConn(allSeg3->getNodalConnectivity()); + nodeConn->rearrange(3); + { + std::vector v(1,2); + midPts=nodeConn->keepSelectedComponents(v); + } + cellSeg3Ids=DataArrayInt::FindPermutationFromFirstToSecond(midPts,_new_pts_ids); + { + std::vector v(2); v[0]=0; v[1]=1; + MCAuto tmp(nodeConn->keepSelectedComponents(v)); + _matrix=tmp->selectByTupleId(cellSeg3Ids->begin(),cellSeg3Ids->end()); + } + nodeConn->rearrange(1); + } + updateData(pmtdToModify); + } +} + +void MEDFileFieldLin2QuadVisitor::updateData(MEDFileFieldPerMeshPerTypePerDisc *pmtd) +{ + pmtd->incrementNbOfVals(_new_pts_ids->getNumberOfTuples()); +} + +void MEDFileFieldLin2QuadVisitor::newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt) +{ + const MEDFileFieldPerMeshPerType *pmpt2(dynamic_cast(pmpt)); + if(!pmpt2) + throw INTERP_KERNEL::Exception("MEDFileFieldLin2QuadVisitor::newPerMeshPerTypeEntry : not managed for structure elements !"); + if(pmpt2->getNumberOfLoc()!=1) + throw INTERP_KERNEL::Exception("MEDFileFieldLin2QuadVisitor::newPerMeshPerTypeEntry : not managed for multi discr per timestep !"); + _gt=pmpt->getGeoType(); +} + +void MEDFileFieldLin2QuadVisitor::newMeshEntry(const MEDFileFieldPerMesh *fpm) +{ + if(fpm->getMeshName()!=_lin->getName()) + throw INTERP_KERNEL::Exception("MEDFileFieldLin2QuadVisitor::newMeshEntry : mismatch into meshName !"); +} + +void MEDFileFieldLin2QuadVisitor::newTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts) +{ + _1ts_update_requested=false; + if(!ts) + return ; + const MEDFileField1TSWithoutSDA *tsd(dynamic_cast(ts)); + if(!tsd) + return ; + MCAuto contentCpy(ts->deepCopy()); + MCAuto contentCpy2(DynamicCastSafe(contentCpy)); + if(contentCpy2.isNull()) + return; + _cur_f1ts=MEDFileField1TS::New(*contentCpy2,true); + _cur_f1ts->shallowCpyGlobs(*_lin_globs); +} + +void MEDFileFieldLin2QuadVisitor::endTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts) +{ + if(_cur_f1ts.isNull()) + return ; + if(_1ts_update_requested) + { + if(!_pfl.empty()) + { + int locId(_cur_f1ts->getLocalizationId(_pfl)); + DataArrayInt *pfl(_cur_f1ts->getProfile(_pfl)); + MCAuto newPfl; + { + std::vector vs(2); + vs[0]=pfl; vs[1]=_new_pts_ids; + newPfl=DataArrayInt::Aggregate(vs); + } + newPfl->setName(_pfl); + { + std::vector locToKill(1,locId); + _cur_f1ts->killLocalizationIds(locToKill); + } + _cur_f1ts->appendProfile(newPfl); + } + DataArrayDouble *arr(_cur_f1ts->getUndergroundDataArray()); + MCAuto res; + { + std::vector v(1,0),v2(1,1); + MCAuto pts0(_matrix->keepSelectedComponents(v)); + MCAuto pts1(_matrix->keepSelectedComponents(v2)); + MCAuto part0(arr->selectByTupleId(*pts0)); + MCAuto part1(arr->selectByTupleId(*pts1)); + res=DataArrayDouble::Add(part0,part1); + res->applyLin(0.5,0.); + } + res=DataArrayDouble::Aggregate(arr,res); + _cur_f1ts->setArray(res); + } + if(_cur_fmts.isNotNull()) + { _cur_fmts->pushBackTimeStep(_cur_f1ts); } + _1ts_update_requested=false; +} + +/*! + * \a newQuad is expected to be the result of MEDFileUMesh::linearToQuadratic of \a oldLin + */ +MCAuto MEDFileFields::linearToQuadratic(const MEDFileMeshes *oldLin, const MEDFileMeshes *newQuad) const +{ + if(!oldLin || !newQuad) + throw INTERP_KERNEL::Exception("MEDFileFields::linearToQuadratic : input meshes must be non NULL !"); + MCAuto ret(MEDFileFields::New()); + for(int i=0;igetNumberOfMeshes();i++) + { + MEDFileMesh *mm(oldLin->getMeshAtPos(i)); + if(!mm) + continue; + MEDFileUMesh *mmu(dynamic_cast(mm)); + if(!mmu) + continue; + MEDFileMesh *mmq(newQuad->getMeshWithName(mmu->getName())); + MEDFileUMesh *mmqu(dynamic_cast(mmq)); + if(!mmqu) + { + std::ostringstream oss; oss << "MEDFileFields::linearToQuadratic : mismatch of name between input meshes for name \"" << mmu->getName() << "\""; + throw INTERP_KERNEL::Exception(oss.str()); + } + MEDFileFieldLin2QuadVisitor vis(mmu,mmqu,this,ret); + accept(vis); + } + return ret; +} + MEDFileAnyTypeFieldMultiTS *MEDFileFields::getFieldAtPos(int i) const { if(i<0 || i>=(int)_fields.size()) diff --git a/src/MEDLoader/MEDFileField.hxx b/src/MEDLoader/MEDFileField.hxx index 724a53321..f173681e8 100644 --- a/src/MEDLoader/MEDFileField.hxx +++ b/src/MEDLoader/MEDFileField.hxx @@ -123,6 +123,7 @@ namespace MEDCoupling MEDLOADER_EXPORT bool changeMeshNames(const std::vector< std::pair >& modifTab); MEDLOADER_EXPORT bool renumberEntitiesLyingOnMesh(const std::string& meshName, const std::vector& oldCode, const std::vector& newCode, const DataArrayInt *renumO2N); MEDLOADER_EXPORT void accept(MEDFileFieldVisitor& visitor) const; + MEDLOADER_EXPORT MCAuto linearToQuadratic(const MEDFileMeshes *oldLin, const MEDFileMeshes *newQuad) const; public: MEDLOADER_EXPORT MEDFileFields *extractPart(const std::map >& extractDef, MEDFileMesh *mm) const; public: diff --git a/src/MEDLoader/MEDFileField.txx b/src/MEDLoader/MEDFileField.txx index e3ba8eeb8..ab4a0d178 100644 --- a/src/MEDLoader/MEDFileField.txx +++ b/src/MEDLoader/MEDFileField.txx @@ -199,6 +199,16 @@ namespace MEDCoupling } } + template + void MEDFileField1TSTemplateWithoutSDA::copyTimeInfoFrom(const typename Traits::FieldType *mcf) + { + if(!mcf) + throw INTERP_KERNEL::Exception("MEDFileField1TSTemplateWithoutSDA::copyTimeInfoFrom : input field is nullptr !"); + int b(0),c(0); + double a(mcf->getTime(b,c)); + setTime(b,c,a); + } + /////////////////////////////////////////////////////// template @@ -397,6 +407,12 @@ namespace MEDCoupling return ReturnSafelyTypedDataArray(arr); } + template + void MEDFileTemplateField1TS::setArray(DataArray *arr) + { + return contentNotNull()->setArray(arr); + } + template typename Traits::ArrayType *MEDFileTemplateField1TS::getUndergroundDataArray() const { @@ -438,6 +454,12 @@ namespace MEDCoupling return ret.retn(); } + template + void MEDFileTemplateField1TS::copyTimeInfoFrom(const typename Traits::FieldType *mcf) + { + contentNotNull()->copyTimeInfoFrom(mcf); + } + /*! * This is the simplest version to fetch a field for MED structure. One drawback : if \a this is a complex field (multi spatial discretization inside a same field) this method will throw exception and more advance * method should be called (getFieldOnMeshAtLevel for example). diff --git a/src/MEDLoader/MEDFileField1TS.hxx b/src/MEDLoader/MEDFileField1TS.hxx index 961fc7506..54f75f887 100644 --- a/src/MEDLoader/MEDFileField1TS.hxx +++ b/src/MEDLoader/MEDFileField1TS.hxx @@ -160,6 +160,7 @@ namespace MEDCoupling 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); + MEDLOADER_EXPORT void copyTimeInfoFrom(const typename Traits::FieldType *mcf); protected: MCAuto< typename Traits::ArrayType > _arr; }; @@ -349,11 +350,13 @@ namespace MEDCoupling public: MEDLOADER_EXPORT static typename Traits::ArrayType *ReturnSafelyTypedDataArray(MCAuto& arr); MEDLOADER_EXPORT typename Traits::ArrayType *getFieldWithProfile(TypeOfField type, int meshDimRelToMax, const MEDFileMesh *mesh, DataArrayInt *&pfl) const; + MEDLOADER_EXPORT void setArray(DataArray *arr); MEDLOADER_EXPORT typename Traits::ArrayType *getUndergroundDataArray() const; MEDLOADER_EXPORT typename Traits::ArrayType *getUndergroundDataArrayExt(std::vector< std::pair,std::pair > >& entries) const; MEDLOADER_EXPORT static MCAuto::FieldType> SetDataArrayInField(MEDCouplingFieldDouble *f, MCAuto& arr); MEDLOADER_EXPORT static MCAuto ToFieldTemplateWithTime(const typename Traits::FieldType *f); public: + MEDLOADER_EXPORT void copyTimeInfoFrom(const typename Traits::FieldType *mcf); MEDLOADER_EXPORT typename Traits::FieldType *field(const MEDFileMesh *mesh) const; MEDLOADER_EXPORT typename Traits::FieldType *getFieldAtLevel(TypeOfField type, int meshDimRelToMax, int renumPol=0) const; MEDLOADER_EXPORT typename Traits::FieldType *getFieldAtTopLevel(TypeOfField type, int renumPol=0) const; diff --git a/src/MEDLoader/MEDFileFieldGlobs.cxx b/src/MEDLoader/MEDFileFieldGlobs.cxx index 11a46ae5d..990cd8801 100644 --- a/src/MEDLoader/MEDFileFieldGlobs.cxx +++ b/src/MEDLoader/MEDFileFieldGlobs.cxx @@ -432,7 +432,7 @@ MEDFileFieldLoc& MEDFileFieldGlobs::getLocalization(const std::string& locName) } /*! - * The returned value is never null. + * The returned value is never null. Borrowed reference returned. */ DataArrayInt *MEDFileFieldGlobs::getProfile(const std::string& pflName) { @@ -1119,7 +1119,7 @@ MEDFileFieldLoc& MEDFileFieldGlobsReal::getLocalization(const std::string& locNa /*! * Returns a profile array, apt for modification, by its name. * \param [in] pflName - the name of the profile of interest. - * \return DataArrayInt * - a non-const pointer to the profile array having the name \a pflName. + * \return DataArrayInt * - Borrowed reference - a non-const pointer to the profile array having the name \a pflName. * \throw If there is no a profile named \a pflName. */ DataArrayInt *MEDFileFieldGlobsReal::getProfile(const std::string& pflName) @@ -1130,7 +1130,7 @@ DataArrayInt *MEDFileFieldGlobsReal::getProfile(const std::string& pflName) /*! * Returns a profile array, apt for modification, by its id. * \param [in] pflId - the id of the profile of interest. - * \return DataArrayInt * - a non-const pointer to the profile array having the id \a pflId. + * \return DataArrayInt * - Borrowed reference - a non-const pointer to the profile array having the id \a pflId. * \throw If there is no a profile with id \a pflId. */ DataArrayInt *MEDFileFieldGlobsReal::getProfileFromId(int pflId) diff --git a/src/MEDLoader/MEDFileFieldInternal.cxx b/src/MEDLoader/MEDFileFieldInternal.cxx index 6df10f498..2a138059b 100644 --- a/src/MEDLoader/MEDFileFieldInternal.cxx +++ b/src/MEDLoader/MEDFileFieldInternal.cxx @@ -807,6 +807,13 @@ int MEDFileFieldPerMeshPerTypePerDisc::getNumberOfTuples() const return _end-_start; } +void MEDFileFieldPerMeshPerTypePerDisc::incrementNbOfVals(int deltaNbVal) +{ + int nbi((_end-_start)/_nval); + _nval+=deltaNbVal; + _end+=nbi*deltaNbVal; +} + DataArray *MEDFileFieldPerMeshPerTypePerDisc::getOrCreateAndGetArray() { return _father->getOrCreateAndGetArray(); @@ -1605,6 +1612,18 @@ const MEDFileFieldPerMeshPerTypePerDisc *MEDFileFieldPerMeshPerTypeCommon::getLe return static_cast(0); } +int MEDFileFieldPerMeshPerTypeCommon::locIdOfLeaf(const MEDFileFieldPerMeshPerTypePerDisc *leaf) const +{ + int ret(0); + for(std::vector< MCAuto >::const_iterator it=_field_pm_pt_pd.begin();it!=_field_pm_pt_pd.end();it++,ret++) + { + const MEDFileFieldPerMeshPerTypePerDisc *cand(*it); + if(cand==leaf) + return ret; + } + throw INTERP_KERNEL::Exception("MEDFileFieldPerMeshPerTypeCommon::locIdOfLeaf : not found such a leaf in this !"); +} + void MEDFileFieldPerMeshPerTypeCommon::fillValues(int& startEntryId, std::vector< std::pair,std::pair > >& entries) const { int i=0; diff --git a/src/MEDLoader/MEDFileFieldInternal.hxx b/src/MEDLoader/MEDFileFieldInternal.hxx index afa3dc9da..8c2e47903 100644 --- a/src/MEDLoader/MEDFileFieldInternal.hxx +++ b/src/MEDLoader/MEDFileFieldInternal.hxx @@ -162,6 +162,7 @@ namespace MEDCoupling int getEnd() const { return _end; } void setEnd(int endd) { _end=endd; } int getNumberOfVals() const { return _nval; } + void incrementNbOfVals(int deltaNbVal); DataArray *getOrCreateAndGetArray(); const DataArray *getOrCreateAndGetArray() const; const std::vector& getInfo() const; @@ -252,6 +253,7 @@ namespace MEDCoupling MEDFileFieldPerMeshPerTypePerDisc *getLeafGivenLocId(int locId); const MEDFileFieldPerMeshPerTypePerDisc *getLeafGivenLocId(int locId) const; int getNumberOfLoc() const { return _field_pm_pt_pd.size(); } + int locIdOfLeaf(const MEDFileFieldPerMeshPerTypePerDisc *leaf) 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); diff --git a/src/MEDLoader/MEDFileFieldVisitor.hxx b/src/MEDLoader/MEDFileFieldVisitor.hxx index 9ad13e9f6..47269e3a9 100644 --- a/src/MEDLoader/MEDFileFieldVisitor.hxx +++ b/src/MEDLoader/MEDFileFieldVisitor.hxx @@ -48,6 +48,7 @@ namespace MEDCoupling virtual void endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt) = 0; // virtual void newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd) = 0; + virtual ~MEDFileFieldVisitor() { } }; } diff --git a/src/MEDLoader/Swig/MEDLoaderCommon.i b/src/MEDLoader/Swig/MEDLoaderCommon.i index e68c6c625..014eade6e 100644 --- a/src/MEDLoader/Swig/MEDLoaderCommon.i +++ b/src/MEDLoader/Swig/MEDLoaderCommon.i @@ -160,6 +160,7 @@ using namespace MEDCoupling; %newobject MEDCoupling::MEDFileFields::partOfThisOnStructureElements; %newobject MEDCoupling::MEDFileFields::__iter__; %newobject MEDCoupling::MEDFileFields::extractPart; +%newobject MEDCoupling::MEDFileFields::linearToQuadratic; %newobject MEDCoupling::MEDFileWritableStandAlone::serialize; %newobject MEDCoupling::MEDFileAnyTypeFieldMultiTS::New; @@ -2221,6 +2222,7 @@ namespace MEDCoupling static MEDFileField1TS *New(DataArrayByte *db) throw(INTERP_KERNEL::Exception); static MEDFileField1TS *New(); MEDCoupling::MEDFileIntField1TS *convertToInt(bool isDeepCpyGlobs=true) const throw(INTERP_KERNEL::Exception); + void copyTimeInfoFrom(MEDCouplingFieldDouble *mcf) throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *field(const MEDFileMesh *mesh) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *getFieldAtLevel(TypeOfField type, int meshDimRelToMax, int renumPol=0) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *getFieldAtTopLevel(TypeOfField type, int renumPol=0) const throw(INTERP_KERNEL::Exception); @@ -2338,6 +2340,7 @@ namespace MEDCoupling // void setFieldNoProfileSBT(const MEDCouplingFieldInt *field) throw(INTERP_KERNEL::Exception); void setFieldProfile(const MEDCouplingFieldInt *field, const MEDFileMesh *mesh, int meshDimRelToMax, const DataArrayInt *profile) throw(INTERP_KERNEL::Exception); + void copyTimeInfoFrom(MEDCouplingFieldInt *mcf) throw(INTERP_KERNEL::Exception); MEDCouplingFieldInt *field(const MEDFileMesh *mesh) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldInt *getFieldAtLevel(TypeOfField type, int meshDimRelToMax, int renumPol=0) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldInt *getFieldAtTopLevel(TypeOfField type, int renumPol=0) const throw(INTERP_KERNEL::Exception); @@ -2408,6 +2411,7 @@ namespace MEDCoupling // void setFieldNoProfileSBT(const MEDCouplingFieldFloat *field) throw(INTERP_KERNEL::Exception); void setFieldProfile(const MEDCouplingFieldFloat *field, const MEDFileMesh *mesh, int meshDimRelToMax, const DataArrayInt *profile) throw(INTERP_KERNEL::Exception); + void copyTimeInfoFrom(MEDCouplingFieldFloat *mcf) throw(INTERP_KERNEL::Exception); MEDCouplingFieldFloat *field(const MEDFileMesh *mesh) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldFloat *getFieldAtLevel(TypeOfField type, int meshDimRelToMax, int renumPol=0) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldFloat *getFieldAtTopLevel(TypeOfField type, int renumPol=0) const throw(INTERP_KERNEL::Exception); @@ -3408,6 +3412,12 @@ namespace MEDCoupling convertToMapIntDataArrayInt(extractDef,extractDefCpp); return self->extractPart(extractDefCpp,mm); } + + MEDFileFields *linearToQuadratic(const MEDFileMeshes *oldLin, const MEDFileMeshes *newQuad) const throw(INTERP_KERNEL::Exception) + { + MCAuto ret(self->linearToQuadratic(oldLin,newQuad)); + return ret.retn(); + } } }; diff --git a/src/MEDLoader/Swig/MEDLoaderTest3.py b/src/MEDLoader/Swig/MEDLoaderTest3.py index e61c9028a..f970090ea 100644 --- a/src/MEDLoader/Swig/MEDLoaderTest3.py +++ b/src/MEDLoader/Swig/MEDLoaderTest3.py @@ -576,6 +576,10 @@ class MEDLoaderTest3(unittest.TestCase): ff1.setTime(3,4,2.3) itt,orr,ti=ff1.getTime() self.assertEqual(3,itt); self.assertEqual(4,orr); self.assertAlmostEqual(2.3,ti,14); + f1.setTime(5.5,7,8) + ff1.copyTimeInfoFrom(f1) + itt,orr,ti=ff1.getTime() + self.assertEqual(7,itt); self.assertEqual(8,orr); self.assertAlmostEqual(5.5,ti,14); da,infos=ff1.getUndergroundDataArrayExt() f2.getArray().setName(da.getName())#da has the same name than f2 self.assertTrue(da.isEqual(f2.getArray(),1e-12)) @@ -6191,6 +6195,44 @@ class MEDLoaderTest3(unittest.TestCase): if os.path.exists(errfname): os.remove(errfname) pass + + def testFieldsLinearToQuadratic(self): + arr=DataArrayDouble([0,1]) + m=MEDCouplingCMesh(); + m.setCoords(arr,arr,arr) + m=m.buildUnstructured() + m2=m.deepCopy() + m2.translate([2,0,0]) + m3=MEDCouplingUMesh.MergeUMeshes([m,m2]) + m3.setName("mesh") + mm=MEDFileUMesh() + mm[0]=m3 + mmq=mm.linearToQuadratic(0) + mms=MEDFileMeshes() ; mms.pushMesh(mm) + mmsq=MEDFileMeshes() ; mmsq.pushMesh(mmq) + # + f=MEDCouplingFieldDouble(ON_NODES) + f.setName("field") + f.setMesh(m3) + f.setTime(3.,1,2) + arr=DataArrayDouble(m3.getNumberOfNodes()) + arr.iota() + f.setArray(arr) + f1ts=MEDFileField1TS() + f1ts.setFieldNoProfileSBT(f) + fmts=MEDFileFieldMultiTS() + fmts.pushBackTimeStep(f1ts) + fs=MEDFileFields() + fs.pushField(fmts) + fs2=fs.linearToQuadratic(mms,mmsq) + # + fToTest=fs2[0][0].field(mmq) + self.assertEqual(fToTest.getTime(),[3.,1,2]) + mTest=MEDCoupling1SGTUMesh(fToTest.getMesh()) + self.assertTrue(mTest.getNodalConnectivity().isEqual(DataArrayInt([1,0,2,3,5,4,6,7,16,17,18,19,20,21,22,23,24,25,26,27,9,8,10,11,13,12,14,15,28,29,30,31,32,33,34,35,36,37,38,39]))) + self.assertTrue(mTest.getCoords().isEqual(DataArrayDouble([0.,0.,0.,1.,0.,0.,0.,1.,0.,1.,1.,0.,0.,0.,1.,1.,0.,1.,0.,1.,1.,1.,1.,1.,2.,0.,0.,3.,0.,0.,2.,1.,0.,3.,1.,0.,2.,0.,1.,3.,0.,1.,2.,1.,1.,3.,1.,1.,0.5, 0.,0.,0.,0.5, 0.,0.5, 1.,0.,1.,0.5, 0.,0.5, 0.,1.,0.,0.5, 1.,0.5, 1.,1.,1.,0.5, 1.,1.,0.,0.5, 0.,0.,0.5, 0.,1.,0.5, 1.,1.,0.5, 2.5, 0.,0.,2.,0.5, 0.,2.5, 1.,0.,3.,0.5, 0.,2.5, 0.,1.,2.,0.5, 1.,2.5, 1.,1.,3.,0.5, 1.,3.,0.,0.5, 2.,0.,0.5, 2.,1.,0.5, 3.,1.,0.5],40,3),1e-12)) + self.assertTrue(fToTest.getArray().isEqual(DataArrayDouble([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,0.5,1,2.5,2,4.5,5,6.5,6,3,2,4,5,8.5,9,10.5,10,12.5,13,14.5,14,11,10,12,13]),1e-12)) + pass pass diff --git a/src/PyWrapping/TestPyWrapGathered_medcoupling.py b/src/PyWrapping/TestPyWrapGathered_medcoupling.py index dd39f8dae..4c91bab5b 100644 --- a/src/PyWrapping/TestPyWrapGathered_medcoupling.py +++ b/src/PyWrapping/TestPyWrapGathered_medcoupling.py @@ -98,6 +98,27 @@ class medcouplingTest(unittest.TestCase): interface=CommInterface() pass + def test5(self): + f=MEDCouplingFieldDouble(ON_NODES) + f.setTime(1.25,3,6) + a,b,c=f.getTime() + self.assertEqual(b,3) ; self.assertEqual(c,6) ; self.assertAlmostEqual(a,1.25,14); + f1ts=MEDFileField1TS() + f1ts.setTime(10,13,10.75) + f.copyTimeInfoFrom(f1ts) + a,b,c=f.getTime() + self.assertEqual(b,10) ; self.assertEqual(c,13) ; self.assertAlmostEqual(a,10.75,14); + f2=MEDCouplingFieldInt(ON_NODES) + f2.copyTimeInfoFrom(f1ts) + a,b,c=f2.getTime() + self.assertEqual(b,10) ; self.assertEqual(c,13) ; self.assertAlmostEqual(a,10.75,14); + f3=MEDCouplingFieldFloat(ON_NODES) + f3.copyTimeInfoFrom(f1ts) + a,b,c=f3.getTime() + self.assertEqual(b,10) ; self.assertEqual(c,13) ; self.assertAlmostEqual(a,10.75,14); + pass + + def partitionerTesterHelper(self,algoSelected): arr=DataArrayDouble(10) ; arr.iota() m=MEDCouplingCMesh() ; m.setCoords(arr,arr) diff --git a/src/PyWrapping/medcoupling_pycode b/src/PyWrapping/medcoupling_pycode index a3dfe1a20..a286b9116 100644 --- a/src/PyWrapping/medcoupling_pycode +++ b/src/PyWrapping/medcoupling_pycode @@ -82,10 +82,19 @@ def MEDCouplingMesh_write(self,fileName): def MEDCouplingField_write(self,fileName): MEDCouplingWriterHelper(self,fileName,WriteField) + +def MEDCouplingFieldT_copyTimeInfoFrom(self,mlf1ts): + assert(isinstance(mlf1ts,MEDFileAnyTypeField1TS)) + a,b,c=mlf1ts.getTime() + self.setTime(c,a,b) + pass MEDCouplingMesh.write=MEDCouplingMesh_write del MEDCouplingMesh_write MEDCouplingField.write=MEDCouplingField_write del MEDCouplingField_write - +MEDCouplingFieldDouble.copyTimeInfoFrom=MEDCouplingFieldT_copyTimeInfoFrom +MEDCouplingFieldInt.copyTimeInfoFrom=MEDCouplingFieldT_copyTimeInfoFrom +MEDCouplingFieldFloat.copyTimeInfoFrom=MEDCouplingFieldT_copyTimeInfoFrom +del MEDCouplingFieldT_copyTimeInfoFrom %} -- 2.30.2