#include <deal.II/base/smartpointer.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/shape_info.h>
-#include <deal.II/matrix_free/mapping_fe_evaluation.h>
+#include <deal.II/matrix_free/mapping_data_on_the_fly.h>
DEAL_II_NAMESPACE_OPEN
+
+
+// forward declarations
namespace parallel
{
namespace distributed
template <typename> class Vector;
}
}
-
-
-
-// forward declarations
namespace internal
{
DeclException0 (ExcAccessToUninitializedField);
}
+template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
+ int n_components_ = 1, typename Number = double > class FEEvaluation;
/**
*/
//@{
/**
- * Initializes the operation pointer to the current cell. Unlike the
- * FEValues::reinit function, where the information related to a particular
+ * Initializes the operation pointer to the current cell. Unlike the reinit
+ * functions taking a cell iterator as argument below and the
+ * FEValues::reinit() methods, where the information related to a particular
* cell is generated in the reinit call, this function is very cheap since
* all data is pre-computed in @p matrix_free, and only a few indices have
* to be set appropriately.
*/
void reinit (const unsigned int cell);
+ /**
+ * Initialize the data to the current cell using a TriaIterator object as
+ * usual in FEValues. The argument is either of type
+ * DoFHandler::active_cell_iterator or DoFHandler::level_cell_iterator. This
+ * option is only available if the FEEvaluation object was created with a
+ * finite element, quadrature formula and correct update flags and
+ * <b>without</b> a MatrixFree object. This initialization method loses the
+ * ability to use vectorization, see also the description of the
+ * FEEvaluation class. When this reinit method is used, FEEvaluation can
+ * also read from vectors (but less efficient than with data coming from
+ * MatrixFree).
+ */
+ template <class DH, bool level_dof_access>
+ void reinit (const TriaIterator<DoFCellAccessor<DH,level_dof_access> > &cell);
+
+ /**
+ * Initialize the data to the current cell using a TriaIterator object as
+ * usual in FEValues. This option is only available if the FEEvaluation
+ * object was created with a finite element, quadrature formula and correct
+ * update flags and <b>without</b> a MatrixFree object. This initialization
+ * method loses the ability to use vectorization, see also the description
+ * of the FEEvaluation class. When this reinit method is used, FEEvaluation
+ * can <b>not</b> read from vectors because no DoFHandler information is
+ * available.
+ */
+ void reinit (const typename Triangulation<dim>::cell_iterator &cell);
+
/**
* For the transformation information stored in MappingInfo, this function
* returns the index which belongs to the current cell as specified in @p
* solution of linear systems, the temporary solution should always have
* homogeneous constraints and this method is the correct one.
*
- * If the class was constructed through a MappingFEEvaluation object, only
- * one single cell is used by this class and this function extracts the
- * values of the underlying components on this cell. This call is slower
- * than the ones done through a MatrixFree object and lead to a structure
- * that does not effectively use vectorization in the evaluate routines
- * based on these values (instead, VectorizedArray<Number>::n_array_elements
- * same copies are worked on).
+ * If this class was constructed without a MatrixFree object and the
+ * information is acquired on the fly through a
+ * DoFHandler<dim>::cell_iterator, only one single cell is used by this
+ * class and this function extracts the values of the underlying components
+ * on the given cell. This call is slower than the ones done through a
+ * MatrixFree object and lead to a structure that does not effectively use
+ * vectorization in the evaluate routines based on these values (instead,
+ * VectorizedArray::n_array_elements same copies are worked on).
*/
template <typename VectorType>
void read_dof_values (const VectorType &src);
* as MatrixFree can only handle homogeneous constraints. Note that if
* vectorization is enabled, the DoF values for several cells are set.
*
- * If the class was constructed through a MappingFEEvaluation object, only
- * one single cell is used by this class and this function extracts the
- * values of the underlying components on this cell. This call is slower
- * than the ones done through a MatrixFree object and lead to a structure
- * that does not effectively use vectorization in the evaluate routines
- * based on these values (instead, VectorizedArray<Number>::n_array_elements
- * same copies are worked on). In that case, no constraints can be
- * processed as these are not available here.
+ * If this class was constructed without a MatrixFree object and the
+ * information is acquired on the fly through a
+ * DoFHandler<dim>::cell_iterator, only one single cell is used by this
+ * class and this function extracts the values of the underlying components
+ * on the given cell. This call is slower than the ones done through a
+ * MatrixFree object and lead to a structure that does not effectively use
+ * vectorization in the evaluate routines based on these values (instead,
+ * VectorizedArray::n_array_elements same copies are worked on).
*/
template <typename VectorType>
void read_dof_values_plain (const VectorType &src);
* function ConstraintMatrix::distribute_local_to_global. If vectorization
* is enabled, the DoF values for several cells are used.
*
- * If the class was constructed through a MappingFEEvaluation object, only
- * one single cell is used by this class and this function extracts the
- * values of the underlying components on this cell. This call is slower
- * than the ones done through a MatrixFree object and lead to a structure
- * that does not effectively use vectorization in the evaluate routines
- * based on these values (instead, VectorizedArray<Number>::n_array_elements
- * same copies are worked on). In that case, no constraints can be
- * processed as these are not available here.
+ * If this class was constructed without a MatrixFree object and the
+ * information is acquired on the fly through a
+ * DoFHandler<dim>::cell_iterator, only one single cell is used by this
+ * class and this function extracts the values of the underlying components
+ * on the given cell. This call is slower than the ones done through a
+ * MatrixFree object and lead to a structure that does not effectively use
+ * vectorization in the evaluate routines based on these values (instead,
+ * VectorizedArray::n_array_elements same copies are worked on).
*/
template<typename VectorType>
void distribute_local_to_global (VectorType &dst) const;
* function ConstraintMatrix::distribute_local_to_global. Note that if
* vectorization is enabled, the DoF values for several cells are used.
*
- * If the class was constructed through a MappingFEEvaluation object, only
- * one single cell is used by this class and this function extracts the
- * values of the underlying components on this cell. This call is slower
- * than the ones done through a MatrixFree object and lead to a structure
- * that does not effectively use vectorization in the evaluate routines
- * based on these values (instead, VectorizedArray<Number>::n_array_elements
- * same copies are worked on). In that case, no constraints can be
- * processed as these are not available here.
+ * If this class was constructed without a MatrixFree object and the
+ * information is acquired on the fly through a
+ * DoFHandler<dim>::cell_iterator, only one single cell is used by this
+ * class and this function extracts the values of the underlying components
+ * on the given cell. This call is slower than the ones done through a
+ * MatrixFree object and lead to a structure that does not effectively use
+ * vectorization in the evaluate routines based on these values (instead,
+ * VectorizedArray::n_array_elements same copies are worked on).
*/
template<typename VectorType>
void set_dof_values (VectorType &dst) const;
*/
VectorizedArray<Number> *begin_hessians ();
+ /**
+ * Returns the numbering of local degrees of freedom within the evaluation
+ * routines of FEEvaluation in terms of the standard numbering on finite
+ * elements.
+ */
+ const std::vector<unsigned int> &
+ get_internal_dof_numbering() const;
+
//@}
protected:
/**
* Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type MappingFEEvaluation
- * and a DoFHandler in order to allow for reading out the finite element
- * data. It uses the data provided by dof_handler.get_fe(). If the element
- * is vector-valued, the optional argument allows to specify the index of
- * the base element (as long as the element is primitive, non-primitive are
- * not supported currently).
+ * FEValues. The arguments are similar to the ones passed to the constructor
+ * of FEValues, with the notable difference that FEEvaluation expects a
+ * one-dimensional quadrature formula, Quadrature<1>, instead of a @p dim
+ * dimensional one. The finite element can be both scalar or vector valued,
+ * but this method always only selects a scalar base element at a time (with
+ * @p n_components copies as specified by the class template argument). For
+ * vector-valued elements, the optional argument @p first_selected_component
+ * allows to specify the index of the base element to be used for
+ * evaluation. Note that the internal data structures always assume that the
+ * base element is primitive, non-primitive are not supported currently.
*
- * With initialization from a FEValues object, no call to a reinit method of
- * this class is necessary. Instead, it is enough if the geometry is
- * initialized to a given cell iterator. It can also read from or write to
- * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
- * types (which is less efficient with MPI since index translation has to be
- * done), but of course only for one cell at a time. Hence, a kernel using
- * this method does not vectorize over several elements (which is most
- * efficient for vector operations), but only possibly within the element if
- * the evaluate/integrate routines are combined (e.g. for computing cell
- * matrices).
- */
- FEEvaluationBase (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component = 0);
-
- /**
- * Copy constructor
+ * As known from FEValues, a call to the reinit method with a
+ * Triangulation::cell_iterator is necessary to make the geometry and
+ * degrees of freedom of the current class known. If the iterator includes
+ * DoFHandler information (i.e., it is a DoFHandler::cell_iterator or
+ * similar), the initialization allows to also read from or write to vectors
+ * in the standard way for DoFHandler::active_cell_iterator types for
+ * one cell at a time. However, this approach is much slower than the path
+ * with MatrixFree with MPI since index translation has to be done. As only
+ * one cell at a time is used, this method does not vectorize over several
+ * elements (which is most efficient for vector operations), but only
+ * possibly within the element if the evaluate/integrate routines are
+ * combined inside user code (e.g. for computing cell matrices).
+ *
+ * The optional FEEvaluationBase object allows several FEEvaluation objects
+ * to share the geometry evaluation, i.e., the underlying mapping and
+ * quadrature points do only need to be evaluated once. This only works if
+ * the quadrature formulas are the same. Otherwise, a new evaluation object
+ * is created. Make sure to not pass an optional object around when you
+ * intend to use the FEEvaluation object in %parallel with another one
+ * because otherwise the intended sharing may create race conditions.
+ */
+ template <int n_components_other>
+ FEEvaluationBase (const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBase<dim,n_components_other,Number> *other);
+
+ /**
+ * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe,
+ * quadrature, and update flags, the underlying geometry evaluation based on
+ * FEValues will be deep-copied in order to allow for using in parallel with
+ * threads.
*/
FEEvaluationBase (const FEEvaluationBase &other);
bool gradients_quad_submitted;
/**
- * Geometry data generated by FEValues on the fly.
+ * Geometry data that can be generated FEValues on the fly with the
+ * respective constructor.
*/
- SmartPointer<const MappingFEEvaluation<dim,Number> > mapped_geometry;
+ std_cxx1x::shared_ptr<internal::MatrixFreeFunctions::MappingDataOnTheFly<dim,Number> > mapped_geometry;
/**
- * A pointer to the underlying DoFHandler.
+ * For use with on-the-fly evaluation, provide a data structure to store the
+ * global dof indices on the current cell from a reinit call.
*/
- const DoFHandler<dim> *dof_handler;
+ std::vector<types::global_dof_index> old_style_dof_indices;
/**
- * For a DoFHandler with more than one finite element, select at which
+ * For a FiniteElement with more than one finite element, select at which
* component this data structure should start.
*/
const unsigned int first_selected_component;
* MatrixFree object was given at initialization.
*/
mutable std::vector<types::global_dof_index> local_dof_indices;
+
+ /**
+ * Make other FEEvaluationBase as well as FEEvaluation objects friends.
+ */
+ template <int, int, typename> friend class FEEvaluationBase;
+ template <int, int, int, int, typename> friend class FEEvaluation;
};
const unsigned int n_q_points);
/**
- * Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type MappingFEEvaluation
- * and a DoFHandler in order to allow for reading out the finite element
- * data. It uses the data provided by dof_handler.get_fe(). If the element
- * is vector-valued, the optional argument allows to specify the index of
- * the base element (as long as the element is primitive, non-primitive are
- * not supported currently).
- *
- * With initialization from a FEValues object, no call to a reinit method of
- * this class is necessary. Instead, it is enough if the geometry is
- * initialized to a given cell iterator. It can also read from or write to
- * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
- * types (which is less efficient with MPI since index translation has to be
- * done), but of course only for one cell at a time. Hence, a kernel using
- * this method does not vectorize over several elements (which is most
- * efficient for vector operations), but only possibly within the element if
- * the evaluate/integrate routines are combined (e.g. for computing cell
- * matrices).
- * With this initialization, no call to a reinit method of this
- * class. Instead, it is enough if the geometry is initialized to a given
- * cell iterator. Moreover, beware that a kernel using this method does not
- * vectorize over several elements (which is most efficient for vector
- * operations), but only possibly within the element if the
- * evaluate/integrate routines are combined (e.g. for matrix assembly).
- */
- FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component = 0);
+ * Constructor with reduced functionality for similar usage of FEEvaluation
+ * as FEValues, including matrix assembly.
+ */
+ template <int n_components_other>
+ FEEvaluationAccess (const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBase<dim,n_components_other,Number> *other);
/**
* Copy constructor
const unsigned int n_q_points);
/**
- * Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type MappingFEEvaluation
- * and a DoFHandler in order to allow for reading out the finite element
- * data. It uses the data provided by dof_handler.get_fe(). If the element
- * is vector-valued, the optional argument allows to specify the index of
- * the base element (as long as the element is primitive, non-primitive are
- * not supported currently).
- *
- * With initialization from a FEValues object, no call to a reinit method of
- * this class is necessary. Instead, it is enough if the geometry is
- * initialized to a given cell iterator. It can also read from or write to
- * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
- * types (which is less efficient with MPI since index translation has to be
- * done), but of course only for one cell at a time. Hence, a kernel using
- * this method does not vectorize over several elements (which is most
- * efficient for vector operations), but only possibly within the element if
- * the evaluate/integrate routines are combined (e.g. for computing cell
- * matrices).
- * With this initialization, no call to a reinit method of this
- * class. Instead, it is enough if the geometry is initialized to a given
- * cell iterator. Moreover, beware that a kernel using this method does not
- * vectorize over several elements (which is most efficient for vector
- * operations), but only possibly within the element if the
- * evaluate/integrate routines are combined (e.g. for matrix assembly).
- */
- FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component = 0);
+ * Constructor with reduced functionality for similar usage of FEEvaluation
+ * as FEValues, including matrix assembly.
+ */
+ template <int n_components_other>
+ FEEvaluationAccess (const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBase<dim,n_components_other,Number> *other);
/**
* Copy constructor
const unsigned int n_q_points);
/**
- * Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type MappingFEEvaluation
- * and a DoFHandler in order to allow for reading out the finite element
- * data. It uses the data provided by dof_handler.get_fe(). If the element
- * is vector-valued, the optional argument allows to specify the index of
- * the base element (as long as the element is primitive, non-primitive are
- * not supported currently).
- *
- * With initialization from a FEValues object, no call to a reinit method of
- * this class is necessary. Instead, it is enough if the geometry is
- * initialized to a given cell iterator. It can also read from or write to
- * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
- * types (which is less efficient with MPI since index translation has to be
- * done), but of course only for one cell at a time. Hence, a kernel using
- * this method does not vectorize over several elements (which is most
- * efficient for vector operations), but only possibly within the element if
- * the evaluate/integrate routines are combined (e.g. for computing cell
- * matrices).
- * With this initialization, no call to a reinit method of this
- * class. Instead, it is enough if the geometry is initialized to a given
- * cell iterator. Moreover, beware that a kernel using this method does not
- * vectorize over several elements (which is most efficient for vector
- * operations), but only possibly within the element if the
- * evaluate/integrate routines are combined (e.g. for matrix assembly).
- */
- FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component = 0);
+ * Constructor with reduced functionality for similar usage of FEEvaluation
+ * as FEValues, including matrix assembly.
+ */
+ template <int n_components_other>
+ FEEvaluationAccess (const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBase<dim,n_components_other,Number> *other);
/**
* Copy constructor
const unsigned int n_q_points);
/**
- * Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type MappingFEEvaluation
- * and a DoFHandler in order to allow for reading out the finite element
- * data. It uses the data provided by dof_handler.get_fe(). If the element
- * is vector-valued, the optional argument allows to specify the index of
- * the base element (as long as the element is primitive, non-primitive are
- * not supported currently).
- *
- * With initialization from a FEValues object, no call to a reinit method of
- * this class is necessary. Instead, it is enough if the geometry is
- * initialized to a given cell iterator. It can also read from or write to
- * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
- * types (which is less efficient with MPI since index translation has to be
- * done), but of course only for one cell at a time. Hence, a kernel using
- * this method does not vectorize over several elements (which is most
- * efficient for vector operations), but only possibly within the element if
- * the evaluate/integrate routines are combined (e.g. for computing cell
- * matrices).
- * With this initialization, no call to a reinit method of this
- * class. Instead, it is enough if the geometry is initialized to a given
- * cell iterator. Moreover, beware that a kernel using this method does not
- * vectorize over several elements (which is most efficient for vector
- * operations), but only possibly within the element if the
- * evaluate/integrate routines are combined (e.g. for matrix assembly).
- */
- FEEvaluationAccess (const MappingFEEvaluation<1,Number> &geometry,
- const DoFHandler<1> &dof_handler,
- const unsigned int first_selected_component = 0);
+ * Constructor with reduced functionality for similar usage of FEEvaluation
+ * as FEValues, including matrix assembly.
+ */
+ template <int n_components_other>
+ FEEvaluationAccess (const Mapping<1> &mapping,
+ const FiniteElement<1> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBase<1,n_components_other,Number> *other);
/**
* Copy constructor
* functions that make it much faster (between 5 and 500, depending on the
* polynomial order).
*
- * This class can be used in two different ways. The first way is to
- * initialize it from a MatrixFree object that caches everything related to
- * the degrees of freedom and the mapping information. This way, it is
- * possible to use vectorization for applying a vector operation for several
- * cells at once. The second form of usage is to initialize it from geometry
- * information generated by FEValues, which is stored in the class
- * MappingFEEvaluation. Here, the operations can only work on a single cell, but
- * possibly be vectorized by combining several operations (e.g. when
- * performing matrix assembly).
+ * <h3>Usage and initialization</h3>
+ *
+ * <h4>Fast usage in combination with MatrixFree</h4>
*
- * This class contains specialized evaluation routines for several elements,
- * including standard FE_Q or FE_DGQ elements and quadrature points symmetric
- * around 0.5 (like Gauss quadrature), FE_DGP elements based on truncated
- * tensor products as well as the faster case of Gauss-Lobatto elements with
- * Gauss-Lobatto quadrature which give diagonal mass matrices and quicker
- * evaluation internally. Note that many of the operations available through
- * this class are inherited from the base class FEEvaluationBase, in
- * particular reading from and writing to vectors. Also, the class inherits
- * from FEEvaluationAccess that implements access to values, gradients and
- * Hessians of the finite element function on quadrature points.
+ * The first and foremost way of usage is to initialize this class from a
+ * MatrixFree object that caches everything related to the degrees of freedom
+ * and the mapping information. This way, it is possible to use vectorization
+ * for applying a vector operation for several cells at once. This setting is
+ * explained in the step-37 and step-48 tutorial programs. For vector-valued
+ * problems, the deal.II test suite includes a few additional examples as
+ * well, e.g. the Stokes operator found at
+ * https://github.com/dealii/dealii/blob/master/tests/matrix_free/matrix_vector_stokes_noflux.cc
+ *
+ * For most operator evaluation tasks, this path provides the most efficient
+ * solution by combining pre-computed data for the mapping (Jacobian
+ * transformations for the geometry description) with on-the-fly evaluation of
+ * basis functions. In other words, the framework provides a trade-off between
+ * memory usage and initialization of objects that is suitable for matrix-free
+ * operator evaluation.
+ *
+ * <h4>Usage without pre-initialized MatrixFree object</h4>
+ *
+ * The second form of usage is to initialize FEEvaluation from geometry
+ * information generated by FEValues. This allows to apply the integration
+ * loops on the fly without prior initialization of MatrixFree objects. This
+ * can be useful when the memory and initialization cost of MatrixFree is not
+ * acceptable, e.g. when a different number of quadrature points should be
+ * used for one single evaluation in error computation. Also, when using the
+ * routines of this class to assemble matrices the trade-off implied by the
+ * MatrixFree class may not be desired. In such a case, the cost to initialize
+ * the necessary geometry data on the fly is comparably low and thus avoiding
+ * a global object MatrixFree can be useful. When used in this way, reinit
+ * methods reminiscent from FEValues with a cell iterator are to be
+ * used. However, note that this model results in working on a single cell at
+ * a time, with geometry data duplicated in all components of the vectorized
+ * array. Thus, vectorization is only useful when it can apply the same
+ * operation on different data, e.g. when performing matrix assembly.
+ *
+ * As an example, consider the following code to assemble the contributions to
+ * the Laplace matrix:
+ *
+ * @code
+ * FEEvaluation<dim,fe_degree> fe_eval (mapping, finite_element,
+ * QGauss<1>(fe_degree+1), flags);
+ * for (cell = dof_handler.begin_active();
+ * cell != dof_handler.end();
+ * ++cell)
+ * {
+ * fe_eval.reinit(cell);
+ * for (unsigned int i=0; i<dofs_per_cell;
+ * i += VectorizedArray<double>::n_array_elements)
+ * {
+ * const unsigned int n_items =
+ * i+VectorizedArray<double>::n_array_elements > dofs_per_cell ?
+ * (dofs_per_cell - i) : VectorizedArray<double>::n_array_elements;
+ *
+ * // Set n_items unit vectors
+ * for (unsigned int j=0; j<dofs_per_cell; ++j)
+ * fe_eval.begin_dof_values()[j] = VectorizedArray<double>();
+ * for (unsigned int v=0; v<n_items; ++v)
+ * fe_eval.begin_dof_values()[i+v][v] = 1.;
+ *
+ * // Apply operator on unit vector
+ * fe_eval.evaluate(true, true);
+ * for (unsigned int q=0; q<n_q_points; ++q)
+ * {
+ * fe_eval.submit_value(10.*fe_eval.get_value(q), q);
+ * fe_eval.submit_gradient(fe_eval.get_gradient(q), q);
+ * }
+ * fe_eval.integrate(true, true);
+ *
+ * // Insert computed entries in matrix
+ * for (unsigned int v=0; v<n_items; ++v)
+ * for (unsigned int j=0; j<dofs_per_cell; ++j)
+ * cell_matrix(fe_eval.get_internal_dof_numbering()[j],
+ * fe_eval.get_internal_dof_numbering()[i+v])
+ * = fe_eval.begin_dof_values()[j][v];
+ * }
+ * ...
+ * }
+ * @endcode
+ *
+ * This code generates the columns of the cell matrix with the loop over @p i
+ * above. The way this is done is the following: FEEvaluation's routines focus
+ * on the evaluation of finite element operators, so the way to get a cell
+ * matrix out of an operator evaluation is to apply it to all the unit vectors
+ * on the cell. This might seem inefficient but the evaluation routines used
+ * here are so quick that they still work much faster than what is possible
+ * with FEValues.
+ *
+ * Due to vectorization, we can actually generate matrix data for several unit
+ * vectors at a time (e.g. 4). The variable @p n_items make sure that we do
+ * the last iteration where the number of cell dofs is not divisible by the
+ * vectorization length correctly. Also note that we need to get the internal
+ * dof numbering applied by fe_eval because FEEvaluation internally uses a
+ * lexicographic numbering of degrees of freedom. This is necessary to
+ * efficiently work with tensor products where all degrees of freedom along a
+ * dimension must be laid out in a regular way.
+ *
+ * <h3>Description of evaluation routines</h3>
+ *
+ * This class contains specialized evaluation routines for several elements
+ * based on tensor-product quadrature formulas and tensor-product-like shape
+ * functions, including standard FE_Q or FE_DGQ elements and quadrature points
+ * symmetric around 0.5 (like Gauss quadrature), FE_DGP elements based on
+ * truncated tensor products as well as the faster case of Gauss-Lobatto
+ * elements with Gauss-Lobatto quadrature which give diagonal mass matrices
+ * and quicker evaluation internally. The main benefit of this class is the
+ * evaluation of all shape functions in all quadrature or integration over all
+ * shape functions in <code>dim (fe_degree+1)<sup>dim+1</sup> </code>
+ * operations instead of the slower <code>
+ * (fe_degree+1)<sup>2*dim</sup></code> complexity in the evaluation routines
+ * of FEValues.
+ *
+ * Note that many of the operations available through this class are inherited
+ * from the base class FEEvaluationBase, in particular reading from and
+ * writing to vectors. Also, the class inherits from FEEvaluationAccess that
+ * implements access to values, gradients and Hessians of the finite element
+ * function on quadrature points.
*
* This class assumes that shape functions of the FiniteElement under
* consideration do <em>not</em> depend on the actual shape of the cells in
* fe_degree+1 degrees of freedom per coordinate direction
*
* @param n_q_points_1d Number of points in the quadrature formula in 1D,
- * usually chosen as fe_degree+1
+ * defaults to fe_degree+1
*
* @param n_components Number of vector components when solving a system of
* PDEs. If the same operation is applied to several
* components of a PDE (e.g. a vector Laplace equation), they
* can be applied simultaneously with one call (and often
- * more efficiently)
+ * more efficiently). Defaults to 1.
*
- * @param Number Number format, usually @p double or @p float
+ * @param Number Number format, usually @p double or @p float.
+ * Defaults to @p double
*
* @author Katharina Kormann and Martin Kronbichler, 2010, 2011
*/
-template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
- int n_components_ = 1, typename Number = double >
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number >
class FEEvaluation : public FEEvaluationAccess<dim,n_components_,Number>
{
public:
const unsigned int fe_no = 0,
const unsigned int quad_no = 0);
- /**
- * Copy constructor
- */
- FEEvaluation (const FEEvaluation &other);
-
/**
* Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type MappingFEEvaluation
- * and a DoFHandler in order to allow for reading out the finite element
- * data. It uses the data provided by dof_handler.get_fe(). If the element
- * is vector-valued, the optional argument allows to specify the index of
- * the base element (as long as the element is primitive, non-primitive are
- * not supported currently).
+ * FEValues. The arguments are similar to the ones passed to the constructor
+ * of FEValues, with the notable difference that FEEvaluation expects a
+ * one-dimensional quadrature formula, Quadrature<1>, instead of a @p dim
+ * dimensional one. The finite element can be both scalar or vector valued,
+ * but this method always only selects a scalar base element at a time (with
+ * @p n_components copies as specified by the class template). For
+ * vector-valued elements, the optional argument @p first_selected_component
+ * allows to specify the index of the base element to be used for
+ * evaluation. Note that the internal data structures always assume that the
+ * base element is primitive, non-primitive are not supported currently.
*
- * With initialization from a FEValues object, no call to a reinit method of
- * this class is necessary. Instead, it is enough if the geometry is
- * initialized to a given cell iterator. It can also read from or write to
- * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
- * types (which is less efficient with MPI since index translation has to be
- * done), but of course only for one cell at a time. Hence, a kernel using
- * this method does not vectorize over several elements (which is most
- * efficient for vector operations), but only possibly within the element if
- * the evaluate/integrate routines are combined (e.g. for computing cell
- * matrices).
- */
- FEEvaluation (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component = 0);
+ * As known from FEValues, a call to the reinit method with a
+ * Triangulation<dim>::cell_iterator is necessary to make the geometry and
+ * degrees of freedom of the current class known. If the iterator includes
+ * DoFHandler information (i.e., it is a DoFHandler<dim>::cell_iterator or
+ * similar), the initialization allows to also read from or write to vectors
+ * in the standard way for DoFHandler<dim>::active_cell_iterator types for
+ * one cell at a time. However, this approach is much slower than the path
+ * with MatrixFree with MPI since index translation has to be done. As only
+ * one cell at a time is used, this method does not vectorize over several
+ * elements (which is most efficient for vector operations), but only
+ * possibly within the element if the evaluate/integrate routines are
+ * combined inside user code (e.g. for computing cell matrices).
+ */
+ FEEvaluation (const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component = 0);
+
+ /**
+ * Constructor for the reduced functionality. Similar as the other
+ * constructor but uses MappingQ1 implicitly.
+ */
+ FEEvaluation (const FiniteElement<dim> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component = 0);
+
+ /**
+ * Constructor for the reduced functionality. Similar to the other
+ * constructor with FiniteElement argument but using another
+ * FEEvaluationBase object to provide info about the geometry. This allows
+ * several FEEvaluation objects to share the geometry evaluation, i.e., the
+ * underlying mapping and quadrature points do only need to be evaluated
+ * once. Make sure to not pass an optional object around when you intend to
+ * use the FEEvaluation object in %parallel to the given one because
+ * otherwise the intended sharing may create race conditions.
+ */
+ template <int n_components_other>
+ FEEvaluation (const FiniteElement<dim> &fe,
+ const FEEvaluationBase<dim,n_components_other,Number> &other,
+ const unsigned int first_selected_component = 0);
+
+ /**
+ * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe,
+ * quadrature, and update flags, the underlying geometry evaluation based on
+ * FEValues will be deep-copied in order to allow for using in parallel with
+ * threads.
+ */
+ FEEvaluation (const FEEvaluation &other);
/**
* Evaluates the function values, the gradients, and the Laplacians of the
:
BaseClass (matrix_free, fe_no, quad_no)
{}
-
- /**
- * Constructor.
- */
- FEEvaluationGeneral (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component = 0) DEAL_II_DEPRECATED
-:
- BaseClass (geometry, dof_handler, first_selected_component)
- {}
};
:
BaseClass (matrix_free, fe_no, quad_no)
{}
-
- /**
- * Constructor.
- */
- FEEvaluationGL (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component = 0) DEAL_II_DEPRECATED
-:
- BaseClass (geometry, dof_handler, first_selected_component)
- {}
};
:
BaseClass (matrix_free, fe_no, quad_no)
{}
-
- /**
- * Constructor.
- */
- FEEvaluationDGP (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component = 0) DEAL_II_DEPRECATED
-:
- BaseClass (geometry, dof_handler, first_selected_component)
- {}
};
cell (numbers::invalid_unsigned_int),
cell_type (internal::MatrixFreeFunctions::undefined),
cell_data_number (0),
- dof_handler (0),
first_selected_component (0)
{
for (unsigned int c=0; c<n_components_; ++c)
template <int dim, int n_components_, typename Number>
+template <int n_components_other>
inline
FEEvaluationBase<dim,n_components_,Number>
-::FEEvaluationBase (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler_in,
- const unsigned int first_selected_component)
+::FEEvaluationBase (const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBase<dim,n_components_other,Number> *other)
:
quad_no (-1),
n_fe_components (n_components_),
matrix_info (0),
dof_info (0),
mapping_info (0),
- stored_shape_info (new internal::MatrixFreeFunctions::ShapeInfo<Number>(geometry.get_quadrature(), dof_handler_in.get_fe(), dof_handler_in.get_fe().component_to_base_index(first_selected_component).first)),
+ // select the correct base element from the given FE component
+ stored_shape_info (new internal::MatrixFreeFunctions::ShapeInfo<Number>(quadrature, fe, fe.component_to_base_index(first_selected_component).first)),
data (stored_shape_info.get()),
cartesian_data (0),
- jacobian (geometry.get_inverse_jacobians().begin()),
- J_value (geometry.get_JxW_values().begin()),
+ jacobian (0),
+ J_value (0),
quadrature_weights (0),
- quadrature_points (geometry.get_quadrature_points().begin()),
+ quadrature_points (0),
jacobian_grad (0),
jacobian_grad_upper(0),
cell (0),
cell_type (internal::MatrixFreeFunctions::general),
cell_data_number (0),
- mapped_geometry (&geometry),
- dof_handler (&dof_handler_in),
- first_selected_component (first_selected_component)
+ // keep the number of the selected component within the current base element
+ // for reading dof values
+ first_selected_component (fe.component_to_base_index(first_selected_component).second)
{
const unsigned int base_element_number =
- dof_handler_in.get_fe().component_to_base_index(first_selected_component).first;
+ fe.component_to_base_index(first_selected_component).first;
for (unsigned int c=0; c<n_components_; ++c)
{
values_dofs[c] = 0;
for (unsigned int d=0; d<(dim*dim+dim)/2; ++d)
hessians_quad[c][d] = 0;
}
- Assert(dof_handler->get_fe().element_multiplicity(base_element_number) == 1 ||
- dof_handler->get_fe().element_multiplicity(base_element_number)-first_selected_component >= n_components_,
+
+ Assert(other == 0 || other->mapped_geometry.get() != 0, ExcInternalError());
+ if (other != 0 &&
+ other->mapped_geometry->get_quadrature() == quadrature)
+ mapped_geometry = other->mapped_geometry;
+ else
+ mapped_geometry.reset(new internal::MatrixFreeFunctions::
+ MappingDataOnTheFly<dim,Number>(mapping, quadrature,
+ update_flags));
+ jacobian = mapped_geometry->get_inverse_jacobians().begin();
+ J_value = mapped_geometry->get_JxW_values().begin();
+ quadrature_points = mapped_geometry->get_quadrature_points().begin();
+
+ Assert(fe.element_multiplicity(base_element_number) == 1 ||
+ fe.element_multiplicity(base_element_number)-first_selected_component >= n_components_,
ExcMessage("The underlying element must at least contain as many "
"components as requested by this class"));
}
cell (other.cell),
cell_type (other.cell_type),
cell_data_number (other.cell_data_number),
- mapped_geometry (other.mapped_geometry),
- dof_handler (other.dof_handler),
first_selected_component (other.first_selected_component)
{
for (unsigned int c=0; c<n_components_; ++c)
for (unsigned int d=0; d<(dim*dim+dim)/2; ++d)
hessians_quad[c][d] = 0;
}
+
+ // Create deep copy of mapped geometry for use in parallel...
+ if (other.mapped_geometry.get() != 0)
+ {
+ mapped_geometry.reset
+ (new internal::MatrixFreeFunctions::
+ MappingDataOnTheFly<dim,Number>(other.mapped_geometry->get_fe_values().get_mapping(),
+ other.mapped_geometry->get_quadrature(),
+ other.mapped_geometry->get_fe_values().get_update_flags()));
+ jacobian = mapped_geometry->get_inverse_jacobians().begin();
+ J_value = mapped_geometry->get_JxW_values().begin();
+ quadrature_points = mapped_geometry->get_quadrature_points().begin();
+ }
}
void
FEEvaluationBase<dim,n_components_,Number>::reinit (const unsigned int cell_in)
{
- Assert (mapped_geometry == 0, ExcMessage("FEEvaluation was initialized without a matrix-free object. Integer indexing is not possible"));
+ Assert (mapped_geometry == 0,
+ ExcMessage("FEEvaluation was initialized without a matrix-free object."
+ " Integer indexing is not possible"));
if (mapped_geometry != 0)
return;
Assert (dof_info != 0, ExcNotInitialized());
+template <int dim, int n_components_, typename Number>
+template <typename DH, bool level_dof_access>
+inline
+void
+FEEvaluationBase<dim,n_components_,Number>
+::reinit (const TriaIterator<DoFCellAccessor<DH,level_dof_access> > &cell)
+{
+ Assert(matrix_info == 0,
+ ExcMessage("Cannot use initialization from cell iterator if "
+ "initialized from MatrixFree object. Use variant for "
+ "on the fly computation with arguments as for FEValues "
+ "instead"));
+ Assert(mapped_geometry.get() != 0, ExcNotInitialized());
+ mapped_geometry->reinit(static_cast<typename Triangulation<dim>::cell_iterator>(cell));
+ local_dof_indices.resize(cell->get_fe().dofs_per_cell);
+ if (level_dof_access)
+ cell->get_mg_dof_indices(local_dof_indices);
+ else
+ cell->get_dof_indices(local_dof_indices);
+}
+
+
+
+template <int dim, int n_components_, typename Number>
+inline
+void
+FEEvaluationBase<dim,n_components_,Number>
+::reinit (const typename Triangulation<dim>::cell_iterator &cell)
+{
+ Assert(matrix_info == 0,
+ ExcMessage("Cannot use initialization from cell iterator if "
+ "initialized from MatrixFree object. Use variant for "
+ "on the fly computation with arguments as for FEValues "
+ "instead"));
+ Assert(mapped_geometry.get() != 0, ExcNotInitialized());
+ mapped_geometry->reinit(cell);
+}
+
+
+
template <int dim, int n_components_, typename Number>
inline
unsigned int
// process constraints and need not care about vectorization
if (matrix_info == 0)
{
- Assert (dof_handler != 0, ExcNotInitialized());
- typename DoFHandler<dim>::cell_iterator cell (&dof_handler->get_tria(),
- mapped_geometry->get_cell()->level(),
- mapped_geometry->get_cell()->index(),
- dof_handler);
- local_dof_indices.resize(dof_handler->get_fe().dofs_per_cell);
- cell->get_dof_indices(local_dof_indices);
+ Assert (!local_dof_indices.empty(), ExcNotInitialized());
unsigned int index = first_selected_component * this->data->dofs_per_cell;
for (unsigned int comp = 0; comp<n_components; ++comp)
/*------------------------------ access to data fields ----------------------*/
+template <int dim, int n_components, typename Number>
+inline
+const std::vector<unsigned int> &
+FEEvaluationBase<dim,n_components,Number>::
+get_internal_dof_numbering() const
+{
+ return data->lexicographic_numbering;
+}
+
+
+
template <int dim, int n_components, typename Number>
inline
template <int dim, int n_components_, typename Number>
+template <int n_components_other>
inline
FEEvaluationAccess<dim,n_components_,Number>
-::FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component)
+::FEEvaluationAccess (const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBase<dim,n_components_other,Number> *other)
:
- FEEvaluationBase <dim,n_components_,Number> (geometry, dof_handler, first_selected_component)
+ FEEvaluationBase <dim,n_components_,Number>(mapping, fe, quadrature, update_flags,
+ first_selected_component, other)
{}
template <int dim, typename Number>
+template <int n_components_other>
inline
FEEvaluationAccess<dim,1,Number>
-::FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component)
+::FEEvaluationAccess (const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBase<dim,n_components_other,Number> *other)
:
- FEEvaluationBase <dim,1,Number> (geometry, dof_handler, first_selected_component)
+ FEEvaluationBase <dim,1,Number> (mapping, fe, quadrature, update_flags,
+ first_selected_component, other)
{}
template <int dim, typename Number>
+template <int n_components_other>
inline
FEEvaluationAccess<dim,dim,Number>
-::FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component)
+::FEEvaluationAccess (const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBase<dim,n_components_other,Number> *other)
:
- FEEvaluationBase <dim,dim,Number> (geometry, dof_handler, first_selected_component)
+ FEEvaluationBase <dim,dim,Number> (mapping, fe, quadrature, update_flags,
+ first_selected_component, other)
{}
template <typename Number>
+template <int n_components_other>
inline
FEEvaluationAccess<1,1,Number>
-::FEEvaluationAccess (const MappingFEEvaluation<1,Number> &geometry,
- const DoFHandler<1> &dof_handler,
- const unsigned int first_selected_component)
+::FEEvaluationAccess (const Mapping<1> &mapping,
+ const FiniteElement<1> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBase<1,n_components_other,Number> *other)
:
- FEEvaluationBase <1,1,Number> (geometry, dof_handler, first_selected_component)
+ FEEvaluationBase <1,1,Number> (mapping, fe, quadrature, update_flags,
+ first_selected_component, other)
{}
typename Number>
inline
FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluation (const MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component)
+::FEEvaluation (const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component)
:
- BaseClass (geometry, dof_handler, first_selected_component),
+ BaseClass (mapping, fe, quadrature, update_flags,
+ first_selected_component,
+ static_cast<FEEvaluationBase<dim,1,Number>*>(0)),
+ dofs_per_cell (this->data->dofs_per_cell)
+{
+ check_template_arguments(numbers::invalid_unsigned_int);
+ set_data_pointers();
+}
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+inline
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluation (const FiniteElement<dim> &fe,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component)
+ :
+ BaseClass (StaticMappingQ1<dim>::mapping, fe, quadrature, update_flags,
+ first_selected_component,
+ static_cast<FEEvaluationBase<dim,1,Number>*>(0)),
+ dofs_per_cell (this->data->dofs_per_cell)
+{
+ check_template_arguments(numbers::invalid_unsigned_int);
+ set_data_pointers();
+}
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+template <int n_components_other>
+inline
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluation (const FiniteElement<dim> &fe,
+ const FEEvaluationBase<dim,n_components_other,Number> &other,
+ const unsigned int first_selected_component)
+ :
+ BaseClass (other.mapped_geometry->get_fe_values().get_mapping(),
+ fe, other.mapped_geometry->get_quadrature(),
+ other.mapped_geometry->get_fe_values().get_update_flags(),
+ first_selected_component, &other),
dofs_per_cell (this->data->dofs_per_cell)
{
check_template_arguments(numbers::invalid_unsigned_int);
message += Utilities::int_to_string(fe_degree) + ",";
message += Utilities::int_to_string(n_q_points_1d);
message += "," + Utilities::int_to_string(n_components);
- message += ",Number>(data, ";
- message += Utilities::int_to_string(fe_no) + ", ";
- message += Utilities::int_to_string(this->quad_no) + ")\n";
+ message += ",Number>(data";
+ if (fe_no != numbers::invalid_unsigned_int)
+ {
+ message += ", " + Utilities::int_to_string(fe_no) + ", ";
+ message += Utilities::int_to_string(this->quad_no);
+ }
+ message += ")\n";
// check whether some other vector component has the correct number of
// points
message += Utilities::int_to_string(fe_degree) + ",";
message += Utilities::int_to_string(n_q_points_1d);
message += "," + Utilities::int_to_string(n_components);
- message += ",Number>(data, ";
- message += Utilities::int_to_string(proposed_dof_comp) + ", ";
- message += Utilities::int_to_string(proposed_quad_comp) + ")?\n";
+ message += ",Number>(data";
+ if (fe_no != numbers::invalid_unsigned_int)
+ {
+ message += ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
+ message += Utilities::int_to_string(proposed_quad_comp);
+ }
+ message += ")?\n";
std::string correct_pos;
if (proposed_dof_comp != fe_no)
correct_pos = " ^ ";
message += Utilities::int_to_string(this->data->fe_degree) + ",";
message += Utilities::int_to_string(proposed_n_q_points_1d);
message += "," + Utilities::int_to_string(n_components);
- message += ",Number>(data, ";
- message += Utilities::int_to_string(fe_no) + ", ";
- message += Utilities::int_to_string(this->quad_no) + ")?\n";
+ message += ",Number>(data";
+ if (fe_no != numbers::invalid_unsigned_int)
+ {
+ message += ", " + Utilities::int_to_string(fe_no) + ", ";
+ message += Utilities::int_to_string(this->quad_no);
+ }
+ message += ")?\n";
std::string correct_pos;
if (this->data->fe_degree != fe_degree)
correct_pos = " ^";
{
Assert (this->dof_values_initialized == true,
internal::ExcAccessToUninitializedField());
+ Assert(this->matrix_info != 0 ||
+ this->mapped_geometry->is_initialized(), ExcNotInitialized());
// Select algorithm matching the element type at run time (the function
// pointer is easy to predict, so negligible in cost)
if (integrate_grad == true)
Assert (this->gradients_quad_submitted == true,
internal::ExcAccessToUninitializedField());
+ Assert(this->matrix_info != 0 ||
+ this->mapped_geometry->is_initialized(), ExcNotInitialized());
// Select algorithm matching the element type at run time (the function
// pointer is easy to predict, so negligible in cost)
projects / dealii.git / commitdiff