std::vector<int> _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)
{
}
- _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);
-
}
bool intersects = true;
for (int idim=0; idim<dim; idim++)
{
- if (bb_ptr[idim*2]-bb[idim*2+1]>-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)
bool intersects = true;
for (int idim=0; idim<dim; idim++)
{
- if (bb_ptr[idim*2]-xx[idim]>1e-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)
*
*
*/
+#include "InterpolationOptions.hxx"
+
namespace INTERP_KERNEL
{
template<class TrueMainInterpolator>
- 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<class MatrixType, class MyMeshType>
void interpolateMeshes(const MyMeshType& mesh1, const MyMeshType& mesh2, MatrixType& result)
{
public:
Interpolation2D() { }
+ Interpolation2D(const InterpolationOptions& io):InterpolationPlanar<Interpolation2D>(io){};
public:
bool doRotate() const { return false; }
double medianPlane() const { return 0.; }
#include "Interpolation.hxx"
#include "NormalizedUnstructuredMesh.hxx"
#include "IntersectorHexa.hxx"
-#include "OptionManager.hxx"
+#include "InterpolationOptions.hxx"
namespace INTERP_KERNEL
{
- class Interpolation3D : public Interpolation<Interpolation3D>,OptionManager
+ class Interpolation3D : public Interpolation<Interpolation3D>
{
public:
Interpolation3D();
+ Interpolation3D(const InterpolationOptions& io);
template<class MatrixType, class MyMeshType>
void interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result);
private:
#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
*/
Interpolation3D::Interpolation3D()
{
- registerOption(&_splitting_policy,"SplittingPolicy",GENERAL_48);
+ }
+ Interpolation3D::Interpolation3D(const InterpolationOptions& io):Interpolation<Interpolation3D>(io)
+ {
}
/**
case 8:
- SplittingPolicy sp;
- getOption("SplittingPolicy",sp);
- intersector = new IntersectorHexa<MyMeshType>(srcMesh, targetMesh, targetIdx,sp);
+ intersector = new IntersectorHexa<MyMeshType>(srcMesh, targetMesh, targetIdx,getSplittingPolicy());
break;
default:
,_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<Interpolation3DSurf>(io)
+ {
+ }
+
+
/**
\brief Function used to set the options for the intersection calculation
\details The following options can be modified:
--- /dev/null
+#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
#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 RealPlanar>
- class InterpolationPlanar : public Interpolation< InterpolationPlanar<RealPlanar> >, public OptionManager
+ class InterpolationPlanar : public Interpolation< InterpolationPlanar<RealPlanar> >
{
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<class MatrixType, class MyMeshType>
#define __INTERPOLATIONPLANAR_TXX__
#include "InterpolationPlanar.hxx"
-#include "OptionManager.hxx"
+#include "InterpolationOptions.hxx"
#include "PlanarIntersector.hxx"
#include "PlanarIntersector.txx"
#include "TriangulationIntersector.hxx"
* two local meshes in two dimensions. Meshes can contain mixed triangular and quadrangular elements.
*/
template<class RealPlanar>
- InterpolationPlanar<RealPlanar>::InterpolationPlanar():_printLevel(0),_precision(DEFAULT_PRECISION),_dimCaracteristic(1),
- _intersectionType(Triangulation),_orientation(1)
+ InterpolationPlanar<RealPlanar>::InterpolationPlanar():_dimCaracteristic(1)
+
{
- registerOption(&_precision, "Precision",DEFAULT_PRECISION);
- registerOption(&_intersectionType, "IntersectionType",Triangulation);
- registerOption(&_printLevel, "PrintLevel",0);
- registerOption(&_orientation, "IntersectionSign",1);
+ }
+ template<class RealPlanar>
+ InterpolationPlanar<RealPlanar>::InterpolationPlanar(const InterpolationOptions& io):Interpolation< InterpolationPlanar<RealPlanar> >(io),_dimCaracteristic(1)
+
+ {
}
+
/**
* \brief Function used to set the options for the intersection calculation
* \details The following options can be modified:
template<class RealPlanar>
void InterpolationPlanar<RealPlanar>::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);
}
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;
PlanarIntersector<MyMeshType>* intersector;
- switch (_intersectionType)
+ switch (InterpolationOptions::getIntersectionType())
{
case Triangulation:
- intersector=new TriangulationIntersector<MyMeshType>(myMesh_P,myMesh_S, _dimCaracteristic,_precision,
- medianPlane(), _printLevel);
+ intersector=new TriangulationIntersector<MyMeshType>(
+ myMesh_P,
+ myMesh_S,
+ _dimCaracteristic,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrintLevel());
break;
case Convex:
- intersector=new ConvexIntersector<MyMeshType>(myMesh_P,myMesh_S, _dimCaracteristic, _precision,
- medianPlane(), doRotate(), _printLevel);
+ intersector=new ConvexIntersector<MyMeshType>(
+ myMesh_P,
+ myMesh_S,
+ _dimCaracteristic,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getDoRotate(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrintLevel());
break;
case Geometric2D:
- intersector=new Geometric2DIntersector<MyMeshType>(myMesh_P, myMesh_S, _dimCaracteristic, _precision);
+ intersector=new Geometric2DIntersector<MyMeshType>(myMesh_P, myMesh_S, _dimCaracteristic, InterpolationOptions::getPrecision());
break;
// case MEDMEM::Generic:
//intersector=new GenericIntersector<SPACEDIM>(myMesh_P,myMesh_S, _DimCaracteristic,_Precision,
// -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<ConnType,numPol>::indFC(i_S),surf));
counter++;
}
/* DEBUG prints */
/***********************************/
- if (_printLevel >=1)
+ if (InterpolationOptions::getPrintLevel() >=1)
{
long end_intersection=clock();
// if (_printLevel >=2)
#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
*
*/
template<class MyMeshType>
- class IntersectorHexa : public TargetIntersector<typename MyMeshType::MyConnType>, public OptionManager
+ class IntersectorHexa : public TargetIntersector<typename MyMeshType::MyConnType>, public InterpolationOptions
{
public:
+++ /dev/null
-#ifndef _OPTION_MANAGER_HXX
-#define _OPTION_MANAGER_HXX
-
-#include <string>
-#include <map>
-
-#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<T*>& 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 <class T> void registerOption( T * var, const std::string& name, T defaut)
- {
- OptionGeneric * newoption = new Option<T>( defaut, var);
- _optionList.insert(make_pair(name, newoption));
- }
-
- template <class T> void getOption(const std::string& name, T& var)
- {
- if(_optionList.find(name) != _optionList.end())
- {
- Option<T>* option_ptr = dynamic_cast<Option<T>*>(_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 <class T> void setOption(const std::string& name, T value)
- {
- if(_optionList.find(name) != _optionList.end())
- {
- Option<T>* option_ptr = dynamic_cast<Option<T>*>(_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<double>* option_double_ptr = dynamic_cast<Option<double>*> (_optionList.find(name)->second);
- if (option_double_ptr!=0)
- setOption(name,(double) value);
- else
- {
- Option<int>* option_ptr =dynamic_cast<Option<int>*>(_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
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)
{
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) \
#include "MultiElement2DTests.hxx"
#include "SingleElementPlanarTests.hxx"
#include "QuadraticPlanarInterpTest.hxx"
-
+#include "InterpolationOptionsTest.hxx"
using namespace INTERP_TEST;
//--- Registers the fixture into the 'registry'
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
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
@{
_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()
#ifndef DEC_HXX_
#define DEC_HXX_
-#include "OptionManager.hxx"
#include <MEDMEM_Field.hxx>
#include <NormalizedUnstructuredMesh.hxx>
+#include "DECOptions.hxx"
namespace ICoCo
{
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<double>* field);
void attachLocalField(const ParaFIELD* field);
ProcessorGroup* _target_group;
const CommInterface* _comm_interface;
- bool _forced_renormalization_flag;
bool _owns_field;
private:
ICoCo::Field* _icoco_field;
--- /dev/null
+#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
{
_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)
_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");
}
{
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];
}
for (int i=0; i<dim; i++)
{
- bbtemp[i*2]=bb1[i*2]-deltamax*_adjustment_eps;
- bbtemp[i*2+1]=bb1[i*2+1]+deltamax*_adjustment_eps;
+ bbtemp[i*2]=bb1[i*2]-deltamax*adjustment_eps;
+ bbtemp[i*2+1]=bb1[i*2+1]+deltamax*adjustment_eps;
}
for (int idim=0; idim < dim; idim++)
#include <vector>
#include <set>
-#include "OptionManager.hxx"
+#include "InterpolationOptions.hxx"
namespace MEDMEM
{
class InterpolationMatrix;
- class ElementLocator: public INTERP_KERNEL::OptionManager
+ class ElementLocator: public INTERP_KERNEL::InterpolationOptions
{
public:
ElementLocator(const ParaMESH& mesh, const ProcessorGroup& distant_group);
const ProcessorGroup& _local_group;
ProcessorGroup* _union_group;
std::vector<int> _distant_proc_ids;
- double _adjustment_eps;
void _computeBoundingBoxes();
bool _intersectsBoundingBox(int irank);
#include "Interpolation3DSurf.txx"
#include "Interpolation3D.txx"
#include "MEDNormalizedUnstructuredMesh.hxx"
+#include "InterpolationOptions.hxx"
/*! \class InterpolationMatrix
\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<nbelems+1; i++)
- _row_offsets[i]=0;
+ InterpolationMatrix::InterpolationMatrix(
+ const ParaMEDMEM::ParaMESH& source_support,
+ const ProcessorGroup& source_group,
+ const ProcessorGroup& target_group,
+ const DECOptions& dec_options,
+ const INTERP_KERNEL::InterpolationOptions& interp_options):
+ _source_support(*source_support.getMesh()),
+ _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.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ _row_offsets.resize(nbelems+1);
+ for (int i=0; i<nbelems+1; i++)
+ _row_offsets[i]=0;
- _coeffs.resize(nbelems);
-}
-
-InterpolationMatrix::~InterpolationMatrix()
-{
-}
-
-/*!
-\brief Adds the contribution of a distant subdomain to the interpolation matrix.
-
-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)
-{
- 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<map<int,double> > 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<map<int,double> > 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()="<<surfaces.size()<<" source_size="<<source_size<<endl;
- throw MEDEXCEPTION("uncoherent number of rows in interpolation matrix");
- }
+ if (surfaces.size() != source_size)
+ {
+ cout<<"surfaces.size()="<<surfaces.size()<<" source_size="<<source_size<<endl;
+ throw MEDEXCEPTION("uncoherent number of rows in interpolation matrix");
+ }
- //computing the vectors containing the source and target element volumes
- MEDMEM::SUPPORT target_support(&distant_support,"all cells", MED_EN::MED_CELL);
- MEDMEM::FIELD<double>* target_triangle_surf = getSupportVolumes(target_support);
- MEDMEM::SUPPORT source_support (const_cast<MEDMEM::MESH*>(&_source_support),"all cells", MED_EN::MED_CELL);
- MEDMEM::FIELD<double>* source_triangle_surf = getSupportVolumes(source_support);
-
- //storing the source volumes
- _source_volume.resize(source_size);
- for (int i=0; i<source_size;i++)
- _source_volume[i]=source_triangle_surf->getValueIJ(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<double>* target_triangle_surf = getSupportVolumes(target_support);
+ MEDMEM::SUPPORT source_support (const_cast<MEDMEM::MESH*>(&_source_support),"all cells", MED_EN::MED_CELL);
+ MEDMEM::FIELD<double>* source_triangle_surf = getSupportVolumes(source_support);
+
+ //storing the source volumes
+ _source_volume.resize(source_size);
+ for (int i=0; i<source_size;i++)
+ _source_volume[i]=source_triangle_surf->getValueIJ(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<int,double>::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<pair<int,int> >::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<int,double>::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<pair<int,int> >::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<double>& field) const
-{
- vector<double> 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<double>& field) const
+ {
+ vector<double> 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; irow<nbrows; irow++)
- for (int icomp=0; icomp< nbcomp; icomp++)
- {
- double coeff_row = field.getValueIJ(irow+1,icomp+1);
- for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
- {
- 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;
- }
- }
+ for (int irow=0; irow<nbrows; irow++)
+ for (int icomp=0; icomp< nbcomp; icomp++)
+ {
+ double coeff_row = field.getValueIJ(irow+1,icomp+1);
+ for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
+ {
+ 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;
+ }
+ }
- // performing VT^(-1).(W.S)
- // where VT^(-1) is the inverse of the diagonal matrix containing
- // the volumes of target cells
+ // 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; icomp<nbcomp; icomp++)
+ target_value[i*nbcomp+icomp] /= _target_volume[i];
+ }
+
+ if (_target_group.containsMyRank())
+ {
+ int nbelems=field.getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
+ int nbcomp = field.getNumberOfComponents();
+ double* value = const_cast<double*> (field.getValue());
+ for (int i=0; i<nbelems*nbcomp;i++)
+ value[i]=0.0;
+ }
+ //on source side : sending T=VT^(-1).(W.S)
+ //on target side :: receiving T and storing it in field
+ _mapping.sendRecv(&target_value[0],field);
+
+ }
+ /*!
+ \brief performs s=WTt, where t is the target field, s is the source field, WT is the transpose matrix from W
- for (int i=0; i<_col_offsets.size();i++)
- for (int icomp=0; icomp<nbcomp; icomp++)
- target_value[i*nbcomp+icomp] /= _target_volume[i];
- }
+ 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.
- if (_target_group.containsMyRank())
+ \param field source field on processors involved on the source side, target field on processors on the target side
+ */
+ void InterpolationMatrix::transposeMultiply(MEDMEM::FIELD<double>& field) const
{
- int nbelems=field.getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
int nbcomp = field.getNumberOfComponents();
- double* value = const_cast<double*> (field.getValue());
- for (int i=0; i<nbelems*nbcomp;i++)
- value[i]=0.0;
- }
- //on source side : sending T=VT^(-1).(W.S)
- //on target side :: receiving T and storing it in field
- _mapping.sendRecv(&target_value[0],field);
-
-}
-/*!
- \brief performs s=WTt, where t is the target field, s is the source field, WT is the transpose matrix from W
+ vector<double> 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<double>& field) const
-{
- int nbcomp = field.getNumberOfComponents();
- vector<double> 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<double> target_value(nbrows*nbcomp,0.0);
+ vector<double> target_value(nbrows*nbcomp,0.0);
- //performing WT.T
- //WT is W transpose
- //T is the target vector
- for (int irow=0; irow<nbrows; irow++)
- for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
- {
- int colid= _coeffs[irow][icol-_row_offsets[irow]].first;
- double value = _coeffs[irow][icol-_row_offsets[irow]].second;
+ //performing WT.T
+ //WT is W transpose
+ //T is the target vector
+ for (int irow=0; irow<nbrows; irow++)
+ for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
+ {
+ int colid= _coeffs[irow][icol-_row_offsets[irow]].first;
+ double value = _coeffs[irow][icol-_row_offsets[irow]].second;
- for (int icomp=0; icomp<nbcomp; icomp++)
- {
- double coeff_row = source_value[(colid-1)*nbcomp+icomp];
- target_value[irow*nbcomp+icomp]+=value*coeff_row;
- }
- }
- //performing VS^(-1).(WT.T)
- //VS^(-1) is the inverse of the diagonal matrix storing
- //volumes of the source cells
- for (int i=0; i<nbrows; i++)
- for (int icomp=0; icomp<nbcomp; icomp++)
- target_value[i*nbcomp+icomp]/=_source_volume[i];
+ for (int icomp=0; icomp<nbcomp; icomp++)
+ {
+ double coeff_row = source_value[(colid-1)*nbcomp+icomp];
+ target_value[irow*nbcomp+icomp]+=value*coeff_row;
+ }
+ }
+ //performing VS^(-1).(WT.T)
+ //VS^(-1) is the inverse of the diagonal matrix storing
+ //volumes of the source cells
+ for (int i=0; i<nbrows; i++)
+ for (int icomp=0; icomp<nbcomp; icomp++)
+ target_value[i*nbcomp+icomp]/=_source_volume[i];
- //storing S= VS^(-1).(WT.T)
- for (int irow=0; irow<nbrows; irow++)
- for (int icomp=0; icomp<nbcomp; icomp++)
- field.setValueIJ(irow+1,icomp+1,target_value[irow*nbcomp+icomp]);
- }
+ //storing S= VS^(-1).(WT.T)
+ for (int irow=0; irow<nbrows; irow++)
+ for (int icomp=0; icomp<nbcomp; icomp++)
+ field.setValueIJ(irow+1,icomp+1,target_value[irow*nbcomp+icomp]);
+ }
-}
+ }
-/*!
-\brief returns the volumes of the cells underlying the field \a field
+ /*!
+ \brief returns the volumes of the cells underlying the field \a field
-For 2D geometries, the returned field contains the areas.
-For 3D geometries, the returned field contains the volumes.
+ For 2D geometries, the returned field contains the areas.
+ For 3D geometries, the returned field contains the volumes.
-\param field field on which cells the volumes are required
-\return field containing the volumes
-*/
- MEDMEM::FIELD<double>* 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<double>* 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");
+ }
+ }
}
#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);
void transposeMultiply(MEDMEM::FIELD<double>&)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:
vector<pair<int,int> > _col_offsets;
const MEDMEM::MESH& _source_support;
MxN_Mapping _mapping;
- string _method;
+
const ProcessorGroup& _source_group;
const ProcessorGroup& _target_group;
vector<double> _target_volume;
#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
@{
*/
- IntersectionDEC::IntersectionDEC()
+IntersectionDEC::IntersectionDEC()
{
}
*/
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()
{
const ParaMEDMEM::ParaMESH* para_mesh = _local_field->getSupport()->getMesh();
//cout <<"size of Interpolation Matrix"<<sizeof(InterpolationMatrix)<<endl;
- _interpolation_matrix = new InterpolationMatrix (*para_mesh, *_source_group,*_target_group,*this);
- _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())
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;
{
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++)
else if (_target_group->containsMyRank())
{
_interpolation_matrix->multiply(*_local_field->getField());
- if (_forced_renormalization_flag)
+ if (getForcedRenormalization())
renormalizeTargetField();
}
{
_interpolation_matrix->multiply(*_local_field->getField());
- if (_forced_renormalization_flag)
+ if (getForcedRenormalization())
renormalizeTargetField();
}
#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();
//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;
-
-
};
}
namespace ParaMEDMEM {
- typedef enum{WithoutTimeInterp,LinearTimeInterp} TimeInterpolationMethod;
+ // typedef enum{WithoutTimeInterp,LinearTimeInterp} TimeInterpolationMethod;
class MPI_AccessDEC {
-#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);
}
//communication phase
- switch (_allToAllMethod) {
+ switch (getAllToAllMethod()) {
case Native:
{
const MPI_Comm* comm = group->getComm();
}
// 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;
// }
//communication phase
- switch (_allToAllMethod) {
+ switch (getAllToAllMethod()) {
case Native:
{
const MPI_Comm* comm = group->getComm();
// 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;
#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<double>& field);
void sendRecv(double* field, MEDMEM::FIELD<double>& field) const ;
void reverseSendRecv(double* field, MEDMEM::FIELD<double>& 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<pair<int,int> > _sending_ids;
std::vector<int> _recv_ids;
# Executables targets
bin_PROGRAMS = TestParaMEDMEM \
TestMPIAccessDEC \
- TestMPIAccess
+ TestMPIAccess\
+ test_bug
dist_TestParaMEDMEM_SOURCES = TestParaMEDMEM.cxx
dist_TestMPIAccessDEC_SOURCES = TestMPIAccessDEC.cxx
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);
field_before_int = parafield->getVolumeIntegral(1);
dec.synchronize();
cout<<"DEC usage"<<endl;
- dec.setOption("ForcedRenormalization",false);
+ dec.setForcedRenormalization(false);
dec.sendData();
paramesh->write(MED_DRIVER,"./sourcesquareb");
filename<<"./sourcesquare_"<<source_group->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");
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_"<<target_group->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;
delete self_group;
delete mesh;
delete paramesh;
- delete parafield;
+ delete parafield;
delete parasupport;
delete [] value;
delete icocofield;
icocofield=new ICoCo::MEDField(paramesh,parafield);
dec.attachLocalField(icocofield);
+
+
}
//loading the geometry for the target group
if (source_group->containsMyRank())
{
cout<<"DEC usage"<<endl;
- dec.setOption("Asynchronous",Asynchronous);
+ dec.setAsynchronous(Asynchronous);
if ( WithInterp ) {
- dec.setOption("TimeInterpolation",LinearTimeInterp);
+ dec.setTimeInterpolationMethod(LinearTimeInterp);
}
if ( WithPointToPoint ) {
- dec.setOption("AllToAllMethod",PointToPoint);
+ dec.setAllToAllMethod(PointToPoint);
}
else {
- dec.setOption("AllToAllMethod",Native);
+ dec.setAllToAllMethod(Native);
}
dec.synchronize();
- dec.setOption("ForcedRenormalization",false);
+ dec.setForcedRenormalization(false);
for (double time=0; time<tmaxA+1e-10; time+=dtA)
{
cout << "testAsynchronousIntersectionDEC_2D" << rank << " time " << time
<< " dtA " << dtA << " tmaxA " << tmaxA << endl ;
- if ( time+dtA < tmaxA+1e-10 ) {
+ if ( time+dtA < tmaxA+1e-7 ) {
dec.sendData( time , dtA );
}
else {
int nb_local=parafield->getField()->getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
for (int i=0; i<nb_local;i++)
value[i]= time+dtA;
+
}
}
if (target_group->containsMyRank())
{
cout<<"DEC usage"<<endl;
- dec.setOption("Asynchronous",Asynchronous);
+ dec.setAsynchronous(Asynchronous);
if ( WithInterp ) {
- dec.setOption("TimeInterpolation",LinearTimeInterp);
+ dec.setTimeInterpolationMethod(LinearTimeInterp);
}
if ( WithPointToPoint ) {
- dec.setOption("AllToAllMethod",PointToPoint);
+ dec.setAllToAllMethod(PointToPoint);
}
else {
- dec.setOption("AllToAllMethod",Native);
+ dec.setAllToAllMethod(Native);
}
dec.synchronize();
- dec.setOption("ForcedRenormalization",false);
+ dec.setForcedRenormalization(false);
vector<double> times;
for (double time=0; time<tmaxB+1e-10; time+=dtB)
{
cout << "testAsynchronousIntersectionDEC_2D" << rank << " time " << time
<< " dtB " << dtB << " tmaxB " << tmaxB << endl ;
dec.recvData( time );
- cout << "testAsynchronousIntersectionDEC_2D" << rank << " time " << time
- << " VolumeIntegral " << parafield->getVolumeIntegral(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);
}
}
