#include <deal.II/base/config.h>
+#include <deal.II/base/table.h>
+#include <deal.II/base/derivative_form.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+
+#include <vector>
+
DEAL_II_NAMESPACE_OPEN
+template <int,int> class FiniteElement;
+
+
/*!@addtogroup feaccess */
/*@{*/
}
+namespace internal
+{
+ namespace FEValues
+ {
+ /**
+ * A class that stores all of the mapping related data used in
+ * dealii::FEValues, dealii::FEFaceValues, and dealii::FESubfaceValues
+ * objects. Objects of this kind will be given
+ * as <i>output</i> argument when dealii::FEValues::reinit()
+ * calls Mapping::fill_fe_values() for a given cell, face, or subface.
+ *
+ * The data herein will then be provided as <i>input</i> argument in
+ * the following call to FiniteElement::fill_fe_values().
+ *
+ * @ingroup feaccess
+ */
+ template <int dim, int spacedim=dim>
+ class MappingRelatedData
+ {
+ public:
+ /**
+ * Initialize all vectors to correct size.
+ */
+ void initialize (const unsigned int n_quadrature_points,
+ const UpdateFlags flags);
+
+ /**
+ * Store an array of weights times the Jacobi determinant at the quadrature
+ * points. This function is reset each time reinit() is called. The Jacobi
+ * determinant is actually the reciprocal value of the Jacobi matrices
+ * stored in this class, see the general documentation of this class for
+ * more information.
+ *
+ * However, if this object refers to an FEFaceValues or FESubfaceValues
+ * object, then the JxW_values correspond to the Jacobian of the
+ * transformation of the face, not the cell, i.e. the dimensionality is that
+ * of a surface measure, not of a volume measure. In this case, it is
+ * computed from the boundary forms, rather than the Jacobian matrix.
+ */
+ std::vector<double> JxW_values;
+
+ /**
+ * Array of the Jacobian matrices at the quadrature points.
+ */
+ std::vector< DerivativeForm<1,dim,spacedim> > jacobians;
+
+ /**
+ * Array of the derivatives of the Jacobian matrices at the quadrature
+ * points.
+ */
+ std::vector<DerivativeForm<2,dim,spacedim> > jacobian_grads;
+
+ /**
+ * Array of the inverse Jacobian matrices at the quadrature points.
+ */
+ std::vector<DerivativeForm<1,spacedim,dim> > inverse_jacobians;
+
+ /**
+ * Array of quadrature points. This array is set up upon calling reinit()
+ * and contains the quadrature points on the real element, rather than on
+ * the reference element.
+ */
+ std::vector<Point<spacedim> > quadrature_points;
+
+ /**
+ * List of outward normal vectors at the quadrature points. This field is
+ * filled in by the finite element class.
+ */
+ std::vector<Point<spacedim> > normal_vectors;
+
+ /**
+ * List of boundary forms at the quadrature points. This field is filled in
+ * by the finite element class.
+ */
+ std::vector<Tensor<1,spacedim> > boundary_forms;
+ };
+
+
+ /**
+ * A class that stores all of the shape function related data used in
+ * dealii::FEValues, dealii::FEFaceValues, and dealii::FESubfaceValues
+ * objects. Objects of this kind will be given
+ * as <i>output</i> argument when dealii::FEValues::reinit()
+ * calls FiniteElement::fill_fe_values().
+ *
+ * @ingroup feaccess
+ */
+ template <int dim, int spacedim=dim>
+ class FiniteElementRelatedData
+ {
+ public:
+ /**
+ * Initialize all vectors to correct size.
+ */
+ void initialize (const unsigned int n_quadrature_points,
+ const FiniteElement<dim,spacedim> &fe,
+ const UpdateFlags flags);
+
+ /**
+ * Storage type for shape values. Each row in the matrix denotes the values
+ * of a single shape function at the different points, columns are for a
+ * single point with the different shape functions.
+ *
+ * If a shape function has more than one non-zero component (in deal.II
+ * diction: it is non-primitive), then we allocate one row per non-zero
+ * component, and shift subsequent rows backward. Lookup of the correct row
+ * for a shape function is thus simple in case the entire finite element is
+ * primitive (i.e. all shape functions are primitive), since then the shape
+ * function number equals the row number. Otherwise, use the
+ * #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 FiniteElement::get_nonzero_components()
+ * function which tells us which components are non-zero and thus have a row
+ * in the array presently under discussion.
+ */
+ typedef dealii::Table<2,double> ShapeVector;
+
+ /**
+ * Storage type for gradients. The layout of data is the same as for the
+ * #ShapeVector data type.
+ */
+ typedef std::vector<std::vector<Tensor<1,spacedim> > > GradientVector;
+
+ /**
+ * Likewise for second order derivatives.
+ */
+ typedef std::vector<std::vector<Tensor<2,spacedim> > > HessianVector;
+
+ /**
+ * Store the values of the shape functions at the quadrature points. See the
+ * description of the data type for the layout of the data in this field.
+ */
+ ShapeVector shape_values;
+
+ /**
+ * Store the gradients of the shape functions at the quadrature points. See
+ * the description of the data type for the layout of the data in this
+ * field.
+ */
+ GradientVector shape_gradients;
+
+ /**
+ * Store the 2nd derivatives of the shape functions at the quadrature
+ * points. See the description of the data type for the layout of the data
+ * in this field.
+ */
+ HessianVector shape_hessians;
+
+ /**
+ * When asked for the value (or gradient, or Hessian) of shape function i's
+ * c-th vector component, we need to look it up in the #shape_values,
+ * #shape_gradients and #shape_hessians arrays. The question is where in
+ * this array does the data for shape function i, component c reside. This
+ * is what this table answers.
+ *
+ * The format of the table is as follows: - It has dofs_per_cell times
+ * n_components entries. - The entry that corresponds to shape function i,
+ * component c is <code>i * n_components + c</code>. - The value stored at
+ * this position indicates the row in #shape_values and the other tables
+ * where the corresponding datum is stored for all the quadrature points.
+ *
+ * In the general, vector-valued context, the number of components is larger
+ * than one, but for a given shape function, not all vector components may
+ * be nonzero (e.g., if a shape function is primitive, then exactly one
+ * vector component is non-zero, while the others are all zero). For such
+ * zero components, #shape_values and friends do not have a row.
+ * Consequently, for vector components for which shape function i is zero,
+ * the entry in the current table is numbers::invalid_unsigned_int.
+ *
+ * On the other hand, the table is guaranteed to have at least one valid
+ * index for each shape function. In particular, for a primitive finite
+ * element, each shape function has exactly one nonzero component and so for
+ * each i, there is exactly one valid index within the range
+ * <code>[i*n_components, (i+1)*n_components)</code>.
+ */
+ std::vector<unsigned int> shape_function_to_row_table;
+ };
+ }
+}
+
+
/*@}*/
//TODO: Add access to mapping values to FEValuesBase
-//TODO: Several FEValuesBase of a system should share Mapping
/**
* A class that contains all data vectors for FEValues, FEFaceValues, and
* FESubfaceValues.
*
- * This class has been extracted from FEValuesBase to encapsulate in one
- * place all of the data, independent of the functions thater later
- * access this data in the public interfaces of the FEValues and related
- * classes. Consequently, this base class is protected in FEValuesBase.
- *
- * The second reason is because in FEValuesBase::reinit, we first need to
- * call Mapping::fill_fe_values() to compute mapping related data, and later
- * call FiniteElement::fill_fe_values() to compute shape function related
- * data. In the first step, Mapping::fill_fe_values() gets a pointer to
- * its own internal data structure and a pointer to the FEValuesData base
- * object of FEValuesBase, and the mapping then places the computed data
- * into the data fields that pertain to the mapping below. In the second
- * step, the finite element receives a pointer to its own internal object,
- * and to the current object, and from both of these computes the shape
- * function related information and, again, places it into the current
- * FEValuesData object.
- *
* More information can be found on the page on
* @ref UpdateFlagsEssay.
*
* @ingroup feaccess
- * @author Guido Kanschat
- * @date 2000
*/
template <int dim, int spacedim=dim>
-class FEValuesData
+class FEValuesData : public internal::FEValues::MappingRelatedData<dim,spacedim>,
+ public internal::FEValues::FiniteElementRelatedData<dim,spacedim>
{
public:
/**
const FiniteElement<dim,spacedim> &fe,
const UpdateFlags flags);
- /**
- * @name Fields filled by the finite element
- * @{
- */
-
- /**
- * Storage type for shape values. Each row in the matrix denotes the values
- * of a single shape function at the different points, columns are for a
- * single point with the different shape functions.
- *
- * If a shape function has more than one non-zero component (in deal.II
- * diction: it is non-primitive), then we allocate one row per non-zero
- * component, and shift subsequent rows backward. Lookup of the correct row
- * for a shape function is thus simple in case the entire finite element is
- * primitive (i.e. all shape functions are primitive), since then the shape
- * function number equals the row number. Otherwise, use the
- * #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 FiniteElement::get_nonzero_components()
- * function which tells us which components are non-zero and thus have a row
- * in the array presently under discussion.
- */
- typedef Table<2,double> ShapeVector;
-
- /**
- * Storage type for gradients. The layout of data is the same as for the
- * #ShapeVector data type.
- */
- typedef std::vector<std::vector<Tensor<1,spacedim> > > GradientVector;
-
- /**
- * Likewise for second order derivatives.
- */
- typedef std::vector<std::vector<Tensor<2,spacedim> > > HessianVector;
-
- /**
- * Store the values of the shape functions at the quadrature points. See the
- * description of the data type for the layout of the data in this field.
- */
- ShapeVector shape_values;
-
- /**
- * Store the gradients of the shape functions at the quadrature points. See
- * the description of the data type for the layout of the data in this
- * field.
- */
- GradientVector shape_gradients;
-
- /**
- * Store the 2nd derivatives of the shape functions at the quadrature
- * points. See the description of the data type for the layout of the data
- * in this field.
- */
- HessianVector shape_hessians;
-
- /**
- * @}
- */
-
- /**
- * @name Fields filled by the mapping
- * @{
- */
-
- /**
- * Store an array of weights times the Jacobi determinant at the quadrature
- * points. This function is reset each time reinit() is called. The Jacobi
- * determinant is actually the reciprocal value of the Jacobi matrices
- * stored in this class, see the general documentation of this class for
- * more information.
- *
- * However, if this object refers to an FEFaceValues or FESubfaceValues
- * object, then the JxW_values correspond to the Jacobian of the
- * transformation of the face, not the cell, i.e. the dimensionality is that
- * of a surface measure, not of a volume measure. In this case, it is
- * computed from the boundary forms, rather than the Jacobian matrix.
- */
- std::vector<double> JxW_values;
-
- /**
- * Array of the Jacobian matrices at the quadrature points.
- */
- std::vector< DerivativeForm<1,dim,spacedim> > jacobians;
-
- /**
- * Array of the derivatives of the Jacobian matrices at the quadrature
- * points.
- */
- std::vector<DerivativeForm<2,dim,spacedim> > jacobian_grads;
-
- /**
- * Array of the inverse Jacobian matrices at the quadrature points.
- */
- std::vector<DerivativeForm<1,spacedim,dim> > inverse_jacobians;
-
- /**
- * Array of quadrature points. This array is set up upon calling reinit()
- * and contains the quadrature points on the real element, rather than on
- * the reference element.
- */
- std::vector<Point<spacedim> > quadrature_points;
-
- /**
- * List of outward normal vectors at the quadrature points. This field is
- * filled in by the finite element class.
- */
- std::vector<Point<spacedim> > normal_vectors;
-
- /**
- * List of boundary forms at the quadrature points. This field is filled in
- * by the finite element class.
- */
- std::vector<Tensor<1,spacedim> > boundary_forms;
-
- /**
- * @}
- */
-
- /**
- * When asked for the value (or gradient, or Hessian) of shape function i's
- * c-th vector component, we need to look it up in the #shape_values,
- * #shape_gradients and #shape_hessians arrays. The question is where in
- * this array does the data for shape function i, component c reside. This
- * is what this table answers.
- *
- * The format of the table is as follows: - It has dofs_per_cell times
- * n_components entries. - The entry that corresponds to shape function i,
- * component c is <code>i * n_components + c</code>. - The value stored at
- * this position indicates the row in #shape_values and the other tables
- * where the corresponding datum is stored for all the quadrature points.
- *
- * In the general, vector-valued context, the number of components is larger
- * than one, but for a given shape function, not all vector components may
- * be nonzero (e.g., if a shape function is primitive, then exactly one
- * vector component is non-zero, while the others are all zero). For such
- * zero components, #shape_values and friends do not have a row.
- * Consequently, for vector components for which shape function i is zero,
- * the entry in the current table is numbers::invalid_unsigned_int.
- *
- * On the other hand, the table is guaranteed to have at least one valid
- * index for each shape function. In particular, for a primitive finite
- * element, each shape function has exactly one nonzero component and so for
- * each i, there is exactly one valid index within the range
- * <code>[i*n_components, (i+1)*n_components)</code>.
- */
- std::vector<unsigned int> shape_function_to_row_table;
-
/**
* Original update flags handed to the constructor of FEValues.
*/
* mapped quadrature points are accessible, as no finite element data is
* actually used).
*/
- const FEValues<dim> &get_fe_values () const;
+ const dealii::FEValues<dim> &get_fe_values () const;
/**
* Return a vector of inverse transpose Jacobians. For compatibility
/**
* An underlying FEValues object that performs the (scalar) evaluation.
*/
- FEValues<dim> fe_values;
+ dealii::FEValues<dim> fe_values;
/**
* Get 1D quadrature formula to be used for reinitializing shape info.
template <int dim, typename Number>
inline
- const FEValues<dim> &
+ const dealii::FEValues<dim> &
MappingDataOnTheFly<dim,Number>::get_fe_values() const
{
return fe_values;
const unsigned int my_q,
CellType (&cell_t_prev)[n_vector_elements],
CellType (&cell_t)[n_vector_elements],
- FEValues<dim,dim> &fe_values,
+ dealii::FEValues<dim,dim> &fe_values,
CellData &cell_data) const;
};
// hp::DoFHandler<dim>::active_cell_iterator, we need to manually
// select the correct finite element, so just hold a vector of
// FEValues
- std::vector<std_cxx11::shared_ptr<FEValues<dim> > >
+ std::vector<std_cxx11::shared_ptr<dealii::FEValues<dim> > >
fe_values (current_data.quadrature.size());
UpdateFlags update_flags_feval =
(update_flags & update_inverse_jacobians ? update_jacobians : update_default) |
const unsigned int n_q_points = current_data.n_q_points[fe_index];
if (fe_values[fe_index].get() == 0)
fe_values[fe_index].reset
- (new FEValues<dim> (mapping, dummy_fe,
- current_data.quadrature[fe_index],
- update_flags_feval));
- FEValues<dim> &fe_val = *fe_values[fe_index];
+ (new dealii::FEValues<dim> (mapping, dummy_fe,
+ current_data.quadrature[fe_index],
+ update_flags_feval));
+ dealii::FEValues<dim> &fe_val = *fe_values[fe_index];
data.resize (n_q_points);
// if the fe index has changed from the previous cell, set the
const unsigned int my_q,
CellType (&cell_t_prev)[n_vector_elements],
CellType (&cell_t)[n_vector_elements],
- FEValues<dim,dim> &fe_val,
+ dealii::FEValues<dim,dim> &fe_val,
CellData &data) const
{
const unsigned int n_q_points = fe_val.n_quadrature_points;
/* --------------------- FEValuesData ----------------- */
+namespace internal
+{
+ namespace FEValues
+ {
+ template <int dim, int spacedim>
+ void
+ MappingRelatedData<dim,spacedim>::initialize (const unsigned int n_quadrature_points,
+ const UpdateFlags flags)
+ {
+ if (flags & update_quadrature_points)
+ this->quadrature_points.resize(n_quadrature_points);
+
+ if (flags & update_JxW_values)
+ this->JxW_values.resize(n_quadrature_points);
+
+ if (flags & update_jacobians)
+ this->jacobians.resize(n_quadrature_points);
+
+ if (flags & update_jacobian_grads)
+ this->jacobian_grads.resize(n_quadrature_points);
+
+ if (flags & update_inverse_jacobians)
+ this->inverse_jacobians.resize(n_quadrature_points);
+
+ if (flags & update_boundary_forms)
+ this->boundary_forms.resize(n_quadrature_points);
+
+ if (flags & update_normal_vectors)
+ this->normal_vectors.resize(n_quadrature_points);
+ }
+
+
+ template <int dim, int spacedim>
+ void
+ FiniteElementRelatedData<dim,spacedim>::initialize (const unsigned int n_quadrature_points,
+ const FiniteElement<dim,spacedim> &fe,
+ const UpdateFlags flags)
+ {
+
+ // initialize the table mapping
+ // from shape function number to
+ // the rows in the tables storing
+ // the data by shape function and
+ // nonzero component
+ this->shape_function_to_row_table
+ = make_shape_function_to_row_table (fe);
+
+ // count the total number of non-zero
+ // components accumulated over all shape
+ // functions
+ unsigned int n_nonzero_shape_components = 0;
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ n_nonzero_shape_components += fe.n_nonzero_components (i);
+ Assert (n_nonzero_shape_components >= fe.dofs_per_cell,
+ ExcInternalError());
+
+ // with the number of rows now
+ // known, initialize those fields
+ // that we will need to their
+ // correct size
+ if (flags & update_values)
+ this->shape_values.reinit(n_nonzero_shape_components,
+ n_quadrature_points);
+
+ if (flags & update_gradients)
+ this->shape_gradients.resize (n_nonzero_shape_components,
+ std::vector<Tensor<1,spacedim> > (n_quadrature_points));
+
+ if (flags & update_hessians)
+ this->shape_hessians.resize (n_nonzero_shape_components,
+ std::vector<Tensor<2,spacedim> > (n_quadrature_points));
+ }
+ }
+}
+
+
template <int dim, int spacedim>
void
FEValuesData<dim,spacedim>::initialize (const unsigned int n_quadrature_points,
const FiniteElement<dim,spacedim> &fe,
const UpdateFlags flags)
{
- this->update_flags = flags;
+ // initialize the base classes
+ internal::FEValues::MappingRelatedData<dim,spacedim>::initialize(n_quadrature_points, flags);
+ internal::FEValues::FiniteElementRelatedData<dim,spacedim>::initialize(n_quadrature_points, fe, flags);
- // initialize the table mapping
- // from shape function number to
- // the rows in the tables storing
- // the data by shape function and
- // nonzero component
- this->shape_function_to_row_table
- = make_shape_function_to_row_table (fe);
-
- // count the total number of non-zero
- // components accumulated over all shape
- // functions
- unsigned int n_nonzero_shape_components = 0;
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- n_nonzero_shape_components += fe.n_nonzero_components (i);
- Assert (n_nonzero_shape_components >= fe.dofs_per_cell,
- ExcInternalError());
-
- // with the number of rows now
- // known, initialize those fields
- // that we will need to their
- // correct size
- if (flags & update_values)
- this->shape_values.reinit(n_nonzero_shape_components,
- n_quadrature_points);
-
- if (flags & update_gradients)
- this->shape_gradients.resize (n_nonzero_shape_components,
- std::vector<Tensor<1,spacedim> > (n_quadrature_points));
-
- if (flags & update_hessians)
- this->shape_hessians.resize (n_nonzero_shape_components,
- std::vector<Tensor<2,spacedim> > (n_quadrature_points));
-
- if (flags & update_quadrature_points)
- this->quadrature_points.resize(n_quadrature_points);
-
- if (flags & update_JxW_values)
- this->JxW_values.resize(n_quadrature_points);
-
- if (flags & update_jacobians)
- this->jacobians.resize(n_quadrature_points);
-
- if (flags & update_jacobian_grads)
- this->jacobian_grads.resize(n_quadrature_points);
-
- if (flags & update_inverse_jacobians)
- this->inverse_jacobians.resize(n_quadrature_points);
-
- if (flags & update_boundary_forms)
- this->boundary_forms.resize(n_quadrature_points);
-
- if (flags & update_normal_vectors)
- this->normal_vectors.resize(n_quadrature_points);
+ this->update_flags = flags;
}