From: vbd Date: Fri, 17 Oct 2008 12:44:14 +0000 (+0000) Subject: modification of option management X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=79efb23e9287f9b6070414645f1f6af84e0d0ef0;p=tools%2Fmedcoupling.git modification of option management --- diff --git a/src/INTERP_KERNEL/BBTree.txx b/src/INTERP_KERNEL/BBTree.txx index 4c9a851b4..5e595560b 100644 --- a/src/INTERP_KERNEL/BBTree.txx +++ b/src/INTERP_KERNEL/BBTree.txx @@ -23,12 +23,30 @@ private: std::vector _elems; bool _terminal; int _nbelems; + double _epsilon; public: - - BBTree(double* bbs, int* elems, int level, int nbelems): - _left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems) + /*! + Constructor of the bounding box tree + \param bbs pointer to the [xmin1 xmax1 ymin1 ymax1 xmin2 xmax2 ...] array containing the bounding boxes that are to be indexed. + \param elems array to the indices of the elements contained in the BBTree + \param level level in the BBTree recursive structure + \param nbelems nb of elements in the BBTree + \param epsilon precision to which points are decided to be coincident + + Parameters \a elems and \a level are used only by BBTree itself for creating trees recursively. A typical use is therefore : +\code + int nbelems=... + double* bbs= new double[2*2*nbelems]; + // filling bbs ... + ... + BBTree<2> tree = new BBTree<2>(elems,0,0,nbelems,1e-12); +\endcode + */ + + BBTree(double* bbs, int* elems, int level, int nbelems, double epsilon=1E-12): + _left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems),_epsilon(epsilon) { if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL) { @@ -86,12 +104,11 @@ public: } - _max_left=max_left+1e-12; - _min_right=min_right-1e-12; - _left=new BBTree(bbs, &(new_elems_left[0]), level+1, new_elems_left.size()); - _right=new BBTree(bbs, &(new_elems_right[0]), level+1, new_elems_right.size()); + _max_left=max_left+_epsilon; + _min_right=min_right-_epsilon; + _left=new BBTree(bbs, &(new_elems_left[0]), level+1, new_elems_left.size(),_epsilon); + _right=new BBTree(bbs, &(new_elems_right[0]), level+1, new_elems_right.size(),_epsilon); - } @@ -121,7 +138,7 @@ public: bool intersects = true; for (int idim=0; idim-1e-12|| bb_ptr[idim*2+1]-bb[idim*2]<1e-12) + if (bb_ptr[idim*2]-bb[idim*2+1]>-_epsilon|| bb_ptr[idim*2+1]-bb[idim*2]<_epsilon) intersects=false; } if (intersects) @@ -166,7 +183,7 @@ public: bool intersects = true; for (int idim=0; idim1e-12|| bb_ptr[idim*2+1]-xx[idim]<-1e-12) + if (bb_ptr[idim*2]-xx[idim]>_epsilon|| bb_ptr[idim*2+1]-xx[idim]<-_epsilon) intersects=false; } if (intersects) diff --git a/src/INTERP_KERNEL/Interpolation.hxx b/src/INTERP_KERNEL/Interpolation.hxx index f5f7b05b1..0a14ba5a2 100644 --- a/src/INTERP_KERNEL/Interpolation.hxx +++ b/src/INTERP_KERNEL/Interpolation.hxx @@ -7,14 +7,24 @@ * * */ +#include "InterpolationOptions.hxx" + namespace INTERP_KERNEL { template - class Interpolation + class Interpolation : public InterpolationOptions { public: - Interpolation() { } - + Interpolation() { + // double InterpolationOptions::_precision=1e-12; +// int InterpolationOptions::_printLevel=0; +// InterpolationOptions::setIntersectionType(Triangulation); +// InterpolationOptions::setMedianPlane(0.5); +// InterpolationOptions::setDoRotate(true); +// InterpolationOptions::setBoundingBoxAdjustment(0.1); +// InterpolationOptions::setSplittingPolicy(GENERAL_48); + } + Interpolation(const InterpolationOptions& io) :InterpolationOptions(io){} //interpolation of two triangular meshes. template void interpolateMeshes(const MyMeshType& mesh1, const MyMeshType& mesh2, MatrixType& result) diff --git a/src/INTERP_KERNEL/Interpolation2D.hxx b/src/INTERP_KERNEL/Interpolation2D.hxx index ce1b7a98d..0749d83ee 100755 --- a/src/INTERP_KERNEL/Interpolation2D.hxx +++ b/src/INTERP_KERNEL/Interpolation2D.hxx @@ -9,6 +9,7 @@ namespace INTERP_KERNEL { public: Interpolation2D() { } + Interpolation2D(const InterpolationOptions& io):InterpolationPlanar(io){}; public: bool doRotate() const { return false; } double medianPlane() const { return 0.; } diff --git a/src/INTERP_KERNEL/Interpolation3D.hxx b/src/INTERP_KERNEL/Interpolation3D.hxx index 5584b898f..f08fdb560 100644 --- a/src/INTERP_KERNEL/Interpolation3D.hxx +++ b/src/INTERP_KERNEL/Interpolation3D.hxx @@ -4,13 +4,14 @@ #include "Interpolation.hxx" #include "NormalizedUnstructuredMesh.hxx" #include "IntersectorHexa.hxx" -#include "OptionManager.hxx" +#include "InterpolationOptions.hxx" namespace INTERP_KERNEL { - class Interpolation3D : public Interpolation,OptionManager + class Interpolation3D : public Interpolation { public: Interpolation3D(); + Interpolation3D(const InterpolationOptions& io); template void interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result); private: diff --git a/src/INTERP_KERNEL/Interpolation3D.txx b/src/INTERP_KERNEL/Interpolation3D.txx index a1654cb69..0ae133e35 100644 --- a/src/INTERP_KERNEL/Interpolation3D.txx +++ b/src/INTERP_KERNEL/Interpolation3D.txx @@ -7,7 +7,6 @@ #include "IntersectorTetra.txx" #include "IntersectorHexa.txx" #include "Log.hxx" -#include "OptionManager.hxx" /// If defined, use recursion to traverse the binary search tree, else use the BBTree class #define USE_RECURSIVE_BBOX_FILTER @@ -36,7 +35,9 @@ namespace INTERP_KERNEL */ Interpolation3D::Interpolation3D() { - registerOption(&_splitting_policy,"SplittingPolicy",GENERAL_48); + } + Interpolation3D::Interpolation3D(const InterpolationOptions& io):Interpolation(io) + { } /** @@ -164,9 +165,7 @@ namespace INTERP_KERNEL case 8: - SplittingPolicy sp; - getOption("SplittingPolicy",sp); - intersector = new IntersectorHexa(srcMesh, targetMesh, targetIdx,sp); + intersector = new IntersectorHexa(srcMesh, targetMesh, targetIdx,getSplittingPolicy()); break; default: diff --git a/src/INTERP_KERNEL/Interpolation3DSurf.txx b/src/INTERP_KERNEL/Interpolation3DSurf.txx index cd0b7189b..1f012f1c4 100644 --- a/src/INTERP_KERNEL/Interpolation3DSurf.txx +++ b/src/INTERP_KERNEL/Interpolation3DSurf.txx @@ -13,11 +13,13 @@ namespace INTERP_KERNEL ,_medianPlane(DFT_MEDIAN_PLANE) ,_surf3DAdjustmentEps(DFT_SURF3D_ADJ_EPS) { - registerOption(&_medianPlane,"MedianPlane",DFT_MEDIAN_PLANE); - registerOption(&_doRotate,"DoRotate",true); - registerOption(&_surf3DAdjustmentEps,"BoundingBoxAdjustment",DFT_SURF3D_ADJ_EPS); } + Interpolation3DSurf::Interpolation3DSurf(const InterpolationOptions& io):InterpolationPlanar(io) + { + } + + /** \brief Function used to set the options for the intersection calculation \details The following options can be modified: diff --git a/src/INTERP_KERNEL/InterpolationOptions.hxx b/src/INTERP_KERNEL/InterpolationOptions.hxx new file mode 100644 index 000000000..fc5a49b8a --- /dev/null +++ b/src/INTERP_KERNEL/InterpolationOptions.hxx @@ -0,0 +1,66 @@ +#ifndef _INTERPOLATIONOPTIONS_HXX_ +#define _INTERPOLATIONOPTIONS_HXX_ + + +namespace INTERP_KERNEL { + typedef enum {Triangulation, Convex, Geometric2D, Generic} 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{ + private : + int _printLevel ; + IntersectionType _intersectionType; + double _precision; + double _medianPlane ; + bool _doRotate ; + double _boundingBoxAdjustment ; + int _orientation ; + SplittingPolicy _splittingPolicy ; + + public: + InterpolationOptions() + {init(); + } + int getPrintLevel() const {return _printLevel;} + void setPrintLevel(int pl) {_printLevel=pl;} + + IntersectionType getIntersectionType() const{return InterpolationOptions::_intersectionType;} + void setIntersectionType(IntersectionType it) {InterpolationOptions::_intersectionType=it;} + + double getPrecision() const {return InterpolationOptions::_precision;} + void setPrecision(double p) {InterpolationOptions::_precision=p;} + + double getMedianPlane() {return InterpolationOptions::_medianPlane;} + void setMedianPlane(double mp) {InterpolationOptions::_medianPlane=mp;} + + bool getDoRotate() {return InterpolationOptions::_doRotate;} + void setDoRotate( bool dr) {InterpolationOptions::_doRotate = dr;} + + double getBoundingBoxAdjustment() {return InterpolationOptions::_boundingBoxAdjustment;} + void setBoundingBoxAdjustment(double bba) {InterpolationOptions::_boundingBoxAdjustment=bba;} + + int getOrientation() {return InterpolationOptions::_orientation;} + void setOrientation(int o) {InterpolationOptions::_orientation=o;} + + SplittingPolicy getSplittingPolicy() {return _splittingPolicy;} + void setSplittingPolicy(SplittingPolicy sp) {_splittingPolicy=sp;} + void init() + { + _printLevel=0; + _intersectionType=Triangulation; + _precision=1e-12;; + _medianPlane =0.5; + _doRotate =true; + _boundingBoxAdjustment =0.1; + _orientation =0; + _splittingPolicy =GENERAL_48; + } + }; + +} +#endif diff --git a/src/INTERP_KERNEL/InterpolationPlanar.hxx b/src/INTERP_KERNEL/InterpolationPlanar.hxx index 94223f375..f68c39332 100755 --- a/src/INTERP_KERNEL/InterpolationPlanar.hxx +++ b/src/INTERP_KERNEL/InterpolationPlanar.hxx @@ -4,29 +4,27 @@ #include "Interpolation.hxx" #include "PlanarIntersector.hxx" #include "NormalizedUnstructuredMesh.hxx" -#include "OptionManager.hxx" +#include "InterpolationOptions.hxx" namespace INTERP_KERNEL { - typedef enum {Triangulation, Convex, Geometric2D, Generic} IntersectionType; template - class InterpolationPlanar : public Interpolation< InterpolationPlanar >, public OptionManager + class InterpolationPlanar : public Interpolation< InterpolationPlanar > { private: - int _printLevel; - double _precision; + double _dimCaracteristic; - IntersectionType _intersectionType; static const double DEFAULT_PRECISION; - int _orientation; + public: InterpolationPlanar(); + InterpolationPlanar(const InterpolationOptions & io); // geometric precision, debug print level, coice of the median plane, intersection etc ... void setOptions(double precision, int printLevel, - IntersectionType intersectionType, int orientation); + IntersectionType intersectionType, int orientation=0); // Main function to interpolate triangular and quadratic meshes template diff --git a/src/INTERP_KERNEL/InterpolationPlanar.txx b/src/INTERP_KERNEL/InterpolationPlanar.txx index 4f3464794..bae20b1d4 100644 --- a/src/INTERP_KERNEL/InterpolationPlanar.txx +++ b/src/INTERP_KERNEL/InterpolationPlanar.txx @@ -2,7 +2,7 @@ #define __INTERPOLATIONPLANAR_TXX__ #include "InterpolationPlanar.hxx" -#include "OptionManager.hxx" +#include "InterpolationOptions.hxx" #include "PlanarIntersector.hxx" #include "PlanarIntersector.txx" #include "TriangulationIntersector.hxx" @@ -29,16 +29,18 @@ namespace INTERP_KERNEL * two local meshes in two dimensions. Meshes can contain mixed triangular and quadrangular elements. */ template - InterpolationPlanar::InterpolationPlanar():_printLevel(0),_precision(DEFAULT_PRECISION),_dimCaracteristic(1), - _intersectionType(Triangulation),_orientation(1) + InterpolationPlanar::InterpolationPlanar():_dimCaracteristic(1) + { - registerOption(&_precision, "Precision",DEFAULT_PRECISION); - registerOption(&_intersectionType, "IntersectionType",Triangulation); - registerOption(&_printLevel, "PrintLevel",0); - registerOption(&_orientation, "IntersectionSign",1); + } + template + InterpolationPlanar::InterpolationPlanar(const InterpolationOptions& io):Interpolation< InterpolationPlanar >(io),_dimCaracteristic(1) + + { } + /** * \brief Function used to set the options for the intersection calculation * \details The following options can be modified: @@ -55,10 +57,10 @@ namespace INTERP_KERNEL template void InterpolationPlanar::setOptions(double precision, int printLevel, IntersectionType intersectionType, int orientation) { - _precision=precision; - _printLevel=printLevel; - _intersectionType=intersectionType; - _orientation = orientation; + InterpolationOptions::setPrecision(precision); + InterpolationOptions::setPrintLevel(printLevel); + InterpolationOptions::setIntersectionType(intersectionType); + InterpolationOptions::setOrientation(orientation); } @@ -111,7 +113,7 @@ namespace INTERP_KERNEL double DimCaracteristic_P=diagonalP/nbMaille_P; _dimCaracteristic=std::min(DimCaracteristic_S, DimCaracteristic_P); - if (_printLevel>=1) + if (InterpolationOptions::getPrintLevel()>=1) { std::cout << " - Characteristic size of the source mesh : " << DimCaracteristic_S << std::endl; std::cout << " - Characteristic size of the target mesh: " << DimCaracteristic_P << std::endl; @@ -120,18 +122,29 @@ namespace INTERP_KERNEL PlanarIntersector* intersector; - switch (_intersectionType) + switch (InterpolationOptions::getIntersectionType()) { case Triangulation: - intersector=new TriangulationIntersector(myMesh_P,myMesh_S, _dimCaracteristic,_precision, - medianPlane(), _printLevel); + intersector=new TriangulationIntersector( + myMesh_P, + myMesh_S, + _dimCaracteristic, + InterpolationOptions::getPrecision(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrintLevel()); break; case Convex: - intersector=new ConvexIntersector(myMesh_P,myMesh_S, _dimCaracteristic, _precision, - medianPlane(), doRotate(), _printLevel); + intersector=new ConvexIntersector( + myMesh_P, + myMesh_S, + _dimCaracteristic, + InterpolationOptions::getPrecision(), + InterpolationOptions::getDoRotate(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrintLevel()); break; case Geometric2D: - intersector=new Geometric2DIntersector(myMesh_P, myMesh_S, _dimCaracteristic, _precision); + intersector=new Geometric2DIntersector(myMesh_P, myMesh_S, _dimCaracteristic, InterpolationOptions::getPrecision()); break; // case MEDMEM::Generic: //intersector=new GenericIntersector(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, @@ -187,7 +200,8 @@ namespace INTERP_KERNEL // -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 && (double)_orientation >=0 ) || ( surf < 0.0 && (double)_orientation <=0 )) + int orientation=InterpolationOptions::getOrientation(); + if (( surf > 0.0 && orientation >=0 ) || ( surf < 0.0 && orientation <=0 )) result[i_P].insert(std::make_pair(OTT::indFC(i_S),surf)); counter++; } @@ -199,7 +213,7 @@ namespace INTERP_KERNEL /* DEBUG prints */ /***********************************/ - if (_printLevel >=1) + if (InterpolationOptions::getPrintLevel() >=1) { long end_intersection=clock(); // if (_printLevel >=2) diff --git a/src/INTERP_KERNEL/IntersectorHexa.hxx b/src/INTERP_KERNEL/IntersectorHexa.hxx index b50e3552c..b2c92d8a5 100644 --- a/src/INTERP_KERNEL/IntersectorHexa.hxx +++ b/src/INTERP_KERNEL/IntersectorHexa.hxx @@ -4,15 +4,11 @@ #include "TargetIntersector.hxx" #include "IntersectorTetra.hxx" #include "NormalizedUnstructuredMesh.hxx" -#include "OptionManager.hxx" +#include "InterpolationOptions.hxx" namespace INTERP_KERNEL { - /// 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. - enum SplittingPolicy { PLANAR_FACE_5 = 5, PLANAR_FACE_6 = 6, GENERAL_24 = 24, GENERAL_48 = 48 }; + /** * \brief Class responsible for calculating intersection between a hexahedron target element and @@ -20,7 +16,7 @@ namespace INTERP_KERNEL * */ template - class IntersectorHexa : public TargetIntersector, public OptionManager + class IntersectorHexa : public TargetIntersector, public InterpolationOptions { public: diff --git a/src/INTERP_KERNEL/OptionManager.hxx b/src/INTERP_KERNEL/OptionManager.hxx deleted file mode 100644 index ed6609b03..000000000 --- a/src/INTERP_KERNEL/OptionManager.hxx +++ /dev/null @@ -1,151 +0,0 @@ -#ifndef _OPTION_MANAGER_HXX -#define _OPTION_MANAGER_HXX - -#include -#include - -#include "InterpolationUtils.hxx" - -namespace INTERP_KERNEL -{ - class OptionGeneric - { - //dummy method to have some polymorphism... - virtual void dummy()=0; - }; - - template < class T> class shared_ptr - { - private : - T* _ptr; - int* _count; - public : - shared_ptr(T* ptr) - { - _ptr=ptr; - _count==new int(1); - } - virtual ~shared_ptr() - { - *_count --; - if (*_count ==0 ) - { - delete _ptr; - delete _count; - } - } - shared_ptr(const shared_ptr& sp) - { - _count = sp._count; - _ptr=sp._ptr; - (*_count)++; - } - - T& operator*() - { - return *_ptr; - } - }; - - template < class T > - class Option : public OptionGeneric - { - private: - T* _var; - public: - Option( T defaut, T * var): - _var(var) - { - *_var= defaut; - } - void setValue(T value) - { - *_var = value; - } - T getValue() - { - return *_var; - } - - //dummy method to have some polymorphism... - void dummy(){} - }; - - class OptionManager - { - private: - std::map< std::string, OptionGeneric* > _optionList; - - public: - - OptionManager() - {} - - OptionManager(const OptionManager& om) - { - std::map< std::string, OptionGeneric*>::const_iterator mi; - for(mi = om._optionList.begin() ; mi != om._optionList.end() ; mi++) - { - _optionList.insert(make_pair(mi->first,mi->second)); - } - } - - ~OptionManager() - { - // std::map< std::string, OptionGeneric*>::iterator mi; - //for(mi = _optionList.begin() ; mi != _optionList.end() ; mi++) {if ((*mi).second!=0) delete (*mi).second; (*mi).second=0;} - } - - template void registerOption( T * var, const std::string& name, T defaut) - { - OptionGeneric * newoption = new Option( defaut, var); - _optionList.insert(make_pair(name, newoption)); - } - - template void getOption(const std::string& name, T& var) - { - if(_optionList.find(name) != _optionList.end()) - { - Option* option_ptr = dynamic_cast*>(_optionList.find(name)->second); - var= option_ptr->getValue(); - } - else throw Exception - ("Option is not listed, please register the option before using the getOption command"); - } - - template void setOption(const std::string& name, T value) - { - if(_optionList.find(name) != _optionList.end()) - { - Option* option_ptr = dynamic_cast*>(_optionList.find(name)->second); - if (option_ptr != 0 ) - option_ptr->setValue(value); - else throw Exception ("Error setOption: Option is registered with a different type"); - } - else throw Exception - ("Option is not listed, please register the option before using the setOption command"); - } - - void setOption(const std::string& name, int value) - { - if(_optionList.find(name) != _optionList.end()) - { - Option* option_double_ptr = dynamic_cast*> (_optionList.find(name)->second); - if (option_double_ptr!=0) - setOption(name,(double) value); - else - { - Option* option_ptr =dynamic_cast*>(_optionList.find(name)->second); - if (option_ptr != 0 ) - option_ptr->setValue(value); - else throw Exception ("Error setOption: Option is registered with a different type"); - } - } - else throw Exception - ("Option is not listed, please register the option before using the setOption command"); - } - - }; -} - -#endif diff --git a/src/INTERP_KERNEL/Remapper.cxx b/src/INTERP_KERNEL/Remapper.cxx index 51bdc13d7..e8ff6740d 100644 --- a/src/INTERP_KERNEL/Remapper.cxx +++ b/src/INTERP_KERNEL/Remapper.cxx @@ -8,6 +8,15 @@ namespace INTERP_KERNEL { + +// int InterpolationOptions::_printLevel=0; +// IntersectionType InterpolationOptions::_intersectionType=Triangulation; +// double InterpolationOptions::_precision=1e-12;; +// double InterpolationOptions::_medianPlane =0.5; +// bool InterpolationOptions::_doRotate =true; +// double InterpolationOptions::_boundingBoxAdjustment =0.1; +// int InterpolationOptions::_orientation =0; +// SplittingPolicy InterpolationOptions::_splittingPolicy =GENERAL_48; Remapper::Remapper():_matrix(0) { diff --git a/src/INTERP_KERNEL/Test/Makefile.am b/src/INTERP_KERNEL/Test/Makefile.am index 356dec9e1..b8ef375b8 100644 --- a/src/INTERP_KERNEL/Test/Makefile.am +++ b/src/INTERP_KERNEL/Test/Makefile.am @@ -54,10 +54,12 @@ dist_libInterpKernelTest_la_SOURCES= \ SingleElementPlanarTests.cxx \ PointLocatorTest.cxx \ MEDMeshMaker.cxx \ + InterpolationOptionsTest.cxx \ QuadraticPlanarInterpTest.cxx \ QuadraticPlanarInterpTest2.cxx \ QuadraticPlanarInterpTest3.cxx \ - QuadraticPlanarInterpTest4.cxx + QuadraticPlanarInterpTest4.cxx + libInterpKernelTest_la_CPPFLAGS= @CPPUNIT_INCLUDES@ $(MED2_INCLUDES) $(HDF5_INCLUDES) \ diff --git a/src/INTERP_KERNEL/Test/TestInterpKernel.cxx b/src/INTERP_KERNEL/Test/TestInterpKernel.cxx index 9e4a1b9e7..8858b7f27 100644 --- a/src/INTERP_KERNEL/Test/TestInterpKernel.cxx +++ b/src/INTERP_KERNEL/Test/TestInterpKernel.cxx @@ -30,7 +30,7 @@ #include "MultiElement2DTests.hxx" #include "SingleElementPlanarTests.hxx" #include "QuadraticPlanarInterpTest.hxx" - +#include "InterpolationOptionsTest.hxx" using namespace INTERP_TEST; //--- Registers the fixture into the 'registry' @@ -45,6 +45,7 @@ CPPUNIT_TEST_SUITE_REGISTRATION( PointLocatorTest); CPPUNIT_TEST_SUITE_REGISTRATION( MultiElement2DTests ); CPPUNIT_TEST_SUITE_REGISTRATION( SingleElementPlanarTests ); CPPUNIT_TEST_SUITE_REGISTRATION( QuadraticPlanarInterpTest ); +CPPUNIT_TEST_SUITE_REGISTRATION( InterpolationOptionsTest ); // --- generic Main program from KERNEL_SRC/src/Basics/Test diff --git a/src/ParaMEDMEM/DEC.cxx b/src/ParaMEDMEM/DEC.cxx index 84e694f9a..2888b815f 100644 --- a/src/ParaMEDMEM/DEC.cxx +++ b/src/ParaMEDMEM/DEC.cxx @@ -48,6 +48,13 @@ The following code excerpt shows how to set options for an object that inherits namespace ParaMEDMEM { + //intitializing DECOptions variables +// string DECOptions::_method="P0"; +// TimeInterpolationMethod DECOptions::_timeInterpolationMethod=WithoutTimeInterp; +// bool DECOptions::_asynchronous = false; +// bool DECOptions::_forcedRenormalization = false; + +// AllToAllMethod DECOptions::_allToAllMethod = Native; /*! \addtogroup dec @{ @@ -56,11 +63,9 @@ namespace ParaMEDMEM _source_group(&source_group), _target_group(&target_group), _owns_field(false), - _forced_renormalization_flag(false), _icoco_field(0) { _union_group = source_group.fuse(target_group); - registerOption(&_forced_renormalization_flag,"ForcedRenormalization",false); } DEC::~DEC() diff --git a/src/ParaMEDMEM/DEC.hxx b/src/ParaMEDMEM/DEC.hxx index 6c627a067..e621e0eb6 100644 --- a/src/ParaMEDMEM/DEC.hxx +++ b/src/ParaMEDMEM/DEC.hxx @@ -1,9 +1,9 @@ #ifndef DEC_HXX_ #define DEC_HXX_ -#include "OptionManager.hxx" #include #include +#include "DECOptions.hxx" namespace ICoCo { @@ -15,10 +15,10 @@ namespace ParaMEDMEM class ProcessorGroup; class ParaFIELD; class CommInterface; - class DEC : public INTERP_KERNEL::OptionManager + class DEC : public DECOptions { public: - DEC():_local_field(0),_forced_renormalization_flag(false){} + DEC():_local_field(0){} DEC(ProcessorGroup& local_group, ProcessorGroup& distant_group); void attachLocalField( MEDMEM::FIELD* field); void attachLocalField(const ParaFIELD* field); @@ -40,7 +40,6 @@ namespace ParaMEDMEM ProcessorGroup* _target_group; const CommInterface* _comm_interface; - bool _forced_renormalization_flag; bool _owns_field; private: ICoCo::Field* _icoco_field; diff --git a/src/ParaMEDMEM/DECOptions.hxx b/src/ParaMEDMEM/DECOptions.hxx new file mode 100644 index 000000000..75c725a0e --- /dev/null +++ b/src/ParaMEDMEM/DECOptions.hxx @@ -0,0 +1,52 @@ +#ifndef _DECOPTIONS_HXX_ +#define _DECOPTIONS_HXX_ + + +namespace ParaMEDMEM { + //Enum describing the allToAll method used in the communication pattern + typedef enum{Native,PointToPoint} AllToAllMethod; + typedef enum{WithoutTimeInterp,LinearTimeInterp} TimeInterpolationMethod; + + class DECOptions + { + private: + string _method; + bool _asynchronous; + TimeInterpolationMethod _timeInterpolationMethod; + AllToAllMethod _allToAllMethod; + bool _forcedRenormalization; + + public: + DECOptions():_method("P0"), + _timeInterpolationMethod(WithoutTimeInterp), + _asynchronous(false), + _forcedRenormalization(false), + _allToAllMethod(Native) + {} + DECOptions(const DECOptions& deco) + { + _method=deco._method; + _timeInterpolationMethod=deco._timeInterpolationMethod; + _asynchronous=deco._asynchronous; + _forcedRenormalization=deco._forcedRenormalization; + _allToAllMethod=deco._allToAllMethod; + } + string getMethod() const {return _method;} + void setMethod(string m) {_method=m;} + + TimeInterpolationMethod getTimeInterpolationMethod() const{return DECOptions::_timeInterpolationMethod;} + void setTimeInterpolationMethod(TimeInterpolationMethod it) {DECOptions::_timeInterpolationMethod=it;} + + bool getForcedRenormalization() const {return DECOptions::_forcedRenormalization;} + void setForcedRenormalization( bool dr) {DECOptions::_forcedRenormalization = dr;} + + bool getAsynchronous() const {return DECOptions::_asynchronous;} + void setAsynchronous( bool dr) {DECOptions::_asynchronous = dr;} + + AllToAllMethod getAllToAllMethod() const {return _allToAllMethod;} + void setAllToAllMethod(AllToAllMethod sp) {_allToAllMethod=sp;} + + }; + +} +#endif diff --git a/src/ParaMEDMEM/ElementLocator.cxx b/src/ParaMEDMEM/ElementLocator.cxx index c91753534..8f80feab3 100644 --- a/src/ParaMEDMEM/ElementLocator.cxx +++ b/src/ParaMEDMEM/ElementLocator.cxx @@ -25,8 +25,6 @@ _local_group(*mesh.getBlockTopology()->getProcGroup()), { _union_group = _local_group.fuse(distant_group); _computeBoundingBoxes(); -//JR registerOption(&_adjustment_eps,"BoundingBoxAdjustement",0.1); - registerOption(&_adjustment_eps,"BoundingBoxAdjustment",0.1); } ElementLocator::ElementLocator(const ParaSUPPORT& support, const ProcessorGroup& distant_group) @@ -34,8 +32,6 @@ ElementLocator::ElementLocator(const ParaSUPPORT& support, const ProcessorGroup& _distant_group(distant_group), _union_group(_local_group.fuse(distant_group)) { -// registerOption(&_adjustment_eps,"BoundingBoxAdjustement",0.1); - registerOption(&_adjustment_eps,"BoundingBoxAdjustment",0.1); throw ("Element Locator SUPPORT constructor not implemented yet"); } @@ -174,6 +170,7 @@ bool ElementLocator::_intersectsBoundingBox(double* bb1, double* bb2, int dim) { double bbtemp[2*dim]; double deltamax=0.0; + double adjustment_eps=getBoundingBoxAdjustment(); for (int i=0; i< dim; i++) { double delta = bb1[2*i+1]-bb1[2*i]; @@ -181,8 +178,8 @@ bool ElementLocator::_intersectsBoundingBox(double* bb1, double* bb2, int dim) } for (int i=0; i #include -#include "OptionManager.hxx" +#include "InterpolationOptions.hxx" namespace MEDMEM { @@ -17,7 +17,7 @@ class ParaSUPPORT; class InterpolationMatrix; - class ElementLocator: public INTERP_KERNEL::OptionManager + class ElementLocator: public INTERP_KERNEL::InterpolationOptions { public: ElementLocator(const ParaMESH& mesh, const ProcessorGroup& distant_group); @@ -33,7 +33,6 @@ private: const ProcessorGroup& _local_group; ProcessorGroup* _union_group; std::vector _distant_proc_ids; - double _adjustment_eps; void _computeBoundingBoxes(); bool _intersectsBoundingBox(int irank); diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx index 0bf7d7c0a..ae09da5c0 100644 --- a/src/ParaMEDMEM/InterpolationMatrix.cxx +++ b/src/ParaMEDMEM/InterpolationMatrix.cxx @@ -8,6 +8,7 @@ #include "Interpolation3DSurf.txx" #include "Interpolation3D.txx" #include "MEDNormalizedUnstructuredMesh.hxx" +#include "InterpolationOptions.hxx" /*! \class InterpolationMatrix @@ -29,308 +30,299 @@ namespace ParaMEDMEM \param target_group processor group containing the distant processors \param method interpolation method */ -InterpolationMatrix::InterpolationMatrix( -const ParaMEDMEM::ParaMESH& source_support, -const ProcessorGroup& source_group, -const ProcessorGroup& target_group, -const OptionManager& option_manager): - OptionManager(option_manager), - _source_support(*source_support.getMesh()), - _mapping(source_group, target_group), - _source_group(source_group), - _target_group(target_group), - _source_volume(0) - -{ - int nbelems= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); - _row_offsets.resize(nbelems+1); - for (int i=0; i > surfaces; - //computation of the intersection volumes between source and target elements - int source_size= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); - + _coeffs.resize(nbelems); + } - - if (distant_support.getMeshDimension()==2 && distant_support.getSpaceDimension()==3) + InterpolationMatrix::~InterpolationMatrix() { - MEDNormalizedUnstructuredMesh<3,2>target_wrapper(&distant_support); - MEDNormalizedUnstructuredMesh<3,2> source_wrapper(&_source_support); - INTERP_KERNEL::Interpolation3DSurf interpolator; - double precision, median_plane; - int print_level; - INTERP_KERNEL::IntersectionType inter_type; - bool do_rotate; - getOption("DoRotate",do_rotate); - getOption("Precision",precision); - getOption("MedianPlane",median_plane); - getOption("PrintLevel",print_level); - getOption("IntersectionType",inter_type); - interpolator.setOptions(precision,print_level,median_plane,inter_type,do_rotate); - interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces); } - else if (distant_support.getMeshDimension()==2 && distant_support.getSpaceDimension()==2) - { - MEDNormalizedUnstructuredMesh<2,2>target_wrapper(&distant_support); - MEDNormalizedUnstructuredMesh<2,2> source_wrapper(&_source_support); - INTERP_KERNEL::Interpolation2D interpolator; - interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces); - } - else if (distant_support.getMeshDimension()==3 && distant_support.getSpaceDimension()==3) - { - MEDNormalizedUnstructuredMesh<3,3>target_wrapper(&distant_support); - MEDNormalizedUnstructuredMesh<3,3> source_wrapper(&_source_support); + /*! + \brief Adds the contribution of a distant subdomain to the interpolation matrix. - INTERP_KERNEL::Interpolation3D interpolator; - interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces); - } - else + The method adds contribution to the interpolation matrix. + For each row of the matrix, elements are addded as + a (column, coeff) pair in the _coeffs array. This column number refers + to an element on the target side via the _col_offsets array. It is made of a series + of (iproc, ielem) pairs. + The number of elements per row is stored in the row_offsets array. + + \param distant_support local representation of the distant subdomain + \param iproc_distant id of the distant subdomain (in the distant group) + \param distant_elems mapping between the local representation of the subdomain and the actual elem ids on the distant subdomain + */ + void InterpolationMatrix::addContribution(MEDMEM::MESH& distant_support, int iproc_distant, int* distant_elems) { - throw MEDMEM::MEDEXCEPTION("no interpolator exists for these mesh and space dimensions "); - } + if (distant_support.getMeshDimension() != _source_support.getMeshDimension() || + distant_support.getMeshDimension() != _source_support.getMeshDimension() ) + throw MEDMEM::MEDEXCEPTION("local and distant meshes do not have the same space and mesh dimensions"); + + //creating the interpolator structure + vector > surfaces; + //computation of the intersection volumes between source and target elements + int source_size= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); + + + + if (distant_support.getMeshDimension()==2 && distant_support.getSpaceDimension()==3) + { + MEDNormalizedUnstructuredMesh<3,2>target_wrapper(&distant_support); + MEDNormalizedUnstructuredMesh<3,2> source_wrapper(&_source_support); + INTERP_KERNEL::Interpolation3DSurf interpolator(*this); + interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces); + } + else if (distant_support.getMeshDimension()==2 && distant_support.getSpaceDimension()==2) + { + MEDNormalizedUnstructuredMesh<2,2>target_wrapper(&distant_support); + MEDNormalizedUnstructuredMesh<2,2> source_wrapper(&_source_support); + + INTERP_KERNEL::Interpolation2D interpolator(*this); + interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces); + } + else if (distant_support.getMeshDimension()==3 && distant_support.getSpaceDimension()==3) + { + MEDNormalizedUnstructuredMesh<3,3>target_wrapper(&distant_support); + MEDNormalizedUnstructuredMesh<3,3> source_wrapper(&_source_support); + + INTERP_KERNEL::Interpolation3D interpolator(*this); + interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces); + } + else + { + throw MEDMEM::MEDEXCEPTION("no interpolator exists for these mesh and space dimensions "); + } - if (surfaces.size() != source_size) - { - cout<<"surfaces.size()="<* target_triangle_surf = getSupportVolumes(target_support); - MEDMEM::SUPPORT source_support (const_cast(&_source_support),"all cells", MED_EN::MED_CELL); - MEDMEM::FIELD* source_triangle_surf = getSupportVolumes(source_support); - - //storing the source volumes - _source_volume.resize(source_size); - for (int i=0; igetValueIJ(i+1,1); + //computing the vectors containing the source and target element volumes + MEDMEM::SUPPORT target_support(&distant_support,"all cells", MED_EN::MED_CELL); + MEDMEM::FIELD* target_triangle_surf = getSupportVolumes(target_support); + MEDMEM::SUPPORT source_support (const_cast(&_source_support),"all cells", MED_EN::MED_CELL); + MEDMEM::FIELD* source_triangle_surf = getSupportVolumes(source_support); + + //storing the source volumes + _source_volume.resize(source_size); + for (int i=0; igetValueIJ(i+1,1); - //loop over the elements to build the interpolation - //matrix structures - for (int ielem=0; ielem < surfaces.size(); ielem++) - { - _row_offsets[ielem+1]+=surfaces[ielem].size(); - // _source_indices.push_back(make_pair(iproc_distant, ielem)); - for (map::const_iterator iter=surfaces[ielem].begin(); - iter!=surfaces[ielem].end(); - iter++) - { - //surface of the target triangle - double surf = target_triangle_surf->getValueIJ(iter->first,1); - - //locating the (iproc, itriangle) pair in the list of columns - vector >::iterator iter2 = find (_col_offsets.begin(), _col_offsets.end(),make_pair(iproc_distant,iter->first)); - int col_id; + //loop over the elements to build the interpolation + //matrix structures + for (int ielem=0; ielem < surfaces.size(); ielem++) + { + _row_offsets[ielem+1]+=surfaces[ielem].size(); + // _source_indices.push_back(make_pair(iproc_distant, ielem)); + for (map::const_iterator iter=surfaces[ielem].begin(); + iter!=surfaces[ielem].end(); + iter++) + { + //surface of the target triangle + double surf = target_triangle_surf->getValueIJ(iter->first,1); + + //locating the (iproc, itriangle) pair in the list of columns + vector >::iterator iter2 = find (_col_offsets.begin(), _col_offsets.end(),make_pair(iproc_distant,iter->first)); + int col_id; - if (iter2==_col_offsets.end()) - { - //(iproc, itriangle) is not registered in the list - //of distant elements - - _col_offsets.push_back(make_pair(iproc_distant,iter->first)); - col_id =_col_offsets.size(); - _mapping.addElementFromSource(iproc_distant,distant_elems[iter->first-1]); - _target_volume.push_back(surf); - } - else - { - col_id=iter2-_col_offsets.begin()+1; - } + if (iter2==_col_offsets.end()) + { + //(iproc, itriangle) is not registered in the list + //of distant elements + + _col_offsets.push_back(make_pair(iproc_distant,iter->first)); + col_id =_col_offsets.size(); + _mapping.addElementFromSource(iproc_distant,distant_elems[iter->first-1]); + _target_volume.push_back(surf); + } + else + { + col_id=iter2-_col_offsets.begin()+1; + } - //the non zero coefficient is stored - //ielem is the row, - //col_id is the number of the column - //iter->second is the value of the coefficient + //the non zero coefficient is stored + //ielem is the row, + //col_id is the number of the column + //iter->second is the value of the coefficient - _coeffs[ielem].push_back(make_pair(col_id,iter->second)); + _coeffs[ielem].push_back(make_pair(col_id,iter->second)); - } - } + } + } delete source_triangle_surf; delete target_triangle_surf; -} + } -/*! The call to this method updates the arrays on the target side -so that they know which amount of data from which processor they -should expect. -That call makes actual interpolations via multiply method -available. -*/ + /*! The call to this method updates the arrays on the target side + so that they know which amount of data from which processor they + should expect. + That call makes actual interpolations via multiply method + available. + */ -void InterpolationMatrix::prepare() -{ - int nbelems = _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); - for (int ielem=0; ielem < nbelems; ielem++) + void InterpolationMatrix::prepare() { - _row_offsets[ielem+1]+=_row_offsets[ielem]; - } - _mapping.prepareSendRecv(); -} + int nbelems = _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); + for (int ielem=0; ielem < nbelems; ielem++) + { + _row_offsets[ielem+1]+=_row_offsets[ielem]; + } + _mapping.prepareSendRecv(); + } -/*! - \brief performs t=Ws, where t is the target field, s is the source field + /*! + \brief performs t=Ws, where t is the target field, s is the source field - The call to this method must be called both on the working side -and on the idle side. On the working side, the vector T=VT^(-1).(W.S) -is computed and sent. On the idle side, no computation is done, but the -result from the working side is received and the field is updated. + The call to this method must be called both on the working side + and on the idle side. On the working side, the vector T=VT^(-1).(W.S) + is computed and sent. On the idle side, no computation is done, but the + result from the working side is received and the field is updated. - \param field source field on processors involved on the source side, target field on processors on the target side - */ -void InterpolationMatrix::multiply(MEDMEM::FIELD& field) const -{ - vector target_value(_col_offsets.size()* field.getNumberOfComponents(),0.0); + \param field source field on processors involved on the source side, target field on processors on the target side + */ + void InterpolationMatrix::multiply(MEDMEM::FIELD& field) const + { + vector target_value(_col_offsets.size()* field.getNumberOfComponents(),0.0); - //computing the matrix multiply on source side - if (_source_group.containsMyRank()) - { - int nbcomp = field.getNumberOfComponents(); - int nbrows= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); + //computing the matrix multiply on source side + if (_source_group.containsMyRank()) + { + int nbcomp = field.getNumberOfComponents(); + int nbrows= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); - // performing W.S - // W is the intersection matrix - // S is the source vector + // performing W.S + // W is the intersection matrix + // S is the source vector - for (int irow=0; irowgetNumberOfElements(MED_EN::MED_ALL_ELEMENTS); + int nbcomp = field.getNumberOfComponents(); + double* value = const_cast (field.getValue()); + for (int i=0; i& field) const { - int nbelems=field.getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS); int nbcomp = field.getNumberOfComponents(); - double* value = const_cast (field.getValue()); - for (int i=0; i source_value(_col_offsets.size()* nbcomp,0.0); + _mapping.reverseSendRecv(&source_value[0],field); - The call to this method must be called both on the working side -and on the idle side. On the working side, the target vector T is received and the - vector S=VS^(-1).(WT.T) is computed to update the field. -On the idle side, no computation is done, but the field is sent. - - \param field source field on processors involved on the source side, target field on processors on the target side - */ -void InterpolationMatrix::transposeMultiply(MEDMEM::FIELD& field) const -{ - int nbcomp = field.getNumberOfComponents(); - vector source_value(_col_offsets.size()* nbcomp,0.0); - _mapping.reverseSendRecv(&source_value[0],field); - - //treatment of the transpose matrix multiply on the source side - if (_source_group.containsMyRank()) - { - int nbrows = _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); + //treatment of the transpose matrix multiply on the source side + if (_source_group.containsMyRank()) + { + int nbrows = _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); - vector target_value(nbrows*nbcomp,0.0); + vector target_value(nbrows*nbcomp,0.0); - //performing WT.T - //WT is W transpose - //T is the target vector - for (int irow=0; irow* InterpolationMatrix::getSupportVolumes(const MEDMEM::SUPPORT& support) - { - const MEDMEM::MESH* mesh=support.getMesh(); - int dim = mesh->getMeshDimension(); - switch (dim) - { - case 2: - return mesh->getArea(&support); - case 3: - return mesh->getVolume(&support); - default: - throw MEDMEM::MEDEXCEPTION("interpolation is not available for this dimension"); - } - } + \param field field on which cells the volumes are required + \return field containing the volumes + */ + MEDMEM::FIELD* InterpolationMatrix::getSupportVolumes(const MEDMEM::SUPPORT& support) + { + const MEDMEM::MESH* mesh=support.getMesh(); + int dim = mesh->getMeshDimension(); + switch (dim) + { + case 2: + return mesh->getArea(&support); + case 3: + return mesh->getVolume(&support); + default: + throw MEDMEM::MEDEXCEPTION("interpolation is not available for this dimension"); + } + } } diff --git a/src/ParaMEDMEM/InterpolationMatrix.hxx b/src/ParaMEDMEM/InterpolationMatrix.hxx index bc55d9a78..ad2610ed6 100644 --- a/src/ParaMEDMEM/InterpolationMatrix.hxx +++ b/src/ParaMEDMEM/InterpolationMatrix.hxx @@ -4,17 +4,21 @@ #include "MEDMEM_Field.hxx" #include "MPI_AccessDEC.hxx" #include "MxN_Mapping.hxx" -#include "OptionManager.hxx" +#include "InterpolationOptions.hxx" +#include "DECOptions.hxx" + namespace ParaMEDMEM { - class InterpolationMatrix : public INTERP_KERNEL::OptionManager + class InterpolationMatrix : public INTERP_KERNEL::InterpolationOptions, public DECOptions { public: InterpolationMatrix(const ParaMEDMEM::ParaMESH& source_support, - const ProcessorGroup& local_group, - const ProcessorGroup& distant_group, - const INTERP_KERNEL::OptionManager& option_manager); + const ProcessorGroup& local_group, + const ProcessorGroup& distant_group, + const DECOptions& dec_opt, + const InterpolationOptions& i_opt); + virtual ~InterpolationMatrix(); void addContribution(MEDMEM::MESH& distant_support, int iproc_distant, int* distant_elems); @@ -22,8 +26,6 @@ namespace ParaMEDMEM void transposeMultiply(MEDMEM::FIELD&)const; void prepare(); int getNbRows() const {return _row_offsets.size();} - void setAllToAllMethod(const AllToAllMethod& method) - { _mapping.setAllToAllMethod(method);} MPI_AccessDEC* getAccessDEC(){return _mapping.getAccessDEC();} private: @@ -35,7 +37,7 @@ namespace ParaMEDMEM vector > _col_offsets; const MEDMEM::MESH& _source_support; MxN_Mapping _mapping; - string _method; + const ProcessorGroup& _source_group; const ProcessorGroup& _target_group; vector _target_volume; diff --git a/src/ParaMEDMEM/IntersectionDEC.cxx b/src/ParaMEDMEM/IntersectionDEC.cxx index 3a88d4c54..27c106bdc 100644 --- a/src/ParaMEDMEM/IntersectionDEC.cxx +++ b/src/ParaMEDMEM/IntersectionDEC.cxx @@ -7,12 +7,11 @@ #include "MPIProcessorGroup.hxx" #include "ParaMESH.hxx" #include "ParaSUPPORT.hxx" -#include "OptionManager.hxx" #include "DEC.hxx" #include "InterpolationMatrix.hxx" #include "IntersectionDEC.hxx" #include "ElementLocator.hxx" -#include "OptionManager.hxx" + namespace ParaMEDMEM @@ -89,7 +88,7 @@ In the P0-P0 case, this matrix is a plain rectangular matrix with coefficients e @{ */ - IntersectionDEC::IntersectionDEC() +IntersectionDEC::IntersectionDEC() { } @@ -103,19 +102,9 @@ The constructor must be called synchronously on all processors of both processor */ IntersectionDEC::IntersectionDEC(ProcessorGroup& local_group, ProcessorGroup& distant_group): - DEC(local_group, distant_group),_interpolation_matrix(0),_asynchronous(false),_timeinterpolationmethod(WithoutTimeInterp),_allToAllMethod(Native) + DEC(local_group, distant_group),_interpolation_matrix(0) { - registerOption(&_method,string("Method"),string("P0")); - registerOption(&_bb_adjustment, "BoundingBoxAdjustment",0.1); - registerOption(&_intersection_type, "IntersectionType",INTERP_KERNEL::Triangulation); - registerOption(&_precision, "Precision", 1e-12); - registerOption(&_median_plane, "MedianPlane",0.5); - registerOption(&_do_rotate, "DoRotate",true); - registerOption(&_asynchronous, "Asynchronous",false); - registerOption(&_splitting_policy, "SplittingPolicy",INTERP_KERNEL::GENERAL_48); - registerOption(&_print_level,"PrintLevel",0); - registerOption(&_timeinterpolationmethod,"TimeInterpolation",WithoutTimeInterp); - registerOption(&_allToAllMethod,"AllToAllMethod",Native); + } IntersectionDEC::~IntersectionDEC() @@ -141,10 +130,7 @@ void IntersectionDEC::synchronize() { const ParaMEDMEM::ParaMESH* para_mesh = _local_field->getSupport()->getMesh(); //cout <<"size of Interpolation Matrix"<setAllToAllMethod(_allToAllMethod); - _interpolation_matrix->getAccessDEC()->Asynchronous( _asynchronous ) ; - _interpolation_matrix->getAccessDEC()->SetTimeInterpolator( _timeinterpolationmethod ) ; + _interpolation_matrix = new InterpolationMatrix (*para_mesh, *_source_group,*_target_group,*this,*this); //setting up the communication DEC on both sides if (_source_group->containsMyRank()) @@ -153,9 +139,7 @@ void IntersectionDEC::synchronize() ElementLocator locator(*para_mesh, *_target_group); //transfering option from IntersectionDEC to ElementLocator - double bb_adj; - getOption("BoundingBoxAdjustment",bb_adj); - locator.setOption("BoundingBoxAdjustment",bb_adj); + locator.setBoundingBoxAdjustment(getBoundingBoxAdjustment()); MESH* distant_mesh=0; int* distant_ids=0; @@ -187,10 +171,9 @@ void IntersectionDEC::synchronize() { ElementLocator locator(*para_mesh, *_source_group); //transfering option from IntersectionDEC to ElementLocator - double bb_adj; - INTERP_KERNEL::OptionManager::getOption("BoundingBoxAdjustment",bb_adj); - locator.setOption("BoundingBoxAdjustment",bb_adj); - + locator.setBoundingBoxAdjustment(getBoundingBoxAdjustment()); + + MESH* distant_mesh=0; int* distant_ids=0; for (int i=0; i<_source_group->size(); i++) @@ -226,7 +209,7 @@ void IntersectionDEC::recvData() else if (_target_group->containsMyRank()) { _interpolation_matrix->multiply(*_local_field->getField()); - if (_forced_renormalization_flag) + if (getForcedRenormalization()) renormalizeTargetField(); } @@ -256,7 +239,7 @@ void IntersectionDEC::sendData() { _interpolation_matrix->multiply(*_local_field->getField()); - if (_forced_renormalization_flag) + if (getForcedRenormalization()) renormalizeTargetField(); } diff --git a/src/ParaMEDMEM/IntersectionDEC.hxx b/src/ParaMEDMEM/IntersectionDEC.hxx index 09ec5067f..a15c2d8d8 100644 --- a/src/ParaMEDMEM/IntersectionDEC.hxx +++ b/src/ParaMEDMEM/IntersectionDEC.hxx @@ -6,12 +6,13 @@ #include "MxN_Mapping.hxx" #include "InterpolationPlanar.hxx" #include "IntersectorHexa.hxx" +#include "InterpolationOptions.hxx" namespace ParaMEDMEM { class InterpolationMatrix; - class IntersectionDEC:public DEC + class IntersectionDEC:public DEC, public INTERP_KERNEL::InterpolationOptions { public: IntersectionDEC(); @@ -41,25 +42,9 @@ namespace ParaMEDMEM //local element number containing the distant points const int* _distant_locations; - - //inerpolation method - string _method; - + InterpolationMatrix* _interpolation_matrix; - //options - double _bb_adjustment; - bool _asynchronous; - TimeInterpolationMethod _timeinterpolationmethod ; - AllToAllMethod _allToAllMethod; - INTERP_KERNEL::IntersectionType _intersection_type; - double _precision; - double _median_plane; - bool _do_rotate; - INTERP_KERNEL::SplittingPolicy _splitting_policy; - int _print_level; - - }; } diff --git a/src/ParaMEDMEM/MPI_AccessDEC.hxx b/src/ParaMEDMEM/MPI_AccessDEC.hxx index 30b4fd5b7..39fc10e37 100644 --- a/src/ParaMEDMEM/MPI_AccessDEC.hxx +++ b/src/ParaMEDMEM/MPI_AccessDEC.hxx @@ -10,7 +10,7 @@ namespace ParaMEDMEM { - typedef enum{WithoutTimeInterp,LinearTimeInterp} TimeInterpolationMethod; + // typedef enum{WithoutTimeInterp,LinearTimeInterp} TimeInterpolationMethod; class MPI_AccessDEC { diff --git a/src/ParaMEDMEM/MxN_Mapping.cxx b/src/ParaMEDMEM/MxN_Mapping.cxx index e2197dcd1..067aa9c19 100644 --- a/src/ParaMEDMEM/MxN_Mapping.cxx +++ b/src/ParaMEDMEM/MxN_Mapping.cxx @@ -1,15 +1,18 @@ -#include "CommInterface.hxx" +#include "CommInterface.hxx" #include "ProcessorGroup.hxx" #include "MPIProcessorGroup.hxx" +#include "MPI_AccessDEC.hxx" #include "MxN_Mapping.hxx" namespace ParaMEDMEM { -MxN_Mapping::MxN_Mapping(const ProcessorGroup& local_group, const ProcessorGroup& distant_group) -: _union_group(local_group.fuse(distant_group)) +MxN_Mapping::MxN_Mapping(const ProcessorGroup& local_group, const ProcessorGroup& distant_group,const DECOptions& dec_options) + : _union_group(local_group.fuse(distant_group)), + DECOptions(dec_options) { - _accessDEC = new MPI_AccessDEC(local_group,distant_group); + _accessDEC = new MPI_AccessDEC(local_group,distant_group,getAsynchronous()); + _accessDEC->SetTimeInterpolator(getTimeInterpolationMethod()); _send_proc_offsets.resize(_union_group->size()+1,0); _recv_proc_offsets.resize(_union_group->size()+1,0); @@ -215,7 +218,7 @@ void MxN_Mapping::sendRecv(double* sendfield, MEDMEM::FIELD& field) cons } //communication phase - switch (_allToAllMethod) { + switch (getAllToAllMethod()) { case Native: { const MPI_Comm* comm = group->getComm(); @@ -247,8 +250,8 @@ void MxN_Mapping::sendRecv(double* sendfield, MEDMEM::FIELD& field) cons } // if (sendbuf!=0) delete[] sendbuf; - if (sendbuf!=0 && _allToAllMethod == Native) delete[] sendbuf; - if (recvbuf!=0) delete[] recvbuf; + if (sendbuf!=0 && getAllToAllMethod()== Native) delete[] sendbuf; + if (recvbuf !=0) delete[] recvbuf; delete[] sendcounts; delete[] recvcounts; delete[] senddispls; @@ -309,7 +312,7 @@ void MxN_Mapping::reverseSendRecv(double* recvfield, MEDMEM::FIELD& fiel // } //communication phase - switch (_allToAllMethod) { + switch (getAllToAllMethod()) { case Native: { const MPI_Comm* comm = group->getComm(); @@ -344,8 +347,8 @@ void MxN_Mapping::reverseSendRecv(double* recvfield, MEDMEM::FIELD& fiel // if (sendbuf!=0) delete[] sendbuf; - if (sendbuf!=0 && _allToAllMethod == Native) delete[] sendbuf; - if (recvbuf!=0) delete[] recvbuf; + if (sendbuf!=0 && getAllToAllMethod() == Native) delete[] sendbuf; + if (recvbuf!=0) delete[] recvbuf; delete[] sendcounts; delete[] recvcounts; delete[] senddispls; diff --git a/src/ParaMEDMEM/MxN_Mapping.hxx b/src/ParaMEDMEM/MxN_Mapping.hxx index edc039c40..694da46a4 100644 --- a/src/ParaMEDMEM/MxN_Mapping.hxx +++ b/src/ParaMEDMEM/MxN_Mapping.hxx @@ -5,33 +5,31 @@ #include "MEDMEM_Field.hxx" #include "MPI_AccessDEC.hxx" - +#include "DECOptions.hxx" namespace ParaMEDMEM { -typedef enum{Native,PointToPoint} AllToAllMethod; class ProcessorGroup; -class MxN_Mapping + class MxN_Mapping : public DECOptions { public: MxN_Mapping(); - MxN_Mapping(const ProcessorGroup& local_group, const ProcessorGroup& distant_group); + MxN_Mapping(const ProcessorGroup& local_group, const ProcessorGroup& distant_group, const DECOptions& dec_options); virtual ~MxN_Mapping(); void addElementFromSource(int distant_proc, int distant_elem); void prepareSendRecv(); void sendRecv(MEDMEM::FIELD& field); void sendRecv(double* field, MEDMEM::FIELD& field) const ; void reverseSendRecv(double* field, MEDMEM::FIELD& field) const ; - void setAllToAllMethod(const AllToAllMethod& method){_allToAllMethod = method ;}; + MPI_AccessDEC* getAccessDEC(){return _accessDEC;} private : // ProcessorGroup& _local_group; // ProcessorGroup& _distant_group; ProcessorGroup* _union_group; MPI_AccessDEC * _accessDEC; - AllToAllMethod _allToAllMethod ; int _nb_comps; std::vector > _sending_ids; std::vector _recv_ids; diff --git a/src/ParaMEDMEM/Test/Makefile.am b/src/ParaMEDMEM/Test/Makefile.am index 30bc30b24..ca46ba3ab 100644 --- a/src/ParaMEDMEM/Test/Makefile.am +++ b/src/ParaMEDMEM/Test/Makefile.am @@ -77,7 +77,8 @@ endif # Executables targets bin_PROGRAMS = TestParaMEDMEM \ TestMPIAccessDEC \ - TestMPIAccess + TestMPIAccess\ + test_bug dist_TestParaMEDMEM_SOURCES = TestParaMEDMEM.cxx dist_TestMPIAccessDEC_SOURCES = TestMPIAccessDEC.cxx diff --git a/src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx b/src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx index 4ebd16f9d..2d1bface1 100644 --- a/src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx +++ b/src/ParaMEDMEM/Test/ParaMEDMEMTest.hxx @@ -38,18 +38,19 @@ class ParaMEDMEMTest : public CppUnit::TestFixture CPPUNIT_TEST(testBlockTopology_constructor); CPPUNIT_TEST(testBlockTopology_serialize); CPPUNIT_TEST(testIntersectionDEC_2D); - CPPUNIT_TEST(testSynchronousEqualIntersectionWithoutInterpNativeDEC_2D); - CPPUNIT_TEST(testSynchronousEqualIntersectionWithoutInterpDEC_2D); - CPPUNIT_TEST(testSynchronousEqualIntersectionDEC_2D); - CPPUNIT_TEST(testSynchronousFasterSourceIntersectionDEC_2D); - CPPUNIT_TEST(testSynchronousSlowerSourceIntersectionDEC_2D); - CPPUNIT_TEST(testSynchronousSlowSourceIntersectionDEC_2D); - CPPUNIT_TEST(testSynchronousFastSourceIntersectionDEC_2D); + + CPPUNIT_TEST(testSynchronousEqualIntersectionWithoutInterpNativeDEC_2D); + CPPUNIT_TEST(testSynchronousEqualIntersectionWithoutInterpDEC_2D); + CPPUNIT_TEST(testSynchronousEqualIntersectionDEC_2D); + CPPUNIT_TEST(testSynchronousFasterSourceIntersectionDEC_2D); + CPPUNIT_TEST(testSynchronousSlowerSourceIntersectionDEC_2D); + CPPUNIT_TEST(testSynchronousSlowSourceIntersectionDEC_2D); + CPPUNIT_TEST(testSynchronousFastSourceIntersectionDEC_2D); CPPUNIT_TEST(testAsynchronousEqualIntersectionDEC_2D); - CPPUNIT_TEST(testAsynchronousFasterSourceIntersectionDEC_2D); - CPPUNIT_TEST(testAsynchronousSlowerSourceIntersectionDEC_2D); - CPPUNIT_TEST(testAsynchronousSlowSourceIntersectionDEC_2D); - CPPUNIT_TEST(testAsynchronousFastSourceIntersectionDEC_2D); + CPPUNIT_TEST(testAsynchronousFasterSourceIntersectionDEC_2D); + CPPUNIT_TEST(testAsynchronousSlowerSourceIntersectionDEC_2D); + CPPUNIT_TEST(testAsynchronousSlowSourceIntersectionDEC_2D); + CPPUNIT_TEST(testAsynchronousFastSourceIntersectionDEC_2D); #ifdef MED_ENABLE_FVM //can be added again after FVM correction for 2D // CPPUNIT_TEST(testNonCoincidentDEC_2D); diff --git a/src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx b/src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx index e74e1a90b..b36a262a1 100644 --- a/src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx +++ b/src/ParaMEDMEM/Test/ParaMEDMEMTest_IntersectionDEC.cxx @@ -183,7 +183,7 @@ void ParaMEDMEMTest::testIntersectionDEC_2D() field_before_int = parafield->getVolumeIntegral(1); dec.synchronize(); cout<<"DEC usage"<write(MED_DRIVER,"./sourcesquareb"); @@ -205,15 +205,20 @@ void ParaMEDMEMTest::testIntersectionDEC_2D() filename<<"./sourcesquare_"<myRank()+1; aRemover.Register(filename.str().c_str()); field_after_int = parafield->getVolumeIntegral(1); - CPPUNIT_ASSERT_DOUBLES_EQUAL(field_before_int, field_after_int, 1e-6); - + + + // MPI_Bcast(&field_before_int,1,MPI_DOUBLE,0,MPI_COMM_WORLD); +// MPI_Bcast(&field_after_int,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + + CPPUNIT_ASSERT_DOUBLES_EQUAL(field_before_int, field_after_int, 1e-6); + } //attaching a DEC to the target group if (target_group->containsMyRank()) { dec.synchronize(); - dec.setOption("ForcedRenormalization",false); + dec.setForcedRenormalization(false); dec.recvData(); paramesh->write(MED_DRIVER, "./targetsquareb"); @@ -226,11 +231,17 @@ void ParaMEDMEMTest::testIntersectionDEC_2D() dec.sendData(); paramesh->write(MED_DRIVER, "./targetsquare"); parafield->write(MED_DRIVER, "./targetsquare", "boundary"); - if (target_group->myRank()==0) + if (target_group->myRank()==0) aRemover.Register("./targetsquareb"); filename<<"./targetsquareb_"<myRank()+1; aRemover.Register(filename.str().c_str()); + // double field_before_int, field_after_int; +// MPI_Bcast(&field_before_int,1,MPI_DOUBLE,0,MPI_COMM_WORLD); +// MPI_Bcast(&field_after_int,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + +// CPPUNIT_ASSERT_DOUBLES_EQUAL(field_before_int, field_after_int, 1e-6); + } delete source_group; @@ -238,7 +249,7 @@ void ParaMEDMEMTest::testIntersectionDEC_2D() delete self_group; delete mesh; delete paramesh; - delete parafield; + delete parafield; delete parasupport; delete [] value; delete icocofield; @@ -397,6 +408,8 @@ void ParaMEDMEMTest::testAsynchronousIntersectionDEC_2D(double dtA, double tmaxA icocofield=new ICoCo::MEDField(paramesh,parafield); dec.attachLocalField(icocofield); + + } //loading the geometry for the target group @@ -436,23 +449,23 @@ void ParaMEDMEMTest::testAsynchronousIntersectionDEC_2D(double dtA, double tmaxA if (source_group->containsMyRank()) { cout<<"DEC usage"<getField()->getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS); for (int i=0; icontainsMyRank()) { cout<<"DEC usage"< times; for (double time=0; timegetVolumeIntegral(1) + double vi = parafield->getVolumeIntegral(1); + cout << "testAsynchronousIntersectionDEC_2D" << rank << " time " << time + << " VolumeIntegral " << vi << " time*10000 " << time*10000 << endl ; - CPPUNIT_ASSERT_DOUBLES_EQUAL(parafield->getVolumeIntegral(1),time*10000,0.001); + + CPPUNIT_ASSERT_DOUBLES_EQUAL(vi,time*10000,0.001); } }