#include <deal.II/base/tensor_product_polynomials.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary.h>
+#include <deal.II/grid/manifold.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/mapping_q.h>
}
-
-template <>
-void
-MappingQ<1,2>::add_line_support_points (const Triangulation<1,2>::cell_iterator &cell,
- std::vector<Point<2> > &a) const
-{
- const unsigned int dim = 1;
- const unsigned int spacedim = 2;
- // Ask for the mid point, if that's the only thing we need.
- if (degree==2)
- {
- const Boundary<dim,spacedim> *const boundary
- = &(cell->get_triangulation().get_boundary(cell->material_id()));
- a.push_back(boundary->get_new_point_on_line(cell));
- }
- else
- // otherwise call the more complicated functions and ask for inner points
- // from the boundary description
- {
- std::vector<Point<spacedim> > line_points (degree-1);
-
- const Boundary<dim,spacedim> *const boundary
- = &(cell->get_triangulation().get_boundary(cell->material_id()));
-
- boundary->get_intermediate_points_on_line (cell, line_points);
- // Append all points
- a.insert (a.end(), line_points.begin(), line_points.end());
- }
-}
-
-
-
-template <>
-void
-MappingQ<1,3>::add_line_support_points (const Triangulation<1,3>::cell_iterator &cell,
- std::vector<Point<3> > &a) const
-{
- const unsigned int dim = 1;
- const unsigned int spacedim = 3;
- // Ask for the mid point, if that's the only thing we need.
- if (degree==2)
- {
- const Boundary<dim,spacedim> *const boundary
- = &(cell->get_triangulation().get_boundary(cell->material_id()));
- a.push_back(boundary->get_new_point_on_line(cell));
- }
- else
- // otherwise call the more complicated functions and ask for inner points
- // from the boundary description
- {
- std::vector<Point<spacedim> > line_points (degree-1);
-
- const Boundary<dim,spacedim> *const boundary
- = &(cell->get_triangulation().get_boundary(cell->material_id()));
-
- boundary->get_intermediate_points_on_line (cell, line_points);
- // Append all points
- a.insert (a.end(), line_points.begin(), line_points.end());
- }
-}
-
-
-
template<int dim, int spacedim>
void
MappingQ<dim,spacedim>::add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
std::vector<Point<spacedim> > &a) const
{
- static const StraightBoundary<dim,spacedim> straight_boundary;
+ std::vector<Point<spacedim> > ps(2);
+ std::vector<double> ws(2, .5);
+
// if we only need the midpoint, then ask for it.
if (degree==2)
{
for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
{
const typename Triangulation<dim,spacedim>::line_iterator line = cell->line(line_no);
- const Boundary<dim,spacedim> *const boundary
- = (line->at_boundary()?
- &line->get_triangulation().get_boundary(line->boundary_indicator()):
- (dim != spacedim) ?
- &line->get_triangulation().get_boundary(cell->material_id()):
- &straight_boundary);
-
- a.push_back(boundary->get_new_point_on_line(line));
+ ps[0] = line->vertex(0);
+ ps[1] = line->vertex(1);
+ a.push_back(line->get_manifold().get_new_point(ps,ws));
};
}
else
for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
{
const typename Triangulation<dim,spacedim>::line_iterator line = cell->line(line_no);
+ ps[0] = line->vertex(0);
+ ps[1] = line->vertex(1);
- const Boundary<dim,spacedim> *const boundary
- = (line->at_boundary()?
- &line->get_triangulation().get_boundary(line->boundary_indicator()) :
- (dim != spacedim) ?
- &line->get_triangulation().get_boundary(cell->material_id()) :
- &straight_boundary);
-
- boundary->get_intermediate_points_on_line (line, line_points);
+ line->get_manifold().get_intermediate_points (line_points, ps);
+
if (dim==3)
{
// in 3D, lines might be in wrong orientation. if so, reverse
lines_per_face = GeometryInfo<3>::lines_per_face,
vertices_per_cell = GeometryInfo<3>::vertices_per_cell;
- static const StraightBoundary<3> straight_boundary;
- // used if face quad at boundary or entirely in the interior of the domain
+ // Used by the new Manifold
+ // interface. This vector collects
+ // the vertices used to compute the
+ // intermediate points.
+ std::vector<Point<3> > vertices(4);
+
std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
// used if only one line of face quad is at boundary
std::vector<Point<3> > b(4*degree);
// points on it
if (face->at_boundary())
{
- face->get_triangulation().get_boundary(face->boundary_indicator())
- .get_intermediate_points_on_quad (face, quad_points);
+ vertices.resize(4);
+ for(unsigned int i=0;i<4; ++i)
+ vertices[i] = face->vertex(i);
+ face->get_manifold().get_intermediate_points(quad_points, vertices);
// 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
}
else
{
- // face is entirely in the interior. get intermediate points
- // from a straight boundary object
- straight_boundary.get_intermediate_points_on_quad (face, quad_points);
+ // face is entirely in the interior.
+ 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);
// 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> > &a) const
{
std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
-
- cell->get_triangulation().get_boundary(cell->material_id())
- .get_intermediate_points_on_quad (cell, quad_points);
+ 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 (quad_points, vertices);
for (unsigned int i=0; i<quad_points.size(); ++i)
a.push_back(quad_points[i]);
}