add_quad_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
std::vector<Point<spacedim> > &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<spacedim> &manifold,
+ const std::vector<Point<spacedim> > &surrounding_points,
+ std::vector<Point<spacedim> > &points) const;
private:
virtual
const typename Triangulation<dim,spacedim>::cell_iterator &cell,
std::vector<Point<spacedim> > &a) const;
+
/**
* Needed by the @p laplace_on_quad function (for <tt>dim==2</tt>). Filled
* by the constructor.
*/
const FE_Q<dim> 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<spacedim> class. This
+ * functionality was originally in the Boundary<dim,spacedim>. We
+ * moved it here after the transition to Manfifold<spacedim>, since
+ * there was an hack in order to make it work there.
+ */
+ QGaussLobatto<1> line_support_points;
+
/**
* Declare other MappingQ classes friends.
*/
get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
const std::vector<double> &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<Point<spacedim> > &surrounding_points,
- std::vector<Point<spacedim> > &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
*/
virtual
Point<spacedim> normal_vector(const Point<spacedim> &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<Point<1> > &
- get_line_support_points (const unsigned int n_intermediate_points) const;
-
-
- private:
-
- /**
- * Point generator for the intermediate points on a boundary.
- */
- mutable std::vector<std_cxx1x::shared_ptr<QGaussLobatto<1> > > points;
-
-
- /**
- * Mutex for protecting the points array.
- */
- mutable Threads::Mutex mutex;
};
// first get the normal vectors at the two vertices of this line
// from the boundary description
std::vector<Point<dim> > 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
// not at boundary
{
static const StraightBoundary<dim> 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());
};
};
n_shape_functions(2),
renumber(0),
use_mapping_q_on_all_cells (false),
- feq(degree)
+ feq(degree),
+ line_support_points(degree+1)
{}
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<>
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.
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<dim> (*mapping.tensor_pols);
laplace_on_quad_vector=mapping.laplace_on_quad_vector;
// get the boundary description and second compute the points on it
for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
{
- const typename Triangulation<dim,spacedim>::line_iterator
- line = (dim == 1
- ?
- static_cast<typename Triangulation<dim,spacedim>::line_iterator>(cell)
- :
- cell->line(line_no));
+ const typename Triangulation<dim,spacedim>::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)
{
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
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
std::vector<Point<3> > 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<quad_points.size(); ++i)
a.push_back(quad_points[i]);
}
+template<int dim, int spacedim>
+void
+MappingQ<dim,spacedim>::get_intermediate_points (const Manifold<spacedim> &manifold,
+ const std::vector<Point<spacedim> > &surrounding_points,
+ std::vector<Point<spacedim> > &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:
+ {
+ // 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<n; ++i)
+ {
+ const double x = line_support_points.point(i+1)[0];
+ w[1] = x; w[0] = (1-x);
+ points[i] = manifold.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());
+
+ // 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<m; ++i)
+ {
+ const double y=line_support_points.point(1+i)[0];
+ for (unsigned int j=0; j<m; ++j)
+ {
+ const double x=line_support_points.point(1+j)[0];
+
+ w[0] = (1-x)*(1-y);
+ w[1] = x*(1-y);
+ w[2] = (1-x)*y ;
+ w[3] = x*y ;
+ points[i*m+j]=manifold.get_new_point(surrounding_points, w);
+ }
+ }
+ }
+
+ break;
+ case 8:
+ Assert(false, ExcNotImplemented());
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ break;
+ }
+}
+
// explicit instantiations
#include "mapping_q.inst"
return point;
}
-template <int spacedim>
-void
-Manifold<spacedim>::get_intermediate_points(const std::vector<Point<spacedim> > &surrounding_points,
- std::vector<Point<spacedim> > &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:
- {
- // Use interior points of QGaussLobatto quadrature formula
- // support points for consistency with MappingQ
- const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
-
- for(unsigned int i=0; i<n; ++i)
- {
- const double x = line_points[i+1][0];
- 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 std::vector<Point<1> > &line_points = this->get_line_support_points(m);
-
- for (unsigned int i=0; i<m; ++i)
- {
- const double y=line_points[1+i][0];
- for (unsigned int j=0; j<m; ++j)
- {
- const double x=line_points[1+j][0];
-
- w[0] = (1-x)*(1-y);
- w[1] = x*(1-y);
- w[2] = (1-x)*y ;
- 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;
- }
-}
-
-
-template <int spacedim>
-const std::vector<Point<1> > &
-Manifold<spacedim>::
-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<QGaussLobatto<1> >
- quadrature (new QGaussLobatto<1>(n_intermediate_points+2));
- points[n_intermediate_points] = quadrature;
- }
- }
- return points[n_intermediate_points]->get_points();
-}
-
-
/* -------------------------- FlatManifold --------------------- */
projects / dealii-svn.git / commitdiff