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
|
// Copyright (C) 2016 EDF
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#include <vector>
#include <memory>
#include <Eigen/Dense>
#include "StOpt/core/grids/LegendreInterpolator.h"
#include "StOpt/core/grids/RegularLegendreGrid.h"
using namespace StOpt ;
using namespace std;
LegendreInterpolator::LegendreInterpolator(const RegularLegendreGrid *p_grid, const Eigen::ArrayXd &p_point)
{
shared_ptr< array< function< double(const double &) >, 11 > > legendre = p_grid->getLegendre();
shared_ptr< vector< Eigen::ArrayXXd > > fInterpol = p_grid->getFInterpol();
Eigen::ArrayXi poly = p_grid->getPoly();
// number of points involved
m_nbWeigth = (1 + poly).prod();
m_weightAndPoints.resize(m_nbWeigth);
// coordinate min
Eigen::ArrayXi coordmin = p_grid->lowerPositionCoord(p_point);
Eigen::ArrayXd meshSize = p_grid->getMeshSize(coordmin);
// get back real coordinates renormalized in [-1,1]
Eigen::ArrayXd xCoord(p_point.size());
Eigen::ArrayXd xCoordMin = p_grid->getCoordinateFromIntCoord(coordmin);
for (int ip = 0; ip < xCoord.size(); ++ip)
xCoord(ip) = std::max(std::min(2.*(p_point(ip) - xCoordMin(ip)) / meshSize(ip) - 1., 1.), -1.);
// calculate weights
Eigen::ArrayXi iCoord(p_point.size()) ; // coordinates in the mesh
// iterate on all points on the mesh
for (int j = 0 ; j < m_nbWeigth ; ++j)
{
int jloc = j;
int nPoint = m_nbWeigth;
for (int id = p_point.size() - 1 ; id >= 0 ; --id)
{
nPoint /= (poly(id) + 1);
iCoord(id) = jloc / nPoint;
jloc = jloc % nPoint;
}
// now iterates
// calculate \f$\sum_i ....\sum_j L_i(\xi_{icoord(0)}) L_i(xCoord) \rho_i \kappa_i ...L_j(\xi_{icoord(nd-1)}) L_j(xCoord(nd-1)) \rho_j (j+0.5) \f$
double weight = 0.;
for (int jj = 0 ; jj < m_nbWeigth ; ++jj)
{
int jjloc = jj;
int nnPoint = m_nbWeigth;
double weightLocal = 1.;
for (int id = p_point.size() - 1 ; id >= 0 ; --id)
{
nnPoint /= (poly(id) + 1);
int iiCoord = jjloc / nnPoint;
weightLocal *= (*fInterpol)[id](iiCoord, iCoord(id)) * (*legendre)[iiCoord](xCoord[id]);
jjloc = jjloc % nnPoint;
}
weight += weightLocal;
}
m_weightAndPoints(j) = make_pair(weight, p_grid->intCoordPerDimToGlobal(iCoord + coordmin));
}
}
|
