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
117
118
119
120
121
122
123
124
125
126
|
// Copyright (C) 2016 EDF
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#include <vector>
#include <array>
#include <iostream>
#include <Eigen/Dense>
#include "StOpt/core/utils/types.h"
#include "StOpt/core/grids/InterpolatorSpectral.h"
using namespace std;
using namespace Eigen;
namespace StOpt
{
// Only upper part is filled in
ArrayXd localConstMatrixCalculation(const ArrayXi &p_simToCell,
const int &p_nbCell)
{
int nbSimul = p_simToCell.size();
// initialization
ArrayXd matReg = ArrayXd::Zero(p_nbCell);
for (int is = 0; is < nbSimul ; ++is)
{
// cell number
int ncell = p_simToCell(is) ;
matReg(ncell) += 1;
}
// check if matrix is null
for (int iReg = 0; iReg < matReg.size(); ++iReg)
if (matReg(iReg) < 0.5) // here use constant because matReg is integer
matReg(iReg) = 1;
// normalization
matReg /= nbSimul ;
return matReg;
}
ArrayXXd localConstSecondMemberCalculation(const ArrayXi &p_simToCell,
const int &p_nbCell,
const ArrayXXd &p_fToRegress)
{
int nbSimul = p_simToCell.size();
// number of function to regress
int iSecMem = p_fToRegress.rows();
ArrayXXd secMember = ArrayXXd::Zero(p_fToRegress.rows(), p_nbCell);
for (int is = 0; is < nbSimul ; ++is)
{
int nCell = p_simToCell(is);
// nest on second members
for (int nsm = 0 ; nsm < iSecMem ; ++nsm)
{
double xtemp = p_fToRegress(nsm, is) ;
// second member of the regression problem
secMember(nsm, nCell) += xtemp ;
}
}
// normalization
secMember /= nbSimul;
return secMember;
}
ArrayXXd localConstReconstruction(const ArrayXi &p_simToCell,
const ArrayXXd &p_foncBasisCoef)
{
int nbSimul = p_simToCell.size();
// initialization
ArrayXXd solution(p_foncBasisCoef.rows(), nbSimul) ;
for (int is = 0; is < nbSimul ; ++is)
{
int nCell = p_simToCell(is) ;
for (int isecMem = 0; isecMem < p_foncBasisCoef.rows(); ++isecMem)
solution(isecMem, is) = p_foncBasisCoef(isecMem, nCell) ;
}
return solution;
}
ArrayXd localConstReconstructionOnePoint(const ArrayXd &p_oneParticle,
const vector< shared_ptr< ArrayXd > > &p_mesh1D,
const ArrayXXd &p_foncBasisCoef)
{
// Values of the functon basis and position of the particle in the mesh
int iCell = 0 ;
int idecCell = 1;
for (int id = 0 ; id < p_oneParticle.size() ; id++)
{
int iMesh = 1 ;
while ((p_oneParticle(id) > (*p_mesh1D[id])(iMesh)) && (iMesh < p_mesh1D[id]->size() - 1)) iMesh++;
iCell += (iMesh - 1) * idecCell;
idecCell *= p_mesh1D[id]->size() - 1;
}
// reconstruction
ArrayXd solution(p_foncBasisCoef.rows());
for (int isecMem = 0; isecMem < solution.size(); ++isecMem)
solution(isecMem) = p_foncBasisCoef(isecMem, iCell) ;
return solution;
}
double localConstReconsOnePointSimStock(const ArrayXd &p_oneParticle, const ArrayXd &p_stock,
const std::vector< std::shared_ptr<InterpolatorSpectral> > &p_interpBaseFunc,
const std::vector< std::shared_ptr< ArrayXd > > &p_mesh1D)
{
// Values of the functon basis and position of the particle in the mesh
int iCell = 0 ;
int idecCell = 1;
for (int id = 0 ; id < p_oneParticle.size() ; id++)
{
int iMesh = 1 ;
while ((p_oneParticle(id) > (*p_mesh1D[id])(iMesh)) && (iMesh < p_mesh1D[id]->size() - 1)) iMesh++;
iCell += (iMesh - 1) * idecCell;
idecCell *= p_mesh1D[id]->size() - 1;
}
// reconstruction
return p_interpBaseFunc[iCell]->apply(p_stock);
}
}
|
