namespace internal
{
DeclException0 (ExcAccessToUninitializedField);
+
+ template <typename FEEval>
+ void do_evaluate (FEEval &, const bool, const bool, const bool);
+ template <typename FEEval>
+ void do_integrate (FEEval &, const bool, const bool);
}
* @author Katharina Kormann and Martin Kronbichler, 2010, 2011
*/
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
class FEEvaluationBase
{
public:
- typedef Tensor<1,n_components,VectorizedArray<Number> > value_type;
- typedef Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > > gradient_type;
+ typedef Number number_type;
+ typedef Tensor<1,n_components_,VectorizedArray<Number> > value_type;
+ typedef Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > gradient_type;
static const unsigned int dimension = dim;
+ static const unsigned int n_components = n_components_;
static const unsigned int dofs_per_cell = dofs_per_cell_;
static const unsigned int n_q_points = n_q_points_;
*/
template <typename VectorType>
void read_dof_values (const std::vector<VectorType> &src,
- const unsigned int first_index=0);
+ const unsigned int first_index=0);
/**
* Reads data from several vectors. Same as
*/
template <typename VectorType>
void read_dof_values (const std::vector<VectorType*> &src,
- const unsigned int first_index=0);
-
- /**
- * For a collection of several vector @p src,
- * read out the values on the degrees of
- * freedom of the current cell for @p
- * n_components (template argument), and store
- * them internally. Similar functionality as
- * the function
- * ConstraintMatrix::read_dof_values. Note
- * that if vectorization is enabled, the DoF
- * values for several cells are set.
- */
- template<typename VectorType>
- void read_dof_values (const VectorType * src_data[]);
+ const unsigned int first_index=0);
/**
* For the vector @p src, read out the values
void read_dof_values_plain (const std::vector<VectorType*> &src,
const unsigned int first_index=0);
- /**
- * For a collection of several vector @p src,
- * read out the values on the degrees of
- * freedom of the current cell for @p
- * n_components (template argument), and store
- * them internally. Similar functionality as
- * the function
- * DoFAccessor::read_dof_values. Note
- * that if vectorization is enabled, the DoF
- * values for several cells are set.
- */
- template<typename VectorType>
- void read_dof_values_plain (const VectorType * src_data[]);
-
/**
* Takes the values stored internally on dof
* values of the current cell and sums them
void distribute_local_to_global (std::vector<VectorType*> &dst,
const unsigned int first_index=0) const;
- /**
- * Takes the values stored internally on dof
- * values of the current cell for a
- * vector-valued problem consisting of @p
- * n_components (template argument) and sums
- * them into the collection of vectors 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.
- */
- template<typename VectorType>
- void distribute_local_to_global (VectorType * dst_data[]) const;
-
/**
* Takes the values stored internally on dof
* values of the current cell and sums them
void set_dof_values (std::vector<VectorType*> &dst,
const unsigned int first_index=0) const;
- /**
- * Takes the values stored internally on dof
- * values of the current cell for a
- * vector-valued problem consisting of @p
- * n_components (template argument) and sums
- * them into the collection of vectors 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.
- */
- template<typename VectorType>
- void set_dof_values (VectorType * dst_data[]) const;
-
//@}
/**
* (n_components == 1) and for the
* vector-valued case (n_components == dim).
*/
- Tensor<1,n_components,Tensor<2,dim,VectorizedArray<Number> > >
+ Tensor<1,n_components_,Tensor<2,dim,VectorizedArray<Number> > >
get_hessian (const unsigned int q_point) const;
/**
//@}
+ /**
+ * @name 4: Access to internal data
+ */
+ //@{
+ /**
+ * Returns a read-only pointer to the first
+ * field of function values on quadrature
+ * points. First come the function values on
+ * all quadrature points for the first
+ * component, then all values for the second
+ * component, and so on. This is related to
+ * the internal data structures used in this
+ * class. The raw data after a call to @p
+ * evaluate only contains unit cell
+ * operations, so possible transformations,
+ * quadrature weights etc. must be applied
+ * manually. In general, it is safer to use
+ * the get_value() function instead, which
+ * does all the transformation internally.
+ */
+ const VectorizedArray<Number> * begin_values () const;
+
+ /**
+ * Returns a read and write pointer to the
+ * first field of function values on
+ * quadrature points. First come the function
+ * values on all quadrature points for the
+ * first component, then all values for the
+ * second component, and so on. This is
+ * related to the internal data structures
+ * used in this class. The raw data after a
+ * call to @p evaluate only contains unit
+ * cell operations, so possible
+ * transformations, quadrature weights
+ * etc. must be applied manually. In general,
+ * it is safer to use the get_value() function
+ * instead, which does all the transformation
+ * internally.
+ */
+ VectorizedArray<Number> * begin_values ();
+
+ /**
+ * Returns a read-only pointer to the first
+ * field of function gradients on quadrature
+ * points. First comes the x-component of the
+ * gradient for the first component on all
+ * quadrature points, then the y-component,
+ * and so on. Next comes the x-component of
+ * the second component, and so on. This is
+ * related to the internal data structures
+ * used in this class. The raw data after a
+ * call to @p evaluate only contains unit
+ * cell operations, so possible
+ * transformations, quadrature weights
+ * etc. must be applied manually. In general,
+ * it is safer to use the get_gradient() function
+ * instead, which does all the transformation
+ * internally.
+ */
+ const VectorizedArray<Number> * begin_gradients () const;
+
+ /**
+ * Returns a read and write pointer to the
+ * first field of function gradients on
+ * quadrature points. First comes the
+ * x-component of the gradient for the first
+ * component on all quadrature points, then
+ * the y-component, and so on. Next comes the
+ * x-component of the second component, and so
+ * on. This is related to the internal data
+ * structures used in this class. The raw data
+ * after a call to @p evaluate only
+ * contains unit cell operations, so possible
+ * transformations, quadrature weights
+ * etc. must be applied manually. In general,
+ * it is safer to use the get_gradient()
+ * function instead, which does all the
+ * transformation internally.
+ */
+ VectorizedArray<Number> * begin_gradients ();
+
+ /**
+ * Returns a read-only pointer to the first
+ * field of function hessians on quadrature
+ * points. First comes the xx-component of the
+ * hessian for the first component on all
+ * quadrature points, then the yy-component,
+ * zz-component in (3D), then the
+ * xy-component, and so on. Next comes the
+ * xx-component of the second component, and
+ * so on. This is related to the internal data
+ * structures used in this class. The raw data
+ * after a call to @p evaluate only
+ * contains unit cell operations, so possible
+ * transformations, quadrature weights
+ * etc. must be applied manually. In general,
+ * it is safer to use the get_laplacian() or
+ * get_hessian() functions instead, which does
+ * all the transformation internally.
+ */
+ const VectorizedArray<Number> * begin_hessians () const;
+
+ /**
+ * Returns a read and write pointer to the
+ * first field of function hessians on
+ * quadrature points. First comes the
+ * xx-component of the hessian for the first
+ * component on all quadrature points, then
+ * the yy-component, zz-component in (3D),
+ * then the xy-component, and so on. Next
+ * comes the xx-component of the second
+ * component, and so on. This is related to
+ * the internal data structures used in this
+ * class. The raw data after a call to @p
+ * evaluate only contains unit cell
+ * operations, so possible transformations,
+ * quadrature weights etc. must be applied
+ * manually. In general, it is safer to use
+ * the get_laplacian() or get_hessian()
+ * functions instead, which does all the
+ * transformation internally.
+ */
+ VectorizedArray<Number> * begin_hessians ();
+
+ //@}
+
protected:
/**
const unsigned int fe_no = 0,
const unsigned int quad_no = 0);
+ /**
+ * A unified function to read from and write
+ * into vectors based on the given template
+ * operation. It can perform the operation for
+ * @p read_dof_values, @p
+ * distribute_local_to_global, and @p
+ * set_dof_values. It performs the operation
+ * for several vectors at a time.
+ */
+ template<typename VectorType, typename VectorOperation>
+ void read_write_operation (const VectorOperation &operation,
+ VectorType *vectors[]) const;
+
+ /**
+ * For a collection of several vector @p src,
+ * read out the values on the degrees of
+ * freedom of the current cell for @p
+ * n_components (template argument), and store
+ * them internally. Similar functionality as
+ * the function
+ * DoFAccessor::read_dof_values. Note
+ * that if vectorization is enabled, the DoF
+ * values for several cells are set.
+ */
+ template<typename VectorType>
+ void read_dof_values_plain (const VectorType * src_data[]);
+
/**
* Internal data fields that store the
* values. Since all array lengths are known
*/
const internal::MatrixFreeFunctions::ShapeInfo<Number> &data;
- /**
- * After a call to reinit(), stores the number
- * of the cell we are currently working with.
- */
- unsigned int cell;
-
- /**
- * Stores the type of the cell we are
- * currently working with after a call to
- * reinit(). Valid values are @p cartesian, @p
- * affine and @p general, which have different
- * implications on how the Jacobian
- * transformations are stored internally in
- * MappingInfo.
- */
- internal::MatrixFreeFunctions::CellType cell_type;
-
- /**
- * The stride to access the correct data in
- * MappingInfo.
- */
- unsigned int cell_data_number;
-
- /**
- * Stores whether the present cell chunk used
- * in vectorization is not completely filled
- * up with physical cells. E.g. if
- * vectorization dictates that four cells
- * should be worked with but only three
- * physical cells are left, this flag will be
- * set to true, otherwise to false. Mainly
- * used for internal checking when reading
- * from vectors or writing to vectors.
- */
- bool at_irregular_cell;
-
- /**
- * If the present cell chunk for vectorization
- * is not completely filled up with data, this
- * field stores how many physical cells are
- * underlying. Is between 1 and
- * VectorizedArray<Number>::n_array_elements-1
- * (inclusive).
- */
- unsigned int n_irreg_components_filled;
-
/**
* A pointer to the Cartesian Jacobian
* information of the present cell. Only set
*/
const Tensor<1,(dim>1?dim*(dim-1)/2:1),Tensor<1,dim,VectorizedArray<Number> > > * jacobian_grad_upper;
+ /**
+ * After a call to reinit(), stores the number
+ * of the cell we are currently working with.
+ */
+ unsigned int cell;
+
+ /**
+ * Stores the type of the cell we are
+ * currently working with after a call to
+ * reinit(). Valid values are @p cartesian, @p
+ * affine and @p general, which have different
+ * implications on how the Jacobian
+ * transformations are stored internally in
+ * MappingInfo.
+ */
+ internal::MatrixFreeFunctions::CellType cell_type;
+
+ /**
+ * The stride to access the correct data in
+ * MappingInfo.
+ */
+ unsigned int cell_data_number;
+
+ /**
+ * If the present cell chunk for vectorization
+ * is not completely filled up with data, this
+ * field stores how many physical cells are
+ * underlying. Is between 1 and
+ * VectorizedArray<Number>::n_array_elements-1
+ * (inclusive).
+ */
+ unsigned int n_irreg_components_filled;
+
+ /**
+ * Stores whether the present cell chunk used
+ * in vectorization is not completely filled
+ * up with physical cells. E.g. if
+ * vectorization dictates that four cells
+ * should be worked with but only three
+ * physical cells are left, this flag will be
+ * set to true, otherwise to false. Mainly
+ * used for internal checking when reading
+ * from vectors or writing to vectors.
+ */
+ bool at_irregular_cell;
+
/**
* Debug information to track whether dof
* values have been initialized before
* accessed. Used to control exceptions when
* uninitialized data is used.
*/
- bool dof_values_initialized;
+ bool dof_values_initialized;
/**
* Debug information to track whether values
* before accessed. Used to control exceptions
* when uninitialized data is used.
*/
- bool values_quad_initialized;
+ bool values_quad_initialized;
/**
* Debug information to track whether
* control exceptions when uninitialized data
* is used.
*/
- bool gradients_quad_initialized;
+ bool gradients_quad_initialized;
/**
* Debug information to track whether
* control exceptions when uninitialized data
* is used.
*/
- bool hessians_quad_initialized;
+ bool hessians_quad_initialized;
/**
* Debug information to track whether values
* actually stared. Used to control exceptions
* when uninitialized data is used.
*/
- bool values_quad_submitted;
+ bool values_quad_submitted;
/**
* Debug information to track whether
* control exceptions when uninitialized data
* is used.
*/
- bool gradients_quad_submitted;
+ bool gradients_quad_submitted;
};
* @author Katharina Kormann and Martin Kronbichler, 2010, 2011
*/
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
class FEEvaluationAccess :
- public FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>
+ public FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
{
public:
- typedef Tensor<1,n_components,VectorizedArray<Number> > value_type;
- typedef Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > > gradient_type;
+ typedef Number number_type;
+ typedef Tensor<1,n_components_,VectorizedArray<Number> > value_type;
+ typedef Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > gradient_type;
static const unsigned int dimension = dim;
+ static const unsigned int n_components = n_components_;
static const unsigned int dofs_per_cell = dofs_per_cell_;
static const unsigned int n_q_points = n_q_points_;
- typedef FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,
+ typedef FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,
Number> BaseClass;
protected:
public FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,1,Number>
{
public:
+ typedef Number number_type;
typedef VectorizedArray<Number> value_type;
typedef Tensor<1,dim,VectorizedArray<Number> > gradient_type;
static const unsigned int dimension = dim;
public FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,dim,Number>
{
public:
- typedef Tensor<1,dim,VectorizedArray<Number> > value_type;
- typedef Tensor<2,dim,VectorizedArray<Number> > gradient_type;
+ typedef Number number_type;
+ typedef Tensor<1,dim,VectorizedArray<Number> > value_type;
+ typedef Tensor<2,dim,VectorizedArray<Number> > gradient_type;
static const unsigned int dimension = dim;
+ static const unsigned int n_components = dim;
static const unsigned int dofs_per_cell = dofs_per_cell_;
static const unsigned int n_q_points = n_q_points_;
typedef FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,dim,Number> BaseClass;
*
* @author Katharina Kormann and Martin Kronbichler, 2010, 2011
*/
-template <int dim, int fe_degree, int n_q_points_1d=fe_degree+1,
- int n_components=1, typename Number=double >
+template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
+ int n_components_ = 1, typename Number = double >
class FEEvaluationGeneral :
public FEEvaluationAccess<dim,
Utilities::fixed_int_power<fe_degree+1,dim>::value,
Utilities::fixed_int_power<n_q_points_1d,dim>::value,
- n_components,Number>
+ n_components_,Number>
{
public:
typedef FEEvaluationAccess<dim,
Utilities::fixed_int_power<fe_degree+1,dim>::value,
Utilities::fixed_int_power<n_q_points_1d,dim>::value,
- n_components, Number> BaseClass;
+ n_components_, Number> BaseClass;
+ typedef Number number_type;
typedef typename BaseClass::value_type value_type;
typedef typename BaseClass::gradient_type gradient_type;
+ static const unsigned int dimension = dim;
+ static const unsigned int n_components = n_components_;
static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell;
static const unsigned int n_q_points = BaseClass::n_q_points;
* @p get_gradient() or @p get_laplacian
* return useful information.
*/
- void evaluate (bool evaluate_val, bool evaluate_grad,
- bool evaluate_hess=false);
+ void evaluate (const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_hess = false);
/**
* This function takes the values and/or
* integrate_grad are used to enable/disable
* some of values or gradients.
*/
- void integrate (bool integrate_val, bool integrate_grad);
+ void integrate (const bool integrate_val,
+ const bool integrate_grad);
/**
* Returns the q-th quadrature point stored in
protected:
/**
- * Internal function that applies the shape
- * function data of the tensor product in a
- * given coordinate direction (first template
+ * Internal function that applies the function
+ * values of the tensor product in a given
+ * coordinate direction (first template
+ * argument), from polynomials to values on
+ * quadrature points (second flag set to true)
+ * or in an integration loop from values on
+ * quadrature points to values tested by
+ * different test function (second flag set to
+ * false), and if the result is to be added to
+ * previous content in the data fields or
+ * not.
+ */
+ template <int direction, bool dof_to_quad, bool add>
+ void apply_values (const VectorizedArray<Number> in [],
+ VectorizedArray<Number> out []);
+
+ /**
+ * Internal function that applies the gradient
+ * operation of the tensor product in a given
+ * coordinate direction (first template
* argument), from polynomials to values on
* quadrature points (second flag set to true)
* or in an integration loop from values on
* quadrature points to values tested by
* different test function (second flag set to
* false), and if the result is to be added to
- * some previous results or not.
+ * previous content in the data fields or
+ * not.
+ */
+ template <int direction, bool dof_to_quad, bool add>
+ void apply_gradients (const VectorizedArray<Number> in [],
+ VectorizedArray<Number> out []);
+
+ /**
+ * Internal function that applies the second
+ * derivative operation (Hessian) of the
+ * tensor product in a given coordinate
+ * direction (first template argument), from
+ * polynomials to values on quadrature points
+ * (second flag set to true) or in an
+ * integration loop from values on quadrature
+ * points to values tested by different test
+ * function (second flag set to false), and if
+ * the result is to be added to previous
+ * content in the data fields or not.
*/
template <int direction, bool dof_to_quad, bool add>
- void apply_tensor_prod (const VectorizedArray<Number> * shape_data,
- const VectorizedArray<Number> in [],
- VectorizedArray<Number> out []);
+ void apply_hessians (const VectorizedArray<Number> in [],
+ VectorizedArray<Number> out []);
+
+ /**
+ * Friend declaration.
+ */
+ template <typename FEEval> friend void
+ internal::do_evaluate (FEEval &, const bool, const bool, const bool);
+ template <typename FEEval> friend void
+ internal::do_integrate (FEEval &, const bool, const bool);
};
*
* @author Katharina Kormann and Martin Kronbichler, 2010, 2011
*/
-template <int dim, int fe_degree, int n_q_points_1d=fe_degree+1,
- int n_components=1, typename Number=double >
+template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
+ int n_components_ = 1, typename Number = double >
class FEEvaluation :
- public FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components,Number>
+ public FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
{
public:
- typedef FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components,Number> BaseClass;
+ typedef FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number> BaseClass;
+ typedef Number number_type;
typedef typename BaseClass::value_type value_type;
typedef typename BaseClass::gradient_type gradient_type;
+ static const unsigned int dimension = dim;
+ static const unsigned int n_components = n_components_;
static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell;
static const unsigned int n_q_points = BaseClass::n_q_points;
* (unless these values have been set
* manually).
*/
- void evaluate (bool evaluate_val, bool evaluate_grad,
- bool evaluate_hess=false);
+ void evaluate (const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_hess = false);
/**
* This function takes the values and/or
* integrate_grad are used to enable/disable
* some of values or gradients.
*/
- void integrate (bool integrate_val, bool integrate_grad);
+ void integrate (const bool integrate_val,
+ const bool integrate_grad);
protected:
+
/**
* Internal function that applies the function
* values of the tensor product in a given
template <int direction, bool dof_to_quad, bool add>
void apply_hessians (const VectorizedArray<Number> in [],
VectorizedArray<Number> out []);
+
+ /**
+ * Friend declarations.
+ */
+ template <typename FEEval> friend void
+ internal::do_evaluate (FEEval &, const bool, const bool, const bool);
+ template <typename FEEval> friend void
+ internal::do_integrate (FEEval &, const bool, const bool);
};
*
* @author Katharina Kormann and Martin Kronbichler, 2010, 2011
*/
-template <int dim, int fe_degree, int n_components=1, typename Number=double >
+template <int dim, int fe_degree, int n_components_ = 1, typename Number = double >
class FEEvaluationGL :
- public FEEvaluation<dim,fe_degree,fe_degree+1,n_components,Number>
+ public FEEvaluation<dim,fe_degree,fe_degree+1,n_components_,Number>
{
public:
- typedef FEEvaluation<dim,fe_degree,fe_degree+1,n_components,Number> BaseClass;
+ typedef FEEvaluation<dim,fe_degree,fe_degree+1,n_components_,Number> BaseClass;
+ typedef Number number_type;
typedef typename BaseClass::value_type value_type;
typedef typename BaseClass::gradient_type gradient_type;
+ static const unsigned int dimension = dim;
+ static const unsigned int n_components = n_components_;
static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell;
static const unsigned int n_q_points = BaseClass::n_q_points;
* (unless these values have been set
* manually).
*/
- void evaluate (bool evaluate_val, bool evaluate_grad,
- bool evaluate_lapl=false);
+ void evaluate (const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_lapl = false);
/**
* This function takes the values and/or
* integrate_grad are used to enable/disable
* some of values or gradients.
*/
- void integrate (bool integrate_val, bool integrate_grad);
+ void integrate (const bool integrate_val,
+ const bool integrate_grad);
protected:
/**
#ifndef DOXYGEN
-/*----------------------- FEEvaluationBase -------------------------------*/
+/*----------------------- FEEvaluationBase ----------------------------------*/
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-FEEvaluationBase (const MatrixFree<dim,Number> &data_in,
- const unsigned int fe_no_in,
- const unsigned int quad_no_in)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::FEEvaluationBase (const MatrixFree<dim,Number> &data_in,
+ const unsigned int fe_no_in,
+ const unsigned int quad_no_in)
:
quad_no (quad_no_in),
n_fe_components (data_in.get_dof_info(fe_no_in).n_components),
data (data_in.get_shape_info
(fe_no_in, quad_no_in, active_fe_index,
active_quad_index)),
- cell (numbers::invalid_unsigned_int),
- cell_type (internal::MatrixFreeFunctions::undefined),
- cell_data_number (0),
- at_irregular_cell (false),
- n_irreg_components_filled (0),
cartesian_data (0),
jacobian (0),
J_value (0),
quadrature_weights[active_quad_index].begin()),
quadrature_points (0),
jacobian_grad (0),
- jacobian_grad_upper(0)
+ jacobian_grad_upper(0),
+ cell (numbers::invalid_unsigned_int),
+ cell_type (internal::MatrixFreeFunctions::undefined),
+ cell_data_number (0),
+ n_irreg_components_filled (0),
+ at_irregular_cell (false)
{
Assert (matrix_info.indices_initialized() == true,
ExcNotInitialized());
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-reinit (const unsigned int cell_in)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::reinit (const unsigned int cell_in)
{
AssertIndexRange (cell_in, dof_info.row_starts.size()-1);
AssertDimension (((dof_info.cell_active_fe_index.size() > 0) ?
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
unsigned int
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_cell_data_number () const
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_cell_data_number () const
{
Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
return cell_data_number;
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
internal::MatrixFreeFunctions::CellType
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_cell_type () const
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_cell_type () const
{
Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
return cell_type;
"Use MatrixFree::initialize_dof_vector to get a "
"compatible vector."));
}
-}
+ // A class to use the same code to read from
+ // and write to vector
+ template <typename Number>
+ struct VectorReader
+ {
+ template <typename VectorType>
+ void process_dof (const unsigned int index,
+ VectorType &vec,
+ Number &res) const
+ {
+ res = vector_access (vec, index);
+ }
+ void pre_constraints (const Number &,
+ Number &res) const
+ {
+ res = Number();
+ }
-template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
-template<typename VectorType>
-inline
-void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values (const VectorType &src)
-{
- AssertDimension (n_components, n_fe_components);
- // only need one component, but to avoid
- // compiler warnings, use n_components copies
- // here (but these will not be used)
- const VectorType * src_data[n_components];
- for (unsigned int d=0; d<n_components; ++d)
- src_data[d] = &src;
- read_dof_values (src_data);
-}
+ template <typename VectorType>
+ void process_constraint (const unsigned int index,
+ const Number weight,
+ VectorType &vec,
+ Number &res) const
+ {
+ res += weight * vector_access (vec, index);
+ }
+ void post_constraints (const Number &sum,
+ Number &write_pos) const
+ {
+ write_pos = sum;
+ }
+ void process_empty (Number &res) const
+ {
+ res = Number();
+ }
+ };
-template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
-template<typename VectorType>
-inline
-void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values (const std::vector<VectorType> &src,
- const unsigned int first_index)
-{
- AssertIndexRange (first_index, src.size());
- Assert (n_fe_components == 1, ExcNotImplemented());
- Assert ((n_fe_components == 1 ?
- (first_index+n_components <= src.size()) : true),
- ExcIndexRange (first_index + n_components, 0, src.size()));
- const VectorType * src_data [n_components];
- for (unsigned int comp=0; comp<n_components; ++comp)
- src_data[comp] = &src[comp+first_index];
- read_dof_values (src_data);
-}
+ // A class to use the same code to read from
+ // and write to vector
+ template <typename Number>
+ struct VectorDistributorLocalToGlobal
+ {
+ template <typename VectorType>
+ void process_dof (const unsigned int index,
+ VectorType &vec,
+ Number &res) const
+ {
+ vector_access (vec, index) += res;
+ }
+ void pre_constraints (const Number &input,
+ Number &res) const
+ {
+ res = input;
+ }
+ template <typename VectorType>
+ void process_constraint (const unsigned int index,
+ const Number weight,
+ VectorType &vec,
+ Number &res) const
+ {
+ vector_access (vec, index) += weight * res;
+ }
+
+ void post_constraints (const Number &,
+ Number &) const
+ {
+ }
+
+ void process_empty (Number &) const
+ {
+ }
+ };
+
+
+ // A class to use the same code to read from
+ // and write to vector
+ template <typename Number>
+ struct VectorSetter
+ {
+ template <typename VectorType>
+ void process_dof (const unsigned int index,
+ VectorType &vec,
+ Number &res) const
+ {
+ vector_access (vec, index) = res;
+ }
+
+ void pre_constraints (const Number &,
+ Number &) const
+ {
+ }
+
+ template <typename VectorType>
+ void process_constraint (const unsigned int,
+ const Number,
+ VectorType&,
+ Number &) const
+ {
+ }
+
+ void post_constraints (const Number &,
+ Number &) const
+ {
+ }
+
+ void process_empty (Number &) const
+ {
+ }
+ };
-template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
-template<typename VectorType>
-inline
-void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values (const std::vector<VectorType*> &src,
- const unsigned int first_index)
-{
- AssertIndexRange (first_index, src.size());
- Assert (n_fe_components == 1, ExcNotImplemented());
- Assert ((n_fe_components == 1 ?
- (first_index+n_components <= src.size()) : true),
- ExcIndexRange (first_index + n_components, 0, src.size()));
- const VectorType * src_data [n_components];
- for (unsigned int comp=0; comp<n_components; ++comp)
- src_data[comp] = src[comp+first_index];
- read_dof_values (src_data);
}
+
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
-template<typename VectorType>
+ int n_components_, typename Number>
+template<typename VectorType, typename VectorOperation>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values (const VectorType * src[])
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_write_operation (const VectorOperation &operation,
+ VectorType *src[]) const
{
+ // This functions processes all the functions
+ // read_dof_values,
+ // distribute_local_to_global, and
+ // set_dof_values with the same code. The
+ // distinction between these three cases is
+ // made by the input VectorOperation that
+ // either reads values from a vector and puts
+ // the data into the local data field or write
+ // local data into the vector. Certain
+ // operations are no-ops for the given use
+ // case.
+
Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
// loop over all local dofs. ind_local holds
VectorizedArray<Number>::n_array_elements * dofs_per_cell;
for (unsigned int comp=0; comp<n_components; ++comp)
internal::check_vector_compatibility (*src[comp], dof_info);
- Number * local_src_number [n_components];
+ Number * local_data [n_components];
for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp] = reinterpret_cast<Number*>(values_dofs[comp]);
+ local_data[comp] =
+ const_cast<Number*>(reinterpret_cast<const Number*>(values_dofs[comp]));
// standard case where there are sufficiently
// many cells to fill all vectors
// run through values up to next constraint
for (unsigned int j=0; j<indicators->first; ++j)
for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][ind_local+j] =
- internal::vector_access (*src[comp], dof_indices[j]);
+ operation.process_dof (dof_indices[j], *src[comp],
+ local_data[comp][ind_local+j]);
+
ind_local += indicators->first;
dof_indices += indicators->first;
// according to constraints
Number value [n_components];
for (unsigned int comp=0; comp<n_components; ++comp)
- value[comp] = 0;
+ operation.pre_constraints (local_data[comp][ind_local],
+ value[comp]);
+
const Number * data_val =
matrix_info.constraint_pool_begin(indicators->second);
const Number * end_pool =
matrix_info.constraint_pool_end(indicators->second);
for ( ; data_val != end_pool; ++data_val, ++dof_indices)
for (unsigned int comp=0; comp<n_components; ++comp)
- value[comp] +=
- (internal::vector_access (*src[comp], *dof_indices) *
- (*data_val));
+ operation.process_constraint (*dof_indices, *data_val,
+ *src[comp], value[comp]);
for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][ind_local] = value[comp];
+ operation.post_constraints (value[comp],
+ local_data[comp][ind_local]);
+
ind_local++;
}
for(; ind_local < n_local_dofs; ++dof_indices, ++ind_local)
{
for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][ind_local] =
- internal::vector_access (*src[comp], *dof_indices);
+ operation.process_dof (*dof_indices, *src[comp],
+ local_data[comp][ind_local]);
}
}
else
static_cast<int>(n_local_dofs));
for (unsigned int j=0; j<n_local_dofs; ++j)
for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][j] =
- internal::vector_access (*src[comp], dof_indices[j]);
+ operation.process_dof (dof_indices[j], *src[comp],
+ local_data[comp][j]);
}
}
// the global vector, src, to the local one,
// local_src.
for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][ind_local] =
- internal::vector_access (*src[comp], dof_indices[j]);
+ operation.process_dof (dof_indices[j], *src[comp],
+ local_data[comp][ind_local]);
// here we jump over all the components that
// are artificial
>= n_irreg_components_filled)
{
for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][ind_local] = 0.;
+ operation.process_empty (local_data[comp][ind_local]);
++ind_local;
}
}
// according to constraint
Number value [n_components];
for (unsigned int comp=0; comp<n_components; ++comp)
- value[comp] = 0.;
+ operation.pre_constraints (local_data[comp][ind_local],
+ value[comp]);
+
const Number * data_val =
matrix_info.constraint_pool_begin(indicators->second);
const Number * end_pool =
matrix_info.constraint_pool_end(indicators->second);
+
for ( ; data_val != end_pool; ++data_val, ++dof_indices)
for (unsigned int comp=0; comp<n_components; ++comp)
- value[comp] +=
- internal::vector_access (*src[comp], *dof_indices) * (*data_val);
+ operation.process_constraint (*dof_indices, *data_val,
+ *src[comp], value[comp]);
+
for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][ind_local] = value[comp];
+ operation.post_constraints (value[comp],
+ local_data[comp][ind_local]);
ind_local++;
while (ind_local % VectorizedArray<Number>::n_array_elements
>= n_irreg_components_filled)
{
for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][ind_local] = 0.;
+ operation.process_empty (local_data[comp][ind_local]);
++ind_local;
}
}
// the global vector, src, to the local one,
// local_dst.
for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][ind_local] =
- internal::vector_access (*src[comp], *dof_indices);
+ operation.process_dof (*dof_indices, *src[comp],
+ local_data[comp][ind_local]);
++ind_local;
while (ind_local % VectorizedArray<Number>::n_array_elements
>= n_irreg_components_filled)
{
for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][ind_local] = 0.;
+ operation.process_empty(local_data[comp][ind_local]);
++ind_local;
}
}
// out the indices.
{
internal::check_vector_compatibility (*src[0], dof_info);
- Assert (n_fe_components == n_components, ExcNotImplemented());
+ Assert (n_fe_components == n_components_, ExcNotImplemented());
const unsigned int n_local_dofs =
dofs_per_cell*VectorizedArray<Number>::n_array_elements * n_components;
- Number * local_src_number = reinterpret_cast<Number*>(values_dofs[0]);
+ Number * local_data =
+ const_cast<Number*>(reinterpret_cast<const Number*>(values_dofs[0]));
if (at_irregular_cell == false)
{
// check whether there is any constraint on
{
// run through values up to next constraint
for (unsigned int j=0; j<indicators->first; ++j)
- local_src_number[ind_local+j] =
- internal::vector_access (*src[0], dof_indices[j]);
+ operation.process_dof (dof_indices[j], *src[0],
+ local_data[ind_local+j]);
ind_local += indicators->first;
dof_indices += indicators->first;
// constrained case: build the local value as
// a linear combination of the global value
// according to constraints
- Number value = 0;
+ Number value;
+ operation.pre_constraints (local_data[ind_local], value);
+
const Number * data_val =
matrix_info.constraint_pool_begin(indicators->second);
const Number * end_pool =
matrix_info.constraint_pool_end(indicators->second);
+
for ( ; data_val != end_pool; ++data_val, ++dof_indices)
- value +=
- (internal::vector_access (*src[0], *dof_indices) *
- (*data_val));
+ operation.process_constraint (*dof_indices, *data_val,
+ *src[0], value);
- local_src_number[ind_local] = value;
+ operation.post_constraints (value, local_data[ind_local]);
ind_local++;
}
// get the dof values past the last
// constraint
for(; ind_local<n_local_dofs; ++dof_indices, ++ind_local)
- local_src_number[ind_local] =
- internal::vector_access (*src[0], *dof_indices);
+ operation.process_dof (*dof_indices, *src[0],
+ local_data[ind_local]);
Assert (dof_indices == dof_info.end_indices(cell),
ExcInternalError());
}
AssertDimension (dof_info.end_indices(cell)-dof_indices,
static_cast<int>(n_local_dofs));
for (unsigned int j=0; j<n_local_dofs; ++j)
- local_src_number[j] =
- internal::vector_access (*src[0], dof_indices[j]);
+ operation.process_dof (dof_indices[j], *src[0],
+ local_data[j]);
}
}
// non-constrained case: copy the data from
// the global vector, src, to the local one,
// local_src.
- local_src_number[ind_local] =
- internal::vector_access (*src[0], dof_indices[j]);
+ operation.process_dof (dof_indices[j], *src[0],
+ local_data[ind_local]);
// here we jump over all the components that
// are artificial
while (ind_local % VectorizedArray<Number>::n_array_elements
>= n_irreg_components_filled)
{
- local_src_number[ind_local] = 0.;
+ operation.process_empty (local_data[ind_local]);
++ind_local;
}
}
// constrained case: build the local value as
// a linear combination of the global value
// according to constraint
- Number value = 0;
+ Number value;
+ operation.pre_constraints (local_data[ind_local], value);
+
const Number * data_val =
matrix_info.constraint_pool_begin(indicators->second);
const Number * end_pool =
matrix_info.constraint_pool_end(indicators->second);
+
for ( ; data_val != end_pool; ++data_val, ++dof_indices)
- value +=
- internal::vector_access (*src[0], *dof_indices) * (*data_val);
- local_src_number[ind_local] = value;
+ operation.process_constraint (*dof_indices, *data_val,
+ *src[0], value);
+
+ operation.post_constraints (value, local_data[ind_local]);
ind_local++;
while (ind_local % VectorizedArray<Number>::n_array_elements
>= n_irreg_components_filled)
{
- local_src_number[ind_local] = 0.;
+ operation.process_empty (local_data[ind_local]);
++ind_local;
}
}
// non-constrained case: copy the data from
// the global vector, src, to the local one,
// local_dst.
- local_src_number[ind_local] =
- internal::vector_access (*src[0], *dof_indices);
+ operation.process_dof (*dof_indices, *src[0],
+ local_data[ind_local]);
++ind_local;
while (ind_local % VectorizedArray<Number>::n_array_elements
>= n_irreg_components_filled)
{
- local_src_number[ind_local] = 0.;
+ operation.process_empty (local_data[ind_local]);
++ind_local;
}
}
}
}
+}
+
+
+
+template <int dim, int dofs_per_cell_, int n_q_points_,
+ int n_components_, typename Number>
+template<typename VectorType>
+inline
+void
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values (const VectorType &src)
+{
+ AssertDimension (n_components_, n_fe_components);
+
+ // only need one component, but to silent
+ // compiler warnings, use n_components copies
+ // here (but these will not be used)
+ VectorType * src_data[n_components];
+ for (unsigned int d=0; d<n_components; ++d)
+ src_data[d] = const_cast<VectorType*>(&src);
+
+ internal::VectorReader<Number> reader;
+ read_write_operation (reader, src_data);
#ifdef DEBUG
dof_values_initialized = true;
+template <int dim, int dofs_per_cell_, int n_q_points_,
+ int n_components_, typename Number>
+template<typename VectorType>
+inline
+void
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values (const std::vector<VectorType> &src,
+ const unsigned int first_index)
+{
+ AssertIndexRange (first_index, src.size());
+ Assert (n_fe_components == 1, ExcNotImplemented());
+ Assert ((n_fe_components == 1 ?
+ (first_index+n_components <= src.size()) : true),
+ ExcIndexRange (first_index + n_components_, 0, src.size()));
+
+ VectorType * src_data [n_components];
+ for (unsigned int comp=0; comp<n_components; ++comp)
+ src_data[comp] = const_cast<VectorType*>(&src[comp+first_index]);
+
+ internal::VectorReader<Number> reader;
+ read_write_operation (reader, src_data);
+
+#ifdef DEBUG
+ dof_values_initialized = true;
+#endif
+}
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values_plain (const VectorType &src)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values (const std::vector<VectorType*> &src,
+ const unsigned int first_index)
+{
+ AssertIndexRange (first_index, src.size());
+ Assert (n_fe_components == 1, ExcNotImplemented());
+ Assert ((n_fe_components == 1 ?
+ (first_index+n_components <= src.size()) : true),
+ ExcIndexRange (first_index + n_components_, 0, src.size()));
+
+ const VectorType * src_data [n_components];
+ for (unsigned int comp=0; comp<n_components; ++comp)
+ src_data[comp] = const_cast<VectorType*>(src[comp+first_index]);
+
+ internal::VectorReader<Number> reader;
+ read_write_operation (reader, src_data);
+
+#ifdef DEBUG
+ dof_values_initialized = true;
+#endif
+}
+
+
+
+template <int dim, int dofs_per_cell_, int n_q_points_,
+ int n_components_, typename Number>
+template<typename VectorType>
+inline
+void
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values_plain (const VectorType &src)
{
- AssertDimension (n_components, n_fe_components);
+ AssertDimension (n_components_, n_fe_components);
// only need one component, but to avoid
// compiler warnings, use n_components copies
// here (but these will not be used)
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values_plain (const std::vector<VectorType> &src,
- const unsigned int first_index)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values_plain (const std::vector<VectorType> &src,
+ const unsigned int first_index)
{
AssertIndexRange (first_index, src.size());
Assert (n_fe_components == 1, ExcNotImplemented());
Assert ((n_fe_components == 1 ?
(first_index+n_components <= src.size()) : true),
- ExcIndexRange (first_index + n_components, 0, src.size()));
+ ExcIndexRange (first_index + n_components_, 0, src.size()));
const VectorType * src_data [n_components];
for (unsigned int comp=0; comp<n_components; ++comp)
src_data[comp] = &src[comp+first_index];
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values_plain (const std::vector<VectorType*> &src,
- const unsigned int first_index)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values_plain (const std::vector<VectorType*> &src,
+ const unsigned int first_index)
{
AssertIndexRange (first_index, src.size());
Assert (n_fe_components == 1, ExcNotImplemented());
Assert ((n_fe_components == 1 ?
(first_index+n_components <= src.size()) : true),
- ExcIndexRange (first_index + n_components, 0, src.size()));
+ ExcIndexRange (first_index + n_components_, 0, src.size()));
const VectorType * src_data [n_components];
for (unsigned int comp=0; comp<n_components; ++comp)
src_data[comp] = src[comp+first_index];
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-read_dof_values_plain (const VectorType * src[])
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::distribute_local_to_global (VectorType &dst) const
{
- Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
- Assert (dof_info.store_plain_indices == true, ExcNotInitialized());
+ AssertDimension (n_components_, n_fe_components);
+ Assert (dof_values_initialized==true,
+ internal::ExcAccessToUninitializedField());
- // loop over all local dofs. ind_local holds
- // local number on cell, index iterates over
- // the elements of index_local_to_global and
- // dof_indices points to the global indices
- // stored in index_local_to_global
- const unsigned int * dof_indices = dof_info.begin_indices_plain(cell);
+ // only need one component, but to avoid
+ // compiler warnings, use n_components copies
+ // here (but these will not be used)
+ VectorType * dst_data [n_components];
+ for (unsigned int d=0; d<n_components; ++d)
+ dst_data[d] = &dst;
- // scalar case (or case when all components
- // have the same degrees of freedom and sit on
- // a different vector each)
- if (n_fe_components == 1)
- {
- const unsigned int n_local_dofs =
- VectorizedArray<Number>::n_array_elements * dofs_per_cell;
- for (unsigned int comp=0; comp<n_components; ++comp)
- internal::check_vector_compatibility (*src[comp], dof_info);
- Number * local_src_number [n_components];
- for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp] = reinterpret_cast<Number*>(values_dofs[comp]);
+ internal::VectorDistributorLocalToGlobal<Number> distributor;
+ read_write_operation (distributor, dst_data);
+}
- // standard case where there are sufficiently
- // many cells to fill all vectors
- if (at_irregular_cell == false)
- {
- for (unsigned int j=0; j<n_local_dofs; ++j)
- for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][j] =
- internal::vector_access (*src[comp], dof_indices[j]);
- }
- // non-standard case: need to fill in zeros
- // for those components that are not present
- // (a bit more expensive), but there is not
- // more than one such cell
- else
- {
- Assert (n_irreg_components_filled > 0, ExcInternalError());
- for(unsigned int ind_local=0; ind_local<n_local_dofs;
- ++dof_indices)
- {
- // non-constrained case: copy the data from
- // the global vector, src, to the local one,
- // local_dst.
- for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][ind_local] =
- internal::vector_access (*src[comp], *dof_indices);
- ++ind_local;
- while (ind_local % VectorizedArray<Number>::n_array_elements >= n_irreg_components_filled)
- {
- for (unsigned int comp=0; comp<n_components; ++comp)
- local_src_number[comp][ind_local] = 0.;
- ++ind_local;
- }
- }
- }
- }
- else
- // case with vector-valued finite elements
- // where all components are included in one
- // single vector. Assumption: first come all
- // entries to the first component, then all
- // entries to the second one, and so on. This
- // is ensured by the way MatrixFree reads
- // out the indices.
- {
- internal::check_vector_compatibility (*src[0], dof_info);
- Assert (n_fe_components == n_components, ExcNotImplemented());
- const unsigned int n_local_dofs =
- dofs_per_cell * VectorizedArray<Number>::n_array_elements * n_components;
- Number * local_src_number = reinterpret_cast<Number*>(values_dofs[0]);
- if (at_irregular_cell == false)
- {
- for (unsigned int j=0; j<n_local_dofs; ++j)
- local_src_number[j] =
- internal::vector_access (*src[0], dof_indices[j]);
- }
- // non-standard case: need to fill in zeros
- // for those components that are not present
- // (a bit more expensive), but there is not
- // more than one such cell
- else
- {
- Assert (n_irreg_components_filled > 0, ExcInternalError());
- for(unsigned int ind_local=0; ind_local<n_local_dofs; ++dof_indices)
- {
- // non-constrained case: copy the data from
- // the global vector, src, to the local one,
- // local_dst.
- local_src_number[ind_local] =
- internal::vector_access (*src[0], *dof_indices);
- ++ind_local;
- while (ind_local % VectorizedArray<Number>::n_array_elements >= n_irreg_components_filled)
- {
- local_src_number[ind_local] = 0.;
- ++ind_local;
- }
- }
- }
- }
+template <int dim, int dofs_per_cell_, int n_q_points_,
+ int n_components_, typename Number>
+template<typename VectorType>
+inline
+void
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::distribute_local_to_global (std::vector<VectorType> &dst,
+ const unsigned int first_index) const
+{
+ AssertIndexRange (first_index, dst.size());
+ Assert (n_fe_components == 1, ExcNotImplemented());
+ Assert ((n_fe_components == 1 ?
+ (first_index+n_components <= dst.size()) : true),
+ ExcIndexRange (first_index + n_components_, 0, dst.size()));
+ Assert (dof_values_initialized==true,
+ internal::ExcAccessToUninitializedField());
-#ifdef DEBUG
- dof_values_initialized = true;
-#endif
+ VectorType * dst_data [n_components];
+ for (unsigned int comp=0; comp<n_components; ++comp)
+ dst_data[comp] = &dst[comp+first_index];
+
+ internal::VectorDistributorLocalToGlobal<Number> distributor;
+ read_write_operation (distributor, dst_data);
}
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-distribute_local_to_global (VectorType &dst) const
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::distribute_local_to_global (std::vector<VectorType*> &dst,
+ const unsigned int first_index) const
{
- AssertDimension (n_components, n_fe_components);
+ AssertIndexRange (first_index, dst.size());
+ Assert (n_fe_components == 1, ExcNotImplemented());
+ Assert ((n_fe_components == 1 ?
+ (first_index+n_components <= dst.size()) : true),
+ ExcIndexRange (first_index + n_components_, 0, dst.size()));
+ Assert (dof_values_initialized==true,
+ internal::ExcAccessToUninitializedField());
+
+ VectorType * dst_data [n_components];
+ for (unsigned int comp=0; comp<n_components; ++comp)
+ dst_data[comp] = dst[comp+first_index];
+
+ internal::VectorDistributorLocalToGlobal<Number> distributor;
+ read_write_operation (distributor, dst_data);
+}
+
+
+
+template <int dim, int dofs_per_cell_, int n_q_points_,
+ int n_components_, typename Number>
+template<typename VectorType>
+inline
+void
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::set_dof_values (VectorType &dst) const
+{
+ AssertDimension (n_components_, n_fe_components);
+ Assert (dof_values_initialized==true,
+ internal::ExcAccessToUninitializedField());
+
// only need one component, but to avoid
// compiler warnings, use n_components copies
// here (but these will not be used)
VectorType * dst_data [n_components];
for (unsigned int d=0; d<n_components; ++d)
dst_data[d] = &dst;
- distribute_local_to_global (dst_data);
+
+ internal::VectorSetter<Number> setter;
+ read_write_operation (setter, dst_data);
}
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-distribute_local_to_global (std::vector<VectorType> &dst,
- const unsigned int first_index) const
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::set_dof_values (std::vector<VectorType> &dst,
+ const unsigned int first_index) const
{
AssertIndexRange (first_index, dst.size());
Assert (n_fe_components == 1, ExcNotImplemented());
Assert ((n_fe_components == 1 ?
(first_index+n_components <= dst.size()) : true),
- ExcIndexRange (first_index + n_components, 0, dst.size()));
+ ExcIndexRange (first_index + n_components_, 0, dst.size()));
+
+ Assert (dof_values_initialized==true,
+ internal::ExcAccessToUninitializedField());
VectorType * dst_data [n_components];
for (unsigned int comp=0; comp<n_components; ++comp)
dst_data[comp] = &dst[comp+first_index];
- distribute_local_to_global (dst_data);
+
+ internal::VectorSetter<Number> setter;
+ read_write_operation (setter, dst_data);
}
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-distribute_local_to_global (std::vector<VectorType*> &dst,
- const unsigned int first_index) const
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::set_dof_values (std::vector<VectorType*> &dst,
+ const unsigned int first_index) const
{
AssertIndexRange (first_index, dst.size());
Assert (n_fe_components == 1, ExcNotImplemented());
Assert ((n_fe_components == 1 ?
(first_index+n_components <= dst.size()) : true),
- ExcIndexRange (first_index + n_components, 0, dst.size()));
+ ExcIndexRange (first_index + n_components_, 0, dst.size()));
+
+ Assert (dof_values_initialized==true,
+ internal::ExcAccessToUninitializedField());
VectorType * dst_data [n_components];
for (unsigned int comp=0; comp<n_components; ++comp)
dst_data[comp] = dst[comp+first_index];
- distribute_local_to_global (dst_data);
+
+ internal::VectorSetter<Number> setter;
+ read_write_operation (setter, dst_data);
}
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-distribute_local_to_global (VectorType * dst[]) const
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::read_dof_values_plain (const VectorType * src[])
{
+ // this is different from the other three
+ // operations because we do not use
+ // constraints here, so this is a separate
+ // function.
Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
- Assert (dof_values_initialized==true,
- internal::ExcAccessToUninitializedField());
+ Assert (dof_info.store_plain_indices == true, ExcNotInitialized());
// loop over all local dofs. ind_local holds
// local number on cell, index iterates over
// the elements of index_local_to_global and
// dof_indices points to the global indices
// stored in index_local_to_global
- const unsigned int * dof_indices = dof_info.begin_indices(cell);
- const std::pair<unsigned short,unsigned short> * indicators =
- dof_info.begin_indicators(cell);
- const std::pair<unsigned short,unsigned short> * indicators_end =
- dof_info.end_indicators(cell);
- unsigned int ind_local = 0;
+ const unsigned int * dof_indices = dof_info.begin_indices_plain(cell);
// scalar case (or case when all components
// have the same degrees of freedom and sit on
const unsigned int n_local_dofs =
VectorizedArray<Number>::n_array_elements * dofs_per_cell;
for (unsigned int comp=0; comp<n_components; ++comp)
- internal::check_vector_compatibility (*dst[comp], dof_info);
-
- const Number * local_dst_number [n_components];
+ internal::check_vector_compatibility (*src[comp], dof_info);
+ Number * local_src_number [n_components];
for (unsigned int comp=0; comp<n_components; ++comp)
- local_dst_number[comp] =
- reinterpret_cast<const Number*>(values_dofs[comp]);
+ local_src_number[comp] = reinterpret_cast<Number*>(values_dofs[comp]);
+
+ // standard case where there are sufficiently
+ // many cells to fill all vectors
if (at_irregular_cell == false)
{
- // check whether there is no constraint at all
- if (indicators != indicators_end)
- {
- // run from one constraint to the next
- for ( ; indicators != indicators_end; ++indicators)
- {
- // distribute values up to the constraint
- // (values not constrained)
- for (unsigned int j=0; j<indicators->first; ++j)
- for (unsigned int comp=0; comp<n_components; ++comp)
- internal::vector_access (*dst[comp], dof_indices[j])
- += local_dst_number[comp][ind_local+j];
- dof_indices += indicators->first;
- ind_local += indicators->first;
-
- // constrained case: build the local value as
- // a linear combination of the global value
- // according to constraint
- const Number * data_val =
- matrix_info.constraint_pool_begin(indicators->second);
- const Number * end_pool =
- matrix_info.constraint_pool_end(indicators->second);
- for ( ; data_val != end_pool; ++data_val, ++dof_indices)
- for (unsigned int comp=0; comp<n_components; ++comp)
- internal::vector_access (*dst[comp], *dof_indices)
- += local_dst_number[comp][ind_local] * (*data_val);
- ++ind_local;
- }
- // distribute values after the last constraint
- // (values not constrained)
- for(; ind_local<n_local_dofs; ++dof_indices, ++ind_local)
- for (unsigned int comp=0; comp<n_components; ++comp)
- internal::vector_access (*dst[comp], *dof_indices)
- += local_dst_number[comp][ind_local];
- }
- // no constraint at all: loop bounds are
- // known, compiler can unroll without checks
- else
- {
- AssertDimension (dof_info.end_indices(cell)-dof_indices,
- static_cast<int>(n_local_dofs));
- for (unsigned int j=0; j<n_local_dofs; ++j)
- for (unsigned int comp=0; comp<n_components; ++comp)
- internal::vector_access (*dst[comp], dof_indices[j])
- += local_dst_number[comp][j];
- }
- return;
+ for (unsigned int j=0; j<n_local_dofs; ++j)
+ for (unsigned int comp=0; comp<n_components; ++comp)
+ local_src_number[comp][j] =
+ internal::vector_access (*src[comp], dof_indices[j]);
}
- // irregular case
- Assert (n_irreg_components_filled > 0, ExcInternalError());
- for ( ; indicators != indicators_end; ++indicators)
+ // non-standard case: need to fill in zeros
+ // for those components that are not present
+ // (a bit more expensive), but there is not
+ // more than one such cell
+ else
{
- for(unsigned int j=0; j<indicators->first; ++j)
+ Assert (n_irreg_components_filled > 0, ExcInternalError());
+ for(unsigned int ind_local=0; ind_local<n_local_dofs;
+ ++dof_indices)
{
// non-constrained case: copy the data from
- // the local vector to the global one.
+ // the global vector, src, to the local one,
+ // local_dst.
for (unsigned int comp=0; comp<n_components; ++comp)
- internal::vector_access (*dst[comp], dof_indices[j])
- += local_dst_number[comp][ind_local];
+ local_src_number[comp][ind_local] =
+ internal::vector_access (*src[comp], *dof_indices);
++ind_local;
- if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
- ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
- }
- dof_indices += indicators->first;
-
- // constrained case: distribute according to
- // the constraint
- const Number * data_val =
- matrix_info.constraint_pool_begin(indicators->second);
- const Number * end_pool =
- matrix_info.constraint_pool_end(indicators->second);
- for ( ; data_val != end_pool; ++data_val, ++dof_indices)
- {
- for (unsigned int comp=0; comp<n_components; ++comp)
- internal::vector_access (*dst[comp], *dof_indices)
- += local_dst_number[comp][ind_local] * (*data_val);
+ while (ind_local % VectorizedArray<Number>::n_array_elements >= n_irreg_components_filled)
+ {
+ for (unsigned int comp=0; comp<n_components; ++comp)
+ local_src_number[comp][ind_local] = 0.;
+ ++ind_local;
+ }
}
- ++ind_local;
- if (ind_local % VectorizedArray<Number>::n_array_elements ==
- n_irreg_components_filled)
- ind_local += VectorizedArray<Number>::n_array_elements-
- n_irreg_components_filled;
- }
- for(; ind_local<n_local_dofs; ++dof_indices)
- {
- Assert (dof_indices != dof_info.end_indices(cell),
- ExcInternalError());
-
- // non-constrained case
- for (unsigned int comp=0; comp<n_components; ++comp)
- internal::vector_access (*dst[comp], *dof_indices)
- += local_dst_number[comp][ind_local];
- ++ind_local;
- if (ind_local % VectorizedArray<Number>::n_array_elements ==
- n_irreg_components_filled)
- ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
}
}
else
// is ensured by the way MatrixFree reads
// out the indices.
{
- internal::check_vector_compatibility (*dst[0], dof_info);
- Assert (n_fe_components == n_components, ExcNotImplemented());
+ internal::check_vector_compatibility (*src[0], dof_info);
+ Assert (n_fe_components == n_components_, ExcNotImplemented());
const unsigned int n_local_dofs =
dofs_per_cell * VectorizedArray<Number>::n_array_elements * n_components;
- const Number * local_dst_number =
- reinterpret_cast<const Number*>(values_dofs[0]);
+ Number * local_src_number = reinterpret_cast<Number*>(values_dofs[0]);
if (at_irregular_cell == false)
{
- // check whether there is no constraint at all
- if (indicators != indicators_end)
- {
- // run from one constraint to the next
- for ( ; indicators != indicators_end; ++indicators)
- {
- // distribute values up to the constraint
- // (values not constrained)
- for (unsigned int j=0; j<indicators->first; ++j)
- internal::vector_access (*dst[0], dof_indices[j])
- += local_dst_number[ind_local+j];
- dof_indices += indicators->first;
- ind_local += indicators->first;
-
- // constrained case: build the local value as
- // a linear combination of the global value
- // according to constraint
- const Number * data_val =
- matrix_info.constraint_pool_begin(indicators->second);
- const Number * end_pool =
- matrix_info.constraint_pool_end(indicators->second);
- for ( ; data_val != end_pool; ++data_val, ++dof_indices)
- internal::vector_access (*dst[0], *dof_indices)
- += local_dst_number[ind_local] * (*data_val);
- ++ind_local;
- }
- // distribute values after the last constraint
- // (values not constrained)
- for(; ind_local<n_local_dofs; ++dof_indices, ++ind_local)
- internal::vector_access (*dst[0], *dof_indices)
- += local_dst_number[ind_local];
- }
- // no constraint at all: loop bounds are
- // known, compiler can unroll without checks
- else
- {
- AssertDimension (dof_info.end_indices(cell)-dof_indices,
- static_cast<int>(n_local_dofs));
- for (unsigned int j=0; j<n_local_dofs; ++j)
- internal::vector_access (*dst[0], dof_indices[j])
- += local_dst_number[j];
- }
- return;
+ for (unsigned int j=0; j<n_local_dofs; ++j)
+ local_src_number[j] =
+ internal::vector_access (*src[0], dof_indices[j]);
}
- // irregular case
- Assert (n_irreg_components_filled > 0, ExcInternalError());
- for ( ; indicators != indicators_end; ++indicators)
+ // non-standard case: need to fill in zeros
+ // for those components that are not present
+ // (a bit more expensive), but there is not
+ // more than one such cell
+ else
{
- for(unsigned int j=0; j<indicators->first; ++j)
+ Assert (n_irreg_components_filled > 0, ExcInternalError());
+ for(unsigned int ind_local=0; ind_local<n_local_dofs; ++dof_indices)
{
// non-constrained case: copy the data from
- // the local vector to the global one.
- internal::vector_access (*dst[0], dof_indices[j])
- += local_dst_number[ind_local];
+ // the global vector, src, to the local one,
+ // local_dst.
+ local_src_number[ind_local] =
+ internal::vector_access (*src[0], *dof_indices);
++ind_local;
- if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
- ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
- }
- dof_indices += indicators->first;
-
- // constrained case: distribute according to
- // the constraint
- const Number * data_val =
- matrix_info.constraint_pool_begin(indicators->second);
- const Number * end_pool =
- matrix_info.constraint_pool_end(indicators->second);
- for ( ; data_val != end_pool; ++data_val, ++dof_indices)
- {
- internal::vector_access (*dst[0], *dof_indices)
- += local_dst_number[ind_local] * (*data_val);
+ while (ind_local % VectorizedArray<Number>::n_array_elements >= n_irreg_components_filled)
+ {
+ local_src_number[ind_local] = 0.;
+ ++ind_local;
+ }
}
- ++ind_local;
- if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
- ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
}
- for(; ind_local<n_local_dofs; ++dof_indices)
- {
- Assert (dof_indices != dof_info.end_indices(cell),
- ExcInternalError());
-
- // non-constrained case
- internal::vector_access (*dst[0], *dof_indices)
- += local_dst_number[ind_local];
- ++ind_local;
- if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
- ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
- }
- Assert (dof_indices == dof_info.end_indices(cell),
- ExcInternalError());
}
+
+#ifdef DEBUG
+ dof_values_initialized = true;
+#endif
}
+
+/*------------------------------ access to data fields ----------------------*/
+
template <int dim, int dofs_per_cell_, int n_q_points_,
int n_components, typename Number>
-template<typename VectorType>
inline
-void
+const VectorizedArray<Number> *
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-set_dof_values (VectorType &dst) const
+begin_values () const
{
- AssertDimension (n_components, n_fe_components);
- // only need one component, but to avoid
- // compiler warnings, use n_components copies
- // here (but these will not be used)
- VectorType * dst_data [n_components];
- for (unsigned int d=0; d<n_components; ++d)
- dst_data[d] = &dst;
- set_dof_values (dst_data);
+ Assert (values_quad_initialized || values_quad_submitted,
+ ExcNotInitialized());
+ return &values_quad[0][0];
}
template <int dim, int dofs_per_cell_, int n_q_points_,
int n_components, typename Number>
-template<typename VectorType>
inline
-void
+VectorizedArray<Number> *
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-set_dof_values (std::vector<VectorType> &dst,
- const unsigned int first_index) const
+begin_values ()
{
- AssertIndexRange (first_index, dst.size());
- Assert (n_fe_components == 1, ExcNotImplemented());
- Assert ((n_fe_components == 1 ?
- (first_index+n_components <= dst.size()) : true),
- ExcIndexRange (first_index + n_components, 0, dst.size()));
-
- VectorType * dst_data [n_components];
- for (unsigned int comp=0; comp<n_components; ++comp)
- dst_data[comp] = &dst[comp+first_index];
- set_dof_values (dst_data);
+#ifdef DEBUG
+ values_quad_submitted = true;
+#endif
+ return &values_quad[0][0];
}
template <int dim, int dofs_per_cell_, int n_q_points_,
int n_components, typename Number>
-template<typename VectorType>
inline
-void
+const VectorizedArray<Number> *
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-set_dof_values (std::vector<VectorType*> &dst,
- const unsigned int first_index) const
+begin_gradients () const
{
- AssertIndexRange (first_index, dst.size());
- Assert (n_fe_components == 1, ExcNotImplemented());
- Assert ((n_fe_components == 1 ?
- (first_index+n_components <= dst.size()) : true),
- ExcIndexRange (first_index + n_components, 0, dst.size()));
-
- VectorType * dst_data [n_components];
- for (unsigned int comp=0; comp<n_components; ++comp)
- dst_data[comp] = dst[comp+first_index];
- set_dof_values (dst_data);
+ Assert (gradients_quad_initialized || gradients_quad_submitted,
+ ExcNotInitialized());
+ return &gradients_quad[0][0][0];
}
template <int dim, int dofs_per_cell_, int n_q_points_,
int n_components, typename Number>
-template<typename VectorType>
inline
-void
+VectorizedArray<Number> *
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-set_dof_values (VectorType * dst[]) const
+begin_gradients ()
{
- Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
- Assert (dof_values_initialized==true,
- internal::ExcAccessToUninitializedField());
+#ifdef DEBUG
+ gradients_quad_submitted = true;
+#endif
+ return &gradients_quad[0][0][0];
+}
- // loop over all local dofs. ind_local holds
- // local number on cell, index iterates over
- // the elements of index_local_to_global and
- // glob_indices points to the global indices
- // stored in index_local_to_global
- const unsigned int * dof_indices = dof_info.begin_indices(cell);
- const std::pair<unsigned short,unsigned short> * indicators =
- dof_info.begin_indicators(cell);
- const std::pair<unsigned short,unsigned short> * indicators_end =
- dof_info.end_indicators(cell);
- unsigned int ind_local = 0;
- if (n_fe_components == 1)
- {
- const unsigned int n_local_dofs =
- VectorizedArray<Number>::n_array_elements * dofs_per_cell;
- for (unsigned int comp=0; comp<n_components; ++comp)
- AssertDimension (dst[comp]->size(),
- dof_info.vector_partitioner->size());
- const Number * local_dst_number [n_components];
- for (unsigned int comp=0; comp<n_components; ++comp)
- local_dst_number[comp] =
- reinterpret_cast<const Number*>(values_dofs[comp]);
- if (at_irregular_cell == false)
- {
- // check whether there is no constraint at all
- if (indicators != indicators_end)
- {
- // run from one constraint to the next
- for ( ; indicators != indicators_end; ++indicators)
- {
- // distribute values up to the constraint
- // (values not constrained)
- for (unsigned int j=0; j<indicators->first; ++j)
- for (unsigned int comp=0; comp<n_components; ++comp)
- internal::vector_access (*dst[comp], dof_indices[j])
- = local_dst_number[comp][ind_local+j];
- dof_indices += indicators->first;
- ind_local += indicators->first;
-
- // jump over constraints
- const unsigned int row_length =
- matrix_info.constraint_pool_end(indicators->second)-
- matrix_info.constraint_pool_begin(indicators->second);
- dof_indices += row_length;
- ++ind_local;
- }
- // distribute values after the last constraint
- // (values not constrained)
- for(; ind_local<n_local_dofs; ++dof_indices, ++ind_local)
- for (unsigned int comp=0; comp<n_components; ++comp)
- internal::vector_access (*dst[comp], *dof_indices)
- = local_dst_number[comp][ind_local];
- }
- // no constraint at all: loop bounds are
- // known, compiler can unroll without checks
- else
- {
- AssertDimension (dof_info.end_indices(cell)-dof_indices,
- n_local_dofs);
- for (unsigned int j=0; j<n_local_dofs; ++j)
- for (unsigned int comp=0; comp<n_components; ++comp)
- internal::vector_access (*dst[comp], dof_indices[j])
- = local_dst_number[comp][j];
- }
- return;
- }
+template <int dim, int dofs_per_cell_, int n_q_points_,
+ int n_components, typename Number>
+inline
+const VectorizedArray<Number> *
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
+begin_hessians () const
+{
+ Assert (hessians_quad_initialized, ExcNotInitialized());
+ return &hessians_quad[0][0][0];
+}
- // irregular case
- Assert (n_irreg_components_filled > 0, ExcInternalError());
- for ( ; indicators != indicators_end; ++indicators)
- {
- for(unsigned int j=0; j<indicators->first; ++j)
- {
- // non-constrained case
- for (unsigned int comp=0; comp<n_components; ++comp)
- internal::vector_access (*dst[comp], dof_indices[j])
- = local_dst_number[comp][ind_local];
- ++ind_local;
- if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
- ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
- }
- dof_indices += indicators->first;
-
- // jump over constraint
- const unsigned int row_length =
- matrix_info.constraint_pool_end(indicators->second)-
- matrix_info.constraint_pool_begin(indicators->second);
- dof_indices += row_length;
- ++ind_local;
- if (ind_local % VectorizedArray<Number>::n_array_elements ==
- n_irreg_components_filled)
- ind_local += VectorizedArray<Number>::n_array_elements -
- n_irreg_components_filled;
- }
- for(; ind_local<n_local_dofs; ++dof_indices)
- {
- Assert (dof_indices != dof_info.end_indices(cell),
- ExcInternalError());
- // non-constrained case
- for (unsigned int comp=0; comp<n_components; ++comp)
- internal::vector_access (*dst[comp], *dof_indices)
- = local_dst_number[comp][ind_local];
- ++ind_local;
- if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
- ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
- }
- }
- else
- // case with vector-valued finite elements
- // where all components are included in one
- // single vector. Assumption: first come all
- // entries to the first component, then all
- // entries to the second one, and so on. This
- // is ensured by the way MatrixFree reads
- // out the indices.
- {
- AssertDimension (dst[0]->size(),
- dof_info.vector_partitioner->size());
- Assert (n_fe_components == n_components, ExcNotImplemented());
- const unsigned int n_local_dofs =
- dofs_per_cell * VectorizedArray<Number>::n_array_elements * n_components;
- const Number * local_dst_number =
- reinterpret_cast<const Number*>(values_dofs[0]);
- if (at_irregular_cell == false)
- {
- // check whether there is no constraint at all
- if (indicators != indicators_end)
- {
- // run from one constraint to the next
- for ( ; indicators != indicators_end; ++indicators)
- {
- // distribute values up to the constraint
- // (values not constrained)
- for (unsigned int j=0; j<indicators->first; ++j)
- internal::vector_access (*dst[0], dof_indices[j])
- = local_dst_number[ind_local+j];
- dof_indices += indicators->first;
- ind_local += indicators->first;
-
- // jump over constraints
- const unsigned int row_length =
- matrix_info.constraint_pool_end(indicators->second) -
- matrix_info.constraint_pool_begin(indicators->second);
- dof_indices += row_length;
- ++ind_local;
- }
- // distribute values after the last constraint
- // (values not constrained)
- for(; ind_local<n_local_dofs; ++dof_indices, ++ind_local)
- internal::vector_access (*dst[0], *dof_indices)
- = local_dst_number[ind_local];
- }
- // no constraint at all: loop bounds are
- // known, compiler can unroll without checks
- else
- {
- AssertDimension (dof_info.end_indices(cell)-dof_indices,
- n_local_dofs);
- for (unsigned int j=0; j<n_local_dofs; ++j)
- internal::vector_access (*dst[0], dof_indices[j])
- = local_dst_number[j];
- }
- return;
- }
- // irregular case
- Assert (n_irreg_components_filled > 0, ExcInternalError());
- for ( ; indicators != indicators_end; ++indicators)
- {
- for(unsigned int j=0; j<indicators->first; ++j)
- {
- // non-constrained case
- internal::vector_access (*dst[0], dof_indices[j])
- = local_dst_number[ind_local];
- ++ind_local;
- if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
- ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
- }
- dof_indices += indicators->first;
-
- // jump over constraint
- const unsigned int row_length =
- matrix_info.constraint_pool_end(indicators->second)-
- matrix_info.constraint_pool_begin(indicators->second);
- dof_indices += row_length;
- ++ind_local;
- if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
- ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
- }
- for(; ind_local<n_local_dofs; ++dof_indices)
- {
- Assert (dof_indices != dof_info.end_indices(cell),
- ExcInternalError());
- // non-constrained case
- internal::vector_access (*dst[0], *dof_indices)
- = local_dst_number[ind_local];
- ++ind_local;
- if (ind_local % VectorizedArray<Number>::n_array_elements == n_irreg_components_filled)
- ind_local += VectorizedArray<Number>::n_array_elements-n_irreg_components_filled;
- }
- Assert (dof_indices == dof_info.end_indices (cell),
- ExcInternalError());
- }
+template <int dim, int dofs_per_cell_, int n_q_points_,
+ int n_components, typename Number>
+inline
+VectorizedArray<Number> *
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
+begin_hessians ()
+{
+ return &hessians_quad[0][0][0];
}
-// ------------------------------ access to data fields ---------------------
-
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
-Tensor<1,n_components,VectorizedArray<Number> >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_dof_value (const unsigned int dof) const
+Tensor<1,n_components_,VectorizedArray<Number> >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_dof_value (const unsigned int dof) const
{
AssertIndexRange (dof, dofs_per_cell);
- Tensor<1,n_components,VectorizedArray<Number> > return_value (false);
+ Tensor<1,n_components_,VectorizedArray<Number> > return_value (false);
for(unsigned int comp=0;comp<n_components;comp++)
return_value[comp] = this->values_dofs[comp][dof];
return return_value;
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
-Tensor<1,n_components,VectorizedArray<Number> >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_value (const unsigned int q_point) const
+Tensor<1,n_components_,VectorizedArray<Number> >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_value (const unsigned int q_point) const
{
Assert (this->values_quad_initialized==true,
internal::ExcAccessToUninitializedField());
AssertIndexRange (q_point, n_q_points);
- Tensor<1,n_components,VectorizedArray<Number> > return_value (false);
+ Tensor<1,n_components_,VectorizedArray<Number> > return_value (false);
for(unsigned int comp=0;comp<n_components;comp++)
return_value[comp] = this->values_quad[comp][q_point];
return return_value;
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
-Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_gradient (const unsigned int q_point) const
+Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_gradient (const unsigned int q_point) const
{
Assert (this->gradients_quad_initialized==true,
internal::ExcAccessToUninitializedField());
AssertIndexRange (q_point, n_q_points);
- Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > > grad_out (false);
+ Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > grad_out (false);
// Cartesian cell
if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
grad_out[comp][d] = (this->gradients_quad[comp][d][q_point] *
cartesian_data[0][d]);
}
- // cell with general Jacobian
- else if (this->cell_type == internal::MatrixFreeFunctions::general)
- {
- for(unsigned int comp=0;comp<n_components;comp++)
- {
- for (unsigned int d=0; d<dim; ++d)
- {
- grad_out[comp][d] = (jacobian[q_point][d][0] *
- this->gradients_quad[comp][0][q_point]);
- for (unsigned e=1; e<dim; ++e)
- grad_out[comp][d] += (jacobian[q_point][d][e] *
- this->gradients_quad[comp][e][q_point]);
- }
- }
- }
- // cell with general Jacobian, but constant
- // within the cell
- else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
+ // cell with general/affine Jacobian
+ else
{
+ const Tensor<2,dim,VectorizedArray<Number> > & jac =
+ this->cell_type == internal::MatrixFreeFunctions::general ?
+ jacobian[q_point] : jacobian[0];
for(unsigned int comp=0;comp<n_components;comp++)
{
for (unsigned int d=0; d<dim; ++d)
{
- grad_out[comp][d] = (jacobian[0][d][0] *
+ grad_out[comp][d] = (jac[d][0] *
this->gradients_quad[comp][0][q_point]);
for (unsigned e=1; e<dim; ++e)
- grad_out[comp][d] += (jacobian[0][d][e] *
+ grad_out[comp][d] += (jac[d][e] *
this->gradients_quad[comp][e][q_point]);
}
}
+namespace internal
+{
+ // compute tmp = hess_unit(u) * J^T. do this
+ // manually because we do not store the lower
+ // diagonal because of symmetry
+ template <int dim, int n_q_points, typename Number>
+ inline
+ void
+ hessian_unit_times_jac (const Tensor<2,dim,VectorizedArray<Number> > &jac,
+ const VectorizedArray<Number> hessians_quad[][n_q_points],
+ const unsigned int q_point,
+ VectorizedArray<Number> tmp[dim][dim])
+ {
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ switch (dim)
+ {
+ case 1:
+ tmp[0][0] = jac[0][0] * hessians_quad[0][q_point];
+ break;
+ case 2:
+ tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
+ jac[d][1] * hessians_quad[2][q_point]);
+ tmp[1][d] = (jac[d][0] * hessians_quad[2][q_point] +
+ jac[d][1] * hessians_quad[1][q_point]);
+ break;
+ case 3:
+ tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
+ jac[d][1] * hessians_quad[3][q_point] +
+ jac[d][2] * hessians_quad[4][q_point]);
+ tmp[1][d] = (jac[d][0] * hessians_quad[3][q_point] +
+ jac[d][1] * hessians_quad[1][q_point] +
+ jac[d][2] * hessians_quad[5][q_point]);
+ tmp[2][d] = (jac[d][0] * hessians_quad[4][q_point] +
+ jac[d][1] * hessians_quad[5][q_point] +
+ jac[d][2] * hessians_quad[2][q_point]);
+ break;
+ default: Assert (false, ExcNotImplemented());
+ }
+ }
+ }
+}
+
+
+
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
-Tensor<1,n_components,Tensor<2,dim,VectorizedArray<Number> > >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_hessian (const unsigned int q_point) const
+Tensor<1,n_components_,Tensor<2,dim,VectorizedArray<Number> > >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_hessian (const unsigned int q_point) const
{
Assert (this->hessians_quad_initialized==true,
internal::ExcAccessToUninitializedField());
// because it needs to access unscaled
// gradient data
VectorizedArray<Number> tmp[dim][dim];
+ internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
+ q_point, tmp);
- // compute tmp = hess_unit(u) * J^T. do this
- // manually because we do not store the lower
- // diagonal because of symmetry
- for (unsigned int d=0; d<dim; ++d)
- {
- switch (dim)
- {
- case 1:
- tmp[0][0] = jac[0][0] * this->hessians_quad[comp][0][q_point];
- break;
- case 2:
- tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
- jac[d][1] * this->hessians_quad[comp][2][q_point]);
- tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][2][q_point] +
- jac[d][1] * this->hessians_quad[comp][1][q_point]);
- break;
- case 3:
- tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
- jac[d][1] * this->hessians_quad[comp][3][q_point] +
- jac[d][2] * this->hessians_quad[comp][4][q_point]);
- tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][3][q_point] +
- jac[d][1] * this->hessians_quad[comp][1][q_point] +
- jac[d][2] * this->hessians_quad[comp][5][q_point]);
- tmp[2][d] = (jac[d][0] * this->hessians_quad[comp][4][q_point] +
- jac[d][1] * this->hessians_quad[comp][5][q_point] +
- jac[d][2] * this->hessians_quad[comp][2][q_point]);
- break;
- default: Assert (false, ExcNotImplemented());
- }
- }
// compute first part of hessian,
// J * tmp = J * hess_unit(u) * J^T
for (unsigned int d=0; d<dim; ++d)
// because it needs to access unscaled
// gradient data
VectorizedArray<Number> tmp[dim][dim];
+ internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
+ q_point, tmp);
- // compute tmp = hess_unit(u) * J^T. do this
- // manually because we do not store the lower
- // diagonal because of symmetry
- for (unsigned int d=0; d<dim; ++d)
- {
- switch (dim)
- {
- case 1:
- tmp[0][0] = jac[0][0] * this->hessians_quad[comp][0][q_point];
- break;
- case 2:
- tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
- jac[d][1] * this->hessians_quad[comp][2][q_point]);
- tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][2][q_point] +
- jac[d][1] * this->hessians_quad[comp][1][q_point]);
- break;
- case 3:
- tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
- jac[d][1] * this->hessians_quad[comp][3][q_point] +
- jac[d][2] * this->hessians_quad[comp][4][q_point]);
- tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][3][q_point] +
- jac[d][1] * this->hessians_quad[comp][1][q_point] +
- jac[d][2] * this->hessians_quad[comp][5][q_point]);
- tmp[2][d] = (jac[d][0] * this->hessians_quad[comp][4][q_point] +
- jac[d][1] * this->hessians_quad[comp][5][q_point] +
- jac[d][2] * this->hessians_quad[comp][2][q_point]);
- break;
- default: Assert (false, ExcNotImplemented());
- }
- }
// compute first part of hessian,
// J * tmp = J * hess_unit(u) * J^T
for (unsigned int d=0; d<dim; ++d)
hessian_out[comp][e][d] = hessian_out[comp][d][e];
}
}
- return Tensor<1,n_components,Tensor<2,dim,VectorizedArray<Number> > >(hessian_out);
+ return Tensor<1,n_components_,Tensor<2,dim,VectorizedArray<Number> > >(hessian_out);
}
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
-Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_hessian_diagonal (const unsigned int q_point) const
+Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_hessian_diagonal (const unsigned int q_point) const
{
Assert (this->hessians_quad_initialized==true,
internal::ExcAccessToUninitializedField());
AssertIndexRange (q_point, n_q_points);
- Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > > hessian_out (false);
+ Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > hessian_out (false);
// Cartesian cell
if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
// because it needs to access unscaled
// gradient data
VectorizedArray<Number> tmp[dim][dim];
+ internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
+ q_point, tmp);
- // compute tmp = hess_unit(u) * J^T. do this
- // manually because we do not store the lower
- // diagonal because of symmetry
- for (unsigned int d=0; d<dim; ++d)
- {
- switch (dim)
- {
- case 1:
- tmp[0][0] = jac[0][0] * this->hessians_quad[comp][0][q_point];
- break;
- case 2:
- tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
- jac[d][1] * this->hessians_quad[comp][2][q_point]);
- tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][2][q_point] +
- jac[d][1] * this->hessians_quad[comp][1][q_point]);
- break;
- case 3:
- tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
- jac[d][1] * this->hessians_quad[comp][3][q_point] +
- jac[d][2] * this->hessians_quad[comp][4][q_point]);
- tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][3][q_point] +
- jac[d][1] * this->hessians_quad[comp][1][q_point] +
- jac[d][2] * this->hessians_quad[comp][5][q_point]);
- tmp[2][d] = (jac[d][0] * this->hessians_quad[comp][4][q_point] +
- jac[d][1] * this->hessians_quad[comp][5][q_point] +
- jac[d][2] * this->hessians_quad[comp][2][q_point]);
- break;
- default: Assert (false, ExcNotImplemented());
- }
- }
// compute only the trace part of hessian,
// J * tmp = J * hess_unit(u) * J^T
for (unsigned int d=0; d<dim; ++d)
// because it needs to access unscaled
// gradient data
VectorizedArray<Number> tmp[dim][dim];
+ internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
+ q_point, tmp);
- // compute tmp = hess_unit(u) * J^T. do this
- // manually because we do not store the lower
- // diagonal because of symmetry
- for (unsigned int d=0; d<dim; ++d)
- {
- switch (dim)
- {
- case 1:
- tmp[0][0] = jac[0][0] * this->hessians_quad[comp][0][q_point];
- break;
- case 2:
- tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
- jac[d][1] * this->hessians_quad[comp][2][q_point]);
- tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][2][q_point] +
- jac[d][1] * this->hessians_quad[comp][1][q_point]);
- break;
- case 3:
- tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] +
- jac[d][1] * this->hessians_quad[comp][3][q_point] +
- jac[d][2] * this->hessians_quad[comp][4][q_point]);
- tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][3][q_point] +
- jac[d][1] * this->hessians_quad[comp][1][q_point] +
- jac[d][2] * this->hessians_quad[comp][5][q_point]);
- tmp[2][d] = (jac[d][0] * this->hessians_quad[comp][4][q_point] +
- jac[d][1] * this->hessians_quad[comp][5][q_point] +
- jac[d][2] * this->hessians_quad[comp][2][q_point]);
- break;
- default: Assert (false, ExcNotImplemented());
- }
- }
// compute only the trace part of hessian,
// J * tmp = J * hess_unit(u) * J^T
for (unsigned int d=0; d<dim; ++d)
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
-Tensor<1,n_components,VectorizedArray<Number> >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-get_laplacian (const unsigned int q_point) const
+Tensor<1,n_components_,VectorizedArray<Number> >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::get_laplacian (const unsigned int q_point) const
{
Assert (this->hessians_quad_initialized==true,
internal::ExcAccessToUninitializedField());
AssertIndexRange (q_point, n_q_points);
- Tensor<1,n_components,VectorizedArray<Number> > laplacian_out (false);
- const Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number> > > hess_diag
+ Tensor<1,n_components_,VectorizedArray<Number> > laplacian_out (false);
+ const Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > hess_diag
= get_hessian_diagonal(q_point);
for (unsigned int comp=0; comp<n_components; ++comp)
{
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-submit_dof_value (const Tensor<1,n_components,VectorizedArray<Number> > val_in,
- const unsigned int dof)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::submit_dof_value (const Tensor<1,n_components_,VectorizedArray<Number> > val_in,
+ const unsigned int dof)
{
#ifdef DEBUG
this->dof_values_initialized = true;
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-submit_value (const Tensor<1,n_components,VectorizedArray<Number> > val_in,
- const unsigned int q_point)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::submit_value (const Tensor<1,n_components_,VectorizedArray<Number> > val_in,
+ const unsigned int q_point)
{
#ifdef DEBUG
Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
void
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-submit_gradient (const Tensor<1,n_components,
- Tensor<1,dim,VectorizedArray<Number> > > grad_in,
- const unsigned int q_point)
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::submit_gradient (const Tensor<1,n_components_,
+ Tensor<1,dim,VectorizedArray<Number> > >grad_in,
+ const unsigned int q_point)
{
#ifdef DEBUG
Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
this->gradients_quad[comp][d][q_point] = (grad_in[comp][d] *
cartesian_data[0][d] * JxW);
}
- else if (this->cell_type == internal::MatrixFreeFunctions::general)
- {
- for (unsigned int comp=0; comp<n_components; ++comp)
- for (unsigned int d=0; d<dim; ++d)
- {
- VectorizedArray<Number> new_val = jacobian[q_point][0][d] * grad_in[comp][0];
- for (unsigned e=1; e<dim; ++e)
- new_val += jacobian[q_point][e][d] * grad_in[comp][e];
- this->gradients_quad[comp][d][q_point] = new_val * J_value[q_point];
- }
- }
- else //if (this->cell_type == internal::MatrixFreeFunctions::affine)
+ else
{
- const VectorizedArray<Number> JxW = J_value[0] * quadrature_weights[q_point];
+ const Tensor<2,dim,VectorizedArray<Number> > &jac =
+ this->cell_type == internal::MatrixFreeFunctions::general ?
+ jacobian[q_point] : jacobian[0];
+ const VectorizedArray<Number> JxW =
+ this->cell_type == internal::MatrixFreeFunctions::general ?
+ J_value[q_point] : J_value[0] * quadrature_weights[q_point];
for (unsigned int comp=0; comp<n_components; ++comp)
for (unsigned int d=0; d<dim; ++d)
{
- VectorizedArray<Number> new_val = jacobian[0][0][d] * grad_in[comp][0];
+ VectorizedArray<Number> new_val = jac[0][d] * grad_in[comp][0];
for (unsigned e=1; e<dim; ++e)
- new_val += jacobian[0][e][d] * grad_in[comp][e];
+ new_val += (jac[e][d] * grad_in[comp][e]);
this->gradients_quad[comp][d][q_point] = new_val * JxW;
}
}
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
-Tensor<1,n_components,VectorizedArray<Number> >
-FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-integrate_value () const
+Tensor<1,n_components_,VectorizedArray<Number> >
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::integrate_value () const
{
#ifdef DEBUG
Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
Assert (this->values_quad_submitted == true,
internal::ExcAccessToUninitializedField());
#endif
- Tensor<1,n_components,VectorizedArray<Number> > return_value (false);
+ Tensor<1,n_components_,VectorizedArray<Number> > return_value (false);
for (unsigned int comp=0; comp<n_components; ++comp)
return_value[comp] = this->values_quad[comp][0];
- for (unsigned int q=0; q<n_q_points; ++q)
+ for (unsigned int q=1; q<n_q_points; ++q)
for (unsigned int comp=0; comp<n_components; ++comp)
return_value[comp] += this->values_quad[comp][q];
return (return_value);
-/*----------------------- FEEvaluationAccess -------------------------------*/
+/*----------------------- FEEvaluationAccess --------------------------------*/
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components_, typename Number>
inline
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
-FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
- const unsigned int fe_no,
- const unsigned int quad_no_in)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
+::FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
+ const unsigned int fe_no,
+ const unsigned int quad_no_in)
:
- FEEvaluationBase <dim,dofs_per_cell_,n_q_points_,n_components,Number>
+ FEEvaluationBase <dim,dofs_per_cell_,n_q_points_,n_components_,Number>
(data_in, fe_no, quad_no_in)
{}
-/*-------------------- FEEvaluationAccess scalar --------------------------*/
+/*-------------------- FEEvaluationAccess scalar ----------------------------*/
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
- const unsigned int fe_no,
- const unsigned int quad_no_in)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
+ const unsigned int fe_no,
+ const unsigned int quad_no_in)
:
FEEvaluationBase <dim,dofs_per_cell_,n_q_points_,1,Number>
(data_in, fe_no, quad_no_in)
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
VectorizedArray<Number>
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-get_dof_value (const unsigned int dof) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::get_dof_value (const unsigned int dof) const
{
AssertIndexRange (dof, dofs_per_cell);
return this->values_dofs[0][dof];
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
VectorizedArray<Number>
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-get_value (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::get_value (const unsigned int q_point) const
{
Assert (this->values_quad_initialized==true,
internal::ExcAccessToUninitializedField());
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
Tensor<1,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-get_gradient (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::get_gradient (const unsigned int q_point) const
{
// could use the base class gradient, but that
// involves too many inefficient
grad_out[d] = (this->gradients_quad[0][d][q_point] *
this->cartesian_data[0][d]);
}
- // cell with general Jacobian
- else if (this->cell_type == internal::MatrixFreeFunctions::general)
- {
- for (unsigned int d=0; d<dim; ++d)
- {
- grad_out[d] = (this->jacobian[q_point][d][0] *
- this->gradients_quad[0][0][q_point]);
- for (unsigned e=1; e<dim; ++e)
- grad_out[d] += (this->jacobian[q_point][d][e] *
- this->gradients_quad[0][e][q_point]);
- }
- }
- // cell with general Jacobian, but constant
- // within the cell
- else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
+ // cell with general/constant Jacobian
+ else
{
+ const Tensor<2,dim,VectorizedArray<Number> > &jac =
+ this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->jacobian[q_point] : this->jacobian[0];
for (unsigned int d=0; d<dim; ++d)
{
- grad_out[d] = (this->jacobian[0][d][0] *
- this->gradients_quad[0][0][q_point]);
+ grad_out[d] = (jac[d][0] * this->gradients_quad[0][0][q_point]);
for (unsigned e=1; e<dim; ++e)
- grad_out[d] += (this->jacobian[0][d][e] *
- this->gradients_quad[0][e][q_point]);
+ grad_out[d] += (jac[d][e] * this->gradients_quad[0][e][q_point]);
}
}
return grad_out;
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
Tensor<2,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-get_hessian (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::get_hessian (const unsigned int q_point) const
{
return BaseClass::get_hessian(q_point)[0];
}
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
Tensor<1,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-get_hessian_diagonal (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::get_hessian_diagonal (const unsigned int q_point) const
{
return BaseClass::get_hessian_diagonal(q_point)[0];
}
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
VectorizedArray<Number>
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-get_laplacian (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::get_laplacian (const unsigned int q_point) const
{
return BaseClass::get_laplacian(q_point)[0];
}
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-submit_dof_value (const VectorizedArray<Number> val_in,
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::submit_dof_value (const VectorizedArray<Number> val_in,
const unsigned int dof)
{
#ifdef DEBUG
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-submit_value (const VectorizedArray<Number> val_in,
- const unsigned int q_point)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::submit_value (const VectorizedArray<Number> val_in,
+ const unsigned int q_point)
{
#ifdef DEBUG
Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-submit_gradient (const Tensor<1,dim,VectorizedArray<Number> > grad_in,
- const unsigned int q_point)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::submit_gradient (const Tensor<1,dim,VectorizedArray<Number> > grad_in,
+ const unsigned int q_point)
{
#ifdef DEBUG
Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
this->cartesian_data[0][d] *
JxW);
}
- else if (this->cell_type == internal::MatrixFreeFunctions::general)
- {
- for (unsigned int d=0; d<dim; ++d)
- {
- VectorizedArray<Number> new_val = this->jacobian[q_point][0][d] * grad_in[0];
- for (unsigned e=1; e<dim; ++e)
- new_val += this->jacobian[q_point][e][d] * grad_in[e];
- this->gradients_quad[0][d][q_point] = new_val * this->J_value[q_point];
- }
- }
- else //if (this->cell_type == internal::MatrixFreeFunctions::affine)
+ // general/affine cell type
+ else
{
- const VectorizedArray<Number> JxW = this->J_value[0] * this->quadrature_weights[q_point];
+ const Tensor<2,dim,VectorizedArray<Number> > & jac =
+ this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->jacobian[q_point] : this->jacobian[0];
+ const VectorizedArray<Number> JxW =
+ this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
for (unsigned int d=0; d<dim; ++d)
{
- VectorizedArray<Number> new_val = this->jacobian[0][0][d] * grad_in[0];
+ VectorizedArray<Number> new_val = jac[0][d] * grad_in[0];
for (unsigned e=1; e<dim; ++e)
- new_val += this->jacobian[0][e][d] * grad_in[e];
+ new_val += jac[e][d] * grad_in[e];
this->gradients_quad[0][d][q_point] = new_val * JxW;
}
}
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
VectorizedArray<Number>
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>::
-integrate_value () const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
+::integrate_value () const
{
return BaseClass::integrate_value()[0];
}
-/*----------------- FEEvaluationAccess vector-valued ----------------------*/
+/*----------------- FEEvaluationAccess vector-valued ------------------------*/
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
- const unsigned int fe_no,
- const unsigned int quad_no_in)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
+ const unsigned int fe_no,
+ const unsigned int quad_no_in)
:
FEEvaluationBase <dim,dofs_per_cell_,n_q_points_,dim,Number>
(data_in, fe_no, quad_no_in)
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
Tensor<2,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-get_gradient (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::get_gradient (const unsigned int q_point) const
{
return BaseClass::get_gradient (q_point);
}
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
VectorizedArray<Number>
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-get_divergence (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::get_divergence (const unsigned int q_point) const
{
Assert (this->gradients_quad_initialized==true,
internal::ExcAccessToUninitializedField());
divergence += (this->gradients_quad[d][d][q_point] *
this->cartesian_data[0][d]);
}
- // cell with general Jacobian
- else if (this->cell_type == internal::MatrixFreeFunctions::general)
- {
- divergence = (this->jacobian[q_point][0][0] *
- this->gradients_quad[0][0][q_point]);
- for (unsigned e=1; e<dim; ++e)
- divergence += (this->jacobian[q_point][0][e] *
- this->gradients_quad[0][e][q_point]);
- for (unsigned int d=1; d<dim; ++d)
- for (unsigned e=0; e<dim; ++e)
- divergence += (this->jacobian[q_point][d][e] *
- this->gradients_quad[d][e][q_point]);
- }
- // cell with general Jacobian, but constant
- // within the cell
- else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
+ // cell with general/constant Jacobian
+ else
{
- divergence = (this->jacobian[0][0][0] *
- this->gradients_quad[0][0][q_point]);
+ const Tensor<2,dim,VectorizedArray<Number> > &jac =
+ this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->jacobian[q_point] : this->jacobian[0];
+ divergence = (jac[0][0] * this->gradients_quad[0][0][q_point]);
for (unsigned e=1; e<dim; ++e)
- divergence += (this->jacobian[0][0][e] *
- this->gradients_quad[0][e][q_point]);
+ divergence += (jac[0][e] * this->gradients_quad[0][e][q_point]);
for (unsigned int d=1; d<dim; ++d)
for (unsigned e=0; e<dim; ++e)
- divergence += (this->jacobian[0][d][e] *
- this->gradients_quad[d][e][q_point]);
+ divergence += (jac[d][e] * this->gradients_quad[d][e][q_point]);
}
return divergence;
}
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
SymmetricTensor<2,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-get_symmetric_gradient (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::get_symmetric_gradient (const unsigned int q_point) const
{
// copy from generic function into
// dim-specialization function
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
Tensor<1,dim==2?1:dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-get_curl (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::get_curl (const unsigned int q_point) const
{
// copy from generic function into
// dim-specialization function
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
Tensor<2,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-get_hessian_diagonal (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::get_hessian_diagonal (const unsigned int q_point) const
{
Assert (this->hessians_quad_initialized==true,
internal::ExcAccessToUninitializedField());
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
Tensor<3,dim,VectorizedArray<Number> >
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-get_hessian (const unsigned int q_point) const
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::get_hessian (const unsigned int q_point) const
{
Assert (this->hessians_quad_initialized==true,
internal::ExcAccessToUninitializedField());
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-submit_gradient (const Tensor<2,dim,VectorizedArray<Number> > grad_in,
- const unsigned int q_point)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::submit_gradient (const Tensor<2,dim,VectorizedArray<Number> > grad_in,
+ const unsigned int q_point)
{
BaseClass::submit_gradient (grad_in, q_point);
}
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-submit_gradient (const Tensor<1,dim,Tensor<1,dim,VectorizedArray<Number> > > grad_in,
- const unsigned int q_point)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::submit_gradient (const Tensor<1,dim,Tensor<1,dim,VectorizedArray<Number> > >
+ grad_in,
+ const unsigned int q_point)
{
BaseClass::submit_gradient(grad_in, q_point);
}
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-submit_symmetric_gradient (const SymmetricTensor<2,dim,VectorizedArray<Number> >
- sym_grad,
- const unsigned int q_point)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::submit_symmetric_gradient(const SymmetricTensor<2,dim,VectorizedArray<Number> >
+ sym_grad,
+ const unsigned int q_point)
{
// could have used base class operator, but
// that involves some overhead which is
this->cartesian_data[0][e]);
}
}
- else if (this->cell_type == internal::MatrixFreeFunctions::general)
+ // general/affine cell type
+ else
{
+ const VectorizedArray<Number> JxW =
+ this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
+ const Tensor<2,dim,VectorizedArray<Number> > &jac =
+ this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->jacobian[q_point] : this->jacobian[0];
VectorizedArray<Number> weighted [dim][dim];
- {
- const VectorizedArray<Number> JxW = this->J_value[q_point];
- for (unsigned int i=0; i<dim; ++i)
- weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
- for (unsigned int i=0, counter=dim; i<dim; ++i)
- for (unsigned int j=i+1; j<dim; ++j, ++counter)
- {
- const VectorizedArray<Number> value = sym_grad.access_raw_entry(counter) * JxW;
- weighted[i][j] = value;
- weighted[j][i] = value;
- }
- }
- for (unsigned int comp=0; comp<dim; ++comp)
- for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<dim; ++i)
+ weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
+ for (unsigned int i=0, counter=dim; i<dim; ++i)
+ for (unsigned int j=i+1; j<dim; ++j, ++counter)
{
- VectorizedArray<Number> new_val = this->jacobian[q_point][0][d] * weighted[comp][0];
- for (unsigned e=1; e<dim; ++e)
- new_val += this->jacobian[q_point][e][d] * weighted[comp][e];
- this->gradients_quad[comp][d][q_point] = new_val;
+ const VectorizedArray<Number> value = sym_grad.access_raw_entry(counter) * JxW;
+ weighted[i][j] = value;
+ weighted[j][i] = value;
}
- }
- else //if (this->cell_type == internal::MatrixFreeFunctions::affine)
- {
- VectorizedArray<Number> weighted [dim][dim];
- {
- const VectorizedArray<Number> JxW = (this->J_value[0] *
- this->quadrature_weights[q_point]);
- for (unsigned int i=0; i<dim; ++i)
- weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
- for (unsigned int i=0, counter=dim; i<dim; ++i)
- for (unsigned int j=i+1; j<dim; ++j, ++counter)
- {
- const VectorizedArray<Number> value = sym_grad.access_raw_entry(counter) * JxW;
- weighted[i][j] = value;
- weighted[j][i] = value;
- }
- }
for (unsigned int comp=0; comp<dim; ++comp)
for (unsigned int d=0; d<dim; ++d)
{
- VectorizedArray<Number> new_val = this->jacobian[q_point][0][d] * weighted[comp][0];
+ VectorizedArray<Number> new_val = jac[0][d] * weighted[comp][0];
for (unsigned e=1; e<dim; ++e)
- new_val += this->jacobian[q_point][e][d] * weighted[comp][e];
+ new_val += jac[e][d] * weighted[comp][e];
this->gradients_quad[comp][d][q_point] = new_val;
}
}
template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
inline
void
-FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>::
-submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray<Number> > curl,
- const unsigned int q_point)
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray<Number> > curl,
+ const unsigned int q_point)
{
Tensor<2,dim,VectorizedArray<Number> > grad;
switch (dim)
/*----------------------- FEEvaluationGeneral -------------------------------*/
-template <int dim, int fe_degree, int n_q_points_1d, int n_components,
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
typename Number>
inline
-FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components,Number>::
-FEEvaluationGeneral (const MatrixFree<dim,Number> &data_in,
- const unsigned int fe_no,
- const unsigned int quad_no_in)
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationGeneral (const MatrixFree<dim,Number> &data_in,
+ const unsigned int fe_no,
+ const unsigned int quad_no_in)
:
BaseClass (data_in, fe_no, quad_no_in)
{
-template <int dim, int fe_degree, int n_q_points_1d, int n_components,
- typename Number>
-inline
-void
-FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components,Number>::
-evaluate (bool evaluate_val, bool evaluate_grad, bool evaluate_lapl)
+namespace internal
{
- Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
- Assert (this->dof_values_initialized == true,
- internal::ExcAccessToUninitializedField());
+ // evaluates the given shape data in 1d-3d
+ // using the tensor product form. does not use
+ // the tensor product form and corresponds to
+ // a usual matrix-matrix product
+ template <int dim, int fe_degree, int n_q_points_1d, typename Number,
+ int direction, bool dof_to_quad, bool add>
+ inline
+ void
+ apply_tensor_product (const VectorizedArray<Number>*shape_data,
+ const VectorizedArray<Number> in [],
+ VectorizedArray<Number> out [])
+ {
+ AssertIndexRange (direction, dim);
+ const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
+ nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
- const VectorizedArray<Number> * val = this->data.shape_values.begin();
- const VectorizedArray<Number> * grad = this->data.shape_gradients.begin();
- const VectorizedArray<Number> * hess = this->data.shape_hessians.begin();
+ const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
+ const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
+ const int stride = Utilities::fixed_int_power<nn,direction>::value;
- for(unsigned int comp=0;comp<n_components;comp++)
+ for (int i2=0; i2<n_blocks2; ++i2)
+ {
+ for (int i1=0; i1<n_blocks1; ++i1)
+ {
+ for (int col=0; col<nn; ++col)
+ {
+ VectorizedArray<Number> val0;
+ if (dof_to_quad == true)
+ val0 = shape_data[col];
+ else
+ val0 = shape_data[col*n_q_points_1d];
+ VectorizedArray<Number> res0 = val0 * in[0];
+ for (int ind=1; ind<mm; ++ind)
+ {
+ if (dof_to_quad == true)
+ val0 = shape_data[ind*n_q_points_1d+col];
+ else
+ val0 = shape_data[col*n_q_points_1d+ind];
+ res0 += val0 * in[stride*ind];
+ }
+ if (add == false)
+ out[stride*col] = res0;
+ else
+ out[stride*col] += res0;
+ }
+
+ // increment: in regular case, just go to the
+ // next point in x-direction. If we are at the
+ // end of one chunk in x-dir, need to jump
+ // over to the next layer in z-direction
+ switch (direction)
+ {
+ case 0:
+ in += mm;
+ out += nn;
+ break;
+ case 1:
+ case 2:
+ ++in;
+ ++out;
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ }
+ if (direction == 1)
+ {
+ in += nn*(mm-1);
+ out += nn*(nn-1);
+ }
+ }
+ }
+
+
+
+ // This performs the evaluation of function
+ // values, gradients and Hessians for
+ // tensor-product finite elements. The
+ // operation is used for both
+ // FEEvaluationGeneral and FEEvaluation, which
+ // provide different functions apply_values,
+ // apply_gradients in the individual
+ // coordinate directions
+ template <typename FEEval>
+ inline
+ void
+ do_evaluate (FEEval &fe_eval,
+ const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_lapl)
+ {
+ Assert (fe_eval.cell != numbers::invalid_unsigned_int,
+ ExcNotInitialized());
+ Assert (fe_eval.dof_values_initialized == true,
+ internal::ExcAccessToUninitializedField());
+
+ const unsigned int temp_size = fe_eval.dofs_per_cell > fe_eval.n_q_points ?
+ fe_eval.dofs_per_cell : fe_eval.n_q_points;
+ const unsigned int n_components = fe_eval.n_components;
+ const unsigned int dim = fe_eval.dimension;
+
+ for(unsigned int c=0; c<n_components; c++)
{
- VectorizedArray<Number> temp1[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
- VectorizedArray<Number> temp2[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
+ VectorizedArray<typename FEEval::number_type> temp1[temp_size];
+ VectorizedArray<typename FEEval::number_type> temp2[temp_size];
- if (dim == 3)
- {
- if (evaluate_grad == true)
- {
- // grad x
- apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], temp1);
- apply_tensor_prod<1,true,false> (val, temp1, temp2);
- apply_tensor_prod<2,true,false> (val, temp2, this->gradients_quad[comp][0]);
- }
+ switch (dim)
+ {
+ case 3:
- if (evaluate_lapl == true)
- {
- // grad xz
- if (evaluate_grad == false)
- {
- apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], temp1);
- apply_tensor_prod<1,true,false> (val, temp1, temp2);
- }
- apply_tensor_prod<2,true,false>(grad, temp2, this->hessians_quad[comp][4]);
+ if (evaluate_grad == true)
+ {
+ // grad x
+ fe_eval.template apply_gradients<0,true,false>
+ (fe_eval.values_dofs[c], temp1);
+ fe_eval.template apply_values<1,true,false>
+ (temp1, temp2);
+ fe_eval.template apply_values<2,true,false>
+ (temp2, fe_eval.gradients_quad[c][0]);
+ }
- // grad xy
- apply_tensor_prod<1,true,false>(grad, temp1, temp2);
- apply_tensor_prod<2,true,false> (val, temp2, this->hessians_quad[comp][3]);
+ if (evaluate_lapl == true)
+ {
+ // grad xz
+ if (evaluate_grad == false)
+ {
+ fe_eval.template apply_gradients<0,true,false>
+ (fe_eval.values_dofs[c], temp1);
+ fe_eval.template apply_values<1,true,false>
+ (temp1, temp2);
+ }
+ fe_eval.template apply_gradients<2,true,false>
+ (temp2, fe_eval.hessians_quad[c][4]);
+
+ // grad xy
+ fe_eval.template apply_gradients<1,true,false>
+ (temp1, temp2);
+ fe_eval.template apply_values<2,true,false>
+ (temp2, fe_eval.hessians_quad[c][3]);
+
+ // grad xx
+ fe_eval.template apply_hessians<0,true,false>
+ (fe_eval.values_dofs[c], temp1);
+ fe_eval.template apply_values<1,true,false>
+ (temp1, temp2);
+ fe_eval.template apply_values<2,true,false>
+ (temp2, fe_eval.hessians_quad[c][0]);
+ }
- // grad xx
- apply_tensor_prod<0,true,false> (hess, this->values_dofs[comp], temp1);
- apply_tensor_prod<1,true,false> (val, temp1, temp2);
- apply_tensor_prod<2,true,false> (val, temp2, this->hessians_quad[comp][0]);
- }
+ // grad y
+ fe_eval.template apply_values<0,true,false>
+ (fe_eval.values_dofs[c], temp1);
+ if (evaluate_grad == true)
+ {
+ fe_eval.template apply_gradients<1,true,false>
+ (temp1, temp2);
+ fe_eval.template apply_values<2,true,false>
+ (temp2, fe_eval.gradients_quad[c][1]);
+ }
- // grad y
- apply_tensor_prod<0,true,false> (val, this->values_dofs[comp], temp1);
- if (evaluate_grad == true)
- {
- apply_tensor_prod<1,true,false> (grad, temp1, temp2);
- apply_tensor_prod<2,true,false> (val, temp2, this->gradients_quad[comp][1]);
- }
+ if (evaluate_lapl == true)
+ {
+ // grad yz
+ if (evaluate_grad == false)
+ fe_eval.template apply_gradients<1,true,false>
+ (temp1, temp2);
+ fe_eval.template apply_gradients<2,true,false>
+ (temp2, fe_eval.hessians_quad[c][5]);
+
+ // grad yy
+ fe_eval.template apply_hessians<1,true,false>
+ (temp1, temp2);
+ fe_eval.template apply_values<2,true,false>
+ (temp2, fe_eval.hessians_quad[c][1]);
+ }
- if (evaluate_lapl == true)
- {
- // grad yz
- if (evaluate_grad == false)
- apply_tensor_prod<1,true,false> (grad, temp1, temp2);
- apply_tensor_prod<2,true,false> (grad, temp2, this->hessians_quad[comp][5]);
+ // grad z: can use the values applied in x direction stored in temp1
+ fe_eval.template apply_values<1,true,false>
+ (temp1, temp2);
+ if (evaluate_grad == true)
+ fe_eval.template apply_gradients<2,true,false>
+ (temp2, fe_eval.gradients_quad[c][2]);
- // grad yy
- apply_tensor_prod<1,true,false> (hess, temp1, temp2);
- apply_tensor_prod<2,true,false> (val, temp2, this->hessians_quad[comp][1]);
- }
+ // grad zz: can use the values applied in x and y direction stored
+ // in temp2
+ if (evaluate_lapl == true)
+ fe_eval.template apply_hessians<2,true,false>
+ (temp2, fe_eval.hessians_quad[c][2]);
- // grad z: can use the values applied in x direction stored in temp1
- apply_tensor_prod<1,true,false> (val, temp1, temp2);
- if (evaluate_grad == true)
- apply_tensor_prod<2,true,false> (grad, temp2, this->gradients_quad[comp][2]);
+ // val: can use the values applied in x & y direction stored in temp2
+ if (evaluate_val == true)
+ fe_eval.template apply_values<2,true,false>
+ (temp2, fe_eval.values_quad[c]);
- // grad zz: can use the values applied in x and y direction stored in temp2
- if (evaluate_lapl == true)
- apply_tensor_prod<2,true,false> (hess, temp2, this->hessians_quad[comp][2]);
+ break;
- // val: can use the values applied in x & y direction stored in temp2
- if (evaluate_val == true)
- apply_tensor_prod<2,true,false> (val, temp2, this->values_quad[comp]);
- }
- else if (dim == 2)
- {
- // grad x
- if (evaluate_grad == true)
- {
- apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], temp1);
- apply_tensor_prod<1,true,false> (val, temp1, this->gradients_quad[comp][0]);
- }
- if (evaluate_lapl == true)
- {
- // grad xy
- if (evaluate_grad == false)
- apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], temp1);
- apply_tensor_prod<1,true,false> (grad, temp1, this->hessians_quad[comp][2]);
+ case 2:
+
+ // grad x
+ if (evaluate_grad == true)
+ {
+ fe_eval.template apply_gradients<0,true,false>
+ (fe_eval.values_dofs[c], temp1);
+ fe_eval.template apply_values<1,true,false>
+ (temp1, fe_eval.gradients_quad[c][0]);
+ }
+ if (evaluate_lapl == true)
+ {
+ // grad xy
+ if (evaluate_grad == false)
+ fe_eval.template apply_gradients<0,true,false>
+ (fe_eval.values_dofs[c], temp1);
+ fe_eval.template apply_gradients<1,true,false>
+ (temp1, fe_eval.hessians_quad[c][2]);
+
+ // grad xx
+ fe_eval.template apply_hessians<0,true,false>
+ (fe_eval.values_dofs[c], temp1);
+ fe_eval.template apply_values<1,true,false>
+ (temp1, fe_eval.hessians_quad[c][0]);
+ }
+
+ // grad y
+ fe_eval.template apply_values<0,true,false>
+ (fe_eval.values_dofs[c], temp1);
+ if (evaluate_grad == true)
+ fe_eval.template apply_gradients<1,true,false>
+ (temp1, fe_eval.gradients_quad[c][1]);
+
+ // grad yy
+ if (evaluate_lapl == true)
+ fe_eval.template apply_hessians<1,true,false>
+ (temp1, fe_eval.hessians_quad[c][1]);
+
+ // val: can use values applied in x
+ if (evaluate_val == true)
+ fe_eval.template apply_values<1,true,false>
+ (temp1, fe_eval.values_quad[c]);
+
+ break;
+
+ case 1:
+ if (evaluate_val == true)
+ fe_eval.template apply_values<0,true,false>
+ (fe_eval.values_dofs[c], fe_eval.values_quad[c]);
+ if (evaluate_grad == true)
+ fe_eval.template apply_gradients<0,true,false>
+ (fe_eval.values_dofs[c], fe_eval.gradients_quad[c][0]);
+ if (evaluate_lapl == true)
+ fe_eval.template apply_hessians<0,true,false>
+ (fe_eval.values_dofs[c], fe_eval.hessians_quad[c][0]);
+ break;
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ }
+
+#ifdef DEBUG
+ if (evaluate_val == true)
+ fe_eval.values_quad_initialized = true;
+ if (evaluate_grad == true)
+ fe_eval.gradients_quad_initialized = true;
+ if (evaluate_lapl == true)
+ fe_eval.hessians_quad_initialized = true;
+#endif
+ }
+
+
+
+ template <typename FEEval>
+ inline
+ void
+ do_integrate (FEEval &fe_eval,
+ const bool integrate_val,
+ const bool integrate_grad)
+ {
+ Assert (fe_eval.cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+ if (integrate_val == true)
+ Assert (fe_eval.values_quad_submitted == true,
+ ExcAccessToUninitializedField());
+ if (integrate_grad == true)
+ Assert (fe_eval.gradients_quad_submitted == true,
+ ExcAccessToUninitializedField());
+
+ const unsigned int temp_size = fe_eval.dofs_per_cell > fe_eval.n_q_points ?
+ fe_eval.dofs_per_cell : fe_eval.n_q_points;
+ const unsigned int n_components = fe_eval.n_components;
+ const unsigned int dim = fe_eval.dimension;
+
+
+ for(unsigned int c=0; c<n_components; c++)
+ {
+ VectorizedArray<typename FEEval::number_type> temp1[temp_size];
+ VectorizedArray<typename FEEval::number_type> temp2[temp_size];
+
+ switch (dim)
+ {
+ case 3:
+
+ if (integrate_val == true)
+ {
+ // val
+ fe_eval.template apply_values<0,false,false>
+ (fe_eval.values_quad[c], temp1);
+ }
+ if (integrate_grad == true)
+ {
+ // grad x: can sum to temporary value in temp1
+ if (integrate_val == true)
+ fe_eval.template apply_gradients<0,false,true>
+ (fe_eval.gradients_quad[c][0], temp1);
+ else
+ fe_eval.template apply_gradients<0,false,false>
+ (fe_eval.gradients_quad[c][0], temp1);
+ }
+ fe_eval.template apply_values<1,false,false>
+ (temp1, temp2);
+ if (integrate_grad == true)
+ {
+ // grad y: can sum to temporary x value in temp2
+ fe_eval.template apply_values<0,false,false>
+ (fe_eval.gradients_quad[c][1], temp1);
+ fe_eval.template apply_gradients<1,false,true>
+ (temp1, temp2);
+ }
+ fe_eval.template apply_values<2,false,false>
+ (temp2, fe_eval.values_dofs[c]);
+ if (integrate_grad == true)
+ {
+ // grad z: can sum to temporary x and y value in output
+ fe_eval.template apply_values<0,false,false>
+ (fe_eval.gradients_quad[c][2], temp1);
+ fe_eval.template apply_values<1,false,false>
+ (temp1, temp2);
+ fe_eval.template apply_gradients<2,false,true>
+ (temp2, fe_eval.values_dofs[c]);
+ }
+
+ break;
+
+ case 2:
+
+ // val
+ if (integrate_val == true)
+ fe_eval.template apply_values<0,false,false>
+ (fe_eval.values_quad[c], temp1);
+ if (integrate_grad == true)
+ {
+ //grad x
+ if (integrate_val == true)
+ fe_eval.template apply_gradients<0,false,true>
+ (fe_eval.gradients_quad[c][0], temp1);
+ else
+ fe_eval.template apply_gradients<0,false,false>
+ (fe_eval.gradients_quad[c][0], temp1);
+ }
+ fe_eval.template apply_values<1,false,false>
+ (temp1, fe_eval.values_dofs[c]);
+ if (integrate_grad == true)
+ {
+ // grad y
+ fe_eval.template apply_values<0,false,false>
+ (fe_eval.gradients_quad[c][1], temp1);
+ fe_eval.template apply_gradients<1,false,true>
+ (temp1, fe_eval.values_dofs[c]);
+ }
- // grad xx
- apply_tensor_prod<0,true,false> (hess, this->values_dofs[comp], temp1);
- apply_tensor_prod<1,true,false> (val, temp1, this->hessians_quad[comp][0]);
- }
+ break;
- // grad y
- apply_tensor_prod<0,true,false> (val, this->values_dofs[comp], temp1);
- if (evaluate_grad == true)
- apply_tensor_prod<1,true,false> (grad, temp1, this->gradients_quad[comp][1]);
+ case 1:
- // grad yy
- if (evaluate_lapl == true)
- apply_tensor_prod<1,true,false> (hess, temp1, this->hessians_quad[comp][1]);
+ if (integrate_grad == true)
+ fe_eval.template apply_gradients<0,false,false>
+ (fe_eval.gradients_quad[c][0], fe_eval.values_dofs[c]);
+ if (integrate_val == true)
+ {
+ if (integrate_grad == true)
+ fe_eval.template apply_values<0,false,true>
+ (fe_eval.values_quad[c], fe_eval.values_dofs[c]);
+ else
+ fe_eval.template apply_values<0,false,false>
+ (fe_eval.values_quad[c], fe_eval.values_dofs[c]);
+ }
+ break;
- // val: can use values applied in x
- if (evaluate_val == true)
- apply_tensor_prod<1,true,false> (val, temp1, this->values_quad[comp]);
- }
- else if (dim == 1)
- {
- if (evaluate_val == true)
- apply_tensor_prod<0,true,false> (val, this->values_dofs[comp],
- this->values_quad[comp]);
- if (evaluate_grad == true)
- apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp],
- this->gradients_quad[comp][0]);
- if (evaluate_lapl == true)
- apply_tensor_prod<0,true,false> (hess, this->values_dofs[comp],
- this->hessians_quad[comp][0]);
- }
+ default:
+ Assert (false, ExcNotImplemented());
+ }
}
#ifdef DEBUG
- if (evaluate_val == true)
- this->values_quad_initialized = true;
- if (evaluate_grad == true)
- this->gradients_quad_initialized = true;
- if (evaluate_lapl == true)
- this->hessians_quad_initialized = true;
+ fe_eval.dof_values_initialized = true;
#endif
+ }
}
-template <int dim, int fe_degree, int n_q_points_1d, int n_components,
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
typename Number>
inline
void
-FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components,Number>::
-integrate (bool integrate_val,bool integrate_grad)
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::evaluate (const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_lapl)
{
-#ifdef DEBUG
- Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
- if (integrate_val == true)
- Assert (this->values_quad_submitted == true,
- internal::ExcAccessToUninitializedField());
- if (integrate_grad == true)
- Assert (this->gradients_quad_submitted == true,
- internal::ExcAccessToUninitializedField());
-#endif
-
- const VectorizedArray<Number> * val = this->data.shape_values.begin();
- const VectorizedArray<Number> * grad = this->data.shape_gradients.begin();
+ internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl);
+}
- for(unsigned int comp=0;comp<n_components;comp++)
- {
- VectorizedArray<Number> temp1[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
- VectorizedArray<Number> temp2[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
- if (dim == 3)
- {
- if (integrate_val == true)
- {
- // val
- apply_tensor_prod<0,false,false> (val, this->values_quad[comp], temp1);
- }
- if (integrate_grad == true)
- {
- // grad x: can sum to temporary value in temp1
- if (integrate_val == true)
- apply_tensor_prod<0,false,true>
- (grad, this->gradients_quad[comp][0],temp1);
- else
- apply_tensor_prod<0,false,false>
- (grad, this->gradients_quad[comp][0],temp1);
- }
- apply_tensor_prod<1,false,false> (val, temp1, temp2);
- if (integrate_grad == true)
- {
- // grad y: can sum to temporary x value in temp2
- apply_tensor_prod<0,false,false> (val, this->gradients_quad[comp][1], temp1);
- apply_tensor_prod<1,false,true> (grad, temp1, temp2);
- }
- apply_tensor_prod<2,false,false> (val, temp2, this->values_dofs[comp]);
- if (integrate_grad == true)
- {
- // grad z: can sum to temporary x and y value in output
- apply_tensor_prod<0,false,false> (val, this->gradients_quad[comp][2], temp1);
- apply_tensor_prod<1,false,false> (val, temp1, temp2);
- apply_tensor_prod<2,false,true> (grad, temp2, this->values_dofs[comp]);
- }
- }
- else if (dim == 2)
- {
- // val
- if (integrate_val == true)
- apply_tensor_prod<0,false,false> (val, this->values_quad[comp], temp1);
- if (integrate_grad == true)
- {
- //grad x
- if (integrate_val == true)
- apply_tensor_prod<0,false,true>
- (grad, this->gradients_quad[comp][0],temp1);
- else
- apply_tensor_prod<0,false,false>
- (grad, this->gradients_quad[comp][0],temp1);
- }
- apply_tensor_prod<1,false,false> (val, temp1, this->values_dofs[comp]);
- if (integrate_grad == true)
- {
- // grad y
- apply_tensor_prod<0,false,false> (grad, this->gradients_quad[comp][1], temp1);
- apply_tensor_prod<1,false,true> (val, temp1, this->values_dofs[comp]);
- }
- }
- else if (dim == 1)
- {
- if (integrate_grad == true)
- apply_tensor_prod<0,false,false> (grad, this->gradients_quad[comp][0],
- this->values_dofs[comp]);
- if (integrate_val == true)
- {
- if (integrate_grad == true)
- apply_tensor_prod<0,false,true> (val, this->values_quad[comp],
- this->values_dofs[comp]);
- else
- apply_tensor_prod<0,false,false> (val, this->values_quad[comp],
- this->values_dofs[comp]);
- }
- }
- }
-#ifdef DEBUG
- this->dof_values_initialized = true;
-#endif
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+inline
+void
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::integrate (const bool integrate_val,
+ const bool integrate_grad)
+{
+ internal::do_integrate (*this, integrate_val, integrate_grad);
}
-template <int dim, int fe_degree, int n_q_points_1d, int n_components,
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
typename Number>
inline
Point<dim,VectorizedArray<Number> >
-FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components,Number>::
-quadrature_point (const unsigned int q) const
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::quadrature_point (const unsigned int q) const
{
Assert (this->mapping_info.quadrature_points_initialized == true,
ExcNotInitialized());
}
- // General tensor product application for up
- // to three spatial dimensions. Does not
- // assume any symmetry in the shape values
- // field
-template <int dim, int fe_degree, int n_q_points_1d, int n_components,
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
typename Number>
template <int direction, bool dof_to_quad, bool add>
inline
void
-FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components,Number>::
-apply_tensor_prod (const VectorizedArray<Number>*shape_data,
- const VectorizedArray<Number> input [],
- VectorizedArray<Number> output [])
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::apply_values(const VectorizedArray<Number> in [],
+ VectorizedArray<Number> out [])
{
- AssertIndexRange (direction, dim);
- const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
- nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
+ internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,Number,
+ direction, dof_to_quad, add>
+ (this->data.shape_values.begin(), in, out);
+}
- const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
- const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
- const int stride = ((direction > 0 ? nn : 1 ) *
- (direction > 1 ? nn : 1));
- const VectorizedArray<Number> * in = &input[0];
- VectorizedArray<Number> * out = &output[0];
- for (int i2=0; i2<n_blocks2; ++i2)
- {
- for (int i1=0; i1<n_blocks1; ++i1)
- {
- for (int col=0; col<nn; ++col)
- {
- VectorizedArray<Number> val0;
- if (dof_to_quad == true)
- val0 = shape_data[col];
- else
- val0 = shape_data[col*n_q_points_1d];
- VectorizedArray<Number> res0 = val0 * in[0];
- for (int ind=1; ind<mm; ++ind)
- {
- if (dof_to_quad == true)
- val0 = shape_data[ind*n_q_points_1d+col];
- else
- val0 = shape_data[col*n_q_points_1d+ind];
- res0 += val0 * in[stride*ind];
- }
- if (add == false)
- out[stride*col] = res0;
- else
- out[stride*col] += res0;
- }
- // increment: in regular case, just go to the
- // next point in x-direction. If we are at the
- // end of one chunk in x-dir, need to jump
- // over to the next layer in z-direction
- switch (direction)
- {
- case 0:
- in += mm;
- out += nn;
- break;
- case 1:
- case 2:
- ++in;
- ++out;
- break;
- default:
- Assert (false, ExcNotImplemented());
- }
- }
- if (direction == 1)
- {
- in += nn*(mm-1);
- out += nn*(nn-1);
- }
- }
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+template <int direction, bool dof_to_quad, bool add>
+inline
+void
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::apply_gradients(const VectorizedArray<Number> in [],
+ VectorizedArray<Number> out [])
+{
+ internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,Number,
+ direction, dof_to_quad, add>
+ (this->data.shape_gradients.begin(), in, out);
+}
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+template <int direction, bool dof_to_quad, bool add>
+inline
+void
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::apply_hessians(const VectorizedArray<Number> in [],
+ VectorizedArray<Number> out [])
+{
+ internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,Number,
+ direction, dof_to_quad, add>
+ (this->data.shape_hessians.begin(), in, out);
}
-/*----------------------- FEEvaluation -------------------------------*/
+/*-------------------------- FEEvaluation -----------------------------------*/
-template <int dim, int fe_degree, int n_q_points_1d, int n_components,
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
typename Number>
inline
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number>::
-FEEvaluation (const MatrixFree<dim,Number> &data_in,
- const unsigned int fe_no,
- const unsigned int quad_no)
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluation (const MatrixFree<dim,Number> &data_in,
+ const unsigned int fe_no,
+ const unsigned int quad_no)
:
BaseClass (data_in, fe_no, quad_no)
{
-template <int dim, int fe_degree, int n_q_points_1d, int n_components,
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
typename Number>
inline
void
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number>::
-evaluate (bool evaluate_val, bool evaluate_grad, bool evaluate_lapl)
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::evaluate (const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_lapl)
{
- Assert (this->cell != numbers::invalid_unsigned_int,
- ExcNotInitialized());
- Assert (this->dof_values_initialized == true,
- internal::ExcAccessToUninitializedField());
-
- for(unsigned int comp=0;comp<n_components;comp++)
- {
- VectorizedArray<Number> temp1[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
- VectorizedArray<Number> temp2[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
-
- if (dim == 3)
- {
- if (evaluate_grad == true)
- {
- // grad x
- apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
- apply_values<1,true,false> (temp1, temp2);
- apply_values<2,true,false> (temp2, this->gradients_quad[comp][0]);
- }
-
- if (evaluate_lapl == true)
- {
- // grad xz
- if (evaluate_grad == false)
- {
- apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
- apply_values<1,true,false> (temp1, temp2);
- }
- apply_gradients<2,true,false>(temp2, this->hessians_quad[comp][4]);
-
- // grad xy
- apply_gradients<1,true,false>(temp1, temp2);
- apply_values<2,true,false> (temp2, this->hessians_quad[comp][3]);
-
- // grad xx
- apply_hessians<0,true,false> (this->values_dofs[comp], temp1);
- apply_values<1,true,false> (temp1, temp2);
- apply_values<2,true,false> (temp2, this->hessians_quad[comp][0]);
- }
-
- // grad y
- apply_values<0,true,false> (this->values_dofs[comp], temp1);
- if (evaluate_grad == true)
- {
- apply_gradients<1,true,false> (temp1, temp2);
- apply_values<2,true,false> (temp2, this->gradients_quad[comp][1]);
- }
-
- if (evaluate_lapl == true)
- {
- // grad yz
- if (evaluate_grad == false)
- apply_gradients<1,true,false> (temp1, temp2);
- apply_gradients<2,true,false> (temp2, this->hessians_quad[comp][5]);
-
- // grad yy
- apply_hessians<1,true,false> (temp1, temp2);
- apply_values<2,true,false> (temp2, this->hessians_quad[comp][1]);
- }
-
- // grad z: can use the values applied in x direction stored in temp1
- apply_values<1,true,false> (temp1, temp2);
- if (evaluate_grad == true)
- apply_gradients<2,true,false> (temp2, this->gradients_quad[comp][2]);
-
- // grad zz: can use the values applied in x and y direction stored in temp2
- if (evaluate_lapl == true)
- apply_hessians<2,true,false> (temp2, this->hessians_quad[comp][2]);
-
- // val: can use the values applied in x & y direction stored in temp2
- if (evaluate_val == true)
- apply_values<2,true,false> (temp2, this->values_quad[comp]);
- }
- else if (dim == 2)
- {
- // grad x
- if (evaluate_grad == true)
- {
- apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
- apply_values<1,true,false> (temp1, this->gradients_quad[comp][0]);
- }
- if (evaluate_lapl == true)
- {
- // grad xy
- if (evaluate_grad == false)
- apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
- apply_gradients<1,true,false> (temp1, this->hessians_quad[comp][2]);
-
- // grad xx
- apply_hessians<0,true,false> (this->values_dofs[comp], temp1);
- apply_values<1,true,false> (temp1, this->hessians_quad[comp][0]);
- }
-
- // grad y
- apply_values<0,true,false> (this->values_dofs[comp], temp1);
- if (evaluate_grad == true)
- apply_gradients<1,true,false> (temp1, this->gradients_quad[comp][1]);
-
- // grad yy
- if (evaluate_lapl == true)
- apply_hessians<1,true,false> (temp1, this->hessians_quad[comp][1]);
-
- // val: can use values applied in x
- if (evaluate_val == true)
- apply_values<1,true,false> (temp1, this->values_quad[comp]);
- }
- else if (dim == 1)
- {
- if (evaluate_val == true)
- apply_values<0,true,false> (this->values_dofs[comp],
- this->values_quad[comp]);
- if (evaluate_grad == true)
- apply_gradients<0,true,false> (this->values_dofs[comp],
- this->gradients_quad[comp][0]);
- if (evaluate_lapl == true)
- apply_hessians<0,true,false> (this->values_dofs[comp],
- this->hessians_quad[comp][0]);
- }
- }
-
-#ifdef DEBUG
- if (evaluate_val == true)
- this->values_quad_initialized = true;
- if (evaluate_grad == true)
- this->gradients_quad_initialized = true;
- if (evaluate_lapl == true)
- this->hessians_quad_initialized = true;
-#endif
+ internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl);
}
-template <int dim, int fe_degree, int n_q_points_1d, int n_components,
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
typename Number>
inline
void
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number>::
-integrate (bool integrate_val,bool integrate_grad)
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::integrate (bool integrate_val,bool integrate_grad)
{
-#ifdef DEBUG
- Assert (this->cell != numbers::invalid_unsigned_int,
- ExcNotInitialized());
- if (integrate_val == true)
- Assert (this->values_quad_submitted == true,
- internal::ExcAccessToUninitializedField());
- if (integrate_grad == true)
- Assert (this->gradients_quad_submitted == true,
- internal::ExcAccessToUninitializedField());
-#endif
-
- for(unsigned int comp=0;comp<n_components;comp++)
- {
- VectorizedArray<Number> temp1[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
- VectorizedArray<Number> temp2[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points];
-
- if (dim == 3)
- {
- if (integrate_val == true)
- {
- // val
- apply_values<0,false,false> (this->values_quad[comp], temp1);
- }
- if (integrate_grad == true)
- {
- // grad x: can sum to temporary value in temp1
- if (integrate_val == true)
- apply_gradients<0,false,true> (this->gradients_quad[comp][0],
- temp1);
- else
- apply_gradients<0,false,false> (this->gradients_quad[comp][0],
- temp1);
- }
- apply_values<1,false,false> (temp1, temp2);
- if (integrate_grad == true)
- {
- // grad y: can sum to temporary x value in temp2
- apply_values<0,false,false> (this->gradients_quad[comp][1], temp1);
- apply_gradients<1,false,true> (temp1, temp2);
- }
- apply_values<2,false,false> (temp2, this->values_dofs[comp]);
- if (integrate_grad == true)
- {
- // grad z: can sum to temporary x and y value in output
- apply_values<0,false,false> (this->gradients_quad[comp][2], temp1);
- apply_values<1,false,false> (temp1, temp2);
- apply_gradients<2,false,true> (temp2, this->values_dofs[comp]);
- }
- }
- else if (dim == 2)
- {
- // val
- if (integrate_val == true)
- apply_values<0,false,false> (this->values_quad[comp], temp1);
- if (integrate_grad == true)
- {
- //grad x
- if (integrate_val == true)
- apply_gradients<0,false,true> (this->gradients_quad[comp][0],
- temp1);
- else
- apply_gradients<0,false,false> (this->gradients_quad[comp][0],
- temp1);
- }
- apply_values<1,false,false> (temp1, this->values_dofs[comp]);
- if (integrate_grad == true)
- {
- // grad y
- apply_values<0,false,false> (this->gradients_quad[comp][1], temp1);
- apply_gradients<1,false,true> (temp1, this->values_dofs[comp]);
- }
- }
- else if (dim == 1)
- {
- if (integrate_grad == true)
- apply_gradients<0,false,false> (this->gradients_quad[comp][0],
- this->values_dofs[comp]);
- if (integrate_val == true)
- {
- if (integrate_grad == true)
- apply_values<0,false,true> (this->values_quad[comp],
- this->values_dofs[comp]);
- else
- apply_values<0,false,false> (this->values_quad[comp],
- this->values_dofs[comp]);
- }
- }
- }
-#ifdef DEBUG
- this->dof_values_initialized = true;
-#endif
+ internal::do_integrate (*this, integrate_val, integrate_grad);
}
-// ----------------- optimized implementation tensor product symmetric case
+/*----------------- optimized implementation tensor product symmetric case --*/
-template <int dim, int fe_degree, int n_q_points_1d, int n_components,
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
typename Number>
template <int direction, bool dof_to_quad, bool add>
inline
void
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number>::
-apply_values (const VectorizedArray<Number> input [],
- VectorizedArray<Number> output [])
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::apply_values (const VectorizedArray<Number> in [],
+ VectorizedArray<Number> out [])
{
AssertIndexRange (direction, dim);
const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
- const int stride = ((direction > 0 ? nn : 1 ) *
- (direction > 1 ? nn : 1));
-
- const VectorizedArray<Number> * in = &input[0];
- VectorizedArray<Number> * out = &output[0];
+ const int stride = Utilities::fixed_int_power<nn,direction>::value;
+
+ // This loop specializes the general
+ // application of tensor-product based
+ // elements for "symmetric" finite elements,
+ // i.e., when the shape functions are
+ // symmetric about 0.5 and the quadrature
+ // points are, too. In that case, the 1D shape
+ // values read (sorted lexicographically, rows
+ // run over 1D dofs, columns over quadrature
+ // points):
+ // Q2 --> [ 0.687 0 -0.087 ]
+ // [ 0.4 1 0.4 ]
+ // [-0.087 0 0.687 ]
+ // Q3 --> [ 0.66 0.003 0.002 0.049 ]
+ // [ 0.521 1.005 -0.01 -0.230 ]
+ // [-0.230 -0.01 1.005 0.521 ]
+ // [ 0.049 0.002 0.003 0.66 ]
+ // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
+ // [ 0.608 1.059 0 0.039 0.176 ]
+ // [-0.409 -0.113 1 -0.113 -0.409 ]
+ // [ 0.176 0.039 0 1.059 0.608 ]
+ // [-0.032 -0.007 0 0.022 0.658 ]
+ //
+ // In these matrices, we want to use avoid
+ // computations involving zeros and ones and
+ // in addition use the symmetry in entries to
+ // reduce the number of read operations.
+ const VectorizedArray<Number> * shape_values = this->data.shape_values.begin();
for (int i2=0; i2<n_blocks2; ++i2)
{
for (int i1=0; i1<n_blocks1; ++i1)
{
for (int col=0; col<n_cols; ++col)
{
- VectorizedArray<Number> val0, val1, res0, res1;
+ VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
if (dof_to_quad == true)
{
- val0 = this->data.shape_values[col];
- val1 = this->data.shape_values[nn-1-col];
+ val0 = shape_values[col];
+ val1 = shape_values[nn-1-col];
}
else
{
- val0 = this->data.shape_values[col*n_q_points_1d];
- val1 = this->data.shape_values[(col+1)*n_q_points_1d-1];
+ val0 = shape_values[col*n_q_points_1d];
+ val1 = shape_values[(col+1)*n_q_points_1d-1];
}
if (mid > 0)
{
- res0 = val0 * in[0];
- res1 = val1 * in[0];
- res0 += val1 * in[stride*(mm-1)];
- res1 += val0 * in[stride*(mm-1)];
+ in0 = in[0];
+ in1 = in[stride*(mm-1)];
+ res0 = val0 * in0;
+ res1 = val1 * in0;
+ res0 += val1 * in1;
+ res1 += val0 * in1;
for (int ind=1; ind<mid; ++ind)
{
if (dof_to_quad == true)
{
- val0 = this->data.shape_values[ind*n_q_points_1d+col];
- val1 = this->data.shape_values[ind*n_q_points_1d+nn-1-col];
+ val0 = shape_values[ind*n_q_points_1d+col];
+ val1 = shape_values[ind*n_q_points_1d+nn-1-col];
}
else
{
- val0 = this->data.shape_values[col*n_q_points_1d+ind];
- val1 = this->data.shape_values[(col+1)*n_q_points_1d-1-ind];
+ val0 = shape_values[col*n_q_points_1d+ind];
+ val1 = shape_values[(col+1)*n_q_points_1d-1-ind];
}
- res0 += val0 * in[stride*ind];
- res1 += val1 * in[stride*ind];
- res0 += val1 * in[stride*(mm-1-ind)];
- res1 += val0 * in[stride*(mm-1-ind)];
+ in0 = in[stride*ind];
+ in1 = in[stride*(mm-1-ind)];
+ res0 += val0 * in0;
+ res1 += val1 * in0;
+ res0 += val1 * in1;
+ res1 += val0 * in1;
}
}
else
{
if (mm % 2 == 1)
{
- val0 = this->data.shape_values[mid*n_q_points_1d+col];
+ val0 = shape_values[mid*n_q_points_1d+col];
val1 = val0 * in[stride*mid];
res0 += val1;
res1 += val1;
{
if (mm % 2 == 1 && nn % 2 == 0)
{
- val0 = this->data.shape_values[col*n_q_points_1d+mid];
+ val0 = shape_values[col*n_q_points_1d+mid];
val1 = val0 * in[stride*mid];
res0 += val1;
res1 += val1;
else if (dof_to_quad == true && nn%2==1)
{
VectorizedArray<Number> res0;
- VectorizedArray<Number> val0 = this->data.shape_values[n_cols];
+ VectorizedArray<Number> val0 = shape_values[n_cols];
if (mid > 0)
{
res0 = in[0] + in[stride*(mm-1)];
res0 *= val0;
for (int ind=1; ind<mid; ++ind)
{
- val0 = this->data.shape_values[ind*n_q_points_1d+n_cols];
+ val0 = shape_values[ind*n_q_points_1d+n_cols];
VectorizedArray<Number> val1 = in[stride*ind] + in[stride*(mm-1-ind)];
val1 *= val0;
res0 += val1;
res0 = VectorizedArray<Number>();
if (mm % 2 == 1)
{
- val0 = this->data.shape_values[mid*n_q_points_1d+n_cols];
+ val0 = shape_values[mid*n_q_points_1d+n_cols];
res0 += val0 * in[stride*mid];
}
if (add == false)
VectorizedArray<Number> res0;
if (mid > 0)
{
- VectorizedArray<Number> val0 = this->data.shape_values[n_cols*n_q_points_1d];
+ VectorizedArray<Number> val0 = shape_values[n_cols*n_q_points_1d];
res0 = in[0] + in[stride*(mm-1)];
res0 *= val0;
for (int ind=1; ind<mid; ++ind)
{
- val0 = this->data.shape_values[n_cols*n_q_points_1d+ind];
+ val0 = shape_values[n_cols*n_q_points_1d+ind];
VectorizedArray<Number> val1 = in[stride*ind] + in[stride*(mm-1-ind)];
val1 *= val0;
res0 += val1;
-template <int dim, int fe_degree, int n_q_points_1d, int n_components,
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
typename Number>
template <int direction, bool dof_to_quad, bool add>
inline
void
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number>::
-apply_gradients (const VectorizedArray<Number> input [],
- VectorizedArray<Number> output [])
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::apply_gradients (const VectorizedArray<Number> in [],
+ VectorizedArray<Number> out [])
{
AssertIndexRange (direction, dim);
const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
- const int stride = ((direction > 0 ? nn : 1 ) *
- (direction > 1 ? nn : 1));
+ const int stride = Utilities::fixed_int_power<nn,direction>::value;
- const VectorizedArray<Number> * in = &input[0];
- VectorizedArray<Number> * out = &output[0];
+ const VectorizedArray<Number> * shape_gradients = this->data.shape_gradients.begin();
for (int i2=0; i2<n_blocks2; ++i2)
{
for (int i1=0; i1<n_blocks1; ++i1)
{
+ // For the specialized loop used for the
+ // gradient computation in here, the 1D shape
+ // values read (sorted lexicographically, rows
+ // run over 1D dofs, columns over quadrature
+ // points):
+ // Q2 --> [-2.549 -1 0.549 ]
+ // [ 3.098 0 -3.098 ]
+ // [-0.549 1 2.549 ]
+ // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
+ // [ 6.07 -1.44 -2.97 2.196 ]
+ // [-2.196 2.97 1.44 -6.07 ]
+ // [ 0.44 -0.5 1.03 4.315 ]
+ // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
+ // [10.111 -2.76 -2.667 2.066 -2.306 ]
+ // [-5.688 5.773 0 -5.773 5.688 ]
+ // [ 2.306 -2.066 2.667 2.76 -10.111 ]
+ // [-0.413 0.353 -0.333 -0.353 0.413 ]
+ //
+ // In these matrices, we want to use avoid
+ // computations involving zeros and ones and
+ // in addition use the symmetry in entries to
+ // reduce the number of read operations.
for (int col=0; col<n_cols; ++col)
{
- VectorizedArray<Number> val0, val1, res0, res1;
+ VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
if (dof_to_quad == true)
{
- val0 = this->data.shape_gradients[col];
- val1 = this->data.shape_gradients[nn-1-col];
+ val0 = shape_gradients[col];
+ val1 = shape_gradients[nn-1-col];
}
else
{
- val0 = this->data.shape_gradients[col*n_q_points_1d];
- val1 = this->data.shape_gradients[(nn-col-1)*n_q_points_1d];
+ val0 = shape_gradients[col*n_q_points_1d];
+ val1 = shape_gradients[(nn-col-1)*n_q_points_1d];
}
if (mid > 0)
{
- res0 = val0 * in[0];
- res1 = val1 * in[0];
- res0 -= val1 * in[stride*(mm-1)];
- res1 -= val0 * in[stride*(mm-1)];
+ in0 = in[0];
+ in1 = in[stride*(mm-1)];
+ res0 = val0 * in0;
+ res1 = val1 * in0;
+ res0 -= val1 * in1;
+ res1 -= val0 * in1;
for (int ind=1; ind<mid; ++ind)
{
if (dof_to_quad == true)
{
- val0 = this->data.shape_gradients[ind*n_q_points_1d+col];
- val1 = this->data.shape_gradients[ind*n_q_points_1d+nn-1-col];
+ val0 = shape_gradients[ind*n_q_points_1d+col];
+ val1 = shape_gradients[ind*n_q_points_1d+nn-1-col];
}
else
{
- val0 = this->data.shape_gradients[col*n_q_points_1d+ind];
- val1 = this->data.shape_gradients[(nn-col-1)*n_q_points_1d+ind];
+ val0 = shape_gradients[col*n_q_points_1d+ind];
+ val1 = shape_gradients[(nn-col-1)*n_q_points_1d+ind];
}
- res0 += val0 * in[stride*ind];
- res1 += val1 * in[stride*ind];
- res0 -= val1 * in[stride*(mm-1-ind)];
- res1 -= val0 * in[stride*(mm-1-ind)];
+ in0 = in[stride*ind];
+ in1 = in[stride*(mm-1-ind)];
+ res0 += val0 * in0;
+ res1 += val1 * in0;
+ res0 -= val1 * in1;
+ res1 -= val0 * in1;
}
}
else
if (mm % 2 == 1)
{
if (dof_to_quad == true)
- val0 = this->data.shape_gradients[mid*n_q_points_1d+col];
+ val0 = shape_gradients[mid*n_q_points_1d+col];
else
- val0 = this->data.shape_gradients[col*n_q_points_1d+mid];
+ val0 = shape_gradients[col*n_q_points_1d+mid];
val1 = val0 * in[stride*mid];
res0 += val1;
res1 -= val1;
{
VectorizedArray<Number> val0, res0;
if (dof_to_quad == true)
- val0 = this->data.shape_gradients[n_cols];
+ val0 = shape_gradients[n_cols];
else
- val0 = this->data.shape_gradients[n_cols*n_q_points_1d];
+ val0 = shape_gradients[n_cols*n_q_points_1d];
res0 = in[0] - in[stride*(mm-1)];
res0 *= val0;
for (int ind=1; ind<mid; ++ind)
{
if (dof_to_quad == true)
- val0 = this->data.shape_gradients[ind*n_q_points_1d+n_cols];
+ val0 = shape_gradients[ind*n_q_points_1d+n_cols];
else
- val0 = this->data.shape_gradients[n_cols*n_q_points_1d+ind];
+ val0 = shape_gradients[n_cols*n_q_points_1d+ind];
VectorizedArray<Number> val1 = in[stride*ind] - in[stride*(mm-1-ind)];
val1 *= val0;
res0 += val1;
// same symmetry relations hold. However, it
// is not possible to omit some values that
// are zero for the values
-template <int dim, int fe_degree, int n_q_points_1d, int n_components,
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
typename Number>
template <int direction, bool dof_to_quad, bool add>
inline
void
-FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number>::
-apply_hessians (const VectorizedArray<Number> input [],
- VectorizedArray<Number> output [])
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::apply_hessians (const VectorizedArray<Number> in [],
+ VectorizedArray<Number> out [])
{
AssertIndexRange (direction, dim);
const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
- const int stride = ((direction > 0 ? nn : 1 ) *
- (direction > 1 ? nn : 1));
+ const int stride = Utilities::fixed_int_power<nn,direction>::value;
- const VectorizedArray<Number> * in = &input[0];
- VectorizedArray<Number> * out = &output[0];
for (int i2=0; i2<n_blocks2; ++i2)
{
for (int i1=0; i1<n_blocks1; ++i1)
{
for (int col=0; col<n_cols; ++col)
{
- VectorizedArray<Number> val0, val1, res0, res1;
+ VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
if (dof_to_quad == true)
{
val0 = this->data.shape_hessians[col];
}
if (mid > 0)
{
- res0 = val0 * in[0];
- res1 = val1 * in[0];
- res0 += val1 * in[stride*(mm-1)];
- res1 += val0 * in[stride*(mm-1)];
+ in0 = in[0];
+ in1 = in[stride*(mm-1)];
+ res0 = val0 * in0;
+ res1 = val1 * in0;
+ res0 += val1 * in1;
+ res1 += val0 * in1;
for (int ind=1; ind<mid; ++ind)
{
if (dof_to_quad == true)
val0 = this->data.shape_hessians[col*n_q_points_1d+ind];
val1 = this->data.shape_hessians[(col+1)*n_q_points_1d-1-ind];
}
- res0 += val0 * in[stride*ind];
- res1 += val1 * in[stride*ind];
- res0 += val1 * in[stride*(mm-1-ind)];
- res1 += val0 * in[stride*(mm-1-ind)];
+ in0 = in[stride*ind];
+ in1 = in[stride*(mm-1-ind)];
+ res0 += val0 * in0;
+ res1 += val1 * in0;
+ res0 += val1 * in1;
+ res1 += val0 * in1;
}
}
else
}
-/*----------------------- FEEvaluationGL -------------------------------*/
+/*------------------------- FEEvaluationGL ----------------------------------*/
-template <int dim, int fe_degree, int n_components, typename Number>
+template <int dim, int fe_degree, int n_components_, typename Number>
inline
-FEEvaluationGL<dim,fe_degree,n_components,Number>::
-FEEvaluationGL (const MatrixFree<dim,Number> &data_in,
+FEEvaluationGL<dim,fe_degree,n_components_,Number>
+::FEEvaluationGL (const MatrixFree<dim,Number> &data_in,
const unsigned int fe_no,
const unsigned int quad_no)
:
std::string error_mess = "FEEvaluationGL not appropriate. It assumes:\n";
error_mess += " - identity operation for shape values\n";
error_mess += " - zero diagonal at interior points for gradients\n";
+ error_mess += " - gradient equal to unity at element boundary\n";
error_mess += "Try FEEvaluation<...> instead!";
const double zero_tol =
for (unsigned int i=1; i<n_points_1d-1; ++i)
Assert (std::fabs(this->data.shape_gradients[i*n_points_1d+i][0])<zero_tol,
ExcMessage (error_mess.c_str()));
+ Assert (std::fabs(this->data.shape_gradients[n_points_1d-1][0]-
+ (n_points_1d%2==0 ? -1. : 1.)) < zero_tol,
+ ExcMessage (error_mess.c_str()));
#endif
}
-template <int dim, int fe_degree, int n_components, typename Number>
+template <int dim, int fe_degree, int n_components_, typename Number>
inline
void
-FEEvaluationGL<dim,fe_degree,n_components,Number>::
-evaluate (bool evaluate_val,bool evaluate_grad,bool evaluate_lapl)
+FEEvaluationGL<dim,fe_degree,n_components_,Number>
+::evaluate (const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_lapl)
{
Assert (this->cell != numbers::invalid_unsigned_int,
ExcNotInitialized());
this->values_quad_initialized = true;
#endif
}
+ // separate implementation here compared to
+ // the general case because the values are an
+ // identity operation
if (evaluate_grad == true)
{
for(unsigned int comp=0;comp<n_components;comp++)
-template <int dim, int fe_degree, int n_components, typename Number>
+template <int dim, int fe_degree, int n_components_, typename Number>
inline
void
-FEEvaluationGL<dim,fe_degree,n_components,Number>::
-integrate (bool integrate_val, bool integrate_grad)
+FEEvaluationGL<dim,fe_degree,n_components_,Number>
+::integrate (const bool integrate_val, const bool integrate_grad)
{
Assert (this->cell != numbers::invalid_unsigned_int,
ExcNotInitialized());
-template <int dim, int fe_degree, int n_components, typename Number>
+template <int dim, int fe_degree, int n_components_, typename Number>
template <int direction, bool dof_to_quad, bool add>
inline
void
-FEEvaluationGL<dim,fe_degree,n_components,Number>::
-apply_gradients (const VectorizedArray<Number> input [],
- VectorizedArray<Number> output [])
+FEEvaluationGL<dim,fe_degree,n_components_,Number>
+::apply_gradients (const VectorizedArray<Number> in [],
+ VectorizedArray<Number> out [])
{
AssertIndexRange (direction, dim);
const int mm = fe_degree+1;
const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
- const int stride = ((direction > 0 ? nn : 1 ) *
- (direction > 1 ? nn : 1));
-
- const VectorizedArray<Number> * in = &input[0];
- VectorizedArray<Number> * out = &output[0];
+ const int stride = Utilities::fixed_int_power<nn,direction>::value;
+
+ // This loop specializes the application of
+ // the tensor product loop for Gauss-Lobatto
+ // elements which are symmetric about 0.5 just
+ // as the general class of elements treated by
+ // FEEvaluation, have diagonal shape matrices
+ // for the values and have the following
+ // gradient matrices (notice the zeros on the
+ // diagonal in the interior points, which is
+ // due to the construction of Legendre
+ // polynomials):
+ // Q2 --> [-3 -1 1 ]
+ // [ 4 0 -4 ]
+ // [-1 1 3 ]
+ // Q3 --> [-6 -1.618 0.618 -1 ]
+ // [ 8.09 0 -2.236 3.09 ]
+ // [-3.09 2.236 0 -8.09 ]
+ // [ 1 -0.618 1.618 6 ]
+ // Q4 --> [-10 -2.482 0.75 -0.518 1 ]
+ // [ 13.51 0 -2.673 1.528 -2.82 ]
+ // [-5.333 3.491 0 -3.491 5.333 ]
+ // [ 2.82 -1.528 2.673 0 -13.51 ]
+ // [-1 0.518 -0.75 2.482 10 ]
for (int i2=0; i2<n_blocks2; ++i2)
{
for (int i1=0; i1<n_blocks1; ++i1)
{
for (int col=0; col<n_cols; ++col)
{
- VectorizedArray<Number> val0, val1, res0, res1;
- if (dof_to_quad == true)
- {
- val0 = this->data.shape_gradients[col];
- val1 = this->data.shape_gradients[nn-1-col];
- }
- else
- {
- val0 = this->data.shape_gradients[col*mm];
- val1 = this->data.shape_gradients[(nn-col-1)*mm];
- }
+ VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
if (mid > 0)
{
- res0 = val0 * in[0];
- res1 = val1 * in[0];
- res0 -= val1 * in[stride*(mm-1)];
- res1 -= val0 * in[stride*(mm-1)];
+ if (dof_to_quad == true)
+ {
+ val0 = this->data.shape_gradients[col];
+ val1 = this->data.shape_gradients[nn-1-col];
+ }
+ else
+ {
+ val0 = this->data.shape_gradients[col*mm];
+ val1 = this->data.shape_gradients[(nn-col-1)*mm];
+ }
+ in0 = in[0];
+ in1 = in[stride*(mm-1)];
+ if (col == 0)
+ {
+ if ((mm+dof_to_quad)%2 == 1)
+ {
+ res0 = val0 * in0;
+ res1 = -in0;
+ res0 += in1;
+ res1 -= val0 * in1;
+ }
+ else
+ {
+ res0 = val0 * in0;
+ res0 -= in1;
+ res1 = in0;
+ res1 -= val0 * in1;
+ }
+ }
+ else
+ {
+ res0 = val0 * in0;
+ res1 = val1 * in0;
+ res0 -= val1 * in1;
+ res1 -= val0 * in1;
+ }
for (int ind=1; ind<mid; ++ind)
{
if (dof_to_quad == true)
// at inner points, the gradient is zero for
// ind==col
+ in0 = in[stride*ind];
+ in1 = in[stride*(mm-1-ind)];
if (ind == col)
{
- res1 += val1 * in[stride*ind];
- res0 -= val1 * in[stride*(mm-1-ind)];
+ res1 += val1 * in0;
+ res0 -= val1 * in1;
}
else
{
- res0 += val0 * in[stride*ind];
- res1 += val1 * in[stride*ind];
- res0 -= val1 * in[stride*(mm-1-ind)];
- res1 -= val0 * in[stride*(mm-1-ind)];
+ res0 += val0 * in0;
+ res1 += val1 * in0;
+ res0 -= val1 * in1;
+ res1 -= val0 * in1;
}
}
}