From: Martin Kronbichler Date: Fri, 17 Mar 2017 14:20:39 +0000 (+0100) Subject: Add new function to query multiple points to Manifold X-Git-Tag: v8.5.0-rc1~17^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=834e865e44ee273361dbdf11bbce986f6be78407;p=dealii.git Add new function to query multiple points to Manifold --- diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h index 8e97d71509..424b331bf4 100644 --- a/include/deal.II/grid/manifold.h +++ b/include/deal.II/grid/manifold.h @@ -29,6 +29,9 @@ DEAL_II_NAMESPACE_OPEN +// forward declaration +template class Table; + /** * We collect here some helper functions used in the Manifold * classes. @@ -396,6 +399,33 @@ public: get_new_point (const std::vector > &surrounding_points, const std::vector &weights) const; + /** + * Compute a new set of points around 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 + * argument @p new_points. After exit of this function, the size of + * @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. + * + * The implementation does not allow for @p surrounding_points and + * @p new_points to point to the same vector, so make sure pass different + * objects into the function. + */ + virtual + void + add_new_points (const std::vector > &surrounding_points, + const Table<2,double> &weights, + std::vector > &new_points) const; + /** * Given a point which lies close to the given manifold, it modifies it and * projects it to manifold itself. @@ -723,6 +753,23 @@ public: get_new_point(const std::vector > &surrounding_points, const std::vector &weights) const; + /** + * Compute a new set of points around 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 + * argument @p new_points. After exit of this function, the size of + * @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. + */ + virtual + void + add_new_points (const std::vector > &surrounding_points, + const Table<2,double> &weights, + std::vector > &new_points) const; /** * Project to FlatManifold. This is the identity function for flat, @@ -934,6 +981,33 @@ public: get_new_point(const std::vector > &surrounding_points, const std::vector &weights) const; + /** + * Compute a new set of points around 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 + * argument @p new_points. After exit of this function, the size of + * @p new_points equals the size at entry plus the number of rows in + * @p weights. + * + * The implementation of this function first transforms the + * @p surrounding_points to the chart 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. + */ + virtual + void + add_new_points (const std::vector > &surrounding_points, + const Table<2,double> &weights, + std::vector > &new_points) const; + /** * Pull back the given point in spacedim to the Euclidean chartdim * dimensional space. diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index f36b8ca3ea..93a949cd00 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -14,6 +14,7 @@ // --------------------------------------------------------------------- #include +#include #include #include #include @@ -50,6 +51,8 @@ template Manifold::~Manifold () {} + + template Point Manifold:: @@ -60,6 +63,8 @@ project_to_manifold (const std::vector > &, return Point(); } + + template Point Manifold:: @@ -73,6 +78,8 @@ get_intermediate_point (const Point &p1, return project_to_manifold(vertices, w * p2 + (1-w)*p1 ); } + + template Point Manifold:: @@ -81,6 +88,8 @@ get_new_point (const Quadrature &quad) const return get_new_point(quad.get_points(),quad.get_weights()); } + + template Point Manifold:: @@ -131,6 +140,29 @@ get_new_point (const std::vector > &surrounding_points, +template +void +Manifold:: +add_new_points (const std::vector > &surrounding_points, + const Table<2,double> &weights, + std::vector > &new_points) const +{ + AssertDimension(surrounding_points.size(), weights.size(1)); + Assert(&surrounding_points != &new_points, + ExcMessage("surrounding_points and new_points cannot be the same " + "array")); + + std::vector local_weights(surrounding_points.size()); + for (unsigned int row=0; row Tensor<1,2> Manifold<2, 2>:: @@ -244,6 +276,7 @@ normal_vector (const Triangulation<3, 3>::face_iterator &face, } + template Tensor<1,spacedim> Manifold:: @@ -275,6 +308,7 @@ get_normals_at_vertices (const Triangulation<2, 2>::face_iterator &face, } + template <> void Manifold<3, 3>:: @@ -307,6 +341,7 @@ get_normals_at_vertices (const Triangulation<3, 3>::face_iterator &face, } + template void Manifold:: @@ -322,7 +357,6 @@ get_normals_at_vertices (const typename Triangulation::face_itera - template Point Manifold:: @@ -344,6 +378,7 @@ get_new_point_on_quad (const typename Triangulation::quad_iterato } + template Point Manifold:: @@ -383,6 +418,7 @@ get_new_point_on_cell (const typename Triangulation::cell_iterator } + template <> Point<1> Manifold<1,1>:: @@ -393,6 +429,7 @@ get_new_point_on_face (const Triangulation<1,1>::face_iterator &) const } + template <> Point<2> Manifold<1,2>:: @@ -414,6 +451,7 @@ get_new_point_on_face (const Triangulation<1,3>::face_iterator &) const } + template <> Point<1> Manifold<1,1>:: @@ -499,6 +537,8 @@ FlatManifold::FlatManifold (const Tensor<1,spacedim> &periodicity, tolerance(tolerance) {} + + template Point FlatManifold:: @@ -518,11 +558,69 @@ get_new_point (const std::vector > &surrounding_points, Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0)-1.0) < 1e-10, ExcMessage("The weights for the individual points should sum to 1!")); - Tensor<1,spacedim> minP = periodicity; + Point p; + // if there is no periodicity, use a shortcut + if (periodicity == Tensor<1,spacedim>()) + { + for (unsigned int i=0; i minP = periodicity; + + for (unsigned int d=0; d 0) + for (unsigned int i=0; i= -tolerance*periodicity[d]), + ExcPeriodicBox(d, surrounding_points[i], periodicity[i])); + } + + // compute the weighted average point, possibly taking into account periodicity + for (unsigned int i=0; i dp; + for (unsigned int d=0; d 0) + dp[d] = ( (surrounding_points[i][d]-minP[d]) > periodicity[d]/2.0 ? + -periodicity[d] : 0.0 ); + + p += (surrounding_points[i]+dp)*weights[i]; + } + + // if necessary, also adjust the weighted point by the periodicity + for (unsigned int d=0; d 0) + if (p[d] < 0) + p[d] += periodicity[d]; + } + + return project_to_manifold(surrounding_points, p); +} + + + +template +void +FlatManifold:: +add_new_points (const std::vector > &surrounding_points, + const Table<2,double> &weights, + std::vector > &new_points) const +{ + AssertDimension(surrounding_points.size(), weights.size(1)); + if (weights.size(0) == 0) + return; + + const unsigned int n_points = surrounding_points.size(); + + Tensor<1,spacedim> minP = periodicity; for (unsigned int d=0; d 0) - for (unsigned int i=0; i > &surrounding_points, ExcPeriodicBox(d, surrounding_points[i], periodicity[i])); } - // compute the weighted average point, possibly taking into account periodicity - Point p; - for (unsigned int i=0; i *surrounding_points_start = &surrounding_points[0]; + std::vector > modified_points; + bool adjust_periodicity = false; + for (unsigned int d=0; d 0) + for (unsigned int i=0; i periodicity[d]/2.0) + { + adjust_periodicity = true; + break; + } + if (adjust_periodicity == true) { - Point dp; + modified_points = surrounding_points; for (unsigned int d=0; d 0) - dp[d] = ( (surrounding_points[i][d]-minP[d]) > periodicity[d]/2.0 ? - -periodicity[d] : 0.0 ); - - p += (surrounding_points[i]+dp)*weights[i]; + for (unsigned int i=0; i periodicity[d]/2.0) + modified_points[i][d] -= periodicity[d]; + surrounding_points_start = &modified_points[0]; } - // if necessary, also adjust the weighted point by the periodicity - for (unsigned int d=0; d 0) - if (p[d] < 0) - p[d] += periodicity[d]; + // Now perform the interpolation + for (unsigned int row=0; row new_point; + for (unsigned int p=0; p 0) + if (new_point[d] < 0) + new_point[d] += periodicity[d]; + + new_points.push_back(project_to_manifold(surrounding_points, new_point)); + } } @@ -642,6 +761,36 @@ get_new_point (const std::vector > &surrounding_points, +template +void +ChartManifold:: +add_new_points (const std::vector > &surrounding_points, + const Table<2,double> &weights, + std::vector > &new_points) const +{ + Assert(weights.size(0) > 0, ExcEmptyObject()); + AssertDimension(surrounding_points.size(), weights.size(1)); + + const unsigned int n_points = surrounding_points.size(); + + std::vector > chart_points(n_points); + std::vector local_weights(n_points); + + for (unsigned int i=0; i p_chart = sub_manifold.get_new_point(chart_points, + local_weights); + new_points.push_back(push_forward(p_chart)); + } +} + + + template DerivativeForm<1,chartdim,spacedim> ChartManifold:: @@ -695,4 +844,3 @@ ChartManifold::get_periodicity() const #include "manifold.inst" DEAL_II_NAMESPACE_CLOSE -