--- /dev/null
+/*---------------------------- fe_lib.system.h ---------------------------*/
+/* $Id$ */
+#ifndef __fe_system_H
+#define __fe_system_H
+/*---------------------------- fe_lib.system.h ---------------------------*/
+
+
+#include <fe/fe.h>
+
+
+
+template <int dim>
+class FESystem : public FiniteElement<dim> {
+ public:
+
+ /**
+ * Constructor. Take a finite element type
+ * and the number of elements you want to
+ * group together using this class.
+ *
+ * In fact, the object #fe# is not used,
+ * apart from getting the number of dofs
+ * per vertex, line, etc for that finite
+ * element class. For this, it would have
+ * been possible to use the #get_fe_data#
+ * function that each element has to
+ * provide. The correct way to write
+ * this constructor would therefore have
+ * been to specify it without the first
+ * argument and let the user specify the
+ * desired finite element by an explicit
+ * template argument list, like this:
+ * #AnyClass::f<int>()#. However, #C++#
+ * does not allow this call sequence for
+ * constructors, so we have to use the
+ * way as shown here, to let the compiler
+ * deduce the template argument itself.
+ *
+ * Obviously, the tenplate finite element
+ * class needs to be of the same dimension
+ * as is this object.
+ */
+ template <typename FE>
+ FESystem (const FE &fe, const unsigned int n_elements);
+
+ /**
+ * Destructor.
+ */
+ virtual ~FESystem ();
+
+ /**
+ * Return the value of the #i#th shape
+ * function at point #p# on the unit cell.
+ *
+ * For an element composed of #N#
+ * subelements, the first #N# shape
+ * functions refer to the zeroth shape
+ * function of the underlying object,
+ * the shape functions #N..2N-1# refer
+ * to the base shape function with
+ * number #1#, and so on. The #i# shape
+ * function therefore equals the
+ * #i/N# the shape function of the
+ * base object.
+ */
+ virtual double shape_value(const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Return the gradient of the #i#th shape
+ * function at point #p# on the unit cell.
+ *
+ * For the ordering of shape functions
+ * refer to the #shape_value# function.
+ */
+ virtual Tensor<1,dim> shape_grad(const unsigned int i,
+ const Point<dim>& p) const;
+
+ /**
+ * Return the tensor of second derivatives
+ * of the #i#th shape function at
+ * point #p# on the unit cell.
+ *
+ * For the ordering of shape functions
+ * refer to the #shape_value# function.
+ */
+ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ *
+ * For the ordering of shape functions
+ * refer to the #shape_value# function.
+ */
+ virtual void get_unit_support_points (vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ *
+ * For the ordering of shape functions
+ * refer to the #shape_value# function.
+ */
+ virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
+ const Boundary<dim> &boundary,
+ vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
+ const Boundary<dim> &boundary,
+ vector<Point<dim> > &support_points) const;
+
+ /**
+ * Fill the local mass matrix. The elements
+ * of this matrix are the integrals
+ * $\int_K \phi_i \phi_j dx$ over a given
+ * cell $K$. However, here only those
+ * elements of the matrix are set for which
+ * the shape functions $\phi_i$ and
+ * $\phi_j$ belong to the same subelement,
+ * i.e. the resulting matrix is a block
+ * matrix where each block is a diagonal
+ * matrix with diagonal values equal to
+ * the respective entry of the local mass
+ * matrix for the underlying finite element
+ * class. This definition of the mass
+ * matrix for systems of finite elements
+ * is consistent with the use of the matrix
+ * for the projection of initial values and
+ * the like, where the subelements are not
+ * coupled to each other. Also in most
+ * other cases you will not want the
+ * coupling terms to appear in the mass
+ * matrix.
+ *
+ * If the shape functions of this element
+ * were numbered such that the first
+ * numbers are for the shape functions of
+ * the first subelement, then those for
+ * the second subelement, and so on, then
+ * the mass matrix generated by this
+ * function would be a block diagonal
+ * matrix with each block being the mass
+ * matrix of the base finite element.
+ *
+ * Refer to the base class for more
+ * information on this function.
+ */
+ virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
+ const Boundary<dim> &boundary,
+ dFMatrix &local_mass_matrix) const;
+
+ /**
+ * Return the value of the #i#th shape
+ * function of the transformation mapping
+ * from unit cell to real cell. Since
+ * the transform functions are not
+ * touched when clustering several finite
+ * element objects together using this
+ * class, this function simply passes down
+ * the call to the respective function of
+ * the underlying element.
+ */
+ virtual double shape_value_transform (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Same as above: return gradient of the
+ * #i#th shape function for the mapping
+ * from unit to real cell.
+ */
+ virtual Tensor<1,dim> shape_grad_transform (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ *
+ * Since this function is only about the
+ * mapping from unit to real cell, it
+ * is not affected by putting several
+ * equal elements together, so this
+ * function simply passes down to the
+ * underlying object.
+ */
+ virtual void get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
+ const Boundary<dim> &boundary,
+ const vector<Point<dim-1> > &unit_points,
+ vector<double> &face_jacobi_determinants) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ *
+ * Since this function is only about the
+ * mapping from unit to real cell, it
+ * is not affected by putting several
+ * equal elements together, so this
+ * function simply passes down to the
+ * underlying object.
+ */
+ virtual void get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
+ const unsigned int subface_no,
+ const vector<Point<dim-1> > &unit_points,
+ vector<double> &face_jacobi_determinants) const;
+
+ /**
+ * Return the normal vectors to the
+ * face with number #face_no# of #cell#.
+ *
+ * Since this function is only about the
+ * mapping from unit to real cell, it
+ * is not affected by putting several
+ * equal elements together, so this
+ * function simply passes down to the
+ * underlying object.
+ *
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Boundary<dim> &boundary,
+ const vector<Point<dim-1> > &unit_points,
+ vector<Point<dim> > &normal_vectors) const;
+
+ /**
+ * Return the normal vectors to the
+ * subface with number #subface_no# of
+ * the face with number #face_no# of #cell#.
+ *
+ * Since this function is only about the
+ * mapping from unit to real cell, it
+ * is not affected by putting several
+ * equal elements together, so this
+ * function simply passes down to the
+ * underlying object.
+ *
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const vector<Point<dim-1> > &unit_points,
+ vector<Point<dim> > &normal_vectors) const;
+
+ private:
+ /**
+ * Pointer to an object of the underlying
+ * finite element class. This object is
+ * created by the constructor.
+ */
+ const FiniteElement<dim> *const base_element;
+
+ /**
+ * Number of subelements of this object.
+ * Since these objects may have
+ * subobjects themselves, this may be
+ * smaller than the total number of finite
+ * elements composed into this structure.
+ */
+ const unsigned int n_sub_elements;
+
+ /**
+ * Helper function used in the constructor:
+ * take a #FiniteElementData# object
+ * and return an object of the same type
+ * with the number of degrees of
+ * freedom per vertex, line, etc.
+ * multiplied by #n#. Don't touch the
+ * number of functions for the
+ * transformation from unit to real
+ * cell.
+ */
+ static FiniteElementData<dim> multiply_dof_numbers (const FiniteElementData<dim> &fe_data,
+ const unsigned int N);
+
+ /**
+ * This function is simply singled out of
+ * the constructor; it sets up the
+ * #restriction# and #prolongation#
+ * matrices. Since the operation of this
+ * function can be done without explicit
+ * knowledge of the data type of the
+ * underlying finite element class, we
+ * don't want to have this function in
+ * the general template definition in
+ * the #.h# file.
+ */
+ void initialize_matrices ();
+};
+
+
+
+
+
+/* ------------------------- template functions ------------------------- */
+
+template <int dim>
+template <typename FE>
+FESystem<dim>::FESystem (const FE &fe, const unsigned int n_elements) :
+ FiniteElement (multiply_dof_numbers(fe, n_elements)),
+ base_element (new FE()),
+ n_sub_elements (n_elements)
+{
+ base_element->subscribe ();
+ initialize_matrices ();
+};
+
+
+
+
+/*---------------------------- fe_lib.system.h ---------------------------*/
+/* end of #ifndef __fe_system_H */
+#endif
+/*---------------------------- fe_lib.system.h ---------------------------*/
--- /dev/null
+/* $Id$ */
+/* Copyright W. Bangerth, University of Heidelberg, 1990 */
+
+
+#include <fe/fe_system.h>
+
+
+
+
+
+template <int dim>
+FESystem<dim>::~FESystem () {
+ base_element->unsubscribe ();
+ delete base_element;
+};
+
+
+
+template <int dim>
+void FESystem<dim>::initialize_matrices () {
+ // distribute the matrices of the base
+ // finite element to the matrices of
+ // this object
+ for (unsigned int i=0; i<base_element->total_dofs; ++i)
+ for (unsigned int j=0; j<base_element->total_dofs; ++j)
+ for (unsigned int n=0; n<n_sub_elements; ++n)
+ // only fill diagonals of the blocks
+ {
+ for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
+ {
+ restriction[child] (i*n_sub_elements + n,
+ j*n_sub_elements + n)
+ = base_element->restriction[child](i,j);
+ prolongation[child] (i*n_sub_elements + n,
+ j*n_sub_elements + n)
+ = base_element->prolongation[child](i,j);
+ };
+
+ interface_constraints (i*n_sub_elements + n,
+ j*n_sub_elements + n)
+ = base_element->interface_constraints(i,j);
+ };
+};
+
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+FiniteElementData<1>
+FESystem<1>::multiply_dof_numbers (const FiniteElementData<1> &fe_data,
+ const unsigned int &N) {
+ return FiniteElementData<1> (fe_data.dofs_per_vertex * N,
+ fe_data.dofs_per_line * N,
+ fe_data.n_transform_functions);
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+FiniteElementData<2>
+FESystem<2>::multiply_dof_numbers (const FiniteElementData<2> &fe_data,
+ const unsigned int N) {
+ return FiniteElementData<2> (fe_data.dofs_per_vertex * N,
+ fe_data.dofs_per_line * N,
+ fe_data.dofs_per_quad * N,
+ fe_data.n_transform_functions);
+};
+
+#endif
+
+
+
+template <int dim>
+double FESystem<dim>::shape_value (const unsigned int i,
+ const Point<dim> &p) const {
+ Assert((i<total_dofs), ExcInvalidIndex(i));
+
+ return base_element->shape_value (i / n_sub_elements, p);
+};
+
+
+
+template <int dim>
+Tensor<1,dim>
+FESystem<dim>::shape_grad (const unsigned int i,
+ const Point<dim> &p) const {
+ Assert((i<total_dofs), ExcInvalidIndex(i));
+
+ return base_element->shape_grad (i / n_sub_elements, p);
+};
+
+
+
+template <int dim>
+Tensor<2,dim>
+FESystem<dim>::shape_grad_grad (const unsigned int i,
+ const Point<dim> &p) const {
+ Assert((i<total_dofs), ExcInvalidIndex(i));
+
+ return base_element->shape_grad_grad (i / n_sub_elements, p);
+};
+
+
+
+template <int dim>
+void FESystem<dim>::get_unit_support_points (vector<Point<dim> > &support_points) const {
+ Assert (support_points.size() == total_dofs,
+ ExcWrongFieldDimension (support_points.size(), total_dofs));
+
+ vector<Point<dim> > base_support_points (base_element->total_dofs);
+ base_element->get_unit_support_points (base_support_points);
+
+ for (unsigned int i=0; i<base_element->total_dofs; ++i)
+ for (unsigned int n=0; n<n_sub_elements; ++n)
+ support_points[i*n_sub_elements+n] = base_support_points[i];
+};
+
+
+
+template <int dim>
+void FESystem<dim>::get_support_points (const DoFHandler<dim>::cell_iterator &cell,
+ const Boundary<dim> &boundary,
+ vector<Point<dim> > &support_points) const {
+ Assert (support_points.size() == total_dofs,
+ ExcWrongFieldDimension (support_points.size(), total_dofs));
+
+ vector<Point<dim> > base_support_points (base_element->total_dofs);
+ base_element->get_support_points (cell, boundary, base_support_points);
+
+ for (unsigned int i=0; i<base_element->total_dofs; ++i)
+ for (unsigned int n=0; n<n_sub_elements; ++n)
+ support_points[i*n_sub_elements+n] = base_support_points[i];
+};
+
+
+
+template <int dim>
+void FESystem<dim>::get_face_support_points (const DoFHandler<dim>::face_iterator &face,
+ const Boundary<dim> &boundary,
+ vector<Point<dim> > &support_points) const {
+ Assert (support_points.size() == dofs_per_face,
+ ExcWrongFieldDimension (support_points.size(), dofs_per_face));
+
+ vector<Point<dim> > base_support_points (base_element->dofs_per_face);
+ base_element->get_face_support_points (face, boundary, base_support_points);
+
+ for (unsigned int i=0; i<base_element->dofs_per_face; ++i)
+ for (unsigned int n=0; n<n_sub_elements; ++n)
+ support_points[i*n_sub_elements+n] = base_support_points[i];
+};
+
+
+
+template <int dim>
+void FESystem<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
+ const Boundary<dim> &boundary,
+ dFMatrix &local_mass_matrix) const {
+ Assert (local_mass_matrix.n() == total_dofs,
+ ExcWrongFieldDimension(local_mass_matrix.n(),total_dofs));
+ Assert (local_mass_matrix.m() == total_dofs,
+ ExcWrongFieldDimension(local_mass_matrix.m(),total_dofs));
+
+ // first get the local mass matrix for
+ // the base object
+ dFMatrix base_mass_matrix (base_element->total_dofs,
+ base_element->total_dofs);
+ base_element->get_local_mass_matrix (cell, boundary, base_mass_matrix);
+
+
+ // now distribute it to the mass matrix
+ // of this object
+ for (unsigned int i=0; i<base_element->total_dofs; ++i)
+ for (unsigned int j=0; j<base_element->total_dofs; ++j)
+ for (unsigned int n=0; n<n_sub_elements; ++n)
+ // only fill diagonals of the blocks
+ local_mass_matrix (i*n_sub_elements + n,
+ j*n_sub_elements + n) = base_mass_matrix (i,j);
+};
+
+
+
+template <int dim>
+double FESystem<dim>::shape_value_transform (const unsigned int i,
+ const Point<dim> &p) const {
+ return base_element->shape_value_transform (i, p);
+};
+
+
+
+template <int dim>
+Tensor<1,dim> FESystem<dim>::shape_grad_transform (const unsigned int i,
+ const Point<dim> &p) const {
+ return base_element->shape_grad_transform (i, p);
+};
+
+
+
+template <int dim>
+void FESystem<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
+ const Boundary<dim> &boundary,
+ const vector<Point<dim-1> > &unit_points,
+ vector<double> &face_jacobi_determinants) const {
+ base_element->get_face_jacobians (face, boundary, unit_points, face_jacobi_determinants);
+};
+
+
+
+template <int dim>
+void FESystem<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
+ const unsigned int subface_no,
+ const vector<Point<dim-1> > &unit_points,
+ vector<double> &face_jacobi_determinants) const {
+ base_element->get_subface_jacobians (face, subface_no, unit_points, face_jacobi_determinants);
+};
+
+
+
+
+template <int dim>
+void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Boundary<dim> &boundary,
+ const vector<Point<dim-1> > &unit_points,
+ vector<Point<dim> > &normal_vectors) const {
+ base_element->get_normal_vectors (cell, face_no, boundary, unit_points, normal_vectors);
+};
+
+
+
+template <int dim>
+void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const vector<Point<dim-1> > &unit_points,
+ vector<Point<dim> > &normal_vectors) const {
+ base_element->get_normal_vectors (cell, face_no, subface_no, unit_points, normal_vectors);
+};
+
+
+
+
+
+// explicit instantiations
+template class FESystem<deal_II_dimension>;
+
+
+
+#include <fe/fe_lib.lagrange.h>
+FELinear<2> fl;
+FESystem<2> fe(fl, 2);
+
+