#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/evaluated_geometry.h>
+#include <deal.II/matrix_free/mapping_fe_evaluation.h>
DEAL_II_NAMESPACE_OPEN
*/
internal::MatrixFreeFunctions::CellType get_cell_type() const;
+ /**
+ * Returns a reference to the ShapeInfo object currently in use.
+ */
+ const internal::MatrixFreeFunctions::ShapeInfo<Number> &
+ get_shape_info() const;
+
+ /**
+ * Fills the JxW values currently used.
+ */
+ void
+ fill_JxW_values(AlignedVector<VectorizedArray<Number> > &JxW_values) const;
+
//@}
/**
* ConstraintMatrix::distribute. When accessing vector entries during the
* 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).
*/
template <typename VectorType>
void read_dof_values (const VectorType &src);
* function is also necessary when inhomogeneous constraints are to be used,
* 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.
*/
template <typename VectorType>
void read_dof_values_plain (const VectorType &src);
* Takes the values stored internally on dof values of the current cell and
* sums them into the vector @p dst. The function also applies constraints
* during the write operation. The functionality is hence similar to the
- * function ConstraintMatrix::distribute_local_to_global. Note that if
- * vectorization is enabled, the DoF values for several cells are used.
+ * 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.
*/
template<typename VectorType>
void distribute_local_to_global (VectorType &dst) const;
* and sums them into the collection of vectors vector @p dst, starting at
* index @p first_index. The function also applies constraints during the
* write operation. The functionality is hence similar to the function
- * ConstraintMatrix::distribute_local_to_global. Note that if vectorization
- * is enabled, the DoF values for several cells are used.
+ * ConstraintMatrix::distribute_local_to_global. If vectorization is
+ * enabled, the DoF values for several cells are used.
*/
template<typename VectorType>
void distribute_local_to_global (std::vector<VectorType> &dst,
* during the write operation. The functionality is hence similar to the
* 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.
*/
template<typename VectorType>
void set_dof_values (VectorType &dst) const;
/**
* Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type EvaluatedGeometry
+ * 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
* operations), but only possibly within the element if the
* evaluate/integrate routines are combined (e.g. for matrix assembly).
*/
- FEEvaluationBase (const EvaluatedGeometry<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int base_element = 0);
+ FEEvaluationBase (const MappingFEEvaluation<dim,Number> &geometry,
+ const DoFHandler<dim> &dof_handler,
+ const unsigned int first_selected_component = 0);
/**
* Copy constructor
const internal::MatrixFreeFunctions::MappingInfo<dim,Number> *mapping_info;
/**
- * Stores a reference to the unit cell data, i.e., values, gradients and
+ * In case the class is initialized from MappingFEEvaluation instead of
+ * MatrixFree, this data structure holds the evaluated shape data.
+ */
+ std_cxx1x::shared_ptr<internal::MatrixFreeFunctions::ShapeInfo<Number> > stored_shape_info;
+
+ /**
+ * Stores a pointer to the unit cell shape data, i.e., values, gradients and
* Hessians in 1D at the quadrature points that constitute the tensor
* product. Also contained in matrix_info, but it simplifies code if we
- * store a reference to it. If the object is initialized without MatrixFree
- * object, the constructor creates this data structure.
+ * store a reference to it.
*/
- std_cxx1x::shared_ptr<const internal::MatrixFreeFunctions::ShapeInfo<Number> > data;
+ const internal::MatrixFreeFunctions::ShapeInfo<Number> *data;
/**
* A pointer to the Cartesian Jacobian information of the present cell. Only
/**
* Geometry data generated by FEValues on the fly.
*/
- SmartPointer<const EvaluatedGeometry<dim,Number> > evaluated_geometry;
+ SmartPointer<const MappingFEEvaluation<dim,Number> > mapped_geometry;
/**
* A pointer to the underlying DoFHandler.
*/
const DoFHandler<dim> *dof_handler;
+
+ /**
+ * For a DoFHandler with more than one finite element, select at which
+ * component this data structure should start.
+ */
+ const unsigned int first_selected_component;
+
+ /**
+ * A temporary data structure necessary to read degrees of freedom when no
+ * MatrixFree object was given at initialization.
+ */
+ mutable std::vector<types::global_dof_index> local_dof_indices;
};
/**
* Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type EvaluatedGeometry
+ * 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
* operations), but only possibly within the element if the
* evaluate/integrate routines are combined (e.g. for matrix assembly).
*/
- FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int base_element = 0);
+ FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
+ const DoFHandler<dim> &dof_handler,
+ const unsigned int first_selected_component = 0);
/**
* Copy constructor
/**
* Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type EvaluatedGeometry
+ * 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
* operations), but only possibly within the element if the
* evaluate/integrate routines are combined (e.g. for matrix assembly).
*/
- FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int base_element = 0);
+ FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
+ const DoFHandler<dim> &dof_handler,
+ const unsigned int first_selected_component = 0);
/**
* Copy constructor
/**
* Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type EvaluatedGeometry
+ * 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
* operations), but only possibly within the element if the
* evaluate/integrate routines are combined (e.g. for matrix assembly).
*/
- FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int base_element = 0);
+ FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
+ const DoFHandler<dim> &dof_handler,
+ const unsigned int first_selected_component = 0);
/**
* Copy constructor
/**
* Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type EvaluatedGeometry
+ * 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
* operations), but only possibly within the element if the
* evaluate/integrate routines are combined (e.g. for matrix assembly).
*/
- FEEvaluationGeneral (const EvaluatedGeometry<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int base_element = 0);
+ FEEvaluationGeneral (const MappingFEEvaluation<dim,Number> &geometry,
+ const DoFHandler<dim> &dof_handler,
+ const unsigned int first_selected_component = 0);
/**
* Copy constructor
* 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
- * EvaluatedGeometry. Here, the operations can only work on a single cell, but
+ * 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).
*
/**
* Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type EvaluatedGeometry
+ * 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
* operations), but only possibly within the element if the
* evaluate/integrate routines are combined (e.g. for matrix assembly).
*/
- FEEvaluation (const EvaluatedGeometry<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int base_element = 0);
+ FEEvaluation (const MappingFEEvaluation<dim,Number> &geometry,
+ const DoFHandler<dim> &dof_handler,
+ const unsigned int first_selected_component = 0);
/**
* Evaluates the function values, the gradients, and the Laplacians of the
/**
* Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type EvaluatedGeometry
+ * 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
* operations), but only possibly within the element if the
* evaluate/integrate routines are combined (e.g. for matrix assembly).
*/
- FEEvaluationGL (const EvaluatedGeometry<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int base_element = 0);
+ FEEvaluationGL (const MappingFEEvaluation<dim,Number> &geometry,
+ const DoFHandler<dim> &dof_handler,
+ const unsigned int first_selected_component = 0);
/**
* Copy constructor
/**
* Constructor that comes with reduced functionality and works similar as
- * FEValues. The user has to provide a structure of type EvaluatedGeometry
+ * 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
* operations), but only possibly within the element if the
* evaluate/integrate routines are combined (e.g. for matrix assembly).
*/
- FEEvaluationDGP (const EvaluatedGeometry<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int base_element = 0);
+ FEEvaluationDGP (const MappingFEEvaluation<dim,Number> &geometry,
+ const DoFHandler<dim> &dof_handler,
+ const unsigned int first_selected_component = 0);
/**
* Copy constructor
#ifndef DOXYGEN
-namespace internal
-{
- namespace MatrixFreeFunctions
- {
- // a small class that gives control over the delete behavior of
- // std::shared_ptr: we need to disable it when we initialize a pointer
- // from another structure.
- template <typename CLASS>
- struct DummyDeleter
- {
- DummyDeleter (const bool do_delete = false)
- :
- do_delete(do_delete)
- {}
-
- void operator () (CLASS *pointer)
- {
- if (do_delete)
- delete pointer;
- }
-
- const bool do_delete;
- };
- }
-}
-
-
/*----------------------- FEEvaluationBase ----------------------------------*/
mapping_info (&data_in.get_mapping_info()),
data (&data_in.get_shape_info
(fe_no_in, quad_no_in, active_fe_index,
- active_quad_index),
- internal::MatrixFreeFunctions::DummyDeleter
- <const internal::MatrixFreeFunctions::ShapeInfo<Number> >(false)),
+ active_quad_index)),
cartesian_data (0),
jacobian (0),
J_value (0),
jacobian_grad_upper(0),
cell (numbers::invalid_unsigned_int),
cell_type (internal::MatrixFreeFunctions::undefined),
- cell_data_number (0)
+ 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>
inline
FEEvaluationBase<dim,n_components_,Number>
-::FEEvaluationBase (const EvaluatedGeometry<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler_in,
- const unsigned int base_element)
+::FEEvaluationBase (const MappingFEEvaluation<dim,Number> &geometry,
+ const DoFHandler<dim> &dof_handler_in,
+ const unsigned int first_selected_component)
:
quad_no (-1),
n_fe_components (n_components_),
matrix_info (0),
dof_info (0),
mapping_info (0),
- data (new internal::MatrixFreeFunctions::ShapeInfo<Number>(geometry.get_quadrature(), dof_handler_in.get_fe(), base_element)),
+ stored_shape_info (new internal::MatrixFreeFunctions::ShapeInfo<Number>(geometry.get_quadrature(), dof_handler_in.get_fe(), 0)),
+ data (stored_shape_info.get()),
cartesian_data (0),
jacobian (geometry.get_inverse_jacobians().begin()),
J_value (geometry.get_JxW_values().begin()),
cell (0),
cell_type (internal::MatrixFreeFunctions::general),
cell_data_number (0),
- evaluated_geometry (&geometry),
- dof_handler (&dof_handler_in)
+ mapped_geometry (&geometry),
+ dof_handler (&dof_handler_in),
+ first_selected_component (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;
}
- Assert(dof_handler->get_fe().element_multiplicity(base_element) == 1 ||
- dof_handler->get_fe().element_multiplicity(base_element) >= n_components_,
+ Assert(dof_handler->get_fe().element_multiplicity(0) == 1 ||
+ dof_handler->get_fe().element_multiplicity(0)-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),
- evaluated_geometry (other.evaluated_geometry),
- dof_handler (other.dof_handler)
+ 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)
{
-
template <int dim, int n_components_, typename Number>
inline
void
-FEEvaluationBase<dim,n_components_,Number>
-::reinit (const unsigned int cell_in)
+FEEvaluationBase<dim,n_components_,Number>::reinit (const unsigned int cell_in)
{
- Assert (evaluated_geometry == 0, ExcMessage("FEEvaluation was initialized without a matrix-free object. Integer indexing is not possible"));
- if (evaluated_geometry != 0)
+ 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());
Assert (mapping_info != 0, ExcNotInitialized());
template <int dim, int n_components_, typename Number>
inline
internal::MatrixFreeFunctions::CellType
-FEEvaluationBase<dim,n_components_,Number>
-::get_cell_type () const
+FEEvaluationBase<dim,n_components_,Number>::get_cell_type () const
{
Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
return cell_type;
+template <int dim, int n_components_, typename Number>
+inline
+const internal::MatrixFreeFunctions::ShapeInfo<Number> &
+FEEvaluationBase<dim,n_components_,Number>::get_shape_info() const
+{
+ Assert(data != 0, ExcInternalError());
+ return *data;
+}
+
+
+
+template <int dim, int n_components_, typename Number>
+inline
+void
+FEEvaluationBase<dim,n_components_,Number>
+::fill_JxW_values(AlignedVector<VectorizedArray<Number> > &JxW_values) const
+{
+ AssertDimension(JxW_values.size(), data->n_q_points);
+ Assert (this->J_value != 0, ExcNotImplemented());
+ if (this->cell_type == internal::MatrixFreeFunctions::cartesian ||
+ this->cell_type == internal::MatrixFreeFunctions::affine)
+ {
+ Assert (this->mapping_info != 0, ExcNotImplemented());
+ VectorizedArray<Number> J = this->J_value[0];
+ for (unsigned int q=0; q<this->data->n_q_points; ++q)
+ JxW_values[q] = J * this->quadrature_weights[q];
+ }
+ else
+ for (unsigned int q=0; q<data->n_q_points; ++q)
+ JxW_values[q] = this->J_value[q];
+}
+
+
+
namespace internal
{
// write access to generic vectors that have operator ().
res = vector_access (const_cast<const VectorType &>(vec), index);
}
+ template <typename VectorType>
+ void process_dof_global (const types::global_dof_index index,
+ VectorType &vec,
+ Number &res) const
+ {
+ res = const_cast<const VectorType &>(vec)(index);
+ }
+
void pre_constraints (const Number &,
Number &res) const
{
vector_access (vec, index) += res;
}
+ template <typename VectorType>
+ void process_dof_global (const types::global_dof_index index,
+ VectorType &vec,
+ Number &res) const
+ {
+ vec(index) += res;
+ }
+
void pre_constraints (const Number &input,
Number &res) const
{
vector_access (vec, index) = res;
}
+ template <typename VectorType>
+ void process_dof_global (const types::global_dof_index index,
+ VectorType &vec,
+ Number &res) const
+ {
+ vec(index) = res;
+ }
+
void pre_constraints (const Number &,
Number &) const
{
// into the local data field or write local data into the vector. Certain
// operations are no-ops for the given use case.
- Assert (matrix_info != 0, ExcNotInitialized());
+ // Case 1: No MatrixFree object given, simple case because we do not need to
+ // 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);
+
+ unsigned int index = first_selected_component * this->data->dofs_per_cell;
+ for (unsigned int comp = 0; comp<n_components; ++comp)
+ {
+ for (unsigned int i=0; i<this->data->dofs_per_cell; ++i, ++index)
+ {
+ operation.process_dof_global(local_dof_indices[this->data->lexicographic_numbering[index]],
+ *src[0], values_dofs[comp][i][0]);
+ for (unsigned int v=1; v<VectorizedArray<Number>::n_array_elements; ++v)
+ operation.process_empty(values_dofs[comp][i][v]);
+ }
+ }
+ return;
+ }
+
Assert (dof_info != 0, ExcNotInitialized());
Assert (matrix_info->indices_initialized() == true,
ExcNotInitialized());
FEEvaluationBase<dim,n_components_,Number>
::read_dof_values_plain (const VectorType *src[])
{
+ // Case without MatrixFree initialization object
+ if (matrix_info == 0)
+ {
+ internal::VectorReader<Number> reader;
+ read_write_operation (reader, src);
+ return;
+ }
+
// this is different from the other three operations because we do not use
// constraints here, so this is a separate function.
- Assert (matrix_info != 0, ExcNotInitialized());
Assert (dof_info != 0, ExcNotInitialized());
Assert (matrix_info->indices_initialized() == true,
ExcNotInitialized());
template <int dim, int n_components_, typename Number>
inline
FEEvaluationAccess<dim,n_components_,Number>
-::FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
const DoFHandler<dim> &dof_handler,
- const unsigned int base_element)
+ const unsigned int first_selected_component)
:
- FEEvaluationBase <dim,n_components_,Number> (geometry, dof_handler, base_element)
+ FEEvaluationBase <dim,n_components_,Number> (geometry, dof_handler, first_selected_component)
{}
template <int dim, typename Number>
inline
FEEvaluationAccess<dim,1,Number>
-::FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
const DoFHandler<dim> &dof_handler,
- const unsigned int base_element)
+ const unsigned int first_selected_component)
:
- FEEvaluationBase <dim,1,Number> (geometry, dof_handler, base_element)
+ FEEvaluationBase <dim,1,Number> (geometry, dof_handler, first_selected_component)
{}
template <int dim, typename Number>
inline
FEEvaluationAccess<dim,dim,Number>
-::FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
const DoFHandler<dim> &dof_handler,
- const unsigned int base_element)
+ const unsigned int first_selected_component)
:
- FEEvaluationBase <dim,dim,Number> (geometry, dof_handler, base_element)
+ FEEvaluationBase <dim,dim,Number> (geometry, dof_handler, first_selected_component)
{}
typename Number>
inline
FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluationGeneral (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluationGeneral (const MappingFEEvaluation<dim,Number> &geometry,
const DoFHandler<dim> &dof_handler,
- const unsigned int base_element)
+ const unsigned int first_selected_component)
:
- BaseClass (geometry, dof_handler, base_element)
+ BaseClass (geometry, dof_handler, first_selected_component)
{
set_data_pointers();
}
typename Number>
inline
FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluation (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluation (const MappingFEEvaluation<dim,Number> &geometry,
const DoFHandler<dim> &dof_handler,
- const unsigned int base_element)
+ const unsigned int first_selected_component)
:
- BaseClass (geometry, dof_handler, base_element)
+ BaseClass (geometry, dof_handler, first_selected_component)
{
compute_even_odd_factors();
}
template <int dim, int fe_degree, int n_components_, typename Number>
inline
FEEvaluationGL<dim,fe_degree,n_components_,Number>
-::FEEvaluationGL (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluationGL (const MappingFEEvaluation<dim,Number> &geometry,
const DoFHandler<dim> &dof_handler,
- const unsigned int base_element)
+ const unsigned int first_selected_component)
:
- BaseClass (geometry, dof_handler, base_element)
+ BaseClass (geometry, dof_handler, first_selected_component)
{}
typename Number>
inline
FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluationDGP (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluationDGP (const MappingFEEvaluation<dim,Number> &geometry,
const DoFHandler<dim> &dof_handler,
- const unsigned int base_element)
+ const unsigned int first_selected_component)
:
- BaseClass (geometry, dof_handler, base_element)
+ BaseClass (geometry, dof_handler, first_selected_component)
{}
// ---------------------------------------------------------------------
-#ifndef __deal2__matrix_free_evaluated_geometry_h
-#define __deal2__matrix_free_evaluated_geometry_h
+#ifndef __deal2__matrix_free_mapping_fe_evaluation_h
+#define __deal2__matrix_free_mapping_fe_evaluation_h
#include <deal.II/base/config.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/mapping_info.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_nothing.h>
DEAL_II_NAMESPACE_OPEN
/**
- * The class makes FEEvaluation with the mapping information generated by
- * FEValues.
+ * This class provides evaluated mapping information using standard deal.II
+ * information in a form that FEEvaluation and friends can use for vectorized
+ * access. Since no vectorization over cells is available with the
+ * DoFHandler/Triangulation cell iterators, the interface to FEEvaluation's
+ * vectorization model is to use @p VectorizedArray::n_array_element copies of
+ * the same element. This interface is thus primarily useful for evaluating
+ * several operators on the same cell, e.g., when assembling cell matrices.
+ *
+ * @author Martin Kronbichler
*/
template <int dim, typename Number=double>
-class EvaluatedGeometry : public Subscriptor
+class MappingFEEvaluation : public Subscriptor
{
public:
/**
- * Constructor, similar to FEValues.
+ * Constructor, similar to FEValues. Since this class only evaluates the
+ * geometry, no finite element has to be specified and the simplest element,
+ * FE_Nothing, is used internally for the underlying FEValues object.
*/
- EvaluatedGeometry (const Mapping<dim> &mapping,
- const FiniteElement<dim> &fe,
+ MappingFEEvaluation (const Mapping<dim> &mapping,
const Quadrature<1> &quadrature,
const UpdateFlags update_flags);
/**
* Constructor. Instead of providing a mapping, use MappingQ1.
*/
- EvaluatedGeometry (const FiniteElement<dim> &fe,
- const Quadrature<1> &quadrature,
+ MappingFEEvaluation (const Quadrature<1> &quadrature,
const UpdateFlags update_flags);
/**
- * Initialize with the given cell iterator.
+ * Initialize with the given cell iterator. This works for all type of
+ * iterators that can be converted into a Triangulation::cell_iterator
+ * (DoFHandler::cell_iterator or Triangulation::cell_iterator).
*/
template <typename ITERATOR>
- void reinit(ITERATOR &cell);
+ void reinit(ITERATOR cell);
+
+ /**
+ * Return a triangulation iterator to the current cell.
+ */
+ typename Triangulation<dim>::cell_iterator get_cell () const;
+
+ /**
+ * Return a reference to the underlying FEValues object that evaluates
+ * certain quantities (only mapping-related ones like Jacobians or mapped
+ * quadrature points are accessible, as no finite element data is actually
+ * used).
+ */
+ const FEValues<dim>& get_fe_values () const;
/**
* Return a vector of inverse transpose Jacobians. For compatibility with
get_quadrature () const;
private:
+ /**
+ * Dummy finite element object necessary for initializing the FEValues
+ * object.
+ */
+ FE_Nothing<dim> fe_dummy;
+
/**
* An underlying FEValues object that performs the (scalar) evaluation.
*/
template <int dim, typename Number>
inline
-EvaluatedGeometry<dim,Number>::EvaluatedGeometry (const Mapping<dim> &mapping,
- const FiniteElement<dim> &fe,
- const Quadrature<1> &quadrature,
- const UpdateFlags update_flags)
+MappingFEEvaluation<dim,Number>::MappingFEEvaluation (const Mapping<dim> &mapping,
+ const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags)
:
- fe_values(mapping, fe, Quadrature<dim>(quadrature),
+ fe_values(mapping, fe_dummy, Quadrature<dim>(quadrature),
internal::MatrixFreeFunctions::MappingInfo<dim,Number>::compute_update_flags(update_flags)),
quadrature_1d(quadrature),
inverse_jacobians(fe_values.get_quadrature().size()),
template <int dim, typename Number>
inline
-EvaluatedGeometry<dim,Number>::EvaluatedGeometry (const FiniteElement<dim> &fe,
- const Quadrature<1> &quadrature,
- const UpdateFlags update_flags)
+MappingFEEvaluation<dim,Number>::MappingFEEvaluation (const Quadrature<1> &quadrature,
+ const UpdateFlags update_flags)
:
- fe_values(fe, Quadrature<dim>(quadrature),
+ fe_values(fe_dummy, Quadrature<dim>(quadrature),
internal::MatrixFreeFunctions::MappingInfo<dim,Number>::compute_update_flags(update_flags)),
quadrature_1d(quadrature),
inverse_jacobians(fe_values.get_quadrature().size()),
template <typename ITERATOR>
inline
void
-EvaluatedGeometry<dim,Number>::reinit(ITERATOR &cell)
+MappingFEEvaluation<dim,Number>::reinit(ITERATOR cell)
{
- fe_values.reinit(cell);
+ fe_values.reinit(static_cast<typename Triangulation<dim>::cell_iterator>(cell));
for (unsigned int q=0; q<fe_values.get_quadrature().size(); ++q)
{
if (fe_values.get_update_flags() & update_inverse_jacobians)
+template <int dim, typename Number>
+inline
+typename Triangulation<dim>::cell_iterator
+MappingFEEvaluation<dim,Number>::get_cell() const
+{
+ return fe_values.get_cell();
+}
+
+
+
+template <int dim, typename Number>
+inline
+const FEValues<dim>&
+MappingFEEvaluation<dim,Number>::get_fe_values() const
+{
+ return fe_values;
+}
+
+
+
template <int dim, typename Number>
inline
const AlignedVector<Tensor<2,dim,VectorizedArray<Number> > >&
-EvaluatedGeometry<dim,Number>::get_inverse_jacobians() const
+MappingFEEvaluation<dim,Number>::get_inverse_jacobians() const
{
return inverse_jacobians;
}
template <int dim, typename Number>
inline
const AlignedVector<Tensor<1,dim,VectorizedArray<Number> > >&
-EvaluatedGeometry<dim,Number>::get_normal_vectors() const
+MappingFEEvaluation<dim,Number>::get_normal_vectors() const
{
return normal_vectors;
}
template <int dim, typename Number>
inline
const AlignedVector<Point<dim,VectorizedArray<Number> > >&
-EvaluatedGeometry<dim,Number>::get_quadrature_points() const
+MappingFEEvaluation<dim,Number>::get_quadrature_points() const
{
return quadrature_points;
}
template <int dim, typename Number>
inline
const AlignedVector<VectorizedArray<Number> >&
-EvaluatedGeometry<dim,Number>::get_JxW_values() const
+MappingFEEvaluation<dim,Number>::get_JxW_values() const
{
return jxw_values;
}
template <int dim, typename Number>
inline
const Quadrature<1>&
-EvaluatedGeometry<dim,Number>::get_quadrature() const
+MappingFEEvaluation<dim,Number>::get_quadrature() const
{
return quadrature_1d;
}
const std::vector<hp::QCollection<1> > &quad,
const typename MatrixFree<dim,Number>::AdditionalData additional_data)
{
+
+ // Reads out the FE information and stores the shape function values,
+ // gradients and Hessians for quadrature points.
+ {
+ const unsigned int n_fe = dof_handler.size();
+ const unsigned int n_quad = quad.size();
+ shape_info.reinit (TableIndices<4>(n_fe, n_quad, 1, 1));
+ for (unsigned int no=0; no<n_fe; no++)
+ for (unsigned int nq =0; nq<n_quad; nq++)
+ {
+ AssertDimension (quad[nq].size(), 1);
+ shape_info(no,nq,0,0).reinit(quad[nq][0], dof_handler[no]->get_fe());
+ }
+ }
+
if (additional_data.initialize_indices == true)
{
clear();
}
}
- // Reads out the FE information and stores the shape function values,
- // gradients and Hessians for quadrature points.
- const unsigned int n_fe = dof_handler.size();
- const unsigned int n_quad = quad.size();
- shape_info.reinit (TableIndices<4>(n_fe, n_quad, 1, 1));
- for (unsigned int no=0; no<n_fe; no++)
- {
- const FiniteElement<dim> &fe = dof_handler[no]->get_fe();
- for (unsigned int nq =0; nq<n_quad; nq++)
- {
- AssertDimension (quad[nq].size(), 1);
- shape_info(no,nq,0,0).reinit(quad[nq][0], fe.base_element(0));
- }
- }
-
// Evaluates transformations from unit to real cell, Jacobian determinants,
// quadrature points in real space, based on the ordering of the cells
// determined in @p extract_local_to_global_indices. The algorithm assumes
const std::vector<hp::QCollection<1> > &quad,
const typename MatrixFree<dim,Number>::AdditionalData additional_data)
{
+ // Reads out the FE information and stores the shape function values,
+ // gradients and Hessians for quadrature points.
+ {
+ const unsigned int n_components = dof_handler.size();
+ const unsigned int n_quad = quad.size();
+ unsigned int n_fe_in_collection = 0;
+ for (unsigned int i=0; i<n_components; ++i)
+ n_fe_in_collection = std::max (n_fe_in_collection,
+ dof_handler[i]->get_fe().size());
+ unsigned int n_quad_in_collection = 0;
+ for (unsigned int q=0; q<n_quad; ++q)
+ n_quad_in_collection = std::max (n_quad_in_collection, quad[q].size());
+ shape_info.reinit (TableIndices<4>(n_components, n_quad,
+ n_fe_in_collection,
+ n_quad_in_collection));
+ for (unsigned int no=0; no<n_components; no++)
+ for (unsigned int fe_no=0; fe_no<dof_handler[no]->get_fe().size(); ++fe_no)
+ for (unsigned int nq =0; nq<n_quad; nq++)
+ for (unsigned int q_no=0; q_no<quad[nq].size(); ++q_no)
+ shape_info(no,nq,fe_no,q_no).reinit (quad[nq][q_no],
+ dof_handler[no]->get_fe()[fe_no]);
+ }
+
if (additional_data.initialize_indices == true)
{
clear();
}
}
- // Reads out the FE information and stores the shape function values,
- // gradients and Hessians for quadrature points.
- const unsigned int n_components = dof_handler.size();
- const unsigned int n_quad = quad.size();
- unsigned int n_fe_in_collection = 0;
- for (unsigned int i=0; i<n_components; ++i)
- n_fe_in_collection = std::max (n_fe_in_collection,
- dof_handler[i]->get_fe().size());
- unsigned int n_quad_in_collection = 0;
- for (unsigned int q=0; q<n_quad; ++q)
- n_quad_in_collection = std::max (n_quad_in_collection, quad[q].size());
- shape_info.reinit (TableIndices<4>(n_components, n_quad,
- n_fe_in_collection,
- n_quad_in_collection));
- for (unsigned int no=0; no<n_components; no++)
- for (unsigned int fe_no=0; fe_no<dof_handler[no]->get_fe().size(); ++fe_no)
- {
- const FiniteElement<dim> &fe = dof_handler[no]->get_fe()[fe_no];
- for (unsigned int nq =0; nq<n_quad; nq++)
- for (unsigned int q_no=0; q_no<quad[nq].size(); ++q_no)
- shape_info(no,nq,fe_no,q_no).reinit (quad[nq][q_no],
- fe.base_element(0));
- }
-
// Evaluates transformations from unit to real cell, Jacobian determinants,
// quadrature points in real space, based on the ordering of the cells
// determined in @p extract_local_to_global_indices.
template <int dim, typename Number>
void MatrixFree<dim,Number>::initialize_indices
(const std::vector<const ConstraintMatrix *> &constraint,
- const std::vector<IndexSet> &locally_owned_set)
+ const std::vector<IndexSet> &locally_owned_set)
{
const unsigned int n_fe = dof_handlers.n_dof_handlers;
const unsigned int n_active_cells = cell_level_index.size();
fes.back()->dofs_per_cell);
}
- lexicographic_inv[no].resize (fes.size());
for (unsigned int fe_index = 0; fe_index<fes.size(); ++fe_index)
{
const FiniteElement<dim> &fe = *fes[fe_index];
Assert (fe.n_base_elements() == 1,
- ExcMessage ("MatrixFree only works for DoFHandler with one base element"));
+ ExcMessage ("MatrixFree currently only works for DoFHandler with one base element"));
const unsigned int n_fe_components = fe.element_multiplicity (0);
// cache number of finite elements and dofs_per_cell
dof_info[no].dimension = dim;
dof_info[no].n_components = n_fe_components;
- // get permutation that gives lexicographic renumbering of the cell
- // dofs renumber (this is necessary for FE_Q, for example, since
- // there the vertex DoFs come first, which is incompatible with the
- // lexicographic ordering necessary to apply tensor products
- // efficiently)
-
- if (n_fe_components == 1)
- {
- lexicographic_inv[no][fe_index] =
- internal::MatrixFreeFunctions::get_lexicographic_numbering_inverse(fe);
- AssertDimension (lexicographic_inv[no][fe_index].size(),
- dof_info[no].dofs_per_cell[fe_index]);
- }
- else
- {
- // ok, we have more than one component, get the inverse
- // permutation, invert it, sort the components one after one,
- // and invert back
- Assert (n_fe_components > 1, ExcInternalError());
- std::vector<unsigned int> scalar_lex =
- Utilities::invert_permutation(internal::MatrixFreeFunctions::get_lexicographic_numbering_inverse(fe.base_element(0)));
- AssertDimension (scalar_lex.size() * n_fe_components,
- dof_info[no].dofs_per_cell[fe_index]);
- std::vector<unsigned int> lexicographic (dof_info[no].dofs_per_cell[fe_index]);
- for (unsigned int comp=0; comp<n_fe_components; ++comp)
- for (unsigned int i=0; i<scalar_lex.size(); ++i)
- lexicographic[fe.component_to_system_index(comp,i)]
- = scalar_lex.size () * comp + scalar_lex[i];
-
- // invert numbering
- lexicographic_inv[no][fe_index] =
- Utilities::invert_permutation(lexicographic);
- }
- AssertDimension (lexicographic_inv[no][fe_index].size(),
+ AssertDimension (shape_info(no,0,fe_index,0).lexicographic_numbering.size(),
dof_info[no].dofs_per_cell[fe_index]);
}
local_dof_indices.resize (dof_info[no].dofs_per_cell[0]);
cell_it->get_dof_indices(local_dof_indices);
dof_info[no].read_dof_indices (local_dof_indices,
- lexicographic_inv[no][0],
+ shape_info(no,0,0,0).lexicographic_numbering,
*constraint[no], counter,
constraint_values,
cell_at_boundary);
local_dof_indices.resize (dof_info[no].dofs_per_cell[0]);
cell_it->get_mg_dof_indices(local_dof_indices);
dof_info[no].read_dof_indices (local_dof_indices,
- lexicographic_inv[no][0],
+ shape_info(no,0,0,0).lexicographic_numbering,
*constraint[no], counter,
constraint_values,
cell_at_boundary);
local_dof_indices.resize (cell_it->get_fe().dofs_per_cell);
cell_it->get_dof_indices(local_dof_indices);
dof_info[no].read_dof_indices (local_dof_indices,
- lexicographic_inv[no][cell_it->active_fe_index()],
+ shape_info(no,0,cell_it->active_fe_index(),0).lexicographic_numbering,
*constraint[no], counter,
constraint_values,
cell_at_boundary);
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2011 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef __deal2__matrix_free_operators_h
+#define __deal2__matrix_free_operators_h
+
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace MatrixFreeOperators
+{
+ /**
+ * This class implements the operation of the action of the inverse of a mass
+ * matrix on an element for the special case of an evaluation object with as
+ * many quadrature points as there are cell degrees of freedom. It uses
+ * algorithms from FEEvaluation and produces the exact mass matrix for DGQ
+ * elements. This algorithm uses tensor products of inverse 1D shape matrices
+ * over quadrature points, so the inverse operation is exactly as expensive as
+ * applying the forward operator on each cell. Of course, for continuous
+ * finite elements this operation does not produce the inverse of a mass
+ * operation as the coupling between the elements cannot be considered by this
+ * operation.
+ *
+ * The equation may contain variable coefficients, so the user is required to
+ * provide an array for the inverse of the local coefficient (this class
+ * provide a helper method 'fill_inverse_JxW_values' to get the inverse of a
+ * constant-coefficient operator).
+ *
+ * @author Martin Kronbichler, 2014
+ */
+ template <int dim, int fe_degree, int n_components = 1, typename Number = double>
+ class CellwiseInverseMassMatrix
+ {
+ public:
+ /**
+ * Constructor. Initializes the shape information from the ShapeInfo field
+ * in the class FEEval.
+ */
+ CellwiseInverseMassMatrix (const FEEvaluationBase<dim,n_components,Number> &fe_eval);
+
+ /**
+ * Applies the inverse mass matrix operation on an input array. It is
+ * assumed that the passed input and output arrays are of correct size,
+ * namely FEEval::dofs_per_cell * n_components long. The inverse of the
+ * local coefficient (also containing the inverse JxW values) must be passed
+ * as first argument. Passing more than one component in the coefficient is
+ * allowed.
+ */
+ void apply(const AlignedVector<VectorizedArray<Number> > &inverse_coefficient,
+ const unsigned int n_actual_components,
+ const VectorizedArray<Number> *in_array,
+ VectorizedArray<Number> *out_array) const;
+
+ /**
+ * Fills the given array with the inverse of the JxW values, i.e., a mass
+ * matrix with coefficient 1. Non-unit coefficients must be multiplied (in
+ * inverse form) to this array.
+ */
+ void fill_inverse_JxW_values(AlignedVector<VectorizedArray<Number> > &inverse_jxw) const;
+
+ private:
+ /**
+ * A reference to the FEEvaluation object for getting the JxW_values.
+ */
+ const FEEvaluationBase<dim,n_components,Number> &fe_eval;
+
+ /**
+ * A structure to hold inverse shape functions, put into an aligned vector
+ * in order to make sure it can be allocated.
+ */
+ struct Inverse1DShape
+ {
+ VectorizedArray<Number> evenodd[fe_degree+1][fe_degree/2+1];
+ };
+ AlignedVector<Inverse1DShape> inverse_shape;
+ };
+
+
+
+ // ------------------------------------ inline functions ---------------------
+
+ template <int dim, int fe_degree, int n_components, typename Number>
+ inline
+ CellwiseInverseMassMatrix<dim,fe_degree,n_components,Number>
+ ::CellwiseInverseMassMatrix (const FEEvaluationBase<dim,n_components,Number> &fe_eval)
+ :
+ fe_eval (fe_eval)
+ {
+ FullMatrix<double> shapes_1d(fe_degree+1, fe_degree+1);
+ for (unsigned int i=0, c=0; i<shapes_1d.m(); ++i)
+ for (unsigned int j=0; j<shapes_1d.n(); ++j, ++c)
+ shapes_1d(i,j) = fe_eval.get_shape_info().shape_values_number[c];
+ shapes_1d.gauss_jordan();
+ inverse_shape.resize(1);
+ for (int i=0; i<(fe_degree+1)/2; ++i)
+ for (unsigned int q=0; q<(fe_degree+2)/2; ++q)
+ {
+ inverse_shape[0].evenodd[i][q] =
+ 0.5 * (shapes_1d(i,q) + shapes_1d(i,fe_degree-q));
+ inverse_shape[0].evenodd[fe_degree-i][q] =
+ 0.5 * (shapes_1d(i,q) - shapes_1d(i,fe_degree-q));
+ }
+ if (fe_degree % 2 == 0)
+ for (unsigned int q=0; q<(fe_degree+2)/2; ++q)
+ inverse_shape[0].evenodd[fe_degree/2][q] = shapes_1d(fe_degree/2,q);
+ }
+
+
+
+ template <int dim, int fe_degree, int n_components, typename Number>
+ inline
+ void
+ CellwiseInverseMassMatrix<dim,fe_degree,n_components,Number>
+ ::fill_inverse_JxW_values(AlignedVector<VectorizedArray<Number> > &inverse_jxw) const
+ {
+ const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
+ Assert(inverse_jxw.size() > 0 &&
+ inverse_jxw.size() % dofs_per_cell == 0,
+ ExcMessage("Expected diagonal to be a multiple of scalar dof per cells"));
+
+ // temporarily reduce size of inverse_jxw to dofs_per_cell to get JxW values
+ // from fe_eval (will not reallocate any memory)
+ const unsigned int previous_size = inverse_jxw.size();
+ inverse_jxw.resize(dofs_per_cell);
+ fe_eval.fill_JxW_values(inverse_jxw);
+
+ // invert
+ inverse_jxw.resize_fast(previous_size);
+ for (unsigned int q=0; q<dofs_per_cell; ++q)
+ inverse_jxw[q] = 1./inverse_jxw[q];
+ // copy values to rest of vector
+ for (unsigned int q=dofs_per_cell; q<previous_size; )
+ for (unsigned int i=0; i<dofs_per_cell; ++i, ++q)
+ inverse_jxw[q] = inverse_jxw[i];
+ }
+
+
+
+ template <int dim, int fe_degree, int n_components, typename Number>
+ inline
+ void
+ CellwiseInverseMassMatrix<dim,fe_degree,n_components,Number>
+ ::apply(const AlignedVector<VectorizedArray<Number> > &inverse_coefficients,
+ const unsigned int n_actual_components,
+ const VectorizedArray<Number> *in_array,
+ VectorizedArray<Number> *out_array) const
+ {
+ const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
+ Assert(inverse_coefficients.size() > 0 &&
+ inverse_coefficients.size() % dofs_per_cell == 0,
+ ExcMessage("Expected diagonal to be a multiple of scalar dof per cells"));
+ if (inverse_coefficients.size() != dofs_per_cell)
+ AssertDimension(n_actual_components * dofs_per_cell, inverse_coefficients.size());
+
+ Assert(dim == 2 || dim == 3, ExcNotImplemented());
+
+ const unsigned int shift_coefficient =
+ inverse_coefficients.size() > dofs_per_cell ? dofs_per_cell : 0;
+ const VectorizedArray<Number> *inv_coefficient = &inverse_coefficients[0];
+ VectorizedArray<Number> temp_data_field[dofs_per_cell];
+ for (unsigned int d=0; d<n_components; ++d)
+ {
+ const VectorizedArray<Number>* in = in_array+d*dofs_per_cell;
+ VectorizedArray<Number>* out = out_array+d*dofs_per_cell;
+ internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
+ VectorizedArray<Number>,0,false,false,2>
+ (inverse_shape[0].evenodd, in, temp_data_field);
+ internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
+ VectorizedArray<Number>,1,false,false,2>
+ (inverse_shape[0].evenodd, temp_data_field, out);
+
+ if (dim == 3)
+ {
+ internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
+ VectorizedArray<Number>,2,false,false,2>
+ (inverse_shape[0].evenodd, out, temp_data_field);
+ for (unsigned int q=0; q<dofs_per_cell; ++q)
+ temp_data_field[q] *= inv_coefficient[q];
+ internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
+ VectorizedArray<Number>,2,true,false,2>
+ (inverse_shape[0].evenodd, temp_data_field, out);
+ }
+ else if (dim == 2)
+ for (unsigned int q=0; q<dofs_per_cell; ++q)
+ out[q] *= inv_coefficient[q];
+
+ internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
+ VectorizedArray<Number>,1,true,false,2>
+ (inverse_shape[0].evenodd, out, temp_data_field);
+ internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
+ VectorizedArray<Number>,0,true,false,2>
+ (inverse_shape[0].evenodd, temp_data_field, out);
+
+ inv_coefficient += shift_coefficient;
+ }
+ }
+
+} // end of namespace MatrixFreeOperators
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
*/
std::vector<Number> shape_gradient_number;
+ /**
+ * Renumbering from deal.II's numbering of cell degrees of freedom to
+ * lexicographic numbering used inside the FEEvaluation schemes of the
+ * underlying element in the DoFHandler. For vector-valued elements, the
+ * renumbering starts with a lexicographic numbering of the first
+ * component, then everything of the second component, and so on.
+ */
+ std::vector<unsigned int> lexicographic_numbering;
+
/**
* Stores the number of quadrature points in
* @p dim dimensions for a cell.
};
- /**
- * Extracts the mapping between the actual numbering inside a scalar
- * element and lexicographic numbering as required by the evaluation
- * routines.
- */
- template <int dim>
- std::vector<unsigned int>
- get_lexicographic_numbering_inverse (const FiniteElement<dim> &fe);
-
-
// ------------------------------------------ inline functions
namespace MatrixFreeFunctions
{
- // helper function
- template <int dim>
- std::vector<unsigned int>
- get_lexicographic_numbering_inverse(const FiniteElement<dim> &fe)
- {
- Assert(fe.n_components() == 1,
- ExcMessage("Expected a scalar element"));
-
- const FE_Poly<TensorProductPolynomials<dim>,dim,dim> *fe_poly =
- dynamic_cast<const FE_Poly<TensorProductPolynomials<dim>,dim,dim>*>(&fe);
-
- const FE_Poly<TensorProductPolynomials<dim,Polynomials::
- PiecewisePolynomial<double> >,dim,dim> *fe_poly_piece =
- dynamic_cast<const FE_Poly<TensorProductPolynomials<dim,
- Polynomials::PiecewisePolynomial<double> >,dim,dim>*> (&fe);
-
- const FE_DGP<dim> *fe_dgp = dynamic_cast<const FE_DGP<dim>*>(&fe);
-
- std::vector<unsigned int> lexicographic;
- if (fe_poly != 0)
- lexicographic = fe_poly->get_poly_space_numbering_inverse();
- else if (fe_poly_piece != 0)
- lexicographic = fe_poly_piece->get_poly_space_numbering_inverse();
- else if (fe_dgp != 0)
- {
- lexicographic.resize(fe_dgp->dofs_per_cell);
- for (unsigned int i=0; i<fe_dgp->dofs_per_cell; ++i)
- lexicographic[i] = i;
- }
- else
- Assert(false, ExcNotImplemented());
-
- return lexicographic;
- }
-
-
-
// ----------------- actual ShapeInfo functions --------------------
template <typename Number>
const unsigned int n_dofs_1d = fe->degree+1,
n_q_points_1d = quad.size();
- std::vector<unsigned int> lexicographic (fe->dofs_per_cell);
// renumber (this is necessary for FE_Q, for example, since there the
// vertex DoFs come first, which is incompatible with the lexicographic
// ordering necessary to apply tensor products efficiently)
+ std::vector<unsigned int> scalar_lexicographic;
{
- lexicographic = get_lexicographic_numbering_inverse(*fe);
+ // find numbering to lexicographic
+ Assert(fe->n_components() == 1,
+ ExcMessage("Expected a scalar element"));
+
+ const FE_Poly<TensorProductPolynomials<dim>,dim,dim> *fe_poly =
+ dynamic_cast<const FE_Poly<TensorProductPolynomials<dim>,dim,dim>*>(fe);
+
+ const FE_Poly<TensorProductPolynomials<dim,Polynomials::
+ PiecewisePolynomial<double> >,dim,dim> *fe_poly_piece =
+ dynamic_cast<const FE_Poly<TensorProductPolynomials<dim,
+ Polynomials::PiecewisePolynomial<double> >,dim,dim>*> (fe);
+
+ const FE_DGP<dim> *fe_dgp = dynamic_cast<const FE_DGP<dim>*>(fe);
+
+ if (fe_poly != 0)
+ scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse();
+ else if (fe_poly_piece != 0)
+ scalar_lexicographic = fe_poly_piece->get_poly_space_numbering_inverse();
+ else if (fe_dgp != 0)
+ {
+ scalar_lexicographic.resize(fe_dgp->dofs_per_cell);
+ for (unsigned int i=0; i<fe_dgp->dofs_per_cell; ++i)
+ scalar_lexicographic[i] = i;
+ }
+ else
+ Assert(false, ExcNotImplemented());
+
+ // Finally store the renumbering into the member variable of this
+ // class
+ if (fe_in.n_components() == 1)
+ lexicographic_numbering = scalar_lexicographic;
+ else
+ {
+ // have more than one component, get the inverse
+ // permutation, invert it, sort the components one after one,
+ // and invert back
+ std::vector<unsigned int> scalar_inv =
+ Utilities::invert_permutation(scalar_lexicographic);
+ std::vector<unsigned int> lexicographic (fe_in.dofs_per_cell);
+ for (unsigned int comp=0; comp<fe_in.n_components(); ++comp)
+ for (unsigned int i=0; i<scalar_inv.size(); ++i)
+ lexicographic[fe_in.component_to_system_index(comp,i)]
+ = scalar_inv.size () * comp + scalar_inv[i];
+
+ // invert numbering again
+ lexicographic_numbering =
+ Utilities::invert_permutation(lexicographic);
+ }
// to evaluate 1D polynomials, evaluate along the line where y=z=0,
// assuming that shape_value(0,Point<dim>()) == 1. otherwise, need
// other entry point (e.g. generating a 1D element by reading the
// name, as done before r29356)
- Assert(std::fabs(fe->shape_value(lexicographic[0], Point<dim>())-1) < 1e-13,
+ Assert(std::fabs(fe->shape_value(scalar_lexicographic[0],
+ Point<dim>())-1) < 1e-13,
ExcInternalError());
}
{
// need to reorder from hierarchical to lexicographic to get the
// DoFs correct
- const unsigned int my_i = lexicographic[i];
+ const unsigned int my_i = scalar_lexicographic[i];
for (unsigned int q=0; q<n_q_points_1d; ++q)
{
// fill both vectors with
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests CellwiseInverseMassMatrix on DG elements by comparing its action on a
+// random vector to a CG solver
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/matrix_free/operators.h>
+
+
+std::ofstream logfile("output");
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+
+ MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+ data (data_in)
+ {};
+
+ void
+ local_mass_operator (const MatrixFree<dim,Number> &data,
+ VECTOR &dst,
+ const VECTOR &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const
+ {
+ FEEvaluation<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
+ const unsigned int n_q_points = fe_eval.n_q_points;
+
+ for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+ {
+ fe_eval.reinit (cell);
+ fe_eval.read_dof_values (src);
+ fe_eval.evaluate (true, false);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ fe_eval.submit_value (fe_eval.get_value(q),q);
+ fe_eval.integrate (true, false);
+ fe_eval.distribute_local_to_global (dst);
+ }
+ }
+
+ void
+ local_inverse_mass_operator (const MatrixFree<dim,Number> &data,
+ VECTOR &dst,
+ const VECTOR &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const
+ {
+ FEEvaluation<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
+ MatrixFreeOperators::CellwiseInverseMassMatrix<dim,fe_degree,1,Number> mass_inv(fe_eval);
+ const unsigned int n_q_points = fe_eval.n_q_points;
+ AlignedVector<VectorizedArray<Number> > inverse_coefficients(n_q_points);
+
+ for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+ {
+ fe_eval.reinit (cell);
+ mass_inv.fill_inverse_JxW_values(inverse_coefficients);
+ fe_eval.read_dof_values (src);
+ mass_inv.apply(inverse_coefficients, 1, fe_eval.begin_dof_values(),
+ fe_eval.begin_dof_values());
+ fe_eval.distribute_local_to_global (dst);
+ }
+ }
+
+ void vmult (VECTOR &dst,
+ const VECTOR &src) const
+ {
+ dst = 0;
+ data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_mass_operator,
+ this, dst, src);
+ };
+
+ void apply_inverse (VECTOR &dst,
+ const VECTOR &src) const
+ {
+ dst = 0;
+ data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_inverse_mass_operator,
+ this, dst, src);
+ };
+
+private:
+ const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof)
+{
+
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+ MappingQ<dim> mapping(4);
+
+ MatrixFree<dim,number> mf_data;
+ {
+ const QGauss<1> quad (fe_degree+1);
+ typename MatrixFree<dim,number>::AdditionalData data;
+ data.tasks_parallel_scheme =
+ MatrixFree<dim,number>::AdditionalData::partition_color;
+ data.tasks_block_size = 3;
+ ConstraintMatrix constraints;
+
+ mf_data.reinit (mapping, dof, constraints, quad, data);
+ }
+
+ MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
+ Vector<number> in (dof.n_dofs()), inverse (dof.n_dofs()), reference(dof.n_dofs());
+
+ for (unsigned int i=0; i<dof.n_dofs(); ++i)
+ {
+ const double entry = Testing::rand()/(double)RAND_MAX;
+ in(i) = entry;
+ }
+
+ mf.apply_inverse (inverse, in);
+
+ SolverControl control(1000, 1e-12);
+ SolverCG<Vector<number> > solver(control);
+ solver.solve (mf, reference, in, PreconditionIdentity());
+
+ inverse -= reference;
+ const double diff_norm = inverse.linfty_norm() / reference.linfty_norm();
+
+ deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_ball (tria);
+ static const HyperBallBoundary<dim> boundary;
+ tria.set_boundary (0, boundary);
+ if (dim < 3 || fe_degree < 2)
+ tria.refine_global(1);
+ tria.begin(tria.n_levels()-1)->set_refine_flag();
+ tria.last()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ typename Triangulation<dim>::active_cell_iterator
+ cell = tria.begin_active (),
+ endc = tria.end();
+ for (; cell!=endc; ++cell)
+ if (cell->center().norm()<1e-8)
+ cell->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+
+ FE_DGQ<dim> fe (fe_degree);
+ DoFHandler<dim> dof (tria);
+ dof.distribute_dofs(fe);
+
+ do_test<dim, fe_degree, double> (dof);
+
+ if (dim == 2)
+ {
+ deallog.push("float");
+ deallog.threshold_double(1.e-6);
+ do_test<dim, fe_degree, float> (dof);
+ deallog.threshold_double(5.e-11);
+ deallog.pop();
+ }
+}
+
+
+
+int main ()
+{
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ deallog << std::setprecision (3);
+
+ {
+ deallog.threshold_double(5.e-11);
+ deallog.push("2d");
+ test<2,1>();
+ test<2,2>();
+ test<2,4>();
+ deallog.pop();
+ deallog.push("3d");
+ test<3,1>();
+ test<3,2>();
+ deallog.pop();
+ }
+}
--- /dev/null
+
+DEAL:2d::Testing FE_DGQ<2>(1)
+DEAL:2d:cg::Starting value 6.32
+DEAL:2d:cg::Convergence step 51 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d:float::Testing FE_DGQ<2>(1)
+DEAL:2d:float:cg::Starting value 6.15
+DEAL:2d:float:cg::Convergence step 66 value 0
+DEAL:2d:float::Norm of difference: 0
+DEAL:2d:float::
+DEAL:2d::Testing FE_DGQ<2>(2)
+DEAL:2d:cg::Starting value 8.62
+DEAL:2d:cg::Convergence step 116 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d:float::Testing FE_DGQ<2>(2)
+DEAL:2d:float:cg::Starting value 9.05
+DEAL:2d:float:cg::Convergence step 158 value 0
+DEAL:2d:float::Norm of difference: 0
+DEAL:2d:float::
+DEAL:2d::Testing FE_DGQ<2>(4)
+DEAL:2d:cg::Starting value 14.6
+DEAL:2d:cg::Convergence step 354 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d:float::Testing FE_DGQ<2>(4)
+DEAL:2d:float:cg::Starting value 14.7
+DEAL:2d:float:cg::Convergence step 437 value 0
+DEAL:2d:float::Norm of difference: 0
+DEAL:2d:float::
+DEAL:3d::Testing FE_DGQ<3>(1)
+DEAL:3d:cg::Starting value 13.7
+DEAL:3d:cg::Convergence step 131 value 0
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_DGQ<3>(2)
+DEAL:3d:cg::Starting value 13.5
+DEAL:3d:cg::Convergence step 255 value 0
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests CellwiseInverseMassMatrix on vector DG elements, otherwise the same
+// as inverse_mass_01.cc
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/matrix_free/operators.h>
+
+
+std::ofstream logfile("output");
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+
+ MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+ data (data_in)
+ {};
+
+ void
+ local_mass_operator (const MatrixFree<dim,Number> &data,
+ VECTOR &dst,
+ const VECTOR &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const
+ {
+ FEEvaluation<dim,fe_degree,fe_degree+1,dim,Number> fe_eval (data);
+ const unsigned int n_q_points = fe_eval.n_q_points;
+
+ for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+ {
+ fe_eval.reinit (cell);
+ fe_eval.read_dof_values (src);
+ fe_eval.evaluate (true, false);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ fe_eval.submit_value (fe_eval.get_value(q),q);
+ fe_eval.integrate (true, false);
+ fe_eval.distribute_local_to_global (dst);
+ }
+ }
+
+ void
+ local_inverse_mass_operator (const MatrixFree<dim,Number> &data,
+ VECTOR &dst,
+ const VECTOR &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const
+ {
+ FEEvaluation<dim,fe_degree,fe_degree+1,dim,Number> fe_eval (data);
+ MatrixFreeOperators::CellwiseInverseMassMatrix<dim,fe_degree,dim,Number> mass_inv(fe_eval);
+ const unsigned int n_q_points = fe_eval.n_q_points;
+ AlignedVector<VectorizedArray<Number> > inverse_coefficients(n_q_points);
+
+ for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+ {
+ fe_eval.reinit (cell);
+ mass_inv.fill_inverse_JxW_values(inverse_coefficients);
+ fe_eval.read_dof_values (src);
+ mass_inv.apply(inverse_coefficients, dim, fe_eval.begin_dof_values(),
+ fe_eval.begin_dof_values());
+ fe_eval.distribute_local_to_global (dst);
+ }
+ }
+
+ void vmult (VECTOR &dst,
+ const VECTOR &src) const
+ {
+ dst = 0;
+ data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_mass_operator,
+ this, dst, src);
+ };
+
+ void apply_inverse (VECTOR &dst,
+ const VECTOR &src) const
+ {
+ dst = 0;
+ data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_inverse_mass_operator,
+ this, dst, src);
+ };
+
+private:
+ const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof)
+{
+
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+ MatrixFree<dim,number> mf_data;
+ {
+ const QGauss<1> quad (fe_degree+1);
+ typename MatrixFree<dim,number>::AdditionalData data;
+ data.tasks_parallel_scheme =
+ MatrixFree<dim,number>::AdditionalData::partition_color;
+ data.tasks_block_size = 3;
+ ConstraintMatrix constraints;
+
+ mf_data.reinit (dof, constraints, quad, data);
+ }
+
+ MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
+ Vector<number> in (dof.n_dofs()), inverse (dof.n_dofs()), reference(dof.n_dofs());
+
+ for (unsigned int i=0; i<dof.n_dofs(); ++i)
+ {
+ const double entry = Testing::rand()/(double)RAND_MAX;
+ in(i) = entry;
+ }
+
+ mf.apply_inverse (inverse, in);
+
+ SolverControl control(1000, 1e-12);
+ SolverCG<Vector<number> > solver(control);
+ solver.solve (mf, reference, in, PreconditionIdentity());
+
+ inverse -= reference;
+ const double diff_norm = inverse.linfty_norm() / reference.linfty_norm();
+
+ deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube (tria, -1, 1);
+ tria.refine_global(1);
+ if (dim < 3 || fe_degree < 2)
+ tria.refine_global(1);
+ tria.begin(tria.n_levels()-1)->set_refine_flag();
+ tria.last()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ typename Triangulation<dim>::active_cell_iterator
+ cell = tria.begin_active (),
+ endc = tria.end();
+ for (; cell!=endc; ++cell)
+ if (cell->center().norm()<1e-8)
+ cell->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+
+ FESystem<dim> fe (FE_DGQ<dim>(fe_degree), dim);
+ DoFHandler<dim> dof (tria);
+ dof.distribute_dofs(fe);
+
+ do_test<dim, fe_degree, double> (dof);
+}
+
+
+
+int main ()
+{
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ deallog << std::setprecision (3);
+
+ {
+ deallog.threshold_double(5.e-11);
+ deallog.push("2d");
+ test<2,1>();
+ test<2,2>();
+ test<2,4>();
+ deallog.pop();
+ deallog.push("3d");
+ test<3,1>();
+ test<3,2>();
+ deallog.pop();
+ }
+}
--- /dev/null
+
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^2]
+DEAL:2d:cg::Starting value 7.90
+DEAL:2d:cg::Convergence step 7 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2]
+DEAL:2d:cg::Starting value 11.6
+DEAL:2d:cg::Convergence step 15 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^2]
+DEAL:2d:cg::Starting value 19.1
+DEAL:2d:cg::Convergence step 46 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3]
+DEAL:3d:cg::Starting value 25.0
+DEAL:3d:cg::Convergence step 10 value 0
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3]
+DEAL:3d:cg::Starting value 24.3
+DEAL:3d:cg::Convergence step 32 value 0
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests CellwiseInverseMassMatrix on vector DG elements, similar test as
+// inverse_mass_02 but using different coefficients on different components
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/matrix_free/operators.h>
+
+
+std::ofstream logfile("output");
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+
+ MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+ data (data_in)
+ {};
+
+ void
+ local_mass_operator (const MatrixFree<dim,Number> &data,
+ VECTOR &dst,
+ const VECTOR &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const
+ {
+ FEEvaluation<dim,fe_degree,fe_degree+1,3,Number> fe_eval (data);
+ const unsigned int n_q_points = fe_eval.n_q_points;
+
+ for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+ {
+ fe_eval.reinit (cell);
+ fe_eval.read_dof_values (src);
+ fe_eval.evaluate (true, false);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ Tensor<1,3,VectorizedArray<double> > val = fe_eval.get_value(q);
+ val[0] *= make_vectorized_array(0.8314159);
+ val[1] *= make_vectorized_array(2.3);
+ val[2] *= make_vectorized_array(1.98);
+ fe_eval.submit_value (val, q);
+ }
+ fe_eval.integrate (true, false);
+ fe_eval.distribute_local_to_global (dst);
+ }
+ }
+
+ void
+ local_inverse_mass_operator (const MatrixFree<dim,Number> &data,
+ VECTOR &dst,
+ const VECTOR &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const
+ {
+ FEEvaluation<dim,fe_degree,fe_degree+1,3,Number> fe_eval (data);
+ MatrixFreeOperators::CellwiseInverseMassMatrix<dim,fe_degree,3,Number> mass_inv(fe_eval);
+ const unsigned int n_q_points = fe_eval.n_q_points;
+ AlignedVector<VectorizedArray<Number> > inverse_coefficients(3*n_q_points);
+
+ for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+ {
+ fe_eval.reinit (cell);
+ mass_inv.fill_inverse_JxW_values(inverse_coefficients);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ inverse_coefficients[q] *= make_vectorized_array(1./0.8314159);
+ inverse_coefficients[n_q_points+q] *= make_vectorized_array(1./2.3);
+ inverse_coefficients[2*n_q_points+q] *= make_vectorized_array(1./1.98);
+ }
+ fe_eval.read_dof_values (src);
+ mass_inv.apply(inverse_coefficients, 3, fe_eval.begin_dof_values(),
+ fe_eval.begin_dof_values());
+ fe_eval.distribute_local_to_global (dst);
+ }
+ }
+
+ void vmult (VECTOR &dst,
+ const VECTOR &src) const
+ {
+ dst = 0;
+ data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_mass_operator,
+ this, dst, src);
+ };
+
+ void apply_inverse (VECTOR &dst,
+ const VECTOR &src) const
+ {
+ dst = 0;
+ data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_inverse_mass_operator,
+ this, dst, src);
+ };
+
+private:
+ const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof)
+{
+
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+ MappingQ<dim> mapping(4);
+ MatrixFree<dim,number> mf_data;
+ {
+ const QGauss<1> quad (fe_degree+1);
+ typename MatrixFree<dim,number>::AdditionalData data;
+ data.tasks_parallel_scheme =
+ MatrixFree<dim,number>::AdditionalData::partition_color;
+ data.tasks_block_size = 3;
+ ConstraintMatrix constraints;
+
+ mf_data.reinit (mapping, dof, constraints, quad, data);
+ }
+
+ MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
+ Vector<number> in (dof.n_dofs()), inverse (dof.n_dofs()), reference(dof.n_dofs());
+
+ for (unsigned int i=0; i<dof.n_dofs(); ++i)
+ {
+ const double entry = Testing::rand()/(double)RAND_MAX;
+ in(i) = entry;
+ }
+
+ mf.apply_inverse (inverse, in);
+
+ SolverControl control(10000, 1e-12);
+ SolverCG<Vector<number> > solver(control);
+ solver.solve (mf, reference, in, PreconditionIdentity());
+
+ inverse -= reference;
+ const double diff_norm = inverse.linfty_norm() / reference.linfty_norm();
+
+ deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_ball (tria);
+ static const HyperBallBoundary<dim> boundary;
+ tria.set_boundary (0, boundary);
+ if (dim < 3 || fe_degree < 2)
+ tria.refine_global(1);
+ tria.begin(tria.n_levels()-1)->set_refine_flag();
+ tria.last()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ typename Triangulation<dim>::active_cell_iterator
+ cell = tria.begin_active (),
+ endc = tria.end();
+ for (; cell!=endc; ++cell)
+ if (cell->center().norm()<1e-8)
+ cell->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+
+ FESystem<dim> fe (FE_DGQ<dim>(fe_degree),3);
+ DoFHandler<dim> dof (tria);
+ dof.distribute_dofs(fe);
+
+ do_test<dim, fe_degree, double> (dof);
+}
+
+
+
+int main ()
+{
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ deallog << std::setprecision (3);
+
+ {
+ deallog.threshold_double(5.e-11);
+ deallog.push("2d");
+ test<2,1>();
+ test<2,2>();
+ test<2,4>();
+ deallog.pop();
+ deallog.push("3d");
+ test<3,1>();
+ test<3,2>();
+ deallog.pop();
+ }
+}
--- /dev/null
+
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^3]
+DEAL:2d:cg::Starting value 10.6
+DEAL:2d:cg::Convergence step 108 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^3]
+DEAL:2d:cg::Starting value 15.1
+DEAL:2d:cg::Convergence step 252 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^3]
+DEAL:2d:cg::Starting value 25.4
+DEAL:2d:cg::Convergence step 699 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3]
+DEAL:3d:cg::Starting value 23.7
+DEAL:3d:cg::Convergence step 412 value 0
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3]
+DEAL:3d:cg::Starting value 23.7
+DEAL:3d:cg::Convergence step 1019 value 0
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
template <int dim, int fe_degree, typename number>
void do_test (const DoFHandler<dim> &dof,
- const ConstraintMatrix &constraints,
- const unsigned int parallel_option = 0)
+ const ConstraintMatrix &constraints)
{
deallog << "Testing " << dof.get_fe().get_name() << std::endl;
- if (parallel_option > 0)
- deallog << "Parallel option: " << parallel_option << std::endl;
- //std::cout << "Number of cells: " << dof.get_tria().n_active_cells() << std::endl;
- //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
- //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl;
MatrixFree<dim,number> mf_data;
{
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this test is similar to matrix_vector_06, but implements the operations on
+// the fly from an FEValues implementation instead of data cached in
+// MatrixFree.
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/lac/vector.h>
+
+
+std::ofstream logfile("output");
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+ MatrixFreeTest(const DoFHandler<dim> &dof_handler,
+ const ConstraintMatrix &constraints)
+ :
+ dof_handler (dof_handler),
+ constraints (constraints)
+ {}
+
+ void vmult (VECTOR &dst,
+ const VECTOR &src) const
+ {
+ MappingFEEvaluation<dim,Number> mapping(QGauss<1>(fe_degree+1),
+ update_gradients | update_values |
+ update_JxW_values);
+ VECTOR src_cpy = src;
+ constraints.distribute(src_cpy);
+ FEEvaluation<dim,fe_degree,fe_degree+1,1,Number> fe_eval(mapping, dof_handler);
+ dst = 0;
+ typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for ( ; cell != endc; ++cell)
+ {
+ mapping.reinit(cell);
+ fe_eval.read_dof_values (src_cpy);
+ fe_eval.evaluate (true, true, false);
+ for (unsigned int q=0; q<fe_eval.n_q_points; ++q)
+ {
+ fe_eval.submit_value (Number(10)*fe_eval.get_value(q),q);
+ fe_eval.submit_gradient (fe_eval.get_gradient(q),q);
+ }
+ fe_eval.integrate (true,true);
+ fe_eval.distribute_local_to_global (dst);
+ }
+ constraints.condense(dst);
+ };
+
+private:
+ const DoFHandler<dim> &dof_handler;
+ const ConstraintMatrix &constraints;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof,
+ const ConstraintMatrix &constraints,
+ const unsigned int parallel_option = 0)
+{
+
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+ if (parallel_option > 0)
+ deallog << "Parallel option: " << parallel_option << std::endl;
+ //std::cout << "Number of cells: " << dof.get_tria().n_active_cells() << std::endl;
+ //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+ //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl;
+
+ MatrixFreeTest<dim,fe_degree,number> mf (dof, constraints);
+ Vector<number> in (dof.n_dofs()), out (dof.n_dofs());
+ Vector<number> in_dist (dof.n_dofs());
+ Vector<number> out_dist (in_dist);
+
+ for (unsigned int i=0; i<dof.n_dofs(); ++i)
+ {
+ if (constraints.is_constrained(i))
+ continue;
+ const double entry = Testing::rand()/(double)RAND_MAX;
+ in(i) = entry;
+ in_dist(i) = entry;
+ }
+
+ mf.vmult (out_dist, in_dist);
+
+
+ // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u)
+ SparsityPattern sparsity;
+ {
+ CompressedSimpleSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+ DoFTools::make_sparsity_pattern (dof, csp, constraints, true);
+ sparsity.copy_from(csp);
+ }
+ SparseMatrix<double> sparse_matrix (sparsity);
+ {
+ QGauss<dim> quadrature_formula(fe_degree+1);
+
+ FEValues<dim> fe_values (dof.get_fe(), quadrature_formula,
+ update_values | update_gradients |
+ update_JxW_values);
+
+ const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof.begin_active(),
+ endc = dof.end();
+ for (; cell!=endc; ++cell)
+ {
+ cell_matrix = 0;
+ fe_values.reinit (cell);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point)
+ +
+ 10. *
+ fe_values.shape_value(i,q_point) *
+ fe_values.shape_value(j,q_point)) *
+ fe_values.JxW(q_point));
+ }
+
+ cell->get_dof_indices(local_dof_indices);
+ constraints.distribute_local_to_global (cell_matrix,
+ local_dof_indices,
+ sparse_matrix);
+ }
+ }
+
+ sparse_matrix.vmult (out, in);
+ out -= out_dist;
+ const double diff_norm = out.linfty_norm() / out_dist.linfty_norm();
+
+ deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_ball (tria);
+ static const HyperBallBoundary<dim> boundary;
+ tria.set_boundary (0, boundary);
+ if (dim < 3 || fe_degree < 2)
+ tria.refine_global(1);
+ tria.begin(tria.n_levels()-1)->set_refine_flag();
+ tria.last()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ typename Triangulation<dim>::active_cell_iterator
+ cell = tria.begin_active (),
+ endc = tria.end();
+ for (; cell!=endc; ++cell)
+ if (cell->center().norm()<1e-8)
+ cell->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+
+ FE_Q<dim> fe (fe_degree);
+ DoFHandler<dim> dof (tria);
+ dof.distribute_dofs(fe);
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints(dof, constraints);
+ VectorTools::interpolate_boundary_values (dof, 0, ZeroFunction<dim>(),
+ constraints);
+ constraints.close();
+
+ do_test<dim, fe_degree, double> (dof, constraints);
+
+ if (dim == 2)
+ {
+ deallog.push("float");
+ deallog.threshold_double(1.e-6);
+ do_test<dim, fe_degree, float> (dof, constraints);
+ deallog.threshold_double(5.e-11);
+ deallog.pop();
+ }
+}
+
+
+
+int main ()
+{
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ deallog << std::setprecision (3);
+
+ {
+ deallog.threshold_double(5.e-11);
+ deallog.push("2d");
+ test<2,1>();
+ test<2,2>();
+ test<2,4>();
+ deallog.pop();
+ deallog.push("3d");
+ test<3,1>();
+ test<3,2>();
+ deallog.pop();
+ }
+}
--- /dev/null
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d:float::Testing FE_Q<2>(1)
+DEAL:2d:float::Norm of difference: 0
+DEAL:2d:float::
+DEAL:2d::Testing FE_Q<2>(2)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d:float::Testing FE_Q<2>(2)
+DEAL:2d:float::Norm of difference: 0
+DEAL:2d:float::
+DEAL:2d::Testing FE_Q<2>(4)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d:float::Testing FE_Q<2>(4)
+DEAL:2d:float::Norm of difference: 0
+DEAL:2d:float::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// evaluating the geometry on the fly as in matrix_vector_15, but using
+// vector-valued finite elements.
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/lac/vector.h>
+
+
+std::ofstream logfile("output");
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+ MatrixFreeTest(const DoFHandler<dim> &dof_handler,
+ const ConstraintMatrix &constraints)
+ :
+ dof_handler (dof_handler),
+ constraints (constraints)
+ {}
+
+ void vmult (VECTOR &dst,
+ const VECTOR &src) const
+ {
+ MappingFEEvaluation<dim,Number> mapping(QGauss<1>(fe_degree+1),
+ update_gradients | update_values |
+ update_JxW_values);
+ VECTOR src_cpy = src;
+ constraints.distribute(src_cpy);
+ FEEvaluation<dim,fe_degree,fe_degree+1,dim,Number> fe_eval(mapping, dof_handler);
+ dst = 0;
+ typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for ( ; cell != endc; ++cell)
+ {
+ mapping.reinit(cell);
+ fe_eval.read_dof_values (src_cpy);
+ fe_eval.evaluate (true, true, false);
+ for (unsigned int q=0; q<fe_eval.n_q_points; ++q)
+ {
+ fe_eval.submit_value (make_vectorized_array<Number>(10.) *
+ fe_eval.get_value(q),q);
+ fe_eval.submit_gradient (fe_eval.get_gradient(q),q);
+ }
+ fe_eval.integrate (true,true);
+ fe_eval.distribute_local_to_global (dst);
+ }
+ constraints.condense(dst);
+ };
+
+private:
+ const DoFHandler<dim> &dof_handler;
+ const ConstraintMatrix &constraints;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof,
+ const ConstraintMatrix &constraints,
+ const unsigned int parallel_option = 0)
+{
+
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+ if (parallel_option > 0)
+ deallog << "Parallel option: " << parallel_option << std::endl;
+ //std::cout << "Number of cells: " << dof.get_tria().n_active_cells() << std::endl;
+ //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+ //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl;
+
+ MatrixFreeTest<dim,fe_degree,number> mf (dof, constraints);
+ Vector<number> in (dof.n_dofs()), out (dof.n_dofs());
+ Vector<number> in_dist (dof.n_dofs());
+ Vector<number> out_dist (in_dist);
+
+ for (unsigned int i=0; i<dof.n_dofs(); ++i)
+ {
+ if (constraints.is_constrained(i))
+ continue;
+ const double entry = Testing::rand()/(double)RAND_MAX;
+ in(i) = entry;
+ in_dist(i) = entry;
+ }
+
+ mf.vmult (out_dist, in_dist);
+
+
+ // assemble sparse matrix with vector-valued form of (\nabla v, \nabla u) +
+ // (v, 10 * u)
+ SparsityPattern sparsity;
+ {
+ CompressedSimpleSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+ DoFTools::make_sparsity_pattern (dof, csp, constraints, true);
+ sparsity.copy_from(csp);
+ }
+ SparseMatrix<double> sparse_matrix (sparsity);
+ {
+ QGauss<dim> quadrature_formula(fe_degree+1);
+
+ FEValues<dim> fe_values (dof.get_fe(), quadrature_formula,
+ update_values | update_gradients |
+ update_JxW_values);
+
+ const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof.begin_active(),
+ endc = dof.end();
+ for (; cell!=endc; ++cell)
+ {
+ cell_matrix = 0;
+ fe_values.reinit (cell);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if (dof.get_fe().system_to_component_index(i).first ==
+ dof.get_fe().system_to_component_index(j).first)
+ cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point)
+ +
+ 10. *
+ fe_values.shape_value(i,q_point) *
+ fe_values.shape_value(j,q_point)) *
+ fe_values.JxW(q_point));
+ }
+
+ cell->get_dof_indices(local_dof_indices);
+ constraints.distribute_local_to_global (cell_matrix,
+ local_dof_indices,
+ sparse_matrix);
+ }
+ }
+
+ sparse_matrix.vmult (out, in);
+ out -= out_dist;
+ const double diff_norm = out.linfty_norm() / out_dist.linfty_norm();
+
+ deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_ball (tria);
+ static const HyperBallBoundary<dim> boundary;
+ tria.set_boundary (0, boundary);
+ if (dim < 3 || fe_degree < 2)
+ tria.refine_global(1);
+ tria.begin(tria.n_levels()-1)->set_refine_flag();
+ tria.last()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ typename Triangulation<dim>::active_cell_iterator
+ cell = tria.begin_active (),
+ endc = tria.end();
+ for (; cell!=endc; ++cell)
+ if (cell->center().norm()<1e-8)
+ cell->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+
+ FESystem<dim> fe(FE_Q<dim>(fe_degree), dim);
+ DoFHandler<dim> dof (tria);
+ dof.distribute_dofs(fe);
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints(dof, constraints);
+ VectorTools::interpolate_boundary_values (dof, 0, ZeroFunction<dim>(dim),
+ constraints);
+ constraints.close();
+
+ do_test<dim, fe_degree, double> (dof, constraints);
+}
+
+
+
+int main ()
+{
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ deallog << std::setprecision (3);
+
+ {
+ deallog.threshold_double(5.e-11);
+ deallog.push("2d");
+ test<2,1>();
+ test<2,2>();
+ test<2,4>();
+ deallog.pop();
+ deallog.push("3d");
+ test<3,1>();
+ test<3,2>();
+ deallog.pop();
+ }
+}
--- /dev/null
+
+DEAL:2d::Testing FESystem<2>[FE_Q<2>(1)^2]
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_Q<2>(2)^2]
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_Q<2>(4)^2]
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Testing FESystem<3>[FE_Q<3>(1)^3]
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FESystem<3>[FE_Q<3>(2)^3]
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
projects / dealii-svn.git / commitdiff