From fd10fe14fcf2d8dee8f9b4a647aa2f36593526ec Mon Sep 17 00:00:00 2001 From: ageay Date: Fri, 3 Apr 2009 15:17:03 +0000 Subject: [PATCH] *** empty log message *** --- .../Bases/NormalizedUnstructuredMesh.hxx | 1 + src/INTERP_KERNEL/CellModel.cxx | 6 + src/INTERP_KERNEL/ConvexIntersector.hxx | 3 +- src/INTERP_KERNEL/ConvexIntersector.txx | 21 + src/INTERP_KERNEL/Geometric2DIntersector.hxx | 2 + src/INTERP_KERNEL/Geometric2DIntersector.txx | 19 + src/INTERP_KERNEL/Interpolation3DSurf.txx | 2 + src/INTERP_KERNEL/InterpolationOptions.hxx | 14 +- src/INTERP_KERNEL/InterpolationPlanar.txx | 29 +- src/INTERP_KERNEL/InterpolationUtils.hxx | 26 ++ src/INTERP_KERNEL/Makefile.am | 6 + src/INTERP_KERNEL/PlanarIntersector.hxx | 3 + src/INTERP_KERNEL/PlanarIntersector.txx | 36 ++ src/INTERP_KERNEL/PlanarIntersectorP1P0.txx | 2 +- src/INTERP_KERNEL/PlanarIntersectorP1P1.hxx | 47 +++ src/INTERP_KERNEL/PlanarIntersectorP1P1.txx | 104 +++++ .../PolyhedronIntersectorP1P0.txx | 2 +- src/INTERP_KERNEL/SplitterTetra.txx | 2 +- .../TriangulationIntersector.hxx | 2 + .../TriangulationIntersector.txx | 29 ++ src/MEDCoupling/MEDCouplingField.cxx | 14 +- src/MEDCoupling/MEDCouplingField.hxx | 14 +- .../MEDCouplingFieldDiscretization.cxx | 118 ++++++ .../MEDCouplingFieldDiscretization.hxx | 73 ++++ src/MEDCoupling/MEDCouplingFieldDouble.cxx | 106 ++--- src/MEDCoupling/MEDCouplingFieldDouble.hxx | 17 +- src/MEDCoupling/MEDCouplingMesh.hxx | 5 + .../MEDCouplingNormalizedUnstructuredMesh.txx | 8 +- src/MEDCoupling/MEDCouplingPointSet.cxx | 142 +++++++ src/MEDCoupling/MEDCouplingPointSet.hxx | 61 +++ ...CouplingSMesh.cxx => MEDCouplingRMesh.cxx} | 41 +- ...CouplingSMesh.hxx => MEDCouplingRMesh.hxx} | 16 +- .../MEDCouplingTimeDiscretization.cxx | 262 ++++++++++++ .../MEDCouplingTimeDiscretization.hxx | 150 +++++++ src/MEDCoupling/MEDCouplingUMesh.cxx | 383 ++++++++++++++---- src/MEDCoupling/MEDCouplingUMesh.hxx | 19 +- src/MEDCoupling/MEDCouplingUMeshDesc.cxx | 230 +++++++++++ src/MEDCoupling/MEDCouplingUMeshDesc.hxx | 64 +++ src/MEDCoupling/Makefile.am | 20 +- src/MEDCoupling/MemArray.cxx | 28 +- src/MEDCoupling/MemArray.hxx | 43 +- src/MEDCoupling/MemArray.txx | 47 ++- src/MEDCoupling/RefCountObject.hxx | 7 + .../Test/MEDCouplingBasicsTest.cxx | 342 +++++++++++++++- .../Test/MEDCouplingBasicsTest.hxx | 13 + src/ParaMEDMEM/BlockTopology.cxx | 4 +- src/ParaMEDMEM/BlockTopology.hxx | 4 +- src/ParaMEDMEM/DEC.cxx | 1 + src/ParaMEDMEM/ElementLocator.cxx | 372 ++--------------- src/ParaMEDMEM/ElementLocator.hxx | 17 +- src/ParaMEDMEM/ICoCoMEDField.cxx | 10 +- src/ParaMEDMEM/InterpolationMatrix.cxx | 230 +---------- src/ParaMEDMEM/InterpolationMatrix.hxx | 11 +- src/ParaMEDMEM/IntersectionDEC.cxx | 4 +- src/ParaMEDMEM/MEDLoader/MEDLoader.cxx | 7 +- src/ParaMEDMEM/ParaFIELD.cxx | 17 +- src/ParaMEDMEM/ParaFIELD.hxx | 2 +- src/ParaMEDMEM/ParaGRID.cxx | 6 +- src/ParaMEDMEM/ParaGRID.hxx | 8 +- src/ParaMEDMEM/ParaMESH.cxx | 4 +- src/ParaMEDMEM/ParaMESH.hxx | 17 +- src/ParaMEDMEM/Test/Makefile.am | 1 + src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx | 3 + .../Test/ParaMEDMEMTest_ICocoTrio.cxx | 251 ++++++++++++ .../Test/ParaMEDMEMTest_IntersectionDEC.cxx | 24 +- ...ParaMEDMEMTest_StructuredCoincidentDEC.cxx | 4 +- src/ParaMEDMEM/Test/test_perf.cxx | 4 +- 67 files changed, 2728 insertions(+), 852 deletions(-) create mode 100644 src/INTERP_KERNEL/PlanarIntersectorP1P1.hxx create mode 100644 src/INTERP_KERNEL/PlanarIntersectorP1P1.txx create mode 100644 src/MEDCoupling/MEDCouplingFieldDiscretization.cxx create mode 100644 src/MEDCoupling/MEDCouplingFieldDiscretization.hxx create mode 100644 src/MEDCoupling/MEDCouplingPointSet.cxx create mode 100644 src/MEDCoupling/MEDCouplingPointSet.hxx rename src/MEDCoupling/{MEDCouplingSMesh.cxx => MEDCouplingRMesh.cxx} (76%) rename src/MEDCoupling/{MEDCouplingSMesh.hxx => MEDCouplingRMesh.hxx} (81%) create mode 100644 src/MEDCoupling/MEDCouplingTimeDiscretization.cxx create mode 100644 src/MEDCoupling/MEDCouplingTimeDiscretization.hxx create mode 100644 src/MEDCoupling/MEDCouplingUMeshDesc.cxx create mode 100644 src/MEDCoupling/MEDCouplingUMeshDesc.hxx create mode 100644 src/ParaMEDMEM/Test/ParaMEDMEMTest_ICocoTrio.cxx diff --git a/src/INTERP_KERNEL/Bases/NormalizedUnstructuredMesh.hxx b/src/INTERP_KERNEL/Bases/NormalizedUnstructuredMesh.hxx index 949e2e655..546c41981 100644 --- a/src/INTERP_KERNEL/Bases/NormalizedUnstructuredMesh.hxx +++ b/src/INTERP_KERNEL/Bases/NormalizedUnstructuredMesh.hxx @@ -30,6 +30,7 @@ namespace INTERP_KERNEL typedef enum { + NORM_POINT0 = 0, NORM_SEG2 = 1, NORM_SEG3 = 2, NORM_TRI3 = 3, diff --git a/src/INTERP_KERNEL/CellModel.cxx b/src/INTERP_KERNEL/CellModel.cxx index db41d916c..c31b7e0e7 100644 --- a/src/INTERP_KERNEL/CellModel.cxx +++ b/src/INTERP_KERNEL/CellModel.cxx @@ -44,6 +44,7 @@ namespace INTERP_KERNEL void CellModel::buildUniqueInstance() { + _map_of_unique_instance.insert(make_pair(NORM_POINT0,CellModel(NORM_POINT0))); _map_of_unique_instance.insert(make_pair(NORM_SEG2,CellModel(NORM_SEG2))); _map_of_unique_instance.insert(make_pair(NORM_SEG3,CellModel(NORM_SEG3))); _map_of_unique_instance.insert(make_pair(NORM_TRI3,CellModel(NORM_TRI3))); @@ -68,6 +69,11 @@ namespace INTERP_KERNEL _dyn=false; switch(type) { + case NORM_POINT0: + { + _nb_of_pts=0; _nb_of_sons=0; _dim=0; + } + break; case NORM_SEG2: { _nb_of_pts=2; _nb_of_sons=0; _dim=1; diff --git a/src/INTERP_KERNEL/ConvexIntersector.hxx b/src/INTERP_KERNEL/ConvexIntersector.hxx index cdbae3d0a..7e1300f96 100644 --- a/src/INTERP_KERNEL/ConvexIntersector.hxx +++ b/src/INTERP_KERNEL/ConvexIntersector.hxx @@ -22,7 +22,7 @@ #include "PlanarIntersectorP0P0.hxx" #include "PlanarIntersectorP0P1.hxx" #include "PlanarIntersectorP1P0.hxx" -#include "InterpolationUtils.hxx" +#include "PlanarIntersectorP1P1.hxx" namespace INTERP_KERNEL { @@ -40,6 +40,7 @@ namespace INTERP_KERNEL bool doRotate, int orientation, int printLevel); double intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS); double intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector& sourceCoords, bool isSourceQuad); + double intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords); private : double _epsilon; }; diff --git a/src/INTERP_KERNEL/ConvexIntersector.txx b/src/INTERP_KERNEL/ConvexIntersector.txx index c15bc3f0e..c61a7ba75 100644 --- a/src/INTERP_KERNEL/ConvexIntersector.txx +++ b/src/INTERP_KERNEL/ConvexIntersector.txx @@ -23,6 +23,7 @@ #include "PlanarIntersectorP0P0.txx" #include "PlanarIntersectorP0P1.txx" #include "PlanarIntersectorP1P0.txx" +#include "PlanarIntersectorP1P1.txx" #include "PolygonAlgorithms.txx" @@ -110,6 +111,26 @@ namespace INTERP_KERNEL return result; } + + template class InterpType> + double ConvexIntersector::intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords) + { + double result = 0; + int nbOfNodesS=sourceCoords.size()/SPACEDIM; + int nbOfNodesT=targetCoords.size()/SPACEDIM; + /*** Compute the intersection area ***/ + INTERP_KERNEL::PolygonAlgorithms P(_epsilon, PlanarIntersector::_precision); + std::deque inter = P.intersectConvexPolygons(&targetCoords[0], &sourceCoords[0], + nbOfNodesT, nbOfNodesS); + double area[SPACEDIM]; + int nb_inter =((int)inter.size())/SPACEDIM; + for(int i = 1; i(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area); + result +=0.5*norm(area); + } + return result; + } } #endif diff --git a/src/INTERP_KERNEL/Geometric2DIntersector.hxx b/src/INTERP_KERNEL/Geometric2DIntersector.hxx index 3bcda434c..6055bad22 100644 --- a/src/INTERP_KERNEL/Geometric2DIntersector.hxx +++ b/src/INTERP_KERNEL/Geometric2DIntersector.hxx @@ -22,6 +22,7 @@ #include "PlanarIntersectorP0P0.hxx" #include "PlanarIntersectorP0P1.hxx" #include "PlanarIntersectorP1P0.hxx" +#include "PlanarIntersectorP1P1.hxx" namespace INTERP_KERNEL { @@ -40,6 +41,7 @@ namespace INTERP_KERNEL double dimCaracteristic, double medianPlane, double precision, int orientation); double intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS); double intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector& sourceCoords, bool isSourceQuad); + double intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords); private: QuadraticPolygon *buildPolygonFrom(const std::vector& coords, NormalizedCellType type); QuadraticPolygon *buildPolygonAFrom(ConnType cell, int nbOfPoints, NormalizedCellType type); diff --git a/src/INTERP_KERNEL/Geometric2DIntersector.txx b/src/INTERP_KERNEL/Geometric2DIntersector.txx index d997fa7be..6763d8a22 100644 --- a/src/INTERP_KERNEL/Geometric2DIntersector.txx +++ b/src/INTERP_KERNEL/Geometric2DIntersector.txx @@ -23,6 +23,7 @@ #include "PlanarIntersectorP0P0.txx" #include "PlanarIntersectorP0P1.txx" #include "PlanarIntersectorP1P0.txx" +#include "PlanarIntersectorP1P1.txx" #include "CellModel.hxx" #include "QuadraticPolygon.hxx" @@ -79,6 +80,24 @@ namespace INTERP_KERNEL return ret; } + template class InterpType> + double Geometric2DIntersector::intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords) + { + int nbOfTargetNodes=targetCoords.size()/SPACEDIM; + std::vector nodes(nbOfTargetNodes); + for(int i=0;i nodes2(nbOfSourceNodes); + for(int i=0;iintersectWith(*p2); + delete p1; delete p2; + return ret; + } + template class InterpType> QuadraticPolygon *Geometric2DIntersector::buildPolygonFrom(const std::vector& coords, NormalizedCellType type) { diff --git a/src/INTERP_KERNEL/Interpolation3DSurf.txx b/src/INTERP_KERNEL/Interpolation3DSurf.txx index cf6eb1d2f..8688ec85b 100644 --- a/src/INTERP_KERNEL/Interpolation3DSurf.txx +++ b/src/INTERP_KERNEL/Interpolation3DSurf.txx @@ -34,6 +34,8 @@ namespace INTERP_KERNEL } Interpolation3DSurf::Interpolation3DSurf(const InterpolationOptions& io):InterpolationPlanar(io) + ,_median_plane(io.getMedianPlane()) + ,_surf_3D_adjustment_eps(io.getBoundingBoxAdjustment()) { } diff --git a/src/INTERP_KERNEL/InterpolationOptions.hxx b/src/INTERP_KERNEL/InterpolationOptions.hxx index d471101e7..34a7d042d 100644 --- a/src/INTERP_KERNEL/InterpolationOptions.hxx +++ b/src/INTERP_KERNEL/InterpolationOptions.hxx @@ -51,28 +51,28 @@ namespace INTERP_KERNEL { double getPrecision() const { return InterpolationOptions::_precision; } void setPrecision(double p) { InterpolationOptions::_precision=p; } - double getMedianPlane() { return InterpolationOptions::_median_plane; } + double getMedianPlane() const { return InterpolationOptions::_median_plane; } void setMedianPlane(double mp) { InterpolationOptions::_median_plane=mp; } - bool getDoRotate() { return InterpolationOptions::_do_rotate; } + bool getDoRotate() const { return InterpolationOptions::_do_rotate; } void setDoRotate( bool dr) { InterpolationOptions::_do_rotate = dr; } - double getBoundingBoxAdjustment() { return InterpolationOptions::_bounding_box_adjustment; } + double getBoundingBoxAdjustment() const { return InterpolationOptions::_bounding_box_adjustment; } void setBoundingBoxAdjustment(double bba) { InterpolationOptions::_bounding_box_adjustment=bba; } - int getOrientation() { return InterpolationOptions::_orientation; } + int getOrientation() const { return InterpolationOptions::_orientation; } void setOrientation(int o) { InterpolationOptions::_orientation=o; } - SplittingPolicy getSplittingPolicy() { return _splitting_policy; } + SplittingPolicy getSplittingPolicy() const { return _splitting_policy; } void setSplittingPolicy(SplittingPolicy sp) { _splitting_policy=sp; } void init() { _print_level=0; _intersection_type=Triangulation; - _precision=1e-12;; + _precision=1e-12; _median_plane=0.5; _do_rotate=true; - _bounding_box_adjustment=0.1; + _bounding_box_adjustment=1e-4; _orientation=0; _splitting_policy=GENERAL_48; } diff --git a/src/INTERP_KERNEL/InterpolationPlanar.txx b/src/INTERP_KERNEL/InterpolationPlanar.txx index decaec360..d48246b49 100644 --- a/src/INTERP_KERNEL/InterpolationPlanar.txx +++ b/src/INTERP_KERNEL/InterpolationPlanar.txx @@ -222,8 +222,35 @@ namespace INTERP_KERNEL break; } } + else if(meth=="P1P1") + { + switch (InterpolationOptions::getIntersectionType()) + { + case Triangulation: + intersector=new TriangulationIntersector(myMeshT,myMeshS,_dim_caracteristic, + InterpolationOptions::getPrecision(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getOrientation(), + InterpolationOptions::getPrintLevel()); + break; + case Convex: + intersector=new ConvexIntersector(myMeshT,myMeshS,_dim_caracteristic, + InterpolationOptions::getPrecision(), + InterpolationOptions::getDoRotate(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getOrientation(), + InterpolationOptions::getPrintLevel()); + break; + case Geometric2D: + intersector=new Geometric2DIntersector(myMeshT, myMeshS, _dim_caracteristic, + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrecision(), + InterpolationOptions::getOrientation()); + break; + } + } else - throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" or \"P1P0\""); + throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\""); /****************************************************************/ /* Create a search tree based on the bounding boxes */ /* Instanciate the intersector and initialise the result vector */ diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index f659d1929..3198b5cc2 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -153,6 +153,32 @@ namespace INTERP_KERNEL std::transform(tmp,tmp+SPACEDIM,quadOut+3*SPACEDIM,std::bind2nd(std::multiplies(),0.5)); } + /*! + * This method builds a potentially non-convex polygon cell built with the first point of 'triIn' the barycenter of two edges starting or ending with + * the first point of 'triIn' and the barycenter of 'triIn'. + * + * @param triIn is a 6 doubles array in full interlace mode, that represents a triangle. + * @param quadOut is a 8 doubles array filled after the following call. + */ + template + inline void fillDualCellOfPolyg(const double *polygIn, int nPtsPolygonIn, double *polygOut) + { + //1st point + std::copy(polygIn,polygIn+SPACEDIM,polygOut); + std::transform(polygIn,polygIn+SPACEDIM,polygIn+SPACEDIM,polygOut+SPACEDIM,std::plus()); + //2nd point + std::transform(polygOut+SPACEDIM,polygOut+2*SPACEDIM,polygOut+SPACEDIM,std::bind2nd(std::multiplies(),0.5)); + double tmp[SPACEDIM]; + // + for(int i=0;i()); + std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+3)*SPACEDIM,std::bind2nd(std::multiplies(),0.5)); + std::transform(polygIn+(i+1)*SPACEDIM,polygIn+(i+2)*SPACEDIM,tmp,tmp,std::plus()); + std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+2)*SPACEDIM,std::bind2nd(std::multiplies(),1./3.)); + } + } + /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ /* calcul les coordonnées du barycentre d'un polygone */ /* le vecteur en entrée est constitué des coordonnées */ diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index 5a5e9f733..45ed3229e 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -67,6 +67,12 @@ PlanarIntersector.hxx \ PlanarIntersector.txx \ PlanarIntersectorP0P0.hxx \ PlanarIntersectorP0P0.txx \ +PlanarIntersectorP0P1.hxx \ +PlanarIntersectorP0P1.txx \ +PlanarIntersectorP1P0.hxx \ +PlanarIntersectorP1P0.txx \ +PlanarIntersectorP1P1.hxx \ +PlanarIntersectorP1P1.txx \ PolygonAlgorithms.hxx \ PolygonAlgorithms.txx \ PolyhedronIntersector.hxx \ diff --git a/src/INTERP_KERNEL/PlanarIntersector.hxx b/src/INTERP_KERNEL/PlanarIntersector.hxx index 9676000bc..3cce25cdc 100644 --- a/src/INTERP_KERNEL/PlanarIntersector.hxx +++ b/src/INTERP_KERNEL/PlanarIntersector.hxx @@ -20,6 +20,7 @@ #define __PLANARINTERSECTOR_HXX__ #include "TargetIntersector.hxx" +#include "NormalizedUnstructuredMesh.hxx" namespace INTERP_KERNEL { @@ -45,6 +46,8 @@ namespace INTERP_KERNEL int projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB); void getRealTargetCoordinates(ConnType icellT, std::vector& coordsT); void getRealSourceCoordinates(ConnType icellS, std::vector& coordsS); + void getRealTargetCoordinatesPermute(ConnType icellT, int offset, std::vector& coordsT); + void getRealSourceCoordinatesPermute(ConnType icellS, int offset, std::vector& coordsS); void getRealCoordinates(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, std::vector& coordsT, std::vector& coordsS, int& orientation); static int projection(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB, double epsilon, double median_plane, bool do_rotate); diff --git a/src/INTERP_KERNEL/PlanarIntersector.txx b/src/INTERP_KERNEL/PlanarIntersector.txx index df6593a01..4493aeb14 100644 --- a/src/INTERP_KERNEL/PlanarIntersector.txx +++ b/src/INTERP_KERNEL/PlanarIntersector.txx @@ -175,6 +175,42 @@ namespace INTERP_KERNEL for(int idim=0; idim::coo2C(_connectS[OTT::conn2C(_connIndexS[OTT::ind2C(icellS)]+iS)])+idim]; } + + /*! + * @param icellT id in target mesh in format of MyMeshType. + * @param offset is a value in C format that indicates the number of circular permutation. + * @param coordsT output val that stores coordinates of the target cell automatically resized to the right length. + */ + template + void PlanarIntersector::getRealTargetCoordinatesPermute(ConnType icellT, int offset, std::vector& coordsT) + { + int nbNodesT=_connIndexT[OTT::ind2C(icellT)+1]-_connIndexT[OTT::ind2C(icellT)]; + coordsT.resize(SPACEDIM*nbNodesT); + for (ConnType iTTmp=0; iTTmp::coo2C(_connectT[OTT::conn2C(_connIndexT[OTT::ind2C(icellT)]+iT)])+idim]; + } + } + + /*! + * @param icellS id in source mesh in format of MyMeshType. + * @param offset is a value in C format that indicates the number of circular permutation. + * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length. + */ + template + void PlanarIntersector::getRealSourceCoordinatesPermute(ConnType icellS, int offset, std::vector& coordsS) + { + int nbNodesS=_connIndexS[OTT::ind2C(icellS)+1]-_connIndexS[OTT::ind2C(icellS)]; + coordsS.resize(SPACEDIM*nbNodesS); + for (ConnType iSTmp=0; iSTmp::coo2C(_connectS[OTT::conn2C(_connIndexS[OTT::ind2C(icellS)]+iS)])+idim]; + } + } /*! * @param icellT id in target mesh in format of MyMeshType. diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P0.txx b/src/INTERP_KERNEL/PlanarIntersectorP1P0.txx index c4750c88b..209ca229c 100644 --- a/src/INTERP_KERNEL/PlanarIntersectorP1P0.txx +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P0.txx @@ -54,7 +54,7 @@ namespace INTERP_KERNEL std::vector targetCellCoords; int orientation=1; PlanarIntersector::getRealTargetCoordinates(OTT::indFC(icellT),targetCellCoords); - NormalizedCellType tT=PlanarIntersector::_meshT.getTypeOfElement(icellT); + NormalizedCellType tT=PlanarIntersector::_meshT.getTypeOfElement(OTT::indFC(icellT)); bool isTargetQuad=CellModel::getCellModel(tT).isQuadratic(); typename MyMatrix::value_type& resRow=res[icellT]; for(typename std::vector::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++) diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P1.hxx b/src/INTERP_KERNEL/PlanarIntersectorP1P1.hxx new file mode 100644 index 000000000..a4012879f --- /dev/null +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P1.hxx @@ -0,0 +1,47 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __PLANARINTERSECTORP1P1_HXX__ +#define __PLANARINTERSECTORP1P1_HXX__ + +#include "PlanarIntersector.hxx" + +namespace INTERP_KERNEL +{ + template + class PlanarIntersectorP1P1 : public PlanarIntersector + { + public: + static const int SPACEDIM=MyMeshType::MY_SPACEDIM; + static const int MESHDIM=MyMeshType::MY_MESHDIM; + typedef typename MyMeshType::MyConnType ConnType; + static const NumberingPolicy numPol=MyMeshType::My_numPol; + protected: + PlanarIntersectorP1P1(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel); + public: + void intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res); + int getNumberOfRowsOfResMatrix() const; + int getNumberOfColsOfResMatrix() const; + + double intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords) { return asLeaf().intersectGeometryGeneral(targetCoords,sourceCoords); } + protected: + ConcreteP1P1Intersector& asLeaf() { return static_cast(*this); } + }; +} + +#endif diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P1.txx b/src/INTERP_KERNEL/PlanarIntersectorP1P1.txx new file mode 100644 index 000000000..e00ec61af --- /dev/null +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P1.txx @@ -0,0 +1,104 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +#ifndef __PLANARINTERSECTORP1P1_TXX__ +#define __PLANARINTERSECTORP1P1_TXX__ + +#include "PlanarIntersectorP1P1.hxx" +#include "InterpolationUtils.hxx" +#include "CellModel.hxx" + +namespace INTERP_KERNEL +{ + template + PlanarIntersectorP1P1::PlanarIntersectorP1P1(const MyMeshType& meshT, const MyMeshType& meshS, + double dimCaracteristic, double precision, double medianPlane, + bool doRotate, int orientation, int printLevel): + PlanarIntersector(meshT,meshS,dimCaracteristic,precision,medianPlane,doRotate,orientation,printLevel) + { + } + + template + int PlanarIntersectorP1P1::getNumberOfRowsOfResMatrix() const + { + return PlanarIntersector::_meshT.getNumberOfNodes(); + } + + template + int PlanarIntersectorP1P1::getNumberOfColsOfResMatrix() const + { + return PlanarIntersector::_meshS.getNumberOfNodes(); + } + + /*! + * This methods split on the fly, into triangles in order to compute dual mesh of target cell (with icellT id in target mesh in C mode). + */ + template + void PlanarIntersectorP1P1::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) + { + int nbNodesT=PlanarIntersector::_connIndexT[icellT+1]-PlanarIntersector::_connIndexT[icellT]; + int orientation=1; + const ConnType *startOfCellNodeConn=PlanarIntersector::_connectT+OTT::conn2C(PlanarIntersector::_connIndexT[icellT]); + std::vector polygT; + PlanarIntersector::getRealTargetCoordinates(OTT::indFC(icellT),polygT); + for(int nodeIdT=0;nodeIdT::coo2C(startOfCellNodeConn[nodeIdT]); + PlanarIntersector::getRealTargetCoordinatesPermute(OTT::indFC(icellT),nodeIdT,polygT); + std::vector polygDualT(SPACEDIM*2*(nbNodesT-1)); + fillDualCellOfPolyg(&polygT[0],polygT.size()/SPACEDIM,&polygDualT[0]); + typename MyMatrix::value_type& resRow=res[curNodeTInCmode]; + for(typename std::vector::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++) + { + int iS=*iter; + int nbNodesS=PlanarIntersector::_connIndexS[iS+1]-PlanarIntersector::_connIndexS[iS]; + const ConnType *startOfCellNodeConnS=PlanarIntersector::_connectS+OTT::conn2C(PlanarIntersector::_connIndexS[iS]); + for(int nodeIdS=0;nodeIdS::coo2C(startOfCellNodeConnS[nodeIdS]); + std::vector polygS; + PlanarIntersector::getRealSourceCoordinatesPermute(OTT::indFC(iS),nodeIdS,polygS); + std::vector polygDualS(SPACEDIM*2*(nbNodesS-1)); + fillDualCellOfPolyg(&polygS[0],polygS.size()/SPACEDIM,&polygDualS[0]); + std::vector polygDualTTmp(polygDualT); + if(SPACEDIM==3) + orientation=PlanarIntersector::projectionThis(&polygDualS[0],&polygDualTTmp[0],polygDualS.size()/SPACEDIM,polygDualT.size()/SPACEDIM); + double surf=orientation*intersectGeometryGeneral(polygDualTTmp,polygDualS); + //filtering out zero surfaces and badly oriented surfaces + // _orientation = -1,0,1 + // -1 : the intersection is taken into account if target and cells have different orientation + // 0 : the intersection is always taken into account + // 1 : the intersection is taken into account if target and cells have the same orientation + if (( surf > 0.0 && PlanarIntersector::_orientation >=0 ) || ( surf < 0.0 && PlanarIntersector::_orientation <=0 )) + { + typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT::indFC(curNodeSInCmode)); + if(iterRes==resRow.end()) + resRow.insert(std::make_pair(OTT::indFC(curNodeSInCmode),surf)); + else + { + double val=(*iterRes).second+surf; + resRow.erase(OTT::indFC(curNodeSInCmode)); + resRow.insert(std::make_pair(OTT::indFC(curNodeSInCmode),val)); + } + } + } + } + } + } +} + +#endif diff --git a/src/INTERP_KERNEL/PolyhedronIntersectorP1P0.txx b/src/INTERP_KERNEL/PolyhedronIntersectorP1P0.txx index 486706873..ab600a106 100644 --- a/src/INTERP_KERNEL/PolyhedronIntersectorP1P0.txx +++ b/src/INTERP_KERNEL/PolyhedronIntersectorP1P0.txx @@ -102,7 +102,7 @@ namespace INTERP_KERNEL else { double val=(*iterRes).second+volume; - resRow.erase(OTT::indFC(*iterCellS)); + resRow.erase(OTT::indFC(sourceNode)); resRow.insert(std::make_pair(OTT::indFC(sourceNode),val)); } } diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx index d470e8b8b..02c20c1cf 100644 --- a/src/INTERP_KERNEL/SplitterTetra.txx +++ b/src/INTERP_KERNEL/SplitterTetra.txx @@ -143,7 +143,7 @@ namespace INTERP_KERNEL const int tab[3][2]={{1,2},{3,2},{1,3}}; const int *curTab=tab[case1]; double pt0[3]; pt0[0]=(tmp[curTab[case2]][0]+tmp[0][0])/2.; pt0[1]=(tmp[curTab[case2]][1]+tmp[0][1])/2.; pt0[2]=(tmp[curTab[case2]][2]+tmp[0][2])/2.; - double pt1[3]; pt1[0]=(tmp[0][0]+tmp[curTab[0]][0]+tmp[curTab[1]][0])/3.; pt1[1]=(tmp[0][1]+tmp[curTab[0]][1]+tmp[curTab[1]][1])/3.; pt1[1]=(tmp[0][2]+tmp[curTab[0]][2]+tmp[curTab[1]][2])/3.; + double pt1[3]; pt1[0]=(tmp[0][0]+tmp[curTab[0]][0]+tmp[curTab[1]][0])/3.; pt1[1]=(tmp[0][1]+tmp[curTab[0]][1]+tmp[curTab[1]][1])/3.; pt1[2]=(tmp[0][2]+tmp[curTab[0]][2]+tmp[curTab[1]][2])/3.; double pt2[3]; pt2[0]=(tmp[0][0]+tmp[1][0]+tmp[2][0]+tmp[3][0])/4.; pt2[1]=(tmp[0][1]+tmp[1][1]+tmp[2][1]+tmp[3][1])/4.; pt2[2]=(tmp[0][2]+tmp[1][2]+tmp[2][2]+tmp[3][2])/4.; std::copy(pt1,pt1+3,output+case2*3); std::copy(pt0,pt0+3,output+(abs(case2-1))*3); diff --git a/src/INTERP_KERNEL/TriangulationIntersector.hxx b/src/INTERP_KERNEL/TriangulationIntersector.hxx index 7ca08511b..e22c430f7 100644 --- a/src/INTERP_KERNEL/TriangulationIntersector.hxx +++ b/src/INTERP_KERNEL/TriangulationIntersector.hxx @@ -22,6 +22,7 @@ #include "PlanarIntersectorP0P0.hxx" #include "PlanarIntersectorP0P1.hxx" #include "PlanarIntersectorP1P0.hxx" +#include "PlanarIntersectorP1P1.hxx" namespace INTERP_KERNEL { @@ -38,6 +39,7 @@ namespace INTERP_KERNEL double dimCaracteristic, double precision, double medianPlane, int orientation, int printLevel); double intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS); double intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector& sourceCoords, bool isSourceQuad); + double intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords); }; } diff --git a/src/INTERP_KERNEL/TriangulationIntersector.txx b/src/INTERP_KERNEL/TriangulationIntersector.txx index e40e4fe9e..dfdf9bc2f 100644 --- a/src/INTERP_KERNEL/TriangulationIntersector.txx +++ b/src/INTERP_KERNEL/TriangulationIntersector.txx @@ -128,6 +128,35 @@ namespace INTERP_KERNEL return result; } + + template class InterpType> + double TriangulationIntersector::intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords) + { + double result = 0.; + ConnType nbNodesS=sourceCoords.size()/SPACEDIM; + ConnType nbNodesT=targetCoords.size()/SPACEDIM; + //Compute the intersection area + double area[SPACEDIM]; + for(ConnType iT = 1; iT inter; + INTERP_KERNEL::intersec_de_triangle(&targetCoords[0],&targetCoords[SPACEDIM*iT],&targetCoords[SPACEDIM*(iT+1)], + &sourceCoords[0],&sourceCoords[SPACEDIM*iS],&sourceCoords[SPACEDIM*(iS+1)], + inter, PlanarIntersector::_dim_caracteristic, + PlanarIntersector::_precision); + ConnType nb_inter=((ConnType)inter.size())/2; + if(nb_inter >3) inter=reconstruct_polygon(inter); + for(ConnType i = 1; i(&inter[0],&inter[2*i],&inter[2*(i+1)],area); + result +=0.5*fabs(area[0]); + } + } + } + return result; + } } #endif diff --git a/src/MEDCoupling/MEDCouplingField.cxx b/src/MEDCoupling/MEDCouplingField.cxx index f36430f33..1a65f4722 100644 --- a/src/MEDCoupling/MEDCouplingField.cxx +++ b/src/MEDCoupling/MEDCouplingField.cxx @@ -18,6 +18,7 @@ // #include "MEDCouplingField.hxx" #include "MEDCouplingMesh.hxx" +#include "MEDCouplingFieldDiscretization.hxx" using namespace ParaMEDMEM; @@ -27,6 +28,11 @@ void MEDCouplingField::updateTime() updateTimeWith(*_mesh); } +TypeOfField MEDCouplingField::getEntity() const +{ + return _type->getEnum(); +} + void MEDCouplingField::setMesh(MEDCouplingMesh *mesh) { if(mesh!=_mesh) @@ -46,11 +52,15 @@ MEDCouplingField::~MEDCouplingField() { if(_mesh) _mesh->decrRef(); + delete _type; +} + +MEDCouplingField::MEDCouplingField(TypeOfField type):_mesh(0),_type(MEDCouplingFieldDiscretization::New(type)) +{ } MEDCouplingField::MEDCouplingField(const MEDCouplingField& other):_name(other._name),_desc(other._name), - _time(other._time),_dt(other._dt),_it(other._it), - _mesh(0),_type(other._type) + _mesh(0),_type(other._type->clone()) { if(other._mesh) { diff --git a/src/MEDCoupling/MEDCouplingField.hxx b/src/MEDCoupling/MEDCouplingField.hxx index d8ab6997c..55c32688f 100644 --- a/src/MEDCoupling/MEDCouplingField.hxx +++ b/src/MEDCoupling/MEDCouplingField.hxx @@ -28,35 +28,29 @@ namespace ParaMEDMEM { class MEDCouplingMesh; + class MEDCouplingFieldDiscretization; class MEDCOUPLING_EXPORT MEDCouplingField : public RefCountObject { public: virtual void checkCoherency() const throw(INTERP_KERNEL::Exception) = 0; void setMesh(ParaMEDMEM::MEDCouplingMesh *mesh); - void setTime(double val) { _time=val; } - double getTime() const { return _time; } - void setDtIt(int dt, int it) { _dt=dt; _it=it; } - void getDtIt(int& dt, int& it) { dt=_dt; it=_it; } ParaMEDMEM::MEDCouplingMesh *getMesh() const { return _mesh; } void setName(const char *name) { _name=name; } void setDescription(const char *desc) { _desc=desc; } const char *getName() const { return _name.c_str(); } - TypeOfField getEntity() const { return _type; } + TypeOfField getEntity() const; protected: void updateTime(); protected: - MEDCouplingField(TypeOfField type):_time(0.),_dt(-1),_it(-1),_mesh(0),_type(type) { } + MEDCouplingField(TypeOfField type); MEDCouplingField(const MEDCouplingField& other); virtual ~MEDCouplingField(); protected: std::string _name; std::string _desc; - double _time; - int _dt; - int _it; MEDCouplingMesh *_mesh; - const TypeOfField _type; + MEDCouplingFieldDiscretization *_type; }; } diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx new file mode 100644 index 000000000..2b45ef279 --- /dev/null +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx @@ -0,0 +1,118 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#include "MEDCouplingFieldDiscretization.hxx" +#include "MEDCouplingMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" + +using namespace ParaMEDMEM; + +const char MEDCouplingFieldDiscretizationP0::REPR[]="P0"; + +const TypeOfField MEDCouplingFieldDiscretizationP0::TYPE=ON_CELLS; + +const char MEDCouplingFieldDiscretizationP1::REPR[]="P1"; + +const TypeOfField MEDCouplingFieldDiscretizationP1::TYPE=ON_NODES; + +MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::New(TypeOfField type) +{ + switch(type) + { + case MEDCouplingFieldDiscretizationP0::TYPE: + return new MEDCouplingFieldDiscretizationP0; + case MEDCouplingFieldDiscretizationP1::TYPE: + return new MEDCouplingFieldDiscretizationP1; + default: + throw INTERP_KERNEL::Exception("Choosen discretization is not implemented yet."); + } +} + +TypeOfField MEDCouplingFieldDiscretizationP0::getEnum() const +{ + return TYPE; +} + +MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationP0::clone() const +{ + return new MEDCouplingFieldDiscretizationP0; +} + +const char *MEDCouplingFieldDiscretizationP0::getStringRepr() const +{ + return REPR; +} + +int MEDCouplingFieldDiscretizationP0::getNumberOfTuples(const MEDCouplingMesh *mesh) const +{ + return mesh->getNumberOfCells(); +} + +void MEDCouplingFieldDiscretizationP0::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) +{ + if(mesh->getNumberOfCells()!=da->getNumberOfTuples()) + { + std::ostringstream message; + message << "Field on cells invalid because there are " << mesh->getNumberOfCells(); + message << " cells in mesh and " << da->getNumberOfTuples() << " tuples in field !"; + throw INTERP_KERNEL::Exception(message.str().c_str()); + } +} + +MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP0::getWeightingField(const MEDCouplingMesh *mesh) const +{ + return mesh->getMeasureField(); +} + +TypeOfField MEDCouplingFieldDiscretizationP1::getEnum() const +{ + return TYPE; +} + +MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationP1::clone() const +{ + return new MEDCouplingFieldDiscretizationP1; +} + +const char *MEDCouplingFieldDiscretizationP1::getStringRepr() const +{ + return REPR; +} + +int MEDCouplingFieldDiscretizationP1::getNumberOfTuples(const MEDCouplingMesh *mesh) const +{ + return mesh->getNumberOfNodes(); +} + +void MEDCouplingFieldDiscretizationP1::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) +{ + if(mesh->getNumberOfNodes()!=da->getNumberOfTuples()) + { + std::ostringstream message; + message << "Field on nodes invalid because there are " << mesh->getNumberOfNodes(); + message << " cells in mesh and " << da->getNumberOfTuples() << " tuples in field !"; + throw INTERP_KERNEL::Exception(message.str().c_str()); + } +} + +MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getWeightingField(const MEDCouplingMesh *mesh) const +{ + //not implemented yet. + //Dual mesh to build + return 0; +} diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx new file mode 100644 index 000000000..3c0c54332 --- /dev/null +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx @@ -0,0 +1,73 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __MEDCOUPLINGFIELDDISCRETIZATION_HXX__ +#define __MEDCOUPLINGFIELDDISCRETIZATION_HXX__ + +#include "MEDCoupling.hxx" +#include "RefCountObject.hxx" +#include "InterpKernelException.hxx" + +namespace ParaMEDMEM +{ + class MEDCouplingMesh; + class DataArrayDouble; + class MEDCouplingFieldDouble; + + class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization + { + public: + static MEDCouplingFieldDiscretization *New(TypeOfField type); + virtual TypeOfField getEnum() const = 0; + virtual MEDCouplingFieldDiscretization *clone() const = 0; + virtual const char *getStringRepr() const = 0; + virtual int getNumberOfTuples(const MEDCouplingMesh *mesh) const = 0; + virtual void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) = 0; + virtual MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const = 0; + }; + + class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationP0 : public MEDCouplingFieldDiscretization + { + public: + TypeOfField getEnum() const; + MEDCouplingFieldDiscretization *clone() const; + const char *getStringRepr() const; + int getNumberOfTuples(const MEDCouplingMesh *mesh) const; + void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception); + MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const; + public: + static const char REPR[]; + static const TypeOfField TYPE; + }; + + class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationP1 : public MEDCouplingFieldDiscretization + { + public: + TypeOfField getEnum() const; + MEDCouplingFieldDiscretization *clone() const; + const char *getStringRepr() const; + int getNumberOfTuples(const MEDCouplingMesh *mesh) const; + void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception); + MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const; + public: + static const char REPR[]; + static const TypeOfField TYPE; + }; +} + +#endif diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx index 77a0ffad6..40cc57247 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx @@ -18,14 +18,16 @@ // #include "MEDCouplingFieldDouble.hxx" #include "MEDCouplingMesh.hxx" +#include "MEDCouplingTimeDiscretization.hxx" +#include "MEDCouplingFieldDiscretization.hxx" #include using namespace ParaMEDMEM; -MEDCouplingFieldDouble *MEDCouplingFieldDouble::New(TypeOfField type) +MEDCouplingFieldDouble *MEDCouplingFieldDouble::New(TypeOfField type, TypeOfTimeDiscretization td) { - return new MEDCouplingFieldDouble(type); + return new MEDCouplingFieldDouble(type,td); } MEDCouplingFieldDouble *MEDCouplingFieldDouble::clone(bool recDeepCpy) const @@ -33,95 +35,95 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::clone(bool recDeepCpy) const return new MEDCouplingFieldDouble(*this,recDeepCpy); } -MEDCouplingFieldDouble::MEDCouplingFieldDouble(TypeOfField type):MEDCouplingField(type),_array(0) +MEDCouplingFieldDouble::MEDCouplingFieldDouble(TypeOfField type, TypeOfTimeDiscretization td):MEDCouplingField(type), + _time_discr(MEDCouplingTimeDiscretization::New(td)) { } -MEDCouplingFieldDouble::MEDCouplingFieldDouble(const MEDCouplingFieldDouble& other, bool deepCpy):MEDCouplingField(other),_array(0) +MEDCouplingFieldDouble::MEDCouplingFieldDouble(const MEDCouplingFieldDouble& other, bool deepCpy):MEDCouplingField(other), + _time_discr(other._time_discr->performCpy(deepCpy)) { - if(other._array) - _array=other._array->performCpy(deepCpy); } MEDCouplingFieldDouble::~MEDCouplingFieldDouble() { - if(_array) - _array->decrRef(); + delete _time_discr; } void MEDCouplingFieldDouble::checkCoherency() const throw(INTERP_KERNEL::Exception) { if(!_mesh) throw INTERP_KERNEL::Exception("Field invalid because no mesh specified !"); - if(!_array) + if(!getArray()) throw INTERP_KERNEL::Exception("Field invalid because no values set !"); - if(_type==ON_CELLS) - { - if(_mesh->getNumberOfCells()!=_array->getNumberOfTuples()) - { - std::ostringstream message; - message << "Field on cells invalid because there are " << _mesh->getNumberOfCells(); - message << " cells in mesh and " << _array->getNumberOfTuples() << " tuples in field !"; - throw INTERP_KERNEL::Exception(message.str().c_str()); - } - } - else if(_type==ON_NODES) - { - if(_mesh->getNumberOfNodes()!=_array->getNumberOfTuples()) - { - std::ostringstream message; - message << "Field on nodes invalid because there are " << _mesh->getNumberOfNodes(); - message << " cells in mesh and " << _array->getNumberOfTuples() << " tuples in field !"; - throw INTERP_KERNEL::Exception(message.str().c_str()); - } - } - else - throw INTERP_KERNEL::Exception("Field of undifined type !!!"); + _type->checkCoherencyBetween(_mesh,getArray()); +} + +double MEDCouplingFieldDouble::accumulate(int compId) const +{ + const double *ptr=getArray()->getConstPointer(); + int nbTuple=getArray()->getNumberOfTuples(); + int nbComps=getArray()->getNumberOfComponents(); + double ret=0.; + for(int i=0;igetWeightingField(_mesh); + const double *ptr=weight->getArray()->getConstPointer(); + int nbOfValues=weight->getArray()->getNbOfElems(); + double ret=0.; + for (int i=0; idecrRef(); + return ret; +} + +void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception) +{ + _time_discr->checkNoTimePresence(); +} + +void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double time, double *res) const throw(INTERP_KERNEL::Exception) +{ + _time_discr->checkTimePresence(time); } void MEDCouplingFieldDouble::applyLin(double a, double b, int compoId) { - double *ptr=_array->getPointer(); + double *ptr=getArray()->getPointer(); ptr+=compoId; - int nbOfComp=_array->getNumberOfComponents(); - int nbOfTuple=_array->getNumberOfTuples(); + int nbOfComp=getArray()->getNumberOfComponents(); + int nbOfTuple=getArray()->getNumberOfTuples(); for(int i=0;igetNumberOfComponents(); + return getArray()->getNumberOfComponents(); } int MEDCouplingFieldDouble::getNumberOfTuples() const throw(INTERP_KERNEL::Exception) { if(!_mesh) throw INTERP_KERNEL::Exception("Impossible to retrieve number of tuples because no mesh specified !"); - if(_type==ON_CELLS) - return _mesh->getNumberOfCells(); - else if(_type==ON_NODES) - return _mesh->getNumberOfNodes(); - else - throw INTERP_KERNEL::Exception("Impossible to retrieve number of tuples because type of entity not implemented yet !"); + return _type->getNumberOfTuples(_mesh); } void MEDCouplingFieldDouble::updateTime() { MEDCouplingField::updateTime(); - if(_array) - updateTimeWith(*_array); + if(getArray()) + updateTimeWith(*getArray()); } void MEDCouplingFieldDouble::setArray(DataArrayDouble *array) { - if(array!=_array) - { - if(_array) - _array->decrRef(); - _array=array; - if(_array) - _array->incrRef(); - declareAsNew(); - } + _time_discr->setArray(array,this); } diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.hxx b/src/MEDCoupling/MEDCouplingFieldDouble.hxx index 57681dffe..8a7d32967 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.hxx @@ -21,6 +21,7 @@ #include "MEDCoupling.hxx" #include "MEDCouplingField.hxx" +#include "MEDCouplingTimeDiscretization.hxx" #include "MemArray.hxx" namespace ParaMEDMEM @@ -28,23 +29,29 @@ namespace ParaMEDMEM class MEDCOUPLING_EXPORT MEDCouplingFieldDouble : public MEDCouplingField { public: - static MEDCouplingFieldDouble *New(TypeOfField type); + static MEDCouplingFieldDouble *New(TypeOfField type, TypeOfTimeDiscretization td=NO_TIME); MEDCouplingFieldDouble *clone(bool recDeepCpy) const; void checkCoherency() const throw(INTERP_KERNEL::Exception); - double getIJ(int tupleId, int compoId) const { return _array->getIJ(tupleId,compoId); } + void setTime(double val, int dt, int it) { _time_discr->setTime(val,dt,it); } + double getTime(int& dt, int& it) const { return _time_discr->getTime(dt,it); } + double getIJ(int tupleId, int compoId) const { return getArray()->getIJ(tupleId,compoId); } void setArray(DataArrayDouble *array); - DataArrayDouble *getArray() const { return _array; } + DataArrayDouble *getArray() const { return _time_discr->getArray(); } + double accumulate(int compId) const; + double measureAccumulate(int compId) const; + void getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception); + void getValueOn(const double *spaceLoc, double time, double *res) const throw(INTERP_KERNEL::Exception); //! \b temporary void applyLin(double a, double b, int compoId); int getNumberOfComponents() const; int getNumberOfTuples() const throw(INTERP_KERNEL::Exception); void updateTime(); private: - MEDCouplingFieldDouble(TypeOfField type); + MEDCouplingFieldDouble(TypeOfField type, TypeOfTimeDiscretization td); MEDCouplingFieldDouble(const MEDCouplingFieldDouble& other, bool deepCpy); ~MEDCouplingFieldDouble(); private: - DataArrayDouble *_array; + MEDCouplingTimeDiscretization *_time_discr; }; } diff --git a/src/MEDCoupling/MEDCouplingMesh.hxx b/src/MEDCoupling/MEDCouplingMesh.hxx index 3102ad515..3b860a226 100644 --- a/src/MEDCoupling/MEDCouplingMesh.hxx +++ b/src/MEDCoupling/MEDCouplingMesh.hxx @@ -25,6 +25,8 @@ namespace ParaMEDMEM { + class MEDCouplingFieldDouble; + class MEDCOUPLING_EXPORT MEDCouplingMesh : public RefCountObject { public: @@ -37,6 +39,9 @@ namespace ParaMEDMEM virtual int getNumberOfNodes() const = 0; virtual int getSpaceDimension() const = 0; virtual int getMeshDimension() const = 0; + // tools + virtual void getBoundingBox(double *bbox) const = 0; + virtual MEDCouplingFieldDouble *getMeasureField() const = 0; protected: MEDCouplingMesh() { } MEDCouplingMesh(const MEDCouplingMesh& other):_name(other._name) { } diff --git a/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.txx b/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.txx index 745137a96..65970cfe9 100644 --- a/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.txx +++ b/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.txx @@ -41,7 +41,7 @@ void MEDCouplingNormalizedUnstructuredMesh::getBoundingBox(dou boundingBox[SPACEDIM+i]=-std::numeric_limits::max(); } ParaMEDMEM::DataArrayDouble *array=_mesh->getCoords(); - const double *ptr=array->getPointer(); + const double *ptr=array->getConstPointer(); int nbOfPts=array->getNbOfElems()/SPACEDIM; for(int j=0;j const double *MEDCouplingNormalizedUnstructuredMesh::getCoordinatesPtr() const { ParaMEDMEM::DataArrayDouble *array=_mesh->getCoords(); - return array->getPointer(); + return array->getConstPointer(); } template @@ -124,8 +124,8 @@ void MEDCouplingNormalizedUnstructuredMesh::prepare() _conn_for_interp=new int[initialConnSize-nbOfCell]; _conn_index_for_interp=new int[nbOfCell+1]; _conn_index_for_interp[0]=0; - const int *work_conn=_mesh->getNodalConnectivity()->getPointer()+1; - const int *work_conn_index=_mesh->getNodalConnectivityIndex()->getPointer(); + const int *work_conn=_mesh->getNodalConnectivity()->getConstPointer()+1; + const int *work_conn_index=_mesh->getNodalConnectivityIndex()->getConstPointer(); int *work_conn_for_interp=_conn_for_interp; int *work_conn_index_for_interp=_conn_index_for_interp; for(int i=0;iperformCpy(deepCpy); +} + +MEDCouplingPointSet::~MEDCouplingPointSet() +{ + if(_coords) + _coords->decrRef(); +} + +int MEDCouplingPointSet::getNumberOfNodes() const +{ + if(_coords) + return _coords->getNumberOfTuples(); + else + throw INTERP_KERNEL::Exception("Unable to get number of nodes because no coordinates specified !"); +} + +int MEDCouplingPointSet::getSpaceDimension() const +{ + if(_coords) + return _coords->getNumberOfComponents(); + else + throw INTERP_KERNEL::Exception("Unable to get space dimension because no coordinates specified !"); +} + +void MEDCouplingPointSet::updateTime() +{ + if(_coords) + { + updateTimeWith(*_coords); + } +} + +bool MEDCouplingPointSet::isStructured() const +{ + return false; +} + +void MEDCouplingPointSet::setCoords(DataArrayDouble *coords) +{ + if( coords != _coords ) + { + if (_coords) + _coords->decrRef(); + _coords=coords; + if(_coords) + _coords->incrRef(); + declareAsNew(); + } +} + +bool MEDCouplingPointSet::areCoordsEqual(const MEDCouplingPointSet& other, double prec) const +{ + return _coords->isEqual(other._coords,prec); +} + +void MEDCouplingPointSet::getBoundingBox(double *bbox) const +{ + int dim=getSpaceDimension(); + for (int idim=0; idim::max(); + bbox[idim*2+1]=-std::numeric_limits::max(); + } + const double *coords=_coords->getConstPointer(); + int nbnodes=getNumberOfNodes(); + for (int i=0; i coords[i*dim+idim] ) + { + bbox[idim*2] = coords[i*dim+idim] ; + } + if ( bbox[idim*2+1] < coords[i*dim+idim] ) + { + bbox[idim*2+1] = coords[i*dim+idim] ; + } + } + } +} + +// ============================================= +// Intersect Bounding Box given 2 Bounding Boxes +// ============================================= +bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps) +{ + double bbtemp[2*dim]; + double deltamax=0.0; + + for (int i=0; i< dim; i++) + { + double delta = bb1[2*i+1]-bb1[2*i]; + if ( delta > deltamax ) + { + deltamax = delta ; + } + } + for (int i=0; i + +namespace ParaMEDMEM +{ + class DataArrayInt; + class DataArrayDouble; + + class MEDCOUPLING_EXPORT MEDCouplingPointSet : public MEDCouplingMesh + { + protected: + MEDCouplingPointSet(); + MEDCouplingPointSet(const MEDCouplingPointSet& other, bool deepCpy); + ~MEDCouplingPointSet(); + public: + void updateTime(); + bool isStructured() const; + int getNumberOfNodes() const; + int getSpaceDimension() const; + void setCoords(DataArrayDouble *coords); + DataArrayDouble *getCoords() const { return _coords; } + bool areCoordsEqual(const MEDCouplingPointSet& other, double prec) const; + void getBoundingBox(double *bbox) const; + virtual MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const = 0; + //! size of returned tinyInfo must be always the same. + virtual void getTinySerializationInformation(std::vector& tinyInfo) const = 0; + virtual void resizeForSerialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2) = 0; + virtual void serialize(DataArrayInt *&a1, DataArrayDouble *&a2) = 0; + virtual MEDCouplingPointSet *buildObjectFromUnserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2) = 0; + virtual void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems) = 0; + protected: + static bool intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps); + protected: + DataArrayDouble *_coords; + }; +} + +#endif diff --git a/src/MEDCoupling/MEDCouplingSMesh.cxx b/src/MEDCoupling/MEDCouplingRMesh.cxx similarity index 76% rename from src/MEDCoupling/MEDCouplingSMesh.cxx rename to src/MEDCoupling/MEDCouplingRMesh.cxx index 4eb480736..70e2fc7da 100644 --- a/src/MEDCoupling/MEDCouplingSMesh.cxx +++ b/src/MEDCoupling/MEDCouplingRMesh.cxx @@ -16,16 +16,16 @@ // // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -#include "MEDCouplingSMesh.hxx" +#include "MEDCouplingRMesh.hxx" #include "MemArray.hxx" using namespace ParaMEDMEM; -MEDCouplingSMesh::MEDCouplingSMesh():_x_array(0),_y_array(0),_z_array(0) +MEDCouplingRMesh::MEDCouplingRMesh():_x_array(0),_y_array(0),_z_array(0) { } -MEDCouplingSMesh::~MEDCouplingSMesh() +MEDCouplingRMesh::~MEDCouplingRMesh() { if(_x_array) _x_array->decrRef(); @@ -35,12 +35,12 @@ MEDCouplingSMesh::~MEDCouplingSMesh() _z_array->decrRef(); } -MEDCouplingSMesh *MEDCouplingSMesh::New() +MEDCouplingRMesh *MEDCouplingRMesh::New() { - return new MEDCouplingSMesh; + return new MEDCouplingRMesh; } -void MEDCouplingSMesh::updateTime() +void MEDCouplingRMesh::updateTime() { if(_x_array) updateTimeWith(*_x_array); @@ -50,15 +50,15 @@ void MEDCouplingSMesh::updateTime() updateTimeWith(*_z_array); } -bool MEDCouplingSMesh::isEqual(const MEDCouplingMesh *other, double prec) const +bool MEDCouplingRMesh::isEqual(const MEDCouplingMesh *other, double prec) const { - const MEDCouplingSMesh *otherC=dynamic_cast(other); + const MEDCouplingRMesh *otherC=dynamic_cast(other); if(!otherC) return false; return true; } -void MEDCouplingSMesh::checkCoherency() const throw(INTERP_KERNEL::Exception) +void MEDCouplingRMesh::checkCoherency() const throw(INTERP_KERNEL::Exception) { const char msg0[]="Invalid "; const char msg1[]=" array ! Must contain more than 1 element."; @@ -82,12 +82,12 @@ void MEDCouplingSMesh::checkCoherency() const throw(INTERP_KERNEL::Exception) } } -bool MEDCouplingSMesh::isStructured() const +bool MEDCouplingRMesh::isStructured() const { return true; } -int MEDCouplingSMesh::getNumberOfCells() const +int MEDCouplingRMesh::getNumberOfCells() const { int ret=1; if(_x_array) @@ -99,7 +99,7 @@ int MEDCouplingSMesh::getNumberOfCells() const return ret; } -int MEDCouplingSMesh::getNumberOfNodes() const +int MEDCouplingRMesh::getNumberOfNodes() const { int ret=1; if(_x_array) @@ -111,7 +111,7 @@ int MEDCouplingSMesh::getNumberOfNodes() const return ret; } -int MEDCouplingSMesh::getSpaceDimension() const +int MEDCouplingRMesh::getSpaceDimension() const { int ret=0; if(_x_array) @@ -123,7 +123,7 @@ int MEDCouplingSMesh::getSpaceDimension() const return ret; } -int MEDCouplingSMesh::getMeshDimension() const +int MEDCouplingRMesh::getMeshDimension() const { int ret=0; if(_x_array) @@ -135,7 +135,7 @@ int MEDCouplingSMesh::getMeshDimension() const return ret; } -DataArrayDouble *MEDCouplingSMesh::getCoordsAt(int i) const throw(INTERP_KERNEL::Exception) +DataArrayDouble *MEDCouplingRMesh::getCoordsAt(int i) const throw(INTERP_KERNEL::Exception) { switch(i) { @@ -150,7 +150,7 @@ DataArrayDouble *MEDCouplingSMesh::getCoordsAt(int i) const throw(INTERP_KERNEL: } } -void MEDCouplingSMesh::setCoords(DataArrayDouble *coordsX, DataArrayDouble *coordsY, DataArrayDouble *coordsZ) +void MEDCouplingRMesh::setCoords(DataArrayDouble *coordsX, DataArrayDouble *coordsY, DataArrayDouble *coordsZ) { if(_x_array) _x_array->decrRef(); @@ -170,3 +170,12 @@ void MEDCouplingSMesh::setCoords(DataArrayDouble *coordsX, DataArrayDouble *coor declareAsNew(); } +void MEDCouplingRMesh::getBoundingBox(double *bbox) const +{ + //not implemented yet ! +} + +MEDCouplingFieldDouble *MEDCouplingRMesh::getMeasureField() const +{ + //not implemented yet ! +} diff --git a/src/MEDCoupling/MEDCouplingSMesh.hxx b/src/MEDCoupling/MEDCouplingRMesh.hxx similarity index 81% rename from src/MEDCoupling/MEDCouplingSMesh.hxx rename to src/MEDCoupling/MEDCouplingRMesh.hxx index fb170beda..e0774596c 100644 --- a/src/MEDCoupling/MEDCouplingSMesh.hxx +++ b/src/MEDCoupling/MEDCouplingRMesh.hxx @@ -16,19 +16,20 @@ // // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -#ifndef __PARAMEDMEM_MEDCOUPLINGSMESH_HXX__ -#define __PARAMEDMEM_MEDCOUPLINGSMESH_HXX__ +#ifndef __PARAMEDMEM_MEDCOUPLINGRMESH_HXX__ +#define __PARAMEDMEM_MEDCOUPLINGRMESH_HXX__ +#include "MEDCoupling.hxx" #include "MEDCouplingMesh.hxx" namespace ParaMEDMEM { class DataArrayDouble; - class MEDCouplingSMesh : public MEDCouplingMesh + class MEDCouplingRMesh : public MEDCouplingMesh { public: - static MEDCouplingSMesh *New(); + static MEDCouplingRMesh *New(); void updateTime(); bool isEqual(const MEDCouplingMesh *other, double prec) const; void checkCoherency() const throw(INTERP_KERNEL::Exception); @@ -41,9 +42,12 @@ namespace ParaMEDMEM void setCoords(DataArrayDouble *coordsX, DataArrayDouble *coordsY=0, DataArrayDouble *coordsZ=0); + // tools + void getBoundingBox(double *bbox) const; + MEDCouplingFieldDouble *getMeasureField() const; private: - MEDCouplingSMesh(); - ~MEDCouplingSMesh(); + MEDCouplingRMesh(); + ~MEDCouplingRMesh(); private: DataArrayDouble *_x_array; DataArrayDouble *_y_array; diff --git a/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx b/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx new file mode 100644 index 000000000..892a52caa --- /dev/null +++ b/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx @@ -0,0 +1,262 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#include "MEDCouplingTimeDiscretization.hxx" +#include "MemArray.hxx" + +#include + +using namespace ParaMEDMEM; + +const char MEDCouplingNoTimeLabel::EXCEPTION_MSG[]="MEDCouplingNoTimeLabel::setTime : no time info attached."; + +const char MEDCouplingWithTimeStep::EXCEPTION_MSG[]="No data on this time."; + +const double MEDCouplingTimeDiscretization::TIME_TOLERANCE_DFT=1.e-12; + +MEDCouplingTimeDiscretization *MEDCouplingTimeDiscretization::New(TypeOfTimeDiscretization type) +{ + switch(type) + { + case MEDCouplingNoTimeLabel::DISCRETIZATION: + return new MEDCouplingNoTimeLabel; + case MEDCouplingWithTimeStep::DISCRETIZATION: + return new MEDCouplingWithTimeStep; + default: + throw INTERP_KERNEL::Exception("Time discretization not implemented yet"); + } +} + +MEDCouplingTimeDiscretization::MEDCouplingTimeDiscretization():_time_tolerance(TIME_TOLERANCE_DFT),_array(0) +{ +} + +MEDCouplingTimeDiscretization::MEDCouplingTimeDiscretization(const MEDCouplingTimeDiscretization& other, bool deepCpy):_time_tolerance(other._time_tolerance) +{ + if(other._array) + _array=other._array->performCpy(deepCpy); + else + _array=0; +} + +MEDCouplingTimeDiscretization::~MEDCouplingTimeDiscretization() +{ + if(_array) + _array->decrRef(); +} + +void MEDCouplingTimeDiscretization::setArray(DataArrayDouble *array, TimeLabel *owner) +{ + if(array!=_array) + { + if(_array) + _array->decrRef(); + _array=array; + if(_array) + _array->incrRef(); + owner->declareAsNew(); + } +} + +void MEDCouplingTimeDiscretization::setArrays(const std::vector& arrays, TimeLabel *owner) throw(INTERP_KERNEL::Exception) +{ + if(arrays.size()!=1) + throw INTERP_KERNEL::Exception("MEDCouplingTimeDiscretization::setArrays : number of arrays must be one."); + setArray(arrays.back(),owner); +} + +void MEDCouplingTimeDiscretization::getArrays(std::vector& arrays) const +{ + arrays.resize(1); + arrays[0]=_array; +} + +bool MEDCouplingTimeDiscretization::isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception) +{ + int dt,it; + double time1=getEndTime(dt,it)-_time_tolerance; + double time2=other->getStartTime(dt,it)+other->getTimeTolerance(); + return time1<=time2; +} + +bool MEDCouplingTimeDiscretization::isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception) +{ + int dt,it; + double time1=getEndTime(dt,it)+_time_tolerance; + double time2=other->getStartTime(dt,it)-other->getTimeTolerance(); + return time1_time_tolerance) + { + std::ostringstream stream; + stream << "The field is defined on time " << _time << " with eps=" << _time_tolerance << " and asking time = " << time << " !"; + throw INTERP_KERNEL::Exception(stream.str().c_str()); + } +} + +DataArrayDouble *MEDCouplingWithTimeStep::getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception) +{ + if(std::fabs(time-_time)<=_time_tolerance) + { + if(_array) + _array->incrRef(); + return _array; + } + else + throw INTERP_KERNEL::Exception(EXCEPTION_MSG); +} + +void MEDCouplingWithTimeStep::getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception) +{ + if(std::fabs(time-_time)<=_time_tolerance) + if(_array) + _array->getTuple(eltId,value); + else + throw INTERP_KERNEL::Exception("No array existing."); + else + throw INTERP_KERNEL::Exception(EXCEPTION_MSG); +} + +void MEDCouplingWithTimeStep::getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception) +{ + if(_dt==dt && _it==it) + if(_array) + _array->getTuple(eltId,value); + else + throw INTERP_KERNEL::Exception("No array existing."); + else + throw INTERP_KERNEL::Exception("No data on this discrete time."); +} + +MEDCouplingTwoTimeSteps::MEDCouplingTwoTimeSteps():_start_time(0.),_end_time(0.),_start_dt(-1),_end_dt(-1),_start_it(-1),_end_it(-1),_end_array(0) +{ +} + +MEDCouplingTwoTimeSteps::~MEDCouplingTwoTimeSteps() +{ + if(_end_array) + _end_array->decrRef(); +} + +void MEDCouplingTwoTimeSteps::checkNoTimePresence() const throw(INTERP_KERNEL::Exception) +{ + throw INTERP_KERNEL::Exception("The field presents a time to be specified in every access !"); +} + +void MEDCouplingTwoTimeSteps::checkTimePresence(double time) const throw(INTERP_KERNEL::Exception) +{ + if(time<_start_time-_time_tolerance || time>_end_time+_time_tolerance) + { + std::ostringstream stream; + stream << "The field is defined between times " << _start_time << " and " << _end_time << " with tolerance "; + stream << _time_tolerance << " and trying to access on time = " << time; + throw INTERP_KERNEL::Exception(stream.str().c_str()); + } +} + +void MEDCouplingTwoTimeSteps::getArrays(std::vector& arrays) const +{ + arrays.resize(2); + arrays[0]=_array; + arrays[1]=_end_array; +} diff --git a/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx b/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx new file mode 100644 index 000000000..9128f0702 --- /dev/null +++ b/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx @@ -0,0 +1,150 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __MEDCOUPLINGTIMEDISCRETIZATION_HXX__ +#define __MEDCOUPLINGTIMEDISCRETIZATION_HXX__ + +#include "MEDCoupling.hxx" +#include "RefCountObject.hxx" +#include "InterpKernelException.hxx" + +#include + +namespace ParaMEDMEM +{ + class DataArrayDouble; + class TimeLabel; + + class MEDCOUPLING_EXPORT MEDCouplingTimeDiscretization + { + protected: + MEDCouplingTimeDiscretization(); + MEDCouplingTimeDiscretization(const MEDCouplingTimeDiscretization& other, bool deepCpy); + public: + static MEDCouplingTimeDiscretization *New(TypeOfTimeDiscretization type); + virtual MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const = 0; + void setTimeTolerance(double val); + double getTimeTolerance() const { return _time_tolerance; } + virtual void checkNoTimePresence() const throw(INTERP_KERNEL::Exception) = 0; + virtual void checkTimePresence(double time) const throw(INTERP_KERNEL::Exception) = 0; + virtual void setArray(DataArrayDouble *array, TimeLabel *owner); + virtual void setArrays(const std::vector& arrays, TimeLabel *owner) throw(INTERP_KERNEL::Exception); + DataArrayDouble *getArray() const { return _array; } + virtual DataArrayDouble *getEndArray() const { return _array; } + //! Warning contrary to getArray method this method returns an object to deal with. + virtual DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception) = 0; + virtual void getArrays(std::vector& arrays) const; + virtual bool isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception); + virtual bool isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception); + double getTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { return getStartTime(dt,it); } + virtual double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) = 0; + virtual double getEndTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) = 0; + void setTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) { setStartTime(time,dt,it); } + virtual void setStartTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) = 0; + virtual void setEndTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) = 0; + virtual void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception) = 0; + virtual void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception) = 0; + virtual ~MEDCouplingTimeDiscretization(); + protected: + double _time_tolerance; + DataArrayDouble *_array; + protected: + static const double TIME_TOLERANCE_DFT; + }; + + class MEDCOUPLING_EXPORT MEDCouplingNoTimeLabel : public MEDCouplingTimeDiscretization + { + public: + MEDCouplingNoTimeLabel(); + MEDCouplingNoTimeLabel(const MEDCouplingTimeDiscretization& other, bool deepCpy); + MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const; + void checkNoTimePresence() const throw(INTERP_KERNEL::Exception) { } + void checkTimePresence(double time) const throw(INTERP_KERNEL::Exception); + DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception); + bool isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception); + bool isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception); + double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception); + double getEndTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception); + void setStartTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception); + void setEndTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception); + void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception); + void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception); + public: + static const TypeOfTimeDiscretization DISCRETIZATION=NO_TIME; + private: + static const char EXCEPTION_MSG[]; + }; + + class MEDCOUPLING_EXPORT MEDCouplingWithTimeStep : public MEDCouplingTimeDiscretization + { + protected: + MEDCouplingWithTimeStep(const MEDCouplingWithTimeStep& other, bool deepCpy); + public: + MEDCouplingWithTimeStep(); + MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const; + void checkNoTimePresence() const throw(INTERP_KERNEL::Exception); + void checkTimePresence(double time) const throw(INTERP_KERNEL::Exception); + void setStartTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) { _time=time; _dt=dt; _it=it; } + void setEndTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) { _time=time; _dt=dt; _it=it; } + double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_dt; it=_it; return _time; } + double getEndTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_dt; it=_it; return _time; } + DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception); + void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception); + void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception); + public: + static const TypeOfTimeDiscretization DISCRETIZATION=ONE_TIME; + private: + static const char EXCEPTION_MSG[]; + protected: + double _time; + int _dt; + int _it; + }; + + class MEDCOUPLING_EXPORT MEDCouplingTwoTimeSteps : public MEDCouplingTimeDiscretization + { + protected: + MEDCouplingTwoTimeSteps(); + ~MEDCouplingTwoTimeSteps(); + public: + void checkNoTimePresence() const throw(INTERP_KERNEL::Exception); + void checkTimePresence(double time) const throw(INTERP_KERNEL::Exception); + void getArrays(std::vector& arrays) const; + DataArrayDouble *getEndArray() const { return _end_array; } + void setStartTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) { _start_time=time; _start_dt=dt; _start_it=it; } + void setEndTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) { _end_time=time; _end_dt=dt; _end_it=it; } + double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_start_dt; it=_start_it; return _start_time; } + double getEndTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_end_dt; it=_end_it; return _end_time; } + protected: + double _start_time; + double _end_time; + int _start_dt; + int _end_dt; + int _start_it; + int _end_it; + DataArrayDouble *_end_array; + }; + + class MEDCOUPLING_EXPORT MEDCouplingLinearTime : public MEDCouplingTwoTimeSteps + { + public: + static const TypeOfTimeDiscretization DISCRETIZATION=LINEAR_TIME; + }; +} + +#endif diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 9c6b4fdb1..d4427d7c4 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -17,7 +17,9 @@ // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // #include "MEDCouplingUMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" #include "CellModel.hxx" +#include "VolSurfFormulae.hxx" #include @@ -37,6 +39,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::clone(bool recDeepCpy) const void MEDCouplingUMesh::updateTime() { + MEDCouplingPointSet::updateTime(); if(_nodal_connec) { updateTimeWith(*_nodal_connec); @@ -45,14 +48,10 @@ void MEDCouplingUMesh::updateTime() { updateTimeWith(*_nodal_connec_index); } - if(_coords) - { - updateTimeWith(*_coords); - } } MEDCouplingUMesh::MEDCouplingUMesh():_iterator(-1),_mesh_dim(-1), - _nodal_connec(0),_nodal_connec_index(0),_coords(0) + _nodal_connec(0),_nodal_connec_index(0) { } @@ -97,19 +96,6 @@ void MEDCouplingUMesh::allocateCells(int nbOfCells) declareAsNew(); } -void MEDCouplingUMesh::setCoords(DataArrayDouble *coords) -{ - if( coords != _coords ) - { - if (_coords) - _coords->decrRef(); - _coords=coords; - if(_coords) - _coords->incrRef(); - declareAsNew(); - } -} - void MEDCouplingUMesh::insertNextCell(INTERP_KERNEL::NormalizedCellType type, int size, const int *nodalConnOfCell) { int *pt=_nodal_connec_index->getPointer(); @@ -122,7 +108,7 @@ void MEDCouplingUMesh::insertNextCell(INTERP_KERNEL::NormalizedCellType type, in void MEDCouplingUMesh::finishInsertingCells() { - int *pt=_nodal_connec_index->getPointer(); + const int *pt=_nodal_connec_index->getConstPointer(); int idx=pt[_iterator]; _nodal_connec->reAlloc(idx); @@ -143,7 +129,7 @@ bool MEDCouplingUMesh::isEqual(const MEDCouplingMesh *other, double prec) const return false; if(_types!=otherC->_types) return false; - if(!_coords->isEqual(otherC->_coords,prec)) + if(!areCoordsEqual(*otherC,prec)) return false; return true; } @@ -159,8 +145,8 @@ void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataA int *revNodalIndxPtr=new int[nbOfNodes+1]; revNodalIndx->useArray(revNodalIndxPtr,true,CPP_DEALLOC,nbOfNodes+1,1); std::fill(revNodalIndxPtr,revNodalIndxPtr+nbOfNodes+1,0); - const int *conn=_nodal_connec->getPointer(); - const int *connIndex=_nodal_connec_index->getPointer(); + const int *conn=_nodal_connec->getConstPointer(); + const int *connIndex=_nodal_connec_index->getConstPointer(); int nbOfCells=getNumberOfCells(); int nbOfEltsInRevNodal=0; for(int eltId=0;eltIduseArray(traducer,true,CPP_DEALLOC,nbOfNodes,1); int nbOfCells=getNumberOfCells(); - const int *connIndex=_nodal_connec_index->getPointer(); + const int *connIndex=_nodal_connec_index->getConstPointer(); int *conn=_nodal_connec->getPointer(); for(int i=0;igetPointer(); + const double *oldCoordsPtr=_coords->getConstPointer(); newCoords->useArray(newCoordsPtr,true,CPP_DEALLOC,newNbOfNodes,spaceDim); int *work=std::find_if(traducer,traducer+nbOfNodes,std::bind2nd(std::not_equal_to(),-1)); for(;work!=traducer+nbOfNodes;work=std::find_if(work,traducer+nbOfNodes,std::bind2nd(std::not_equal_to(),-1))) @@ -244,7 +230,7 @@ DataArrayInt *MEDCouplingUMesh::zipCoordsTraducer() return ret; } -MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const +MEDCouplingPointSet *MEDCouplingUMesh::buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const { MEDCouplingUMesh *ret=buildPartOfMySelfKeepCoords(start,end); if(!keepCoords) @@ -252,10 +238,55 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelf(const int *start, const in return ret; } +/*! + * Given a boundary box 'bbox' returns elements 'elems' contained in this 'bbox'. + * Warning 'elems' is incremented during the call so if elems is not empty before call returned elements will be + * added in 'elems' parameter. + */ +void MEDCouplingUMesh::giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems) +{ + int dim=getSpaceDimension(); + double* elem_bb=new double[2*dim]; + const int* conn = getNodalConnectivity()->getConstPointer(); + const int* conn_index= getNodalConnectivityIndex()->getConstPointer(); + const double* coords = getCoords()->getConstPointer(); + int nbOfCells=getNumberOfCells(); + for ( int ielem=0; ielem::max(); + elem_bb[i*2+1]=-std::numeric_limits::max(); + } + + for (int inode=conn_index[ielem]+1; inode elem_bb[idim*2+1] ) + { + elem_bb[idim*2+1] = coords[node*dim+idim] ; + } + } + } + if (intersectsBoundingBox(elem_bb, bbox, dim, eps)) + { + elems.push_back(ielem); + } + } + delete [] elem_bb; +} + INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::getTypeOfCell(int cellId) const { - int *ptI=_nodal_connec_index->getPointer(); - int *pt=_nodal_connec->getPointer(); + const int *ptI=_nodal_connec_index->getConstPointer(); + const int *pt=_nodal_connec->getConstPointer(); return (INTERP_KERNEL::NormalizedCellType) pt[ptI[cellId]]; } @@ -267,36 +298,20 @@ int MEDCouplingUMesh::getNumberOfNodesInCell(int cellId) const void MEDCouplingUMesh::setConnectivity(DataArrayInt *conn, DataArrayInt *connIndex, bool isComputingTypes) { - if(_nodal_connec!=conn) - { - if(_nodal_connec) - _nodal_connec->decrRef(); - _nodal_connec=conn; - if(_nodal_connec) - _nodal_connec->incrRef(); - } - if(_nodal_connec_index!=connIndex) - { - if(_nodal_connec_index) - _nodal_connec_index->decrRef(); - _nodal_connec_index=connIndex; - if(_nodal_connec_index) - _nodal_connec_index->incrRef(); - } + DataArrayInt::setArrayIn(conn,_nodal_connec); + DataArrayInt::setArrayIn(connIndex,_nodal_connec_index); if(isComputingTypes) computeTypes(); } -MEDCouplingUMesh::MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCpy):MEDCouplingMesh(other),_iterator(-1),_mesh_dim(other._mesh_dim), - _nodal_connec(0),_nodal_connec_index(0),_coords(0), +MEDCouplingUMesh::MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCpy):MEDCouplingPointSet(other,deepCpy),_iterator(-1),_mesh_dim(other._mesh_dim), + _nodal_connec(0),_nodal_connec_index(0), _types(other._types) { if(other._nodal_connec) _nodal_connec=other._nodal_connec->performCpy(deepCpy); if(other._nodal_connec_index) _nodal_connec_index=other._nodal_connec_index->performCpy(deepCpy); - if(other._coords) - _coords=other._coords->performCpy(deepCpy); } MEDCouplingUMesh::~MEDCouplingUMesh() @@ -305,8 +320,6 @@ MEDCouplingUMesh::~MEDCouplingUMesh() _nodal_connec->decrRef(); if(_nodal_connec_index) _nodal_connec_index->decrRef(); - if(_coords) - _coords->decrRef(); } void MEDCouplingUMesh::computeTypes() @@ -314,8 +327,8 @@ void MEDCouplingUMesh::computeTypes() if(_nodal_connec && _nodal_connec_index) { _types.clear(); - const int *conn=_nodal_connec->getPointer(); - const int *connIndex=_nodal_connec_index->getPointer(); + const int *conn=_nodal_connec->getConstPointer(); + const int *connIndex=_nodal_connec_index->getConstPointer(); int nbOfElem=_nodal_connec_index->getNbOfElems()-1; for(const int *pt=connIndex;pt!=connIndex+nbOfElem;pt++) _types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]); @@ -331,11 +344,6 @@ void MEDCouplingUMesh::checkFullyDefined() const throw(INTERP_KERNEL::Exception) throw INTERP_KERNEL::Exception("Reverse nodal connectivity computation requires full connectivity and coordinates set in unstructured mesh."); } -bool MEDCouplingUMesh::isStructured() const -{ - return false; -} - int MEDCouplingUMesh::getNumberOfCells() const { if(_nodal_connec_index) @@ -344,28 +352,68 @@ int MEDCouplingUMesh::getNumberOfCells() const else return _iterator; else - throw INTERP_KERNEL::Exception("Unable to get number of cells because no coordinates specified !"); + throw INTERP_KERNEL::Exception("Unable to get number of cells because no connectivity specified !"); } -int MEDCouplingUMesh::getNumberOfNodes() const +int MEDCouplingUMesh::getMeshLength() const { - if(_coords) - return _coords->getNumberOfTuples(); - else - throw INTERP_KERNEL::Exception("Unable to get number of nodes because no coordinates specified !"); + return _nodal_connec->getNbOfElems(); } -int MEDCouplingUMesh::getSpaceDimension() const +void MEDCouplingUMesh::getTinySerializationInformation(std::vector& tinyInfo) const { - if(_coords) - return _coords->getNumberOfComponents(); - else - throw INTERP_KERNEL::Exception("Unable to get space dimension because no coordinates specified !"); + tinyInfo.resize(5); + tinyInfo[0] = getSpaceDimension(); + tinyInfo[1] = getMeshDimension(); + tinyInfo[2] = getNumberOfNodes(); + tinyInfo[3] = getNumberOfCells(); + tinyInfo[4] = getMeshLength(); } -int MEDCouplingUMesh::getMeshLength() const +/*! + * @param tinyInfo must be equal to the result given by getTinySerializationInformation method. + */ +void MEDCouplingUMesh::resizeForSerialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2) { - return _nodal_connec->getNbOfElems(); + a1->alloc(tinyInfo[4]+tinyInfo[3]+1,1); + a2->alloc(tinyInfo[2],tinyInfo[0]); +} + +void MEDCouplingUMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) +{ + a1=DataArrayInt::New(); + a1->alloc(getMeshLength()+getNumberOfCells()+1,1); + int *ptA1=a1->getPointer(); + const int *conn=getNodalConnectivity()->getConstPointer(); + const int *index=getNodalConnectivityIndex()->getConstPointer(); + ptA1=std::copy(index,index+getNumberOfCells()+1,ptA1); + std::copy(conn,conn+getMeshLength(),ptA1); + a2=getCoords(); + a2->incrRef(); +} + +/*! + * @param tinyInfo must be equal to the result given by getTinySerializationInformation method. + */ +MEDCouplingPointSet *MEDCouplingUMesh::buildObjectFromUnserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2) +{ + MEDCouplingUMesh* meshing = MEDCouplingUMesh::New() ; + // Coordinates + meshing->setCoords(a2); + // Connectivity + const int *recvBuffer=a1->getConstPointer(); + DataArrayInt* myConnecIndex=DataArrayInt::New(); + myConnecIndex->alloc(tinyInfo[3]+1,1); + std::copy(recvBuffer,recvBuffer+tinyInfo[3]+1,myConnecIndex->getPointer()); + DataArrayInt* myConnec=DataArrayInt::New(); + myConnec->alloc(tinyInfo[4],1); + std::copy(recvBuffer+tinyInfo[3]+1,recvBuffer+tinyInfo[3]+1+tinyInfo[4],myConnec->getPointer()); + meshing->setConnectivity(myConnec, myConnecIndex) ; + myConnec->decrRef(); + myConnecIndex->decrRef(); + // + meshing->setMeshDimension(tinyInfo[1]); + return meshing; } MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoords(const int *start, const int *end) const @@ -379,8 +427,8 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoords(const int *start int nbOfElemsRet=end-start; int *connIndexRet=new int[nbOfElemsRet+1]; connIndexRet[0]=0; - const int *conn=_nodal_connec->getPointer(); - const int *connIndex=_nodal_connec_index->getPointer(); + const int *conn=_nodal_connec->getConstPointer(); + const int *connIndex=_nodal_connec_index->getConstPointer(); int newNbring=0; for(const int *work=start;work!=end;work++,newNbring++) connIndexRet[newNbring+1]=connIndexRet[newNbring]+connIndex[*work+1]-connIndex[*work]; @@ -403,3 +451,196 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoords(const int *start return ret; } +/*! + * brief returns the volumes of the cells underlying the field \a field + * + * For 2D geometries, the returned field contains the areas. + * For 3D geometries, the returned field contains the volumes. + * + * param field field on which cells the volumes are required + * return field containing the volumes, area or length depending the meshdimension. + */ +MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField() const +{ + int ipt, type ; + int nbelem = getNumberOfCells() ; + int dim_mesh = getMeshDimension(); + int dim_space = getSpaceDimension() ; + const double *coords = getCoords()->getConstPointer() ; + const int *connec = getNodalConnectivity()->getConstPointer() ; + const int *connec_index = getNodalConnectivityIndex()->getConstPointer() ; + + + MEDCouplingFieldDouble* field = MEDCouplingFieldDouble::New(ON_CELLS); + DataArrayDouble* array = DataArrayDouble::New() ; + array->alloc(nbelem, 1) ; + double *area_vol = array->getPointer() ; + + switch (dim_mesh) + { + case 2: // getting the areas + for ( int iel=0 ; ielsetArray(array) ; + array->decrRef(); + field->setMesh(const_cast(this)); + return field ; +} diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 94a9f896b..714fb865b 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -20,14 +20,14 @@ #define __PARAMEDMEM_MEDCOUPLINGUMESH_HXX__ #include "MEDCoupling.hxx" -#include "MEDCouplingMesh.hxx" +#include "MEDCouplingPointSet.hxx" #include "MemArray.hxx" #include namespace ParaMEDMEM { - class MEDCOUPLING_EXPORT MEDCouplingUMesh : public MEDCouplingMesh + class MEDCOUPLING_EXPORT MEDCouplingUMesh : public MEDCouplingPointSet { public: static MEDCouplingUMesh *New(); @@ -37,8 +37,6 @@ namespace ParaMEDMEM void checkCoherency() const throw(INTERP_KERNEL::Exception); void setMeshDimension(unsigned meshDim); void allocateCells(int nbOfCells); - void setCoords(DataArrayDouble *coords); - DataArrayDouble *getCoords() const { return _coords; } void insertNextCell(INTERP_KERNEL::NormalizedCellType type, int size, const int *nodalConnOfCell); void finishInsertingCells(); const std::set& getAllTypes() const { return _types; } @@ -47,17 +45,21 @@ namespace ParaMEDMEM DataArrayInt *getNodalConnectivityIndex() const { return _nodal_connec_index; } INTERP_KERNEL::NormalizedCellType getTypeOfCell(int cellId) const; int getNumberOfNodesInCell(int cellId) const; - bool isStructured() const; int getNumberOfCells() const; - int getNumberOfNodes() const; - int getSpaceDimension() const; int getMeshDimension() const { return _mesh_dim; } int getMeshLength() const; + //! size of returned tinyInfo must be always the same. + void getTinySerializationInformation(std::vector& tinyInfo) const; + void resizeForSerialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2); + void serialize(DataArrayInt *&a1, DataArrayDouble *&a2); + MEDCouplingPointSet *buildObjectFromUnserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2); //tools void zipCoords(); DataArrayInt *zipCoordsTraducer(); void getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const; - MEDCouplingUMesh *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const; + MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const; + void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems); + MEDCouplingFieldDouble *getMeasureField() const; private: MEDCouplingUMesh(); MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCpy); @@ -72,7 +74,6 @@ namespace ParaMEDMEM unsigned _mesh_dim; DataArrayInt *_nodal_connec; DataArrayInt *_nodal_connec_index; - DataArrayDouble *_coords; std::set _types; private: static const char PART_OF_NAME[]; diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx new file mode 100644 index 000000000..71d4d1faa --- /dev/null +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx @@ -0,0 +1,230 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#include "MEDCouplingUMeshDesc.hxx" +#include "CellModel.hxx" +#include "MemArray.hxx" + +using namespace ParaMEDMEM; + +MEDCouplingUMeshDesc::MEDCouplingUMeshDesc():_mesh_dim(-1),_desc_connec(0),_desc_connec_index(0), + _nodal_connec_face(0),_nodal_connec_face_index(0) +{ +} + +MEDCouplingUMeshDesc::~MEDCouplingUMeshDesc() +{ + if(_desc_connec) + _desc_connec->decrRef(); + if(_desc_connec_index) + _desc_connec_index->decrRef(); + if(_nodal_connec_face) + _nodal_connec_face->decrRef(); + if(_nodal_connec_face_index) + _nodal_connec_face_index->decrRef(); +} + +MEDCouplingUMeshDesc *MEDCouplingUMeshDesc::New() +{ + return new MEDCouplingUMeshDesc; +} + +void MEDCouplingUMeshDesc::checkCoherency() const throw(INTERP_KERNEL::Exception) +{ + for(std::set::const_iterator iter=_types.begin();iter!=_types.end();iter++) + { + if(INTERP_KERNEL::CellModel::getCellModel(*iter).getDimension()!=_mesh_dim) + { + std::ostringstream message; + message << "MeshDesc invalid because dimension is " << _mesh_dim << " and there is presence of cell(s) with type " << (*iter); + throw INTERP_KERNEL::Exception(message.str().c_str()); + } + } +} + +void MEDCouplingUMeshDesc::setMeshDimension(unsigned meshDim) +{ + _mesh_dim=meshDim; + declareAsNew(); +} + +int MEDCouplingUMeshDesc::getNumberOfCells() const +{ + if(_desc_connec_index) + return _desc_connec_index->getNumberOfTuples()-1; + else + throw INTERP_KERNEL::Exception("Unable to get number of cells because no connectivity specified !"); +} + +int MEDCouplingUMeshDesc::getNumberOfFaces() const +{ + if(_nodal_connec_face_index) + return _nodal_connec_face_index->getNumberOfTuples()-1; + else + throw INTERP_KERNEL::Exception("Unable to get number of faces because no connectivity specified !"); +} + +int MEDCouplingUMeshDesc::getCellMeshLength() const +{ + return _desc_connec->getNbOfElems(); +} + +int MEDCouplingUMeshDesc::getFaceMeshLength() const +{ + return _nodal_connec_face->getNbOfElems(); +} + +void MEDCouplingUMeshDesc::setConnectivity(DataArrayInt *descConn, DataArrayInt *descConnIndex, DataArrayInt *nodalFaceConn, DataArrayInt *nodalFaceConnIndx) +{ + DataArrayInt::setArrayIn(descConn,_desc_connec); + DataArrayInt::setArrayIn(descConnIndex,_desc_connec_index); + DataArrayInt::setArrayIn(nodalFaceConn,_nodal_connec_face); + DataArrayInt::setArrayIn(nodalFaceConnIndx,_nodal_connec_face_index); + computeTypes(); +} + +void MEDCouplingUMeshDesc::getTinySerializationInformation(std::vector& tinyInfo) const +{ + tinyInfo.resize(7); + tinyInfo[0]=getSpaceDimension(); + tinyInfo[1]=getMeshDimension(); + tinyInfo[2]=getNumberOfNodes(); + tinyInfo[3]=getNumberOfCells(); + tinyInfo[4]=getCellMeshLength(); + tinyInfo[5]=getNumberOfFaces(); + tinyInfo[6]=getFaceMeshLength(); +} + +void MEDCouplingUMeshDesc::resizeForSerialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2) +{ + a1->alloc(tinyInfo[4]+tinyInfo[3]+1+tinyInfo[6]+tinyInfo[5]+1,1); + a2->alloc(tinyInfo[2],tinyInfo[0]); +} + +void MEDCouplingUMeshDesc::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) +{ + a1=DataArrayInt::New(); + a1->alloc(getCellMeshLength()+getNumberOfCells()+1+getFaceMeshLength()+getNumberOfFaces()+1,1); + int *ptA1=a1->getPointer(); + const int *descConn=_desc_connec->getConstPointer(); + const int *descConnIndex=_desc_connec_index->getConstPointer(); + const int *faceConn=_nodal_connec_face->getConstPointer(); + const int *faceConnIndex=_nodal_connec_face_index->getConstPointer(); + ptA1=std::copy(descConn,descConn+getCellMeshLength(),ptA1); + ptA1=std::copy(descConnIndex,descConnIndex+getNumberOfCells()+1,ptA1); + ptA1=std::copy(faceConn,faceConn+getFaceMeshLength(),ptA1); + std::copy(faceConnIndex,faceConnIndex+getNumberOfFaces()+1,ptA1); + a2=getCoords(); + a2->incrRef(); +} + +MEDCouplingPointSet *MEDCouplingUMeshDesc::buildObjectFromUnserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2) +{ + MEDCouplingUMeshDesc *meshing=MEDCouplingUMeshDesc::New(); + meshing->setCoords(a2); + const int *recvBuffer=a1->getConstPointer(); + DataArrayInt *descConn=DataArrayInt::New(); + descConn->alloc(tinyInfo[4],1); + std::copy(recvBuffer,recvBuffer+tinyInfo[4],descConn->getPointer()); + DataArrayInt *descConnIndex=DataArrayInt::New(); + descConnIndex->alloc(tinyInfo[3]+1,1); + std::copy(recvBuffer+tinyInfo[4],recvBuffer+tinyInfo[4]+tinyInfo[3]+1,descConnIndex->getPointer()); + DataArrayInt *faceConn=DataArrayInt::New(); + faceConn->alloc(tinyInfo[6],1); + std::copy(recvBuffer+tinyInfo[4]+tinyInfo[3]+1,recvBuffer+tinyInfo[4]+tinyInfo[3]+1+tinyInfo[6],faceConn->getPointer()); + DataArrayInt *faceConnIndex=DataArrayInt::New(); + faceConnIndex->alloc(tinyInfo[5]+1,1); + std::copy(recvBuffer+tinyInfo[4]+tinyInfo[3]+1+tinyInfo[6], + recvBuffer+tinyInfo[4]+tinyInfo[3]+1+tinyInfo[6]+tinyInfo[5]+1,faceConnIndex->getPointer()); + meshing->setConnectivity(descConn,descConnIndex,faceConn,faceConnIndex); + descConn->decrRef(); + descConnIndex->decrRef(); + faceConn->decrRef(); + faceConnIndex->decrRef(); + meshing->setMeshDimension(tinyInfo[1]); + return meshing; +} + +void MEDCouplingUMeshDesc::giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems) +{ + int dim=getSpaceDimension(); + double* elem_bb=new double[2*dim]; + const int* conn = _desc_connec->getConstPointer(); + const int* conn_index= _desc_connec_index->getConstPointer(); + const int* face = _nodal_connec_face->getConstPointer(); + const int* face_index= _nodal_connec_face_index->getConstPointer(); + const double* coords = getCoords()->getConstPointer(); + int nbOfCells=getNumberOfCells(); + for ( int ielem=0; ielem::max(); + elem_bb[i*2+1]=-std::numeric_limits::max(); + } + + for (int iface=conn_index[ielem]+1; iface elem_bb[idim*2+1] ) + { + elem_bb[idim*2+1] = coords[node*dim+idim] ; + } + } + } + } + if (intersectsBoundingBox(elem_bb, bbox, dim, eps)) + { + elems.push_back(ielem); + } + } + delete [] elem_bb; +} + +MEDCouplingPointSet *MEDCouplingUMeshDesc::buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const +{ + //not implemented yet. + return 0; +} + +MEDCouplingFieldDouble *MEDCouplingUMeshDesc::getMeasureField() const +{ + //not implemented yet. + return 0; +} + +void MEDCouplingUMeshDesc::computeTypes() +{ + if(_desc_connec && _desc_connec_index) + { + _types.clear(); + const int *conn=_desc_connec->getConstPointer(); + const int *connIndex=_desc_connec_index->getConstPointer(); + int nbOfElem=_desc_connec_index->getNbOfElems()-1; + for(const int *pt=connIndex;pt!=connIndex+nbOfElem;pt++) + _types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]); + } +} diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx new file mode 100644 index 000000000..04267f63a --- /dev/null +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx @@ -0,0 +1,64 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __PARAMEDMEM_MEDCOUPLINGUMESHDESC_HXX__ +#define __PARAMEDMEM_MEDCOUPLINGUMESHDESC_HXX__ + +#include "MEDCouplingPointSet.hxx" +#include "MEDCoupling.hxx" +#include "NormalizedUnstructuredMesh.hxx" + +#include + +namespace ParaMEDMEM +{ + class MEDCOUPLING_EXPORT MEDCouplingUMeshDesc : public MEDCouplingPointSet + { + public: + static MEDCouplingUMeshDesc *New(); + void checkCoherency() const throw(INTERP_KERNEL::Exception); + void setMeshDimension(unsigned meshDim); + int getNumberOfCells() const; + int getNumberOfFaces() const; + int getCellMeshLength() const; + int getFaceMeshLength() const; + int getMeshDimension() const { return _mesh_dim; } + void setConnectivity(DataArrayInt *descConn, DataArrayInt *descConnIndex, DataArrayInt *nodalFaceConn, DataArrayInt *nodalFaceConnIndx); + //tools to overload + void getTinySerializationInformation(std::vector& tinyInfo) const; + void resizeForSerialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2); + void serialize(DataArrayInt *&a1, DataArrayDouble *&a2); + MEDCouplingPointSet *buildObjectFromUnserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2); + void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems); + MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const; + MEDCouplingFieldDouble *getMeasureField() const; + private: + MEDCouplingUMeshDesc(); + ~MEDCouplingUMeshDesc(); + void computeTypes(); + private: + unsigned _mesh_dim; + DataArrayInt *_desc_connec; + DataArrayInt *_desc_connec_index; + DataArrayInt *_nodal_connec_face; + DataArrayInt *_nodal_connec_face_index; + std::set _types; + }; +} + +#endif diff --git a/src/MEDCoupling/Makefile.am b/src/MEDCoupling/Makefile.am index 27d0d69fc..60e83b9c2 100644 --- a/src/MEDCoupling/Makefile.am +++ b/src/MEDCoupling/Makefile.am @@ -31,14 +31,18 @@ salomeinclude_HEADERS = MEDCoupling.hxx \ MEDCouplingFieldDouble.hxx MEDCouplingMesh.hxx MEDCouplingUMesh.hxx TimeLabel.hxx \ MEDCouplingField.hxx MEDCouplingNormalizedUnstructuredMesh.hxx MemArray.hxx \ MEDCouplingNormalizedUnstructuredMesh.txx MemArray.txx RefCountObject.hxx \ -MEDCouplingSMesh.hxx +MEDCouplingRMesh.hxx MEDCouplingTimeDiscretization.hxx \ +MEDCouplingFieldDiscretization.hxx MEDCouplingPointSet.hxx \ +MEDCouplingUMeshDesc.hxx # Libraries targets dist_libmedcoupling_la_SOURCES = \ - MEDCouplingField.cxx MEDCouplingFieldDouble.cxx \ - MEDCouplingUMesh.cxx MemArray.cxx TimeLabel.cxx \ - MEDCouplingSMesh.cxx + MEDCouplingField.cxx MEDCouplingFieldDouble.cxx \ + MEDCouplingUMesh.cxx MemArray.cxx TimeLabel.cxx \ + MEDCouplingRMesh.cxx MEDCouplingTimeDiscretization.cxx \ + MEDCouplingFieldDiscretization.cxx \ + MEDCouplingPointSet.cxx MEDCouplingUMeshDesc.cxx libmedcoupling_la_LDFLAGS= @@ -56,13 +60,17 @@ endif EXTRA_DIST += \ MEDCouplingFieldDouble.hxx \ + MEDCouplingTimeDiscretization.hxx \ + MEDCouplingFieldDiscretization.hxx \ MEDCouplingMesh.hxx \ MEDCouplingUMesh.hxx \ - MEDCouplingSMesh.hxx \ + MEDCouplingRMesh.hxx \ TimeLabel.hxx \ MEDCouplingField.hxx \ MEDCouplingNormalizedUnstructuredMesh.hxx \ MemArray.hxx \ MEDCouplingNormalizedUnstructuredMesh.txx \ MemArray.txx \ - RefCountObject.hxx + RefCountObject.hxx \ + MEDCouplingPointSet.hxx \ + MEDCouplingUMeshDesc.hxx diff --git a/src/MEDCoupling/MemArray.cxx b/src/MEDCoupling/MemArray.cxx index d6741d1ab..6240b21f0 100644 --- a/src/MEDCoupling/MemArray.cxx +++ b/src/MEDCoupling/MemArray.cxx @@ -66,7 +66,19 @@ void DataArrayDouble::reAlloc(int nbOfTuples) declareAsNew(); } -void DataArrayDouble::useArray(double *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo) +void DataArrayDouble::setArrayIn(DataArrayDouble *newArray, DataArrayDouble* &arrayToSet) +{ + if(newArray!=arrayToSet) + { + if(arrayToSet) + arrayToSet->decrRef(); + arrayToSet=newArray; + if(arrayToSet) + arrayToSet->incrRef(); + } +} + +void DataArrayDouble::useArray(const double *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo) { _nb_of_tuples=nbOfTuple; _info_on_compo.resize(nbOfCompo); @@ -103,7 +115,7 @@ void DataArrayInt::alloc(int nbOfTuple, int nbOfCompo) declareAsNew(); } -void DataArrayInt::useArray(int *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo) +void DataArrayInt::useArray(const int *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo) { _nb_of_tuples=nbOfTuple; _info_on_compo.resize(nbOfCompo); @@ -117,3 +129,15 @@ void DataArrayInt::reAlloc(int nbOfTuples) _nb_of_tuples=nbOfTuples; declareAsNew(); } + +void DataArrayInt::setArrayIn(DataArrayInt *newArray, DataArrayInt* &arrayToSet) +{ + if(newArray!=arrayToSet) + { + if(arrayToSet) + arrayToSet->decrRef(); + arrayToSet=newArray; + if(arrayToSet) + arrayToSet->incrRef(); + } +} diff --git a/src/MEDCoupling/MemArray.hxx b/src/MEDCoupling/MemArray.hxx index c35e50a9f..fcfbd0292 100644 --- a/src/MEDCoupling/MemArray.hxx +++ b/src/MEDCoupling/MemArray.hxx @@ -21,25 +21,45 @@ #include "MEDCoupling.hxx" #include "RefCountObject.hxx" +#include "InterpKernelException.hxx" #include #include namespace ParaMEDMEM { + template + class MEDCouplingPointer + { + public: + MEDCouplingPointer():_internal(0),_external(0) { } + void null() { _internal=0; _external=0; } + bool isNull() const { return _internal==0 && _external==0; } + void setInternal(T *pointer); + void setExternal(const T *pointer); + const T *getConstPointer() const { if(_internal) return _internal; else return _external; } + const T *getConstPointerLoc(int offset) const { if(_internal) return _internal+offset; else return _external+offset; } + T *getPointer() const { if(_internal) return _internal; throw INTERP_KERNEL::Exception("Trying to write on an external pointer."); } + private: + T *_internal; + const T *_external; + }; + template class MemArray { public: - MemArray():_nb_of_elem(-1),_ownership(false),_pointer(0),_dealloc(CPP_DEALLOC) { } + MemArray():_nb_of_elem(-1),_ownership(false),_dealloc(CPP_DEALLOC) { } MemArray(const MemArray& other); - T *getPointer() const { return _pointer; } + const T *getConstPointerLoc(int offset) const { return _pointer.getConstPointerLoc(offset); } + const T *getConstPointer() const { return _pointer.getConstPointer(); } + T *getPointer() const { return _pointer.getPointer(); } MemArray &operator=(const MemArray& other); - T operator[](int id) const { return _pointer[id]; } - T& operator[](int id) { return _pointer[id]; } + T operator[](int id) const { return _pointer.getConstPointer()[id]; } + T& operator[](int id) { return _pointer.getPointer()[id]; } void alloc(int nbOfElements); void reAlloc(int newNbOfElements); - void useArray(void *array, bool ownership, DeallocType type, int nbOfElem); + void useArray(const T *array, bool ownership, DeallocType type, int nbOfElem); void writeOnPlace(int id, T element0, const T *others, int sizeOfOthers); ~MemArray() { destroy(); } private: @@ -48,7 +68,8 @@ namespace ParaMEDMEM private: int _nb_of_elem; bool _ownership; - T *_pointer; + MEDCouplingPointer _pointer; + //T *_pointer; DeallocType _dealloc; }; @@ -85,10 +106,13 @@ namespace ParaMEDMEM bool isEqual(DataArrayDouble *other, double prec) const; //!alloc or useArray should have been called before. void reAlloc(int nbOfTuples); + void getTuple(int tupleId, double *res) const { std::copy(_mem.getConstPointerLoc(tupleId*_info_on_compo.size()),_mem.getConstPointerLoc((tupleId+1)*_info_on_compo.size()),res); } double getIJ(int tupleId, int compoId) const { return _mem[tupleId*_info_on_compo.size()+compoId]; } void setIJ(int tupleId, int compoId, double newVal) { _mem[tupleId*_info_on_compo.size()+compoId]=newVal; } double *getPointer() const { return _mem.getPointer(); } - void useArray(double *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo); + static void setArrayIn(DataArrayDouble *newArray, DataArrayDouble* &arrayToSet); + const double *getConstPointer() const { return _mem.getConstPointer(); } + void useArray(const double *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo); void writeOnPlace(int id, double element0, const double *others, int sizeOfOthers) { _mem.writeOnPlace(id,element0,others,sizeOfOthers); } //! nothing to do here because this class does not aggregate any TimeLabel instance. void updateTime() { } @@ -107,10 +131,13 @@ namespace ParaMEDMEM void alloc(int nbOfTuple, int nbOfCompo); //!alloc or useArray should have been called before. void reAlloc(int nbOfTuples); + void getTuple(int tupleId, int *res) const { std::copy(_mem.getConstPointerLoc(tupleId*_info_on_compo.size()),_mem.getConstPointerLoc((tupleId+1)*_info_on_compo.size()),res); } int getIJ(int tupleId, int compoId) const { return _mem[tupleId*_info_on_compo.size()+compoId]; } void setIJ(int tupleId, int compoId, int newVal) { _mem[tupleId*_info_on_compo.size()+compoId]=newVal; } int *getPointer() const { return _mem.getPointer(); } - void useArray(int *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo); + static void setArrayIn(DataArrayInt *newArray, DataArrayInt* &arrayToSet); + const int *getConstPointer() const { return _mem.getConstPointer(); } + void useArray(const int *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo); void writeOnPlace(int id, int element0, const int *others, int sizeOfOthers) { _mem.writeOnPlace(id,element0,others,sizeOfOthers); } //! nothing to do here because this class does not aggregate any TimeLabel instance. void updateTime() { } diff --git a/src/MEDCoupling/MemArray.txx b/src/MEDCoupling/MemArray.txx index 92076ca90..c32654327 100644 --- a/src/MEDCoupling/MemArray.txx +++ b/src/MEDCoupling/MemArray.txx @@ -29,22 +29,39 @@ namespace ParaMEDMEM { template - MemArray::MemArray(const MemArray& other):_nb_of_elem(-1),_ownership(false),_pointer(0),_dealloc(CPP_DEALLOC) + void MEDCouplingPointer::setInternal(T *pointer) { - if(other._pointer) + _internal=pointer; + _external=0; + } + + template + void MEDCouplingPointer::setExternal(const T *pointer) + { + _external=pointer; + _internal=0; + } + + template + MemArray::MemArray(const MemArray& other):_nb_of_elem(-1),_ownership(false),_dealloc(CPP_DEALLOC) + { + if(!other._pointer.isNull()) { T *pointer=new T[other._nb_of_elem]; - std::copy(other._pointer,other._pointer+other._nb_of_elem,pointer); + std::copy(other._pointer.getConstPointer(),other._pointer.getConstPointer()+other._nb_of_elem,pointer); useArray(pointer,true,CPP_DEALLOC,other._nb_of_elem); } } template - void MemArray::useArray(void *array, bool ownership, DeallocType type, int nbOfElem) + void MemArray::useArray(const T *array, bool ownership, DeallocType type, int nbOfElem) { _nb_of_elem=nbOfElem; destroy(); - _pointer=(T *)array; + if(ownership) + _pointer.setInternal((T *)array); + else + _pointer.setExternal(array); _ownership=ownership; _dealloc=type; } @@ -54,8 +71,9 @@ namespace ParaMEDMEM { if(id+sizeOfOthers>=_nb_of_elem) reAlloc(2*_nb_of_elem+sizeOfOthers+1); - _pointer[id]=element0; - memcpy(_pointer+id+1,others,sizeOfOthers*sizeof(T)); + T *pointer=_pointer.getPointer(); + pointer[id]=element0; + std::copy(others,others+sizeOfOthers,pointer+id+1); } template @@ -63,7 +81,7 @@ namespace ParaMEDMEM { destroy(); _nb_of_elem=nbOfElements; - _pointer=new T[_nb_of_elem]; + _pointer.setInternal(new T[_nb_of_elem]); _ownership=true; _dealloc=CPP_DEALLOC; } @@ -72,9 +90,10 @@ namespace ParaMEDMEM void MemArray::reAlloc(int newNbOfElements) { T *pointer=new T[newNbOfElements]; - memcpy(pointer,_pointer,std::min(_nb_of_elem,newNbOfElements)*sizeof(int)); - destroyPointer(_pointer,_dealloc); - _pointer=pointer; + std::copy(_pointer.getConstPointer(),_pointer.getConstPointer()+std::min(_nb_of_elem,newNbOfElements),pointer); + if(_ownership) + destroyPointer((T *)_pointer.getConstPointer(),_dealloc); + _pointer.setInternal(pointer); _nb_of_elem=newNbOfElements; _ownership=true; _dealloc=CPP_DEALLOC; @@ -106,8 +125,8 @@ namespace ParaMEDMEM void MemArray::destroy() { if(_ownership) - destroyPointer(_pointer,_dealloc); - _pointer=0; + destroyPointer((T *)_pointer.getConstPointer(),_dealloc); + _pointer.null(); _ownership=false; } @@ -115,7 +134,7 @@ namespace ParaMEDMEM MemArray &MemArray::operator=(const MemArray& other) { alloc(other._nb_of_elem); - memcpy(_pointer,other._pointer,_nb_of_elem*sizeof(T)); + std::copy(other._pointer.getConstPointer(),other._pointer.getConstPointer()+_nb_of_elem,_pointer.getPointer()); return *this; } } diff --git a/src/MEDCoupling/RefCountObject.hxx b/src/MEDCoupling/RefCountObject.hxx index 3afc8b98c..16a2f0b34 100644 --- a/src/MEDCoupling/RefCountObject.hxx +++ b/src/MEDCoupling/RefCountObject.hxx @@ -35,6 +35,13 @@ namespace ParaMEDMEM ON_NODES = 1 } TypeOfField; + typedef enum + { + NO_TIME = 4, + ONE_TIME = 5, + LINEAR_TIME = 6 + } TypeOfTimeDiscretization; + class RefCountObject : public TimeLabel { protected: diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx index 21b4c6407..15087df5e 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx @@ -22,6 +22,7 @@ #include "MemArray.hxx" #include "Interpolation2D.txx" #include "Interpolation3DSurf.txx" +#include "Interpolation3D.txx" #include "MEDCouplingNormalizedUnstructuredMesh.txx" @@ -30,6 +31,36 @@ using namespace std; using namespace ParaMEDMEM; +void MEDCouplingBasicsTest::testArray() +{ + int tmp1[6]={7,6,5,4,3,2}; + const int tmp2[3]={8,9,10}; + { + MemArray mem; + mem.useArray(tmp1,false,CPP_DEALLOC,6); + CPPUNIT_ASSERT(tmp1==mem.getConstPointer()); + CPPUNIT_ASSERT_THROW(mem.getPointer(),INTERP_KERNEL::Exception); + CPPUNIT_ASSERT_THROW(mem[2]=7,INTERP_KERNEL::Exception); + CPPUNIT_ASSERT_THROW(mem.writeOnPlace(0,12,tmp2,3),INTERP_KERNEL::Exception); + mem.writeOnPlace(4,12,tmp2,3); + } + { + int *tmp3=new int[6]; + std::copy(tmp1,tmp1+6,tmp3); + MemArray mem2; + mem2.useArray(tmp3,true,CPP_DEALLOC,6); + CPPUNIT_ASSERT(tmp3==mem2.getConstPointer()); + CPPUNIT_ASSERT(tmp3==mem2.getPointer()); + CPPUNIT_ASSERT_EQUAL(5,mem2[2]); + mem2[2]=7; + CPPUNIT_ASSERT_EQUAL(7,mem2[2]); + mem2.writeOnPlace(0,12,tmp2,3); + CPPUNIT_ASSERT_EQUAL(9,mem2[2]); + CPPUNIT_ASSERT_EQUAL(12,mem2[0]); + mem2.writeOnPlace(4,12,tmp2,3); + } +} + void MEDCouplingBasicsTest::testMesh() { const int nbOfCells=6; @@ -164,6 +195,37 @@ void MEDCouplingBasicsTest::testMesh() mesh->decrRef(); } +void MEDCouplingBasicsTest::testMeshPointsCloud() +{ + double targetCoords[27]={-0.3,-0.3,0.5, 0.2,-0.3,1., 0.7,-0.3,1.5, -0.3,0.2,0.5, 0.2,0.2,1., 0.7,0.2,1.5, -0.3,0.7,0.5, 0.2,0.7,1., 0.7,0.7,1.5}; + int *targetConn=0; + MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New(); + targetMesh->setMeshDimension(0); + targetMesh->allocateCells(8); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn); + targetMesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(9,3); + std::copy(targetCoords,targetCoords+27,myCoords->getPointer()); + targetMesh->setCoords(myCoords); + myCoords->decrRef(); + // + targetMesh->checkCoherency(); + CPPUNIT_ASSERT_EQUAL(3,targetMesh->getSpaceDimension()); + CPPUNIT_ASSERT_EQUAL(8,targetMesh->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(9,targetMesh->getNumberOfNodes()); + CPPUNIT_ASSERT_EQUAL(0,targetMesh->getMeshDimension()); + // + targetMesh->decrRef(); +} + void MEDCouplingBasicsTest::testDeepCopy() { DataArrayDouble *array=DataArrayDouble::New(); @@ -217,7 +279,9 @@ void MEDCouplingBasicsTest::testBuildPartOfMySelf() const int tab1[2]={0,4}; const int tab2[3]={0,2,3}; // - MEDCouplingUMesh *subMesh=mesh->buildPartOfMySelf(tab1,tab1+2,true); + MEDCouplingPointSet *subMeshSimple=mesh->buildPartOfMySelf(tab1,tab1+2,true); + MEDCouplingUMesh *subMesh=dynamic_cast(subMeshSimple); + CPPUNIT_ASSERT(subMesh); std::string name(subMesh->getName()); CPPUNIT_ASSERT_EQUAL(2,(int)mesh->getAllTypes().size()); CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,*mesh->getAllTypes().begin()); @@ -235,7 +299,9 @@ void MEDCouplingBasicsTest::testBuildPartOfMySelf() CPPUNIT_ASSERT(std::equal(subConnIndex,subConnIndex+3,subMesh->getNodalConnectivityIndex()->getPointer())); subMesh->decrRef(); // - subMesh=mesh->buildPartOfMySelf(tab2,tab2+3,true); + subMeshSimple=mesh->buildPartOfMySelf(tab2,tab2+3,true); + subMesh=dynamic_cast(subMeshSimple); + CPPUNIT_ASSERT(subMesh); name=subMesh->getName(); CPPUNIT_ASSERT_EQUAL(2,(int)subMesh->getAllTypes().size()); CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,*subMesh->getAllTypes().begin()); @@ -279,7 +345,9 @@ void MEDCouplingBasicsTest::testZipCoords() oldCoords->decrRef(); // const int tab1[2]={0,4}; - MEDCouplingUMesh *subMesh=mesh->buildPartOfMySelf(tab1,tab1+2,true); + MEDCouplingPointSet *subMeshPtSet=mesh->buildPartOfMySelf(tab1,tab1+2,true); + MEDCouplingUMesh *subMesh=dynamic_cast(subMeshPtSet); + CPPUNIT_ASSERT(subMesh); DataArrayInt *traducer=subMesh->zipCoordsTraducer(); const int expectedTraducer[9]={0,1,-1,2,3,4,-1,5,6}; CPPUNIT_ASSERT(std::equal(expectedTraducer,expectedTraducer+9,traducer->getPointer())); @@ -295,7 +363,9 @@ void MEDCouplingBasicsTest::testZipCoords() CPPUNIT_ASSERT(std::equal(subConnIndex,subConnIndex+3,subMesh->getNodalConnectivityIndex()->getPointer())); subMesh->decrRef(); // - subMesh=mesh->buildPartOfMySelf(tab1,tab1+2,false); + subMeshPtSet=mesh->buildPartOfMySelf(tab1,tab1+2,false); + subMesh=dynamic_cast(subMeshPtSet); + CPPUNIT_ASSERT(subMesh); CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_QUAD4,*subMesh->getAllTypes().begin()); CPPUNIT_ASSERT_EQUAL(2,subMesh->getNumberOfCells()); CPPUNIT_ASSERT_EQUAL(7,subMesh->getNumberOfNodes()); @@ -419,6 +489,45 @@ void MEDCouplingBasicsTest::test2DInterpP1P0_1() targetMesh->decrRef(); } +void MEDCouplingBasicsTest::test2DInterpP1P1_1() +{ + MEDCouplingUMesh *sourceMesh=build2DSourceMesh_1(); + MEDCouplingUMesh *targetMesh=build2DTargetMesh_2(); + // + MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(sourceMesh); + MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(targetMesh); + INTERP_KERNEL::Interpolation2D myInterpolator; + vector > res; + INTERP_KERNEL::IntersectionType types[2]={INTERP_KERNEL::Triangulation, INTERP_KERNEL::Geometric2D}; + for(int i=0;i<2;i++) + { + myInterpolator.setPrecision(1e-12); + myInterpolator.setIntersectionType(types[i]); + myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P1P1"); + CPPUNIT_ASSERT_EQUAL(9,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.08333333333333334,res[0][0],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.05416666666666665,res[1][0],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666666,res[1][1],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.08333333333333334,res[2][1],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.05416666666666665,res[3][0],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666668,res[3][2],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.1416666666666666,res[4][0],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02499999999999999,res[4][1],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02499999999999999,res[4][2],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09999999999999999,res[4][3],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666666,res[5][1],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09583333333333333,res[5][3],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.08333333333333333,res[6][2],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666667,res[7][2],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09583333333333331,res[7][3],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.04166666666666668,res[8][3],1.e-12); + res.clear(); + } + //clean up + sourceMesh->decrRef(); + targetMesh->decrRef(); +} + void MEDCouplingBasicsTest::test3DSurfInterpP0P0_1() { MEDCouplingUMesh *sourceMesh=build3DSurfSourceMesh_1(); @@ -522,8 +631,162 @@ void MEDCouplingBasicsTest::test3DSurfInterpP1P0_1() void MEDCouplingBasicsTest::test3DInterpP0P0_1() { MEDCouplingUMesh *sourceMesh=build3DSourceMesh_1(); + MEDCouplingUMesh *targetMesh=build3DTargetMesh_1(); + // + MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(sourceMesh); + MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh); + INTERP_KERNEL::Interpolation3D myInterpolator; + vector > res; + myInterpolator.setPrecision(1e-12); + myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0"); + CPPUNIT_ASSERT_EQUAL(8,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(8.e6,sumAll(res),1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[0][0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(41666.66666666667,res[0][6],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[0][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[0][8],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[0][10],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(41666.66666666667,res[1][2],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[1][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[1][8],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[2][0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[2][5],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[2][6],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[2][9],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[2][11],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(395833.3333333333,res[3][0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[3][2],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333331,res[3][3],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[3][5],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(395833.3333333333,res[3][8],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[4][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[4][4],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[4][6],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[4][9],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[4][10],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[5][2],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333331,res[5][3],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[5][4],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(395833.3333333333,res[5][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(395833.3333333333,res[5][10],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[6][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(250000,res[6][6],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(541666.6666666667,res[6][9],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[6][11],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(83333.33333333331,res[7][0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(479166.6666666667,res[7][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(333333.3333333333,res[7][2],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(624999.9999999997,res[7][3],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(479166.6666666667,res[7][4],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(479166.6666666667,res[7][5],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(83333.33333333333,res[7][6],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(83333.33333333331,res[7][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(83333.33333333333,res[7][8],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(83333.33333333333,res[7][9],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(83333.33333333331,res[7][10],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(479166.6666666667,res[7][11],1e-7); + //clean up + sourceMesh->decrRef(); + targetMesh->decrRef(); +} +void MEDCouplingBasicsTest::test3DInterpP0P1_1() +{ + MEDCouplingUMesh *sourceMesh=build3DTargetMesh_1(); + MEDCouplingUMesh *targetMesh=build3DSourceMesh_1(); + // + MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(sourceMesh); + MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh); + INTERP_KERNEL::Interpolation3D myInterpolator; + vector > res; + myInterpolator.setPrecision(1e-12); + myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P1"); + CPPUNIT_ASSERT_EQUAL(9,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(244444.4444444445,res[0][4],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[0][5],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(291666.6666666666,res[0][6],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[0][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(125000,res[1][0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(140277.7777777778,res[1][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(119444.4444444444,res[1][2],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[1][3],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(119444.4444444444,res[1][4],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[1][5],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(26388.88888888889,res[1][6],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(348611.1111111111,res[2][6],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888888,res[2][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(244444.4444444444,res[3][2],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333334,res[3][3],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(291666.6666666666,res[3][6],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[3][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(536111.111111111,res[4][5],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(297222.2222222221,res[4][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(223611.1111111111,res[5][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(125000,res[5][3],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(125000,res[5][5],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(26388.88888888892,res[5][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[6][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(536111.1111111109,res[7][3],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(297222.2222222221,res[7][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(11111.1111111111,res[8][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(11111.11111111111,res[8][2],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666666,res[8][3],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(11111.11111111111,res[8][4],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[8][5],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[8][6],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1466666.666666668,res[8][7],1e-7); + //clean up + sourceMesh->decrRef(); + targetMesh->decrRef(); +} + +void MEDCouplingBasicsTest::test3DInterpP1P0_1() +{ + MEDCouplingUMesh *sourceMesh=build3DSourceMesh_1(); + MEDCouplingUMesh *targetMesh=build3DTargetMesh_1(); + // + MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(sourceMesh); + MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh); + INTERP_KERNEL::Interpolation3D myInterpolator; + vector > res; + myInterpolator.setPrecision(1e-12); + myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P1P0"); + CPPUNIT_ASSERT_EQUAL(8,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(125000,res[0][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(140277.7777777778,res[1][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(223611.1111111111,res[1][5],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(11111.1111111111,res[1][8],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(119444.4444444444,res[2][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(244444.4444444445,res[2][3],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(11111.11111111111,res[2][8],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[3][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[3][3],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(125000,res[3][5],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(536111.1111111109,res[3][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[3][8],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(244444.4444444445,res[4][0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(119444.4444444445,res[4][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(11111.11111111111,res[4][8],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[5][0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[5][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(536111.1111111109,res[5][4],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(125000,res[5][5],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666666,res[5][8],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(291666.6666666666,res[6][0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(26388.88888888889,res[6][1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(348611.1111111112,res[6][2],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(291666.6666666667,res[6][3],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666666,res[6][8],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[7][0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[7][2],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[7][3],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(297222.2222222221,res[7][4],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(26388.88888888892,res[7][5],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[7][6],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(297222.2222222222,res[7][7],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1466666.666666668,res[7][8],1e-7); //clean up sourceMesh->decrRef(); + targetMesh->decrRef(); } MEDCouplingUMesh *MEDCouplingBasicsTest::build2DSourceMesh_1() @@ -565,6 +828,30 @@ MEDCouplingUMesh *MEDCouplingBasicsTest::build2DTargetMesh_1() return targetMesh; } +MEDCouplingUMesh *MEDCouplingBasicsTest::build2DTargetMesh_2() +{ + double targetCoords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 }; + int targetConn[24]={0,3,4, 0,4,1, 1,4,2, 4,5,2, 3,6,4, 6,7,4, 4,7,5, 7,8,5 }; + MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New(); + targetMesh->setMeshDimension(2); + targetMesh->allocateCells(8); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+3); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+6); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+9); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+12); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+15); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+18); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+21); + targetMesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(9,2); + std::copy(targetCoords,targetCoords+18,myCoords->getPointer()); + targetMesh->setCoords(myCoords); + myCoords->decrRef(); + return targetMesh; +} + MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfSourceMesh_1() { double sourceCoords[12]={-0.3,-0.3,0.5, 0.7,-0.3,1.5, -0.3,0.7,0.5, 0.7,0.7,1.5}; @@ -604,6 +891,30 @@ MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfTargetMesh_1() return targetMesh; } +MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfTargetMesh_2() +{ + double targetCoords[27]={-0.3,-0.3,0.5, 0.2,-0.3,1., 0.7,-0.3,1.5, -0.3,0.2,0.5, 0.2,0.2,1., 0.7,0.2,1.5, -0.3,0.7,0.5, 0.2,0.7,1., 0.7,0.7,1.5}; + int targetConn[24]={0,3,4, 0,4,1, 1,4,2, 4,5,2, 3,6,4, 6,7,4, 4,7,5, 7,8,5 }; + MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New(); + targetMesh->setMeshDimension(2); + targetMesh->allocateCells(8); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+3); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+6); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+9); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+12); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+15); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+18); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+21); + targetMesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(9,3); + std::copy(targetCoords,targetCoords+27,myCoords->getPointer()); + targetMesh->setCoords(myCoords); + myCoords->decrRef(); + return targetMesh; +} + MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSourceMesh_1() { double sourceCoords[27]={ 0.0, 0.0, 200.0, 0.0, 0.0, 0.0, 0.0, 200.0, 200.0, 0.0, 200.0, 0.0, 200.0, 0.0, 200.0, @@ -627,12 +938,33 @@ MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSourceMesh_1() sourceMesh->finishInsertingCells(); DataArrayDouble *myCoords=DataArrayDouble::New(); myCoords->alloc(9,3); - std::copy(sourceCoords,sourceCoords+12,myCoords->getPointer()); + std::copy(sourceCoords,sourceCoords+27,myCoords->getPointer()); sourceMesh->setCoords(myCoords); myCoords->decrRef(); return sourceMesh; } +MEDCouplingUMesh *MEDCouplingBasicsTest::build3DTargetMesh_1() +{ + double targetCoords[81]={ 0., 0., 0., 50., 0., 0. , 200., 0., 0. , 0., 50., 0., 50., 50., 0. , 200., 50., 0., 0., 200., 0., 50., 200., 0. , 200., 200., 0. , + 0., 0., 50., 50., 0., 50. , 200., 0., 50. , 0., 50., 50., 50., 50., 50. , 200., 50., 50., 0., 200., 50., 50., 200., 50. , 200., 200., 50. , + 0., 0., 200., 50., 0., 200. , 200., 0., 200. , 0., 50., 200., 50., 50., 200. , 200., 50., 200., 0., 200., 200., 50., 200., 200. , 200., 200., 200. }; + int targetConn[64]={0,1,4,3,9,10,13,12, 1,2,5,4,10,11,14,13, 3,4,7,6,12,13,16,15, 4,5,8,7,13,14,17,16, + 9,10,13,12,18,19,22,21, 10,11,14,13,19,20,23,22, 12,13,16,15,21,22,25,24, 13,14,17,16,22,23,26,25}; + MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New(); + targetMesh->setMeshDimension(3); + targetMesh->allocateCells(12); + for(int i=0;i<8;i++) + targetMesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,targetConn+8*i); + targetMesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(27,3); + std::copy(targetCoords,targetCoords+81,myCoords->getPointer()); + targetMesh->setCoords(myCoords); + myCoords->decrRef(); + return targetMesh; +} + double MEDCouplingBasicsTest::sumAll(const std::vector< std::map >& matrix) { double ret=0.; diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx index eccf4d090..85a903541 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx @@ -31,7 +31,9 @@ namespace ParaMEDMEM class MEDCouplingBasicsTest : public CppUnit::TestFixture { CPPUNIT_TEST_SUITE(MEDCouplingBasicsTest); + CPPUNIT_TEST( testArray ); CPPUNIT_TEST( testMesh ); + CPPUNIT_TEST( testMeshPointsCloud ); CPPUNIT_TEST( testDeepCopy ); CPPUNIT_TEST( testRevNodal ); CPPUNIT_TEST( testBuildPartOfMySelf ); @@ -40,13 +42,18 @@ namespace ParaMEDMEM CPPUNIT_TEST( test2DInterpP0P0_1 ); CPPUNIT_TEST( test2DInterpP0P1_1 ); CPPUNIT_TEST( test2DInterpP1P0_1 ); + CPPUNIT_TEST( test2DInterpP1P1_1 ); CPPUNIT_TEST( test3DSurfInterpP0P0_1 ); CPPUNIT_TEST( test3DSurfInterpP0P1_1 ); CPPUNIT_TEST( test3DSurfInterpP1P0_1 ); CPPUNIT_TEST( test3DInterpP0P0_1 ); + CPPUNIT_TEST( test3DInterpP0P1_1 ); + CPPUNIT_TEST( test3DInterpP1P0_1 ); CPPUNIT_TEST_SUITE_END(); public: + void testArray(); void testMesh(); + void testMeshPointsCloud(); void testDeepCopy(); void testRevNodal(); void testBuildPartOfMySelf(); @@ -55,16 +62,22 @@ namespace ParaMEDMEM void test2DInterpP0P0_1(); void test2DInterpP0P1_1(); void test2DInterpP1P0_1(); + void test2DInterpP1P1_1(); void test3DSurfInterpP0P0_1(); void test3DSurfInterpP0P1_1(); void test3DSurfInterpP1P0_1(); void test3DInterpP0P0_1(); + void test3DInterpP0P1_1(); + void test3DInterpP1P0_1(); private: MEDCouplingUMesh *build2DSourceMesh_1(); MEDCouplingUMesh *build2DTargetMesh_1(); + MEDCouplingUMesh *build2DTargetMesh_2(); MEDCouplingUMesh *build3DSurfSourceMesh_1(); MEDCouplingUMesh *build3DSurfTargetMesh_1(); + MEDCouplingUMesh *build3DSurfTargetMesh_2(); MEDCouplingUMesh *build3DSourceMesh_1(); + MEDCouplingUMesh *build3DTargetMesh_1(); double sumAll(const std::vector< std::map >& matrix); }; } diff --git a/src/ParaMEDMEM/BlockTopology.cxx b/src/ParaMEDMEM/BlockTopology.cxx index e2157cf6b..c8bb7c4b1 100644 --- a/src/ParaMEDMEM/BlockTopology.cxx +++ b/src/ParaMEDMEM/BlockTopology.cxx @@ -18,7 +18,7 @@ // #include "BlockTopology.hxx" #include "MemArray.hxx" -#include "MEDCouplingSMesh.hxx" +#include "MEDCouplingRMesh.hxx" #include "CommInterface.hxx" #include "ProcessorGroup.hxx" #include "MPIProcessorGroup.hxx" @@ -123,7 +123,7 @@ namespace ParaMEDMEM * instead of making the best choice with respect to the * values of the different axes. */ - BlockTopology::BlockTopology(const ProcessorGroup& group, MEDCouplingSMesh *grid): + BlockTopology::BlockTopology(const ProcessorGroup& group, MEDCouplingRMesh *grid): _proc_group(&group), _dimension(grid->getSpaceDimension()), _owns_processor_group(false) { vector axis_length(_dimension); diff --git a/src/ParaMEDMEM/BlockTopology.hxx b/src/ParaMEDMEM/BlockTopology.hxx index 4a9f3a1ad..5ce27bc13 100644 --- a/src/ParaMEDMEM/BlockTopology.hxx +++ b/src/ParaMEDMEM/BlockTopology.hxx @@ -27,7 +27,7 @@ namespace ParaMEDMEM { class ComponentTopology; - class MEDCouplingSMesh; + class MEDCouplingRMesh; typedef enum{Block,Cycle} CYCLE_TYPE; @@ -35,7 +35,7 @@ namespace ParaMEDMEM { public: BlockTopology() { } - BlockTopology(const ProcessorGroup& group, MEDCouplingSMesh *grid); + BlockTopology(const ProcessorGroup& group, MEDCouplingRMesh *grid); BlockTopology(const BlockTopology& geom_topo, const ComponentTopology& comp_topo); BlockTopology(const ProcessorGroup& group, int nb_elem); virtual ~BlockTopology(); diff --git a/src/ParaMEDMEM/DEC.cxx b/src/ParaMEDMEM/DEC.cxx index a1620c4d1..27b9aa79d 100644 --- a/src/ParaMEDMEM/DEC.cxx +++ b/src/ParaMEDMEM/DEC.cxx @@ -154,6 +154,7 @@ namespace ParaMEDMEM localgroup=_source_group; else localgroup=_target_group; + //delete _icoco_field; _icoco_field=new ICoCo::MEDField(*const_cast(triofield), *localgroup); attachLocalField(_icoco_field); return; diff --git a/src/ParaMEDMEM/ElementLocator.cxx b/src/ParaMEDMEM/ElementLocator.cxx index 2069f2805..11ba8a0ee 100644 --- a/src/ParaMEDMEM/ElementLocator.cxx +++ b/src/ParaMEDMEM/ElementLocator.cxx @@ -60,7 +60,7 @@ namespace ParaMEDMEM // between the distant ids and the ids of the local reconstruction // ========================================================================== void ElementLocator::exchangeMesh(int idistantrank, - MEDCouplingUMesh*& distant_mesh, + MEDCouplingPointSet*& distant_mesh, int*& distant_ids) { int dim = _local_cell_mesh->getSpaceDimension(); @@ -71,47 +71,11 @@ namespace ParaMEDMEM return; } - set elems; + vector elems; double* distant_bb = _domain_bounding_boxes+rank*2*dim; - double* elem_bb=new double[2*dim]; - - //defining pointers to med - const int* conn = _local_cell_mesh->getNodalConnectivity()->getPointer() ; - const int* conn_index= _local_cell_mesh->getNodalConnectivityIndex()->getPointer(); - const double* coords = _local_cell_mesh->getCoords()->getPointer() ; - - for ( int ielem=0; ielem<_local_cell_mesh->getNumberOfCells() ; ielem++) - { - for (int i=0; i::max(); - elem_bb[i*2+1]=-std::numeric_limits::max(); - } - - for (int inode=conn_index[ielem]+1; inode elem_bb[idim*2+1] ) - { - elem_bb[idim*2+1] = coords[node*dim+idim] ; - } - } - } - if (_intersectsBoundingBox(elem_bb, distant_bb, dim)) - { - elems.insert(ielem); - } - } + _local_cell_mesh->giveElemsInBoundingBox(distant_bb,getBoundingBoxAdjustment(),elems); //send_mesh contains null pointer if elems is empty - MEDCouplingUMesh* send_mesh = _meshFromElems(elems); - + MEDCouplingPointSet* send_mesh = _local_cell_mesh->buildPartOfMySelf(&elems[0],&elems[elems.size()],false); // Constituting an array containing the ids of the elements that are // going to be sent to the distant subdomain. // This array enables the correct redistribution of the data when the @@ -122,7 +86,7 @@ namespace ParaMEDMEM { distant_ids_send = new int[elems.size()]; int index=0; - for (std::set::const_iterator iter = elems.begin(); iter!= elems.end(); iter++) + for (std::vector::const_iterator iter = elems.begin(); iter!= elems.end(); iter++) { distant_ids_send[index]=*iter; index++; @@ -130,7 +94,6 @@ namespace ParaMEDMEM } _exchangeMesh(send_mesh, distant_mesh, idistantrank, distant_ids_send, distant_ids); delete[] distant_ids_send; - delete[] elem_bb; if(send_mesh) send_mesh->decrRef(); } @@ -163,30 +126,8 @@ namespace ParaMEDMEM CommInterface comm_interface =_union_group->getCommInterface(); int dim = _local_cell_mesh->getSpaceDimension(); _domain_bounding_boxes = new double[2*dim*_union_group->size()]; - const double* coords = _local_cell_mesh->getCoords()->getPointer() ; - - int nbnodes = _local_cell_mesh->getNumberOfNodes(); double * minmax=new double [2*dim]; - for (int idim=0; idim::max(); - minmax[idim*2+1]=-std::numeric_limits::max(); - } - - for (int i=0; i coords[i*dim+idim] ) - { - minmax[idim*2] = coords[i*dim+idim] ; - } - if ( minmax[idim*2+1] < coords[i*dim+idim] ) - { - minmax[idim*2+1] = coords[i*dim+idim] ; - } - } - } + _local_cell_mesh->getBoundingBox(minmax); MPIProcessorGroup* group=static_cast (_union_group); const MPI_Comm* comm = group->getComm(); @@ -226,45 +167,11 @@ namespace ParaMEDMEM return true; } - // ============================================= - // Intersect Bounding Box given 2 Bounding Boxes - // ============================================= - bool ElementLocator::_intersectsBoundingBox(double* bb1, double* bb2, int dim) - { - double bbtemp[2*dim]; - double deltamax=0.0; - double adjustment_eps=getBoundingBoxAdjustment(); - - for (int i=0; i< dim; i++) - { - double delta = bb1[2*i+1]-bb1[2*i]; - if ( delta > deltamax ) - { - deltamax = delta ; - } - // deltamax = (delta>deltamax)?delta:deltamax; - } - for (int i=0; igetNumberOfCells(); - nbconn = local_mesh->getMeshLength(); - send_buffer[0] = local_mesh->getSpaceDimension(); - send_buffer[1] = local_mesh->getMeshDimension(); - send_buffer[2] = local_mesh->getNumberOfNodes(); - send_buffer[3] = nbelems; - send_buffer[4] = nbconn; - } - else - { - for (int i=0; i<5; i++) - { - send_buffer[i]=0; - } - } - + vector tinyInfoLocal,tinyInfoDistant; + local_mesh->getTinySerializationInformation(tinyInfoLocal); + tinyInfoLocal.push_back(local_mesh->getNumberOfCells()); + tinyInfoDistant.resize(tinyInfoLocal.size()); + std::fill(tinyInfoDistant.begin(),tinyInfoDistant.end(),0); MPIProcessorGroup* group=static_cast (_union_group); - const MPI_Comm* comm=(group->getComm()); + const MPI_Comm* comm=group->getComm(); MPI_Status status; - + // iproc_distant is the number of proc in distant group // it must be converted to union numbering before communication int iprocdistant_in_union = group->translateRank(&_distant_group, iproc_distant); - - comm_interface.sendRecv(send_buffer, 5, MPI_INT, iprocdistant_in_union, 1112, - recv_buffer, 5, MPI_INT,iprocdistant_in_union,1112, + + comm_interface.sendRecv(&tinyInfoLocal[0], tinyInfoLocal.size(), MPI_INT, iprocdistant_in_union, 1112, + &tinyInfoDistant[0], tinyInfoDistant.size(), MPI_INT,iprocdistant_in_union,1112, *comm, &status); - - int distant_space_dim = recv_buffer[0]; - int distant_mesh_dim = recv_buffer[1]; - int distant_nnodes = recv_buffer[2]; - int distant_nb_elems = recv_buffer[3]; - int distant_nb_conn = recv_buffer[4]; - - delete[] send_buffer; - delete[] recv_buffer; - - // Second stage : exchanging connectivity buffers - // ---------------------------------------------- - - int nb_integers = nbconn + 2*nbelems + 1; - send_buffer = new int[nb_integers]; - const int* conn = 0; - const int* global_numbering=0; - int* ptr_buffer = send_buffer; - - if (local_mesh != 0) - { - conn = local_mesh->getNodalConnectivity()->getPointer(); - - global_numbering = local_mesh->getNodalConnectivityIndex()->getPointer() ; - - //copying the data in the integer buffer - - memcpy(ptr_buffer, global_numbering, (nbelems+1)*sizeof(int)); - ptr_buffer += nbelems+1; - memcpy(ptr_buffer,conn, nbconn*sizeof(int)); - ptr_buffer += nbconn; - memcpy(ptr_buffer, distant_ids_send, nbelems*sizeof(int)); - } - - // Preparing the recv buffers - int nb_recv_integers = distant_nb_conn + 2*distant_nb_elems + 1 ; - recv_buffer=new int[nb_recv_integers]; - - // Exchanging integer buffer - comm_interface.sendRecv(send_buffer, nb_integers, MPI_INT, + DataArrayInt *v1Local=0; + DataArrayDouble *v2Local=0; + DataArrayInt *v1Distant=DataArrayInt::New(); + DataArrayDouble *v2Distant=DataArrayDouble::New(); + local_mesh->serialize(v1Local,v2Local); + local_mesh->resizeForSerialization(tinyInfoDistant,v1Distant,v2Distant); + comm_interface.sendRecv(v1Local->getPointer(), v1Local->getNbOfElems(), MPI_INT, iprocdistant_in_union, 1111, - recv_buffer, nb_recv_integers, MPI_INT, + v1Distant->getPointer(), v1Distant->getNbOfElems(), MPI_INT, iprocdistant_in_union,1111, *comm, &status); - - if ( nb_integers>0 ) - { - delete[] send_buffer; - } - - // Third stage : exchanging coordinates - // ------------------------------------ - - int nb_recv_floats = distant_space_dim*distant_nnodes; - int nb_send_floats = 0; - double* coords=0; - - if ( local_mesh!=0 ) - { - nb_send_floats = local_mesh->getSpaceDimension() - * local_mesh->getNumberOfNodes(); - coords = local_mesh->getCoords()->getPointer(); - } - - DataArrayDouble* myCoords=DataArrayDouble::New(); - myCoords->alloc(distant_nnodes,distant_space_dim); - - comm_interface.sendRecv(coords, nb_send_floats, MPI_DOUBLE, + comm_interface.sendRecv(v2Local->getPointer(), v2Local->getNbOfElems(), MPI_DOUBLE, iprocdistant_in_union, 1112, - myCoords->getPointer(), nb_recv_floats, MPI_DOUBLE, + v2Distant->getPointer(), v2Distant->getNbOfElems(), MPI_DOUBLE, iprocdistant_in_union, 1112, - *group->getComm(), &status); - - - // Reconstructing an image of the distant mesh locally - - if ( nb_recv_integers>0 && distant_space_dim !=0 ) - { - MEDCouplingUMesh* meshing = MEDCouplingUMesh::New() ; - - // Coordinates - meshing->setCoords(myCoords) ; - myCoords->decrRef(); - // Connectivity - - int *work=recv_buffer; - DataArrayInt* myConnecIndex=DataArrayInt::New(); - myConnecIndex->alloc(distant_nb_elems+1,1); - memcpy(myConnecIndex->getPointer(), work, (distant_nb_elems+1)*sizeof(int)); - work += distant_nb_elems + 1 ; - - DataArrayInt* myConnec=DataArrayInt::New(); - myConnec->alloc(distant_nb_conn,1); - memcpy(myConnec->getPointer(), work, (distant_nb_conn)*sizeof(int)); - work+=distant_nb_conn; - meshing->setConnectivity(myConnec, myConnecIndex) ; - myConnec->decrRef(); - myConnecIndex->decrRef(); - - // correspondence between the distant ids and the ids of - // the local reconstruction - - distant_ids_recv=new int [distant_nb_elems]; - for (int i=0; isetMeshDimension(distant_mesh_dim); - - distant_mesh=meshing; - delete[] recv_buffer; - } - - } - - - // ============== - // _meshFromElems - // ============== - - MEDCouplingUMesh* ElementLocator::_meshFromElems(set& elems) - { - //returns null pointer if there are no elems in the mesh - if ( elems.size()==0 ) return 0; - - // Defining pointers - const int* conn_mesh = - const_cast (_local_cell_mesh->getNodalConnectivity()->getPointer()); - - const int* conn_index = - const_cast (_local_cell_mesh->getNodalConnectivityIndex()->getPointer()); - - const double* coords = - const_cast ( _local_cell_mesh->getCoords()->getPointer()); - - set nodes; - int nbconn=0; - for (set::const_iterator iter=elems.begin(); iter!=elems.end(); iter++) - { - // Conn_index : C-like Addresses - for (int inode=conn_index[*iter]+1; inode big2small; - int i=0; - for (set::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++) - { - big2small[*iter]=i; - i++; - } - - // Memory allocate - DataArrayInt *conn = DataArrayInt::New() ; - conn->alloc(nbconn+elems.size(),1) ; - int *connPtr=conn->getPointer(); - - DataArrayInt * connIndex = DataArrayInt::New() ; - connIndex->alloc(elems.size()+1,1) ; - int* connIndexPtr=connIndex->getPointer(); - - DataArrayDouble *new_coords = DataArrayDouble::New() ; - new_coords->alloc(nodes.size(), _local_cell_mesh->getSpaceDimension()) ; - double *new_coords_ptr = new_coords->getPointer(); - - // New connectivity table - int index=0; - int mainIndex=0; - for (set::const_iterator iter=elems.begin(); iter!=elems.end(); iter++,mainIndex++) - { - connIndexPtr[mainIndex]=index; - connPtr[index++]=conn_mesh[conn_index[*iter]]; - for (int inode = conn_index[*iter]+1; inode < conn_index[*iter+1]; inode++) - { - connPtr[index]=big2small[conn_mesh[inode]] ; // C-like number - index++; - } - } - connIndexPtr[mainIndex]=index; - // Coordinates - index=0; - for (set::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++) + *comm, &status); + if(v1Distant->getNbOfElems()>0) { - int dim = _local_cell_mesh->getSpaceDimension(); - for (int i=0; ibuildObjectFromUnserialization(tinyInfoDistant,v1Distant,v2Distant); } - - // Initialize - MEDCouplingUMesh* meshing = MEDCouplingUMesh::New() ; - meshing->setCoords(new_coords) ; - new_coords->decrRef(); - meshing->setConnectivity(conn, connIndex) ; - conn->decrRef(); - connIndex->decrRef(); - meshing->setMeshDimension(_local_cell_mesh->getMeshDimension()); - - return meshing; + distant_ids_recv=new int[tinyInfoDistant.back()]; + comm_interface.sendRecv((void *)distant_ids_send,tinyInfoLocal.back(), MPI_INT, + iprocdistant_in_union, 1113, + distant_ids_recv,tinyInfoDistant.back(), MPI_INT, + iprocdistant_in_union,1113, + *comm, &status); + v1Local->decrRef(); + v2Local->decrRef(); + v1Distant->decrRef(); + v2Distant->decrRef(); } } diff --git a/src/ParaMEDMEM/ElementLocator.hxx b/src/ParaMEDMEM/ElementLocator.hxx index eb05f12d9..82d84c234 100644 --- a/src/ParaMEDMEM/ElementLocator.hxx +++ b/src/ParaMEDMEM/ElementLocator.hxx @@ -20,7 +20,6 @@ #define __ELEMENTLOCATOR_HXX__ #include "InterpolationOptions.hxx" -#include "MEDCouplingUMesh.hxx" #include #include @@ -31,7 +30,7 @@ namespace ParaMEDMEM class ProcessorGroup; class ParaSUPPORT; class InterpolationMatrix; - + class MEDCouplingPointSet; class ElementLocator : public INTERP_KERNEL::InterpolationOptions { @@ -40,15 +39,15 @@ namespace ParaMEDMEM virtual ~ElementLocator(); void exchangeMesh(int idistantrank, - MEDCouplingUMesh*& target_mesh, + MEDCouplingPointSet*& target_mesh, int*& distant_ids); void exchangeMethod(const std::string& sourceMeth, int idistantrank, std::string& targetMeth); private: const ParaMESH& _local_para_mesh ; - MEDCouplingUMesh* _local_cell_mesh; - MEDCouplingUMesh* _local_face_mesh; - std::vector _distant_cell_meshes; - std::vector _distant_face_meshes; + MEDCouplingPointSet* _local_cell_mesh; + MEDCouplingPointSet* _local_face_mesh; + std::vector _distant_cell_meshes; + std::vector _distant_face_meshes; double* _domain_bounding_boxes; const ProcessorGroup& _distant_group; const ProcessorGroup& _local_group; @@ -57,11 +56,9 @@ namespace ParaMEDMEM void _computeBoundingBoxes(); bool _intersectsBoundingBox(int irank); - bool _intersectsBoundingBox(double* bb1, double* bb2, int dim); - void _exchangeMesh(MEDCouplingUMesh* local_mesh, MEDCouplingUMesh*& distant_mesh, + void _exchangeMesh(MEDCouplingPointSet* local_mesh, MEDCouplingPointSet*& distant_mesh, int iproc_distant, const int* distant_ids_send, int*& distant_ids_recv); - MEDCouplingUMesh* _meshFromElems(std::set& elems); }; } diff --git a/src/ParaMEDMEM/ICoCoMEDField.cxx b/src/ParaMEDMEM/ICoCoMEDField.cxx index c82a59cd1..487139e93 100644 --- a/src/ParaMEDMEM/ICoCoMEDField.cxx +++ b/src/ParaMEDMEM/ICoCoMEDField.cxx @@ -107,11 +107,10 @@ namespace ICoCo //field on the sending end if (triofield._field!=0) { - _field = new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,_mesh, *_comp_topology ); + _field = new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,ParaMEDMEM::ONE_TIME,_mesh, *_comp_topology ); ParaMEDMEM::DataArrayDouble *fieldArr=_field->getField()->getArray(); _field->getField()->setName(triofield.getName().c_str()); - _field->getField()->setTime(triofield._time1); - _field->getField()->setDtIt(0,triofield._itnumber); + _field->getField()->setTime(triofield._time1,0,triofield._itnumber); for (int i =0; igetField()->setName(triofield.getName().c_str()); - _field->getField()->setTime(triofield._time1); - _field->getField()->setDtIt(0,triofield._itnumber); + _field->getField()->setTime(triofield._time1,0,triofield._itnumber); // the trio field points to the pointer inside the MED field triofield._field=const_cast (_field->getField()->getArray()->getPointer()); for (int i=0; i > surfaces; int colSize=0; //computation of the intersection volumes between source and target elements - + MEDCouplingUMesh *distant_supportC=dynamic_cast(&distant_support); + MEDCouplingUMesh *source_supportC=dynamic_cast(_source_support); if ( distant_support.getMeshDimension() == 2 && distant_support.getSpaceDimension() == 3 ) { - MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(&distant_support); - MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(_source_support); + MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(distant_supportC); + MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(source_supportC); INTERP_KERNEL::Interpolation3DSurf interpolator (*this); colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str()); @@ -127,8 +128,8 @@ namespace ParaMEDMEM else if ( distant_support.getMeshDimension() == 2 && distant_support.getSpaceDimension() == 2) { - MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(&distant_support); - MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(_source_support); + MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC); + MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC); INTERP_KERNEL::Interpolation2D interpolator (*this); colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str()); @@ -138,8 +139,8 @@ namespace ParaMEDMEM else if ( distant_support.getMeshDimension() == 3 && distant_support.getSpaceDimension() ==3 ) { - MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(&distant_support); - MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(_source_support); + MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC); + MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC); INTERP_KERNEL::Interpolation3D interpolator (*this); colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str()); @@ -153,10 +154,8 @@ namespace ParaMEDMEM int source_size=surfaces.size(); - MEDCouplingFieldDouble *target_triangle_surf = - getSupportVolumes(&distant_support); - MEDCouplingFieldDouble *source_triangle_surf = - getSupportVolumes(_source_support) ; + MEDCouplingFieldDouble *target_triangle_surf = distant_support.getMeasureField(); + MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(); // Storing the source volumes _source_volume.resize(source_size); @@ -364,209 +363,4 @@ namespace ParaMEDMEM } } - - MEDCouplingFieldDouble* InterpolationMatrix::getSupportVolumes(MEDCouplingMesh * mesh) - { - if(!mesh->isStructured()) - return getSupportUnstructuredVolumes((MEDCouplingUMesh *)mesh); - else - throw INTERP_KERNEL::Exception("Not implemented yet !!!"); - } - - // ==================================================================== - // brief returns the volumes of the cells underlying the field \a field - - // For 2D geometries, the returned field contains the areas. - // For 3D geometries, the returned field contains the volumes. - - // param field field on which cells the volumes are required - // return field containing the volumes - // ==================================================================== - - MEDCouplingFieldDouble* InterpolationMatrix::getSupportUnstructuredVolumes(MEDCouplingUMesh * mesh) - { - int ipt, type ; - int nbelem = mesh->getNumberOfCells() ; - int dim_mesh = mesh->getMeshDimension(); - int dim_space = mesh->getSpaceDimension() ; - double *coords = mesh->getCoords()->getPointer() ; - int *connec = mesh->getNodalConnectivity()->getPointer() ; - int *connec_index = mesh->getNodalConnectivityIndex()->getPointer() ; - - - MEDCouplingFieldDouble* field = MEDCouplingFieldDouble::New(ON_CELLS); - DataArrayDouble* array = DataArrayDouble::New() ; - array->alloc(nbelem, 1) ; - double *area_vol = array->getPointer() ; - - switch (dim_mesh) - { - case 2: // getting the areas - for ( int iel=0 ; ielsetArray(array) ; - array->decrRef(); - field->setMesh(mesh) ; - - return field ; - } } diff --git a/src/ParaMEDMEM/InterpolationMatrix.hxx b/src/ParaMEDMEM/InterpolationMatrix.hxx index 1c8906984..f3ae94344 100644 --- a/src/ParaMEDMEM/InterpolationMatrix.hxx +++ b/src/ParaMEDMEM/InterpolationMatrix.hxx @@ -39,20 +39,17 @@ namespace ParaMEDMEM virtual ~InterpolationMatrix(); - void addContribution(MEDCouplingUMesh& distant_support, int iproc_distant, + void addContribution(MEDCouplingPointSet& distant_support, int iproc_distant, int* distant_elems, const std::string& srcMeth, const std::string& targetMeth); void multiply(MEDCouplingFieldDouble& field) const; void transposeMultiply(MEDCouplingFieldDouble& field)const; void prepare(); - int getNbRows() const {return _row_offsets.size();} - MPIAccessDEC* getAccessDEC(){return _mapping.getAccessDEC();} - - static MEDCouplingFieldDouble *getSupportVolumes(MEDCouplingMesh *field); - static MEDCouplingFieldDouble *getSupportUnstructuredVolumes(MEDCouplingUMesh *field); + int getNbRows() const { return _row_offsets.size(); } + MPIAccessDEC* getAccessDEC() { return _mapping.getAccessDEC(); } private: std::vector _row_offsets; std::map, int > _col_offsets; - MEDCouplingUMesh *_source_support; + MEDCouplingPointSet *_source_support; MxN_Mapping _mapping; const ProcessorGroup& _source_group; diff --git a/src/ParaMEDMEM/IntersectionDEC.cxx b/src/ParaMEDMEM/IntersectionDEC.cxx index fc209b47e..6cfa69ef2 100644 --- a/src/ParaMEDMEM/IntersectionDEC.cxx +++ b/src/ParaMEDMEM/IntersectionDEC.cxx @@ -157,7 +157,7 @@ namespace ParaMEDMEM //transfering option from IntersectionDEC to ElementLocator locator.setBoundingBoxAdjustment(getBoundingBoxAdjustment()); - MEDCouplingUMesh* distant_mesh=0; + MEDCouplingPointSet* distant_mesh=0; int* distant_ids=0; for (int i=0; i<_target_group->size(); i++) { @@ -189,7 +189,7 @@ namespace ParaMEDMEM //transfering option from IntersectionDEC to ElementLocator locator.setBoundingBoxAdjustment(getBoundingBoxAdjustment()); - MEDCouplingUMesh* distant_mesh=0; + MEDCouplingPointSet* distant_mesh=0; int* distant_ids=0; for (int i=0; i<_source_group->size(); i++) { diff --git a/src/ParaMEDMEM/MEDLoader/MEDLoader.cxx b/src/ParaMEDMEM/MEDLoader/MEDLoader.cxx index e6f2b5178..da1c1d01a 100644 --- a/src/ParaMEDMEM/MEDLoader/MEDLoader.cxx +++ b/src/ParaMEDMEM/MEDLoader/MEDLoader.cxx @@ -960,10 +960,9 @@ namespace MEDLoader ParaMEDMEM::MEDCouplingUMesh *mesh=readUMeshFromFileLev1(fileName,meshName,meshDimRelToMax,familiesToKeep,typesToKeep,meshDim); if(typeOfOutField==ON_CELLS) keepSpecifiedMeshDim(fieldPerCellType,meshDim); - ParaMEDMEM::MEDCouplingFieldDouble *ret=ParaMEDMEM::MEDCouplingFieldDouble::New(typeOfOutField); - ret->setDtIt(iteration,order); + ParaMEDMEM::MEDCouplingFieldDouble *ret=ParaMEDMEM::MEDCouplingFieldDouble::New(typeOfOutField,ONE_TIME); ret->setName(fieldName); - ret->setTime(time); + ret->setTime(time,iteration,order); ret->setMesh(mesh); mesh->decrRef(); ParaMEDMEM::DataArrayDouble *arr=buildArrayFromRawData(fieldPerCellType); @@ -1089,7 +1088,7 @@ void MEDLoader::writeParaMesh(const char *fileName, ParaMEDMEM::ParaMESH *mesh) } if(myRank==0) writeMasterFile(fileName,fileNames,mesh->getCellMesh()->getName()); - writeUMesh(fileNames[myRank].c_str(),mesh->getCellMesh()); + writeUMesh(fileNames[myRank].c_str(),dynamic_cast(mesh->getCellMesh())); } void MEDLoader::writeParaField(const char *fileName, const char *meshName, ParaMEDMEM::ParaFIELD *f) diff --git a/src/ParaMEDMEM/ParaFIELD.cxx b/src/ParaMEDMEM/ParaFIELD.cxx index a59eb9125..75eda8e0e 100644 --- a/src/ParaMEDMEM/ParaFIELD.cxx +++ b/src/ParaMEDMEM/ParaFIELD.cxx @@ -62,7 +62,7 @@ namespace ParaMEDMEM \endverbatim */ - ParaFIELD::ParaFIELD(TypeOfField type, ParaMESH* para_support, const ComponentTopology& component_topology) + ParaFIELD::ParaFIELD(TypeOfField type, TypeOfTimeDiscretization td, ParaMESH* para_support, const ComponentTopology& component_topology) :_support(para_support), _component_topology(component_topology),_topology(0), _field(0), @@ -88,7 +88,7 @@ namespace ParaMEDMEM int nb_components = component_topology.nbLocalComponents(); if (nb_components!=0) { - _field=MEDCouplingFieldDouble::New(type); + _field=MEDCouplingFieldDouble::New(type,td); _field->setMesh(_support->getCellMesh()); DataArrayDouble *array=DataArrayDouble::New(); array->alloc(_field->getNumberOfTuples(),nb_components); @@ -99,8 +99,6 @@ namespace ParaMEDMEM _field->setName("Default ParaFIELD name"); _field->setDescription("Default ParaFIELD description"); - _field->setDtIt(0,0); - _field->setTime(0.0); } /*! \brief Constructor creating the ParaFIELD @@ -172,15 +170,7 @@ namespace ParaMEDMEM double ParaFIELD::getVolumeIntegral(int icomp) const { CommInterface comm_interface = _topology->getProcGroup()->getCommInterface(); - MEDCouplingFieldDouble *volume=InterpolationMatrix::getSupportVolumes(_field->getMesh()); - double *ptr=volume->getArray()->getPointer(); - int nbOfValues=volume->getArray()->getNbOfElems(); - double integral=0.; - for (int i=0; igetIJ(i,icomp)*ptr[i]; - - volume->decrRef(); - + double integral=_field->measureAccumulate(icomp); double total=0.; const MPI_Comm* comm = (dynamic_cast(_topology->getProcGroup()))->getComm(); comm_interface.allReduce(&integral, &total, 1, MPI_DOUBLE, MPI_SUM, *comm); @@ -188,4 +178,3 @@ namespace ParaMEDMEM return total; } } - diff --git a/src/ParaMEDMEM/ParaFIELD.hxx b/src/ParaMEDMEM/ParaFIELD.hxx index 5b8944dda..1b0db272b 100644 --- a/src/ParaMEDMEM/ParaFIELD.hxx +++ b/src/ParaMEDMEM/ParaFIELD.hxx @@ -34,7 +34,7 @@ namespace ParaMEDMEM { public: - ParaFIELD(TypeOfField type, ParaMESH* mesh, const ComponentTopology& component_topology); + ParaFIELD(TypeOfField type, TypeOfTimeDiscretization td, ParaMESH* mesh, const ComponentTopology& component_topology); ParaFIELD(MEDCouplingFieldDouble* field, const ProcessorGroup& group); diff --git a/src/ParaMEDMEM/ParaGRID.cxx b/src/ParaMEDMEM/ParaGRID.cxx index 3abe033e6..ce36ce316 100644 --- a/src/ParaMEDMEM/ParaGRID.cxx +++ b/src/ParaMEDMEM/ParaGRID.cxx @@ -20,7 +20,7 @@ #include "Topology.hxx" #include "BlockTopology.hxx" #include "MemArray.hxx" -#include "MEDCouplingSMesh.hxx" +#include "MEDCouplingRMesh.hxx" #include "InterpKernelUtilities.hxx" #include @@ -30,7 +30,7 @@ using namespace std; namespace ParaMEDMEM { - ParaGRID::ParaGRID(MEDCouplingSMesh* global_grid, Topology* topology) throw(INTERP_KERNEL::Exception) + ParaGRID::ParaGRID(MEDCouplingRMesh* global_grid, Topology* topology) throw(INTERP_KERNEL::Exception) { _block_topology = dynamic_cast(topology); @@ -59,7 +59,7 @@ namespace ParaMEDMEM coordinates_names.push_back(array->getName()); coordinates_units.push_back(array->getInfoOnComponentAt(0)); } - _grid=MEDCouplingSMesh::New(); + _grid=MEDCouplingRMesh::New(); _grid->set(xyz_array, coordinates_names,coordinates_units); _grid->setName(global_grid->getName()); _grid->setDescription(global_grid->getDescription());*/ diff --git a/src/ParaMEDMEM/ParaGRID.hxx b/src/ParaMEDMEM/ParaGRID.hxx index f7abab109..8bb5a5f5f 100644 --- a/src/ParaMEDMEM/ParaGRID.hxx +++ b/src/ParaMEDMEM/ParaGRID.hxx @@ -27,17 +27,17 @@ namespace ParaMEDMEM { class Topology; class BlockTopology; - class MEDCouplingSMesh; + class MEDCouplingRMesh; class ParaGRID { public: - ParaGRID(MEDCouplingSMesh* global_grid, Topology* topology) throw(INTERP_KERNEL::Exception); + ParaGRID(MEDCouplingRMesh* global_grid, Topology* topology) throw(INTERP_KERNEL::Exception); BlockTopology * getBlockTopology() const { return _block_topology; } virtual ~ParaGRID(); - MEDCouplingSMesh* getGrid() const { return _grid; } + MEDCouplingRMesh* getGrid() const { return _grid; } private: - MEDCouplingSMesh* _grid; + MEDCouplingRMesh* _grid; // structured grid topology ParaMEDMEM::BlockTopology* _block_topology; // stores the x,y,z axes on the global grid diff --git a/src/ParaMEDMEM/ParaMESH.cxx b/src/ParaMEDMEM/ParaMESH.cxx index 438aa3894..755ab3bfb 100644 --- a/src/ParaMEDMEM/ParaMESH.cxx +++ b/src/ParaMEDMEM/ParaMESH.cxx @@ -31,7 +31,7 @@ using namespace std; namespace ParaMEDMEM { - ParaMESH::ParaMESH( MEDCouplingUMesh *subdomain_mesh, MEDCouplingUMesh *subdomain_face, + ParaMESH::ParaMESH( MEDCouplingPointSet *subdomain_mesh, MEDCouplingPointSet *subdomain_face, DataArrayInt *CorrespElt_local2global, DataArrayInt *CorrespFace_local2global, DataArrayInt *CorrespNod_local2global, const ProcessorGroup& proc_group ): _cell_mesh(subdomain_mesh), @@ -55,7 +55,7 @@ namespace ParaMEDMEM CorrespNod_local2global->incrRef(); } - ParaMESH::ParaMESH( MEDCouplingUMesh *mesh, const ProcessorGroup& proc_group, const std::string& name): + ParaMESH::ParaMESH( MEDCouplingPointSet *mesh, const ProcessorGroup& proc_group, const std::string& name): _cell_mesh(mesh), _face_mesh(0), _my_domain_id(proc_group.myRank()), diff --git a/src/ParaMEDMEM/ParaMESH.hxx b/src/ParaMEDMEM/ParaMESH.hxx index 7e82ed3c4..bf3ef4345 100644 --- a/src/ParaMEDMEM/ParaMESH.hxx +++ b/src/ParaMEDMEM/ParaMESH.hxx @@ -19,8 +19,9 @@ #ifndef __PARAMESH_HXX__ #define __PARAMESH_HXX__ -#include "MEDCouplingUMesh.hxx" +#include "MEDCouplingPointSet.hxx" #include "ProcessorGroup.hxx" +#include "MemArray.hxx" #include #include @@ -34,20 +35,20 @@ namespace ParaMEDMEM class ParaMESH { public: - ParaMESH( MEDCouplingUMesh *subdomain_mesh, - MEDCouplingUMesh *subdomain_face, + ParaMESH( MEDCouplingPointSet *subdomain_mesh, + MEDCouplingPointSet *subdomain_face, DataArrayInt *CorrespElt_local2global, DataArrayInt *CorrespFace_local2global, DataArrayInt *CorrespNod_local2global, const ProcessorGroup& proc_group ) ; - ParaMESH( MEDCouplingUMesh *mesh, + ParaMESH( MEDCouplingPointSet *mesh, const ProcessorGroup& proc_group, const std::string& name); virtual ~ParaMESH(); Topology* getTopology() const { return _explicit_topology; } bool isStructured() const { return _cell_mesh->isStructured(); } - MEDCouplingUMesh *getCellMesh() const { return _cell_mesh; } - MEDCouplingUMesh *getFaceMesh() const { return _face_mesh; } + MEDCouplingPointSet *getCellMesh() const { return _cell_mesh; } + MEDCouplingPointSet *getFaceMesh() const { return _face_mesh; } BlockTopology* getBlockTopology() const { return _block_topology; } const int* getGlobalNumberingNode() const { return _node_global->getPointer(); } @@ -56,8 +57,8 @@ namespace ParaMEDMEM private: //mesh object underlying the ParaMESH object - MEDCouplingUMesh *_cell_mesh ; - MEDCouplingUMesh *_face_mesh ; + MEDCouplingPointSet *_cell_mesh ; + MEDCouplingPointSet *_face_mesh ; //id of the local grid int _my_domain_id; diff --git a/src/ParaMEDMEM/Test/Makefile.am b/src/ParaMEDMEM/Test/Makefile.am index b00a9498c..305298642 100644 --- a/src/ParaMEDMEM/Test/Makefile.am +++ b/src/ParaMEDMEM/Test/Makefile.am @@ -34,6 +34,7 @@ dist_libParaMEDMEMTest_la_SOURCES= \ ParaMEDMEMTest_IntersectionDEC.cxx \ ParaMEDMEMTest_StructuredCoincidentDEC.cxx \ ParaMEDMEMTest_MEDLoader.cxx \ + ParaMEDMEMTest_ICocoTrio.cxx \ MPIAccessDECTest.cxx \ test_AllToAllDEC.cxx \ test_AllToAllvDEC.cxx \ diff --git a/src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx b/src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx index 014652395..35b7e8593 100644 --- a/src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx +++ b/src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx @@ -58,6 +58,7 @@ class ParaMEDMEMTest : public CppUnit::TestFixture #endif CPPUNIT_TEST(testStructuredCoincidentDEC); CPPUNIT_TEST(testStructuredCoincidentDEC); + CPPUNIT_TEST(testICocoTrio1); CPPUNIT_TEST(testMEDLoaderRead1); CPPUNIT_TEST(testMEDLoaderPolygonRead); CPPUNIT_TEST(testMEDLoaderPolyhedronRead); @@ -99,6 +100,8 @@ public: void testAsynchronousSlowSourceIntersectionDEC_2D(); void testAsynchronousFastSourceIntersectionDEC_2D(); // + void testICocoTrio1(); + // void testMEDLoaderRead1(); void testMEDLoaderPolygonRead(); void testMEDLoaderPolyhedronRead(); diff --git a/src/ParaMEDMEM/Test/ParaMEDMEMTest_ICocoTrio.cxx b/src/ParaMEDMEM/Test/ParaMEDMEMTest_ICocoTrio.cxx new file mode 100644 index 000000000..27998113f --- /dev/null +++ b/src/ParaMEDMEM/Test/ParaMEDMEMTest_ICocoTrio.cxx @@ -0,0 +1,251 @@ +#include "ParaMEDMEMTest.hxx" +#include +#include "CommInterface.hxx" +#include "ProcessorGroup.hxx" +#include "MPIProcessorGroup.hxx" +#include "DEC.hxx" +#include "IntersectionDEC.hxx" +#include +#include +#include "ICoCoTrioField.hxx" +#include +#include + +using namespace std; +using namespace ParaMEDMEM; +using namespace ICoCo; + +typedef enum {sync_and,sync_or} synctype; +void synchronize_bool(bool& stop, synctype s) +{ + int my_stop; + int my_stop_temp = stop?1:0; + + if (s==sync_and) + MPI_Allreduce(&my_stop_temp,&my_stop,1,MPI_INTEGER,MPI_MIN,MPI_COMM_WORLD); + else if (s==sync_or) + MPI_Allreduce(&my_stop_temp,&my_stop,1,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD); + stop =(my_stop==1); +} + +void synchronize_dt(double& dt) +{ + double dttemp=dt; + MPI_Allreduce(&dttemp,&dt,1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD); +} + + +void affiche( const TrioField& field) +{ + cout < emetteur_ids; + set recepteur_ids; + emetteur_ids.insert(0); + recepteur_ids.insert(1); + + MPIProcessorGroup recepteur_group(comm,recepteur_ids); + MPIProcessorGroup emetteur_group(comm,emetteur_ids); + + + string cas; + if (recepteur_group.containsMyRank()) + { + cas="recepteur"; + + } + else + cas="emetteur"; + + IntersectionDEC dec_emetteur(emetteur_group, recepteur_group); + + TrioField champ_emetteur, champ_recepteur; + + init_triangle(champ_emetteur); + //init_triangle(champ_emetteur); + init_quad(champ_recepteur); + //init_emetteur(champ_recepteur); + + if (cas=="emetteur") + { + champ_emetteur._field=new double[champ_emetteur._nb_elems]; + for (int ele=0;elegetNumberOfCells(); @@ -170,9 +170,9 @@ void ParaMEDMEMTest::testIntersectionDEC_2D_(const char *srcMeth, const char *ta // ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT(support,*target_group); ParaMEDMEM::ComponentTopology comptopo; if(targetM=="P0") - parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo); + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); else - parafield = new ParaFIELD(ON_NODES,paramesh, comptopo); + parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); int nb_local; if(targetM=="P0") nb_local=mesh->getNumberOfCells(); @@ -341,9 +341,9 @@ void ParaMEDMEMTest::testIntersectionDEC_3D_(const char *srcMeth, const char *ta // ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT( support,*source_group); ParaMEDMEM::ComponentTopology comptopo; if(srcM=="P0") - parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo); + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); else - parafield = new ParaFIELD(ON_NODES,paramesh, comptopo); + parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); int nb_local; if(srcM=="P0") nb_local=mesh->getNumberOfCells(); @@ -374,9 +374,9 @@ void ParaMEDMEMTest::testIntersectionDEC_3D_(const char *srcMeth, const char *ta // ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT(support,*target_group); ParaMEDMEM::ComponentTopology comptopo; if(targetM=="P0") - parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo); + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); else - parafield = new ParaFIELD(ON_NODES,paramesh, comptopo); + parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); int nb_local; if(targetM=="P0") nb_local=mesh->getNumberOfCells(); @@ -597,9 +597,9 @@ void ParaMEDMEMTest::testAsynchronousIntersectionDEC_2D(double dtA, double tmaxA // ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT( support,*source_group); ParaMEDMEM::ComponentTopology comptopo; if(srcM=="P0") - parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo); + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); else - parafield = new ParaFIELD(ON_NODES,paramesh, comptopo); + parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); int nb_local; if(srcM=="P0") @@ -634,9 +634,9 @@ void ParaMEDMEMTest::testAsynchronousIntersectionDEC_2D(double dtA, double tmaxA // ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT(support,*target_group); ParaMEDMEM::ComponentTopology comptopo; if(targetM=="P0") - parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo); + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); else - parafield = new ParaFIELD(ON_NODES,paramesh, comptopo); + parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); int nb_local; if(targetM=="P0") diff --git a/src/ParaMEDMEM/Test/ParaMEDMEMTest_StructuredCoincidentDEC.cxx b/src/ParaMEDMEM/Test/ParaMEDMEMTest_StructuredCoincidentDEC.cxx index c0abcf871..e2b7d25e5 100644 --- a/src/ParaMEDMEM/Test/ParaMEDMEMTest_StructuredCoincidentDEC.cxx +++ b/src/ParaMEDMEM/Test/ParaMEDMEMTest_StructuredCoincidentDEC.cxx @@ -105,7 +105,7 @@ void ParaMEDMEMTest::testStructuredCoincidentDEC() { paramesh=new ParaMESH (mesh,source_group,"source mesh"); ParaMEDMEM::ComponentTopology comptopo(6); - parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo); + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); int nb_local=mesh->getNumberOfCells(); const int* global_numbering = paramesh->getGlobalNumberingCell(); @@ -132,7 +132,7 @@ void ParaMEDMEMTest::testStructuredCoincidentDEC() { paramesh=new ParaMESH (mesh,self_group,"target mesh"); ParaMEDMEM::ComponentTopology comptopo(6, &target_group); - parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo); + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); int nb_local=mesh->getNumberOfCells(); double *value=parafield->getField()->getArray()->getPointer(); diff --git a/src/ParaMEDMEM/Test/test_perf.cxx b/src/ParaMEDMEM/Test/test_perf.cxx index 6720d1173..db5a2c33c 100644 --- a/src/ParaMEDMEM/Test/test_perf.cxx +++ b/src/ParaMEDMEM/Test/test_perf.cxx @@ -159,7 +159,7 @@ void testIntersectionDEC_2D(const string& filename_xml1, const string& meshname1 paramesh=new ParaMESH (mesh,*source_group,"source mesh"); ParaMEDMEM::ComponentTopology comptopo; - parafield = new ParaFIELD(ON_CELLS, paramesh, comptopo); + parafield = new ParaFIELD(ON_CELLS, NO_TIME, paramesh, comptopo); int nb_local=mesh->getNumberOfCells(); double *value=parafield->getField()->getArray()->getPointer(); @@ -192,7 +192,7 @@ void testIntersectionDEC_2D(const string& filename_xml1, const string& meshname1 paramesh=new ParaMESH (mesh,*target_group,"target mesh"); ParaMEDMEM::ComponentTopology comptopo; - parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo); + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); int nb_local=mesh->getNumberOfCells(); double *value=parafield->getField()->getArray()->getPointer(); -- 2.39.2