]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
GeometryInfo moved; declaration of generalized support points and interpolation in...
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 19 Jul 2005 09:43:00 +0000 (09:43 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 19 Jul 2005 09:43:00 +0000 (09:43 +0000)
git-svn-id: https://svn.dealii.org/trunk@11174 0785d39b-7218-0410-832d-ea1e28bc413d

21 files changed:
deal.II/base/include/base/data_out_base.h
deal.II/base/include/base/geometry_info.h [new file with mode: 0644]
deal.II/base/source/geometry_info.cc [moved from deal.II/deal.II/source/grid/geometry_info.cc with 78% similarity]
deal.II/base/source/quadrature.cc
deal.II/deal.II/include/fe/fe_base.h
deal.II/deal.II/include/fe/fe_nedelec.h
deal.II/deal.II/include/fe/fe_raviart_thomas.h
deal.II/deal.II/include/grid/geometry_info.h
deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/include/grid/tria_accessor.templates.h
deal.II/deal.II/source/dofs/dof_handler.cc
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_data.cc
deal.II/deal.II/source/grid/tria.cc
deal.II/deal.II/source/grid/tria_accessor.cc
deal.II/deal.II/source/multigrid/mg_dof_handler.cc
deal.II/deal.II/source/numerics/error_estimator.cc
deal.II/deal.II/source/numerics/matrices.cc
tests/bits/geometry_info_1.cc
tests/bits/geometry_info_2.cc
tests/deal.II/wave-test-3.cc

index fd8ac652e5d915a742cd7c5cd3ee0e655a6c2ea0..8e5a897daefda3d7030b0159e73b46d6c683b782 100644 (file)
@@ -17,7 +17,7 @@
 #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>
diff --git a/deal.II/base/include/base/geometry_info.h b/deal.II/base/include/base/geometry_info.h
new file mode 100644 (file)
index 0000000..30c6deb
--- /dev/null
@@ -0,0 +1,659 @@
+//---------------------------------------------------------------------------
+//    $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
similarity index 78%
rename from deal.II/deal.II/source/grid/geometry_info.cc
rename to deal.II/base/source/geometry_info.cc
index 7280afa11d42bfd5bd4d71af28689715cf303331..c9b1806caa74d8571f3a7e10074f3012228f88af 100644 (file)
@@ -12,7 +12,7 @@
 //---------------------------------------------------------------------------
 
 
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
 
 
 
@@ -27,6 +27,9 @@ template <int dim> const unsigned int GeometryInfo<dim>::lines_per_cell;
 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
@@ -41,6 +44,7 @@ namespace internal
 
   template void define_variables<2> ();
   template void define_variables<3> ();
+  template void define_variables<4> ();
 }
 
 
@@ -50,32 +54,59 @@ template <>
 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
@@ -104,6 +135,16 @@ const 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 <>
@@ -162,6 +203,18 @@ GeometryInfo<3>::child_cell_on_face (const unsigned int face,
 
 
 
+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,
@@ -209,8 +262,16 @@ GeometryInfo<3>::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,
@@ -219,8 +280,10 @@ 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;
 }
 
 
@@ -234,6 +297,7 @@ GeometryInfo<2>::face_to_cell_lines (const unsigned int face,
   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;
 }
 
