From 60ba06a28bbc9587e9a16944c2ea90151c96415a Mon Sep 17 00:00:00 2001 From: heltai Date: Sun, 25 Aug 2013 11:35:07 +0000 Subject: [PATCH] Added get_intermediate_points interface to Manifold. git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@30478 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/grid/manifold.h | 65 ++++++++++++++----------- deal.II/source/grid/manifold.cc | 55 +++++++++++++++++++++ 2 files changed, 92 insertions(+), 28 deletions(-) diff --git a/deal.II/include/deal.II/grid/manifold.h b/deal.II/include/deal.II/grid/manifold.h index 77941d28cd..c55367e249 100644 --- a/deal.II/include/deal.II/grid/manifold.h +++ b/deal.II/include/deal.II/grid/manifold.h @@ -84,37 +84,46 @@ public: get_new_point(const std::vector > &surrounding_points, const std::vector &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 > &points, + const std::vector > &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 &candidate) const; diff --git a/deal.II/source/grid/manifold.cc b/deal.II/source/grid/manifold.cc index 10e3a51409..77e93e2d5d 100644 --- a/deal.II/source/grid/manifold.cc +++ b/deal.II/source/grid/manifold.cc @@ -39,6 +39,61 @@ Manifold::project_to_manifold (Point &candidate) const Assert (false, ExcPureFunctionCalled()); } +template +void +Manifold::get_intermediate_points(std::vector > &points, + const std::vector > &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 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= 2, ExcImpossibleInDim(spacedim)); + const unsigned m= + static_cast(std::sqrt(static_cast(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