From: heltai Date: Sat, 18 Jan 2014 12:38:37 +0000 (+0000) Subject: Moved get_intermediate_points from Manifold to Mapping X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d12121cc8fa4d5c67ae7c94ac6d12c9a0313c8d7;p=dealii-svn.git Moved get_intermediate_points from Manifold to Mapping git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32244 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/fe/mapping_q.h b/deal.II/include/deal.II/fe/mapping_q.h index e43960e015..b5b7fccef7 100644 --- a/deal.II/include/deal.II/fe/mapping_q.h +++ b/deal.II/include/deal.II/fe/mapping_q.h @@ -265,6 +265,15 @@ protected: add_quad_support_points(const typename Triangulation::cell_iterator &cell, std::vector > &a) const; + /** + * Ask the manifold descriptor to return intermediate points on + * lines or faces. According to the dimension of the + * surrounding_points vector, intermediate points on lines or on + * quads are computed. + */ + void get_intermediate_points(const Manifold &manifold, + const std::vector > &surrounding_points, + std::vector > &points) const; private: virtual @@ -356,6 +365,7 @@ private: const typename Triangulation::cell_iterator &cell, std::vector > &a) const; + /** * Needed by the @p laplace_on_quad function (for dim==2). Filled * by the constructor. @@ -434,6 +444,16 @@ private: */ const FE_Q feq; + /* + * The default line support points. These are used when computing + * the location in real space of the support points on line and + * quads, which are queried to the Manifold class. This + * functionality was originally in the Boundary. We + * moved it here after the transition to Manfifold, since + * there was an hack in order to make it work there. + */ + QGaussLobatto<1> line_support_points; + /** * Declare other MappingQ classes friends. */ diff --git a/deal.II/include/deal.II/grid/manifold.h b/deal.II/include/deal.II/grid/manifold.h index 357b6bd560..e26dd7e4e5 100644 --- a/deal.II/include/deal.II/grid/manifold.h +++ b/deal.II/include/deal.II/grid/manifold.h @@ -84,34 +84,6 @@ public: get_new_point(const std::vector > &surrounding_points, const std::vector &weights) const = 0; - /** - * Return equally spaced intermediate points 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 eight: an exception is thrown if spacedim - * is not three, and the number of requested points is not the cube - * of an integer number. - * - * 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, in general, to overload it. It is - * made virtual so that the user can overload it if the default - * implementation is too slow. - */ - virtual - void - get_intermediate_points(const std::vector > &surrounding_points, - std::vector > &points) const; - /** * Given a point which lies close to the given manifold, it modifies * it and projects it to manifold itself. This function is used in @@ -147,34 +119,6 @@ public: */ virtual Point normal_vector(const Point &point) const; - - protected: - - /** - * Returns the support points of the Gauss-Lobatto quadrature - * formula used for intermediate points. - * - * @note Since the manifold description is closely tied to the unit - * cell support points of MappingQ, new boundary descriptions need - * to explicitly use these Gauss-Lobatto points and not equidistant - * points. - */ - const std::vector > & - get_line_support_points (const unsigned int n_intermediate_points) const; - - - private: - - /** - * Point generator for the intermediate points on a boundary. - */ - mutable std::vector > > points; - - - /** - * Mutex for protecting the points array. - */ - mutable Threads::Mutex mutex; }; diff --git a/deal.II/source/fe/mapping_c1.cc b/deal.II/source/fe/mapping_c1.cc index 9097bb47fd..b7e0ec446c 100644 --- a/deal.II/source/fe/mapping_c1.cc +++ b/deal.II/source/fe/mapping_c1.cc @@ -68,7 +68,7 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce // first get the normal vectors at the two vertices of this line // from the boundary description std::vector > face_vertex_normals(2); - line->get_boundary().get_normals_at_points (face_vertex_normals, points); + line->get_boundary().get_normals_at_points (points, face_vertex_normals); // then transform them into interpolation points for a cubic // polynomial @@ -137,7 +137,7 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce // not at boundary { static const StraightBoundary straight_boundary; - straight_boundary.get_intermediate_points (line_points, points); + get_intermediate_points (straight_boundary, points, line_points); a.insert (a.end(), line_points.begin(), line_points.end()); }; }; diff --git a/deal.II/source/fe/mapping_q.cc b/deal.II/source/fe/mapping_q.cc index 91e743e957..200d5b887b 100644 --- a/deal.II/source/fe/mapping_q.cc +++ b/deal.II/source/fe/mapping_q.cc @@ -71,7 +71,8 @@ MappingQ<1>::MappingQ (const unsigned int, n_shape_functions(2), renumber(0), use_mapping_q_on_all_cells (false), - feq(degree) + feq(degree), + line_support_points(degree+1) {} @@ -85,7 +86,8 @@ MappingQ<1>::MappingQ (const MappingQ<1> &m): n_shape_functions(2), renumber(0), use_mapping_q_on_all_cells (m.use_mapping_q_on_all_cells), - feq(degree) + feq(degree), + line_support_points(degree+1) {} template<> @@ -128,7 +130,8 @@ MappingQ::MappingQ (const unsigned int p, degree))), use_mapping_q_on_all_cells (use_mapping_q_on_all_cells || (dim != spacedim)), - feq(degree) + feq(degree), + line_support_points(degree+1) { // Construct the tensor product polynomials used as shape functions for the // Qp mapping of cells at the boundary. @@ -161,7 +164,8 @@ MappingQ::MappingQ (const MappingQ &mapping) n_shape_functions(mapping.n_shape_functions), renumber(mapping.renumber), use_mapping_q_on_all_cells (mapping.use_mapping_q_on_all_cells), - feq(degree) + feq(degree), + line_support_points(degree+1) { tensor_pols=new TensorProductPolynomials (*mapping.tensor_pols); laplace_on_quad_vector=mapping.laplace_on_quad_vector; @@ -785,16 +789,11 @@ MappingQ::add_line_support_points (const typename Triangulation::lines_per_cell; ++line_no) { - const typename Triangulation::line_iterator - line = (dim == 1 - ? - static_cast::line_iterator>(cell) - : - cell->line(line_no)); + const typename Triangulation::line_iterator line = cell->line(line_no); ps[0] = line->vertex(0); ps[1] = line->vertex(1); - line->get_manifold().get_intermediate_points (ps, line_points); + get_intermediate_points (line->get_manifold(), ps, line_points); if (dim==3) { @@ -873,7 +872,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, vertices.resize(4); for(unsigned int i=0;i<4; ++i) vertices[i] = face->vertex(i); - face->get_manifold().get_intermediate_points(vertices, quad_points); + get_intermediate_points(face->get_manifold(), vertices, quad_points); // in 3D, the orientation, flip and rotation of the face might not // match what we expect here, namely the standard orientation. thus // reorder points accordingly. since a Mapping uses the same shape @@ -937,7 +936,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, vertices.resize(4); for(unsigned int i=0;i<4; ++i) vertices[i] = face->vertex(i); - face->get_manifold().get_intermediate_points (vertices, quad_points); + get_intermediate_points (face->get_manifold(), vertices, quad_points); // in 3D, the orientation, flip and rotation of the face might // not match what we expect here, namely the standard // orientation. thus reorder points accordingly. since a Mapping @@ -965,8 +964,7 @@ add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell, std::vector > vertices(4); for(unsigned int i=0; i<4; ++i) vertices[i] = cell->vertex(i); - - cell->get_manifold().get_intermediate_points (vertices, quad_points); + get_intermediate_points (cell->get_manifold(), vertices, quad_points); for (unsigned int i=0; i::clone () const +template +void +MappingQ::get_intermediate_points (const Manifold &manifold, + const std::vector > &surrounding_points, + std::vector > &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: + { + // If two points are passed, these are the two vertices, and + // we can only compute degree-1 intermediate points. + AssertDimension(n, degree-1); + 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()); + + // If four points are passed, these are the two vertices, and + // we can only compute (degree-1)*(degree-1) intermediate + // points. + AssertDimension(m, degree-1); + + for (unsigned int i=0; i::normal_vector(const Point &point) const return point; } -template -void -Manifold::get_intermediate_points(const std::vector > &surrounding_points, - std::vector > &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: - { - // Use interior points of QGaussLobatto quadrature formula - // support points for consistency with MappingQ - const std::vector > &line_points = this->get_line_support_points(n); - - 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 std::vector > &line_points = this->get_line_support_points(m); - - for (unsigned int i=0; i -const std::vector > & -Manifold:: -get_line_support_points (const unsigned int n_intermediate_points) const -{ - if (points.size() <= n_intermediate_points || - points[n_intermediate_points].get() == 0) - { - Threads::Mutex::ScopedLock lock(mutex); - if (points.size() <= n_intermediate_points) - points.resize(n_intermediate_points+1); - - // another thread might have created points in the meantime - if (points[n_intermediate_points].get() == 0) - { - std_cxx1x::shared_ptr > - quadrature (new QGaussLobatto<1>(n_intermediate_points+2)); - points[n_intermediate_points] = quadrature; - } - } - return points[n_intermediate_points]->get_points(); -} - - /* -------------------------- FlatManifold --------------------- */