Updated tests, reordering of components in face evaluation, and storing jacobian for affine face evaluation
template <typename EvalType>
static EvalType
create_evaluator_tensor_product(
- const MatrixFreeFunctions::UnivariateShapeData<Number> &data,
- const unsigned int subface_index,
- const unsigned int direction)
+ const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const unsigned int subface_index,
+ const unsigned int direction)
{
+ const MatrixFreeFunctions::UnivariateShapeData<Number> &data =
+ (std::is_same<EvalType, EvalNormal>::value) ? shape_info.data.front() :
+ shape_info.data.back();
if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
return EvalType(data.shape_values,
data.shape_gradients,
template <bool integrate>
static void
evaluate_or_integrate_in_face(
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
- Number * values_dofs,
- FEEvaluationData<dim, Number, true> & fe_eval,
- Number * scratch_data,
- const unsigned int subface_index,
- const unsigned int face_no)
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ Number * values_dofs,
+ FEEvaluationData<dim, Number, true> & fe_eval,
+ Number * scratch_data,
+ const unsigned int subface_index,
+ const unsigned int face_no)
{
- // TODO. Make sure hanging nodes also are supported.
- // The following part probably needs a rethink.
- EvalNormal eval_normal =
- create_evaluator_tensor_product<EvalNormal>(shape_info.data.front(),
- subface_index,
- 0);
- EvalTangent eval_tangent =
- create_evaluator_tensor_product<EvalTangent>(shape_info.data.back(),
- subface_index,
- 1);
-
- // Used for normal faces which are isotropic
- EvalGeneral eval_general =
- create_evaluator_tensor_product<EvalGeneral>(shape_info.data.back(),
- subface_index,
- 0);
-
- // Note, n_dofs on tangent face
- const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
- const std::size_t n_dofs_normal =
- n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
-
const unsigned int face_direction = face_no / 2;
- if (face_direction == 0)
- {
- evaluate_in_face_apply<-1, 0>(
- eval_general,
- eval_general,
- values_dofs,
- fe_eval,
- scratch_data,
- evaluation_flag,
- n_dofs_normal,
- std::integral_constant<bool, integrate>());
-
- values_dofs += 3 * n_dofs_normal;
-
- evaluate_in_face_apply<0, 1>(
- eval_normal,
- eval_tangent,
- values_dofs,
- fe_eval,
- scratch_data,
- evaluation_flag,
- n_dofs_tangent,
- std::integral_constant<bool, integrate>());
-
- values_dofs += 3 * n_dofs_tangent;
-
- if (dim == 3)
- {
- evaluate_in_face_apply<1, 2>(
- eval_tangent,
- eval_normal,
- values_dofs,
- fe_eval,
- scratch_data,
- evaluation_flag,
- n_dofs_tangent,
- std::integral_constant<bool, integrate>());
- }
- }
- else if (face_direction == 1)
- {
- // NOTE. Take zx-coordinates into account for dim == 3
- if (dim == 3)
- evaluate_in_face_apply<1, 0>(
- eval_tangent,
- eval_normal,
- values_dofs,
- fe_eval,
- scratch_data,
- evaluation_flag,
- n_dofs_tangent,
- std::integral_constant<bool, integrate>());
- else
- evaluate_in_face_apply<0, 0>(
- eval_normal,
- eval_tangent,
- values_dofs,
- fe_eval,
- scratch_data,
- evaluation_flag,
- n_dofs_tangent,
- std::integral_constant<bool, integrate>());
-
- values_dofs += 3 * n_dofs_tangent;
+ evaluate_in_face_apply<0, EvalNormal, EvalTangent>(
+ values_dofs,
+ fe_eval,
+ scratch_data,
+ evaluation_flag,
+ face_direction,
+ subface_index,
+ std::integral_constant<bool, integrate>());
- evaluate_in_face_apply<-1, 1>(
- eval_general,
- eval_general,
- values_dofs,
- fe_eval,
- scratch_data,
- evaluation_flag,
- n_dofs_normal,
- std::integral_constant<bool, integrate>());
-
- values_dofs += 3 * n_dofs_normal;
-
- if (dim == 3)
- {
- // NOTE. Take zx-coordinates into account
- evaluate_in_face_apply<0, 2>(
- eval_normal,
- eval_tangent,
- values_dofs,
- fe_eval,
- scratch_data,
- evaluation_flag,
- n_dofs_tangent,
- std::integral_constant<bool, integrate>());
- }
- }
- else
- {
- evaluate_in_face_apply<0, 0>(
- eval_normal,
- eval_tangent,
- values_dofs,
- fe_eval,
- scratch_data,
- evaluation_flag,
- n_dofs_tangent,
- std::integral_constant<bool, integrate>());
-
- values_dofs += 3 * n_dofs_tangent;
-
- evaluate_in_face_apply<1, 1>(
- eval_tangent,
- eval_normal,
- values_dofs,
- fe_eval,
- scratch_data,
- evaluation_flag,
- n_dofs_tangent,
- std::integral_constant<bool, integrate>());
-
- values_dofs += 3 * n_dofs_tangent;
+ if (dim == 3)
+ evaluate_in_face_apply<1, EvalTangent, EvalNormal>(
+ values_dofs,
+ fe_eval,
+ scratch_data,
+ evaluation_flag,
+ face_direction,
+ subface_index,
+ std::integral_constant<bool, integrate>());
- if (dim == 3)
- {
- evaluate_in_face_apply<-1, 2>(
- eval_general,
- eval_general,
- values_dofs,
- fe_eval,
- scratch_data,
- evaluation_flag,
- n_dofs_normal,
- std::integral_constant<bool, integrate>());
- }
- }
+ evaluate_in_face_apply<2, EvalGeneral, EvalGeneral>(
+ values_dofs,
+ fe_eval,
+ scratch_data,
+ evaluation_flag,
+ face_direction,
+ subface_index,
+ std::integral_constant<bool, integrate>());
}
/*
* Helper function which applies the 1D kernels for on one
* component in a face. normal_dir indicates the direction of the continuous
* component of the RT space. std::integral_constant<bool, false> is the
- * evaluation path, and std::integral_constant<bool, true> bellow is the
- * integration path.
+ * evaluation path, and std::integral_constant<bool, true> below is the
+ * integration path. These two functions can be fused together since all
+ * offsets and pointers are the exact same.
*/
- template <int normal_dir, int component, typename Eval0, typename Eval1>
+ template <int normal_dir, typename Eval0, typename Eval1>
static inline void
evaluate_in_face_apply(
- const Eval0 & eval0,
- const Eval1 & eval1,
Number * values_dofs,
FEEvaluationData<dim, Number, true> & fe_eval,
Number * scratch_data,
const EvaluationFlags::EvaluationFlags evaluation_flag,
- const std::size_t dofs_stride,
+ const unsigned int face_direction,
+ const unsigned int subface_index,
std::integral_constant<bool, false>)
{
- constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
+ // TODO. Need to make sure hanging nodes work, i.e call
+ // create_eval_tensor_product with correct direction.
+ Eval0 eval0 =
+ create_evaluator_tensor_product<Eval0>(fe_eval.get_shape_info(),
+ subface_index,
+ 0);
+ Eval1 eval1 =
+ create_evaluator_tensor_product<Eval1>(fe_eval.get_shape_info(),
+ subface_index,
+ 0);
- Number *values_quad = fe_eval.begin_values() + n_q_points * component;
+ constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
+ const std::size_t n_dofs_tangent =
+ fe_eval.get_shape_info().dofs_per_component_on_face;
+ const std::size_t n_dofs_normal =
+ n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
+ const std::size_t dofs_stride =
+ (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
+ n_dofs_tangent;
+
+ const unsigned int component_table[3][3] = {{1, 2, 0},
+ {2, 0, 1},
+ {0, 1, 2}};
+ const unsigned int component =
+ (dim == 2 && normal_dir == 0 && face_direction == 1) ?
+ 0 :
+ component_table[face_direction][normal_dir];
+
+ // Initial offsets
+ values_dofs +=
+ 3 * ((component == 0) ?
+ 0 :
+ ((component == 1) ?
+ ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
+ ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
+ n_dofs_normal + n_dofs_tangent)));
+ const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
+ Number *values_quad = fe_eval.begin_values() + n_q_points * shift;
Number *gradients_quad =
- fe_eval.begin_gradients() + dim * n_q_points * component;
+ fe_eval.begin_gradients() + dim * n_q_points * shift;
Number *hessians_quad =
- fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * component;
+ fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift;
// Evaluation path
-
if ((evaluation_flag & EvaluationFlags::values) &&
!(evaluation_flag & EvaluationFlags::gradients))
{
}
}
- template <int normal_dir, int component, typename Eval0, typename Eval1>
+ template <int normal_dir, typename Eval0, typename Eval1>
static inline void
evaluate_in_face_apply(
- const Eval0 & eval0,
- const Eval1 & eval1,
Number * values_dofs,
FEEvaluationData<dim, Number, true> & fe_eval,
Number * scratch_data,
const EvaluationFlags::EvaluationFlags evaluation_flag,
- const std::size_t dofs_stride,
+ const unsigned int face_direction,
+ const unsigned int subface_index,
std::integral_constant<bool, true>)
{
- constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
+ // TODO. Need to make sure hanging nodes work, i.e call
+ // create_eval_tensor_product with correct direction.
+ Eval0 eval0 =
+ create_evaluator_tensor_product<Eval0>(fe_eval.get_shape_info(),
+ subface_index,
+ 0);
+ Eval1 eval1 =
+ create_evaluator_tensor_product<Eval1>(fe_eval.get_shape_info(),
+ subface_index,
+ 0);
- Number *values_quad = fe_eval.begin_values() + n_q_points * component;
+ constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
+ const std::size_t n_dofs_tangent =
+ fe_eval.get_shape_info().dofs_per_component_on_face;
+ const std::size_t n_dofs_normal =
+ n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
+ const std::size_t dofs_stride =
+ (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
+ n_dofs_tangent;
+
+ const unsigned int component_table[3][3] = {{1, 2, 0},
+ {2, 0, 1},
+ {0, 1, 2}};
+ const unsigned int component =
+ (dim == 2 && normal_dir == 0 && face_direction == 1) ?
+ 0 :
+ component_table[face_direction][normal_dir];
+
+ // Initial offsets
+ values_dofs +=
+ 3 * ((component == 0) ?
+ 0 :
+ ((component == 1) ?
+ ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
+ ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
+ n_dofs_normal + n_dofs_tangent)));
+ const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
+ Number *values_quad = fe_eval.begin_values() + n_q_points * shift;
Number *gradients_quad =
- fe_eval.begin_gradients() + dim * n_q_points * component;
+ fe_eval.begin_gradients() + dim * n_q_points * shift;
Number *hessians_quad =
- fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * component;
+ fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift;
// Integration path
if ((evaluation_flag & EvaluationFlags::values) &&
(n_q_points_1d < 1) ? 1 :
n_q_points_1d,
Number>::
- template evaluate_or_integrate_in_face<false>(
- evaluation_flag,
- shape_info,
- temp,
- fe_eval,
- scratch_data,
- subface_index,
- fe_eval.get_face_no());
+ template evaluate_or_integrate_in_face<false>(evaluation_flag,
+ temp,
+ fe_eval,
+ scratch_data,
+ subface_index,
+ fe_eval.get_face_no());
}
else if (fe_degree > -1 &&
subface_index >= GeometryInfo<dim>::max_children_per_cell &&
n_q_points_1d,
Number>::
template evaluate_or_integrate_in_face<true>(integration_flag,
- shape_info,
temp,
fe_eval,
scratch_data,
using BaseClass =
FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>;
+ /**
+ * @copydoc FEEvaluationBase<dim,dim,Number,is_face>::get_value()
+ */
+ value_type
+ get_value(const unsigned int q_point) const;
+
/**
* @copydoc FEEvaluationBase<dim,dim,Number,is_face>::get_gradient()
*/
gradient_type
get_hessian_diagonal(const unsigned int q_point) const;
+ /**
+ * @copydoc FEEvaluationBase<dim,dim,Number,is_face>::submit_value()
+ */
+ void
+ submit_value(const Tensor<1, dim, VectorizedArrayType> val_in,
+ const unsigned int q_point);
+
/**
* @copydoc FEEvaluationBase<dim,dim,Number,is_face>::submit_gradient()
*/
}
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
+FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::get_value(
+ const unsigned int q_point) const
+{
+ // Check if Piola transform is required
+ if (this->data->element_type ==
+ internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
+ {
+# ifdef DEBUG
+ Assert(this->values_quad_initialized == true,
+ internal::ExcAccessToUninitializedField());
+# endif
+
+ AssertIndexRange(q_point, this->n_quadrature_points);
+ Assert(this->J_value != nullptr,
+ internal::ExcMatrixFreeAccessToUninitializedMappingField(
+ "update_values"));
+ const std::size_t nqp = this->n_quadrature_points;
+ Tensor<1, dim, VectorizedArrayType> value_out;
+
+ // Cartesian cell
+ if (!is_face &&
+ this->cell_type == internal::MatrixFreeFunctions::cartesian)
+ {
+ const Tensor<2, dim, dealii::VectorizedArray<Number>> jac =
+ this->jacobian[1];
+ const VectorizedArrayType inv_det = determinant(this->jacobian[0]);
+
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ value_out[comp] = this->values_quad[comp * nqp + q_point] *
+ jac[comp][comp] *
+ inv_det; // / this->jacobian[0][comp][comp];
+ }
+
+ // Affine or general cell
+ else
+ {
+ const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
+ (this->cell_type > internal::MatrixFreeFunctions::affine) ?
+ this->jacobian[q_point] :
+ this->jacobian[0];
+ const Tensor<2, dim, VectorizedArrayType> &jac =
+ (this->cell_type > internal::MatrixFreeFunctions::affine) ?
+ transpose(invert(inv_t_jac)) :
+ this->jacobian[1];
+
+ // Derivatives are reordered for faces. Need to take this into account
+ const VectorizedArrayType inv_det =
+ (is_face && dim == 2 && this->get_face_no() < 2) ?
+ -determinant(inv_t_jac) :
+ determinant(inv_t_jac);
+ // J * u * det(J^-1)
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ {
+ value_out[comp] =
+ this->values_quad[q_point] * jac[comp][0] * inv_det;
+ for (unsigned int e = 1; e < dim; ++e)
+ value_out[comp] +=
+ this->values_quad[e * nqp + q_point] * jac[comp][e] * inv_det;
+ }
+ }
+ return value_out;
+ }
+ else
+ {
+ return BaseClass::get_value(q_point);
+ }
+}
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
get_gradient(const unsigned int q_point) const
{
- return BaseClass::get_gradient(q_point);
+ // Check if Piola transform is required
+ if (this->data->element_type ==
+ internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
+ {
+# ifdef DEBUG
+ Assert(this->gradients_quad_initialized == true,
+ internal::ExcAccessToUninitializedField());
+# endif
+
+ AssertIndexRange(q_point, this->n_quadrature_points);
+ Assert(this->jacobian != nullptr,
+ internal::ExcMatrixFreeAccessToUninitializedMappingField(
+ "update_gradients"));
+ const std::size_t nqp = this->n_quadrature_points;
+ Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_out;
+
+ // Cartesian cell
+ if (!is_face &&
+ this->cell_type == internal::MatrixFreeFunctions::cartesian)
+ {
+ const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
+ this->jacobian[0];
+ const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
+ const VectorizedArrayType inv_det = determinant(inv_t_jac);
+
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ grad_out[comp][d] =
+ this->gradients_quad[(comp * dim + d) * nqp + q_point] *
+ inv_t_jac[d][d] * jac[comp][comp] * inv_det;
+ }
+ // Affine cell
+ else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+ {
+ const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
+ this->jacobian[0];
+ const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
+
+ // Derivatives are reordered for faces. Need to take this into account
+ const VectorizedArrayType inv_det =
+ (is_face && dim == 2 && this->get_face_no() < 2) ?
+ -determinant(inv_t_jac) :
+ determinant(inv_t_jac);
+
+ VectorizedArrayType tmp;
+ // J * grad_quad * J^-1 * det(J^-1)
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ tmp = 0;
+ for (unsigned int f = 0; f < dim; ++f)
+ for (unsigned int e = 0; e < dim; ++e)
+ tmp += jac[comp][f] * inv_t_jac[d][e] * inv_det *
+ this->gradients_quad[(f * dim + e) * nqp + q_point];
+
+ grad_out[comp][d] = tmp;
+ }
+ }
+ // General cell TODO
+ else
+ {
+ // Here we need the jacobian gradient and not the inverse which is
+ // stored in this->jacobian_gradients
+ AssertThrow(false, ExcNotImplemented());
+ }
+ return grad_out;
+ }
+ else
+ {
+ return BaseClass::get_gradient(q_point);
+ }
}
VectorizedArrayType divergence;
const std::size_t nqp = this->n_quadrature_points;
- // Cartesian cell
- if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
+ if (this->data->element_type ==
+ internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
{
- divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
- for (unsigned int d = 1; d < dim; ++d)
- divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
- this->jacobian[0][d][d];
+ // Affine cell
+ if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+ {
+ // Derivatives are reordered for faces. Need to take this into account
+ const VectorizedArrayType inv_det =
+ (is_face && dim == 2 && this->get_face_no() < 2) ?
+ -determinant(this->jacobian[0]) :
+ determinant(this->jacobian[0]);
+
+ divergence = this->gradients_quad[q_point] * inv_det;
+ for (unsigned int d = 1; d < dim; ++d)
+ divergence +=
+ this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
+ }
+ // General cell TODO
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ }
}
- // cell with general/constant Jacobian
else
{
- const Tensor<2, dim, VectorizedArrayType> &jac =
- this->cell_type == internal::MatrixFreeFunctions::general ?
- this->jacobian[q_point] :
- this->jacobian[0];
- divergence = jac[0][0] * this->gradients_quad[q_point];
- for (unsigned int e = 1; e < dim; ++e)
- divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
- for (unsigned int d = 1; d < dim; ++d)
- for (unsigned int e = 0; e < dim; ++e)
- divergence +=
- jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
+ // Cartesian cell
+ if (!is_face &&
+ this->cell_type == internal::MatrixFreeFunctions::cartesian)
+ {
+ divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
+ for (unsigned int d = 1; d < dim; ++d)
+ divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
+ this->jacobian[0][d][d];
+ }
+ // cell with general/constant Jacobian
+ else
+ {
+ const Tensor<2, dim, VectorizedArrayType> &jac =
+ this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->jacobian[q_point] :
+ this->jacobian[0];
+ divergence = jac[0][0] * this->gradients_quad[q_point];
+ for (unsigned int e = 1; e < dim; ++e)
+ divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
+ for (unsigned int d = 1; d < dim; ++d)
+ for (unsigned int e = 0; e < dim; ++e)
+ divergence +=
+ jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
+ }
}
return divergence;
}
}
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline DEAL_II_ALWAYS_INLINE void
+FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
+ submit_value(const Tensor<1, dim, VectorizedArrayType> val_in,
+ const unsigned int q_point)
+{
+ // Check if Piola transform is required
+ if (this->data->element_type ==
+ internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
+ {
+ AssertIndexRange(q_point, this->n_quadrature_points);
+
+ // This is not needed, but might be good to check anyway?
+ Assert(this->J_value != nullptr,
+ internal::ExcMatrixFreeAccessToUninitializedMappingField(
+ "update_value"));
+# ifdef DEBUG
+ Assert(this->is_reinitialized, ExcNotInitialized());
+ this->values_quad_submitted = true;
+# endif
+
+ const std::size_t nqp = this->n_quadrature_points;
+ if (!is_face &&
+ this->cell_type == internal::MatrixFreeFunctions::cartesian)
+ {
+ const Tensor<2, dim, dealii::VectorizedArray<Number>> jac =
+ this->jacobian[1];
+ const VectorizedArrayType weight = this->quadrature_weights[q_point];
+
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ this->values_quad[comp * nqp + q_point] =
+ val_in[comp] * weight * jac[comp][comp];
+ }
+ // Affine or general cell
+ else
+ {
+ const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
+ (this->cell_type > internal::MatrixFreeFunctions::affine) ?
+ this->jacobian[q_point] :
+ this->jacobian[0];
+ const Tensor<2, dim, VectorizedArrayType> &jac =
+ (this->cell_type > internal::MatrixFreeFunctions::affine) ?
+ invert(inv_t_jac) :
+ this->jacobian[1];
+
+ // Derivatives are reordered for faces. Need to take this into account
+ // and 1/inv_det != J_value for faces
+ const VectorizedArrayType fac =
+ (!is_face) ?
+ this->quadrature_weights[q_point] :
+ (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
+ this->J_value[q_point] :
+ this->J_value[0] * this->quadrature_weights[q_point]) *
+ ((dim == 2 && this->get_face_no() < 2) ?
+ -determinant(inv_t_jac) :
+ determinant(inv_t_jac)));
+
+ // J^T * u * w
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ {
+ this->values_quad[comp * nqp + q_point] =
+ val_in[0] * jac[0][comp] * fac;
+ for (unsigned int e = 1; e < dim; ++e)
+ this->values_quad[comp * nqp + q_point] +=
+ val_in[e] * jac[e][comp] * fac;
+ }
+ }
+ }
+ else
+ {
+ BaseClass::submit_value(val_in, q_point);
+ }
+}
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE void
submit_gradient(const Tensor<2, dim, VectorizedArrayType> grad_in,
const unsigned int q_point)
{
- BaseClass::submit_gradient(grad_in, q_point);
+ // Check if Piola transform is required
+ if (this->data->element_type ==
+ internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
+ {
+# ifdef DEBUG
+ Assert(this->is_reinitialized, ExcNotInitialized());
+# endif
+ AssertIndexRange(q_point, this->n_quadrature_points);
+ Assert(this->J_value != nullptr,
+ internal::ExcMatrixFreeAccessToUninitializedMappingField(
+ "update_gradients"));
+ Assert(this->jacobian != nullptr,
+ internal::ExcMatrixFreeAccessToUninitializedMappingField(
+ "update_gradients"));
+# ifdef DEBUG
+ this->gradients_quad_submitted = true;
+# endif
+
+ const std::size_t nqp = this->n_quadrature_points;
+ // Cartesian cell
+ if (!is_face &&
+ this->cell_type == internal::MatrixFreeFunctions::cartesian)
+ {
+ const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
+ this->jacobian[0];
+ const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
+ const VectorizedArrayType weight = this->quadrature_weights[q_point];
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ this->gradients_quad[(comp * dim + d) * nqp + q_point] =
+ grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight;
+ }
+ // Affine cell
+ else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+ {
+ const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
+ this->jacobian[0];
+ const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
+
+ // Derivatives are reordered for faces. Need to take this into account
+ // and 1/inv_det != J_value for faces
+ const VectorizedArrayType fac =
+ (!is_face) ? this->quadrature_weights[q_point] :
+ this->J_value[0] * this->quadrature_weights[q_point] *
+ ((dim == 2 && this->get_face_no() < 2) ?
+ -determinant(inv_t_jac) :
+ determinant(inv_t_jac));
+
+ // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ VectorizedArrayType tmp = 0;
+ for (unsigned int f = 0; f < dim; ++f)
+ for (unsigned int e = 0; e < dim; ++e)
+ tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e] * fac;
+
+ this->gradients_quad[(comp * dim + d) * nqp + q_point] = tmp;
+ }
+ }
+ // General cell TODO
+ else
+ {
+ AssertThrow(false, ExcNotImplemented());
+ }
+ }
+ else
+ {
+ BaseClass::submit_gradient(grad_in, q_point);
+ }
}
const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
const unsigned int q_point)
{
- BaseClass::submit_gradient(grad_in, q_point);
+ // Check if Piola transform is required
+ if (this->data->element_type ==
+ internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
+ {
+# ifdef DEBUG
+ Assert(this->is_reinitialized, ExcNotInitialized());
+# endif
+ AssertIndexRange(q_point, this->n_quadrature_points);
+ Assert(this->J_value != nullptr,
+ internal::ExcMatrixFreeAccessToUninitializedMappingField(
+ "update_gradients"));
+ Assert(this->jacobian != nullptr,
+ internal::ExcMatrixFreeAccessToUninitializedMappingField(
+ "update_gradients"));
+# ifdef DEBUG
+ this->gradients_quad_submitted = true;
+# endif
+
+ const std::size_t nqp = this->n_quadrature_points;
+ // Cartesian cell
+ if (!is_face &&
+ this->cell_type == internal::MatrixFreeFunctions::cartesian)
+ {
+ const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
+ this->jacobian[0];
+ const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
+ const VectorizedArrayType weight = this->quadrature_weights[q_point];
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ this->gradients_quad[(comp * dim + d) * nqp + q_point] =
+ grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight;
+ }
+ // Affine cell
+ else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+ {
+ const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
+ this->jacobian[0];
+ const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
+
+ // Derivatives are reordered for faces. Need to take this into account
+ // and 1/inv_det != J_value for faces
+ const VectorizedArrayType fac =
+ (!is_face) ? this->quadrature_weights[q_point] :
+ this->J_value[0] * this->quadrature_weights[q_point] *
+ ((dim == 2 && this->get_face_no() < 2) ?
+ -determinant(inv_t_jac) :
+ determinant(inv_t_jac));
+
+ // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ VectorizedArrayType tmp = 0;
+ for (unsigned int f = 0; f < dim; ++f)
+ for (unsigned int e = 0; e < dim; ++e)
+ tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e] * fac;
+
+ this->gradients_quad[(comp * dim + d) * nqp + q_point] = tmp;
+ }
+ }
+ // General cell TODO
+ else
+ {
+ AssertThrow(false, ExcNotImplemented());
+ }
+ }
+ else
+ {
+ BaseClass::submit_gradient(grad_in, q_point);
+ }
}
# endif
const std::size_t nqp = this->n_quadrature_points;
- if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
+ if (this->data->element_type ==
+ internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
{
- const VectorizedArrayType fac =
- this->J_value[0] * this->quadrature_weights[q_point] * div_in;
- for (unsigned int d = 0; d < dim; ++d)
+ // Affine cell
+ if (this->cell_type <= internal::MatrixFreeFunctions::affine)
{
- this->gradients_quad[(d * dim + d) * nqp + q_point] =
- (fac * this->jacobian[0][d][d]);
- for (unsigned int e = d + 1; e < dim; ++e)
+ // Derivatives are reordered for faces. Need to take this into account
+ // and 1/inv_det != J_value for faces
+ const VectorizedArrayType fac =
+ ((!is_face) ?
+ 1 :
+ this->J_value[0] * ((dim == 2 && this->get_face_no() < 2) ?
+ -determinant(this->jacobian[0]) :
+ determinant(this->jacobian[0]))) *
+ this->quadrature_weights[q_point] * div_in;
+
+ for (unsigned int d = 0; d < dim; ++d)
{
- this->gradients_quad[(d * dim + e) * nqp + q_point] =
- VectorizedArrayType();
- this->gradients_quad[(e * dim + d) * nqp + q_point] =
- VectorizedArrayType();
+ this->gradients_quad[(dim * d + d) * nqp + q_point] = fac;
+ for (unsigned int e = d + 1; e < dim; ++e)
+ {
+ this->gradients_quad[(dim * d + e) * nqp + q_point] =
+ VectorizedArrayType();
+ this->gradients_quad[(dim * e + d) * nqp + q_point] =
+ VectorizedArrayType();
+ }
}
}
+ // General cell TODO
+ else
+ {
+ AssertThrow(false, ExcNotImplemented());
+ }
}
else
{
- const Tensor<2, dim, VectorizedArrayType> &jac =
- this->cell_type == internal::MatrixFreeFunctions::general ?
- this->jacobian[q_point] :
- this->jacobian[0];
- const VectorizedArrayType fac =
- (this->cell_type == internal::MatrixFreeFunctions::general ?
- this->J_value[q_point] :
- this->J_value[0] * this->quadrature_weights[q_point]) *
- div_in;
- for (unsigned int d = 0; d < dim; ++d)
+ if (!is_face &&
+ this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
- for (unsigned int e = 0; e < dim; ++e)
- this->gradients_quad[(d * dim + e) * nqp + q_point] =
- jac[d][e] * fac;
+ const VectorizedArrayType fac =
+ this->J_value[0] * this->quadrature_weights[q_point] * div_in;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ this->gradients_quad[(d * dim + d) * nqp + q_point] =
+ (fac * this->jacobian[0][d][d]);
+ for (unsigned int e = d + 1; e < dim; ++e)
+ {
+ this->gradients_quad[(d * dim + e) * nqp + q_point] =
+ VectorizedArrayType();
+ this->gradients_quad[(e * dim + d) * nqp + q_point] =
+ VectorizedArrayType();
+ }
+ }
+ }
+ else
+ {
+ const Tensor<2, dim, VectorizedArrayType> &jac =
+ this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->jacobian[q_point] :
+ this->jacobian[0];
+ const VectorizedArrayType fac =
+ (this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->J_value[q_point] :
+ this->J_value[0] * this->quadrature_weights[q_point]) *
+ div_in;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ for (unsigned int e = 0; e < dim; ++e)
+ this->gradients_quad[(d * dim + e) * nqp + q_point] =
+ jac[d][e] * fac;
+ }
}
}
}
const SymmetricTensor<2, dim, VectorizedArrayType> sym_grad,
const unsigned int q_point)
{
+ // TODO
+ AssertThrow(
+ this->data->element_type !=
+ internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas,
+ ExcNotImplemented());
+
// could have used base class operator, but that involves some overhead
// which is inefficient. it is nice to have the symmetric tensor because
// that saves some operations
const Point<dim, Number> *quadrature_points;
/**
- * A pointer to the Jacobian information of the present cell. Only set to a
- * useful value if on a non-Cartesian cell.
+ * A pointer to the inverse transpose Jacobian information of the present
+ * cell. Only set to a useful value if on a non-Cartesian cell. If the cell is
+ * Cartesian/affine then the Jacobian is stored at index 1. For faces on
+ * hypercube elements, the derivatives are reorder s.t the derivative
+ * orthogonal to the face is stored last, i.e for dim = 3 and face_no = 0 or
+ * 1, the derivatives are ordered as [dy, dz, dx], face_no = 2 or 3: [dz, dx,
+ * dy], and face_no = 5 or 6: [dx, dy, dz]. If the Jacobian also is stored,
+ * the components are instead reordered in the same way.
*/
const Tensor<2, dim, Number> *jacobian;
vv,
my_data.jacobians[0][offset + q][d][e]);
}
+ if (face_type[face] <= affine)
+ for (unsigned int e = 0; e < dim; ++e)
+ {
+ const unsigned int ee =
+ ExtractFaceHelper::reorder_face_derivative_indices<
+ dim>(interior_face_no, e);
+ for (unsigned int d = 0; d < dim; ++d)
+ store_vectorized_array(
+ jac[d][ee],
+ vv,
+ my_data.jacobians[0][offset + q + 1][d][e]);
+ }
if (update_flags_faces & update_jacobian_grads)
{
vv,
my_data.jacobians[1][offset + q][d][e]);
}
+ if (face_type[face] <= affine)
+ for (unsigned int e = 0; e < dim; ++e)
+ {
+ const unsigned int ee = ExtractFaceHelper::
+ reorder_face_derivative_indices<dim>(
+ exterior_face_no, e);
+ for (unsigned int d = 0; d < dim; ++d)
+ store_vectorized_array(
+ jac[d][ee],
+ vv,
+ my_data.jacobians[1][offset + q + 1][d][e]);
+ }
if (update_flags_faces & update_jacobian_grads)
{
max_size =
std::max(max_size,
my_data.data_index_offsets[face] +
- (face_type[face] <= affine ? 1 : n_q_points));
+ (face_type[face] <= affine ? 2 : n_q_points));
}
const UpdateFlags update_flags_common =
// ---------------------------------------------------------------------
// This function tests the correctness of the matrix-free implementation
-// of the FE_RaviartThomasNodal element by evaluating a simple fe operator
-// and comparing the result with FEVaules which is considered the
-// reference. The mesh is a hypercube mesh with no hanging nodes and no other
-// constraints
+// of the FE_RaviartThomasNodal element by evaluating values + gradients
+// as well as the divergence and comparing the result with FEVaules which
+// is considered the reference. The mesh is a hypercube mesh with no
+// hanging nodes and no other constraints.
#include "../tests.h"
AffineConstraints<double> constraints;
constraints.close();
- do_test<dim, fe_degree, double>(dof, constraints);
+
+ deallog << "Using " << dof.get_fe().get_name() << std::endl;
+ deallog << "Number of cells: " << dof.get_triangulation().n_active_cells()
+ << std::endl;
+ deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl
+ << std::endl;
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::values);
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::gradients);
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::divergence);
}
-DEAL:2d::Testing FE_RaviartThomasNodal<2>(1)
+DEAL:2d::Using FE_RaviartThomasNodal<2>(1)
DEAL:2d::Number of cells: 16
DEAL:2d::Number of degrees of freedom: 144
DEAL:2d::
-DEAL:2d::Norm of difference: 6.20418e-16
+DEAL:2d::Testing Values
+DEAL:2d::Norm of difference: 6.05894e-16
DEAL:2d::
-DEAL:2d::Testing FE_RaviartThomasNodal<2>(2)
+DEAL:2d::Testing Gradients
+DEAL:2d::Norm of difference: 6.33970e-16
+DEAL:2d::
+DEAL:2d::Testing Divergence
+DEAL:2d::Norm of difference: 5.24419e-16
+DEAL:2d::
+DEAL:2d::Using FE_RaviartThomasNodal<2>(2)
DEAL:2d::Number of cells: 16
DEAL:2d::Number of degrees of freedom: 312
DEAL:2d::
-DEAL:2d::Norm of difference: 7.03559e-16
+DEAL:2d::Testing Values
+DEAL:2d::Norm of difference: 6.67799e-16
+DEAL:2d::
+DEAL:2d::Testing Gradients
+DEAL:2d::Norm of difference: 1.06191e-15
+DEAL:2d::
+DEAL:2d::Testing Divergence
+DEAL:2d::Norm of difference: 5.38383e-16
DEAL:2d::
-DEAL:3d::Testing FE_RaviartThomasNodal<3>(1)
+DEAL:3d::Using FE_RaviartThomasNodal<3>(1)
DEAL:3d::Number of cells: 64
DEAL:3d::Number of degrees of freedom: 1728
DEAL:3d::
-DEAL:3d::Norm of difference: 1.04798e-15
+DEAL:3d::Testing Values
+DEAL:3d::Norm of difference: 6.43004e-16
DEAL:3d::
-DEAL:3d::Testing FE_RaviartThomasNodal<3>(2)
+DEAL:3d::Testing Gradients
+DEAL:3d::Norm of difference: 8.38689e-16
+DEAL:3d::
+DEAL:3d::Testing Divergence
+DEAL:3d::Norm of difference: 8.84868e-16
+DEAL:3d::
+DEAL:3d::Using FE_RaviartThomasNodal<3>(2)
DEAL:3d::Number of cells: 64
DEAL:3d::Number of degrees of freedom: 5616
DEAL:3d::
-DEAL:3d::Norm of difference: 1.49180e-15
+DEAL:3d::Testing Values
+DEAL:3d::Norm of difference: 1.10927e-15
+DEAL:3d::
+DEAL:3d::Testing Gradients
+DEAL:3d::Norm of difference: 1.50806e-15
+DEAL:3d::
+DEAL:3d::Testing Divergence
+DEAL:3d::Norm of difference: 1.72044e-15
DEAL:3d::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// This test it the same as matrix_vector_rt_01.cc however with
+// non-Cartesian (but still affine) cells.
+
+#include <deal.II/grid/grid_tools.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_rt_common.h"
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ const unsigned int n_subdivisions = 2;
+ Point<dim> corners[dim];
+ corners[0] = (dim == 2) ? Point<dim>(1, 0) : Point<dim>(1, 0, 0);
+ corners[1] = (dim == 2) ? Point<dim>(0.5, 0.5) : Point<dim>(0.5, 1, 0);
+ if (dim == 3)
+ corners[2] = Point<dim>(0.5, 0, 1);
+ GridGenerator::subdivided_parallelepiped(tria, n_subdivisions, corners);
+
+ FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+
+ AffineConstraints<double> constraints;
+ constraints.close();
+
+ deallog << "Using " << dof.get_fe().get_name() << std::endl;
+ deallog << "Number of cells: " << dof.get_triangulation().n_active_cells()
+ << std::endl;
+ deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl
+ << std::endl;
+
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::values_gradients);
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::divergence);
+}
--- /dev/null
+
+DEAL:2d::Using FE_RaviartThomasNodal<2>(1)
+DEAL:2d::Number of cells: 4
+DEAL:2d::Number of degrees of freedom: 40
+DEAL:2d::
+DEAL:2d::Testing Values and Gradients
+DEAL:2d::Norm of difference: 3.77395e-16
+DEAL:2d::
+DEAL:2d::Testing Divergence
+DEAL:2d::Norm of difference: 7.00345e-16
+DEAL:2d::
+DEAL:2d::Using FE_RaviartThomasNodal<2>(2)
+DEAL:2d::Number of cells: 4
+DEAL:2d::Number of degrees of freedom: 84
+DEAL:2d::
+DEAL:2d::Testing Values and Gradients
+DEAL:2d::Norm of difference: 7.63491e-16
+DEAL:2d::
+DEAL:2d::Testing Divergence
+DEAL:2d::Norm of difference: 7.46533e-16
+DEAL:2d::
+DEAL:3d::Using FE_RaviartThomasNodal<3>(1)
+DEAL:3d::Number of cells: 8
+DEAL:3d::Number of degrees of freedom: 240
+DEAL:3d::
+DEAL:3d::Testing Values and Gradients
+DEAL:3d::Norm of difference: 1.22886e-15
+DEAL:3d::
+DEAL:3d::Testing Divergence
+DEAL:3d::Norm of difference: 1.16405e-15
+DEAL:3d::
+DEAL:3d::Using FE_RaviartThomasNodal<3>(2)
+DEAL:3d::Number of cells: 8
+DEAL:3d::Number of degrees of freedom: 756
+DEAL:3d::
+DEAL:3d::Testing Values and Gradients
+DEAL:3d::Norm of difference: 2.06049e-15
+DEAL:3d::
+DEAL:3d::Testing Divergence
+DEAL:3d::Norm of difference: 1.23490e-15
+DEAL:3d::
void
test();
+enum TestType : unsigned char
+{
+ values = 0,
+ values_gradients = 1,
+ gradients = 2,
+ divergence = 3
+};
+std::string
+enum_to_string(TestType const enum_type)
+{
+ std::string string_type;
+ switch (enum_type)
+ {
+ case TestType::values:
+ string_type = "Values ";
+ break;
+ case TestType::gradients:
+ string_type = "Gradients ";
+ break;
+ case TestType::values_gradients:
+ string_type = "Values and Gradients ";
+ break;
+ case TestType::divergence:
+ string_type = "Divergence ";
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ break;
+ }
+ return string_type;
+}
template <int dim,
int fe_degree,
class MatrixFreeTest
{
public:
- MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
- : data(data_in){};
+ MatrixFreeTest(const MatrixFree<dim, Number> &data_in,
+ const TestType test_type)
+ : data(data_in)
+ , test_type(test_type)
+ {
+ evaluation_flag =
+ (test_type == TestType::values) ?
+ EvaluationFlags::values :
+ ((test_type == TestType::gradients) ?
+ EvaluationFlags::gradients :
+ ((test_type == TestType::values_gradients) ?
+ EvaluationFlags::values | EvaluationFlags::gradients :
+ EvaluationFlags::gradients));
+ };
virtual ~MatrixFreeTest(){};
{
FEEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_eval(data);
- // OBS! This will need to be modified once the Piola transform is
- // implemented
- unsigned int n_cells =
- data.get_dof_handler().get_triangulation().n_active_cells();
- Number piola =
- (dim == 2) ? n_cells : Utilities::pow((int)std::cbrt(n_cells), 4);
-
for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
{
fe_eval.reinit(cell);
- fe_eval.gather_evaluate(src,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
+ fe_eval.gather_evaluate(src, evaluation_flag);
for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
{
- fe_eval.submit_value(Number(10 * piola) * fe_eval.get_value(q), q);
- fe_eval.submit_gradient(Number(piola) * fe_eval.get_gradient(q), q);
+ if (test_type < TestType::gradients)
+ fe_eval.submit_value(Number(10) * fe_eval.get_value(q), q);
+ if (test_type == TestType::gradients ||
+ test_type == TestType::values_gradients)
+ fe_eval.submit_gradient(fe_eval.get_gradient(q), q);
+ else if (test_type == TestType::divergence)
+ fe_eval.submit_divergence(fe_eval.get_divergence(q), q);
}
- fe_eval.integrate_scatter(EvaluationFlags::values |
- EvaluationFlags::gradients,
- dst);
+ fe_eval.integrate_scatter(evaluation_flag, dst);
}
};
};
protected:
- const MatrixFree<dim, Number> &data;
+ const MatrixFree<dim, Number> & data;
+ EvaluationFlags::EvaluationFlags evaluation_flag;
+ const TestType test_type;
};
template <int dim, int fe_degree, typename Number>
void
do_test(const DoFHandler<dim> & dof,
- const AffineConstraints<double> &constraints)
+ const AffineConstraints<double> &constraints,
+ const TestType test_type)
{
- deallog << "Testing " << dof.get_fe().get_name() << std::endl;
- deallog << "Number of cells: " << dof.get_triangulation().n_active_cells()
- << std::endl;
- deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl
- << std::endl;
-
+ deallog << "Testing " << enum_to_string(test_type) << std::endl;
// constraints.distribute(solution);
MatrixFree<dim, Number> mf_data;
}
// MatrixFree solution
- MatrixFreeTest<dim, fe_degree, fe_degree + 2, Number> mf(mf_data);
+ MatrixFreeTest<dim, fe_degree, fe_degree + 2, Number> mf(mf_data, test_type);
mf.test_functions(solution, initial_condition);
{
const Tensor<1, dim> phi_i = fe_val[velocities].value(i, q) * 10.;
const Tensor<2, dim> grad_phi_i = fe_val[velocities].gradient(i, q);
+ const Number div_phi_i = fe_val[velocities].divergence(i, q);
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
const Tensor<1, dim> phi_j = fe_val[velocities].value(j, q);
const Tensor<2, dim> grad_phi_j =
fe_val[velocities].gradient(j, q);
-
- local_matrix(i, j) +=
- (phi_j * phi_i + scalar_product(grad_phi_i, grad_phi_j)) *
- fe_val.JxW(q);
+ const Number div_phi_j = fe_val[velocities].divergence(j, q);
+
+ if (test_type < TestType::gradients)
+ local_matrix(i, j) += phi_j * phi_i * fe_val.JxW(q);
+ if (test_type == TestType::gradients ||
+ test_type == TestType::values_gradients)
+ local_matrix(i, j) +=
+ scalar_product(grad_phi_i, grad_phi_j) * fe_val.JxW(q);
+ else if (test_type == TestType::divergence)
+ local_matrix(i, j) += div_phi_i * div_phi_j * fe_val.JxW(q);
}
}
cell->get_dof_indices(local_dof_indices);
// ---------------------------------------------------------------------
// This function tests the correctness of the matrix-free implementation
-// of the FE_RaviartThomasNodal element by evaluating a face operator
-// and comparing the result with FEVaules which is considered the
-// reference. The mesh is a hypercube mesh with no hanging nodes and no other
-// constraints
+// of the FE_RaviartThomasNodal element by evaluating values + gradients
+// as well as the divergence on faces and comparing the result with
+// FEFaceVaules which is considered the reference. The mesh is a hypercube
+// mesh with no hanging nodes and no other constraints.
#include "../tests.h"
AffineConstraints<double> constraints;
constraints.close();
- do_test<dim, fe_degree, double>(dof, constraints);
+ deallog << "Using " << dof.get_fe().get_name() << std::endl;
+ deallog << "Number of cells: " << dof.get_triangulation().n_active_cells()
+ << std::endl;
+ deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl
+ << std::endl;
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::values_gradients);
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::divergence);
}
-DEAL:2d::Testing FE_RaviartThomasNodal<2>(1)
+DEAL:2d::Using FE_RaviartThomasNodal<2>(1)
DEAL:2d::Number of cells: 4
DEAL:2d::Number of degrees of freedom: 40
DEAL:2d::
-DEAL:2d::Norm of difference: 3.44596e-16
+DEAL:2d::Testing Values and Gradients
+DEAL:2d::Norm of difference: 2.36910e-16
DEAL:2d::
-DEAL:2d::Testing FE_RaviartThomasNodal<2>(2)
+DEAL:2d::Testing Divergence
+DEAL:2d::Norm of difference: 4.52039e-16
+DEAL:2d::
+DEAL:2d::Using FE_RaviartThomasNodal<2>(2)
DEAL:2d::Number of cells: 4
DEAL:2d::Number of degrees of freedom: 84
DEAL:2d::
-DEAL:2d::Norm of difference: 3.47432e-15
+DEAL:2d::Testing Values and Gradients
+DEAL:2d::Norm of difference: 7.31531e-15
+DEAL:2d::
+DEAL:2d::Testing Divergence
+DEAL:2d::Norm of difference: 2.97155e-15
DEAL:2d::
-DEAL:3d::Testing FE_RaviartThomasNodal<3>(1)
+DEAL:3d::Using FE_RaviartThomasNodal<3>(1)
DEAL:3d::Number of cells: 8
DEAL:3d::Number of degrees of freedom: 240
DEAL:3d::
-DEAL:3d::Norm of difference: 2.28531e-15
+DEAL:3d::Testing Values and Gradients
+DEAL:3d::Norm of difference: 3.73456e-15
+DEAL:3d::
+DEAL:3d::Testing Divergence
+DEAL:3d::Norm of difference: 4.51598e-15
DEAL:3d::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// This test it the same as matrix_vector_rt_face_01.cc however with
+// non-Cartesian (but still affine) cells.
+
+#include "../tests.h"
+
+#include "matrix_vector_rt_face_common.h"
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ const unsigned int n_subdivisions = 2;
+ Point<dim> corners[dim];
+ corners[0] = (dim == 2) ? Point<dim>(1, 0) : Point<dim>(1, 0, 0);
+ corners[1] = (dim == 2) ? Point<dim>(0.5, 0.5) : Point<dim>(0.5, 1, 0.25);
+ if (dim == 3)
+ corners[2] = Point<dim>(0.5, 0, 1);
+ GridGenerator::subdivided_parallelepiped(tria, n_subdivisions, corners);
+
+ FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+
+ AffineConstraints<double> constraints;
+ constraints.close();
+
+ deallog << "Using " << dof.get_fe().get_name() << std::endl;
+ deallog << "Number of cells: " << dof.get_triangulation().n_active_cells()
+ << std::endl;
+ deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl
+ << std::endl;
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::values_gradients);
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::divergence);
+}
--- /dev/null
+
+DEAL:2d::Using FE_RaviartThomasNodal<2>(1)
+DEAL:2d::Number of cells: 4
+DEAL:2d::Number of degrees of freedom: 40
+DEAL:2d::
+DEAL:2d::Testing Values and Gradients
+DEAL:2d::Norm of difference: 6.92435e-16
+DEAL:2d::
+DEAL:2d::Testing Divergence
+DEAL:2d::Norm of difference: 2.84791e-16
+DEAL:2d::
+DEAL:2d::Using FE_RaviartThomasNodal<2>(2)
+DEAL:2d::Number of cells: 4
+DEAL:2d::Number of degrees of freedom: 84
+DEAL:2d::
+DEAL:2d::Testing Values and Gradients
+DEAL:2d::Norm of difference: 4.97871e-15
+DEAL:2d::
+DEAL:2d::Testing Divergence
+DEAL:2d::Norm of difference: 4.05232e-15
+DEAL:2d::
+DEAL:3d::Using FE_RaviartThomasNodal<3>(1)
+DEAL:3d::Number of cells: 8
+DEAL:3d::Number of degrees of freedom: 240
+DEAL:3d::
+DEAL:3d::Testing Values and Gradients
+DEAL:3d::Norm of difference: 2.96549e-15
+DEAL:3d::
+DEAL:3d::Testing Divergence
+DEAL:3d::Norm of difference: 4.91752e-15
+DEAL:3d::
void
test();
+enum TestType : unsigned char
+{
+ values = 0,
+ values_gradients = 1,
+ gradients = 2,
+ divergence = 3
+};
+std::string
+enum_to_string(TestType const enum_type)
+{
+ std::string string_type;
+ switch (enum_type)
+ {
+ case TestType::values:
+ string_type = "Values ";
+ break;
+ case TestType::gradients:
+ string_type = "Gradients ";
+ break;
+ case TestType::values_gradients:
+ string_type = "Values and Gradients ";
+ break;
+ case TestType::divergence:
+ string_type = "Divergence ";
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ break;
+ }
+ return string_type;
+}
template <int dim,
int fe_degree,
class MatrixFreeTest
{
public:
- MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
- : data(data_in){};
+ MatrixFreeTest(const MatrixFree<dim, Number> &data_in,
+ const TestType test_type)
+ : data(data_in)
+ , test_type(test_type)
+ {
+ evaluation_flag =
+ (test_type == TestType::values) ?
+ EvaluationFlags::values :
+ ((test_type == TestType::gradients) ?
+ EvaluationFlags::gradients :
+ ((test_type == TestType::values_gradients) ?
+ EvaluationFlags::values | EvaluationFlags::gradients :
+ EvaluationFlags::gradients));
+ };
virtual ~MatrixFreeTest(){};
FEFaceEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_eval_n(
data, false);
- // Note that this will need to be modified once the Piola transform is
- // implemented
- const unsigned int n_cells =
- data.get_dof_handler().get_triangulation().n_active_cells();
- const Number piola =
- (dim == 2) ? n_cells : Utilities::pow((int)std::cbrt(n_cells), 4);
-
for (unsigned int face = face_range.first; face < face_range.second; ++face)
{
fe_eval.reinit(face);
fe_eval_n.reinit(face);
- fe_eval.gather_evaluate(src,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
- fe_eval_n.gather_evaluate(src,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
-
+ fe_eval.gather_evaluate(src, evaluation_flag);
+ fe_eval_n.gather_evaluate(src, evaluation_flag);
for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
{
- fe_eval.submit_value(Number(10. * piola) * fe_eval.get_value(q), q);
- fe_eval_n.submit_value(Number(10. * piola) * fe_eval_n.get_value(q),
- q);
-
- fe_eval.submit_gradient(Number(piola) * fe_eval.get_gradient(q), q);
- fe_eval_n.submit_gradient(Number(piola) * fe_eval_n.get_gradient(q),
- q);
+ if (test_type < TestType::gradients)
+ {
+ fe_eval.submit_value(10. * fe_eval.get_value(q), q);
+ fe_eval_n.submit_value(10. * fe_eval_n.get_value(q), q);
+ }
+ if (test_type == TestType::gradients ||
+ test_type == TestType::values_gradients)
+ {
+ fe_eval.submit_gradient(fe_eval.get_gradient(q), q);
+ fe_eval_n.submit_gradient(fe_eval_n.get_gradient(q), q);
+ }
+ else if (test_type == TestType::divergence)
+ {
+ fe_eval.submit_divergence(fe_eval.get_divergence(q), q);
+ fe_eval_n.submit_divergence(fe_eval_n.get_divergence(q), q);
+ }
}
-
- fe_eval.integrate_scatter(EvaluationFlags::values |
- EvaluationFlags::gradients,
- dst);
- fe_eval_n.integrate_scatter(EvaluationFlags::values |
- EvaluationFlags::gradients,
- dst);
+ fe_eval.integrate_scatter(evaluation_flag, dst);
+ fe_eval_n.integrate_scatter(evaluation_flag, dst);
}
};
FEFaceEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_eval(data,
true);
- // Note that this will need to be modified once the Piola transform is
- // implemented
- const unsigned int n_cells =
- data.get_dof_handler().get_triangulation().n_active_cells();
- const Number piola =
- (dim == 2) ? n_cells : Utilities::pow((int)std::cbrt(n_cells), 4);
-
for (unsigned int face = face_range.first; face < face_range.second; ++face)
{
fe_eval.reinit(face);
- fe_eval.gather_evaluate(src,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
+ fe_eval.gather_evaluate(src, evaluation_flag);
for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
{
- fe_eval.submit_value(Number(10. * piola) * fe_eval.get_value(q), q);
- fe_eval.submit_gradient(Number(piola) * fe_eval.get_gradient(q), q);
+ if (test_type < TestType::gradients)
+ fe_eval.submit_value(10. * fe_eval.get_value(q), q);
+ if (test_type == TestType::gradients ||
+ test_type == TestType::values_gradients)
+ fe_eval.submit_gradient(fe_eval.get_gradient(q), q);
+ else if (test_type == TestType::divergence)
+ fe_eval.submit_divergence(fe_eval.get_divergence(q), q);
}
- fe_eval.integrate_scatter(EvaluationFlags::values |
- EvaluationFlags::gradients,
- dst);
+ fe_eval.integrate_scatter(evaluation_flag, dst);
}
};
};
protected:
- const MatrixFree<dim, Number> &data;
+ const MatrixFree<dim, Number> & data;
+ EvaluationFlags::EvaluationFlags evaluation_flag;
+ const TestType test_type;
};
template <int dim, int fe_degree, typename Number>
void
do_test(const DoFHandler<dim> & dof,
- const AffineConstraints<double> &constraints)
+ const AffineConstraints<double> &constraints,
+ const TestType test_type)
{
- deallog << "Testing " << dof.get_fe().get_name() << std::endl;
- deallog << "Number of cells: " << dof.get_triangulation().n_active_cells()
- << std::endl;
- deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl
- << std::endl;
-
+ deallog << "Testing " << enum_to_string(test_type) << std::endl;
- // constraints.distribute(solution);
MatrixFree<dim, Number> mf_data;
{
const QGaussLobatto<1> quad(fe_degree + 2);
initial_condition[i] = random_value<Number>();
}
- MatrixFreeTest<dim, fe_degree, fe_degree + 2, Number> mf(mf_data);
+ MatrixFreeTest<dim, fe_degree, fe_degree + 2, Number> mf(mf_data, test_type);
mf.test_functions(solution, initial_condition);
const Tensor<1, dim> phi_i = fe_val[velocities].value(i, q) * 10.;
const Tensor<2, dim> grad_phi_i =
fe_val[velocities].gradient(i, q);
+ const Number div_phi_i = fe_val[velocities].divergence(i, q);
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
const Tensor<1, dim> phi_j = fe_val[velocities].value(j, q);
const Tensor<2, dim> grad_phi_j =
fe_val[velocities].gradient(j, q);
- local_matrix(i, j) +=
- (phi_j * phi_i + scalar_product(grad_phi_i, grad_phi_j)) *
- fe_val.JxW(q);
+ const Number div_phi_j = fe_val[velocities].divergence(j, q);
+
+ if (test_type < TestType::gradients)
+ local_matrix(i, j) += phi_j * phi_i * fe_val.JxW(q);
+ if (test_type == TestType::gradients ||
+ test_type == TestType::values_gradients)
+ local_matrix(i, j) +=
+ scalar_product(grad_phi_i, grad_phi_j) * fe_val.JxW(q);
+ else if (test_type == TestType::divergence)
+ local_matrix(i, j) += div_phi_i * div_phi_j * fe_val.JxW(q);
}
}
cell->get_dof_indices(local_dof_indices);
projects / dealii.git / commitdiff