1 # ifndef INTERPOLATION_HXX
2 # define INTERPOLATION_HXX
4 //template < class T> class FIELD;
5 //template < int > class Wrapper_Nuage;
6 //template < int > class Wrapper_Noeud;
7 //template <class ,class ,int ,int > class dTree;
10 #include "utilities.h"
11 #include "MEDMEM_Exception.hxx"
12 #include "MEDMEM_define.hxx"
14 #include "MEDMEM_InterpolationHighLevelObjects.hxx"
15 #include "MEDMEM_Mesh.hxx"
16 #include "MEDMEM_Field.hxx"
21 template <int DIMENSION=3> class INTERPOLATION
25 FIELD<double> * _fromField;
26 FIELD<double> * _toField;
30 Meta_Wrapper<DIMENSION> * _fromWrapper;
31 Meta_Wrapper<DIMENSION> * _toWrapper;
33 Meta_Mapping <DIMENSION> * _mapping;
35 // only used when multi timestep are interpolated
36 // but always coherent
38 int _isConvexFromMesh;
44 // Initialize INTERPOLATION in order to get :
45 // 1- the node number in the MESH <fromMesh> which
46 // is the nearest from a given one ( use method : getNearestNode( double * node ) );
47 // 2- the cell number (if exists) in the MESH <fromMesh> which
48 // contains a specified node ( use method : getContainingCell ( double * node) )
49 INTERPOLATION(const MESH & fromMesh );
50 // Initialize INTERPOLATION in order to get :
51 // 1- the complete mapping ( use methode : getMapping() )
52 // 2- the functionalities above
53 INTERPOLATION(const MESH & fromMesh,const MESH & toMesh );
54 // Initialize INTERPOLATION in order to get the interpolation of <field> on <toMesh>
55 // Moreover, all the others functionalities are so available
56 INTERPOLATION(const FIELD<double> & fromField, const MESH & toMesh);
60 // Get the node number in the MESH <fromMesh> which is the nearest from a given one
61 int getNearestNode ( double * node );
62 // Get the cell number (if exists) in the MESH <fromMesh> which contains a specified node
63 int getContainingCell ( double * node , int beginingCell=0, int flagIsConvexMesh=0 );
64 // Get the complete mapping, defaultly, fromMesh is supposed to be non-convex, if it is false, set flagIsConvexMesh to 1
65 vector<int> getMapping ( int flagIsConvexMesh=0 );
66 // Get the interpolated field toField
67 FIELD<double> * interpolate( /*med_interpolation_type*/ int itype,int flagIsConvexFromMesh=0);
68 // reset the <from> parameters in order not to redo the mapping (if the mesh are identical)
69 // and then get the interpoated field toField
70 // this method is specifictly used on multi-timestep (or order number) fields
71 // it has only to be used after the first step, the interpolation paramaters are the same for every step
72 FIELD<double> * interpolateNextStep(const FIELD <double> &nextFromField ,int & flagNewMapping);
76 template <int DIMENSION> void INTERPOLATION<DIMENSION>::init() {
78 const char * LOC = "INTERPOLATION::init(): ";
81 _fromField = ( FIELD<double> * ) NULL;
82 _toField = ( FIELD<double> * ) NULL;
83 _fromMesh = ( MESH * ) NULL;
84 _toMesh = ( MESH * ) NULL;
85 _fromWrapper = ( Meta_Wrapper<DIMENSION> * ) NULL;
86 _toWrapper = ( Meta_Wrapper<DIMENSION> * ) NULL;
87 _mapping = ( Meta_Mapping<DIMENSION> * ) NULL;
89 _isConvexFromMesh = UNDEFINED ;
94 template <int DIMENSION> INTERPOLATION<DIMENSION>::INTERPOLATION(const MESH & fromMesh ) {
96 const char * LOC = "INTERPOLATION::INTERPOLATION(MESH * fromMesh ) : ";
101 _fromMesh=const_cast<MESH * > (&fromMesh);
103 if (! _fromMesh ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromMesh is a NULL pointer !")) ;
105 int spaceDimension = _fromMesh->getSpaceDimension();
106 if (spaceDimension != DIMENSION ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The spaceDimension of mesh |" << _fromMesh->getName() << "| is |" << spaceDimension << "| and should be |" << DIMENSION << "|" << endl)) ;
108 _fromWrapper = new Meta_Wrapper<DIMENSION>(_fromMesh->getNumberOfNodes(),
109 const_cast<double *> (_fromMesh->getCoordinates(MED_FULL_INTERLACE)),
110 const_cast<CONNECTIVITY *> (_fromMesh->getConnectivityptr())
113 _mapping = new Meta_Mapping<DIMENSION> (_fromWrapper);
118 template <int DIMENSION> INTERPOLATION<DIMENSION>::INTERPOLATION(const MESH & fromMesh,const MESH & toMesh ) {
120 const char * LOC = "INTERPOLATION::INTERPOLATION(MESH * fromMesh,,const MESH & toMesh) : ";
125 _fromMesh = const_cast<MESH * > ( &fromMesh );
126 _toMesh = const_cast<MESH * > ( &toMesh );
128 if (! _fromMesh ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromMesh is a NULL pointer !")) ;
129 if (! _toMesh ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toMesh is a NULL pointer !")) ;
131 int fromSpaceDimension = _fromMesh->getSpaceDimension();
132 int toSpaceDimension = _toMesh->getSpaceDimension();
134 if (fromSpaceDimension != DIMENSION ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The spaceDimension of mesh |" << _fromMesh->getName() << "| is |" << spaceDimension << "| and should be |" << DIMENSION << "|" << endl)) ;
135 if ( toSpaceDimension != DIMENSION ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The spaceDimension of mesh |" << _toMesh->getName() << "| is |" << spaceDimension << "| and should be |" << DIMENSION << "|" << endl)) ;
137 _fromWrapper = new Meta_Wrapper<DIMENSION>(_fromMesh->getNumberOfNodes(),
138 const_cast<double *> (_fromMesh->getCoordinates(MED_FULL_INTERLACE)),
139 const_cast<CONNECTIVITY *> (_fromMesh->getConnectivityptr())
142 _toWrapper = new Meta_Wrapper<DIMENSION>(_toMesh->getNumberOfNodes(),
143 const_cast<double *> (_toMesh->getCoordinates(MED_FULL_INTERLACE))
146 _mapping = new Meta_Mapping<DIMENSION> (_fromWrapper);
151 template <int DIMENSION> INTERPOLATION<DIMENSION>::INTERPOLATION(const FIELD<double> & fromField,const MESH & toMesh) {
153 const char * LOC = "INTERPOLATION(const FIELD<double> & field,const MESH & toMesh) : ";
158 _toMesh = const_cast<MESH *>(&toMesh);
159 _fromField = const_cast<FIELD<double> *>(&fromField);
161 if ( ! _toMesh ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toMesh is a NULL pointer !")) ;
162 if ( ! _fromField ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"field is a NULL pointer !")) ;
164 _fromMesh = _fromField->getSupport()->getMesh();
166 if ( ! _fromMesh ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromMesh is a NULL pointer !")) ;
168 _fromWrapper = new Meta_Wrapper<DIMENSION>(_fromMesh->getNumberOfNodes(),
169 const_cast<double *> (_fromMesh->getCoordinates(MED_FULL_INTERLACE)),
170 const_cast<CONNECTIVITY *> (_fromMesh->getConnectivityptr()),
171 const_cast<FIELD<double> *>(_fromField)
175 _toWrapper = new Meta_Wrapper<DIMENSION>(_toMesh->getNumberOfNodes(),
176 const_cast<double *> (_toMesh->getCoordinates(MED_FULL_INTERLACE))
180 _mapping = new Meta_Mapping<DIMENSION> (_fromWrapper);
186 template <int DIMENSION> INTERPOLATION<DIMENSION>::~INTERPOLATION()
188 if ( _fromWrapper ) delete _fromWrapper ;
189 if ( _toWrapper ) delete _toWrapper ;
190 if ( _mapping ) delete _mapping ;
193 template <int DIMENSION> int INTERPOLATION<DIMENSION>::getNearestNode( double * node ) {
195 const char * LOC = "INTERPOLATION::getNearestNode( double * node ) ";
199 if ( ! _mapping ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer !")) ;
201 return _mapping->Donne_dTree()->trouve_plus_proche_point(Wrapper_Noeud<DIMENSION > (node) );
207 template <int DIMENSION> int INTERPOLATION<DIMENSION>::getContainingCell ( double * node , int beginingCell, int flagIsConvexMesh ) {
209 const char * LOC = "INTERPOLATION::getContainingCell( double * node ) ";
213 if ( ! _mapping ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer !")) ;
215 _isConvexFromMesh=flagIsConvexMesh;
217 _mapping->Cree_Mapping(_toWrapper,_isConvexFromMesh);
219 return _mapping->Trouve_Maille_Contenant_Noeud(node,beginingCell,flagIsConvexMesh);
225 template <int DIMENSION> vector<int> INTERPOLATION<DIMENSION>::getMapping ( int flagIsConvexMesh ) {
227 const char * LOC = "INTERPOLATION::getMapping( ) ";
231 if ( ! _mapping ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer !")) ;
232 if ( ! _toWrapper ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toWrapper is a NULL pointer !")) ;
234 _isConvexFromMesh=flagIsConvexMesh;
236 _mapping->Cree_Mapping(_toWrapper,_isConvexFromMesh);
238 return _mapping->Get_Mapping();
244 template <int DIMENSION> FIELD<double> * INTERPOLATION<DIMENSION>::interpolate(int itype,int flagIsConvexFromMesh) {
246 const char * LOC = "INTERPOLATION::interpolate(int itype,int flagIsConvexFromMesh) ";
250 if ( ! _mapping ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer !")) ;
251 if ( ! _toWrapper ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toWrapper is a NULL pointer !")) ;
252 if ( ! _fromField ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromField is a NULL pointer !")) ;
256 _isConvexFromMesh=flagIsConvexFromMesh;
258 _mapping->Cree_Mapping(_toWrapper,_isConvexFromMesh);
260 Wrapper_Nuage_Noeud<DIMENSION> * toNodes = _toWrapper->Get_Nuage_Noeuds();
262 Wrapper_MED_Field resultat;
265 cout<<"Mapping"<<endl;
267 cout<<"Mailles"<<endl;
268 _fromWrapper->Get_Maillage()->DONNE_POINTEUR_NUAGEMAILLE()->affiche();
269 cout<<"Noeuds"<<endl;
270 _fromWrapper->Get_Nuage_Noeuds()->affiche();
275 case 0 : // INTERPOLATION P0
276 cout<<"Avant ="<<endl;
277 resultat=Meta_Interpolateur< Meta_Calcul_Interpolation_P0<DIMENSION>,DIMENSION >(_mapping,_fromWrapper).Perform_Interpolation(toNodes);
279 case 1 : // INTERPOLATION P-Hybride (Interpole avec la fonction d'interpolation naturelle de la maille contenant le point)
280 resultat=Meta_Interpolateur< Meta_Calcul_Interpolation_Hybride<DIMENSION>,DIMENSION >(_mapping,_fromWrapper).Perform_Interpolation(toNodes);
282 case 2 : // INTERPOLATION (P/Q) 1 forcée (Interpole avec la fonction élément fini de la maille de degré 1 -meme si la maille est de degré supérieur-)
283 resultat=Meta_Interpolateur< Meta_Calcul_Interpolation_Hybride_P1<DIMENSION>, DIMENSION >(_mapping,_fromWrapper).Perform_Interpolation(toNodes);
286 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Interpolation type "<<itype<<" not yet implemented !")) ;
289 cout<<"On a fini l'interpolation"<<endl;
291 _toField = new FIELD<double>;
293 _toField->setName ( _fromField->getName() );
294 _toField->setDescription ( _fromField->getDescription() );
295 _toField->setNumberOfComponents ( _fromField->getNumberOfComponents() );
296 _toField->setNumberOfValues ( _toMesh ->getNumberOfNodes() );
297 _toField->setComponentsNames ( _fromField->getComponentsNames() );
298 _toField->setComponentsDescriptions ( _fromField->getComponentsDescriptions() );
299 _toField->setMEDComponentsUnits ( _fromField->getMEDComponentsUnits() );
300 _toField->setIterationNumber ( _fromField->getIterationNumber() );
301 _toField->setTime ( _fromField->getTime() );
302 _toField->setOrderNumber ( _fromField->getOrderNumber() );
303 _toField->setValueType ( MED_EN::MED_REEL64 );
305 SUPPORT * mySupport(new SUPPORT(_toMesh,"support",MED_NODE));
306 _toField->setSupport(mySupport);
308 _toField->allocValue(_toField->getNumberOfComponents(),_toField->getNumberOfValues());
310 _toField->setValue(MED_FULL_INTERLACE,resultat.Get_Valeurs());
312 _toWrapper->Construit_Wrapper_Champ(_toField);
320 template <int DIMENSION> FIELD<double> * INTERPOLATION<DIMENSION>::interpolateNextStep(const FIELD <double> & nextFromField, int & flagNewMapping) {
322 const char * LOC = "INTERPOLATION::interpolateNextStep(int itype,int flagIsConvexFromMesh) ";
326 if ( ! _mapping ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer !" ));
327 if ( ! _toWrapper ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toWrapper is a NULL pointer !" )) ;
328 if ( ! _fromWrapper ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromWrapper is a NULL pointer !" )) ;
329 if ( ! _fromField ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromField is a NULL pointer !" )) ;
332 if ( ! _toField ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toField is a NULL pointer, wrong use of interpolateNextStep" )) ;
333 if ( _iType==UNDEFINED ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"_iType is not defined, wrong use of interpolateNextStep" )) ;
334 if ( _isConvexFromMesh==UNDEFINED ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"_isConvexFromMesh is not defined, wrong use of interpolateNextStep" )) ;
336 // delete _toField; ????????????????????????????
338 // if the mesh are identical, the mapping is the same, if not, the mapping has to be re-calculated
339 if (nextFromField.getSupport()->getMesh()->getName()!=_fromMesh->getName())
347 _fromField = const_cast<FIELD<double> *>(&nextFromField);
349 _fromMesh = _fromField->getSupport()->getMesh();
351 _fromWrapper = new Meta_Wrapper<DIMENSION>(_fromMesh->getNumberOfNodes(),
352 const_cast<double *> (_fromMesh->getCoordinates(MED_FULL_INTERLACE)),
353 const_cast<CONNECTIVITY *> (_fromMesh->getConnectivityptr()),
354 const_cast<FIELD<double> *>(_fromField)
356 _mapping = new Meta_Mapping<DIMENSION> (_fromWrapper);
358 _mapping->Cree_Mapping(_toWrapper,_isConvexFromMesh);
366 _fromField = const_cast<FIELD<double> *>(&nextFromField);
367 _fromWrapper->Change_Champ(const_cast<FIELD<double> *>(_fromField));
370 return interpolate(_iType,_isConvexFromMesh);