@@ -261,6 +325,18 @@ GeometryInfo<3>::face_to_cell_lines (const unsigned int 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
@@ -277,4 +353,4 @@ GeometryInfo<dim>::face_to_cell_vertices (const unsigned int face,
 template class GeometryInfo<1>;
 template class GeometryInfo<2>;
 template class GeometryInfo<3>;
-
+template class GeometryInfo<4>;
index e1d9f01ee7fa27bde4d0da6901db2bda593ea229..08edb028e13f2cd14c1a23339804d2f9509bdf60 100644 (file)
@@ -12,7 +12,7 @@
 //---------------------------------------------------------------------------
 
 
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
 #include <base/quadrature.h>
 #include <base/qprojector.h>
 #include <base/memory_consumption.h>
index e755bc5547556a0bdbbf8c16a097efbb86d305a5..4296c5b1cfa76a1a76ba82c37ae0b08b96d5ff46 100644 (file)
 #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;
 
@@ -1231,6 +1233,11 @@ class FiniteElementBase : public Subscriptor,
                                      */
     bool restriction_is_additive (const unsigned int index) const;
 
+                                    /**
+                                     * @name Support points and interpolation
+                                     * @{
+                                     */
+    
                                     /**
                                      * Return the support points of
                                      * the trial functions on the
@@ -1398,7 +1405,7 @@ class FiniteElementBase : public Subscriptor,
 
                                      /**
                                       * The function corresponding to
-                                      * the @p unit_support_point
+                                      * the unit_support_point()
                                       * function, but for faces. See
                                       * there for more information.
                                       */
@@ -1406,6 +1413,86 @@ class FiniteElementBase : public Subscriptor,
     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
@@ -1518,6 +1605,8 @@ class FiniteElementBase : public Subscriptor,
 
                                     /**
                                      * Exception
+                                     *
+                                     * @ingroup Exceptions
                                      */
     DeclException1 (ExcShapeFunctionNotPrimitive,
                    int,
@@ -1529,35 +1618,61 @@ class FiniteElementBase : public Subscriptor,
                    << "_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,
@@ -1566,6 +1681,7 @@ class FiniteElementBase : public Subscriptor,
                    << ", which is not reasonable in the present dimension.");
                                     /**
                                      * Exception
+                                     * @ingroup Exceptions
                                      */
     DeclException2 (ExcComponentIndexInvalid,
                    int, int,
@@ -1573,13 +1689,14 @@ class FiniteElementBase : public Subscriptor,
                    << ") 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
@@ -1847,7 +1964,21 @@ class FiniteElementBase : public Subscriptor,
                                      * 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
index 26b1e3f67ee407ac6fe651adafe4d7745f31eb55..644e4d414e5100d31b2ae6dee65ca4e8d9f1cec8 100644 (file)
@@ -14,7 +14,7 @@
 #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;
index 9891ee9130462e0a67875cf4eeccb5f249fa7c40..6c50fa1d29bbeb0830c2c868a3416cfa3f82f359 100644 (file)
@@ -16,7 +16,7 @@
 #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>
 
index ad92073bfc7cbc2123f5d935c208bb8c17e07ab7..2d8eb3e9ee87cbd0b1834f47b925ee0f2324109e 100644 (file)
 //    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>
index dd27943f9925cdc6eb4a1cd7ef0f725d73b182ae..b7bad5a0758bbc2868de2e80eb6dbcdcd640b139 100644 (file)
@@ -18,7 +18,7 @@
 #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>
@@ -1489,8 +1489,8 @@ struct TriaNumberCache<3> : public TriaNumberCache<2>
  *   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
index af2a9b799c3db42e84afe2757afd9b95aa7798fe..94cabd7d9a9e4b67184431750a2b57554cee3b37 100644 (file)
@@ -20,7 +20,7 @@
 #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>
 
index 1d92f1b4ae9dd6be09c9783c69a044afc3cbb0fe..a541cab9b3515ae5a361ee6a575160ed681ea16c 100644 (file)
@@ -20,7 +20,7 @@
 #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>
index 90a794414ce10549a9a7213cfb9150c69d74e4cd..1463402d32e22107ae02899c12f2820f6c1204ce 100644 (file)
@@ -437,6 +437,31 @@ FiniteElementBase<dim>::has_support_points () const
 
 
 
+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
@@ -475,6 +500,31 @@ FiniteElementBase<dim>::has_face_support_points () 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
@@ -488,6 +538,79 @@ 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
index 06c21aa2846d23d0a361868ab26fafd7fddcb52b..c57c9856105030b122917bd28e1ed844fe41749d 100644 (file)
@@ -12,7 +12,7 @@
 //---------------------------------------------------------------------------
 
 
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
 #include <fe/fe.h>
 
 
index e6d34d7f4296d2bffae037ba657535daa3a5cf34..888ea42df701f1ba2c0f109ea913fcd9464e902f 100644 (file)
@@ -18,7 +18,7 @@
 #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>
index 36a8ef0fce8f49f765363fba3dded2201e058cd0..dc5ce426034c7e7a7f38c21bffbc199b37ce9177 100644 (file)
@@ -18,7 +18,7 @@
 #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>
 
index 381b3f9879caa70deb3e142aecfcf2f00f0ab1cc..8fc92c3e9fd072ab6e34c539d2b85e738d98ba42 100644 (file)
@@ -20,7 +20,7 @@
 #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>
index 16d2abf23a5d85328d4bc11c8d455657140dc3d0..05dbb1f59ca5be213930d9fb894e96e897d625df 100644 (file)
@@ -20,7 +20,7 @@
 #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>
index 5660fb286aee4b55b72e65686ba87d2fb39aa1d8..7282cc0305929039b4a623e574ed423212165300 100644 (file)
@@ -18,7 +18,7 @@
 #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>
index 9bdfe55124aa2349bf51b47f272d578703b86323..70492b4b42a6e2267944f0bc0316f2db51c15843 100644 (file)
@@ -16,7 +16,7 @@
 
 #include "../tests.h"
 #include <base/logstream.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
 
 #include <fstream>
 #include <cstdlib>
index e00737ff25ef3be793e2c8227719f60f41c2109e..0d2dfd2993704eab97134461a4d1fd4f61996590 100644 (file)
@@ -16,7 +16,7 @@
 
 #include "../tests.h"
 #include <base/logstream.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
 
 #include <fstream>
 #include <cstdlib>
index 3073e50b835ed63e7e9090679724a4a9163aa93c..907f1b21aa96ca9eede2e00897bdb026104ee9f4 100644 (file)
@@ -2287,7 +2287,7 @@ template class HuyghensWave<2>;
 #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>
@@ -2807,7 +2807,7 @@ template class EvaluateHuyghensWave<2>;
 #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>
@@ -4463,7 +4463,7 @@ void SweepInfo::write_summary (const std::list<EvaluationBase<2>*> &eval_list,
 #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>
@@ -6516,7 +6516,7 @@ template class TimeStep<2>;
 #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>

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.