template <int dim> class Quadrature;
+//TODO: Add access to mapping values to FEValuesBase
/*!@addtogroup febase */
/*@{*/
/**
- * Contains all data vectors for @p FEValues.
- * This class has been extracted from @p FEValuesBase to be handed
- * over to the fill functions of @p Mapping and
- * @p FiniteElement.
+ * Contains all data vectors for FEValues.
+ * This class has been extracted from FEValuesBase to be handed
+ * over to the fill functions of Mapping and
+ * FiniteElement.
*
* @note All data fields are public, but this is not
- * critical, because access to this object is private in @p FEValues.
+ * critical, because access to this object is private in FEValues.
*
* @author Guido Kanschat, 2000
*/
* shape function number equals
* the row number. Otherwise, use
* the
- * @p shape_function_to_row_table
+ * #shape_function_to_row_table
* array to get at the first row
* that belongs to this
* particular shape function, and
* navigate among all the rows
* for this shape function using
* the
- * @p FiniteElement::get_nonzero_components
+ * FiniteElement::get_nonzero_components()
* function which tells us which
* components are non-zero and
* thus have a row in the array
* Storage type for
* gradients. The layout of data
* is the same as for the
- * ShapeVector data type.
+ * #ShapeVector data type.
*/
typedef std::vector<std::vector<Tensor<1,dim> > > GradientVector;
* times the Jacobi determinant
* at the quadrature points. This
* function is reset each time
- * @p reinit is called. The
+ * reinit() is called. The
* Jacobi determinant is actually
* the reciprocal value of the
* Jacobi matrices stored in this
/**
* Array of quadrature points. This array
- * is set up upon calling @p reinit and
+ * is set up upon calling reinit() and
* contains the quadrature points on the
* real element, rather than on the
* reference element.
/**
* Indicate the first row which a
* given shape function occupies
- * in the @p ShapeVector,
- * @p ShapeGradient, etc
+ * in the #shape_values,
+ * #shape_gradients and
+ * #shape_2nd_derivatives
* arrays. If all shape functions
* are primitive, then this is
* the identity mapping. If, on
* follows: we allocate one row
* for each non-zero component as
* indicated by the
- * <tt>FiniteElement::get_nonzero_components()</tt>
+ * FiniteElement::get_nonzero_components()
* function, and the rows are in
* ascending order exactly those
* non-zero components.
/**
* Original update flags handed
* to the constructor of
- * @p FEValues.
+ * FEValues.
*/
UpdateFlags update_flags;
};
/**
* @brief Common features of <tt>FEValues*</tt> classes.
*
- * <tt>FEValues*</tt> objects are programming interfaces to finite element
- * and mapping classes on the one hand side, to cells and quadrature
- * rules on the other side. The reason for their existence is possible
- * optimization. Depending on the type of finite element and mapping,
- * some values can be computed once on the unit cell. Others must be
- * computed on each cell, but maybe computation of several values at
- * the same time offers ways for optimization. Since this interlay may
- * be complex and depends on the actual finite element, it cannot be
- * left to the applications programmer.
+ * FEValues, FEFaceValues and FESubfaceValues objects are programming
+ * interfaces to finite element and mapping classes on the one hand
+ * side, to cells and quadrature rules on the other side. The reason
+ * for their existence is possible optimization. Depending on the type
+ * of finite element and mapping, some values can be computed once on
+ * the unit cell. Others must be computed on each cell, but maybe
+ * computation of several values at the same time offers ways for
+ * optimization. Since this interlay may be complex and depends on the
+ * actual finite element, it cannot be left to the applications
+ * programmer.
*
- * <tt>FEValues*</tt> provides only data handling: computations are left to
- * objects of type Mapping and FiniteElement. These
- * provide functions <tt>get_*_data</tt> and <tt>fill_*_values</tt> which are
- * called by the constructor and <tt>reinit</tt> functions of
- * <tt>FEValues*</tt>, respectively.
+ * FEValues, FEFaceValues and FESubfaceValues provide only data
+ * handling: computations are left to objects of type Mapping and
+ * FiniteElement. These provide functions <tt>get_*_data</tt> and
+ * <tt>fill_*_values</tt> which are called by the constructor and
+ * <tt>reinit</tt> functions of <tt>FEValues*</tt>, respectively.
*
* @sect3{FEValuesBaseGeneral General usage}
*
void get_function_values (const InputVector &fe_function,
std::vector<Vector<number> > &values) const;
+ /**
+ * Generate function values from
+ * an arbitrary vector.
+ *
+ * This function offers the
+ * possibility to extract
+ * function values in quadrature
+ * points from vectors not
+ * corresponding to a whole
+ * discretization.
+ *
+ * You may want to use this
+ * function, if you want to
+ * access just a single block
+ * from a BlockVector, if you
+ * have a multi-level vector or
+ * if you already have a local
+ * representation of your finite
+ * element data.
+ */
+ template <class InputVector, typename number>
+ void get_function_values (const InputVector& fe_function,
+ const std::vector<unsigned int>& indices,
+ std::vector<number>& values) const;
+
+ /**
+ * Generate vector function
+ * values from an arbitrary
+ * vector.
+ *
+ * This function offers the
+ * possibility to extract
+ * function values in quadrature
+ * points from vectors not
+ * corresponding to a whole
+ * discretization.
+ *
+ * The length of the vector
+ * <tt>indices</tt> may even be a
+ * multiple of the number of dofs
+ * per cell. Then, the vectors in
+ * <tt>value</tt> should allow
+ * for the same multiple of the
+ * components of the finite
+ * element.
+ *
+ * You may want to use this
+ * function, if you want to
+ * access just a single block
+ * from a BlockVector, if you
+ * have a multi-level vector or
+ * if you already have a local
+ * representation of your finite
+ * element data.
+ *
+ * Since this function allows for
+ * fairly general combinations of
+ * argument sizes, be aware that
+ * the checks on the arguments
+ * may not detect errors.
+ */
+ template <class InputVector, typename number>
+ void get_function_values (const InputVector& fe_function,
+ const std::vector<unsigned int>& indices,
+ std::vector<Vector<number> >& values) const;
+
/**
* Compute the gradient of the
* @p ith shape function at the
#include <grid/tria_accessor.h>
#include <grid/tria_boundary.h>
#include <dofs/dof_accessor.h>
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+#include <lac/petsc_vector.h>
#include <fe/mapping_q1.h>
#include <fe/fe_values.h>
#include <fe/fe.h>
template <int dim>
template <class InputVector, typename number>
-void FEValuesBase<dim>::get_function_values (const InputVector &fe_function,
- std::vector<number> &values) const
+void FEValuesBase<dim>::get_function_values (
+ const InputVector &fe_function,
+ std::vector<number> &values) const
{
Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
Assert (fe->n_components() == 1,
template <int dim>
template <class InputVector, typename number>
-void FEValuesBase<dim>::get_function_values (const InputVector &fe_function,
- std::vector<Vector<number> > &values) const
+void FEValuesBase<dim>::get_function_values (
+ const InputVector& fe_function,
+ const std::vector<unsigned int>& indices,
+ std::vector<number> &values) const
{
- Assert (n_quadrature_points == values.size(),
+ Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
+ // This function fills a single
+ // component only
+ Assert (fe->n_components() == 1,
+ ExcWrongNoOfComponents());
+ // One index for each dof
+ Assert (indices.size() == dofs_per_cell,
+ ExcDimensionMismatch(indices.size(), dofs_per_cell));
+ // This vector has one entry for
+ // each quadrature point
+ Assert (values.size() == n_quadrature_points,
+ ExcWrongVectorSize(values.size(), n_quadrature_points));
+
+ // initialize with zero
+ std::fill_n (values.begin(), n_quadrature_points, 0);
+
+ // add up contributions of trial
+ // functions. note that here we
+ // deal with scalar finite
+ // elements, so no need to check
+ // for non-primitivity of shape
+ // functions
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+ values[point] += (fe_function(indices[shape_func]) *
+ this->shape_value(shape_func, point));
+}
+
+
+
+template <int dim>
+template <class InputVector, typename number>
+void FEValuesBase<dim>::get_function_values (
+ const InputVector& fe_function,
+ std::vector<Vector<number> >& values) const
+{
+//TODO: Find out how to do this assertion.
+ // This vector must correspond to a
+ // complete discretization
+// Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(),
+// ExcWrongVectorSize(fe_function.size(),
+// present_cell->get_dof_handler().n_dofs()));
+ // One entry per quadrature point
+ Assert (values.size() == n_quadrature_points,
ExcWrongVectorSize(values.size(), n_quadrature_points));
const unsigned int n_components = fe->n_components();
+ // Assert that we can write all
+ // components into the result
+ // vectors
for (unsigned i=0;i<values.size();++i)
Assert (values[i].size() == n_components, ExcWrongNoOfComponents());
-
+
Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(),
ExcWrongVectorSize(fe_function.size(), present_cell->n_dofs_for_dof_handler()));
+template <int dim>
+template <class InputVector, typename number>
+void FEValuesBase<dim>::get_function_values (
+ const InputVector& fe_function,
+ const std::vector<unsigned int>& indices,
+ std::vector<Vector<number> >& values) const
+{
+ // One value per quadrature point
+ Assert (n_quadrature_points == values.size(),
+ ExcWrongVectorSize(values.size(), n_quadrature_points));
+
+ const unsigned int n_components = fe->n_components();
+
+ // Size of indices must be a
+ // multiple of dofs_per_cell such
+ // that an integer number of
+ // function values is generated in
+ // each point.
+ Assert (indices.size() % dofs_per_cell == 0,
+ ExcNotMultiple(indices.size(), dofs_per_cell));
+
+ // The number of components of the
+ // result may be a multiple of the
+ // number of components of the
+ // finite element
+ const unsigned int result_components = indices.size() / dofs_per_cell;
+
+ for (unsigned i=0;i<values.size();++i)
+ Assert (values[i].size() == result_components,
+ ExcDimensionMismatch(values[i].size(), result_components));
+
+ // If the result has more
+ // components than the finite
+ // element, we need this number for
+ // loops filling all components
+ const unsigned int component_multiple = result_components / n_components;
+
+ Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
+
+ // initialize with zero
+ for (unsigned i=0;i<values.size();++i)
+ std::fill_n (values[i].begin(), values[i].size(), 0);
+
+ // add up contributions of trial
+ // functions. now check whether the
+ // shape function is primitive or
+ // not. if it is, then set its only
+ // non-zero component, otherwise
+ // loop over components
+ for (unsigned int mc = 0; mc < component_multiple; ++mc)
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+ if (fe->is_primitive(shape_func))
+ values[point](fe->system_to_component_index(shape_func).first
+ +mc * n_components)
+ += (fe_function(indices[shape_func+mc*dofs_per_cell])
+ * shape_value(shape_func, point));
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ values[point](c) += (fe_function(indices[shape_func]) *
+ shape_value_component(shape_func, point, c));
+}
+
+
+
template <int dim>
const std::vector<Point<dim> > &
FEValuesBase<dim>::get_quadrature_points () const
// passed to the constructor and
// used by the DoFHandler used by
// this cell, are the same
- Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
- static_cast<const FiniteElementData<dim>&>(cell->get_fe()),
- typename FEValuesBase<dim>::ExcFEDontMatch());
+
+//TODO: This was documented out ith the repository. Why?
+
+// Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
+// static_cast<const FiniteElementData<dim>&>(cell->get_fe()),
+// typename FEValuesBase<dim>::ExcFEDontMatch());
// set new cell. auto_ptr will take
// care that old object gets
// passed to the constructor and
// used by the DoFHandler used by
// this cell, are the same
- Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
- static_cast<const FiniteElementData<dim>&>(cell->get_dof_handler().get_fe()),
- typename FEValuesBase<dim>::ExcFEDontMatch());
+// Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
+// static_cast<const FiniteElementData<dim>&>(cell->get_dof_handler().get_fe()),
+// typename FEValuesBase<dim>::ExcFEDontMatch());
Assert (face_no < GeometryInfo<dim>::faces_per_cell,
ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
// passed to the constructor and
// used by the DoFHandler used by
// this cell, are the same
- Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
- static_cast<const FiniteElementData<dim>&>(cell->get_dof_handler().get_fe()),
- typename FEValuesBase<dim>::ExcFEDontMatch());
+// Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
+// static_cast<const FiniteElementData<dim>&>(cell->get_dof_handler().get_fe()),
+// typename FEValuesBase<dim>::ExcFEDontMatch());
Assert (face_no < GeometryInfo<dim>::faces_per_cell,
ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
Assert (subface_no < GeometryInfo<dim>::subfaces_per_face,
//-----------------------------------------------------------------------------
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<Vector<double> >
-(const Vector<double>&,
- std::vector<double>&) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<Vector<float> >
-(const Vector<float>&,
- std::vector<double>&) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<Vector<float> >
-(const Vector<float>&,
- std::vector<float>&) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<BlockVector<double> >
-(const BlockVector<double>&,
- std::vector<double>&) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<BlockVector<float> >
-(const BlockVector<float>&,
- std::vector<double>&) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<BlockVector<float> >
-(const BlockVector<float>&,
- std::vector<float>&) const;
-
-#ifdef DEAL_II_USE_PETSC
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<PETScWrappers::Vector>
-(const PETScWrappers::Vector &,
- std::vector<double>&) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<PETScWrappers::Vector>
-(const PETScWrappers::Vector &,
- std::vector<float>&) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<PETScWrappers::BlockVector>
-(const PETScWrappers::BlockVector &,
- std::vector<double>&) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<PETScWrappers::BlockVector>
-(const PETScWrappers::BlockVector &,
- std::vector<float>&) const;
-#endif
-
-//-----------------------------------------------------------------------------
-
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<Vector<double> >
-(const Vector<double> &,
- std::vector<Vector<double> > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<Vector<float> >
-(const Vector<float> &,
- std::vector<Vector<double> > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<Vector<float> >
-(const Vector<float> &,
- std::vector<Vector<float> > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<BlockVector<double> >
-(const BlockVector<double> &,
- std::vector<Vector<double> > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<BlockVector<float> >
-(const BlockVector<float> &,
- std::vector<Vector<double> > &) const;
-
-#ifdef DEAL_II_USE_PETSC
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<PETScWrappers::Vector>
-(const PETScWrappers::Vector &,
- std::vector<Vector<double> > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_values<PETScWrappers::BlockVector>
-(const PETScWrappers::BlockVector &,
- std::vector<Vector<double> > &) const;
-#endif
-
-//-----------------------------------------------------------------------------
-
-template
-void FEValuesBase<deal_II_dimension>::get_function_grads<Vector<double> >
-(const Vector<double> &,
- std::vector<Tensor<1,deal_II_dimension> > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_grads<Vector<float> >
-(const Vector<float> &,
- std::vector<Tensor<1,deal_II_dimension> > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_grads<BlockVector<double> >
-(const BlockVector<double> &,
- std::vector<Tensor<1,deal_II_dimension> > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_grads<BlockVector<float> >
-(const BlockVector<float> &,
- std::vector<Tensor<1,deal_II_dimension> > &) const;
-
-#ifdef DEAL_II_USE_PETSC
-template
-void FEValuesBase<deal_II_dimension>::get_function_grads<PETScWrappers::Vector>
-(const PETScWrappers::Vector &,
- std::vector<Tensor<1,deal_II_dimension> > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_grads<PETScWrappers::BlockVector>
-(const PETScWrappers::BlockVector &,
- std::vector<Tensor<1,deal_II_dimension> > &) const;
-#endif
-
-//-----------------------------------------------------------------------------
-
-template
-void FEValuesBase<deal_II_dimension>::get_function_grads<Vector<double> >
-(const Vector<double> &,
- std::vector<std::vector<Tensor<1,deal_II_dimension> > > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_grads<Vector<float> >
-(const Vector<float> &,
- std::vector<std::vector<Tensor<1,deal_II_dimension> > > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_grads<BlockVector<double> >
-(const BlockVector<double> &,
- std::vector<std::vector<Tensor<1,deal_II_dimension> > > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_grads<BlockVector<float> >
-(const BlockVector<float> &,
- std::vector<std::vector<Tensor<1,deal_II_dimension> > > &) const;
-
-#ifdef DEAL_II_USE_PETSC
-template
-void FEValuesBase<deal_II_dimension>::get_function_grads<PETScWrappers::Vector>
-(const PETScWrappers::Vector &,
- std::vector<std::vector<Tensor<1,deal_II_dimension> > > &) const;
-
-template
-void FEValuesBase<deal_II_dimension>::get_function_grads<PETScWrappers::BlockVector>
-(const PETScWrappers::BlockVector &,
- std::vector<std::vector<Tensor<1,deal_II_dimension> > > &) const;
-#endif
-
-//-----------------------------------------------------------------------------
+// Instantiations are in a different file using the macro IN for the vector type.
+// This way, we avoid code reduplication
-template
-void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<Vector<double> >
-(const Vector<double> &,
- std::vector<Tensor<2,deal_II_dimension> > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<Vector<float> >
-(const Vector<float> &,
- std::vector<Tensor<2,deal_II_dimension> > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<BlockVector<double> >
-(const BlockVector<double> &,
- std::vector<Tensor<2,deal_II_dimension> > &) const;
+#define IN Vector<double>
+#include "fe_values.instance.h"
+#undef IN
-#ifdef DEAL_II_USE_PETSC
-template
-void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<PETScWrappers::Vector>
-(const PETScWrappers::Vector &,
- std::vector<Tensor<2,deal_II_dimension> > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<PETScWrappers::BlockVector>
-(const PETScWrappers::BlockVector &,
- std::vector<Tensor<2,deal_II_dimension> > &) const;
-#endif
+#define IN BlockVector<double>
+#include "fe_values.instance.h"
+#undef IN
-//-----------------------------------------------------------------------------
+#define IN Vector<float>
+#include "fe_values.instance.h"
+#undef IN
-template
-void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<Vector<double> >
-(const Vector<double> &,
- std::vector<std::vector<Tensor<2,deal_II_dimension> > > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<Vector<float> >
-(const Vector<float> &,
- std::vector<std::vector<Tensor<2,deal_II_dimension> > > &) const;
-template
-void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<BlockVector<double> >
-(const BlockVector<double> &,
- std::vector<std::vector<Tensor<2,deal_II_dimension> > > &) const;
+#define IN BlockVector<float>
+#include "fe_values.instance.h"
+#undef IN
#ifdef DEAL_II_USE_PETSC
-template
-void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<PETScWrappers::Vector>
-(const PETScWrappers::Vector &,
- std::vector<std::vector<Tensor<2,deal_II_dimension> > > &) const;
-
-template
-void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<PETScWrappers::BlockVector>
-(const PETScWrappers::BlockVector &,
- std::vector<std::vector<Tensor<2,deal_II_dimension> > > &) const;
+#define IN PETScWrappers::Vector
+#include "fe_values.instance.h"
+#undef IN
#endif
projects / dealii-svn.git / commitdiff