From: guido Date: Tue, 19 Jul 2005 09:43:00 +0000 (+0000) Subject: GeometryInfo moved; declaration of generalized support points and interpolation in... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=36bca6e07b4420c9980b3dfd9b636b9363a04185;p=dealii-svn.git GeometryInfo moved; declaration of generalized support points and interpolation in FiniteElementBase git-svn-id: https://svn.dealii.org/trunk@11174 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/data_out_base.h b/deal.II/base/include/base/data_out_base.h index fd8ac652e5..8e5a897dae 100644 --- a/deal.II/base/include/base/data_out_base.h +++ b/deal.II/base/include/base/data_out_base.h @@ -17,7 +17,7 @@ #include #include #include -#include +#include #include #include diff --git a/deal.II/base/include/base/geometry_info.h b/deal.II/base/include/base/geometry_info.h new file mode 100644 index 0000000000..30c6deb71a --- /dev/null +++ b/deal.II/base/include/base/geometry_info.h @@ -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 +#include +#include + + +template 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 + * for (i=0; i, + * at least if @p i is an unsigned int. + */ + 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 +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::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::vertices_per_cell; + + /** + * Number of lines on each face. + */ + static const unsigned int lines_per_face + = GeometryInfo::lines_per_cell; + + /** + * Number of quads on each face. + */ + static const unsigned int quads_per_face + = GeometryInfo::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::lines_per_cell + + GeometryInfo::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::quads_per_cell + + GeometryInfo::lines_per_cell); + + /** + * Number of hexahedra of a + * cell. + */ + static const unsigned int hexes_per_cell + = (2*GeometryInfo::hexes_per_cell + + GeometryInfo::quads_per_cell); + + /** + * List of numbers which denotes which + * face is opposite to a given face. In + * 1d, this list is {1,0}, in 2d + * {2, 3, 0, 1}, in 3d {1, + * 0, 4, 5, 2, 3}. + */ + 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 + * deal.II 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 + * 1 or -1, + * 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 face_orientation + * argument this function handles + * faces oriented in both, the + * standard and non-standard + * orientation. + * face_orientation + * defaults to true + * (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 + * vertexth vertex of + * line line, e.g. + * GeometryInfo<2>::line_to_cell_vertices(2,0)=3. + * + * 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 dim=2 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 + * vertexth vertex of + * face face, e.g. + * GeometryInfo<2>::face_to_cell_vertices(2,0)=3. + * + * Through the + * face_orientation + * argument this function handles + * faces oriented in both, the + * standard and non-standard + * orientation. + * face_orientation + * defaults to true + * (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 + * lineth line of face + * face, e.g. + * GeometryInfo<3>::face_to_cell_lines(3,1)=5. + * + * Through the + * face_orientation + * argument this function handles + * faces oriented in both, the + * standard and non-standard + * orientation. + * face_orientation + * defaults to true + * (standard orientation). + * + * This function is useful and + * implemented for + * dim=3, 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 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 + * GeometryInfo::children_per_cell. + * + * The order of child cells is + * described the documentation of + * the Triangulation class. + */ + static unsigned int child_cell_from_point (const Point &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 cell_to_child_coordinates (const Point &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 child_to_cell_coordinates (const Point &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 &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 +inline +Point +GeometryInfo::unit_cell_vertex (const unsigned int) +{ + Assert(false, ExcNotImplemented()); + + return Point (); +} + + + +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 +inline +unsigned int +GeometryInfo::child_cell_from_point (const Point &) +{ + Assert(false, ExcNotImplemented()); + + return 0; +} + + + +template +inline +Point +GeometryInfo::cell_to_child_coordinates (const Point &p, + const unsigned int child_index) +{ + Assert (child_index < GeometryInfo::children_per_cell, + ExcIndexRange (child_index, 0, GeometryInfo::children_per_cell)); + + return 2*p - unit_cell_vertex(child_index); +} + + + +template +inline +Point +GeometryInfo::child_to_cell_coordinates (const Point &p, + const unsigned int child_index) +{ + Assert (child_index < GeometryInfo::children_per_cell, + ExcIndexRange (child_index, 0, GeometryInfo::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 diff --git a/deal.II/deal.II/source/grid/geometry_info.cc b/deal.II/base/source/geometry_info.cc similarity index 78% rename from deal.II/deal.II/source/grid/geometry_info.cc rename to deal.II/base/source/geometry_info.cc index 7280afa11d..c9b1806caa 100644 --- a/deal.II/deal.II/source/grid/geometry_info.cc +++ b/deal.II/base/source/geometry_info.cc @@ -12,7 +12,7 @@ //--------------------------------------------------------------------------- -#include +#include @@ -27,6 +27,9 @@ template const unsigned int GeometryInfo::lines_per_cell; template const unsigned int GeometryInfo::quads_per_cell; template const unsigned int GeometryInfo::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::face_to_cell_lines (const unsigned int face, Assert (face::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 inline unsigned int @@ -277,4 +353,4 @@ GeometryInfo::face_to_cell_vertices (const unsigned int face, template class GeometryInfo<1>; template class GeometryInfo<2>; template class GeometryInfo<3>; - +template class GeometryInfo<4>; diff --git a/deal.II/base/source/quadrature.cc b/deal.II/base/source/quadrature.cc index e1d9f01ee7..08edb028e1 100644 --- a/deal.II/base/source/quadrature.cc +++ b/deal.II/base/source/quadrature.cc @@ -12,7 +12,7 @@ //--------------------------------------------------------------------------- -#include +#include #include #include #include diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index e755bc5547..4296c5b1cf 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -19,12 +19,14 @@ #include #include #include -#include +#include +#include #include #include #include #include +#include template 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 unit_face_support_point (const unsigned int index) const; + /** + * Return a support point vector + * for generalized interpolation. + */ + const std::vector > & + get_generalized_support_points () const; + + /** + * + */ + bool has_generalized_support_points () const; + + /** + * + */ + const std::vector > & + 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& local_dofs, + const std::vector& 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, offset 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& local_dofs, + const std::vector >& values, + unsigned int offset = 0) const; + + /** + * Interpolate a set of vector + * values, computed in the + * generalized support points. + */ + virtual void interpolate( + std::vector& local_dofs, + const VectorSlice > >& 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 > unit_face_support_points; - + + /** + * Support points used for + * interpolation functions of + * non-Lagrangian elements. + */ + std::vector > generalized_support_points; + + /** + * Face support points used for + * interpolation functions of + * non-Lagrangian elements. + */ + std::vector > generalized_face_support_points; + /** * For each shape function, give * a vector of bools (with size diff --git a/deal.II/deal.II/include/fe/fe_nedelec.h b/deal.II/deal.II/include/fe/fe_nedelec.h index 26b1e3f67e..644e4d414e 100644 --- a/deal.II/deal.II/include/fe/fe_nedelec.h +++ b/deal.II/deal.II/include/fe/fe_nedelec.h @@ -14,7 +14,7 @@ #define __deal2__fe_nedelec_h #include -#include +#include #include template class MappingQ; diff --git a/deal.II/deal.II/include/fe/fe_raviart_thomas.h b/deal.II/deal.II/include/fe/fe_raviart_thomas.h index 9891ee9130..6c50fa1d29 100644 --- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -16,7 +16,7 @@ #include #include #include -#include +#include #include #include diff --git a/deal.II/deal.II/include/grid/geometry_info.h b/deal.II/deal.II/include/grid/geometry_info.h index ad92073bfc..2d8eb3e9ee 100644 --- a/deal.II/deal.II/include/grid/geometry_info.h +++ b/deal.II/deal.II/include/grid/geometry_info.h @@ -10,731 +10,8 @@ // 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 -#include -#include - - -template 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 - * for (i=0; i, - * at least if @p i is an unsigned int. - */ - 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 -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::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::vertices_per_cell; - - /** - * Number of lines on each face. - */ - static const unsigned int lines_per_face - = GeometryInfo::lines_per_cell; - - /** - * Number of quads on each face. - */ - static const unsigned int quads_per_face - = GeometryInfo::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::lines_per_cell + - GeometryInfo::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::quads_per_cell + - GeometryInfo::lines_per_cell); - - /** - * Number of hexahedra of a - * cell. - */ - static const unsigned int hexes_per_cell - = (2*GeometryInfo::hexes_per_cell + - GeometryInfo::quads_per_cell); - - /** - * List of numbers which denotes which - * face is opposite to a given face. In - * 1d, this list is {1,0}, in 2d - * {2, 3, 0, 1}, in 3d {1, - * 0, 4, 5, 2, 3}. - */ - 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 - * deal.II 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 - * 1 or -1, - * 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 face_orientation - * argument this function handles - * faces oriented in both, the - * standard and non-standard - * orientation. - * face_orientation - * defaults to true - * (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 - * vertexth vertex of - * line line, e.g. - * GeometryInfo<2>::line_to_cell_vertices(2,0)=3. - * - * 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 dim=2 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 - * vertexth vertex of - * face face, e.g. - * GeometryInfo<2>::face_to_cell_vertices(2,0)=3. - * - * Through the - * face_orientation - * argument this function handles - * faces oriented in both, the - * standard and non-standard - * orientation. - * face_orientation - * defaults to true - * (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 - * lineth line of face - * face, e.g. - * GeometryInfo<3>::face_to_cell_lines(3,1)=5. - * - * Through the - * face_orientation - * argument this function handles - * faces oriented in both, the - * standard and non-standard - * orientation. - * face_orientation - * defaults to true - * (standard orientation). - * - * This function is useful and - * implemented for - * dim=3, 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 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 - * GeometryInfo::children_per_cell. - * - * The order of child cells is - * described the documentation of - * the Triangulation class. - */ - static unsigned int child_cell_from_point (const Point &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 cell_to_child_coordinates (const Point &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 child_to_cell_coordinates (const Point &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 &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 -inline -Point -GeometryInfo::unit_cell_vertex (const unsigned int) -{ - Assert(false, ExcNotImplemented()); - - return Point (); -} - - - -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 -inline -unsigned int -GeometryInfo::child_cell_from_point (const Point &) -{ - Assert(false, ExcNotImplemented()); - - return 0; -} - - - -template -inline -Point -GeometryInfo::cell_to_child_coordinates (const Point &p, - const unsigned int child_index) -{ - Assert (child_index < GeometryInfo::children_per_cell, - ExcIndexRange (child_index, 0, GeometryInfo::children_per_cell)); - - return 2*p - unit_cell_vertex(child_index); -} - - - -template -inline -Point -GeometryInfo::child_to_cell_coordinates (const Point &p, - const unsigned int child_index) -{ - Assert (child_index < GeometryInfo::children_per_cell, - ExcIndexRange (child_index, 0, GeometryInfo::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 diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index dd27943f99..b7bad5a075 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -18,7 +18,7 @@ #include #include #include -#include +#include #include #include @@ -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 class Triangulation : public Subscriptor diff --git a/deal.II/deal.II/include/grid/tria_accessor.templates.h b/deal.II/deal.II/include/grid/tria_accessor.templates.h index af2a9b799c..94cabd7d9a 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.templates.h +++ b/deal.II/deal.II/include/grid/tria_accessor.templates.h @@ -20,7 +20,7 @@ #include #include #include -#include +#include #include diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc index 1d92f1b4ae..a541cab9b3 100644 --- a/deal.II/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/deal.II/source/dofs/dof_handler.cc @@ -20,7 +20,7 @@ #include #include #include -#include +#include #include #include diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 90a794414c..1463402d32 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -437,6 +437,31 @@ FiniteElementBase::has_support_points () const +template +const std::vector > & +FiniteElementBase::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 +bool +FiniteElementBase::has_generalized_support_points () const +{ + return (generalized_support_points.size() != 0); +} + + + template Point FiniteElementBase::unit_support_point (const unsigned index) const @@ -475,6 +500,31 @@ FiniteElementBase::has_face_support_points () const +template +const std::vector > & +FiniteElementBase::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 +bool +FiniteElementBase::has_generalized_face_support_points () const +{ + return (generalized_face_support_points.size() != 0); +} + + + template Point FiniteElementBase::unit_face_support_point (const unsigned index) const @@ -488,6 +538,79 @@ FiniteElementBase::unit_face_support_point (const unsigned index) const +template +void +FiniteElementBase::interpolate( + std::vector& local_dofs, + const std::vector& 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 +void +FiniteElementBase::interpolate( + std::vector& local_dofs, + const std::vector >& 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;idofs_per_cell;++i) + { + const std::pair index + = this->system_to_component_index(i); + local_dofs[i] = values[index.second](offset+index.first); + } +} + + + + +template +void +FiniteElementBase::interpolate( + std::vector& local_dofs, + const VectorSlice > >& 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;idofs_per_cell;++i) + { + const std::pair index + = this->system_to_component_index(i); + local_dofs[i] = values[index.first][index.second]; + } +} + + + + template unsigned int FiniteElementBase::memory_consumption () const diff --git a/deal.II/deal.II/source/fe/fe_data.cc b/deal.II/deal.II/source/fe/fe_data.cc index 06c21aa284..c57c985610 100644 --- a/deal.II/deal.II/source/fe/fe_data.cc +++ b/deal.II/deal.II/source/fe/fe_data.cc @@ -12,7 +12,7 @@ //--------------------------------------------------------------------------- -#include +#include #include diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index e6d34d7f42..888ea42df7 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -18,7 +18,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc index 36a8ef0fce..dc5ce42603 100644 --- a/deal.II/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/deal.II/source/grid/tria_accessor.cc @@ -18,7 +18,7 @@ #include #include #include -#include +#include #include #include diff --git a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc index 381b3f9879..8fc92c3e9f 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc @@ -20,7 +20,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index 16d2abf23a..05dbb1f59c 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -20,7 +20,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index 5660fb286a..7282cc0305 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -18,7 +18,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/tests/bits/geometry_info_1.cc b/tests/bits/geometry_info_1.cc index 9bdfe55124..70492b4b42 100644 --- a/tests/bits/geometry_info_1.cc +++ b/tests/bits/geometry_info_1.cc @@ -16,7 +16,7 @@ #include "../tests.h" #include -#include +#include #include #include diff --git a/tests/bits/geometry_info_2.cc b/tests/bits/geometry_info_2.cc index e00737ff25..0d2dfd2993 100644 --- a/tests/bits/geometry_info_2.cc +++ b/tests/bits/geometry_info_2.cc @@ -16,7 +16,7 @@ #include "../tests.h" #include -#include +#include #include #include diff --git a/tests/deal.II/wave-test-3.cc b/tests/deal.II/wave-test-3.cc index 3073e50b83..907f1b21aa 100644 --- a/tests/deal.II/wave-test-3.cc +++ b/tests/deal.II/wave-test-3.cc @@ -2287,7 +2287,7 @@ template class HuyghensWave<2>; #include #include #include -#include +#include #include #include #include @@ -2807,7 +2807,7 @@ template class EvaluateHuyghensWave<2>; #include #include #include -#include +#include #include #include @@ -4463,7 +4463,7 @@ void SweepInfo::write_summary (const std::list*> &eval_list, #include #include #include -#include +#include #include #include #include @@ -6516,7 +6516,7 @@ template class TimeStep<2>; #include #include #include -#include +#include #include #include