33 #include <gsl/gsl_sf_pow_int.h>
34 #include <gsl/gsl_sf_gamma.h>
46 namespace interpolation {
51 std::vector<double> getDeltaPos(
Mesh* mesh) {
54 for (
size_t i = 0; i < it.
getState().size(); ++i)
57 for (
size_t i = 0; i < pos2.size(); ++i)
64 std::vector<std::vector<double> > values,
67 : polyPatchOrder_m(polyPatchOrder),
68 smoothingOrder_m(smoothingOrder),
83 "PPSolveFactory::Solve",
84 "NULL Mesh input to Solve");
88 ss <<
"Mismatch between Mesh and values in Solve. Mesh size was "
89 << points->
end().
toInteger() <<
" but field values size was "
92 "PPSolveFactory::Solve",
95 if (smoothingOrder < polyPatchOrder) {
97 "PPSolveFactory::Solve",
98 "Polynomial order must be <= smoothing order in Solve");
103 "PPSolveFactory::Solve",
104 "Dual of Mesh was NULL");
108 edgePoints_m = std::vector< std::vector< std::vector<int> > >(posDim+1);
109 for (
int i = 1; i <= posDim; ++i)
111 (i, 0, smoothingOrder-polyPatchOrder);
119 std::vector<int> state = outOfBoundsIt.
getState();
120 std::vector<int> delta = std::vector<int>(state.size(), 2);
124 for (
size_t i = 0; i < state.size(); ++i)
127 }
else if (state[i] > end[i]) {
132 for (
size_t i = 0; i < state.size(); ++i)
134 position[i] = (pos2[i] - pos1[i])*(state[i]-1) + position[i];
135 }
else if (state[i] > end[i]) {
136 position[i] = (pos2[i] - pos1[i])*(state[i]-end[i]) + position[i];
143 std::vector< std::vector<int> > dataPoints =
145 thisPoints_m = std::vector< std::vector<double> >(dataPoints.size());
146 std::vector<double> deltaPos = getDeltaPos(
points_m);
149 for (
size_t i = 0; i < dataPoints.size(); ++i) {
151 for (
int j = 0; j < posDim; ++j)
152 thisPoints_m[i][j] = (0.5-dataPoints[i][j])*deltaPos[j];
159 std::vector< std::vector<int> > dataPoints =
161 size_t dataPointsSize = dataPoints.size();
165 for (
size_t i = 0; i < dataPointsSize; ++i) {
166 std::vector<int> itState = it.
getState();
167 for (
int j = 0; j < posDim; ++j)
170 bool outOfBounds =
false;
171 for (
int j = 0; j < posDim; ++j) {
172 itCurrent[j] -= dataPoints[i][j];
173 outOfBounds = outOfBounds ||
175 itCurrent[j] > end[j];
195 std::vector<double> deltaPos = getDeltaPos(
points_m);
200 std::vector<int> equalAxes;
203 equalAxes.push_back(j);
204 std::vector<std::vector<int> > edgePoints =
207 for (
size_t j = 0; j < edgePoints.size(); ++j) {
208 derivIndices_m.push_back(std::vector<int>(posDim));
215 for (
size_t k = 0; k < edgePoints[j].size(); ++k) {
216 derivIndices_m.back()[equalAxes[k]] = edgePoints[j][k];
223 std::vector<int> polyIndex =
225 for (
int k = 0; k < posDim && isEqual; ++k) {
226 isEqual = polyIndex[k] == derivIndices_m.back()[k];
242 for (
int j = 0; j < nPolyCoeffs; ++j) {
243 std::vector<int> pVec =
245 std::vector<int> qVec = derivIndices_m[i];
246 for (
int l = 0; l < posDim; ++l) {
247 int pow = pVec[l]-qVec[l];
249 derivPolyVec_m.back()(j+1) = 0.;
252 derivPolyVec_m.back()(j+1) *= gsl_sf_fact(pVec[l])/gsl_sf_fact(pow)*gsl_sf_pow_int(-deltaPos[l]/2., pow);
269 std::vector<double> point = std::vector<double>(posDim, 0.);
271 for (
int j = 0; j < posDim; ++j) {
274 for (
int j = 0; j < posDim; ++j) {
276 if (nearest[j] > end[j]) {
288 for(
int j = 0; j < posDim; ++j) {
299 polynomials_m = std::vector<SquarePolynomialVector*>(meshSize, NULL);
311 double oldPercentage = 0.;
315 double newPercentage = (total-it.toInteger())/
double(total)*100.;
316 if (newPercentage - oldPercentage > 10.) {
317 *gmsg <<
" Done " << newPercentage <<
" %" <<
endl;
318 oldPercentage = newPercentage;
331 std::vector<int> check,
334 std::vector<std::vector<int> >& nearbyPoints) {
335 check[checkIndex] = polyPower;
336 nearbyPoints.push_back(check);
337 if (checkIndex+1 == check.size())
339 for (
int i = 1; i < int(polyPower); ++i)
341 for (
int i = 0; i <= int(polyPower); ++i) {
342 check[checkIndex] = i;
350 int polyOrderUpper) {
353 "Point dimension must be > 0"));
354 if (polyOrderLower > polyOrderUpper)
356 "Polynomial lower bound must be <= polynomial upper bound"));
358 if (polyOrderUpper == polyOrderLower || polyOrderUpper < 0)
359 return std::vector<std::vector<int> >();
361 std::vector<std::vector<int> > nearbyPoints(1, std::vector<int>(pointDim, 0));
363 for (
size_t thisPolyOrder = 0;
364 thisPolyOrder < size_t(polyOrderUpper);
367 if (polyOrderLower < 0)
369 int nPointsLower = 1;
370 for (
int dim = 0; dim < pointDim; ++dim)
371 nPointsLower *= polyOrderLower+1;
372 nearbyPoints.erase(nearbyPoints.begin(), nearbyPoints.begin()+nPointsLower);
static std::vector< int > IndexByPower(int index, int nInputVariables)
std::vector< double > outOfBoundsPosition(Mesh::Iterator outOfBoundsIt)
static void nearbyPointsRecursive(std::vector< int > check, size_t checkIndex, size_t polyPower, std::vector< std::vector< int > > &nearbyPoints)
std::vector< std::vector< double > > thisPoints_m
PolynomialCoefficient represents a coefficient in a multi-dimensional polynomial. ...
PolynomialCoefficient Coeff
virtual Mesh::Iterator end() const =0
virtual Mesh * dual() const =0
void getDerivs(Mesh::Iterator it)
void getValues(Mesh::Iterator it)
PPSolveFactory(Mesh *points, std::vector< std::vector< double > > values, int polyPatchOrder, int smoothingOrder)
std::vector< std::vector< int > > derivIndices_m
PolynomialPatch * solve()
std::vector< int > derivIndexByPower_m
SolveFactory is a factory class for solving a set of linear equations to generate a polynomial based ...
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
virtual Mesh::Iterator begin() const =0
Base class for meshing routines.
std::vector< std::vector< double > > values_m
std::vector< std::vector< double > > derivValues_m
static std::vector< std::vector< int > > getNearbyPointsSquares(int pointDim, int polyOrderLower, int polyOrderUpper)
std::vector< SquarePolynomialVector * > polynomials_m
virtual void getPosition(double *point) const
std::vector< int > getState() const
Patches together many SquarePolynomialVectors to make a multidimensional polynomial spline...
std::vector< std::vector< double > > derivPoints_m
std::vector< std::vector< double > > thisValues_m
std::vector< std::vector< int > > smoothingPoints_m
std::vector< MVector< double > > derivPolyVec_m
const Mesh * getMesh() const
virtual int getPositionDimension() const =0
std::vector< std::vector< std::vector< int > > > edgePoints_m
Inform & endl(Inform &inf)
static unsigned int NumberOfPolynomialCoefficients(int pointDimension, int order)
std::vector< std::vector< int > > derivOrigins_m