33 template <
class triangulationType,
class xDataType,
class yDataType>
35 const triangulationType *triangulation,
36 std::vector<std::tuple<int, double, double>> &surfaceValues,
37 std::array<const long, 2> &surfaceDim,
38 const xDataType *
const tableXValues,
39 const yDataType *
const tableYValues,
40 const size_t nTableValues,
41 std::vector<std::vector<double>> &inputPoints) {
42 unsigned int noPoints = nTableValues;
43 inputPoints = std::vector<std::vector<double>>(
44 noPoints, std::vector<double>(3, 0.0));
50 std::vector<std::array<int, 4>> quadPoints(noPoints);
51 for(
unsigned int i = 0; i < noPoints; ++i) {
53 for(
unsigned int j = 0; j < surfaceDim[0]; ++j)
54 if(std::get<1>(surfaceValues[j * surfaceDim[1]]) < tableXValues[i])
56 for(
unsigned int j = 0; j < surfaceDim[1]; ++j)
57 if(std::get<2>(surfaceValues[j]) < tableYValues[i])
59 quadPoints[i][0] = i0 * surfaceDim[1] + i1 + 1;
60 quadPoints[i][1] = i0 * surfaceDim[1] + i1;
61 quadPoints[i][2] = (i0 + 1) * surfaceDim[1] + i1 + 1;
62 quadPoints[i][3] = (i0 + 1) * surfaceDim[1] + i1;
66 std::vector<std::array<double, 4>> coef(noPoints);
67 for(
unsigned int i = 0; i < noPoints; ++i) {
68 std::vector<std::array<double, 2>> points(4);
69 for(
unsigned int j = 0; j < 4; ++j) {
70 points[j][0] = std::get<1>(surfaceValues[quadPoints[i][j]]);
71 points[j][1] = std::get<2>(surfaceValues[quadPoints[i][j]]);
75 std::array<double, 3> tableValues{
static_cast<double>(tableXValues[i]),
76 static_cast<double>(tableYValues[i]),
78 std::array<std::array<double, 3>, 3> trianglePoints;
79 std::array<int, 3> indexes;
80 std::array<double, 2> mid{(points[3][0] - points[1][0]) / 2.0,
81 (points[0][1] - points[1][1]) / 2.0};
82 if(tableXValues[i] > points[1][0] + mid[0]) {
85 if(tableYValues[i] > points[1][1] + mid[1])
92 if(tableYValues[i] > points[1][1] + mid[1])
98 for(
unsigned int j = 0; j < 3; ++j) {
99 for(
unsigned int k = 0; k < 2; ++k)
100 trianglePoints[j][k] = points[indexes[j]][k];
101 trianglePoints[j][2] = 0;
105 tableValues.data(), trianglePoints[0].data(),
106 trianglePoints[1].data(), coef[i][indexes[2]]);
108 tableValues.data(), trianglePoints[0].data(),
109 trianglePoints[2].data(), coef[i][indexes[1]]);
111 tableValues.data(), trianglePoints[1].data(),
112 trianglePoints[2].data(), coef[i][indexes[0]]);
115 = coef[i][indexes[2]] + coef[i][indexes[1]] + coef[i][indexes[0]];
116 for(
unsigned int j = 0; j < 3; ++j)
117 coef[i][indexes[j]] /= sumArea;
120 for(
unsigned int k = 0; k < 4; ++k) {
121 float surfacePoint[3];
122 triangulation->getVertexPoint(
123 std::get<0>(surfaceValues[quadPoints[i][k]]), surfacePoint[0],
124 surfacePoint[1], surfacePoint[2]);
125 for(
unsigned int j = 0; j < 3; ++j)
126 inputPoints[i][j] += coef[i][k] * surfacePoint[j];
Minimalist debugging class.
void computeInputPoints(const triangulationType *triangulation, std::vector< std::tuple< int, double, double > > &surfaceValues, std::array< const long, 2 > &surfaceDim, const xDataType *const tableXValues, const yDataType *const tableYValues, const size_t nTableValues, std::vector< std::vector< double > > &inputPoints)
int computeTriangleArea(const T *p0, const T *p1, const T *p2, T &area)