*/
/// @{
+ /**
+ * Return an intermediate point between two given
+ * points. Overloading this function allows the default recursive
+ * implementation of the method get_new_point() that takes a
+ * Quadrature object as input to work properly.
+ *
+ * An implementation of this function should returns a parametric
+ * curve on the manifold, joining the points `p1` and `p2`, with
+ * parameter `w` in the interval [0,1]. In particular
+ * `get_new_point(p1, p2, 0.0)` should return `p1` and
+ * `get_new_point(p1, p2, 1.0)` should return `p2`.
+ *
+ * In its default implementation, this function calls the
+ * project_to_manifold() method with the convex combination of `p1`
+ * and `p2`. User classes can get away by simply implementing the
+ * project_to_manifold() method.
+ */
+ virtual
+ Point<spacedim>
+ get_new_point( const Point<spacedim> &p1,
+ const Point<spacedim> &p2,
+ const double w) const;
+
/**
* Return the point which shall become the new vertex surrounded by the
* given points which make up the quadrature. We use a quadrature object,
* which should be filled with the surrounding points together with
* appropriate weights.
*
- * In its default implementation it calls internally the function
- * project_to_manifold. User classes can get away by simply implementing
+ * In its default implementation it calls recursively the function
+ * get_new_point(). User classes can get away by simply implementing
* that method.
*/
virtual
Manifold<dim, spacedim>::~Manifold ()
{}
-
-
template <int dim, int spacedim>
Point<spacedim>
Manifold<dim, spacedim>::
return Point<spacedim>();
}
-
+template <int dim, int spacedim>
+Point<spacedim>
+Manifold<dim, spacedim>::
+get_new_point (const Point<spacedim> &p1,
+ const Point<spacedim> &p2,
+ const double w) const
+{
+ std::vector<Point<spacedim> > vertices;
+ vertices.push_back(p1);
+ vertices.push_back(p2);
+ return project_to_manifold(vertices, w * p1 + (1-w)*p2 );
+}
template <int dim, int spacedim>
Point<spacedim>
Manifold<dim, spacedim>::
get_new_point (const Quadrature<spacedim> &quad) const
{
+ Assert(quad.size() > 0,
+ ExcMessage("Quadrature should have at least one point."));
+
Assert(std::abs(std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.0)-1.0) < 1e-10,
ExcMessage("The weights for the individual points should sum to 1!"));
- Point<spacedim> p;
- for (unsigned int i=0; i<quad.size(); ++i)
- p += quad.weight(i) * quad.point(i);
+ QSorted<spacedim> sorted_quad(quad);
+ Point<spacedim> p = sorted_quad.point(0);
+ double w = sorted_quad.weight(0);
+
+ for (unsigned int i=1; i<sorted_quad.size(); ++i)
+ {
+ if ( w != 0 )
+ p = get_new_point(p, sorted_quad.point(i) , w/(sorted_quad.weight(i) + w) );
+ w += sorted_quad.weight(i);
+ }
- return project_to_manifold(quad.get_points(), p);
+ return p;
}