Make class templated on generic type Number that can be VectorizedArray<double>.
unsigned int size_shape_values = n_dofs_1d*n_q_points_1d*sizeof(Number);
cudaError_t cuda_error = cudaMemcpyToSymbol(internal::global_shape_values,
- &shape_info.shape_values_number[0],
+ &shape_info.shape_values[0],
size_shape_values,
0,
cudaMemcpyHostToDevice);
if (update_flags & update_gradients)
{
cuda_error = cudaMemcpyToSymbol(internal::global_shape_gradients,
- &shape_info.shape_gradient_number[0],
+ &shape_info.shape_gradient[0],
size_shape_values,
0,
cudaMemcpyHostToDevice);
struct FEEvaluationImpl
{
static
- void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ void evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
VectorizedArray<Number> *values_dofs_actual[],
VectorizedArray<Number> *values_quad[],
VectorizedArray<Number> *gradients_quad[][dim],
const bool evaluate_hessians);
static
- void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ void integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
VectorizedArray<Number> *values_dofs_actual[],
VectorizedArray<Number> *values_quad[],
VectorizedArray<Number> *gradients_quad[][dim],
inline
void
FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
- ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ ::evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
VectorizedArray<Number> *values_dofs_actual[],
VectorizedArray<Number> *values_quad[],
VectorizedArray<Number> *gradients_quad[][dim],
inline
void
FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
- ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ ::integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
VectorizedArray<Number> *values_dofs_actual[],
VectorizedArray<Number> *values_quad[],
VectorizedArray<Number> *gradients_quad[][dim],
struct FEEvaluationImplTransformToCollocation
{
static
- void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ void evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
VectorizedArray<Number> *values_dofs[],
VectorizedArray<Number> *values_quad[],
VectorizedArray<Number> *gradients_quad[][dim],
const bool evaluate_hessians);
static
- void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ void integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
VectorizedArray<Number> *values_dofs[],
VectorizedArray<Number> *values_quad[],
VectorizedArray<Number> *gradients_quad[][dim],
inline
void
FEEvaluationImplTransformToCollocation<dim, fe_degree, n_components, Number>
- ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ ::evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
VectorizedArray<Number> *values_dofs[],
VectorizedArray<Number> *values_quad[],
VectorizedArray<Number> *gradients_quad[][dim],
inline
void
FEEvaluationImplTransformToCollocation<dim, fe_degree, n_components, Number>
- ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ ::integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
VectorizedArray<Number> *values_dofs[],
VectorizedArray<Number> *values_quad[],
VectorizedArray<Number> *gradients_quad[][dim],
struct FEEvaluationImplCollocation
{
static
- void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ void evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
VectorizedArray<Number> *values_dofs[],
VectorizedArray<Number> *values_quad[],
VectorizedArray<Number> *gradients_quad[][dim],
const bool evaluate_hessians);
static
- void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ void integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
VectorizedArray<Number> *values_dofs[],
VectorizedArray<Number> *values_quad[],
VectorizedArray<Number> *gradients_quad[][dim],
inline
void
FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
- ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ ::evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
VectorizedArray<Number> *values_dofs[],
VectorizedArray<Number> *values_quad[],
VectorizedArray<Number> *gradients_quad[][dim],
inline
void
FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
- ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ ::integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
VectorizedArray<Number> *values_dofs[],
VectorizedArray<Number> *values_quad[],
VectorizedArray<Number> *gradients_quad[][dim],
/**
* Return a reference to the ShapeInfo object currently in use.
*/
- const internal::MatrixFreeFunctions::ShapeInfo<Number> &
- get_shape_info() const;
+ const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
+ get_shape_info() const;
/**
* Fills the JxW values currently used.
* product. Also contained in matrix_info, but it simplifies code if we
* store a reference to it.
*/
- const internal::MatrixFreeFunctions::ShapeInfo<Number> *data;
+ const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> *data;
/**
* A pointer to the Cartesian Jacobian information of the present cell. Only
dof_info (nullptr),
mapping_info (nullptr),
// select the correct base element from the given FE component
- data (new internal::MatrixFreeFunctions::ShapeInfo<Number>(quadrature, fe, fe.component_to_base_index(first_selected_component).first)),
+ data (new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(quadrature, fe, fe.component_to_base_index(first_selected_component).first)),
cartesian_data (nullptr),
jacobian (nullptr),
J_value (nullptr),
dof_info (other.dof_info),
mapping_info (other.mapping_info),
data (other.matrix_info == nullptr ?
- new internal::MatrixFreeFunctions::ShapeInfo<Number>(*other.data) :
+ new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(*other.data) :
other.data),
cartesian_data (nullptr),
jacobian (nullptr),
mapping_info = other.mapping_info;
if (other.matrix_info == nullptr)
{
- data = new internal::MatrixFreeFunctions::ShapeInfo<Number>(*other.data);
+ data = new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(*other.data);
scratch_data_array = new AlignedVector<VectorizedArray<Number> >();
}
else
template <int dim, int n_components_, typename Number>
inline
-const internal::MatrixFreeFunctions::ShapeInfo<Number> &
-FEEvaluationBase<dim,n_components_,Number>::get_shape_info() const
+const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
+ FEEvaluationBase<dim,n_components_,Number>::get_shape_info() const
{
Assert(data != nullptr, ExcInternalError());
return *data;
/**
* Return the unit cell information for given hp index.
*/
- const internal::MatrixFreeFunctions::ShapeInfo<Number> &
- get_shape_info (const unsigned int fe_component = 0,
- const unsigned int quad_index = 0,
- const unsigned int hp_active_fe_index = 0,
- const unsigned int hp_active_quad_index = 0) const;
+ const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
+ get_shape_info (const unsigned int fe_component = 0,
+ const unsigned int quad_index = 0,
+ const unsigned int hp_active_fe_index = 0,
+ const unsigned int hp_active_quad_index = 0) const;
/**
* Obtains a scratch data object for internal use. Make sure to release it
/**
* Contains shape value information on the unit cell.
*/
- Table<4,internal::MatrixFreeFunctions::ShapeInfo<Number> > shape_info;
+ Table<4,internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>> shape_info;
/**
* Describes how the cells are gone through. With the cell level (first
template <int dim, typename Number>
inline
-const internal::MatrixFreeFunctions::ShapeInfo<Number> &
-MatrixFree<dim,Number>::get_shape_info (const unsigned int index_fe,
- const unsigned int index_quad,
- const unsigned int active_fe_index,
- const unsigned int active_quad_index) const
+const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
+ MatrixFree<dim,Number>::get_shape_info (const unsigned int index_fe,
+ const unsigned int index_quad,
+ const unsigned int active_fe_index,
+ const unsigned int active_quad_index) const
{
AssertIndexRange (index_fe, shape_info.size(0));
AssertIndexRange (index_quad, shape_info.size(1));
#include <deal.II/base/exceptions.h>
#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/vectorization.h>
#include <deal.II/base/aligned_vector.h>
#include <deal.II/fe/fe.h>
* this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
* points are the index running fastest.
*/
- AlignedVector<VectorizedArray<Number> > shape_values;
+ AlignedVector<Number> shape_values;
/**
* Stores the shape gradients of the 1D finite element evaluated on all
* this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
* points are the index running fastest.
*/
- AlignedVector<VectorizedArray<Number> > shape_gradients;
+ AlignedVector<Number> shape_gradients;
/**
* Stores the shape Hessians of the 1D finite element evaluated on all
* this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
* points are the index running fastest.
*/
- AlignedVector<VectorizedArray<Number> > shape_hessians;
+ AlignedVector<Number> shape_hessians;
/**
* Stores the shape values in a different format, namely the so-called
* even-odd scheme where the symmetries in shape_values are used for
* faster evaluation.
*/
- AlignedVector<VectorizedArray<Number> > shape_values_eo;
+ AlignedVector<Number> shape_values_eo;
/**
* Stores the shape gradients in a different format, namely the so-
* called even-odd scheme where the symmetries in shape_gradients are
* used for faster evaluation.
*/
- AlignedVector<VectorizedArray<Number> > shape_gradients_eo;
+ AlignedVector<Number> shape_gradients_eo;
/**
* Stores the shape second derivatives in a different format, namely the
* so-called even-odd scheme where the symmetries in shape_hessians are
* used for faster evaluation.
*/
- AlignedVector<VectorizedArray<Number> > shape_hessians_eo;
+ AlignedVector<Number> shape_hessians_eo;
/**
* Stores the shape gradients of the shape function space associated to
* the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). For
* faster evaluation only the even-odd format is necessary.
*/
- AlignedVector<VectorizedArray<Number> > shape_gradients_collocation_eo;
+ AlignedVector<Number> shape_gradients_collocation_eo;
/**
* Stores the shape hessians of the shape function space associated to
* the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). For
* faster evaluation only the even-odd format is necessary.
*/
- AlignedVector<VectorizedArray<Number> > shape_hessians_collocation_eo;
+ AlignedVector<Number> shape_hessians_collocation_eo;
/**
* Collects all data of 1D shape values evaluated at the point 0 and 1
* (the vertices) in one data structure. Sorting is first the values,
* then gradients, then second derivatives.
*/
- AlignedVector<VectorizedArray<Number> > shape_data_on_face[2];
+ AlignedVector<Number> shape_data_on_face[2];
/**
* Stores one-dimensional values of shape functions on subface. Since
* there are two subfaces, store two variants.
*/
- AlignedVector<VectorizedArray<Number> > values_within_subface[2];
+ AlignedVector<Number> values_within_subface[2];
/**
* Stores one-dimensional gradients of shape functions on subface. Since
* there are two subfaces, store two variants.
*/
- AlignedVector<VectorizedArray<Number> > gradients_within_subface[2];
+ AlignedVector<Number> gradients_within_subface[2];
/**
* Stores one-dimensional gradients of shape functions on subface. Since
* there are two subfaces, store two variants.
*/
- AlignedVector<VectorizedArray<Number> > hessians_within_subface[2];
+ AlignedVector<Number> hessians_within_subface[2];
/**
* Renumbering from deal.II's numbering of cell degrees of freedom to
bool nodal_at_cell_boundaries;
/**
- * For nodal cells, we might get around by simply loading the indices to
- * the degrees of freedom that act on a particular face, rather than the
- * whole set of indices that is then interpolated down to the
- * element. This array stores this indirect addressing.
+ * For nodal basis functions with nodes located at the boundary of the
+ * unit cell, face integrals that involve only the values of the shape
+ * functions (approximations of first derivatives in DG) do not need to
+ * load all degrees of freedom of the cell but rather only the degrees
+ * of freedom located on the face. While it would also be possible to
+ * compute these indices on the fly, we choose to simplify the code and
+ * store the indirect addressing in this class.
*
* The first table index runs through the faces of a cell, and the
* second runs through the nodal degrees of freedom of the face, using
* @p dofs_per_face entries.
*
+ * The indices stored in this member variable are as follows. Consider
+ * for example a 2D element of degree 3 with the following degrees of
+ * freedom in lexicographic numbering:
+ * @code
+ * 12 13 14 15
+ * 8 9 10 11
+ * 4 5 6 7
+ * 0 1 2 3
+ * @endcode
+ *
+ * The first row stores the indices on the face with index 0, i.e., the
+ * numbers <code>0, 4, 8, 12</code>, the second row holds the indices
+ * <code>3, 7, 11, 15</code> for face 1, the third row holds the indices
+ * <code>0, 1, 2, 3</code> for face 2, and the last (fourth) row holds
+ * the indices <code>12, 13, 14, 15</code>. Similarly, the indices are
+ * stored in 3D. (Note that the y faces in 3D use indices reversed in
+ * terms of the lexicographic numbers due to the orientation of the
+ * coordinate system.)
+ *
* @note This object is only filled in case @p nodal_at_cell_boundaries
* evaluates to @p true.
*/
dealii::Table<2,unsigned int> face_to_cell_index_nodal;
/**
- * For Hermite-type basis functions, the @p face_to_cell_index_nodal for
- * the values on a face of the cell is used together with a respective
- * slot in the derivative. In the lexicographic ordering, this index is
- * in the next "layer" of degrees of freedom. This array stores the
- * indirect addressing of both the values and the gradient.
+ * The @p face_to_cell_index_nodal provides a shortcut for the
+ * evaluation of values on the faces. For Hermite-type basis functions,
+ * one can go one step further and also use shortcuts to get derivatives
+ * more cheaply where only two layers of degrees of freedom contribute
+ * to the derivative on the face. In the lexicographic ordering, the
+ * respective indices is in the next "layer" of degrees of freedom as
+ * compared to the nodal values. This array stores the indirect
+ * addressing of both the values and the gradient.
*
* The first table index runs through the faces of a cell, and the
* second runs through the pairs of the nodal degrees of freedom of the
* face and the derivatives, using <code>2*dofs_per_face</code> entries.
*
+ * The indices stored in this member variable are as follows. Consider
+ * for example a 2D element of degree 3 with the following degrees of
+ * freedom in lexicographic numbering:
+ * @code
+ * 20 21 22 23 24
+ * 15 16 17 18 19
+ * 10 11 12 13 14
+ * 5 6 7 8 9
+ * 0 1 2 3 4
+ * @endcode
+ *
+ * The first row stores the indices for values and gradients on the face
+ * with index 0, i.e., the numbers <code>0, 1, 5, 6, 10, 11, 15, 16, 20,
+ * 21</code>, the second row holds the indices <code>4, 3, 9, 8, 14, 13,
+ * 19, 18, 24, 23</code> for face 1, the third row holds the indices
+ * <code>0, 5, 1, 6, 2, 7, 3, 8, 4, 9</code> for face 2, and the last
+ * (fourth) row holds the indices <code>20, 15, 21, 16, 22, 17, 23, 18,
+ * 24, 19</code>. Similarly, the indices are stored in 3D. (Note that
+ * the y faces in 3D use indices reversed in terms of the lexicographic
+ * numbers due to the orientation of the coordinate system.)
+ *
* @note This object is only filled in case @p element_type evaluates to
* @p tensor_symmetric_hermite.
*/
// ----------------- actual ShapeInfo functions --------------------
+ namespace
+ {
+ template <typename Number>
+ Number get_first_array_element(const Number a)
+ {
+ return a;
+ }
+
+ template <typename Number>
+ Number get_first_array_element(const VectorizedArray<Number> a)
+ {
+ return a[0];
+ }
+ }
+
template <typename Number>
ShapeInfo<Number>::ShapeInfo ()
:
if (fe->has_support_points())
unit_point = fe->get_unit_support_points()[scalar_lexicographic[0]];
Assert(fe->dofs_per_cell == 0 ||
- std::fabs(fe->shape_value(scalar_lexicographic[0],
- unit_point)-1) < 1e-13,
+ std::abs(fe->shape_value(scalar_lexicographic[0],
+ unit_point)-1) < 1e-13,
ExcInternalError("Could not decode 1D shape functions for the "
"element " + fe->get_name()));
}
const unsigned int my_i = scalar_lexicographic[i];
for (unsigned int q=0; q<n_q_points_1d; ++q)
{
- // fill both vectors with
- // VectorizedArray<Number>::n_array_elements copies for the
- // shape information and non-vectorized fields
Point<dim> q_point = unit_point;
q_point[0] = quad.get_points()[q][0];
// check if we are a Hermite type
element_type = tensor_symmetric_hermite;
for (unsigned int i=1; i<n_dofs_1d; ++i)
- if (std::abs(this->shape_data_on_face[0][i][0]) > 1e-12)
+ if (std::abs(get_first_array_element(shape_data_on_face[0][i])) > 1e-12)
element_type = tensor_symmetric;
for (unsigned int i=2; i<n_dofs_1d; ++i)
- if (std::abs(this->shape_data_on_face[0][n_dofs_1d+i][0]) > 1e-12)
+ if (std::abs(get_first_array_element(shape_data_on_face[0][n_dofs_1d+i])) > 1e-12)
element_type = tensor_symmetric;
}
}
nodal_at_cell_boundaries = true;
for (unsigned int i=1; i<n_dofs_1d; ++i)
- if (std::abs(this->shape_data_on_face[0][i][0]) > 1e-13 ||
- std::abs(this->shape_data_on_face[1][i-1][0]) > 1e-13)
+ if (std::abs(get_first_array_element(shape_data_on_face[0][i])) > 1e-13 ||
+ std::abs(get_first_array_element(shape_data_on_face[1][i-1])) > 1e-13)
nodal_at_cell_boundaries = false;
if (nodal_at_cell_boundaries == true)
const unsigned int n_dofs_1d = fe_degree + 1;
for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
for (unsigned int j=0; j<n_q_points_1d; ++j)
- if (std::fabs(shape_values[i*n_q_points_1d+j][0] -
- shape_values[(n_dofs_1d-i)*n_q_points_1d
- -j-1][0]) >
+ if (std::abs(get_first_array_element(shape_values[i*n_q_points_1d+j] -
+ shape_values[(n_dofs_1d-i)*n_q_points_1d-j-1])) >
std::max(zero_tol, zero_tol*
- std::abs(shape_values[i*n_q_points_1d+j][0])))
+ std::abs(get_first_array_element(shape_values[i*n_q_points_1d+j]))))
return false;
// shape values should be zero at x=0.5 for all basis functions except
if (n_q_points_1d%2 == 1 && n_dofs_1d%2 == 1)
{
for (unsigned int i=0; i<n_dofs_1d/2; ++i)
- if (std::fabs(shape_values[i*n_q_points_1d+
- n_q_points_1d/2][0]) > zero_tol)
+ if (std::abs(get_first_array_element(shape_values[i*n_q_points_1d+
+ n_q_points_1d/2])) > zero_tol)
return false;
- if (std::fabs(shape_values[(n_dofs_1d/2)*n_q_points_1d+
- n_q_points_1d/2][0]-1.)> zero_tol)
+ if (std::abs(get_first_array_element(shape_values[(n_dofs_1d/2)*n_q_points_1d+
+ n_q_points_1d/2])-1.)> zero_tol)
return false;
}
const double zero_tol_gradient = zero_tol * std::sqrt(fe_degree+1.)*(fe_degree+1);
for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
for (unsigned int j=0; j<n_q_points_1d; ++j)
- if (std::fabs(shape_gradients[i*n_q_points_1d+j][0] +
- shape_gradients[(n_dofs_1d-i)*n_q_points_1d-
- j-1][0]) > zero_tol_gradient)
+ if (std::abs(get_first_array_element(shape_gradients[i*n_q_points_1d+j] +
+ shape_gradients[(n_dofs_1d-i)*n_q_points_1d-
+ j-1])) > zero_tol_gradient)
return false;
if (n_dofs_1d%2 == 1 && n_q_points_1d%2 == 1)
- if (std::fabs(shape_gradients[(n_dofs_1d/2)*n_q_points_1d+
- (n_q_points_1d/2)][0]) > zero_tol_gradient)
+ if (std::abs(get_first_array_element(shape_gradients[(n_dofs_1d/2)*n_q_points_1d+
+ (n_q_points_1d/2)]))
+ > zero_tol_gradient)
return false;
// symmetry for Hessian. Multiply tolerance by degree^3 of the element
const double zero_tol_hessian = zero_tol * (fe_degree+1)*(fe_degree+1)*(fe_degree+1);
for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
for (unsigned int j=0; j<n_q_points_1d; ++j)
- if (std::fabs(shape_hessians[i*n_q_points_1d+j][0] -
- shape_hessians[(n_dofs_1d-i)*n_q_points_1d-
- j-1][0]) > zero_tol_hessian)
+ if (std::abs(get_first_array_element(shape_hessians[i*n_q_points_1d+j] -
+ shape_hessians[(n_dofs_1d-i)*n_q_points_1d-
+ j-1])) > zero_tol_hessian)
return false;
const unsigned int stride = (n_q_points_1d+1)/2;
for (unsigned int q=0; q<stride; ++q)
{
shape_values_eo[i*stride+q] =
- Number(0.5) * (shape_values[i*n_q_points_1d+q] +
- shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
+ 0.5 * (shape_values[i*n_q_points_1d+q] +
+ shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
shape_values_eo[(fe_degree-i)*stride+q] =
- Number(0.5) * (shape_values[i*n_q_points_1d+q] -
- shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
+ 0.5 * (shape_values[i*n_q_points_1d+q] -
+ shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
shape_gradients_eo[i*stride+q] =
- Number(0.5) * (shape_gradients[i*n_q_points_1d+q] +
- shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
+ 0.5 * (shape_gradients[i*n_q_points_1d+q] +
+ shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
shape_gradients_eo[(fe_degree-i)*stride+q] =
- Number(0.5) * (shape_gradients[i*n_q_points_1d+q] -
- shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
+ 0.5 * (shape_gradients[i*n_q_points_1d+q] -
+ shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
shape_hessians_eo[i*stride+q] =
- Number(0.5) * (shape_hessians[i*n_q_points_1d+q] +
- shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
+ 0.5 * (shape_hessians[i*n_q_points_1d+q] +
+ shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
shape_hessians_eo[(fe_degree-i)*stride+q] =
- Number(0.5) * (shape_hessians[i*n_q_points_1d+q] -
- shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
+ 0.5 * (shape_hessians[i*n_q_points_1d+q] -
+ shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
}
if (fe_degree % 2 == 0)
for (unsigned int q=0; q<stride; ++q)
for (unsigned int j=0; j<n_points_1d; ++j)
if (i!=j)
{
- if (std::fabs(shape_values[i*n_points_1d+j][0])>zero_tol)
+ if (std::abs(get_first_array_element(shape_values[i*n_points_1d+j]))>zero_tol)
return false;
}
else
{
- if (std::fabs(shape_values[i*n_points_1d+
- j][0]-1.)>zero_tol)
+ if (std::abs(get_first_array_element(shape_values[i*n_points_1d+j])-1.)>zero_tol)
return false;
}
return true;
#include <deal.II/matrix_free/matrix_free.templates.h>
#include <deal.II/base/utilities.h>
+#include <deal.II/base/vectorization.h>
#include <deal.II/base/conditional_ostream.h>
#include <iostream>
template struct internal::MatrixFreeFunctions::ShapeInfo<double>;
template struct internal::MatrixFreeFunctions::ShapeInfo<float>;
+template struct internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<double>>;
+template struct internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<float>>;
DEAL_II_NAMESPACE_CLOSE
projects / dealii.git / commitdiff