#include <base/config.h>
#include <base/point.h>
#include <base/table.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <vector>
#include <string>
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__geometry_info_h
+#define __deal2__geometry_info_h
+
+
+#include <base/config.h>
+#include <base/exceptions.h>
+#include <base/point.h>
+
+
+template <int dim> class GeometryInfo;
+
+
+/**
+ * Topological description of zero dimensional cells,
+ * i.e. points. This class might not look too useful but often is if
+ * in a certain dimension we would like to enquire information about
+ * objects with dimension one lower than the present, e.g. about
+ * faces.
+ *
+ * This class contains as static members information on vertices and
+ * faces of a @p dim-dimensional grid cell. The interface is the same
+ * for all dimensions. If a value is of no use in a low dimensional
+ * cell, it is (correctly) set to zero, e.g. @p subfaces_per_cell in
+ * 1d.
+ *
+ * This information should always replace hard-coded numbers of
+ * vertices, neighbors and so on, since it can be used dimension
+ * independently.
+ *
+ * @ingroup grid
+ * @author Wolfgang Bangerth, 1998
+ */
+template <>
+struct GeometryInfo<0>
+{
+
+ /**
+ * Number of children a cell has.
+ */
+ static const unsigned int children_per_cell = 1;
+
+ /**
+ * Number of faces a cell has.
+ */
+ static const unsigned int faces_per_cell = 0;
+
+ /**
+ * Number of children each face has
+ * when the adjacent cell is refined.
+ */
+ static const unsigned int subfaces_per_face = 0;
+
+ /**
+ * Number of vertices a cell has.
+ */
+ static const unsigned int vertices_per_cell = 1;
+
+ /**
+ * Number of vertices each face has.
+ * Since this is not useful in one
+ * dimension, we provide a useless
+ * number (in the hope that a compiler
+ * may warn when it sees constructs like
+ * <tt>for (i=0; i<vertices_per_face; ++i)</tt>,
+ * at least if @p i is an <tt>unsigned int</tt>.
+ */
+ static const unsigned int vertices_per_face = 0;
+
+ /**
+ * Number of lines each face has.
+ */
+ static const unsigned int lines_per_face = 0;
+
+ /**
+ * Number of quads on each face.
+ */
+ static const unsigned int quads_per_face = 0;
+
+ /**
+ * Number of lines of a cell.
+ */
+ static const unsigned int lines_per_cell = 0;
+
+ /**
+ * Number of quadrilaterals of a
+ * cell.
+ */
+ static const unsigned int quads_per_cell = 0;
+
+ /**
+ * Number of hexahedra of a
+ * cell.
+ */
+ static const unsigned int hexes_per_cell = 0;
+};
+
+
+
+/**
+ * This template specifies the interface to all topological structure
+ * of the mesh cells. This template class can be instantiated for
+ * dim=1,2 and 3.
+ *
+ * @ingroup grid
+ * @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2005
+ */
+template <int dim>
+struct GeometryInfo
+{
+
+ /**
+ * Number of children of a refined cell.
+ */
+ static const unsigned int children_per_cell = 1 << dim;
+
+ /**
+ * Number of faces of a cell.
+ */
+ static const unsigned int faces_per_cell = 2 * dim;
+
+ /**
+ * Number of children each face has
+ * when the adjacent cell is refined.
+ */
+ static const unsigned int subfaces_per_face = GeometryInfo<dim-1>::children_per_cell;
+
+ /**
+ * Number of vertices of a cell.
+ */
+ static const unsigned int vertices_per_cell = 1 << dim;
+
+ /**
+ * Number of vertices on each
+ * face.
+ */
+ static const unsigned int vertices_per_face = GeometryInfo<dim-1>::vertices_per_cell;
+
+ /**
+ * Number of lines on each face.
+ */
+ static const unsigned int lines_per_face
+ = GeometryInfo<dim-1>::lines_per_cell;
+
+ /**
+ * Number of quads on each face.
+ */
+ static const unsigned int quads_per_face
+ = GeometryInfo<dim-1>::quads_per_cell;
+
+ /**
+ * Number of lines of a cell.
+ *
+ * The formula to compute this makes use
+ * of the fact that when going from one
+ * dimension to the next, the object of
+ * the lower dimension is copied once
+ * (thus twice the old number of lines)
+ * and then a new line is inserted
+ * between each vertex of the old object
+ * and the corresponding one in the copy.
+ */
+ static const unsigned int lines_per_cell
+ = (2*GeometryInfo<dim-1>::lines_per_cell +
+ GeometryInfo<dim-1>::vertices_per_cell);
+
+ /**
+ * Number of quadrilaterals of a
+ * cell.
+ *
+ * This number is computed recursively
+ * just as the previous one, with the
+ * exception that new quads result from
+ * connecting an original line and its
+ * copy.
+ */
+ static const unsigned int quads_per_cell
+ = (2*GeometryInfo<dim-1>::quads_per_cell +
+ GeometryInfo<dim-1>::lines_per_cell);
+
+ /**
+ * Number of hexahedra of a
+ * cell.
+ */
+ static const unsigned int hexes_per_cell
+ = (2*GeometryInfo<dim-1>::hexes_per_cell +
+ GeometryInfo<dim-1>::quads_per_cell);
+
+ /**
+ * List of numbers which denotes which
+ * face is opposite to a given face. In
+ * 1d, this list is <tt>{1,0}</tt>, in 2d
+ * <tt>{2, 3, 0, 1}</tt>, in 3d <tt>{1,
+ * 0, 4, 5, 2, 3}</tt>.
+ */
+ static const unsigned int opposite_face[faces_per_cell];
+
+
+ /**
+ * Rearrange vertices for OpenDX
+ * output. For a cell being
+ * written in OpenDX format, each
+ * entry in this field contains
+ * the number of a vertex in
+ * <tt>deal.II</tt> that corresponds
+ * to the DX numbering at this
+ * location.
+ *
+ * Typical example: write a cell
+ * and arrange the vertices, such
+ * that OpenDX understands them.
+ *
+ * \begin{verbatim}
+ * for (i=0; i< n_vertices; ++i)
+ * out << cell->vertex(dx_to_deal[i]);
+ * \end{verbatim}
+ */
+ static const unsigned int dx_to_deal[vertices_per_cell];
+
+ /**
+ * For each face of the reference
+ * cell, this field stores the
+ * coordinate direction in which
+ * its normal vector points.
+ *
+ * Remark that this is only the
+ * coordinate number. The acual
+ * direction of the normal vector
+ * is obtained by multiplying the
+ * unit vector in this direction
+ * with #unit_normal_orientation.
+ */
+ static const unsigned int unit_normal_direction[faces_per_cell];
+
+ /**
+ * Orientation of the unit normal
+ * vector of a face of the
+ * reference cell.
+ *
+ * Each value is either
+ * <tt>1</tt> or <tt>-1</tt>,
+ * corresponding to a normal
+ * vector pointing in the
+ * positive or negative
+ * coordinate direction,
+ * respectively.
+ */
+ static const int unit_normal_orientation[faces_per_cell];
+
+ /**
+ * This field stores which child
+ * cells are adjacent to a
+ * certain face of the mother
+ * cell.
+ *
+ * For example, in 2D the layout of
+ * a cell is as follows:
+ * @verbatim
+ * . 2
+ * . 3-->--2
+ * . | |
+ * . 3 ^ ^ 1
+ * . | |
+ * . 0-->--1
+ * . 0
+ * @endverbatim
+ * Vertices and faces are indicated
+ * with their numbers, faces also with
+ * their directions.
+ *
+ * Now, when refined, the layout is
+ * like this:
+ * @verbatim
+ * *--*--*
+ * | 3|2 |
+ * *--*--*
+ * | 0|1 |
+ * *--*--*
+ * @endverbatim
+ *
+ * Thus, the child cells on face zero
+ * are (ordered in the direction of the
+ * face) 0 and 1, on face 2 they are
+ * 3 and 2, etc.
+ *
+ * For three spatial dimensions,
+ * the exact order of the
+ * children is laid down in the
+ * documentation of the
+ * Triangulation class. Through
+ * the <tt>face_orientation</tt>
+ * argument this function handles
+ * faces oriented in both, the
+ * standard and non-standard
+ * orientation.
+ * <tt>face_orientation</tt>
+ * defaults to <tt>true</tt>
+ * (standard orientation) and has
+ * no effect in 2d.
+ */
+ static unsigned int child_cell_on_face (const unsigned int face,
+ const unsigned int subface,
+ const bool face_orientation = true);
+
+ /**
+ * Map line vertex number to cell
+ * vertex number, i.e. give the
+ * cell vertex number of the
+ * <tt>vertex</tt>th vertex of
+ * line <tt>line</tt>, e.g.
+ * <tt>GeometryInfo<2>::line_to_cell_vertices(2,0)=3</tt>.
+ *
+ * The order of the lines, as
+ * well as their direction (which
+ * in turn determines which is
+ * the first and which the second
+ * vertex on a line) is the
+ * canonical one in deal.II, as
+ * described in the documentation
+ * of the Triangulation
+ * class.
+ *
+ * For <tt>dim=2</tt> this call
+ * is simply passed down to the
+ * face_to_cell_vertices()
+ * function.
+ */
+ static unsigned int line_to_cell_vertices (const unsigned int line,
+ const unsigned int vertex);
+
+ /**
+ * Map face vertex number to cell
+ * vertex number, i.e. give the
+ * cell vertex number of the
+ * <tt>vertex</tt>th vertex of
+ * face <tt>face</tt>, e.g.
+ * <tt>GeometryInfo<2>::face_to_cell_vertices(2,0)=3</tt>.
+ *
+ * Through the
+ * <tt>face_orientation</tt>
+ * argument this function handles
+ * faces oriented in both, the
+ * standard and non-standard
+ * orientation.
+ * <tt>face_orientation</tt>
+ * defaults to <tt>true</tt>
+ * (standard orientation) and has
+ * no effect in 2d.
+ *
+ * As the children of a cell are
+ * ordered according to the
+ * vertices of the cell, this
+ * call is passed down to the
+ * child_cell_on_face() function.
+ * Hence this function is simply
+ * a wrapper of
+ * child_cell_on_face() giving it
+ * a suggestive name.
+ */
+ static unsigned int face_to_cell_vertices (const unsigned int face,
+ const unsigned int vertex,
+ const bool face_orientation = true);
+
+ /**
+ * Map face line number to cell
+ * line number, i.e. give the
+ * cell line number of the
+ * <tt>line</tt>th line of face
+ * <tt>face</tt>, e.g.
+ * <tt>GeometryInfo<3>::face_to_cell_lines(3,1)=5</tt>.
+ *
+ * Through the
+ * <tt>face_orientation</tt>
+ * argument this function handles
+ * faces oriented in both, the
+ * standard and non-standard
+ * orientation.
+ * <tt>face_orientation</tt>
+ * defaults to <tt>true</tt>
+ * (standard orientation).
+ *
+ * This function is useful and
+ * implemented for
+ * <tt>dim=3</tt>, only.
+ */
+ static unsigned int face_to_cell_lines (const unsigned int face,
+ const unsigned int line,
+ const bool face_orientation = true);
+
+ /**
+ * Return the position of the
+ * @p ith vertex on the unit
+ * cell. The order of vertices is
+ * the canonical one in deal.II,
+ * as described in the
+ * documentation of the
+ * Triangulation class.
+ */
+ static Point<dim> unit_cell_vertex (const unsigned int vertex);
+
+ /**
+ * Given a point @p p in unit
+ * coordinates, return the number
+ * of the child cell in which it
+ * would lie in. If the point
+ * lies on the interface of two
+ * children, return any one of
+ * their indices. The result is
+ * always less than
+ * <tt>GeometryInfo<dimension>::children_per_cell</tt>.
+ *
+ * The order of child cells is
+ * described the documentation of
+ * the Triangulation class.
+ */
+ static unsigned int child_cell_from_point (const Point<dim> &p);
+
+ /**
+ * Given coordinates @p p on the
+ * unit cell, return the values
+ * of the coordinates of this
+ * point in the coordinate system
+ * of the given child. Neither
+ * original nor returned
+ * coordinates need actually be
+ * inside the cell, we simply
+ * perform a scale-and-shift
+ * operation with a shift that
+ * depends on the number of the
+ * child.
+ */
+ static Point<dim> cell_to_child_coordinates (const Point<dim> &p,
+ const unsigned int child_index);
+
+ /**
+ * The reverse function to the
+ * one above: take a point in the
+ * coordinate system of the
+ * child, and transform it to the
+ * coordinate system of the
+ * mother cell.
+ */
+ static Point<dim> child_to_cell_coordinates (const Point<dim> &p,
+ const unsigned int child_index);
+
+ /**
+ * Return true if the given point
+ * is inside the unit cell of the
+ * present space dimension.
+ */
+ static bool is_inside_unit_cell (const Point<dim> &p);
+
+ /**
+ * Exception
+ */
+ DeclException1 (ExcInvalidCoordinate,
+ double,
+ << "The coordinates must satisfy 0 <= x_i <= 1, "
+ << "but here we have x_i=" << arg1);
+};
+
+
+/* -------------- declaration of explicit specializations ------------- */
+
+/// @if NoDoc
+
+
+
+template <>
+inline
+Point<1>
+GeometryInfo<1>::unit_cell_vertex (const unsigned int vertex)
+{
+ Assert (vertex < vertices_per_cell,
+ ExcIndexRange (vertex, 0, vertices_per_cell));
+
+ static const Point<1> vertices[vertices_per_cell] =
+ { Point<1>(0.), Point<1>(1.) };
+ return vertices[vertex];
+}
+
+
+
+template <>
+inline
+Point<2>
+GeometryInfo<2>::unit_cell_vertex (const unsigned int vertex)
+{
+ Assert (vertex < vertices_per_cell,
+ ExcIndexRange (vertex, 0, vertices_per_cell));
+
+ static const Point<2> vertices[vertices_per_cell] =
+ { Point<2>(0., 0.), Point<2>(1., 0.),
+ Point<2>(1.,1.), Point<2>(0.,1.) };
+ return vertices[vertex];
+}
+
+
+
+template <>
+inline
+Point<3>
+GeometryInfo<3>::unit_cell_vertex (const unsigned int vertex)
+{
+ Assert (vertex < vertices_per_cell,
+ ExcIndexRange (vertex, 0, vertices_per_cell));
+
+ static const Point<3> vertices[vertices_per_cell] =
+ { Point<3>(0., 0., 0.), Point<3>(1., 0., 0.),
+ Point<3>(1., 0., 1.), Point<3>(0., 0., 1.),
+ Point<3>(0., 1., 0.), Point<3>(1., 1., 0.),
+ Point<3>(1., 1., 1.), Point<3>(0., 1., 1.) };
+ return vertices[vertex];
+}
+
+
+template <int dim>
+inline
+Point<dim>
+GeometryInfo<dim>::unit_cell_vertex (const unsigned int)
+{
+ Assert(false, ExcNotImplemented());
+
+ return Point<dim> ();
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<1>::child_cell_from_point (const Point<1> &p)
+{
+ Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
+
+ return (p[0] <= 0.5 ? 0 : 1);
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<2>::child_cell_from_point (const Point<2> &p)
+{
+ Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
+ Assert ((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
+
+ return (p[0] <= 0.5 ?
+ (p[1] <= 0.5 ? 0 : 3) :
+ (p[1] <= 0.5 ? 1 : 2));
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<3>::child_cell_from_point (const Point<3> &p)
+{
+ Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
+ Assert ((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
+ Assert ((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
+
+ return (p[0] <= 0.5 ?
+ (p[1] <= 0.5 ?
+ (p[2] <= 0.5 ? 0 : 3) :
+ (p[2] <= 0.5 ? 4 : 7)) :
+ (p[1] <= 0.5 ?
+ (p[2] <= 0.5 ? 1 : 2) :
+ (p[2] <= 0.5 ? 5 : 6)));
+}
+
+
+template <int dim>
+inline
+unsigned int
+GeometryInfo<dim>::child_cell_from_point (const Point<dim> &)
+{
+ Assert(false, ExcNotImplemented());
+
+ return 0;
+}
+
+
+
+template <int dim>
+inline
+Point<dim>
+GeometryInfo<dim>::cell_to_child_coordinates (const Point<dim> &p,
+ const unsigned int child_index)
+{
+ Assert (child_index < GeometryInfo<dim>::children_per_cell,
+ ExcIndexRange (child_index, 0, GeometryInfo<dim>::children_per_cell));
+
+ return 2*p - unit_cell_vertex(child_index);
+}
+
+
+
+template <int dim>
+inline
+Point<dim>
+GeometryInfo<dim>::child_to_cell_coordinates (const Point<dim> &p,
+ const unsigned int child_index)
+{
+ Assert (child_index < GeometryInfo<dim>::children_per_cell,
+ ExcIndexRange (child_index, 0, GeometryInfo<dim>::children_per_cell));
+
+ return (p + unit_cell_vertex(child_index))/2;
+}
+
+
+template <>
+inline
+bool
+GeometryInfo<1>::is_inside_unit_cell (const Point<1> &p)
+{
+ return (p[0] >= 0.) && (p[0] <= 1.);
+}
+
+
+
+template <>
+inline
+bool
+GeometryInfo<2>::is_inside_unit_cell (const Point<2> &p)
+{
+ return (p[0] >= 0.) && (p[0] <= 1.) &&
+ (p[1] >= 0.) && (p[1] <= 1.);
+}
+
+
+
+template <>
+inline
+bool
+GeometryInfo<3>::is_inside_unit_cell (const Point<3> &p)
+{
+ return (p[0] >= 0.) && (p[0] <= 1.) &&
+ (p[1] >= 0.) && (p[1] <= 1.) &&
+ (p[2] >= 0.) && (p[2] <= 1.);
+}
+
+
+/// @endif
+
+#endif
//---------------------------------------------------------------------------
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
template <int dim> const unsigned int GeometryInfo<dim>::quads_per_cell;
template <int dim> const unsigned int GeometryInfo<dim>::hexes_per_cell;
+
+using namespace deal_II_numbers;
+
// make sure that also the icc compiler defines (and not only declares)
// these variables
namespace internal
template void define_variables<2> ();
template void define_variables<3> ();
+ template void define_variables<4> ();
}
const unsigned int GeometryInfo<1>::opposite_face[GeometryInfo<1>::faces_per_cell]
= { 0, 1 };
-
template <>
const unsigned int GeometryInfo<2>::opposite_face[GeometryInfo<2>::faces_per_cell]
= { 2, 3, 0, 1 };
-
template <>
const unsigned int GeometryInfo<3>::opposite_face[GeometryInfo<3>::faces_per_cell]
= { 1, 0, 4, 5, 2, 3 };
+template <>
+const unsigned int GeometryInfo<4>::opposite_face[GeometryInfo<4>::faces_per_cell]
+= { invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int };
+
template <>
const unsigned int GeometryInfo<1>::dx_to_deal[GeometryInfo<1>::vertices_per_cell]
= { 0, 1};
-
template <>
const unsigned int GeometryInfo<2>::dx_to_deal[GeometryInfo<2>::vertices_per_cell]
= { 0, 3, 1, 2};
-
template <>
const unsigned int GeometryInfo<3>::dx_to_deal[GeometryInfo<3>::vertices_per_cell]
= { 0, 3, 4, 7, 1, 2, 5, 6};
+template <>
+const unsigned int GeometryInfo<4>::dx_to_deal[GeometryInfo<4>::vertices_per_cell]
+= { invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int,
+ invalid_unsigned_int };
+
+
//TODO: Use these values in mappings.
template <>
const unsigned int
GeometryInfo<3>::unit_normal_orientation[GeometryInfo<3>::faces_per_cell]
= { -1, 1, -1, 1, 1, -1 };
+// These are already for the new, regular numbering
+template <>
+const unsigned int
+GeometryInfo<4>::unit_normal_direction[GeometryInfo<4>::faces_per_cell]
+= { 0, 0, 1, 1, 2, 2, 3, 3 };
+template <>
+const int
+GeometryInfo<4>::unit_normal_orientation[GeometryInfo<4>::faces_per_cell]
+= { -1, 1, -1, 1, -1, 1, -1, 1 };
+
template <>
+template <>
+unsigned int
+GeometryInfo<4>::child_cell_on_face (const unsigned int,
+ const unsigned int,
+ const bool)
+{
+ Assert(false, ExcNotImplemented());
+ return invalid_unsigned_int;
+}
+
+
+
template <>
unsigned int
GeometryInfo<1>::line_to_cell_vertices (const unsigned int line,
}
+template <>
+unsigned int
+GeometryInfo<4>::line_to_cell_vertices (const unsigned int,
+ const unsigned int)
+{
+ Assert(false, ExcNotImplemented());
+ return invalid_unsigned_int;
+}
+
-//TODO:[GK] Is this the correct notion in 1D?
template <>
unsigned int
GeometryInfo<1>::face_to_cell_lines (const unsigned int face,
{
Assert (face+1<faces_per_cell+1, ExcIndexRange(face, 0, faces_per_cell));
Assert (line+1<lines_per_face+1, ExcIndexRange(line, 0, lines_per_face));
-
- return face;
+
+ // There is only a single line, so
+ // it must be this.
+ return 0;
}
Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
Assert (line<lines_per_face, ExcIndexRange(line, 0, lines_per_face));
+ // The face is a line itself.
return face;
}
+template <>
+unsigned int
+GeometryInfo<4>::face_to_cell_lines (const unsigned int,
+ const unsigned int,
+ const bool)
+{
+ Assert(false, ExcNotImplemented());
+ return invalid_unsigned_int;
+}
+
+
+
template <int dim>
inline
unsigned int
template class GeometryInfo<1>;
template class GeometryInfo<2>;
template class GeometryInfo<3>;
-
+template class GeometryInfo<4>;
//---------------------------------------------------------------------------
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <base/quadrature.h>
#include <base/qprojector.h>
#include <base/memory_consumption.h>
#include <base/point.h>
#include <base/tensor.h>
#include <base/table.h>
-#include <grid/geometry_info.h>
+#include <base/vector_slice.h>
+#include <base/geometry_info.h>
#include <lac/full_matrix.h>
#include <fe/fe_update_flags.h>
#include <fe/mapping.h>
#include <string>
+#include <vector>
template<int dim> class FESystem;
*/
bool restriction_is_additive (const unsigned int index) const;
+ /**
+ * @name Support points and interpolation
+ * @{
+ */
+
/**
* Return the support points of
* the trial functions on the
/**
* The function corresponding to
- * the @p unit_support_point
+ * the unit_support_point()
* function, but for faces. See
* there for more information.
*/
Point<dim-1>
unit_face_support_point (const unsigned int index) const;
+ /**
+ * Return a support point vector
+ * for generalized interpolation.
+ */
+ const std::vector<Point<dim> > &
+ get_generalized_support_points () const;
+
+ /**
+ *
+ */
+ bool has_generalized_support_points () const;
+
+ /**
+ *
+ */
+ const std::vector<Point<dim-1> > &
+ get_generalized_face_support_points () const;
+
+ /**
+ * Return whether a finite
+ * element has defined support
+ * points on faces. If the result
+ * is true, then a call to the
+ * @p get_unit_support_points
+ * yields a non-empty array.
+ *
+ * For more information, see the
+ * documentation for the
+ * has_support_points()
+ * function.
+ */
+ bool has_generalized_face_support_points () const;
+
+ /**
+ * Interpolate a set of scalar
+ * values, computed in the
+ * generalized support points.
+ *
+ * @note This function is
+ * implemented in
+ * FiniteElementBase for the case
+ * that the element has support
+ * points. In this case, the
+ * resulting coefficients are
+ * just the values in the suport
+ * points. All other elements
+ * must reimplement it.
+ */
+ virtual void interpolate(std::vector<double>& local_dofs,
+ const std::vector<double>& values) const;
+
+ /**
+ * Interpolate a set of vector
+ * values, computed in the
+ * generalized support points.
+ *
+ * Since a finite element often
+ * only interpolates part of a
+ * vector, <tt>offset</tt> is
+ * used to determine the first
+ * component of the vector to be
+ * interpolated. Maybe consider
+ * changing your data structures
+ * to use the next function.
+ */
+ virtual void interpolate(std::vector<double>& local_dofs,
+ const std::vector<Vector<double> >& values,
+ unsigned int offset = 0) const;
+
+ /**
+ * Interpolate a set of vector
+ * values, computed in the
+ * generalized support points.
+ */
+ virtual void interpolate(
+ std::vector<double>& local_dofs,
+ const VectorSlice<const std::vector<std::vector<double> > >& values) const;
+
+ //@}
+
/**
* Return in which of the vector
* components of this finite
/**
* Exception
+ *
+ * @ingroup Exceptions
*/
DeclException1 (ExcShapeFunctionNotPrimitive,
int,
<< "_component suffix?");
/**
* Exception
+ *
+ * @ingroup Exceptions
*/
DeclException0 (ExcFENotPrimitive);
/**
* Exception
+ *
+ * @ingroup Exceptions
*/
DeclException0 (ExcUnitShapeValuesDoNotExist);
/**
- * Exception
+ * Attempt to access support
+ * points of a finite element
+ * which is not Lagrangian.
+ *
+ * @ingroup Exceptions
*/
DeclException0 (ExcFEHasNoSupportPoints);
/**
- * Exception
+ * Attempt to access embedding
+ * matrices of a finite element
+ * which did not implement these
+ * matrices.
+ *
+ * @ingroup Exceptions
*/
DeclException0 (ExcEmbeddingVoid);
/**
+ * Attempt to access restriction
+ * matrices of a finite element
+ * which did not implement these
+ * matrices.
+ *
* Exception
+ * @ingroup Exceptions
*/
DeclException0 (ExcProjectionVoid);
/**
+ * Attempt to access constraint
+ * matrices of a finite element
+ * which did not implement these
+ * matrices.
+ *
* Exception
+ * @ingroup Exceptions
*/
DeclException0 (ExcConstraintsVoid);
/**
* Exception
+ * @ingroup Exceptions
*/
DeclException2 (ExcWrongInterfaceMatrixSize,
int, int,
<< ", which is not reasonable in the present dimension.");
/**
* Exception
+ * @ingroup Exceptions
*/
DeclException2 (ExcComponentIndexInvalid,
int, int,
<< ") is invalid, i.e. non-existent");
/**
* Exception
+ * @ingroup Exceptions
*/
DeclException0 (ExcInterpolationNotImplemented);
protected:
/**
* Array of projection matrices. See
- * @p get_restriction_matrix above.
+ * get_restriction_matrix() above.
*
* Matrices in this array are
* automatically initialized to
* support point.
*/
std::vector<Point<dim-1> > unit_face_support_points;
-
+
+ /**
+ * Support points used for
+ * interpolation functions of
+ * non-Lagrangian elements.
+ */
+ std::vector<Point<dim> > generalized_support_points;
+
+ /**
+ * Face support points used for
+ * interpolation functions of
+ * non-Lagrangian elements.
+ */
+ std::vector<Point<dim-1> > generalized_face_support_points;
+
/**
* For each shape function, give
* a vector of bools (with size
#define __deal2__fe_nedelec_h
#include <base/config.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <fe/fe.h>
template <int dim> class MappingQ;
#include <base/config.h>
#include <base/polynomials_raviart_thomas.h>
#include <base/tensor_product_polynomials.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <fe/fe.h>
#include <fe/fe_poly_tensor.h>
// further information on this license.
//
//---------------------------------------------------------------------------
-#ifndef __deal2__geometry_info_h
-#define __deal2__geometry_info_h
+// This class was moved to the base library. Therefore, we provide
+// this file for compatibility reasons only.
-#include <base/config.h>
-#include <base/exceptions.h>
-#include <base/point.h>
-
-
-template <int dim> class GeometryInfo;
-
-
-/**
- * Topological description of zero dimensional cells,
- * i.e. points. This class might not look too useful but often is if
- * in a certain dimension we would like to enquire information about
- * objects with dimension one lower than the present, e.g. about
- * faces.
- *
- * This class contains as static members information on vertices and
- * faces of a @p dim-dimensional grid cell. The interface is the same
- * for all dimensions. If a value is of no use in a low dimensional
- * cell, it is (correctly) set to zero, e.g. @p subfaces_per_cell in
- * 1d.
- *
- * This information should always replace hard-coded numbers of
- * vertices, neighbors and so on, since it can be used dimension
- * independently.
- *
- * @ingroup grid
- * @author Wolfgang Bangerth, 1998
- */
-template <>
-struct GeometryInfo<0>
-{
-
- /**
- * Number of children a cell has.
- */
- static const unsigned int children_per_cell = 1;
-
- /**
- * Number of faces a cell has.
- */
- static const unsigned int faces_per_cell = 0;
-
- /**
- * Number of children each face has
- * when the adjacent cell is refined.
- */
- static const unsigned int subfaces_per_face = 0;
-
- /**
- * Number of vertices a cell has.
- */
- static const unsigned int vertices_per_cell = 1;
-
- /**
- * Number of vertices each face has.
- * Since this is not useful in one
- * dimension, we provide a useless
- * number (in the hope that a compiler
- * may warn when it sees constructs like
- * <tt>for (i=0; i<vertices_per_face; ++i)</tt>,
- * at least if @p i is an <tt>unsigned int</tt>.
- */
- static const unsigned int vertices_per_face = 0;
-
- /**
- * Number of lines each face has.
- */
- static const unsigned int lines_per_face = 0;
-
- /**
- * Number of quads on each face.
- */
- static const unsigned int quads_per_face = 0;
-
- /**
- * Number of lines of a cell.
- */
- static const unsigned int lines_per_cell = 0;
-
- /**
- * Number of quadrilaterals of a
- * cell.
- */
- static const unsigned int quads_per_cell = 0;
-
- /**
- * Number of hexahedra of a
- * cell.
- */
- static const unsigned int hexes_per_cell = 0;
-};
-
-
-
-/**
- * Topological description of four dimensional cells. This class is
- * required in some exotic cases where we compute information in a
- * one-larger dimension than the present, and do so also in 3d (for
- * example, stacking the solutions of a d-dimensional time dependent
- * equation timestep after timestep in a (d+1)-dimensional space).
- *
- * This class contains as static members information on vertices and
- * faces of a @p dim-dimensional grid cell. The interface is the same
- * for all dimensions. If a value is of no use in a low dimensional
- * cell, it is (correctly) set to zero, e.g. @p subfaces_per_cell in
- * 1d.
- *
- * This information should always replace hard-coded numbers of
- * vertices, neighbors and so on, since it can be used dimension
- * independently.
- *
- * @ingroup grid
- * @author Wolfgang Bangerth, 1998
- */
-template <>
-struct GeometryInfo<4>
-{
-
- /**
- * Number of children a cell has.
- */
- static const unsigned int children_per_cell = 16;
-
- /**
- * Number of faces a cell has.
- */
- static const unsigned int faces_per_cell = 8;
-
- /**
- * Number of children each face has
- * when the adjacent cell is refined.
- */
- static const unsigned int subfaces_per_face = 8;
-
- /**
- * Number of vertices a cell has.
- */
- static const unsigned int vertices_per_cell = 16;
-
- /**
- * Number of vertices each face
- * has.
- */
- static const unsigned int vertices_per_face = 8;
-
- /**
- * Number of lines each face has.
- */
- static const unsigned int lines_per_face = 12;
-
- /**
- * Number of quads on each face.
- */
- static const unsigned int quads_per_face = 6;
-
- /**
- * Number of lines of a cell.
- */
- static const unsigned int lines_per_cell = 32;
-
- /**
- * Number of quadrilaterals of a
- * cell.
- */
- static const unsigned int quads_per_cell = 24;
-
- /**
- * Number of hexahedra of a
- * cell.
- */
- static const unsigned int hexes_per_cell = 8;
-};
-
-
-
-/**
- * This template specifies the interface to all topological structure
- * of the mesh cells. This template class can be instantiated for
- * dim=1,2 and 3.
- *
- * @ingroup grid
- * @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2005
- */
-template <int dim>
-struct GeometryInfo
-{
-
- /**
- * Number of children of a refined cell.
- */
- static const unsigned int children_per_cell = 1 << dim;
-
- /**
- * Number of faces of a cell.
- */
- static const unsigned int faces_per_cell = 2 * dim;
-
- /**
- * Number of children each face has
- * when the adjacent cell is refined.
- */
- static const unsigned int subfaces_per_face = GeometryInfo<dim-1>::children_per_cell;
-
- /**
- * Number of vertices of a cell.
- */
- static const unsigned int vertices_per_cell = 1 << dim;
-
- /**
- * Number of vertices on each
- * face.
- */
- static const unsigned int vertices_per_face = GeometryInfo<dim-1>::vertices_per_cell;
-
- /**
- * Number of lines on each face.
- */
- static const unsigned int lines_per_face
- = GeometryInfo<dim-1>::lines_per_cell;
-
- /**
- * Number of quads on each face.
- */
- static const unsigned int quads_per_face
- = GeometryInfo<dim-1>::quads_per_cell;
-
- /**
- * Number of lines of a cell.
- *
- * The formula to compute this makes use
- * of the fact that when going from one
- * dimension to the next, the object of
- * the lower dimension is copied once
- * (thus twice the old number of lines)
- * and then a new line is inserted
- * between each vertex of the old object
- * and the corresponding one in the copy.
- */
- static const unsigned int lines_per_cell
- = (2*GeometryInfo<dim-1>::lines_per_cell +
- GeometryInfo<dim-1>::vertices_per_cell);
-
- /**
- * Number of quadrilaterals of a
- * cell.
- *
- * This number is computed recursively
- * just as the previous one, with the
- * exception that new quads result from
- * connecting an original line and its
- * copy.
- */
- static const unsigned int quads_per_cell
- = (2*GeometryInfo<dim-1>::quads_per_cell +
- GeometryInfo<dim-1>::lines_per_cell);
-
- /**
- * Number of hexahedra of a
- * cell.
- */
- static const unsigned int hexes_per_cell
- = (2*GeometryInfo<dim-1>::hexes_per_cell +
- GeometryInfo<dim-1>::quads_per_cell);
-
- /**
- * List of numbers which denotes which
- * face is opposite to a given face. In
- * 1d, this list is <tt>{1,0}</tt>, in 2d
- * <tt>{2, 3, 0, 1}</tt>, in 3d <tt>{1,
- * 0, 4, 5, 2, 3}</tt>.
- */
- static const unsigned int opposite_face[faces_per_cell];
-
-
- /**
- * Rearrange vertices for OpenDX
- * output. For a cell being
- * written in OpenDX format, each
- * entry in this field contains
- * the number of a vertex in
- * <tt>deal.II</tt> that corresponds
- * to the DX numbering at this
- * location.
- *
- * Typical example: write a cell
- * and arrange the vertices, such
- * that OpenDX understands them.
- *
- * \begin{verbatim}
- * for (i=0; i< n_vertices; ++i)
- * out << cell->vertex(dx_to_deal[i]);
- * \end{verbatim}
- */
- static const unsigned int dx_to_deal[vertices_per_cell];
-
- /**
- * For each face of the reference
- * cell, this field stores the
- * coordinate direction in which
- * its normal vector points.
- *
- * Remark that this is only the
- * coordinate number. The acual
- * direction of the normal vector
- * is obtained by multiplying the
- * unit vector in this direction
- * with #unit_normal_orientation.
- */
- static const unsigned int unit_normal_direction[faces_per_cell];
-
- /**
- * Orientation of the unit normal
- * vector of a face of the
- * reference cell.
- *
- * Each value is either
- * <tt>1</tt> or <tt>-1</tt>,
- * corresponding to a normal
- * vector pointing in the
- * positive or negative
- * coordinate direction,
- * respectively.
- */
- static const int unit_normal_orientation[faces_per_cell];
-
- /**
- * This field stores which child
- * cells are adjacent to a
- * certain face of the mother
- * cell.
- *
- * For example, in 2D the layout of
- * a cell is as follows:
- * @verbatim
- * . 2
- * . 3-->--2
- * . | |
- * . 3 ^ ^ 1
- * . | |
- * . 0-->--1
- * . 0
- * @endverbatim
- * Vertices and faces are indicated
- * with their numbers, faces also with
- * their directions.
- *
- * Now, when refined, the layout is
- * like this:
- * @verbatim
- * *--*--*
- * | 3|2 |
- * *--*--*
- * | 0|1 |
- * *--*--*
- * @endverbatim
- *
- * Thus, the child cells on face zero
- * are (ordered in the direction of the
- * face) 0 and 1, on face 2 they are
- * 3 and 2, etc.
- *
- * For three spatial dimensions,
- * the exact order of the
- * children is laid down in the
- * documentation of the
- * Triangulation class. Through
- * the <tt>face_orientation</tt>
- * argument this function handles
- * faces oriented in both, the
- * standard and non-standard
- * orientation.
- * <tt>face_orientation</tt>
- * defaults to <tt>true</tt>
- * (standard orientation) and has
- * no effect in 2d.
- */
- static unsigned int child_cell_on_face (const unsigned int face,
- const unsigned int subface,
- const bool face_orientation = true);
-
- /**
- * Map line vertex number to cell
- * vertex number, i.e. give the
- * cell vertex number of the
- * <tt>vertex</tt>th vertex of
- * line <tt>line</tt>, e.g.
- * <tt>GeometryInfo<2>::line_to_cell_vertices(2,0)=3</tt>.
- *
- * The order of the lines, as
- * well as their direction (which
- * in turn determines which is
- * the first and which the second
- * vertex on a line) is the
- * canonical one in deal.II, as
- * described in the documentation
- * of the Triangulation
- * class.
- *
- * For <tt>dim=2</tt> this call
- * is simply passed down to the
- * face_to_cell_vertices()
- * function.
- */
- static unsigned int line_to_cell_vertices (const unsigned int line,
- const unsigned int vertex);
-
- /**
- * Map face vertex number to cell
- * vertex number, i.e. give the
- * cell vertex number of the
- * <tt>vertex</tt>th vertex of
- * face <tt>face</tt>, e.g.
- * <tt>GeometryInfo<2>::face_to_cell_vertices(2,0)=3</tt>.
- *
- * Through the
- * <tt>face_orientation</tt>
- * argument this function handles
- * faces oriented in both, the
- * standard and non-standard
- * orientation.
- * <tt>face_orientation</tt>
- * defaults to <tt>true</tt>
- * (standard orientation) and has
- * no effect in 2d.
- *
- * As the children of a cell are
- * ordered according to the
- * vertices of the cell, this
- * call is passed down to the
- * child_cell_on_face() function.
- * Hence this function is simply
- * a wrapper of
- * child_cell_on_face() giving it
- * a suggestive name.
- */
- static unsigned int face_to_cell_vertices (const unsigned int face,
- const unsigned int vertex,
- const bool face_orientation = true);
-
- /**
- * Map face line number to cell
- * line number, i.e. give the
- * cell line number of the
- * <tt>line</tt>th line of face
- * <tt>face</tt>, e.g.
- * <tt>GeometryInfo<3>::face_to_cell_lines(3,1)=5</tt>.
- *
- * Through the
- * <tt>face_orientation</tt>
- * argument this function handles
- * faces oriented in both, the
- * standard and non-standard
- * orientation.
- * <tt>face_orientation</tt>
- * defaults to <tt>true</tt>
- * (standard orientation).
- *
- * This function is useful and
- * implemented for
- * <tt>dim=3</tt>, only.
- */
- static unsigned int face_to_cell_lines (const unsigned int face,
- const unsigned int line,
- const bool face_orientation = true);
-
- /**
- * Return the position of the
- * @p ith vertex on the unit
- * cell. The order of vertices is
- * the canonical one in deal.II,
- * as described in the
- * documentation of the
- * Triangulation class.
- */
- static Point<dim> unit_cell_vertex (const unsigned int vertex);
-
- /**
- * Given a point @p p in unit
- * coordinates, return the number
- * of the child cell in which it
- * would lie in. If the point
- * lies on the interface of two
- * children, return any one of
- * their indices. The result is
- * always less than
- * <tt>GeometryInfo<dimension>::children_per_cell</tt>.
- *
- * The order of child cells is
- * described the documentation of
- * the Triangulation class.
- */
- static unsigned int child_cell_from_point (const Point<dim> &p);
-
- /**
- * Given coordinates @p p on the
- * unit cell, return the values
- * of the coordinates of this
- * point in the coordinate system
- * of the given child. Neither
- * original nor returned
- * coordinates need actually be
- * inside the cell, we simply
- * perform a scale-and-shift
- * operation with a shift that
- * depends on the number of the
- * child.
- */
- static Point<dim> cell_to_child_coordinates (const Point<dim> &p,
- const unsigned int child_index);
-
- /**
- * The reverse function to the
- * one above: take a point in the
- * coordinate system of the
- * child, and transform it to the
- * coordinate system of the
- * mother cell.
- */
- static Point<dim> child_to_cell_coordinates (const Point<dim> &p,
- const unsigned int child_index);
-
- /**
- * Return true if the given point
- * is inside the unit cell of the
- * present space dimension.
- */
- static bool is_inside_unit_cell (const Point<dim> &p);
-
- /**
- * Exception
- */
- DeclException1 (ExcInvalidCoordinate,
- double,
- << "The coordinates must satisfy 0 <= x_i <= 1, "
- << "but here we have x_i=" << arg1);
-};
-
-
-/* -------------- declaration of explicit specializations ------------- */
-
-/// @if NoDoc
-
-
-
-template <>
-inline
-Point<1>
-GeometryInfo<1>::unit_cell_vertex (const unsigned int vertex)
-{
- Assert (vertex < vertices_per_cell,
- ExcIndexRange (vertex, 0, vertices_per_cell));
-
- static const Point<1> vertices[vertices_per_cell] =
- { Point<1>(0.), Point<1>(1.) };
- return vertices[vertex];
-}
-
-
-
-template <>
-inline
-Point<2>
-GeometryInfo<2>::unit_cell_vertex (const unsigned int vertex)
-{
- Assert (vertex < vertices_per_cell,
- ExcIndexRange (vertex, 0, vertices_per_cell));
-
- static const Point<2> vertices[vertices_per_cell] =
- { Point<2>(0., 0.), Point<2>(1., 0.),
- Point<2>(1.,1.), Point<2>(0.,1.) };
- return vertices[vertex];
-}
-
-
-
-template <>
-inline
-Point<3>
-GeometryInfo<3>::unit_cell_vertex (const unsigned int vertex)
-{
- Assert (vertex < vertices_per_cell,
- ExcIndexRange (vertex, 0, vertices_per_cell));
-
- static const Point<3> vertices[vertices_per_cell] =
- { Point<3>(0., 0., 0.), Point<3>(1., 0., 0.),
- Point<3>(1., 0., 1.), Point<3>(0., 0., 1.),
- Point<3>(0., 1., 0.), Point<3>(1., 1., 0.),
- Point<3>(1., 1., 1.), Point<3>(0., 1., 1.) };
- return vertices[vertex];
-}
-
-
-template <int dim>
-inline
-Point<dim>
-GeometryInfo<dim>::unit_cell_vertex (const unsigned int)
-{
- Assert(false, ExcNotImplemented());
-
- return Point<dim> ();
-}
-
-
-
-template <>
-inline
-unsigned int
-GeometryInfo<1>::child_cell_from_point (const Point<1> &p)
-{
- Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
-
- return (p[0] <= 0.5 ? 0 : 1);
-}
-
-
-
-template <>
-inline
-unsigned int
-GeometryInfo<2>::child_cell_from_point (const Point<2> &p)
-{
- Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
- Assert ((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
-
- return (p[0] <= 0.5 ?
- (p[1] <= 0.5 ? 0 : 3) :
- (p[1] <= 0.5 ? 1 : 2));
-}
-
-
-
-template <>
-inline
-unsigned int
-GeometryInfo<3>::child_cell_from_point (const Point<3> &p)
-{
- Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
- Assert ((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
- Assert ((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
-
- return (p[0] <= 0.5 ?
- (p[1] <= 0.5 ?
- (p[2] <= 0.5 ? 0 : 3) :
- (p[2] <= 0.5 ? 4 : 7)) :
- (p[1] <= 0.5 ?
- (p[2] <= 0.5 ? 1 : 2) :
- (p[2] <= 0.5 ? 5 : 6)));
-}
-
-
-template <int dim>
-inline
-unsigned int
-GeometryInfo<dim>::child_cell_from_point (const Point<dim> &)
-{
- Assert(false, ExcNotImplemented());
-
- return 0;
-}
-
-
-
-template <int dim>
-inline
-Point<dim>
-GeometryInfo<dim>::cell_to_child_coordinates (const Point<dim> &p,
- const unsigned int child_index)
-{
- Assert (child_index < GeometryInfo<dim>::children_per_cell,
- ExcIndexRange (child_index, 0, GeometryInfo<dim>::children_per_cell));
-
- return 2*p - unit_cell_vertex(child_index);
-}
-
-
-
-template <int dim>
-inline
-Point<dim>
-GeometryInfo<dim>::child_to_cell_coordinates (const Point<dim> &p,
- const unsigned int child_index)
-{
- Assert (child_index < GeometryInfo<dim>::children_per_cell,
- ExcIndexRange (child_index, 0, GeometryInfo<dim>::children_per_cell));
-
- return (p + unit_cell_vertex(child_index))/2;
-}
-
-
-template <>
-inline
-bool
-GeometryInfo<1>::is_inside_unit_cell (const Point<1> &p)
-{
- return (p[0] >= 0.) && (p[0] <= 1.);
-}
-
-
-
-template <>
-inline
-bool
-GeometryInfo<2>::is_inside_unit_cell (const Point<2> &p)
-{
- return (p[0] >= 0.) && (p[0] <= 1.) &&
- (p[1] >= 0.) && (p[1] <= 1.);
-}
-
-
-
-template <>
-inline
-bool
-GeometryInfo<3>::is_inside_unit_cell (const Point<3> &p)
-{
- return (p[0] >= 0.) && (p[0] <= 1.) &&
- (p[1] >= 0.) && (p[1] <= 1.) &&
- (p[2] >= 0.) && (p[2] <= 1.);
-}
-
-
-/// @endif
-
-#endif
+#include <base/geometry_info.h>
#include <base/point.h>
#include <base/subscriptor.h>
#include <base/smartpointer.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <grid/tria_iterator_selector.h>
#include <vector>
* object, you should be well aware that you might involuntarily alter the
* data stored in the triangulation.
*
- * @ref TriaRawIterator
- * @author Wolfgang Bangerth, 1998
+ * @ingroup grid
+ * @author Wolfgang Bangerth, 1998
*/
template <int dim>
class Triangulation : public Subscriptor
#include <grid/tria_iterator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.templates.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <cmath>
#include <grid/tria_iterator.h>
#include <grid/tria_levels.h>
#include <grid/tria.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <fe/fe.h>
#include <set>
+template <int dim>
+const std::vector<Point<dim> > &
+FiniteElementBase<dim>::get_generalized_support_points () const
+{
+ // a finite element may define
+ // support points, but only if
+ // there are as many as there are
+ // degrees of freedom
+ Assert ((generalized_support_points.size() == 0) ||
+ (generalized_support_points.size() == this->dofs_per_cell),
+ ExcInternalError());
+ return generalized_support_points;
+}
+
+
+
+template <int dim>
+bool
+FiniteElementBase<dim>::has_generalized_support_points () const
+{
+ return (generalized_support_points.size() != 0);
+}
+
+
+
template <int dim>
Point<dim>
FiniteElementBase<dim>::unit_support_point (const unsigned index) const
+template <int dim>
+const std::vector<Point<dim-1> > &
+FiniteElementBase<dim>::get_generalized_face_support_points () const
+{
+ // a finite element may define
+ // support points, but only if
+ // there are as many as there are
+ // degrees of freedom on a face
+ Assert ((generalized_face_support_points.size() == 0) ||
+ (generalized_face_support_points.size() == this->dofs_per_face),
+ ExcInternalError());
+ return generalized_face_support_points;
+}
+
+
+
+template <int dim>
+bool
+FiniteElementBase<dim>::has_generalized_face_support_points () const
+{
+ return (generalized_face_support_points.size() != 0);
+}
+
+
+
template <int dim>
Point<dim-1>
FiniteElementBase<dim>::unit_face_support_point (const unsigned index) const
+template <int dim>
+void
+FiniteElementBase<dim>::interpolate(
+ std::vector<double>& local_dofs,
+ const std::vector<double>& values) const
+{
+ Assert (has_support_points(), ExcFEHasNoSupportPoints());
+ Assert (values.size() == unit_support_points.size(),
+ ExcDimensionMismatch(values.size(), unit_support_points.size()));
+ Assert (local_dofs.size() == this->dofs_per_cell,
+ ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+ Assert (this->n_components() == 1,
+ ExcDimensionMismatch(this->n_components(), 1));
+
+ Assert(false, ExcNotImplemented());
+
+// std::fill(values.begin(), values.end(), local_dofs.begin());
+}
+
+
+
+
+template <int dim>
+void
+FiniteElementBase<dim>::interpolate(
+ std::vector<double>& local_dofs,
+ const std::vector<Vector<double> >& values,
+ unsigned int offset) const
+{
+ Assert (has_support_points(), ExcFEHasNoSupportPoints());
+ Assert (values.size() == unit_support_points.size(),
+ ExcDimensionMismatch(values.size(), unit_support_points.size()));
+ Assert (local_dofs.size() == this->dofs_per_cell,
+ ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+ Assert (values[0].size() >= offset+this->n_components(),
+ ExcDimensionMismatch(values[0].size(),offset+this->n_components()));
+
+ for (unsigned int i=0;i<this->dofs_per_cell;++i)
+ {
+ const std::pair<unsigned int, unsigned int> index
+ = this->system_to_component_index(i);
+ local_dofs[i] = values[index.second](offset+index.first);
+ }
+}
+
+
+
+
+template <int dim>
+void
+FiniteElementBase<dim>::interpolate(
+ std::vector<double>& local_dofs,
+ const VectorSlice<const std::vector<std::vector<double> > >& values) const
+{
+ Assert (has_support_points(), ExcFEHasNoSupportPoints());
+ Assert (values[0].size() == unit_support_points.size(),
+ ExcDimensionMismatch(values.size(), unit_support_points.size()));
+ Assert (local_dofs.size() == this->dofs_per_cell,
+ ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+ Assert (values.size() == this->n_components(),
+ ExcDimensionMismatch(values.size(), this->n_components()));
+
+ for (unsigned int i=0;i<this->dofs_per_cell;++i)
+ {
+ const std::pair<unsigned int, unsigned int> index
+ = this->system_to_component_index(i);
+ local_dofs[i] = values[index.first][index.second];
+ }
+}
+
+
+
+
template <int dim>
unsigned int
FiniteElementBase<dim>::memory_consumption () const
//---------------------------------------------------------------------------
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <fe/fe.h>
#include <grid/tria_boundary.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <grid/grid_tools.h>
#include <grid/magic_numbers.h>
#include <lac/vector.h>
#include <grid/tria_accessor.h>
#include <grid/tria_accessor.templates.h>
#include <grid/tria_iterator.templates.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <grid/grid_tools.h>
#include <fe/mapping_q1.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/tria.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <fe/fe.h>
#include <lac/sparse_matrix.h>
#include <base/exceptions.h>
#include <lac/petsc_vector.h>
#include <lac/petsc_block_vector.h>
#include <grid/tria_iterator.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <fe/fe.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <grid/tria_iterator.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <base/quadrature.h>
#include <fe/fe.h>
#include <fe/fe_values.h>
#include "../tests.h"
#include <base/logstream.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <fstream>
#include <cstdlib>
#include "../tests.h"
#include <base/logstream.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <fstream>
#include <cstdlib>
#include <dofs/dof_accessor.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_constraints.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <fe/fe.h>
#include <fe/fe_values.h>
#include <base/quadrature.h>
#include <lac/vector.h>
#include <dofs/dof_constraints.h>
#include <dofs/dof_handler.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <cmath>
#include <algorithm>
#include <lac/solver_cg.h>
#include <lac/vector_memory.h>
#include <lac/precondition.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <dofs/dof_constraints.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_accessor.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
#include <numerics/data_out_stack.h>
#include <fstream>
projects / dealii-svn.git / commitdiff