const std::vector<double> &weights) const;
/**
- * Compute a new set of points around the given points
+ * Compute a new set of points that interpolate between the given points
* @p surrounding_points. @p weights is a table with as many columns as
* @p surrounding_points.size(). The number of rows in @p weights determines
* how many new points will be computed and appended to the last input
* @p new_points equals the size at entry plus the number of rows in
* @p weights.
*
- * In its default implementation, this function simply calls
- * get_new_point() on each row of @weights and appends those points to the
- * output vector @p new_points. However, this function is more efficient if
- * multiple new points need to be generated like in MappingQGeneric and the
- * manifold does expensive transformations between a chart space and the
- * physical space, such as ChartManifold. If efficiency is not important,
- * you may get away by implementing only the get_new_point() function.
+ * In its default implementation, this function simply calls get_new_point()
+ * on each row of @weights and appends those points to the output vector
+ * @p new_points. However, this function is more efficient if multiple new
+ * points need to be generated like in MappingQGeneric and the manifold does
+ * expensive transformations between a chart space and the physical space,
+ * such as ChartManifold. For this function, the surrounding points need to
+ * be transformed back to the chart sparse only once, rather than for every
+ * call to get_new_point(). If efficiency is not important, you may get away
+ * by implementing only the get_new_point() function.
*
* The implementation does not allow for @p surrounding_points and
- * @p new_points to point to the same vector, so make sure pass different
+ * @p new_points to point to the same vector, so make sure to pass different
* objects into the function.
*/
virtual
const std::vector<double> &weights) const;
/**
- * Compute a new set of points around the given points
+ * Compute a new set of points that interpolate between the given points
* @p surrounding_points. @p weights is a table with as many columns as
* @p surrounding_points.size(). The number of rows in @p weights determines
* how many new points will be computed and appended to the last input
* @p new_points equals the size at entry plus the number of rows in
* @p weights.
*
- * For this particular implementation, an interpolation of the
- * @p surrounding_points according to the @p weights is performed.
+ * For this particular implementation, the interpolation of the
+ * @p surrounding_points according to the @p weights is simply performed in
+ * Cartesian space.
*/
virtual
void
const std::vector<double> &weights) const;
/**
- * Compute a new set of points around the given points
+ * Compute a new set of points that interpolate between the given points
* @p surrounding_points. @p weights is a table with as many columns as
* @p surrounding_points.size(). The number of rows in @p weights determines
* how many new points will be computed and appended to the last input
* @p weights.
*
* The implementation of this function first transforms the
- * @p surrounding_points to the chart by calling pull_back(). Then, new
+ * @p surrounding_points to the chart space by calling pull_back(). Then, new
* points are computed on the chart by usual interpolation according to the
* given @p weights, which are finally transformed to the image space by
* push_forward().
*
* This implementation can be much more efficient for computing multiple new
* points from the same surrounding points than separate calls to
- * get_new_point() in case the pull_back() operation is expensive. Often,
- * this is indeed the case because pull_back() might involve Newton
- * iterations or something similar in non-trivial manifolds.
+ * get_new_point() in case the pull_back() operation is expensive. This is
+ * because pull_back() is only called once for the surrounding points and
+ * the interpolation is done for all given weights using this set of
+ * points. Often, pull_back() is also more expensive than push_forward()
+ * because the former might involve some kind of Newton iteration in
+ * non-trivial manifolds.
*/
virtual
void
public:
CompareWeights(const std::vector<double> &weights)
:
- compare_weights(&weights)
+ compare_weights(weights)
{}
bool operator() (unsigned int a, unsigned int b) const
{
- return (*compare_weights)[a] < (*compare_weights)[b];
+ return compare_weights[a] < compare_weights[b];
}
private:
- const std::vector<double> *compare_weights;
+ const std::vector<double> &compare_weights;
};
/* -------------------------- Manifold --------------------- */
// First sort points in the order of their weights. This is done to
// produce unique points even if get_intermediate_points is not
// associative (as for the SphericalManifold).
- std::vector<unsigned int> permutation(n_points);
+ unsigned int permutation_short[30];
+ std::vector<unsigned int> permutation_long;
+ unsigned int *permutation;
+ if (n_points > 30)
+ {
+ permutation_long.resize(n_points);
+ permutation = &permutation_long[0];
+ }
+ else
+ permutation = &permutation_short[0];
+
for (unsigned int i=0; i<n_points; ++i)
permutation[i] = i;
- std::sort(permutation.begin(),
- permutation.end(),
+ std::sort(permutation,
+ permutation + n_points,
CompareWeights(weights));
// Now loop over points in the order of their associated weight
else
weight = w/(weights[permutation[i]] + w);
- p = get_intermediate_point(p, surrounding_points[permutation[i]],1.0 - weight );
+ if (std::abs(weight) > 1e-14)
+ p = get_intermediate_point(p, surrounding_points[permutation[i]],1.0 - weight );
w += weights[permutation[i]];
}
ExcMessage("surrounding_points and new_points cannot be the same "
"array"));
- std::vector<double> local_weights(surrounding_points.size());
+ const unsigned int n_points = surrounding_points.size();
+ std::vector<double> local_weights(n_points);
for (unsigned int row=0; row<weights.size(0); ++row)
{
- for (unsigned int p=0; p<surrounding_points.size(); ++p)
- local_weights[p] = weights[row][p];
+ for (unsigned int i=0; i<n_points; ++i)
+ local_weights[i] = weights(row,i);
new_points.push_back(get_new_point(surrounding_points, local_weights));
}
}
const unsigned int n_points = surrounding_points.size();
std::vector<Point<chartdim> > chart_points(n_points);
- std::vector<double> local_weights(n_points);
-
for (unsigned int i=0; i<n_points; ++i)
chart_points[i] = pull_back(surrounding_points[i]);
+ std::vector<Point<chartdim> > new_points_on_chart;
+ new_points_on_chart.reserve(weights.size(0));
+ sub_manifold.add_new_points(chart_points, weights, new_points_on_chart);
+
for (unsigned int row=0; row<weights.size(0); ++row)
- {
- for (unsigned int p=0; p<n_points; ++p)
- local_weights[p] = weights[row][p];
- const Point<chartdim> p_chart = sub_manifold.get_new_point(chart_points,
- local_weights);
- new_points.push_back(push_forward(p_chart));
- }
+ new_points.push_back(push_forward(new_points_on_chart[row]));
}
projects / dealii.git / commitdiff