From ca0bb8b3046c36c8b640885f8f40a15fe48be9f6 Mon Sep 17 00:00:00 2001 From: ageay Date: Tue, 19 May 2009 10:23:50 +0000 Subject: [PATCH] *** empty log message *** --- src/INTERP_KERNEL/ConvexIntersector.hxx | 2 +- src/INTERP_KERNEL/ConvexIntersector.txx | 6 +- src/INTERP_KERNEL/Geometric2DIntersector.hxx | 2 +- src/INTERP_KERNEL/Geometric2DIntersector.txx | 4 +- src/INTERP_KERNEL/Interpolation3DSurf.hxx | 15 +- src/INTERP_KERNEL/Interpolation3DSurf.txx | 13 +- src/INTERP_KERNEL/InterpolationOptions.cxx | 25 + src/INTERP_KERNEL/InterpolationOptions.hxx | 52 +- src/INTERP_KERNEL/InterpolationPlanar.txx | 12 + src/INTERP_KERNEL/Makefile.am | 3 +- src/INTERP_KERNEL/PlanarIntersector.hxx | 10 +- src/INTERP_KERNEL/PlanarIntersector.txx | 112 +++-- src/INTERP_KERNEL/PlanarIntersectorP0P0.hxx | 2 +- src/INTERP_KERNEL/PlanarIntersectorP0P0.txx | 12 +- src/INTERP_KERNEL/PlanarIntersectorP0P1.hxx | 2 +- src/INTERP_KERNEL/PlanarIntersectorP0P1.txx | 14 +- src/INTERP_KERNEL/PlanarIntersectorP1P0.hxx | 2 +- src/INTERP_KERNEL/PlanarIntersectorP1P0.txx | 9 +- src/INTERP_KERNEL/PlanarIntersectorP1P1.hxx | 2 +- src/INTERP_KERNEL/PlanarIntersectorP1P1.txx | 14 +- .../TriangulationIntersector.hxx | 2 +- .../TriangulationIntersector.txx | 4 +- src/MEDCoupling/MEDCouplingField.cxx | 6 + src/MEDCoupling/MEDCouplingField.hxx | 2 + .../MEDCouplingFieldDiscretization.cxx | 42 ++ .../MEDCouplingFieldDiscretization.hxx | 5 + src/MEDCoupling/MEDCouplingPointSet.cxx | 95 +++- src/MEDCoupling/MEDCouplingPointSet.hxx | 9 + src/MEDCoupling/MEDCouplingUMesh.cxx | 87 +++- src/MEDCoupling/MEDCouplingUMesh.hxx | 3 +- src/MEDCoupling/MEDCouplingUMeshDesc.cxx | 24 + src/MEDCoupling/MEDCouplingUMeshDesc.hxx | 4 + .../Test/MEDCouplingBasicsTest.cxx | 287 +++++++++++ .../Test/MEDCouplingBasicsTest.hxx | 10 + src/ParaMEDMEM/DEC.cxx | 11 +- src/ParaMEDMEM/DEC.hxx | 3 +- src/ParaMEDMEM/ElementLocator.cxx | 43 +- src/ParaMEDMEM/ElementLocator.hxx | 22 +- src/ParaMEDMEM/GlobalizerMesh.cxx | 143 ++++++ src/ParaMEDMEM/GlobalizerMesh.hxx | 72 +++ src/ParaMEDMEM/ICoCoMEDField.cxx | 2 +- src/ParaMEDMEM/InterpolationMatrix.cxx | 296 +++++++++--- src/ParaMEDMEM/InterpolationMatrix.hxx | 32 +- src/ParaMEDMEM/IntersectionDEC.cxx | 29 +- src/ParaMEDMEM/Makefile.am | 2 + src/ParaMEDMEM/MxN_Mapping.hxx | 3 + src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx | 4 + .../Test/ParaMEDMEMTest_ICocoTrio.cxx | 1 + .../Test/ParaMEDMEMTest_IntersectionDEC.cxx | 448 +++++++++++++++++- src/ParaMEDMEM_Swig/test_IntersectionDEC.py | 2 + 50 files changed, 1712 insertions(+), 294 deletions(-) create mode 100644 src/INTERP_KERNEL/InterpolationOptions.cxx create mode 100644 src/ParaMEDMEM/GlobalizerMesh.cxx create mode 100644 src/ParaMEDMEM/GlobalizerMesh.hxx diff --git a/src/INTERP_KERNEL/ConvexIntersector.hxx b/src/INTERP_KERNEL/ConvexIntersector.hxx index 7e1300f96..094b4c76f 100644 --- a/src/INTERP_KERNEL/ConvexIntersector.hxx +++ b/src/INTERP_KERNEL/ConvexIntersector.hxx @@ -36,7 +36,7 @@ namespace INTERP_KERNEL static const NumberingPolicy numPol=MyMeshType::My_numPol; public: ConvexIntersector(const MyMeshType& meshT, const MyMeshType& meshS, - double dimCaracteristic, double precision, double medianPlane, + double dimCaracteristic, double precision, double md3DSurf, double medianPlane, 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); diff --git a/src/INTERP_KERNEL/ConvexIntersector.txx b/src/INTERP_KERNEL/ConvexIntersector.txx index c61a7ba75..70576f39f 100644 --- a/src/INTERP_KERNEL/ConvexIntersector.txx +++ b/src/INTERP_KERNEL/ConvexIntersector.txx @@ -33,9 +33,9 @@ namespace INTERP_KERNEL { template class InterpType> ConvexIntersector::ConvexIntersector(const MyMeshType& meshT, const MyMeshType& meshS, - double dimCaracteristic, double precision, - double medianPlane, bool doRotate , int oriantation, int printLevel) - :InterpType >(meshT,meshS,dimCaracteristic, precision, medianPlane, doRotate, oriantation, printLevel), + double dimCaracteristic, double precision, double md3DSurf, + double medianPlane, bool doRotate , int oriantation, int printLevel) + :InterpType >(meshT,meshS,dimCaracteristic, precision, md3DSurf, medianPlane, doRotate, oriantation, printLevel), _epsilon(precision*dimCaracteristic) { if(PlanarIntersector::_print_level >= 1) diff --git a/src/INTERP_KERNEL/Geometric2DIntersector.hxx b/src/INTERP_KERNEL/Geometric2DIntersector.hxx index 6055bad22..257d6516b 100644 --- a/src/INTERP_KERNEL/Geometric2DIntersector.hxx +++ b/src/INTERP_KERNEL/Geometric2DIntersector.hxx @@ -38,7 +38,7 @@ namespace INTERP_KERNEL static const NumberingPolicy numPol=MyMeshType::My_numPol; public: Geometric2DIntersector(const MyMeshType& meshT, const MyMeshType& meshS, - double dimCaracteristic, double medianPlane, double precision, int orientation); + double dimCaracteristic, double md3DSurf, 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); diff --git a/src/INTERP_KERNEL/Geometric2DIntersector.txx b/src/INTERP_KERNEL/Geometric2DIntersector.txx index 6763d8a22..3592c35b2 100644 --- a/src/INTERP_KERNEL/Geometric2DIntersector.txx +++ b/src/INTERP_KERNEL/Geometric2DIntersector.txx @@ -35,8 +35,8 @@ namespace INTERP_KERNEL { template class InterpType> Geometric2DIntersector::Geometric2DIntersector(const MyMeshType& meshT, const MyMeshType& meshS, - double dimCaracteristic, double medianPlane, double precision, int orientation): - InterpType >(meshT,meshS,dimCaracteristic, precision, medianPlane, true, orientation, 0) + double dimCaracteristic, double md3DSurf, double medianPlane, double precision, int orientation): + InterpType >(meshT,meshS,dimCaracteristic, precision, md3DSurf, medianPlane, true, orientation, 0) { QUADRATIC_PLANAR::_precision=dimCaracteristic*precision; } diff --git a/src/INTERP_KERNEL/Interpolation3DSurf.hxx b/src/INTERP_KERNEL/Interpolation3DSurf.hxx index 4fff636cf..7aed59254 100644 --- a/src/INTERP_KERNEL/Interpolation3DSurf.hxx +++ b/src/INTERP_KERNEL/Interpolation3DSurf.hxx @@ -31,21 +31,10 @@ namespace INTERP_KERNEL Interpolation3DSurf(const InterpolationOptions& io); void setOptions(double precision, int printLevel, double medianPlane, IntersectionType intersectionType, bool doRotate, int orientation=0); - public: - bool doRotate() const { return _do_rotate; } - double medianPlane() const { return _median_plane; } - double surf3DAdjustmentEps() const { return _surf_3D_adjustment_eps; } - void setSurf3DAdjustmentEps(double val) { _surf_3D_adjustment_eps=val; } template - void performAdjustmentOfBB(PlanarIntersector* intersector, std::vector& bbox) const - { intersector->adjustBoundingBoxes(bbox,_surf_3D_adjustment_eps); } - protected: - bool _do_rotate; - double _median_plane; - double _surf_3D_adjustment_eps; - static const double DFT_MEDIAN_PLANE; - static const double DFT_SURF3D_ADJ_EPS; + void performAdjustmentOfBB(PlanarIntersector* intersector, std::vector& bbox) const + { intersector->adjustBoundingBoxes(bbox,InterpolationPlanar::getBoundingBoxAdjustment(),InterpolationPlanar::getBoundingBoxAdjustmentAbs()); } }; } diff --git a/src/INTERP_KERNEL/Interpolation3DSurf.txx b/src/INTERP_KERNEL/Interpolation3DSurf.txx index 8688ec85b..343067382 100644 --- a/src/INTERP_KERNEL/Interpolation3DSurf.txx +++ b/src/INTERP_KERNEL/Interpolation3DSurf.txx @@ -24,18 +24,11 @@ namespace INTERP_KERNEL { - const double Interpolation3DSurf::DFT_MEDIAN_PLANE=0.5; - const double Interpolation3DSurf::DFT_SURF3D_ADJ_EPS=1e-4; - - Interpolation3DSurf::Interpolation3DSurf():_do_rotate(true) - ,_median_plane(DFT_MEDIAN_PLANE) - ,_surf_3D_adjustment_eps(DFT_SURF3D_ADJ_EPS) + Interpolation3DSurf::Interpolation3DSurf() { } Interpolation3DSurf::Interpolation3DSurf(const InterpolationOptions& io):InterpolationPlanar(io) - ,_median_plane(io.getMedianPlane()) - ,_surf_3D_adjustment_eps(io.getBoundingBoxAdjustment()) { } @@ -63,8 +56,8 @@ namespace INTERP_KERNEL IntersectionType intersectionType, bool doRotate, int orientation) { InterpolationPlanar::setOptions(precision,printLevel,intersectionType, orientation); - _do_rotate=doRotate; - _median_plane=medianPlane; + InterpolationPlanar::setDoRotate(doRotate); + InterpolationPlanar::setMedianPlane(medianPlane); } } diff --git a/src/INTERP_KERNEL/InterpolationOptions.cxx b/src/INTERP_KERNEL/InterpolationOptions.cxx new file mode 100644 index 000000000..a23422564 --- /dev/null +++ b/src/INTERP_KERNEL/InterpolationOptions.cxx @@ -0,0 +1,25 @@ +// 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 "InterpolationOptions.hxx" + +const double INTERP_KERNEL::InterpolationOptions::DFT_MEDIAN_PLANE=0.5; + +const double INTERP_KERNEL::InterpolationOptions::DFT_SURF3D_ADJ_EPS=1.e-4; + +const double INTERP_KERNEL::InterpolationOptions::DFT_MAX_DIST_3DSURF_INTERSECT=-1.; diff --git a/src/INTERP_KERNEL/InterpolationOptions.hxx b/src/INTERP_KERNEL/InterpolationOptions.hxx index 34a7d042d..a18f76d44 100644 --- a/src/INTERP_KERNEL/InterpolationOptions.hxx +++ b/src/INTERP_KERNEL/InterpolationOptions.hxx @@ -19,24 +19,28 @@ #ifndef __INTERPOLATIONOPTIONS_HXX__ #define __INTERPOLATIONOPTIONS_HXX__ - -namespace INTERP_KERNEL { +namespace INTERP_KERNEL +{ typedef enum { Triangulation, Convex, Geometric2D } IntersectionType; /// Type describing the different ways in which the hexahedron can be split into tetrahedra. /// The PLANAR_* policies persume that each face is to be considered planar, while the general /// policies make no such hypothesis. The integer at the end gives the number of tetrahedra /// that result from the split. typedef enum { PLANAR_FACE_5 = 5, PLANAR_FACE_6 = 6, GENERAL_24 = 24, GENERAL_48 = 48 } SplittingPolicy; - - class InterpolationOptions{ + class InterpolationOptions + { private : int _print_level ; IntersectionType _intersection_type; double _precision; double _median_plane ; bool _do_rotate ; + //! this measure is relative to the caracteristic dimension double _bounding_box_adjustment ; + //! this measure is absolute \b not relative to the cell size + double _bounding_box_adjustment_abs ; + double _max_distance_for_3Dsurf_intersect; int _orientation ; SplittingPolicy _splitting_policy ; @@ -45,23 +49,29 @@ namespace INTERP_KERNEL { int getPrintLevel() const { return _print_level; } void setPrintLevel(int pl) { _print_level=pl; } - IntersectionType getIntersectionType() const { return InterpolationOptions::_intersection_type; } - void setIntersectionType(IntersectionType it) { InterpolationOptions::_intersection_type=it; } + IntersectionType getIntersectionType() const { return _intersection_type; } + void setIntersectionType(IntersectionType it) { _intersection_type=it; } - double getPrecision() const { return InterpolationOptions::_precision; } - void setPrecision(double p) { InterpolationOptions::_precision=p; } + double getPrecision() const { return _precision; } + void setPrecision(double p) { _precision=p; } - double getMedianPlane() const { return InterpolationOptions::_median_plane; } - void setMedianPlane(double mp) { InterpolationOptions::_median_plane=mp; } + double getMedianPlane() const { return _median_plane; } + void setMedianPlane(double mp) { _median_plane=mp; } - bool getDoRotate() const { return InterpolationOptions::_do_rotate; } - void setDoRotate( bool dr) { InterpolationOptions::_do_rotate = dr; } + bool getDoRotate() const { return _do_rotate; } + void setDoRotate( bool dr) { _do_rotate = dr; } - double getBoundingBoxAdjustment() const { return InterpolationOptions::_bounding_box_adjustment; } - void setBoundingBoxAdjustment(double bba) { InterpolationOptions::_bounding_box_adjustment=bba; } + double getBoundingBoxAdjustment() const { return _bounding_box_adjustment; } + void setBoundingBoxAdjustment(double bba) { _bounding_box_adjustment=bba; } + + double getBoundingBoxAdjustmentAbs() const { return _bounding_box_adjustment_abs; } + void setBoundingBoxAdjustmentAbs(double bba) { _bounding_box_adjustment_abs=bba; } - int getOrientation() const { return InterpolationOptions::_orientation; } - void setOrientation(int o) { InterpolationOptions::_orientation=o; } + double getMaxDistance3DSurfIntersect() const { return _max_distance_for_3Dsurf_intersect; } + void setMaxDistance3DSurfIntersect(double bba) { _max_distance_for_3Dsurf_intersect=bba; } + + int getOrientation() const { return _orientation; } + void setOrientation(int o) { _orientation=o; } SplittingPolicy getSplittingPolicy() const { return _splitting_policy; } void setSplittingPolicy(SplittingPolicy sp) { _splitting_policy=sp; } @@ -70,12 +80,18 @@ namespace INTERP_KERNEL { _print_level=0; _intersection_type=Triangulation; _precision=1e-12; - _median_plane=0.5; + _median_plane=DFT_MEDIAN_PLANE; _do_rotate=true; - _bounding_box_adjustment=1e-4; + _bounding_box_adjustment=DFT_SURF3D_ADJ_EPS; + _bounding_box_adjustment_abs=0.; + _max_distance_for_3Dsurf_intersect=DFT_MAX_DIST_3DSURF_INTERSECT; _orientation=0; _splitting_policy=GENERAL_48; } + private: + static const double DFT_MEDIAN_PLANE; + static const double DFT_SURF3D_ADJ_EPS; + static const double DFT_MAX_DIST_3DSURF_INTERSECT; }; } diff --git a/src/INTERP_KERNEL/InterpolationPlanar.txx b/src/INTERP_KERNEL/InterpolationPlanar.txx index d48246b49..1e720fd25 100644 --- a/src/INTERP_KERNEL/InterpolationPlanar.txx +++ b/src/INTERP_KERNEL/InterpolationPlanar.txx @@ -148,6 +148,7 @@ namespace INTERP_KERNEL case Triangulation: intersector=new TriangulationIntersector(myMeshT,myMeshS,_dim_caracteristic, InterpolationOptions::getPrecision(), + InterpolationOptions::getMaxDistance3DSurfIntersect(), InterpolationOptions::getMedianPlane(), InterpolationOptions::getOrientation(), InterpolationOptions::getPrintLevel()); @@ -155,6 +156,7 @@ namespace INTERP_KERNEL case Convex: intersector=new ConvexIntersector(myMeshT,myMeshS,_dim_caracteristic, InterpolationOptions::getPrecision(), + InterpolationOptions::getMaxDistance3DSurfIntersect(), InterpolationOptions::getDoRotate(), InterpolationOptions::getMedianPlane(), InterpolationOptions::getOrientation(), @@ -162,6 +164,7 @@ namespace INTERP_KERNEL break; case Geometric2D: intersector=new Geometric2DIntersector(myMeshT, myMeshS, _dim_caracteristic, + InterpolationOptions::getMaxDistance3DSurfIntersect(), InterpolationOptions::getMedianPlane(), InterpolationOptions::getPrecision(), InterpolationOptions::getOrientation()); @@ -175,6 +178,7 @@ namespace INTERP_KERNEL case Triangulation: intersector=new TriangulationIntersector(myMeshT,myMeshS,_dim_caracteristic, InterpolationOptions::getPrecision(), + InterpolationOptions::getMaxDistance3DSurfIntersect(), InterpolationOptions::getMedianPlane(), InterpolationOptions::getOrientation(), InterpolationOptions::getPrintLevel()); @@ -182,6 +186,7 @@ namespace INTERP_KERNEL case Convex: intersector=new ConvexIntersector(myMeshT,myMeshS,_dim_caracteristic, InterpolationOptions::getPrecision(), + InterpolationOptions::getMaxDistance3DSurfIntersect(), InterpolationOptions::getDoRotate(), InterpolationOptions::getMedianPlane(), InterpolationOptions::getOrientation(), @@ -189,6 +194,7 @@ namespace INTERP_KERNEL break; case Geometric2D: intersector=new Geometric2DIntersector(myMeshT, myMeshS, _dim_caracteristic, + InterpolationOptions::getMaxDistance3DSurfIntersect(), InterpolationOptions::getMedianPlane(), InterpolationOptions::getPrecision(), InterpolationOptions::getOrientation()); @@ -202,6 +208,7 @@ namespace INTERP_KERNEL case Triangulation: intersector=new TriangulationIntersector(myMeshT,myMeshS,_dim_caracteristic, InterpolationOptions::getPrecision(), + InterpolationOptions::getMaxDistance3DSurfIntersect(), InterpolationOptions::getMedianPlane(), InterpolationOptions::getOrientation(), InterpolationOptions::getPrintLevel()); @@ -209,6 +216,7 @@ namespace INTERP_KERNEL case Convex: intersector=new ConvexIntersector(myMeshT,myMeshS,_dim_caracteristic, InterpolationOptions::getPrecision(), + InterpolationOptions::getMaxDistance3DSurfIntersect(), InterpolationOptions::getDoRotate(), InterpolationOptions::getMedianPlane(), InterpolationOptions::getOrientation(), @@ -216,6 +224,7 @@ namespace INTERP_KERNEL break; case Geometric2D: intersector=new Geometric2DIntersector(myMeshT, myMeshS, _dim_caracteristic, + InterpolationOptions::getMaxDistance3DSurfIntersect(), InterpolationOptions::getMedianPlane(), InterpolationOptions::getPrecision(), InterpolationOptions::getOrientation()); @@ -229,6 +238,7 @@ namespace INTERP_KERNEL case Triangulation: intersector=new TriangulationIntersector(myMeshT,myMeshS,_dim_caracteristic, InterpolationOptions::getPrecision(), + InterpolationOptions::getMaxDistance3DSurfIntersect(), InterpolationOptions::getMedianPlane(), InterpolationOptions::getOrientation(), InterpolationOptions::getPrintLevel()); @@ -236,6 +246,7 @@ namespace INTERP_KERNEL case Convex: intersector=new ConvexIntersector(myMeshT,myMeshS,_dim_caracteristic, InterpolationOptions::getPrecision(), + InterpolationOptions::getMaxDistance3DSurfIntersect(), InterpolationOptions::getDoRotate(), InterpolationOptions::getMedianPlane(), InterpolationOptions::getOrientation(), @@ -243,6 +254,7 @@ namespace INTERP_KERNEL break; case Geometric2D: intersector=new Geometric2DIntersector(myMeshT, myMeshS, _dim_caracteristic, + InterpolationOptions::getMaxDistance3DSurfIntersect(), InterpolationOptions::getMedianPlane(), InterpolationOptions::getPrecision(), InterpolationOptions::getOrientation()); diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index 45ed3229e..7a348f546 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -115,7 +115,8 @@ dist_libinterpkernel_la_SOURCES = \ BoundingBox.cxx \ TetraAffineTransform.cxx\ CellModel.cxx\ - UnitTetraIntersectionBary.cxx + UnitTetraIntersectionBary.cxx \ + InterpolationOptions.cxx libinterpkernel_la_CPPFLAGS=-I$(srcdir)/Geometric2D -I$(srcdir)/Bases diff --git a/src/INTERP_KERNEL/PlanarIntersector.hxx b/src/INTERP_KERNEL/PlanarIntersector.hxx index 3cce25cdc..d78414b66 100644 --- a/src/INTERP_KERNEL/PlanarIntersector.hxx +++ b/src/INTERP_KERNEL/PlanarIntersector.hxx @@ -36,11 +36,11 @@ namespace INTERP_KERNEL static const NumberingPolicy numPol=MyMeshType::My_numPol; public: //! \addtogroup InterpKerGrpIntPlan @{ - PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel); + PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel); //! @} virtual ~PlanarIntersector(); void createBoundingBoxes(const MyMeshType& mesh, std::vector& bbox); - void adjustBoundingBoxes(std::vector& bbox, double Surf3DAdjustmentEps); + void adjustBoundingBoxes(std::vector& bbox, double surf3DAdjustmentEps, double surf3DAdjustmentEpsAbs); inline void getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes); protected : int projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB); @@ -49,8 +49,9 @@ namespace INTERP_KERNEL 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); + double getValueRegardingOption(double val) const; + static int projection(double *Coords_A, double *Coords_B, + int nb_NodesA, int nb_NodesB, double epsilon, double md3DSurf, double median_plane, bool do_rotate); static void rotate3DTriangle( double* PP1, double*PP2, double*PP3, TranslationRotationMatrix& rotation_matrix); protected: @@ -63,6 +64,7 @@ namespace INTERP_KERNEL const MyMeshType& _meshT; const MyMeshType& _meshS; double _dim_caracteristic; + double _max_distance_3Dsurf_intersect; double _precision; double _median_plane; bool _do_rotate; diff --git a/src/INTERP_KERNEL/PlanarIntersector.txx b/src/INTERP_KERNEL/PlanarIntersector.txx index 4493aeb14..2a0234587 100644 --- a/src/INTERP_KERNEL/PlanarIntersector.txx +++ b/src/INTERP_KERNEL/PlanarIntersector.txx @@ -29,9 +29,9 @@ namespace INTERP_KERNEL { template - PlanarIntersector::PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel): + PlanarIntersector::PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel): _meshT(meshT),_meshS(meshS), - _dim_caracteristic(dimCaracteristic),_precision(precision),_median_plane(medianPlane), + _dim_caracteristic(dimCaracteristic),_max_distance_3Dsurf_intersect(md3DSurf),_precision(precision),_median_plane(medianPlane), _do_rotate(doRotate),_orientation(orientation),_print_level(printLevel) { _connectT=meshT.getConnectivityPtr(); @@ -127,7 +127,7 @@ namespace INTERP_KERNEL \param bbox vector containing the bounding boxes */ template - void PlanarIntersector::adjustBoundingBoxes(std::vector& bbox, double Surf3DAdjustmentEps) + void PlanarIntersector::adjustBoundingBoxes(std::vector& bbox, double surf3DAdjustmentEps, double surf3DAdjustmentEpsAbs) { /* We build the segment tree for locating possible matching intersections*/ @@ -142,8 +142,8 @@ namespace INTERP_KERNEL } for(int idim=0; idim void PlanarIntersector::getRealCoordinates(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, std::vector& coordsT, std::vector& coordsS, int& orientation) { - /*double epsilon=_precision*_dim_caracteristic; - coordsT.resize(SPACEDIM*nbNodesT); - coordsS.resize(SPACEDIM*nbNodesS); - int nb_dist_NodesT=nbNodesT; - int nb_dist_NodesS=nbNodesS; - int i_last = nbNodesT - 1; - const double * Pi_last=_coordsT +_connectT[OTT::conn2C(_connIndexT[OTT::ind2C(icellT)]+i_last)]; - - for (int iT=0; iT::coo2C(_connectT[OTT::conn2C(_connIndexT[OTT::ind2C(icellT)]+iT)]); - if(distance2(Pi_last, PiT)>epsilon) - { - for (int idim=0; idim::coo2C(_connectS[OTT::conn2C(_connIndexS[OTT::ind2C(icellS)]+i_last)]); - for (int iS=0; iS::coo2C(_connectS[OTT::conn2C(_connIndexS[OTT::ind2C(icellS)]+iS)]); - if(distance2(Pi_last, PiS)>epsilon) - { - for (int idim=0; idim= 3) - { - std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl; - std::cout << std::endl << "icellT= " << icellT << ", nb nodes T= " << nbNodesT << std::endl; - for(int iT =0; iT< nbNodesT; iT++) - {for (int idim=0; idim + double PlanarIntersector::getValueRegardingOption(double val) const + { + if(_orientation==0) + return val; + if(_orientation==2) + return fabs(val); + if (( val > 0.0 && _orientation==1) || ( val < 0.0 && _orientation==-1 )) + return val; + return 0.; + } template int PlanarIntersector::projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB) { - return projection(Coords_A,Coords_B,nb_NodesA,nb_NodesB,_dim_caracteristic*_precision,_median_plane,_do_rotate); + return projection(Coords_A,Coords_B,nb_NodesA,nb_NodesB,_dim_caracteristic*_precision,_max_distance_3Dsurf_intersect,_median_plane,_do_rotate); } template int PlanarIntersector::projection(double *Coords_A, double *Coords_B, - int nb_NodesA, int nb_NodesB, double epsilon, double median_plane, bool do_rotate) + int nb_NodesA, int nb_NodesB, double epsilon, double md3DSurf, double median_plane, bool do_rotate) { double normal_A[3]={0,0,0}; double normal_B[3]={0,0,0}; @@ -344,6 +312,32 @@ namespace INTERP_KERNEL normB = sqrt(dotprod(normal_B,normal_B)); } + //fabien option + if(md3DSurf>0.) + { + double coords_GA[3]; + for (int i=0;i<3;i++) + { + coords_GA[i]=0.; + for (int j=0;jmd3DSurf) + return 0; + } if(i_A2 PlanarIntersectorP0P0::PlanarIntersectorP0P0(const MyMeshType& meshT, const MyMeshType& meshS, - double dimCaracteristic, double precision, double medianPlane, + double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel): - PlanarIntersector(meshT,meshS,dimCaracteristic,precision,medianPlane,doRotate,orientation,printLevel) + PlanarIntersector(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,doRotate,orientation,printLevel) { } @@ -52,12 +52,8 @@ namespace INTERP_KERNEL int iS=*iter; int nbNodesS=PlanarIntersector::_connIndexS[iS+1]-PlanarIntersector::_connIndexS[iS]; double surf=intersectGeometry(OTT::indFC(icellT),OTT::indFC(iS),nbNodesT,nbNodesS); - //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 )) + surf=PlanarIntersector::getValueRegardingOption(surf); + if(surf!=0.) resRow.insert(std::make_pair(OTT::indFC(iS),surf)); } } diff --git a/src/INTERP_KERNEL/PlanarIntersectorP0P1.hxx b/src/INTERP_KERNEL/PlanarIntersectorP0P1.hxx index 964133cec..a5a530216 100644 --- a/src/INTERP_KERNEL/PlanarIntersectorP0P1.hxx +++ b/src/INTERP_KERNEL/PlanarIntersectorP0P1.hxx @@ -32,7 +32,7 @@ namespace INTERP_KERNEL typedef typename MyMeshType::MyConnType ConnType; static const NumberingPolicy numPol=MyMeshType::My_numPol; protected: - PlanarIntersectorP0P1(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel); + PlanarIntersectorP0P1(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel); public: void intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res); int getNumberOfRowsOfResMatrix() const; diff --git a/src/INTERP_KERNEL/PlanarIntersectorP0P1.txx b/src/INTERP_KERNEL/PlanarIntersectorP0P1.txx index ac1969fce..f0f3a5a36 100644 --- a/src/INTERP_KERNEL/PlanarIntersectorP0P1.txx +++ b/src/INTERP_KERNEL/PlanarIntersectorP0P1.txx @@ -26,9 +26,9 @@ namespace INTERP_KERNEL { template PlanarIntersectorP0P1::PlanarIntersectorP0P1(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) + double dimCaracteristic, double precision, double md3DSurf, double medianPlane, + bool doRotate, int orientation, int printLevel): + PlanarIntersector(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,doRotate,orientation,printLevel) { } @@ -80,12 +80,8 @@ namespace INTERP_KERNEL orientation=PlanarIntersector::projectionThis(&sourceCellCoordsTmp[0],quadrangle,sourceCellCoords.size()/SPACEDIM,4); NormalizedCellType tS=PlanarIntersector::_meshS.getTypeOfElement(iS); double surf=orientation*intersectGeometryWithQuadrangle(quadrangle,sourceCellCoordsTmp,CellModel::getCellModel(tS).isQuadratic()); - //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 )) + surf=PlanarIntersector::getValueRegardingOption(surf); + if(surf!=0.) { typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT::indFC(iS)); if(iterRes==resRow.end()) diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P0.hxx b/src/INTERP_KERNEL/PlanarIntersectorP1P0.hxx index 5c2f0934c..c7b5aaf23 100644 --- a/src/INTERP_KERNEL/PlanarIntersectorP1P0.hxx +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P0.hxx @@ -32,7 +32,7 @@ namespace INTERP_KERNEL typedef typename MyMeshType::MyConnType ConnType; static const NumberingPolicy numPol=MyMeshType::My_numPol; protected: - PlanarIntersectorP1P0(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel); + PlanarIntersectorP1P0(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel); public: void intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res); int getNumberOfRowsOfResMatrix() const; diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P0.txx b/src/INTERP_KERNEL/PlanarIntersectorP1P0.txx index 209ca229c..577cb1693 100644 --- a/src/INTERP_KERNEL/PlanarIntersectorP1P0.txx +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P0.txx @@ -25,9 +25,9 @@ namespace INTERP_KERNEL { template PlanarIntersectorP1P0::PlanarIntersectorP1P0(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) + double dimCaracteristic, double precision, double md3DSurf, double medianPlane, + bool doRotate, int orientation, int printLevel): + PlanarIntersector(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,doRotate,orientation,printLevel) { } @@ -80,7 +80,8 @@ namespace INTERP_KERNEL if(SPACEDIM==3) orientation=PlanarIntersector::projectionThis(&targetCellCoordsTmp[0],quadrangle,targetCellCoords.size()/SPACEDIM,4); double surf=orientation*intersectGeometryWithQuadrangle(quadrangle,targetCellCoordsTmp,isTargetQuad); - if (( surf > 0.0 && PlanarIntersector::_orientation >=0 ) || ( surf < 0.0 && PlanarIntersector::_orientation <=0 )) + surf=PlanarIntersector::getValueRegardingOption(surf); + if(surf!=0.) { typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT::indFC(curNodeSInCmode)); if(iterRes==resRow.end()) diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P1.hxx b/src/INTERP_KERNEL/PlanarIntersectorP1P1.hxx index a4012879f..44854c8ef 100644 --- a/src/INTERP_KERNEL/PlanarIntersectorP1P1.hxx +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P1.hxx @@ -32,7 +32,7 @@ namespace INTERP_KERNEL 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); + PlanarIntersectorP1P1(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel); public: void intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res); int getNumberOfRowsOfResMatrix() const; diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P1.txx b/src/INTERP_KERNEL/PlanarIntersectorP1P1.txx index e00ec61af..66ab5e59d 100644 --- a/src/INTERP_KERNEL/PlanarIntersectorP1P1.txx +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P1.txx @@ -26,9 +26,9 @@ 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) + double dimCaracteristic, double precision, double md3DSurf, double medianPlane, + bool doRotate, int orientation, int printLevel): + PlanarIntersector(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,doRotate,orientation,printLevel) { } @@ -78,12 +78,8 @@ namespace INTERP_KERNEL 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 )) + surf=PlanarIntersector::getValueRegardingOption(surf); + if(surf!=0.) { typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT::indFC(curNodeSInCmode)); if(iterRes==resRow.end()) diff --git a/src/INTERP_KERNEL/TriangulationIntersector.hxx b/src/INTERP_KERNEL/TriangulationIntersector.hxx index e22c430f7..2f93f9c40 100644 --- a/src/INTERP_KERNEL/TriangulationIntersector.hxx +++ b/src/INTERP_KERNEL/TriangulationIntersector.hxx @@ -36,7 +36,7 @@ namespace INTERP_KERNEL static const NumberingPolicy numPol=MyMeshType::My_numPol; public: TriangulationIntersector(const MyMeshType& meshT, const MyMeshType& meshS, - double dimCaracteristic, double precision, double medianPlane, int orientation, int printLevel); + double dimCaracteristic, double precision, double md3DSurf, 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 dfdf9bc2f..ad7e91efb 100644 --- a/src/INTERP_KERNEL/TriangulationIntersector.txx +++ b/src/INTERP_KERNEL/TriangulationIntersector.txx @@ -33,9 +33,9 @@ namespace INTERP_KERNEL { template class InterpType> TriangulationIntersector::TriangulationIntersector(const MyMeshType& meshT, const MyMeshType& meshS, - double DimCaracteristic, double Precision, + double DimCaracteristic, double Precision, double md3DSurf, double MedianPlane, int orientation, int PrintLevel) - :InterpType >(meshT,meshS,DimCaracteristic, Precision, MedianPlane, true, orientation, PrintLevel) + :InterpType >(meshT,meshS,DimCaracteristic, Precision, md3DSurf, MedianPlane, true, orientation, PrintLevel) { if(PlanarIntersector::_print_level >= 1) { diff --git a/src/MEDCoupling/MEDCouplingField.cxx b/src/MEDCoupling/MEDCouplingField.cxx index c5b53a5f8..1443a82a8 100644 --- a/src/MEDCoupling/MEDCouplingField.cxx +++ b/src/MEDCoupling/MEDCouplingField.cxx @@ -85,3 +85,9 @@ MEDCouplingField::MEDCouplingField(const MEDCouplingField& other):_name(other._n _mesh->incrRef(); } } + + +MEDCouplingMesh *MEDCouplingField::buildSubMeshData(const int *start, const int *end, DataArrayInt *&di) const +{ + return _type->buildSubMeshData(start,end,_mesh,di); +} diff --git a/src/MEDCoupling/MEDCouplingField.hxx b/src/MEDCoupling/MEDCouplingField.hxx index 5f59bccfe..40a1b3857 100644 --- a/src/MEDCoupling/MEDCouplingField.hxx +++ b/src/MEDCoupling/MEDCouplingField.hxx @@ -28,6 +28,7 @@ namespace ParaMEDMEM { + class DataArrayInt; class MEDCouplingMesh; class MEDCouplingFieldDiscretization; @@ -43,6 +44,7 @@ namespace ParaMEDMEM void setDescription(const char *desc) { _desc=desc; } const char *getName() const { return _name.c_str(); } TypeOfField getEntity() const; + MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, DataArrayInt *&di) const; protected: void updateTime(); protected: diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx index 844f67eee..6987c1bff 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx @@ -20,6 +20,8 @@ #include "MEDCouplingMesh.hxx" #include "MEDCouplingFieldDouble.hxx" +#include "MEDCouplingPointSet.hxx" + using namespace ParaMEDMEM; const char MEDCouplingFieldDiscretizationP0::REPR[]="P0"; @@ -89,6 +91,21 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP0::getWeightingField(cons return mesh->getMeasureField(); } +/*! + * This method returns a submesh of 'mesh' instance constituting cell ids contained in array defined as an interval [start;end). + * @ param di is an array returned that specifies entity ids (here cells ids) in mesh 'mesh' of entity in returned submesh. + * Example : The first cell id of returned mesh has the (*di)[0] id in 'mesh' + */ +MEDCouplingMesh *MEDCouplingFieldDiscretizationP0::buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const +{ + MEDCouplingPointSet* ret=((const MEDCouplingPointSet *) mesh)->buildPartOfMySelf(start,end,false); + di=DataArrayInt::New(); + di->alloc(end-start,1); + int *pt=di->getPointer(); + std::copy(start,end,pt); + return ret; +} + TypeOfField MEDCouplingFieldDiscretizationP1::getEnum() const { return TYPE; @@ -138,3 +155,28 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getWeightingField(cons //Dual mesh to build return 0; } + +/*! + * This method invert array 'di' that is a conversion map from Old to New node numbering to New to Old node numbering. + */ +DataArrayInt *MEDCouplingFieldDiscretizationP1::invertArrayO2N2N2O(const MEDCouplingMesh *mesh, const DataArrayInt *di) +{ + DataArrayInt *ret=DataArrayInt::New(); + ret->alloc(mesh->getNumberOfNodes(),1); + int nbOfOldNodes=di->getNumberOfTuples(); + const int *old2New=di->getConstPointer(); + int *pt=ret->getPointer(); + for(int i=0;i!=nbOfOldNodes;i++) + if(old2New[i]!=-1) + pt[old2New[i]]=i; + return ret; +} + +MEDCouplingMesh *MEDCouplingFieldDiscretizationP1::buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const +{ + MEDCouplingPointSet* ret=((const MEDCouplingPointSet *) mesh)->buildPartOfMySelf(start,end,true); + DataArrayInt *diInv=ret->zipCoordsTraducer(); + di=invertArrayO2N2N2O(ret,diInv); + diInv->decrRef(); + return ret; +} diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx index 35b83f809..6061bf20f 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx @@ -26,6 +26,7 @@ namespace ParaMEDMEM { + class DataArrayInt; class MEDCouplingMesh; class DataArrayDouble; class MEDCouplingFieldDouble; @@ -42,6 +43,7 @@ namespace ParaMEDMEM virtual void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception) = 0; virtual void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) = 0; virtual MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const = 0; + virtual MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const = 0; }; class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationP0 : public MEDCouplingFieldDiscretization @@ -55,6 +57,7 @@ namespace ParaMEDMEM void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception); void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const; + MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const; public: static const char REPR[]; static const TypeOfField TYPE; @@ -71,6 +74,8 @@ namespace ParaMEDMEM void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception); void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const; + MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const; + static DataArrayInt *invertArrayO2N2N2O(const MEDCouplingMesh *mesh, const DataArrayInt *di); public: static const char REPR[]; static const TypeOfField TYPE; diff --git a/src/MEDCoupling/MEDCouplingPointSet.cxx b/src/MEDCoupling/MEDCouplingPointSet.cxx index abbc4cfd0..d4b975eab 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.cxx +++ b/src/MEDCoupling/MEDCouplingPointSet.cxx @@ -21,6 +21,10 @@ #include "MEDCouplingUMeshDesc.hxx" #include "MEDCouplingMemArray.hxx" +#include +#include +#include + using namespace ParaMEDMEM; MEDCouplingPointSet::MEDCouplingPointSet():_coords(0) @@ -118,6 +122,38 @@ void MEDCouplingPointSet::getBoundingBox(double *bbox) const } } +void MEDCouplingPointSet::zipCoords() +{ + checkFullyDefined(); + DataArrayInt *traducer=zipCoordsTraducer(); + traducer->decrRef(); +} + +void MEDCouplingPointSet::rotate(const double *center, const double *vector, double angle) +{ + int spaceDim=getSpaceDimension(); + if(spaceDim==3) + rotate3D(center,vector,angle); + else if(spaceDim==2) + rotate2D(center,angle); + else + throw INTERP_KERNEL::Exception("MEDCouplingPointSet::rotate : invalid space dim for rotation must be 2 or 3"); + _coords->declareAsNew(); + updateTime(); +} + +void MEDCouplingPointSet::translate(const double *vector) +{ + double *coords=_coords->getPointer(); + int nbNodes=getNumberOfNodes(); + int dim=getSpaceDimension(); + for(int i=0; ideclareAsNew(); + updateTime(); +} + MEDCouplingPointSet *MEDCouplingPointSet::buildInstanceFromMeshType(MEDCouplingMeshType type) { switch(type) @@ -219,8 +255,65 @@ bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* { bool intersects = (bbtemp[idim*2](),1/norm)); + //rotation matrix computation + matrix[0]=cosa; matrix[1]=0.; matrix[2]=0.; matrix[3]=0.; matrix[4]=cosa; matrix[5]=0.; matrix[6]=0.; matrix[7]=0.; matrix[8]=cosa; + matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2]; + matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2]; + matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2]; + std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies(),1-cosa)); + std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus()); + matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1]; + matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0]; + matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.; + std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies(),sina)); + std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus()); + //rotation matrix computed. + double *coords=_coords->getPointer(); + int nbNodes=getNumberOfNodes(); + double tmp[3]; + for(int i=0; i()); + coords[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0]; + coords[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1]; + coords[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2]; + } +} + +/*! + * 'This' is expected to be of spaceDim==2. Idem for 'center' and 'vect' + */ +void MEDCouplingPointSet::rotate2D(const double *center, double angle) +{ + double cosa=cos(angle); + double sina=sin(angle); + double matrix[4]; + matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa; + double *coords=_coords->getPointer(); + int nbNodes=getNumberOfNodes(); + double tmp[2]; + for(int i=0; i()); + coords[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0]; + coords[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1]; + } +} diff --git a/src/MEDCoupling/MEDCouplingPointSet.hxx b/src/MEDCoupling/MEDCouplingPointSet.hxx index 26694f413..e79f7b20b 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.hxx +++ b/src/MEDCoupling/MEDCouplingPointSet.hxx @@ -44,8 +44,13 @@ namespace ParaMEDMEM DataArrayDouble *getCoords() const { return _coords; } bool areCoordsEqual(const MEDCouplingPointSet& other, double prec) const; void getBoundingBox(double *bbox) const; + void zipCoords(); + void rotate(const double *center, const double *vector, double angle); + void translate(const double *vector); static MEDCouplingPointSet *buildInstanceFromMeshType(MEDCouplingMeshType type); virtual MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const = 0; + virtual MEDCouplingPointSet *buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const = 0; + virtual void renumberConnectivity(const int *newNodeNumbers) = 0; //! size of returned tinyInfo must be always the same. virtual void getTinySerializationInformation(std::vector& tinyInfo, std::vector& littleStrings) const; virtual void resizeForUnserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector& littleStrings); @@ -53,8 +58,12 @@ namespace ParaMEDMEM virtual void unserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, const std::vector& littleStrings); virtual void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems) = 0; + virtual DataArrayInt *zipCoordsTraducer() = 0; protected: + virtual void checkFullyDefined() const throw(INTERP_KERNEL::Exception) = 0; static bool intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps); + void rotate2D(const double *center, double angle); + void rotate3D(const double *center, const double *vect, double angle); protected: DataArrayDouble *_coords; }; diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 552a96d02..29a27c0ab 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -22,6 +22,7 @@ #include "VolSurfFormulae.hxx" #include +#include using namespace ParaMEDMEM; @@ -198,13 +199,6 @@ void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataA } } -void MEDCouplingUMesh::zipCoords() -{ - checkFullyDefined(); - DataArrayInt *traducer=zipCoordsTraducer(); - traducer->decrRef(); -} - struct MEDCouplingAccVisit { MEDCouplingAccVisit():_new_nb_of_nodes(0) { } @@ -262,6 +256,49 @@ MEDCouplingPointSet *MEDCouplingUMesh::buildPartOfMySelf(const int *start, const return ret; } +/*! + * Keeps from 'this' only cells which constituing point id is in the ids specified by 'start' and 'end'. + * The return newly allocated mesh will share the same coordinates as 'this'. + * Parameter 'fullyIn' specifies if a cell that has part of its nodes in ids array is kept or not. + * If 'fullyIn' is true only cells whose ids are \b fully contained in ['start','end') tab will be kept. + */ +MEDCouplingPointSet *MEDCouplingUMesh::buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const +{ + std::set fastFinder(start,end); + const int *conn=getNodalConnectivity()->getConstPointer(); + const int *connIndex=getNodalConnectivityIndex()->getConstPointer(); + int nbOfCells=getNumberOfCells(); + std::vector cellIdsKept; + for(int i=0;i connOfCell(conn+connIndex[i]+1,conn+connIndex[i+1]); + connOfCell.erase(-1);//polyhedron separator + int refLgth=std::min(connOfCell.size(),fastFinder.size()); + std::set locMerge; + std::insert_iterator< std::set > it(locMerge,locMerge.begin()); + std::set_intersection(connOfCell.begin(),connOfCell.end(),fastFinder.begin(),fastFinder.end(),it); + if(locMerge.size()==refLgth && fullyIn || locMerge.size()!=0 && !fullyIn) + cellIdsKept.push_back(i); + } + return buildPartOfMySelf(&cellIdsKept[0],&cellIdsKept[cellIdsKept.size()],true); +} + +void MEDCouplingUMesh::renumberConnectivity(const int *newNodeNumbers) +{ + int *conn=getNodalConnectivity()->getPointer(); + const int *connIndex=getNodalConnectivityIndex()->getConstPointer(); + int nbOfCells=getNumberOfCells(); + for(int i=0;i=0)//avoid polyhedron separator + { + node=newNodeNumbers[node]; + } + } +} + /*! * 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 @@ -286,16 +323,18 @@ void MEDCouplingUMesh::giveElemsInBoundingBox(const double *bbox, double eps, st for (int inode=conn_index[ielem]+1; inode=0)//avoid polyhedron separator { - if ( coords[node*dim+idim] < elem_bb[idim*2] ) - { - elem_bb[idim*2] = coords[node*dim+idim] ; - } - if ( coords[node*dim+idim] > elem_bb[idim*2+1] ) + for (int idim=0; idim elem_bb[idim*2+1] ) + { + elem_bb[idim*2+1] = coords[node*dim+idim] ; + } } } } @@ -532,10 +571,10 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField() const int N2 = connec[ipt+2]; int N3 = connec[ipt+3]; - area_vol[iel]=INTERP_KERNEL::calculateAreaForTria(coords+(dim_space*N1), - coords+(dim_space*N2), - coords+(dim_space*N3), - dim_space); + area_vol[iel]=fabs(INTERP_KERNEL::calculateAreaForTria(coords+(dim_space*N1), + coords+(dim_space*N2), + coords+(dim_space*N3), + dim_space)); } break ; @@ -547,11 +586,11 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField() const int N3 = connec[ipt+3]; int N4 = connec[ipt+4]; - area_vol[iel]=INTERP_KERNEL::calculateAreaForQuad(coords+dim_space*N1, - coords+dim_space*N2, - coords+dim_space*N3, - coords+dim_space*N4, - dim_space) ; + area_vol[iel]=fabs(INTERP_KERNEL::calculateAreaForQuad(coords+dim_space*N1, + coords+dim_space*N2, + coords+dim_space*N3, + coords+dim_space*N4, + dim_space)); } break ; diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 1c065e2d1..f45a780a8 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -57,10 +57,11 @@ namespace ParaMEDMEM void unserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, const std::vector& littleStrings); //tools - void zipCoords(); DataArrayInt *zipCoordsTraducer(); void getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const; MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const; + MEDCouplingPointSet *buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const; + void renumberConnectivity(const int *newNodeNumbers); void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems); MEDCouplingFieldDouble *getMeasureField() const; private: diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx index 995e3f77b..16b07423d 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx @@ -20,6 +20,8 @@ #include "CellModel.hxx" #include "MEDCouplingMemArray.hxx" +#include + using namespace ParaMEDMEM; MEDCouplingUMeshDesc::MEDCouplingUMeshDesc():_mesh_dim(-1),_desc_connec(0),_desc_connec_index(0), @@ -219,12 +221,28 @@ MEDCouplingPointSet *MEDCouplingUMeshDesc::buildPartOfMySelf(const int *start, c return 0; } +MEDCouplingPointSet *MEDCouplingUMeshDesc::buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const +{ + //not implemented yet + return 0; +} + +void MEDCouplingUMeshDesc::renumberConnectivity(const int *newNodeNumbers) +{ +} + MEDCouplingFieldDouble *MEDCouplingUMeshDesc::getMeasureField() const { //not implemented yet. return 0; } +DataArrayInt *MEDCouplingUMeshDesc::zipCoordsTraducer() +{ + //not implemented yet. + return 0; +} + void MEDCouplingUMeshDesc::computeTypes() { if(_desc_connec && _desc_connec_index) @@ -237,3 +255,9 @@ void MEDCouplingUMeshDesc::computeTypes() _types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]); } } + +void MEDCouplingUMeshDesc::checkFullyDefined() const throw(INTERP_KERNEL::Exception) +{ + if(!_desc_connec || !_desc_connec_index || !_nodal_connec_face || !_nodal_connec_face_index || !_coords) + throw INTERP_KERNEL::Exception("full connectivity and coordinates not set in unstructured mesh."); +} diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx index 79641fe67..ceced5640 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx @@ -49,11 +49,15 @@ namespace ParaMEDMEM const std::vector& littleStrings); void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems); MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const; + MEDCouplingPointSet *buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const; + void renumberConnectivity(const int *newNodeNumbers); MEDCouplingFieldDouble *getMeasureField() const; + DataArrayInt *zipCoordsTraducer(); private: MEDCouplingUMeshDesc(); ~MEDCouplingUMeshDesc(); void computeTypes(); + void checkFullyDefined() const throw(INTERP_KERNEL::Exception); private: unsigned _mesh_dim; DataArrayInt *_desc_connec; diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx index 031acfd26..a3ef04485 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx @@ -352,6 +352,54 @@ void MEDCouplingBasicsTest::testBuildPartOfMySelf() mesh->decrRef(); } +void MEDCouplingBasicsTest::testBuildPartOfMySelfNode() +{ + MEDCouplingUMesh *mesh=build2DTargetMesh_1(); + const int tab1[2]={5,7}; + MEDCouplingPointSet *subMeshSimple=mesh->buildPartOfMySelfNode(tab1,tab1+2,true); + MEDCouplingUMesh *subMesh=dynamic_cast(subMeshSimple); + CPPUNIT_ASSERT(subMesh); + CPPUNIT_ASSERT_EQUAL(1,(int)subMesh->getAllTypes().size()); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_QUAD4,*subMesh->getAllTypes().begin()); + CPPUNIT_ASSERT_EQUAL(1,subMesh->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(5,subMesh->getNodalConnectivity()->getNbOfElems()); + CPPUNIT_ASSERT_EQUAL(2,subMesh->getNodalConnectivityIndex()->getNbOfElems()); + const int subConn[5]={4,7,8,5,4}; + const int subConnIndex[3]={0,5}; + CPPUNIT_ASSERT(std::equal(subConn,subConn+5,subMesh->getNodalConnectivity()->getPointer())); + CPPUNIT_ASSERT(std::equal(subConnIndex,subConnIndex+2,subMesh->getNodalConnectivityIndex()->getPointer())); + CPPUNIT_ASSERT(subMesh->getCoords()==mesh->getCoords()); + subMeshSimple->decrRef(); + // + subMeshSimple=mesh->buildPartOfMySelfNode(tab1,tab1+2,false); + subMesh=dynamic_cast(subMeshSimple); + CPPUNIT_ASSERT(subMesh); + CPPUNIT_ASSERT_EQUAL(2,(int)subMesh->getAllTypes().size()); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,*subMesh->getAllTypes().begin()); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_QUAD4,*(++subMesh->getAllTypes().begin())); + CPPUNIT_ASSERT_EQUAL(3,subMesh->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(14,subMesh->getNodalConnectivity()->getNbOfElems()); + CPPUNIT_ASSERT_EQUAL(4,subMesh->getNodalConnectivityIndex()->getNbOfElems()); + const int subConn2[14]={3,4,5,2,4,6,7,4,3,4,7,8,5,4}; + const int subConnIndex2[4]={0,4,9,14}; + CPPUNIT_ASSERT(std::equal(subConn2,subConn2+14,subMesh->getNodalConnectivity()->getPointer())); + CPPUNIT_ASSERT(std::equal(subConnIndex2,subConnIndex2+4,subMesh->getNodalConnectivityIndex()->getPointer())); + CPPUNIT_ASSERT(subMesh->getCoords()==mesh->getCoords()); + subMeshSimple->decrRef(); + //testing the case where length of tab2 is greater than max number of node per cell. + const int tab2[7]={0,3,2,1,4,5,6}; + subMeshSimple=mesh->buildPartOfMySelfNode(tab2,tab2+7,true); + subMesh=dynamic_cast(subMeshSimple); + CPPUNIT_ASSERT(subMesh); + CPPUNIT_ASSERT_EQUAL(2,(int)subMesh->getAllTypes().size()); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,*subMesh->getAllTypes().begin()); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_QUAD4,*(++subMesh->getAllTypes().begin())); + CPPUNIT_ASSERT_EQUAL(3,subMesh->getNumberOfCells()); + subMeshSimple->decrRef(); + // + mesh->decrRef(); +} + void MEDCouplingBasicsTest::testZipCoords() { MEDCouplingUMesh *mesh=build2DTargetMesh_1(); @@ -808,6 +856,203 @@ void MEDCouplingBasicsTest::test3DSurfInterpP1P0_1() targetMesh->decrRef(); } +void MEDCouplingBasicsTest::test3DSurfInterpP1P1_1() +{ + MEDCouplingUMesh *sourceMesh=build3DSurfSourceMesh_1(); + MEDCouplingUMesh *targetMesh=build3DSurfTargetMesh_2(); + // + MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper(sourceMesh); + MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper(targetMesh); + INTERP_KERNEL::Interpolation3DSurf 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*sqrt(2.),res[0][0],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.05416666666666665*sqrt(2.),res[1][0],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666666*sqrt(2.),res[1][1],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.08333333333333334*sqrt(2.),res[2][1],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.05416666666666665*sqrt(2.),res[3][0],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666668*sqrt(2.),res[3][2],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.1416666666666666*sqrt(2.),res[4][0],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02499999999999999*sqrt(2.),res[4][1],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02499999999999999*sqrt(2.),res[4][2],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09999999999999999*sqrt(2.),res[4][3],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666666*sqrt(2.),res[5][1],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09583333333333333*sqrt(2.),res[5][3],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.08333333333333333*sqrt(2.),res[6][2],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666667*sqrt(2.),res[7][2],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09583333333333331*sqrt(2.),res[7][3],1.e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.04166666666666668*sqrt(2.),res[8][3],1.e-12); + res.clear(); + } + // + sourceMesh->decrRef(); + targetMesh->decrRef(); +} + +void MEDCouplingBasicsTest::test3DSurfInterpP0P0_2() +{ + MEDCouplingUMesh *sourceMesh=build3DSurfSourceMesh_1(); + MEDCouplingUMesh *targetMesh=build3DSurfTargetMeshPerm_1(); + // + MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper(sourceMesh); + MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper(targetMesh); + INTERP_KERNEL::Interpolation3DSurf myInterpolator; + vector > res; + myInterpolator.setPrecision(1e-12); + myInterpolator.setIntersectionType(INTERP_KERNEL::Triangulation); + { + myInterpolator.setOrientation(2); + myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0"); + CPPUNIT_ASSERT_EQUAL(5,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[1][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[2][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[3][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[4][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[4][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.*sqrt(2.),sumAll(res),1e-12); + res.clear(); + } + { + myInterpolator.setOrientation(0); + myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0"); + CPPUNIT_ASSERT_EQUAL(5,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[1][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.125*sqrt(2.),res[2][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[3][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[4][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[4][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.75*sqrt(2.),sumAll(res),1e-12); + res.clear(); + } + { + myInterpolator.setOrientation(1); + myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0"); + CPPUNIT_ASSERT_EQUAL(5,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[1][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[3][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[4][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[4][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.875*sqrt(2.),sumAll(res),1e-12); + res.clear(); + } + { + myInterpolator.setOrientation(-1); + myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0"); + CPPUNIT_ASSERT_EQUAL(5,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.125*sqrt(2.),res[2][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.125*sqrt(2.),sumAll(res),1e-12); + res.clear(); + } + //clean up + sourceMesh->decrRef(); + targetMesh->decrRef(); +} + +/*! + * Test of precision option implemented by Fabien that represents distance of "barycenter" to the other cell. + */ +void MEDCouplingBasicsTest::test3DSurfInterpP0P0_3() +{ + INTERP_KERNEL::Interpolation3DSurf myInterpolator; + vector > res; + double vecTrans[3]={0.,0.,1.e-10}; + double vec[3]={0.,-1.,0.}; + double pt[3]={-0.3,-0.3,5.e-11}; + const int N=32; + const double deltaA=M_PI/N; + myInterpolator.setPrecision(1e-12); + myInterpolator.setIntersectionType(INTERP_KERNEL::Triangulation); + myInterpolator.setMaxDistance3DSurfIntersect(1e-9); + for(int i=0;irotate(pt,vec,i*deltaA); + MEDCouplingUMesh *targetMesh=build3DSurfSourceMesh_2(); + targetMesh->translate(vecTrans); + targetMesh->rotate(pt,vec,i*deltaA); + MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper(sourceMesh); + MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper(targetMesh); + myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0"); + CPPUNIT_ASSERT_EQUAL(2,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,sumAll(res),1e-12); + sourceMesh->decrRef(); + targetMesh->decrRef(); + } + // + myInterpolator.setMaxDistance3DSurfIntersect(1e-11); + for(int i=0;irotate(pt,vec,i*deltaA); + MEDCouplingUMesh *targetMesh=build3DSurfSourceMesh_2(); + targetMesh->translate(vecTrans); + targetMesh->rotate(pt,vec,i*deltaA); + MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper(sourceMesh); + MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper(targetMesh); + myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0"); + CPPUNIT_ASSERT_EQUAL(2,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,sumAll(res),1e-12); + sourceMesh->decrRef(); + targetMesh->decrRef(); + } + // + res.clear(); + myInterpolator.setMaxDistance3DSurfIntersect(-1.);//unactivate fabien lookup + MEDCouplingUMesh *sourceMesh=build3DSurfSourceMesh_2(); + MEDCouplingUMesh *targetMesh=build3DSurfSourceMesh_2(); + targetMesh->translate(vecTrans); + myInterpolator.setBoundingBoxAdjustment(1e-11); + MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper0(sourceMesh); + MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper0(targetMesh); + myInterpolator.interpolateMeshes(sourceWrapper0,targetWrapper0,res,"P0P0"); + CPPUNIT_ASSERT_EQUAL(2,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,sumAll(res),1e-12); + sourceMesh->decrRef(); + targetMesh->decrRef(); + // + res.clear(); + sourceMesh=build3DSurfSourceMesh_2(); + targetMesh=build3DSurfSourceMesh_2(); + targetMesh->translate(vecTrans); + myInterpolator.setBoundingBoxAdjustment(1e-9); + MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper1(sourceMesh); + MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper1(targetMesh); + myInterpolator.interpolateMeshes(sourceWrapper1,targetWrapper1,res,"P0P0"); + CPPUNIT_ASSERT_EQUAL(2,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,sumAll(res),1e-12); + sourceMesh->decrRef(); + targetMesh->decrRef(); + //keeping the same bbox adj == 1.e-11 but trying rotation + res.clear(); + sourceMesh=build3DSurfSourceMesh_2(); + sourceMesh->rotate(pt,vec,M_PI/4.); + targetMesh=build3DSurfSourceMesh_2(); + targetMesh->translate(vecTrans); + targetMesh->rotate(pt,vec,M_PI/4.); + myInterpolator.setBoundingBoxAdjustment(1e-11); + MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper2(sourceMesh); + MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper2(targetMesh); + myInterpolator.interpolateMeshes(sourceWrapper2,targetWrapper2,res,"P0P0"); + CPPUNIT_ASSERT_EQUAL(2,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,sumAll(res),1e-12); + sourceMesh->decrRef(); + targetMesh->decrRef(); +} + void MEDCouplingBasicsTest::test3DInterpP0P0_1() { MEDCouplingUMesh *sourceMesh=build3DSourceMesh_1(); @@ -1049,6 +1294,24 @@ MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfSourceMesh_1() return sourceMesh; } +MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfSourceMesh_2() +{ + double sourceCoords[12]={-0.3,-0.3,0., 0.7,-0.3,0., -0.3,0.7,0., 0.7,0.7,0.}; + int sourceConn[6]={0,3,1,0,2,3}; + MEDCouplingUMesh *sourceMesh=MEDCouplingUMesh::New(); + sourceMesh->setMeshDimension(2); + sourceMesh->allocateCells(2); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,sourceConn); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,sourceConn+3); + sourceMesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(4,3); + std::copy(sourceCoords,sourceCoords+12,myCoords->getPointer()); + sourceMesh->setCoords(myCoords); + myCoords->decrRef(); + return sourceMesh; +} + MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfTargetMesh_1() { 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}; @@ -1070,6 +1333,30 @@ MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfTargetMesh_1() return targetMesh; } +/*! + * Idem build3DSurfTargetMesh_1 except that cell id 2 is not correctly numbered. + */ +MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfTargetMeshPerm_1() +{ + 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[18]={0,3,4,1, 1,4,2, 4,2,5, 6,7,4,3, 7,8,5,4}; + MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New(); + targetMesh->setMeshDimension(2); + targetMesh->allocateCells(5); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,targetConn); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+4); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+7); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,targetConn+10); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,targetConn+14); + 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::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}; diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx index 167e64e22..e3f38c6ec 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx @@ -38,6 +38,7 @@ namespace ParaMEDMEM CPPUNIT_TEST( testDeepCopy ); CPPUNIT_TEST( testRevNodal ); CPPUNIT_TEST( testBuildPartOfMySelf ); + CPPUNIT_TEST( testBuildPartOfMySelfNode ); CPPUNIT_TEST( testZipCoords ); CPPUNIT_TEST( testEqualMesh ); CPPUNIT_TEST( testEqualFieldDouble ); @@ -49,6 +50,9 @@ namespace ParaMEDMEM CPPUNIT_TEST( test3DSurfInterpP0P0_1 ); CPPUNIT_TEST( test3DSurfInterpP0P1_1 ); CPPUNIT_TEST( test3DSurfInterpP1P0_1 ); + CPPUNIT_TEST( test3DSurfInterpP1P1_1 ); + CPPUNIT_TEST( test3DSurfInterpP0P0_2 ); + CPPUNIT_TEST( test3DSurfInterpP0P0_3 ); CPPUNIT_TEST( test3DInterpP0P0_1 ); CPPUNIT_TEST( test3DInterpP0P1_1 ); CPPUNIT_TEST( test3DInterpP1P0_1 ); @@ -61,6 +65,7 @@ namespace ParaMEDMEM void testDeepCopy(); void testRevNodal(); void testBuildPartOfMySelf(); + void testBuildPartOfMySelfNode(); void testZipCoords(); void testEqualMesh(); void testEqualFieldDouble(); @@ -72,6 +77,9 @@ namespace ParaMEDMEM void test3DSurfInterpP0P0_1(); void test3DSurfInterpP0P1_1(); void test3DSurfInterpP1P0_1(); + void test3DSurfInterpP1P1_1(); + void test3DSurfInterpP0P0_2(); + void test3DSurfInterpP0P0_3(); void test3DInterpP0P0_1(); void test3DInterpP0P1_1(); void test3DInterpP1P0_1(); @@ -80,7 +88,9 @@ namespace ParaMEDMEM MEDCouplingUMesh *build2DTargetMesh_1(); MEDCouplingUMesh *build2DTargetMesh_2(); MEDCouplingUMesh *build3DSurfSourceMesh_1(); + MEDCouplingUMesh *build3DSurfSourceMesh_2(); MEDCouplingUMesh *build3DSurfTargetMesh_1(); + MEDCouplingUMesh *build3DSurfTargetMeshPerm_1(); MEDCouplingUMesh *build3DSurfTargetMesh_2(); MEDCouplingUMesh *build3DSourceMesh_1(); MEDCouplingUMesh *build3DTargetMesh_1(); diff --git a/src/ParaMEDMEM/DEC.cxx b/src/ParaMEDMEM/DEC.cxx index 27b9aa79d..b17616c6d 100644 --- a/src/ParaMEDMEM/DEC.cxx +++ b/src/ParaMEDMEM/DEC.cxx @@ -89,7 +89,14 @@ namespace ParaMEDMEM delete _local_field; delete _icoco_field; delete _union_group; - } + } + + void DEC::setNature(NatureOfField nature) + { + if(_local_field) + _local_field->getField()->setNature(nature); + } + /*! Attaches a local field to a DEC. If the processor is on the receiving end of the DEC, the field will be updated by a recvData() call. @@ -154,7 +161,7 @@ namespace ParaMEDMEM localgroup=_source_group; else localgroup=_target_group; - //delete _icoco_field; + //delete _icoco_field; _icoco_field=new ICoCo::MEDField(*const_cast(triofield), *localgroup); attachLocalField(_icoco_field); return; diff --git a/src/ParaMEDMEM/DEC.hxx b/src/ParaMEDMEM/DEC.hxx index 50b8dd1c2..0f4c545bb 100644 --- a/src/ParaMEDMEM/DEC.hxx +++ b/src/ParaMEDMEM/DEC.hxx @@ -37,7 +37,8 @@ namespace ParaMEDMEM { public: DEC():_local_field(0) { } - DEC(ProcessorGroup& source_group, ProcessorGroup& target_group); + DEC(ProcessorGroup& source_group, ProcessorGroup& target_group); + void setNature(NatureOfField nature); void attachLocalField( MEDCouplingFieldDouble* field); void attachLocalField(const ParaFIELD* field); void attachLocalField(const ICoCo::Field* field); diff --git a/src/ParaMEDMEM/ElementLocator.cxx b/src/ParaMEDMEM/ElementLocator.cxx index 2d243a4c2..773aa5086 100644 --- a/src/ParaMEDMEM/ElementLocator.cxx +++ b/src/ParaMEDMEM/ElementLocator.cxx @@ -21,9 +21,11 @@ #include "ElementLocator.hxx" #include "Topology.hxx" #include "BlockTopology.hxx" +#include "ParaFIELD.hxx" #include "ParaMESH.hxx" #include "ProcessorGroup.hxx" #include "MPIProcessorGroup.hxx" +#include "MEDCouplingFieldDouble.hxx" #include #include @@ -33,12 +35,12 @@ using namespace std; namespace ParaMEDMEM { - ElementLocator::ElementLocator(const ParaMESH& sourceMesh, + ElementLocator::ElementLocator(const ParaFIELD& sourceField, const ProcessorGroup& distant_group) - : _local_para_mesh(sourceMesh), - _local_cell_mesh(sourceMesh.getCellMesh()), - _local_face_mesh(sourceMesh.getFaceMesh()), - _local_group(*sourceMesh.getBlockTopology()->getProcGroup()), + : _local_para_field(sourceField), + _local_cell_mesh(sourceField.getSupport()->getCellMesh()), + _local_face_mesh(sourceField.getSupport()->getFaceMesh()), + _local_group(*sourceField.getSupport()->getBlockTopology()->getProcGroup()), _distant_group(distant_group) { _union_group = _local_group.fuse(distant_group); @@ -51,6 +53,12 @@ namespace ParaMEDMEM delete [] _domain_bounding_boxes; } + const MPI_Comm *ElementLocator::getCommunicator() const + { + MPIProcessorGroup* group=static_cast (_union_group); + return group->getComm(); + } + // ========================================================================== // Procedure for exchanging mesh between a distant proc and a local processor // param idistantrank proc id on distant group @@ -74,6 +82,16 @@ namespace ParaMEDMEM vector elems; double* distant_bb = _domain_bounding_boxes+rank*2*dim; _local_cell_mesh->giveElemsInBoundingBox(distant_bb,getBoundingBoxAdjustment(),elems); + + DataArrayInt *distant_ids_send; + MEDCouplingPointSet *send_mesh = (MEDCouplingPointSet *)_local_para_field.getField()->buildSubMeshData(&elems[0],&elems[elems.size()],distant_ids_send); + _exchangeMesh(send_mesh, distant_mesh, idistantrank, distant_ids_send, distant_ids); + distant_ids_send->decrRef(); + + if(send_mesh) + send_mesh->decrRef(); +#if 0 + int* distant_ids_send=0; //send_mesh contains null pointer if elems is empty MEDCouplingPointSet* send_mesh = _local_cell_mesh->buildPartOfMySelf(&elems[0],&elems[elems.size()],false); // Constituting an array containing the ids of the elements that are @@ -81,21 +99,16 @@ namespace ParaMEDMEM // This array enables the correct redistribution of the data when the // interpolated field is transmitted to the target array - int* distant_ids_send=0; if (elems.size()>0) { distant_ids_send = new int[elems.size()]; - int index=0; - for (std::vector::const_iterator iter = elems.begin(); iter!= elems.end(); iter++) - { - distant_ids_send[index]=*iter; - index++; - } + std::copy(elems.begin(),elems.end(),distant_ids_send); } _exchangeMesh(send_mesh, distant_mesh, idistantrank, distant_ids_send, distant_ids); delete[] distant_ids_send; if(send_mesh) send_mesh->decrRef(); +#endif } void ElementLocator::exchangeMethod(const std::string& sourceMeth, int idistantrank, std::string& targetMeth) @@ -173,7 +186,7 @@ namespace ParaMEDMEM void ElementLocator::_exchangeMesh( MEDCouplingPointSet* local_mesh, MEDCouplingPointSet*& distant_mesh, int iproc_distant, - const int* distant_ids_send, + const DataArrayInt* distant_ids_send, int*& distant_ids_recv) { CommInterface comm_interface=_union_group->getCommInterface(); @@ -185,7 +198,7 @@ namespace ParaMEDMEM //Getting tiny info of local mesh to allow the distant proc to initialize and allocate //the transmitted mesh. local_mesh->getTinySerializationInformation(tinyInfoLocal,tinyInfoLocalS); - tinyInfoLocal.push_back(local_mesh->getNumberOfCells()); + tinyInfoLocal.push_back(distant_ids_send->getNumberOfTuples()); tinyInfoDistant.resize(tinyInfoLocal.size()); std::fill(tinyInfoDistant.begin(),tinyInfoDistant.end(),0); MPIProcessorGroup* group=static_cast (_union_group); @@ -229,7 +242,7 @@ namespace ParaMEDMEM else distant_mesh_tmp->decrRef(); distant_ids_recv=new int[tinyInfoDistant.back()]; - comm_interface.sendRecv((void *)distant_ids_send,tinyInfoLocal.back(), MPI_INT, + comm_interface.sendRecv((void *)distant_ids_send->getConstPointer(),tinyInfoLocal.back(), MPI_INT, iprocdistant_in_union, 1113, distant_ids_recv,tinyInfoDistant.back(), MPI_INT, iprocdistant_in_union,1113, diff --git a/src/ParaMEDMEM/ElementLocator.hxx b/src/ParaMEDMEM/ElementLocator.hxx index 82d84c234..38a49bd0c 100644 --- a/src/ParaMEDMEM/ElementLocator.hxx +++ b/src/ParaMEDMEM/ElementLocator.hxx @@ -21,29 +21,39 @@ #include "InterpolationOptions.hxx" +#include #include #include namespace ParaMEDMEM { - class ParaMESH; + class ParaFIELD; class ProcessorGroup; class ParaSUPPORT; class InterpolationMatrix; class MEDCouplingPointSet; + class DataArrayInt; class ElementLocator : public INTERP_KERNEL::InterpolationOptions { public: - ElementLocator(const ParaMESH& sourceMesh, const ProcessorGroup& distant_group); + ElementLocator(const ParaFIELD& sourceField, const ProcessorGroup& distant_group); virtual ~ElementLocator(); void exchangeMesh(int idistantrank, MEDCouplingPointSet*& target_mesh, int*& distant_ids); void exchangeMethod(const std::string& sourceMeth, int idistantrank, std::string& targetMeth); + const std::vector& getDistantProcIds() const { return _distant_proc_ids; } + const MPI_Comm *getCommunicator() const; private: - const ParaMESH& _local_para_mesh ; + void _computeBoundingBoxes(); + bool _intersectsBoundingBox(int irank); + void _exchangeMesh(MEDCouplingPointSet* local_mesh, MEDCouplingPointSet*& distant_mesh, + int iproc_distant, const DataArrayInt* distant_ids_send, + int*& distant_ids_recv); + private: + const ParaFIELD& _local_para_field ; MEDCouplingPointSet* _local_cell_mesh; MEDCouplingPointSet* _local_face_mesh; std::vector _distant_cell_meshes; @@ -53,12 +63,6 @@ namespace ParaMEDMEM const ProcessorGroup& _local_group; ProcessorGroup* _union_group; std::vector _distant_proc_ids; - - void _computeBoundingBoxes(); - bool _intersectsBoundingBox(int irank); - void _exchangeMesh(MEDCouplingPointSet* local_mesh, MEDCouplingPointSet*& distant_mesh, - int iproc_distant, const int* distant_ids_send, - int*& distant_ids_recv); }; } diff --git a/src/ParaMEDMEM/GlobalizerMesh.cxx b/src/ParaMEDMEM/GlobalizerMesh.cxx new file mode 100644 index 000000000..909215846 --- /dev/null +++ b/src/ParaMEDMEM/GlobalizerMesh.cxx @@ -0,0 +1,143 @@ +// 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 "GlobalizerMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" +#include "CommInterface.hxx" + +using namespace std; + +namespace ParaMEDMEM +{ + GlobalizerMesh::GlobalizerMesh(const MPI_Comm *comm, MEDCouplingFieldDouble *localField):_comm(comm),_local_field(localField) + { + if(_local_field) + _local_field->incrRef(); + } + + GlobalizerMesh::~GlobalizerMesh() + { + if(_local_field) + _local_field->decrRef(); + } + + NatureOfField GlobalizerMesh::getLocalNature() const + { + return _local_field->getNature(); + } + + GlobalizerMeshWorkingSide::GlobalizerMeshWorkingSide(const MPI_Comm *comm, MEDCouplingFieldDouble *localField, + const std::string& distantMeth, const std::vector& lazyProcs):GlobalizerMesh(comm,localField),_distant_method(distantMeth),_lazy_procs(lazyProcs) + { + } + + GlobalizerMeshWorkingSide::~GlobalizerMeshWorkingSide() + { + } + + const std::vector& GlobalizerMeshWorkingSide::getProcIdsInInteraction() const + { + return _lazy_procs; + } + + /*! + * connected with GlobalizerMeshLazySide::recvFromWorkingSide + */ + void GlobalizerMeshWorkingSide::sendSumToLazySide(const std::vector< std::vector >& distantLocEltIds, const std::vector< std::vector >& partialSumRelToDistantIds) + { + int procId=0; + CommInterface comm; + for(vector::const_iterator iter=_lazy_procs.begin();iter!=_lazy_procs.end();iter++,procId++) + { + const vector& eltIds=distantLocEltIds[procId]; + const vector& valued=partialSumRelToDistantIds[procId]; + int lgth=eltIds.size(); + comm.send(&lgth,1,MPI_INT,*iter,1114,*_comm); + comm.send((void *)&eltIds[0],lgth,MPI_INT,*iter,1115,*_comm); + comm.send((void *)&valued[0],lgth,MPI_DOUBLE,*iter,1116,*_comm); + } + } + + /*! + * connected with GlobalizerMeshLazySide::sendToWorkingSide + */ + void GlobalizerMeshWorkingSide::recvSumFromLazySide(std::vector< std::vector >& globalSumRelToDistantIds) + { + int procId=0; + CommInterface comm; + MPI_Status status; + for(vector::const_iterator iter=_lazy_procs.begin();iter!=_lazy_procs.end();iter++,procId++) + { + std::vector& vec=globalSumRelToDistantIds[procId]; + comm.recv(&vec[0],vec.size(),MPI_DOUBLE,*iter,1117,*_comm,&status); + } + } + + GlobalizerMeshLazySide::GlobalizerMeshLazySide(const MPI_Comm *comm, MEDCouplingFieldDouble *localField, const std::vector& computeProcs):GlobalizerMesh(comm,localField),_compute_procs(computeProcs) + { + } + + GlobalizerMeshLazySide::~GlobalizerMeshLazySide() + { + } + + /*! + * connected with GlobalizerMeshWorkingSide::sendSumToLazySide + */ + void GlobalizerMeshLazySide::recvFromWorkingSide() + { + _values_added.resize(_local_field->getNumberOfTuples()); + int procId=0; + CommInterface comm; + _ids_per_working_proc.resize(_compute_procs.size()); + MPI_Status status; + for(vector::const_iterator iter=_compute_procs.begin();iter!=_compute_procs.end();iter++,procId++) + { + int lgth; + comm.recv(&lgth,1,MPI_INT,*iter,1114,*_comm,&status); + vector& ids=_ids_per_working_proc[procId]; + ids.resize(lgth); + vector values(lgth); + comm.recv(&ids[0],lgth,MPI_INT,*iter,1115,*_comm,&status); + comm.recv(&values[0],lgth,MPI_DOUBLE,*iter,1116,*_comm,&status); + for(int i=0;i::const_iterator iter=_compute_procs.begin();iter!=_compute_procs.end();iter++,procId++) + { + vector& ids=_ids_per_working_proc[procId]; + vector valsToSend(ids.size()); + vector::iterator iter3=valsToSend.begin(); + for(vector::const_iterator iter2=ids.begin();iter2!=ids.end();iter2++,iter3++) + *iter3=_values_added[*iter2]; + comm.send(&valsToSend[0],ids.size(),MPI_DOUBLE,*iter,1117,*_comm); + ids.clear(); + } + _ids_per_working_proc.clear(); + } +} + diff --git a/src/ParaMEDMEM/GlobalizerMesh.hxx b/src/ParaMEDMEM/GlobalizerMesh.hxx new file mode 100644 index 000000000..3436545a8 --- /dev/null +++ b/src/ParaMEDMEM/GlobalizerMesh.hxx @@ -0,0 +1,72 @@ +// 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 __GLOBALIZERMESH_HXX__ +#define __GLOBALIZERMESH_HXX__ + +#include "MEDCouplingNatureOfField.hxx" + +#include +#include +#include + +namespace ParaMEDMEM +{ + class MEDCouplingFieldDouble; + + class GlobalizerMesh + { + protected: + GlobalizerMesh(const MPI_Comm *comm, MEDCouplingFieldDouble *localField); + public: + NatureOfField getLocalNature() const; + virtual ~GlobalizerMesh(); + protected: + const MPI_Comm *_comm; + MEDCouplingFieldDouble *_local_field; + }; + + class GlobalizerMeshWorkingSide : public GlobalizerMesh + { + public: + GlobalizerMeshWorkingSide(const MPI_Comm *comm, MEDCouplingFieldDouble *localField, const std::string& distantMeth, const std::vector& lazyProcs); + ~GlobalizerMeshWorkingSide(); + // + const std::vector& getProcIdsInInteraction() const; + void sendSumToLazySide(const std::vector< std::vector >& distantLocEltIds, const std::vector< std::vector >& partialSumRelToDistantIds); + void recvSumFromLazySide(std::vector< std::vector >& globalSumRelToDistantIds); + private: + std::string _distant_method; + std::vector _lazy_procs; + }; + + class GlobalizerMeshLazySide : public GlobalizerMesh + { + public: + GlobalizerMeshLazySide(const MPI_Comm *comm, MEDCouplingFieldDouble *localField, const std::vector& computeProcs); + ~GlobalizerMeshLazySide(); + void recvFromWorkingSide(); + void sendToWorkingSide(); + private: + std::vector _compute_procs; + std::vector _values_added; + std::vector< std::vector > _ids_per_working_proc; + }; +} + +#endif diff --git a/src/ParaMEDMEM/ICoCoMEDField.cxx b/src/ParaMEDMEM/ICoCoMEDField.cxx index 9ba7ed5eb..f2e687091 100644 --- a/src/ParaMEDMEM/ICoCoMEDField.cxx +++ b/src/ParaMEDMEM/ICoCoMEDField.cxx @@ -94,7 +94,7 @@ namespace ICoCo for(int j=0;jinsertNextCell(elemtype,triofield._nodes_per_elem,conn); } delete[] conn; diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx index 80e44c9cb..d75d5428c 100644 --- a/src/ParaMEDMEM/InterpolationMatrix.cxx +++ b/src/ParaMEDMEM/InterpolationMatrix.cxx @@ -17,6 +17,7 @@ // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // #include "ParaMESH.hxx" +#include "ParaFIELD.hxx" #include "ProcessorGroup.hxx" #include "MxN_Mapping.hxx" #include "InterpolationMatrix.hxx" @@ -29,6 +30,7 @@ #include "MEDCouplingNormalizedUnstructuredMesh.txx" #include "InterpolationOptions.hxx" #include "NormalizedUnstructuredMesh.hxx" +#include "GlobalizerMesh.hxx" // class InterpolationMatrix // This class enables the storage of an interpolation matrix Wij mapping @@ -51,28 +53,23 @@ namespace ParaMEDMEM // param method interpolation method // ==================================================================== - InterpolationMatrix::InterpolationMatrix(ParaMEDMEM::ParaMESH *source_support, + InterpolationMatrix::InterpolationMatrix(const ParaMEDMEM::ParaFIELD *source_field, const ProcessorGroup& source_group, const ProcessorGroup& target_group, const DECOptions& dec_options, const INTERP_KERNEL::InterpolationOptions& interp_options): - _source_support(source_support->getCellMesh()), + _source_field(source_field), + _source_support(source_field->getSupport()->getCellMesh()), _mapping(source_group, target_group, dec_options), _source_group(source_group), _target_group(target_group), - _source_volume(0), DECOptions(dec_options), INTERP_KERNEL::InterpolationOptions(interp_options) { - int nbelems = _source_support->getNumberOfCells(); - + int nbelems = source_field->getField()->getNumberOfTuples(); _row_offsets.resize(nbelems+1); - for (int i=0; i > surfaces; int colSize=0; @@ -137,7 +134,7 @@ namespace ParaMEDMEM source_wrapper.ReleaseTempArrays(); } else if ( distant_support.getMeshDimension() == 3 - && distant_support.getSpaceDimension() ==3 ) + && distant_support.getSpaceDimension() == 3 ) { MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC); MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC); @@ -153,22 +150,15 @@ namespace ParaMEDMEM } int source_size=surfaces.size(); + bool needTargetSurf=isSurfaceComputationNeeded(targetMeth); - MEDCouplingFieldDouble *target_triangle_surf = distant_support.getMeasureField(); - MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(); - - // Storing the source volumes - _source_volume.resize(source_size); - for (int i=0; igetIJ(i,0); - } - - source_triangle_surf->decrRef(); + MEDCouplingFieldDouble *target_triangle_surf; + if(needTargetSurf) + target_triangle_surf = distant_support.getMeasureField(); //loop over the elements to build the interpolation //matrix structures - for (int ielem=0; ielem < surfaces.size(); ielem++) + for (int ielem=0; ielem < source_size; ielem++) { _row_offsets[ielem+1] += surfaces[ielem].size(); //_source_indices.push_back(make_pair(iproc_distant, ielem)); @@ -178,7 +168,9 @@ namespace ParaMEDMEM iter++) { //surface of the target triangle - double surf = target_triangle_surf->getIJ(iter->first,0); + double surf; + if(needTargetSurf) + surf = target_triangle_surf->getIJ(iter->first,0); //locating the (iproc, itriangle) pair in the list of columns map,int >::iterator iter2 = _col_offsets.find(make_pair(iproc_distant,iter->first)); @@ -189,11 +181,11 @@ namespace ParaMEDMEM //(iproc, itriangle) is not registered in the list //of distant elements - col_id =_col_offsets.size()+1; + col_id =_col_offsets.size(); _col_offsets.insert(make_pair(make_pair(iproc_distant,iter->first),col_id)); _mapping.addElementFromSource(iproc_distant, distant_elems[iter->first]); - _target_volume.push_back(surf); + //target_volume.push_back(surf); } else { @@ -204,13 +196,208 @@ namespace ParaMEDMEM //ielem is the row, //col_id is the number of the column //iter->second is the value of the coefficient - + if(needTargetSurf) + _target_volume[ielem].push_back(surf); _coeffs[ielem].push_back(make_pair(col_id,iter->second)); } } - target_triangle_surf->decrRef(); + if(needTargetSurf) + target_triangle_surf->decrRef(); + } + + void InterpolationMatrix::finishContributionW(GlobalizerMeshWorkingSide& globalizer) + { + NatureOfField nature=globalizer.getLocalNature(); + switch(nature) + { + case ConservativeVolumic: + computeConservVolDenoW(globalizer); + break; + case Integral: + computeIntegralDenoW(globalizer); + break; + case IntegralGlobConstraint: + computeGlobConstraintDenoW(globalizer); + break; + default: + throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field."); + break; + } + /*for(map,int>::iterator iter=_col_offsets.begin();iter!=_col_offsets.end();iter++) + if((*iter).second==9) + cout << (*iter).first.first << " ** " << (*iter).first.second << endl; + cout << "--> " << _col_offsets[make_pair(4,74)] << endl; + for(vector > >::iterator iter3=_coeffs.begin();iter3!=_coeffs.end();iter3++) + for(vector >::iterator iter4=(*iter3).begin();iter4!=(*iter3).end();iter4++) + if((*iter4).first==569) + cout << " __ " << iter3-_coeffs.begin() << "___" << (*iter4).second << endl; + ostringstream st; st << "deno_" << _deno_multiply[0][0]; + ofstream fs(st.str().c_str()); + for(vector >::iterator iter1=_deno_multiply.begin();iter1!=_deno_multiply.end();iter1++) + { + for(vector::iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++) + fs << *iter2 << " "; + fs << endl; + }*/ + } + + void InterpolationMatrix::finishContributionL(GlobalizerMeshLazySide& globalizer) + { + NatureOfField nature=globalizer.getLocalNature(); + switch(nature) + { + case ConservativeVolumic: + computeConservVolDenoL(globalizer); + break; + case Integral: + computeIntegralDenoL(globalizer); + break; + case IntegralGlobConstraint: + computeGlobConstraintDenoL(globalizer); + break; + default: + throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field."); + break; + } + } + + void InterpolationMatrix::computeConservVolDenoW(GlobalizerMeshWorkingSide& globalizer) + { + computeGlobalRowSum(globalizer,_deno_multiply); + computeGlobalColSum(_deno_reverse_multiply); } + void InterpolationMatrix::computeConservVolDenoL(GlobalizerMeshLazySide& globalizer) + { + globalizer.recvFromWorkingSide(); + globalizer.sendToWorkingSide(); + } + + void InterpolationMatrix::computeIntegralDenoW(GlobalizerMeshWorkingSide& globalizer) + { + MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(); + _deno_multiply.resize(_coeffs.size()); + vector >::iterator iter6=_deno_multiply.begin(); + const double *values=source_triangle_surf->getArray()->getConstPointer(); + for(vector > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++) + { + (*iter6).resize((*iter4).size()); + std::fill((*iter6).begin(),(*iter6).end(),*values); + } + source_triangle_surf->decrRef(); + _deno_reverse_multiply=_target_volume; + } + + /*! + * Nothing to do because surface computation is on working side. + */ + void InterpolationMatrix::computeIntegralDenoL(GlobalizerMeshLazySide& globalizer) + { + } + + void InterpolationMatrix::computeGlobConstraintDenoW(GlobalizerMeshWorkingSide& globalizer) + { + computeGlobalColSum(_deno_multiply); + computeGlobalRowSum(globalizer,_deno_reverse_multiply); + } + + void InterpolationMatrix::computeGlobConstraintDenoL(GlobalizerMeshLazySide& globalizer) + { + globalizer.recvFromWorkingSide(); + globalizer.sendToWorkingSide(); + } + + void InterpolationMatrix::computeGlobalRowSum(GlobalizerMeshWorkingSide& globalizer, std::vector >& denoStrorage) + { + //stores id in distant procs sorted by lazy procs connected with + vector< vector > rowsPartialSumI; + //stores the corresponding values. + vector< vector > rowsPartialSumD; + computeLocalRowSum(globalizer.getProcIdsInInteraction(),rowsPartialSumI,rowsPartialSumD); + globalizer.sendSumToLazySide(rowsPartialSumI,rowsPartialSumD); + globalizer.recvSumFromLazySide(rowsPartialSumD); + divideByGlobalRowSum(globalizer.getProcIdsInInteraction(),rowsPartialSumI,rowsPartialSumD,denoStrorage); + } + + /*! + * @param distantProcs in parameter that indicates which lazy procs are concerned. + * @param resPerProcI out parameter that must be cleared before calling this method. The size of 1st dimension is equal to the size of 'distantProcs'. + * It contains the element ids (2nd dimension) of the corresponding lazy proc. + * @param resPerProcD out parameter with the same format than 'resPerProcI'. It contains corresponding sum values. + */ + void InterpolationMatrix::computeLocalRowSum(const std::vector& distantProcs, std::vector >& resPerProcI, + std::vector >& resPerProcD) const + { + resPerProcI.resize(distantProcs.size()); + resPerProcD.resize(distantProcs.size()); + std::vector res(_col_offsets.size()); + for(vector > >::const_iterator iter=_coeffs.begin();iter!=_coeffs.end();iter++) + for(vector >::const_iterator iter3=(*iter).begin();iter3!=(*iter).end();iter3++) + res[(*iter3).first]+=(*iter3).second; + set procsSet; + int id; + const vector >& mapping=_mapping.getSendingIds(); + for(vector >::const_iterator iter2=mapping.begin();iter2!=mapping.end();iter2++) + { + std::pair::iterator,bool> isIns=procsSet.insert((*iter2).first); + if(isIns.second) + id=std::find(distantProcs.begin(),distantProcs.end(),(*iter2).first)-distantProcs.begin(); + resPerProcI[id].push_back((*iter2).second); + resPerProcD[id].push_back(res[iter2-mapping.begin()]); + } + /* + for(map, int >::const_iterator iter2=_col_offsets.begin();iter2!=_col_offsets.end();iter2++) + { + std::pair::iterator,bool> isIns=procsSet.insert((*iter2).first.first); + if(isIns.second) + id=std::find(distantProcs.begin(),distantProcs.end(),(*iter2).first.first)-distantProcs.begin(); + resPerProcI[id].push_back((*iter2).first.second); + resPerProcD[id].push_back(res[(*iter2).second]); + }*/ + } + + void InterpolationMatrix::divideByGlobalRowSum(const std::vector& distantProcs, const std::vector >& resPerProcI, + const std::vector >& resPerProcD, std::vector >& deno) + { + map fastSums; + int procId=0; + const vector >& mapping=_mapping.getSendingIds(); + map< int, map > distIdToLocId; + for(map< pair,int >::iterator iter8=_col_offsets.begin();iter8!=_col_offsets.end();iter8++) + distIdToLocId[(*iter8).first.first][mapping[(*iter8).second].second]=(*iter8).first.second; + + for(vector::const_iterator iter1=distantProcs.begin();iter1!=distantProcs.end();iter1++,procId++) + { + const std::vector& currentProcI=resPerProcI[procId]; + const std::vector& currentProcD=resPerProcD[procId]; + vector::const_iterator iter3=currentProcD.begin(); + for(vector::const_iterator iter2=currentProcI.begin();iter2!=currentProcI.end();iter2++,iter3++) + fastSums[_col_offsets[std::make_pair(*iter1,distIdToLocId[*iter1][*iter2])]]=*iter3; + } + deno.resize(_coeffs.size()); + vector >::iterator iter6=deno.begin(); + for(vector > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++) + { + (*iter6).resize((*iter4).size()); + vector::iterator iter7=(*iter6).begin(); + for(vector >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++,iter7++) + *iter7=fastSums[(*iter5).first]; + } + } + + void InterpolationMatrix::computeGlobalColSum(std::vector >& denoStrorage) + { + denoStrorage.resize(_coeffs.size()); + vector >::iterator iter2=denoStrorage.begin(); + for(vector > >::const_iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter2++) + { + (*iter2).resize((*iter1).size()); + double sumOfCurrentRow=0.; + for(vector >::const_iterator iter3=(*iter1).begin();iter3!=(*iter1).end();iter3++) + sumOfCurrentRow+=(*iter3).second; + std::fill((*iter2).begin(),(*iter2).end(),sumOfCurrentRow); + } + } // ================================================================== // The call to this method updates the arrays on the target side @@ -222,11 +409,11 @@ namespace ParaMEDMEM void InterpolationMatrix::prepare() { - int nbelems = _source_support->getNumberOfCells(); + int nbelems = _source_field->getField()->getNumberOfTuples(); for (int ielem=0; ielem < nbelems; ielem++) { _row_offsets[ielem+1]+=_row_offsets[ielem]; - } + } _mapping.prepareSendRecv(); } @@ -251,7 +438,7 @@ namespace ParaMEDMEM //computing the matrix multiply on source side if (_source_group.containsMyRank()) { - int nbrows = _source_support->getNumberOfCells(); + int nbrows = _coeffs.size(); // performing W.S // W is the intersection matrix @@ -266,23 +453,11 @@ namespace ParaMEDMEM { int colid= _coeffs[irow][icol-_row_offsets[irow]].first; double value = _coeffs[irow][icol-_row_offsets[irow]].second; - target_value[(colid-1)*nbcomp+icomp]+=value*coeff_row; + double deno = _deno_multiply[irow][icol-_row_offsets[irow]]; + target_value[colid*nbcomp+icomp]+=value*coeff_row/deno; } } } - - // performing VT^(-1).(W.S) - // where VT^(-1) is the inverse of the diagonal matrix containing - // the volumes of target cells - - for (int i=0; i<_col_offsets.size();i++) - { - for (int icomp=0; icompgetNumberOfComponents(); vector source_value(_col_offsets.size()* nbcomp,0.0); _mapping.reverseSendRecv(&source_value[0],field); @@ -325,7 +499,7 @@ namespace ParaMEDMEM //treatment of the transpose matrix multiply on the source side if (_source_group.containsMyRank()) { - int nbrows = _source_support->getNumberOfCells(); + int nbrows = _coeffs.size(); double *array = field.getArray()->getPointer() ; // Initialization @@ -340,27 +514,19 @@ namespace ParaMEDMEM { int colid = _coeffs[irow][icol-_row_offsets[irow]].first; double value = _coeffs[irow][icol-_row_offsets[irow]].second; - + double deno = _deno_reverse_multiply[irow][icol-_row_offsets[irow]]; for (int icomp=0; icomp& res) const; + void computeLocalRowSum(const std::vector& distantProcs, std::vector >& resPerProcI, + std::vector >& resPerProcD) const; + void computeGlobalRowSum(GlobalizerMeshWorkingSide& globalizer, std::vector >& denoStrorage); + void computeGlobalColSum(std::vector >& denoStrorage); + void divideByGlobalRowSum(const std::vector& distantProcs, const std::vector >& resPerProcI, + const std::vector >& resPerProcD, std::vector >& deno); + private: + bool isSurfaceComputationNeeded(const std::string& method) const; + private: + const ParaMEDMEM::ParaFIELD *_source_field; std::vector _row_offsets; std::map, int > _col_offsets; - MEDCouplingPointSet *_source_support; + MEDCouplingPointSet *_source_support; MxN_Mapping _mapping; const ProcessorGroup& _source_group; const ProcessorGroup& _target_group; - std::vector _target_volume; - std::vector _source_volume; + std::vector< std::vector > _target_volume; std::vector > > _coeffs; + std::vector > _deno_multiply; + std::vector > _deno_reverse_multiply; }; } diff --git a/src/ParaMEDMEM/IntersectionDEC.cxx b/src/ParaMEDMEM/IntersectionDEC.cxx index 6cfa69ef2..8a97efdac 100644 --- a/src/ParaMEDMEM/IntersectionDEC.cxx +++ b/src/ParaMEDMEM/IntersectionDEC.cxx @@ -28,8 +28,7 @@ #include "InterpolationMatrix.hxx" #include "IntersectionDEC.hxx" #include "ElementLocator.hxx" - - +#include "GlobalizerMesh.hxx" namespace ParaMEDMEM { @@ -144,21 +143,21 @@ namespace ParaMEDMEM */ void IntersectionDEC::synchronize() { - ParaMEDMEM::ParaMESH* para_mesh = _local_field->getSupport(); - //cout <<"size of Interpolation Matrix"<containsMyRank()) { //locate the distant meshes - ElementLocator locator(*para_mesh, *_target_group); + ElementLocator locator(*_local_field, *_target_group); //transfering option from IntersectionDEC to ElementLocator locator.setBoundingBoxAdjustment(getBoundingBoxAdjustment()); MEDCouplingPointSet* distant_mesh=0; int* distant_ids=0; + std::string distantMeth; for (int i=0; i<_target_group->size(); i++) { // int idistant_proc = (i+_source_group->myRank())%_target_group->size(); @@ -166,26 +165,26 @@ namespace ParaMEDMEM //gathers pieces of the target meshes that can intersect the local mesh locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids); - std::string distantMeth; - locator.exchangeMethod(_method,idistant_proc,distantMeth); if (distant_mesh !=0) { + locator.exchangeMethod(_method,idistant_proc,distantMeth); //adds the contribution of the distant mesh on the local one int idistant_proc_in_union=_union_group->translateRank(_target_group,idistant_proc); std::cout <<"add contribution from proc "<myRank()<addContribution(*distant_mesh,idistant_proc_in_union,distant_ids,_method,distantMeth); - distant_mesh->decrRef(); delete[] distant_ids; distant_mesh=0; distant_ids=0; } - } + } + GlobalizerMeshWorkingSide globalizer(locator.getCommunicator(),_local_field->getField(),distantMeth,locator.getDistantProcIds()); + _interpolation_matrix->finishContributionW(globalizer); } if (_target_group->containsMyRank()) { - ElementLocator locator(*para_mesh, *_source_group); + ElementLocator locator(*_local_field, *_source_group); //transfering option from IntersectionDEC to ElementLocator locator.setBoundingBoxAdjustment(getBoundingBoxAdjustment()); @@ -198,16 +197,18 @@ namespace ParaMEDMEM //gathers pieces of the target meshes that can intersect the local mesh locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids); std::cout << " Data sent from "<<_union_group->myRank()<<" to source proc "<< idistant_proc<decrRef(); delete[] distant_ids; distant_mesh=0; distant_ids=0; } - } + } + GlobalizerMeshLazySide globalizer(locator.getCommunicator(),_local_field->getField(),locator.getDistantProcIds()); + _interpolation_matrix->finishContributionL(globalizer); } _interpolation_matrix->prepare(); } diff --git a/src/ParaMEDMEM/Makefile.am b/src/ParaMEDMEM/Makefile.am index b2d0e0bd5..f164a3c7e 100644 --- a/src/ParaMEDMEM/Makefile.am +++ b/src/ParaMEDMEM/Makefile.am @@ -54,6 +54,7 @@ InterpolationMatrix.hxx\ IntersectionDEC.hxx\ ExplicitCoincidentDEC.hxx\ ElementLocator.hxx\ +GlobalizerMesh.hxx\ ExplicitMapping.hxx\ ICoCoField.hxx \ ICoCoMEDField.hxx \ @@ -74,6 +75,7 @@ StructuredCoincidentDEC.cxx\ ExplicitCoincidentDEC.cxx\ IntersectionDEC.cxx\ ElementLocator.cxx\ +GlobalizerMesh.cxx\ MPIAccessDEC.cxx \ TimeInterpolator.cxx \ LinearTimeInterpolator.cxx\ diff --git a/src/ParaMEDMEM/MxN_Mapping.hxx b/src/ParaMEDMEM/MxN_Mapping.hxx index 9e019efe9..e87a05297 100644 --- a/src/ParaMEDMEM/MxN_Mapping.hxx +++ b/src/ParaMEDMEM/MxN_Mapping.hxx @@ -42,6 +42,9 @@ namespace ParaMEDMEM void sendRecv(double* field, MEDCouplingFieldDouble& field) const ; void reverseSendRecv(double* field, MEDCouplingFieldDouble& field) const ; + // + const std::vector >& getSendingIds() const { return _sending_ids; }//tmp + MPIAccessDEC* getAccessDEC(){ return _access_DEC; } private : ProcessorGroup* _union_group; diff --git a/src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx b/src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx index 35b7e8593..77d95ecfc 100644 --- a/src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx +++ b/src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx @@ -38,6 +38,8 @@ class ParaMEDMEMTest : public CppUnit::TestFixture CPPUNIT_TEST(testIntersectionDEC_2D); CPPUNIT_TEST(testIntersectionDEC_2DP0P1); CPPUNIT_TEST(testIntersectionDEC_3D); + CPPUNIT_TEST(testIntersectionDECNonOverlapp_2D_P0P0); + CPPUNIT_TEST(testIntersectionDECNonOverlapp_2D_P0P1P1P0); CPPUNIT_TEST(testSynchronousEqualIntersectionWithoutInterpNativeDEC_2D); CPPUNIT_TEST(testSynchronousEqualIntersectionWithoutInterpDEC_2D); @@ -81,6 +83,8 @@ public: void testIntersectionDEC_2D(); void testIntersectionDEC_2DP0P1(); void testIntersectionDEC_3D(); + void testIntersectionDECNonOverlapp_2D_P0P0(); + void testIntersectionDECNonOverlapp_2D_P0P1P1P0(); #ifdef MED_ENABLE_FVM void testNonCoincidentDEC_2D(); void testNonCoincidentDEC_3D(); diff --git a/src/ParaMEDMEM/Test/ParaMEDMEMTest_ICocoTrio.cxx b/src/ParaMEDMEM/Test/ParaMEDMEMTest_ICocoTrio.cxx index 27998113f..bb377556d 100644 --- a/src/ParaMEDMEM/Test/ParaMEDMEMTest_ICocoTrio.cxx +++ b/src/ParaMEDMEM/Test/ParaMEDMEMTest_ICocoTrio.cxx @@ -229,6 +229,7 @@ void ParaMEDMEMTest::testICocoTrio1() else dec_emetteur.attachLocalField((ICoCo::Field*) &champ_recepteur); + dec_emetteur.setNature(ConservativeVolumic); if(init) dec_emetteur.synchronize(); diff --git a/src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx b/src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx index 0ca5504d4..ffc018393 100644 --- a/src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx +++ b/src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx @@ -138,7 +138,10 @@ void ParaMEDMEMTest::testIntersectionDEC_2D_(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,NO_TIME,paramesh, comptopo); + { + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + parafield->getField()->setNature(ConservativeVolumic); + } else parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); int nb_local; @@ -171,7 +174,10 @@ 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,NO_TIME,paramesh, comptopo); + { + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + parafield->getField()->setNature(ConservativeVolumic); + } else parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); int nb_local; @@ -342,7 +348,10 @@ 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,NO_TIME,paramesh, comptopo); + { + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + parafield->getField()->setNature(ConservativeVolumic); + } else parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); int nb_local; @@ -375,7 +384,10 @@ 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,NO_TIME,paramesh, comptopo); + { + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + parafield->getField()->setNature(ConservativeVolumic); + } else parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); int nb_local; @@ -521,6 +533,420 @@ void ParaMEDMEMTest::testAsynchronousFastSourceIntersectionDEC_2D() testAsynchronousIntersectionDEC_2D(0.01,1,0.11,1,true,true,true,"P0","P0"); } +void ParaMEDMEMTest::testIntersectionDECNonOverlapp_2D_P0P0() +{ + // + const double sourceCoordsAll[2][8]={{0.4,0.5,0.4,1.5,1.6,1.5,1.6,0.5}, + {0.3,-0.5,1.6,-0.5,1.6,-1.5,0.3,-1.5}}; + const double targetCoordsAll[3][16]={{0.7,1.45,0.7,1.65,0.9,1.65,0.9,1.45, 1.1,1.4,1.1,1.6,1.3,1.6,1.3,1.4}, + {0.7,-0.6,0.7,0.7,0.9,0.7,0.9,-0.6, 1.1,-0.7,1.1,0.6,1.3,0.6,1.3,-0.7}, + {0.7,-1.55,0.7,-1.35,0.9,-1.35,0.9,-1.55, 1.1,-1.65,1.1,-1.45,1.3,-1.45,1.3,-1.65}}; + int conn4All[8]={0,1,2,3,4,5,6,7}; + double targetResults[3][2]={{34.,34.},{38.333333333333336,42.666666666666664},{47.,47.}}; + double targetResults2[3][2]={{0.28333333333333344,0.56666666666666687},{1.8564102564102569,2.0128205128205132},{1.0846153846153845,0.36153846153846159}}; + double targetResults3[3][2]={{3.7777777777777781,7.5555555555555562},{24.511111111111113,26.355555555555558},{14.1,4.7}}; + // + int size; + int rank; + MPI_Comm_size(MPI_COMM_WORLD,&size); + MPI_Comm_rank(MPI_COMM_WORLD,&rank); + // + if(size!=5) + return ; + int nproc_source = 2; + set self_procs; + set procs_source; + set procs_target; + + for (int i=0; icontainsMyRank()) + { + std::ostringstream stream; stream << "sourcemesh2D proc " << rank; + mesh=MEDCouplingUMesh::New(stream.str().c_str(),2); + mesh->allocateCells(2); + mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn4All); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(4,2); + const double *sourceCoords=sourceCoordsAll[rank]; + std::copy(sourceCoords,sourceCoords+8,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + paramesh=new ParaMESH(mesh,*source_group,"source mesh"); + ParaMEDMEM::ComponentTopology comptopo; + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + double *value=parafield->getField()->getArray()->getPointer(); + value[0]=34+13*((double)rank); + } + else + { + std::ostringstream stream; stream << "targetmesh2D proc " << rank-nproc_source; + mesh=MEDCouplingUMesh::New(stream.str().c_str(),2); + mesh->allocateCells(2); + mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn4All); + mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn4All+4); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(8,2); + const double *targetCoords=targetCoordsAll[rank-nproc_source]; + std::copy(targetCoords,targetCoords+16,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + paramesh=new ParaMESH (mesh,*target_group,"target mesh"); + ParaMEDMEM::ComponentTopology comptopo; + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + } + //test 1 - Conservative volumic + ParaMEDMEM::IntersectionDEC dec(*source_group,*target_group); + parafield->getField()->setNature(ConservativeVolumic); + if (source_group->containsMyRank()) + { + dec.setMethod("P0"); + dec.attachLocalField(parafield); + dec.synchronize(); + dec.setForcedRenormalization(false); + dec.sendData(); + } + else + { + dec.setMethod("P0"); + dec.attachLocalField(parafield); + dec.synchronize(); + dec.setForcedRenormalization(false); + dec.recvData(); + const double *res=parafield->getField()->getArray()->getConstPointer(); + const double *expected=targetResults[rank-nproc_source]; + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[0],res[0],1e-13); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[1],res[1],1e-13); + } + //test 2 - Integral + ParaMEDMEM::IntersectionDEC dec2(*source_group,*target_group); + parafield->getField()->setNature(Integral); + if (source_group->containsMyRank()) + { + dec2.setMethod("P0"); + dec2.attachLocalField(parafield); + dec2.synchronize(); + dec2.setForcedRenormalization(false); + dec2.sendData(); + } + else + { + dec2.setMethod("P0"); + dec2.attachLocalField(parafield); + dec2.synchronize(); + dec2.setForcedRenormalization(false); + dec2.recvData(); + const double *res=parafield->getField()->getArray()->getConstPointer(); + const double *expected=targetResults2[rank-nproc_source]; + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[0],res[0],1e-13); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[1],res[1],1e-13); + } + //test 3 - Integral with global constraint + ParaMEDMEM::IntersectionDEC dec3(*source_group,*target_group); + parafield->getField()->setNature(IntegralGlobConstraint); + if (source_group->containsMyRank()) + { + dec3.setMethod("P0"); + dec3.attachLocalField(parafield); + dec3.synchronize(); + dec3.setForcedRenormalization(false); + dec3.sendData(); + } + else + { + dec3.setMethod("P0"); + dec3.attachLocalField(parafield); + dec3.synchronize(); + dec3.setForcedRenormalization(false); + dec3.recvData(); + const double *res=parafield->getField()->getArray()->getConstPointer(); + const double *expected=targetResults3[rank-nproc_source]; + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[0],res[0],1e-13); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[1],res[1],1e-13); + } + //test 4 - Conservative volumic reversed + ParaMEDMEM::IntersectionDEC dec4(*source_group,*target_group); + parafield->getField()->setNature(ConservativeVolumic); + if (source_group->containsMyRank()) + { + dec4.setMethod("P0"); + dec4.attachLocalField(parafield); + dec4.synchronize(); + dec4.setForcedRenormalization(false); + dec4.recvData(); + const double *res=parafield->getField()->getArray()->getConstPointer(); + CPPUNIT_ASSERT_EQUAL(1,parafield->getField()->getNumberOfTuples()); + const double expected[]={37.8518518518519,43.5333333333333}; + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[rank],res[0],1e-13); + } + else + { + dec4.setMethod("P0"); + dec4.attachLocalField(parafield); + dec4.synchronize(); + dec4.setForcedRenormalization(false); + double *res=parafield->getField()->getArray()->getPointer(); + const double *toSet=targetResults[rank-nproc_source]; + res[0]=toSet[0]; + res[1]=toSet[1]; + dec4.sendData(); + } + //test 5 - Integral reversed + ParaMEDMEM::IntersectionDEC dec5(*source_group,*target_group); + parafield->getField()->setNature(Integral); + if (source_group->containsMyRank()) + { + dec5.setMethod("P0"); + dec5.attachLocalField(parafield); + dec5.synchronize(); + dec5.setForcedRenormalization(false); + dec5.recvData(); + const double *res=parafield->getField()->getArray()->getConstPointer(); + CPPUNIT_ASSERT_EQUAL(1,parafield->getField()->getNumberOfTuples()); + const double expected[]={0.794600591715977,1.35631163708087}; + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[rank],res[0],1e-13); + } + else + { + dec5.setMethod("P0"); + dec5.attachLocalField(parafield); + dec5.synchronize(); + dec5.setForcedRenormalization(false); + double *res=parafield->getField()->getArray()->getPointer(); + const double *toSet=targetResults2[rank-nproc_source]; + res[0]=toSet[0]; + res[1]=toSet[1]; + dec5.sendData(); + } + //test 6 - Integral with global constraint reversed + ParaMEDMEM::IntersectionDEC dec6(*source_group,*target_group); + parafield->getField()->setNature(IntegralGlobConstraint); + if (source_group->containsMyRank()) + { + dec6.setMethod("P0"); + dec6.attachLocalField(parafield); + dec6.synchronize(); + dec6.setForcedRenormalization(false); + dec6.recvData(); + const double *res=parafield->getField()->getArray()->getConstPointer(); + CPPUNIT_ASSERT_EQUAL(1,parafield->getField()->getNumberOfTuples()); + const double expected[]={36.4592592592593,44.5407407407407}; + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[rank],res[0],1e-13); + } + else + { + dec6.setMethod("P0"); + dec6.attachLocalField(parafield); + dec6.synchronize(); + dec6.setForcedRenormalization(false); + double *res=parafield->getField()->getArray()->getPointer(); + const double *toSet=targetResults3[rank-nproc_source]; + res[0]=toSet[0]; + res[1]=toSet[1]; + dec6.sendData(); + } + // + delete parafield; + mesh->decrRef(); + delete paramesh; + delete self_group; + delete target_group; + delete source_group; + // + MPI_Barrier(MPI_COMM_WORLD); +} + +void ParaMEDMEMTest::testIntersectionDECNonOverlapp_2D_P0P1P1P0() +{ + int size; + int rank; + MPI_Comm_size(MPI_COMM_WORLD,&size); + MPI_Comm_rank(MPI_COMM_WORLD,&rank); + // + if(size!=5) + return ; + int nproc_source = 2; + set self_procs; + set procs_source; + set procs_target; + + for (int i=0; icontainsMyRank()) + { + if(rank==0) + { + double coords[6]={-0.3,-0.3, 0.7,0.7, 0.7,-0.3}; + int conn[3]={0,1,2}; + int globalNode[3]={1,2,0}; + mesh=MEDCouplingUMesh::New("Source mesh Proc0",2); + mesh->allocateCells(1); + mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(3,2); + std::copy(coords,coords+6,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + } + if(rank==1) + { + double coords[6]={-0.3,-0.3, -0.3,0.7, 0.7,0.7}; + int conn[3]={0,1,2}; + int globalNode[3]={1,3,2}; + mesh=MEDCouplingUMesh::New("Source mesh Proc1",2); + mesh->allocateCells(1); + mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(3,2); + std::copy(coords,coords+6,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + } + paramesh=new ParaMESH(mesh,*source_group,"source mesh"); + ParaMEDMEM::ComponentTopology comptopo; + parafieldP0 = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + parafieldP1 = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); + double *valueP0=parafieldP0->getField()->getArray()->getPointer(); + double *valueP1=parafieldP1->getField()->getArray()->getPointer(); + parafieldP0->getField()->setNature(ConservativeVolumic); + parafieldP1->getField()->setNature(ConservativeVolumic); + if(rank==0) + { + valueP0[0]=31.; + valueP1[0]=34.; valueP1[1]=77.; valueP1[2]=53.; + } + if(rank==1) + { + valueP0[0]=47.; + valueP1[0]=34.; valueP1[1]=57.; valueP1[2]=77.; + } + } + else + { + if(rank==2) + { + double coords[10]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2 }; + int conn[7]={0,3,4,1, 1,4,2}; + int globalNode[5]={4,3,0,2,1}; + mesh=MEDCouplingUMesh::New("Target mesh Proc2",2); + mesh->allocateCells(2); + mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn); + mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn+4); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(5,2); + std::copy(coords,coords+10,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + } + if(rank==3) + { + double coords[6]={0.2,0.2, 0.7,-0.3, 0.7,0.2}; + int conn[3]={0,2,1}; + int globalNode[3]={1,0,5}; + mesh=MEDCouplingUMesh::New("Target mesh Proc3",2); + mesh->allocateCells(1); + mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(3,2); + std::copy(coords,coords+6,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + } + if(rank==4) + { + double coords[12]={-0.3,0.2, -0.3,0.7, 0.2,0.7, 0.2,0.2, 0.7,0.7, 0.7,0.2}; + int conn[8]={0,1,2,3, 3,2,4,5}; + int globalNode[6]={2,6,7,1,8,5}; + mesh=MEDCouplingUMesh::New("Target mesh Proc4",2); + mesh->allocateCells(2); + mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn); + mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+4); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(6,2); + std::copy(coords,coords+12,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + } + ParaMEDMEM::ComponentTopology comptopo; + paramesh=new ParaMESH(mesh,*target_group,"target mesh"); + parafieldP0 = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + parafieldP1 = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); + parafieldP0->getField()->setNature(ConservativeVolumic); + parafieldP1->getField()->setNature(ConservativeVolumic); + } + // test 1 - P0 P1 + ParaMEDMEM::IntersectionDEC dec(*source_group,*target_group); + if (source_group->containsMyRank()) + { + dec.setMethod("P0"); + dec.attachLocalField(parafieldP0); + dec.synchronize(); + dec.setForcedRenormalization(false); + dec.sendData(); + } + else + { + dec.setMethod("P1"); + dec.attachLocalField(parafieldP1); + dec.synchronize(); + dec.setForcedRenormalization(false); + dec.recvData(); + /*const double *res=parafield->getField()->getArray()->getConstPointer(); + const double *expected=targetResults[rank-nproc_source]; + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[0],res[0],1e-13); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[1],res[1],1e-13);*/ + } + // + delete parafieldP0; + delete parafieldP1; + mesh->decrRef(); + delete paramesh; + delete self_group; + delete target_group; + delete source_group; + // + MPI_Barrier(MPI_COMM_WORLD); +} + /*! * Tests an asynchronous exchange between two codes * one sends data with dtA as an interval, the max time being tmaxA @@ -592,13 +1018,16 @@ void ParaMEDMEMTest::testAsynchronousIntersectionDEC_2D(double dtA, double tmaxA meshname<< "Mesh_2_"<< rank+1; mesh=MEDLoader::ReadUMeshFromFile(strstream.str().c_str(),meshname.str().c_str(),0); - + paramesh=new ParaMESH (mesh,*source_group,"source mesh"); // ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT( support,*source_group); ParaMEDMEM::ComponentTopology comptopo; if(srcM=="P0") - parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + { + parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + parafield->getField()->setNature(ConservativeVolumic);//InvertIntegral);//ConservativeVolumic); + } else parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); @@ -630,12 +1059,15 @@ void ParaMEDMEMTest::testAsynchronousIntersectionDEC_2D(double dtA, double tmaxA meshname<< "Mesh_3_"<getField()->setNature(ConservativeVolumic);//InvertIntegral);//ConservativeVolumic); + } else parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo); diff --git a/src/ParaMEDMEM_Swig/test_IntersectionDEC.py b/src/ParaMEDMEM_Swig/test_IntersectionDEC.py index 31957b346..581f7ac2a 100755 --- a/src/ParaMEDMEM_Swig/test_IntersectionDEC.py +++ b/src/ParaMEDMEM_Swig/test_IntersectionDEC.py @@ -63,6 +63,7 @@ class ParaMEDMEMBasicsTest(unittest.TestCase): paramesh=ParaMESH(mesh,source_group,"source mesh") comptopo = ComponentTopology() parafield = ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo) + parafield.getField().setNature(ConservativeVolumic) nb_local=mesh.getNumberOfCells() value = [1.0]*nb_local parafield.getField().setValues(value) @@ -77,6 +78,7 @@ class ParaMEDMEMBasicsTest(unittest.TestCase): paramesh=ParaMESH(mesh,target_group,"target mesh") comptopo = ComponentTopology() parafield = ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo) + parafield.getField().setNature(ConservativeVolumic) nb_local=mesh.getNumberOfCells() value = [0.0]*nb_local parafield.getField().setValues(value) -- 2.39.2