get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
const std::vector<double> &weights) const = 0;
-
/**
- * Given a candidate point and a
- * line segment characterized by
- * the iterator, return a point
- * that lies on the surface
- * described by this object. This
- * function is used in some mesh
- * smoothing algorithms that try
- * to move around points in order
- * to improve the mesh quality
- * but need to ensure that points
- * that were on the boundary
- * remain on the boundary.
+ * Return equally spaced intermediate points on between the given
+ * surrounding points.
+ *
+ * The number of points requested is given by the size of the vector
+ * @p points. It is the task of the derived classes to arrange the
+ * points in approximately equal distances. If the surrounding
+ * points are only two, then the returned points are distributed
+ * along a line. If they are 4, and spacedim >= 2, they are
+ * distributed equally on the patch identified by the four points
+ * (which are thought to be the vertices of a quad), and an
+ * exception is thrown if the number of requested points is not the
+ * square of an integer number. The same is true if the number of
+ * surrounding points is 8: an exception is thrown if spacedim is
+ * not three, and the number of requested points is not the cube of
+ * an integer number.
*
- * If spacedim==1, then the line
- * represented by the line
- * iterator is the entire space
- * (i.e. it is a cell, not a part
- * of the boundary), and the
- * returned point equals the
- * given input point.
+ * This function is called by the @p MappingQ class, and it calls
+ * internally the get_new_point class with appropriate weights, so
+ * that the user does not need to overload anything else. It is made
+ * virtual so that the user can overload it if the default
+ * implementation is too slow.
+ */
+ virtual
+ void
+ get_intermediate_points(std::vector<Point<spacedim> > &points,
+ const std::vector<Point<spacedim> > &surrounding_points) const;
+
+ /**
+ * Given a candidate point and a line segment characterized by the
+ * iterator, return a point that lies on the surface described by
+ * this object. This function is used in some mesh smoothing
+ * algorithms that try to move around points in order to improve the
+ * mesh quality but need to ensure that points that were on the
+ * boundary remain on the boundary.
*
- * Derived classes do not need to
- * implement this function unless
- * mesh smoothing algorithms are
- * used with a particular
- * boundary object. The default
- * implementation of this
- * function throws an exception
- * of type ExcPureFunctionCalled.
+ * Derived classes do not need to implement this function unless
+ * mesh smoothing algorithms are used with a particular boundary
+ * object. The default implementation of this function throws an
+ * exception of type ExcPureFunctionCalled.
*/
virtual
void project_to_manifold (Point<spacedim> &candidate) const;
Assert (false, ExcPureFunctionCalled());
}
+template <int spacedim>
+void
+Manifold<spacedim>::get_intermediate_points(std::vector<Point<spacedim> > &points,
+ const std::vector<Point<spacedim> > &surrounding_points) const
+{
+ Assert(surrounding_points.size() >= 2, ExcMessage("At least 2 surrounding points are required"));
+ const unsigned int n=points.size();
+ Assert(n>0, ExcMessage("You can't ask for 0 intermediate points."));
+ std::vector<double> w(surrounding_points.size());
+
+ switch(surrounding_points.size())
+ {
+ case 2:
+ {
+ const double dx=1./(n+1);
+ double x = dx;
+ for(unsigned int i=0; i<n; ++i, x+=dx)
+ {
+ w[1] = x; w[0] = (1-x);
+ points[i] = get_new_point(surrounding_points, w);
+ }
+ }
+
+ break;
+ case 4:
+ {
+ Assert(spacedim >= 2, ExcImpossibleInDim(spacedim));
+ const unsigned m=
+ static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
+ // is n a square number
+ Assert(m*m==n, ExcInternalError());
+ const double ds=1./(m+1);
+ double y=ds;
+ for (unsigned int i=0; i<m; ++i, y+=ds)
+ {
+ double x=ds;
+ for (unsigned int j=0; j<m; ++j, x+=ds)
+ {
+ w[0] = (1-x); w[1] = x*(1-y);
+ w[2] = y*(1-x); w[3] = x*y;
+ points[i*m+j]=get_new_point(surrounding_points, w);
+ }
+ }
+ }
+
+ break;
+ case 8:
+ Assert(false, ExcNotImplemented());
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ break;
+ }
+}
+
/* -------------------------- FlatManifold --------------------- */