*
* <h3>Notes on the implementation of derived classes</h3>
*
+ * The following sections list the information to be provided by
+ * derived classes, depending on the space dimension. They are
+ * followed by a list of functions helping to generate these values.
+ *
* <h4>Finite elements in one dimension</h4>
*
* Finite elements in one dimension need only set the #restriction
* constraints that are entered more than once (as is necessary for the case
* above), it insists that the weights are exactly the same.
*
+ * <h4>Helper functions</h4>
+ *
+ * Construction of a finite element and computation of the matrices
+ * described above may be a tedious task, in particular if it has to
+ * be performed for several space dimensions. Therefore, some
+ * functions in FETools have been provided to help with these tasks.
+ *
+ * First, aready the basis of the shape function space may be
+ * difficult to implement for arbitrary order and dimension. On the
+ * other hand, if the @ref GlossNodes "node values" are given, then
+ * the duality relation between node functionals and basis functions
+ * defines the basis. As a result, the shape function space may be
+ * defined with arbitrary "raw" basis functions, such that the actual
+ * finite element basis is computed from linear combinations of
+ * them. The coefficients of these combinations are determined by the
+ * duality of node values.
+ *
+ * Using this matrix allows the construction of the basis of shape
+ * functions in two steps.
+ * <ol>
+ *
+ * <li>Define the space of shape functions using an arbitrary basis
+ * <i>w<sub>j</sub></i> and compute the matrix <i>M</i> of node
+ * functionals <i>N<sub>i</sub></i> applied to these basis functions,
+ * such that its entries are <i>m<sub>ij</sub> =
+ * N<sub>i</sub>(w<sub>j</sub>)</i>.
+ *
+ * <li>Compute the basis <i>v<sub>j</sub></i> of the finite element
+ * shape function space by applying <i>M<sup>-1</sup></i> to the basis
+ * <i>w<sub>j</sub></i>.
+ * </ol>
+ *
+ * The function computing the matrix <i>M</i> for you is
+ * FETools::compute_node_matrix(). It relies on the existence of
+ * #generalized_support_points and implementation of interpolate()
+ * with VectorSlice argument.
+ *
* @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1998, 2000, 2001, 2005
*/
template <int dim>
// basis functions
initialize_unit_support_points (deg);
initialize_node_matrix();
-
+
for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
this->prolongation[i].reinit (this->dofs_per_cell,
this->dofs_per_cell);
}
+template <int dim>
+void
+FE_RaviartThomasNodal<dim>::interpolate(
+ std::vector<double>&,
+ const std::vector<double>&) const
+{
+ Assert(false, ExcNotImplemented());
+}
+
+
+template <int dim>
+void
+FE_RaviartThomasNodal<dim>::interpolate(
+ std::vector<double>& local_dofs,
+ const std::vector<Vector<double> >& values,
+ unsigned int offset) const
+{
+ Assert (values.size() == generalized_support_points.size(),
+ ExcDimensionMismatch(values.size(), generalized_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()));
+
+ // First do interpolation on
+ // faces. There, the component
+ // evaluated depends on the face
+ // direction and orientation.
+ unsigned int fbase = 0;
+ unsigned int f=0;
+ for (;f<GeometryInfo<dim>::faces_per_cell;
+ ++f, fbase+=this->dofs_per_face)
+ {
+ for (unsigned int i=0;i<this->dofs_per_face;++i)
+ {
+ local_dofs[fbase+i] = values[fbase+i](offset+GeometryInfo<dim>::unit_normal_direction[f]);
+ }
+ }
+
+ // The remaining points form dim
+ // chunks, one for each component.
+ const unsigned int istep = (this->dofs_per_cell - fbase) / dim;
+ Assert ((this->dofs_per_cell - fbase) % dim == 0, ExcInternalError());
+
+ f = 0;
+ while (fbase < this->dofs_per_cell)
+ {
+ for (unsigned int i=0;i<istep;++i)
+ {
+ local_dofs[fbase+i] = values[fbase+i](offset+f);
+ }
+ fbase+=istep;
+ ++f;
+ }
+ Assert (fbase == this->dofs_per_cell, ExcInternalError());
+}
+
+
+
+template <int dim>
+void
+FE_RaviartThomasNodal<dim>::interpolate(
+ std::vector<double>& local_dofs,
+ const VectorSlice<const std::vector<std::vector<double> > >& values) const
+{
+ Assert (values.size() == this->n_components(),
+ ExcDimensionMismatch(values.size(), this->n_components()));
+ Assert (values[0].size() == generalized_support_points.size(),
+ ExcDimensionMismatch(values.size(), generalized_support_points.size()));
+ Assert (local_dofs.size() == this->dofs_per_cell,
+ ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+ // First do interpolation on
+ // faces. There, the component
+ // evaluated depends on the face
+ // direction and orientation.
+ unsigned int fbase = 0;
+ unsigned int f=0;
+ for (;f<GeometryInfo<dim>::faces_per_cell;
+ ++f, fbase+=this->dofs_per_face)
+ {
+ for (unsigned int i=0;i<this->dofs_per_face;++i)
+ {
+ local_dofs[fbase+i] = values[GeometryInfo<dim>::unit_normal_direction[f]][fbase+i];
+ }
+ }
+ // The remaining points form dim
+ // chunks, one for each component.
+ const unsigned int istep = (this->dofs_per_cell - fbase) / dim;
+ Assert ((this->dofs_per_cell - fbase) % dim == 0, ExcInternalError());
+
+ f = 0;
+ while (fbase < this->dofs_per_cell)
+ {
+ for (unsigned int i=0;i<istep;++i)
+ {
+ local_dofs[fbase+i] = values[f][fbase+i];
+ }
+ fbase+=istep;
+ ++f;
+ }
+ Assert (fbase == this->dofs_per_cell, ExcInternalError());
+}
+
+
+
template <int dim>
FiniteElement<dim> *
void
FE_RaviartThomasNodal<dim>::initialize_unit_support_points (const unsigned int deg)
{
- this->unit_support_points.resize (this->dofs_per_cell);
- this->unit_face_support_points.resize (this->dofs_per_face);
+ this->generalized_support_points.resize (this->dofs_per_cell);
+ this->generalized_face_support_points.resize (this->dofs_per_face);
unsigned int current = 0;
Assert (face_points.n_quadrature_points == this->dofs_per_face,
ExcInternalError());
for (unsigned int k=0;k<this->dofs_per_face;++k)
- this->unit_face_support_points[k] = face_points.point(k);
+ this->generalized_face_support_points[k] = face_points.point(k);
Quadrature<dim> faces = QProjector<dim>::project_to_all_faces(face_points);
for (unsigned int k=0;k<//faces.n_quadrature_points
this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;++k)
- this->unit_support_points[k] = faces.point(k);
+ this->generalized_support_points[k] = faces.point(k);
current = this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
}
((d==1) ? low : high),
((d==2) ? low : high));
for (unsigned int k=0;k<quadrature->n_quadrature_points;++k)
- this->unit_support_points[current++] = quadrature->point(k);
+ this->generalized_support_points[current++] = quadrature->point(k);
delete quadrature;
}
Assert (current == this->dofs_per_cell, ExcInternalError());
FE_RaviartThomasNodal<dim>::initialize_node_matrix ()
{
const unsigned int n_dofs = this->dofs_per_cell;
-
// We use an auxiliary matrix in
// this function. Therefore,
// inverse_node_matrix is still
// empty and shape_value_component
// returns the 'raw' shape values.
- FullMatrix<double> N(n_dofs, n_dofs);
-
- // The curent node functional index
- unsigned int current = 0;
-
- // For each face and all quadrature
- // points on it, the node value is
- // the normal component of the
- // shape function, possibly
- // pointing in negative direction.
- for (unsigned int face = 0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- for (unsigned int k=0;k<this->dofs_per_face;++k)
- {
- for (unsigned int i=0;i<n_dofs;++i)
- N(current,i) = this->shape_value_component(
- i, this->unit_support_points[current],
- GeometryInfo< dim >::unit_normal_direction[face]);
- ++current;
- }
- // Interior degrees of freedom in each direction
- const unsigned int n_cell = (n_dofs - current) / dim;
-
- for (unsigned int d=0;d<dim;++d)
- for (unsigned int k=0;k<n_cell;++k)
- {
- for (unsigned int i=0;i<n_dofs;++i)
- N(current,i) = this->shape_value_component(i, this->unit_support_points[current], d);
- ++current;
- }
- Assert (current == n_dofs, ExcInternalError());
+ FullMatrix<double> M(n_dofs, n_dofs);
+ FETools::compute_node_matrix(M, *this);
this->inverse_node_matrix.reinit(n_dofs, n_dofs);
- this->inverse_node_matrix.invert(N);
+ this->inverse_node_matrix.invert(M);
}