1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
|
// Copyright (C) 2016 EDF
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#ifndef SPARSEGRIDHIERARONEPOINTLINNOBOUND_H
#define SPARSEGRIDHIERARONEPOINTLINNOBOUND_H
#include <Eigen/Dense>
#include "StOpt/core/sparse/sparseGridTypes.h"
#include "StOpt/core/sparse/sparseGridUtils.h"
#include "StOpt/core/sparse/SparseGridHierarOnePointNoBound.h"
/** \file SparseGridHierarOnePointLinNoBound.h
* \brief Class to Hierarchize a single point with data structure without boundary points
* \author Xavier Warin
*/
namespace StOpt
{
/// \class SparseGridHierarOnePointLinNoBound SparseGridHierarOnePointLinNoBound.h
/// It permits to hierarchize a single point performing Hierarchization in the successive
/// dimensions : linear case, no boundary points, approximation of boundary conditions
template< class T, class TT>
class SparseGridHierarOnePointLinNoBound : public SparseGridHierarOnePointNoBound
{
/// \brief Recursive hierarchization
/// \param p_levelCurrent index level for the current point
/// \param p_positionCurrent position of the current point
/// \param p_iterLevel iterator on the level of the current point in the data structure
/// \param p_idim current working dimension
/// \param p_dataSet sparse grid data structure
/// \param p_nodalValues Array of nodal values depending on point number. Depending on TT, can be a two dimension array. Then each column corresponds to a point.
T recursive1DHierarchization(Eigen::ArrayXc &p_levelCurrent,
Eigen::ArrayXui &p_positionCurrent,
const SparseSet::const_iterator &p_iterLevel,
const int &p_idim,
const SparseSet &p_dataSet,
const TT &p_nodalValues)
{
// position
SparseLevel::const_iterator iterPosition = p_iterLevel->second.find(p_positionCurrent);
if (p_idim == -1)
{
// position
int iposPoint = iterPosition->second;
return DoubleOrArray()(p_nodalValues, iposPoint);
}
else
{
T ret = recursive1DHierarchization(p_levelCurrent, p_positionCurrent, p_iterLevel, p_idim - 1, p_dataSet, p_nodalValues);
if (p_levelCurrent(p_idim) > 1)
{
// store current configuration
char currentLevel = p_levelCurrent(p_idim);
unsigned int currentPosition = p_positionCurrent(p_idim);
// switch to direct father
GetDirectFatherNoBound()(p_levelCurrent(p_idim), p_positionCurrent(p_idim));
// find level iterator
SparseSet::const_iterator iterLevelDirectFather = p_dataSet.find(p_levelCurrent);
T hierarDirect = recursive1DHierarchization(p_levelCurrent, p_positionCurrent, iterLevelDirectFather, p_idim - 1, p_dataSet, p_nodalValues);
if (currentLevel == 2)
ret -= hierarDirect;
else
{
// not boundary point
if ((currentPosition != 0) && (currentPosition != lastNode[currentLevel - 1]))
{
// non direct father level and position
p_levelCurrent(p_idim) = currentLevel;
p_positionCurrent(p_idim) = currentPosition;
GetNonDirectFatherNoBound()(p_levelCurrent(p_idim), p_positionCurrent(p_idim));
// find level iterator
SparseSet::const_iterator iterLevelNonDirectFather = p_dataSet.find(p_levelCurrent);
T hierarNonDirect = recursive1DHierarchization(p_levelCurrent, p_positionCurrent, iterLevelNonDirectFather, p_idim - 1, p_dataSet, p_nodalValues);
ret -= 0.5 * (hierarDirect + hierarNonDirect);
}
else
{
// use extrapolation with grandfather
GetDirectFatherNoBound()(p_levelCurrent(p_idim), p_positionCurrent(p_idim));
SparseSet::const_iterator iterLevelDirectGrandFather = p_dataSet.find(p_levelCurrent);
T hierarGrandDirect = recursive1DHierarchization(p_levelCurrent, p_positionCurrent, iterLevelDirectGrandFather, p_idim - 1, p_dataSet, p_nodalValues);
T hierarNonDirect = 2 * hierarDirect - hierarGrandDirect ;
ret -= 0.5 * (hierarDirect + hierarNonDirect);
}
}
// go back to initial state
p_levelCurrent(p_idim) = currentLevel;
p_positionCurrent(p_idim) = currentPosition;
}
return ret ;
}
}
public :
/// \brief Hierarchize a single point (existing in the data structure)
/// \param p_levelCurrent index level for the current point
/// \param p_positionCurrent position of the current point
/// \param p_dataSet sparse grid data structure
/// \param p_nodalValues Array of nodal values depending on point number. Depending on TT, can be a two dimension array. Then each column corresponds to a point.
T operator()(const Eigen::ArrayXc &p_levelCurrent,
const Eigen::ArrayXui &p_positionCurrent,
const SparseSet &p_dataSet,
const TT &p_nodalValues)
{
// get iterator for the level in the data Structure
SparseSet::const_iterator iterLevel = p_dataSet.find(p_levelCurrent);
// copy level and position
Eigen::ArrayXc levelCurrent = p_levelCurrent;
Eigen::ArrayXui positionCurrent = p_positionCurrent;
return recursive1DHierarchization(levelCurrent, positionCurrent, iterLevel, p_levelCurrent.size() - 1, p_dataSet, p_nodalValues);
}
};
}
#endif /* SPARSEGRIDHIERARONEPOINTLINNOBOUND_H */
|
