mkdir build
mkdir c:/project
cd build
- cmake -DCMAKE_BUILD_TYPE=Debug -DCMAKE_INSTALL_PREFIX=c:/project -DDEAL_II_WITH_ZLIB=off -DDEAL_II_CXX_FLAGS="-WX" -T host=x64 -A x64 ..
+ cmake -DCMAKE_BUILD_TYPE=Debug -DCMAKE_INSTALL_PREFIX=c:/project -DDEAL_II_WITH_ZLIB=off -DDEAL_II_CXX_FLAGS="-WX /D FE_EVAL_FACTORY_DEGREE_MAX=2" -T host=x64 -A x64 ..
- name: archive logs
uses: actions/upload-artifact@v1
with:
* FiniteElementData<dim>::degree is higher by one than the constructor
* argument!
*/
+
+namespace internal
+{
+ template <int dim>
+ std::vector<unsigned int>
+ get_lexicographic_numbering_rt_nodal(const unsigned int degree);
+} // namespace internal
+
template <int dim>
class FE_RaviartThomasNodal : public FE_PolyTensor<dim>
{
static const EvaluatorVariant variant = evaluate_evenodd;
};
+ template <bool is_long>
+ struct EvaluatorSelector<MatrixFreeFunctions::tensor_raviart_thomas, is_long>
+ {
+ static const EvaluatorVariant variant = evaluate_raviart_thomas;
+ };
+
/**
const bool add_into_values_array);
};
+ /**
+ * Specialization for MatrixFreeFunctions::tensor_raviart_thomas, which use
+ * specific sum-factorization kernels and with normal/tangential shape_data
+ */
+ template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ struct FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
+ dim,
+ fe_degree,
+ n_q_points_1d,
+ Number>
+ {
+ using EvalNormal =
+ EvaluatorTensorProduct<evaluate_raviart_thomas,
+ dim,
+ (fe_degree == -1) ? 1 : fe_degree + 1,
+ n_q_points_1d,
+ Number>;
+
+ using EvalTangent =
+ EvaluatorTensorProduct<evaluate_raviart_thomas,
+ dim,
+ (fe_degree == -1) ? 1 : fe_degree,
+ n_q_points_1d,
+ Number>;
+ template <bool integrate>
+ static void
+ evaluate_or_integrate(
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ Number * values_dofs_actual,
+ FEEvaluationData<dim, Number, false> & fe_eval,
+ const bool add_into_values_array = false);
+
+ private:
+ template <typename EvalType>
+ static EvalType
+ create_evaluator_tensor_product(
+ const MatrixFreeFunctions::ShapeInfo<Number> &shape_info)
+ {
+ // Shape_data is of size 2 with normal direction in [0] and tangential
+ // direction in [1]. See shape_info.template.h
+ const MatrixFreeFunctions::UnivariateShapeData<Number> *shape_data;
+ if (std::is_same<EvalType, EvalNormal>::value)
+ {
+ shape_data = &shape_info.data[0];
+ }
+ else
+ {
+ shape_data = &shape_info.data[1];
+ }
+ return EvalType(shape_data->shape_values,
+ shape_data->shape_gradients,
+ shape_data->shape_hessians);
+ }
+
+ template <int normal_dir>
+ static void
+ evaluate_tensor_product_per_component(
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ Number * values_dofs_actual,
+ FEEvaluationData<dim, Number, false> & fe_eval,
+ const bool add_into_values_array,
+ std::integral_constant<bool, false>);
+
+ template <int normal_dir>
+ static void
+ evaluate_tensor_product_per_component(
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ Number * values_dofs_actual,
+ FEEvaluationData<dim, Number, false> & fe_eval,
+ const bool add_into_values_array,
+ std::integral_constant<bool, true>);
+ };
+
template <MatrixFreeFunctions::ElementType type,
}
+ template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ template <bool integrate>
+ inline void
+ FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
+ dim,
+ fe_degree,
+ n_q_points_1d,
+ Number>::
+ evaluate_or_integrate(
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ Number * values_dofs_actual,
+ FEEvaluationData<dim, Number, false> & fe_eval,
+ const bool add_into_values_array)
+ {
+ if (evaluation_flag == EvaluationFlags::nothing)
+ return;
+
+ AssertDimension(fe_eval.get_shape_info().data.size(), 2);
+ // First component:
+ evaluate_tensor_product_per_component<0>(
+ evaluation_flag,
+ values_dofs_actual,
+ fe_eval,
+ add_into_values_array,
+ std::integral_constant<bool, integrate>());
+ // Second component :
+ evaluate_tensor_product_per_component<1>(
+ evaluation_flag,
+ values_dofs_actual,
+ fe_eval,
+ add_into_values_array,
+ std::integral_constant<bool, integrate>());
+ if (dim == 3)
+ {
+ // Third component
+ evaluate_tensor_product_per_component<2>(
+ evaluation_flag,
+ values_dofs_actual,
+ fe_eval,
+ add_into_values_array,
+ std::integral_constant<bool, integrate>());
+ }
+ }
+
+ // Helper function that applies the 1d evaluation kernals.
+ // std::integral_constant<bool, false> is the interpolation path, and
+ // std::integral_constant<bool, true> bellow is the integration path.
+ template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ template <int normal_dir>
+ inline void
+ FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
+ dim,
+ fe_degree,
+ n_q_points_1d,
+ Number>::
+ evaluate_tensor_product_per_component(
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ Number * values_dofs_actual,
+ FEEvaluationData<dim, Number, false> & fe_eval,
+ const bool add_into_values_array,
+ std::integral_constant<bool, false>)
+ {
+ (void)add_into_values_array;
+ using Eval0 =
+ typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
+ using Eval1 =
+ typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
+ using Eval2 =
+ typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
+
+ Eval0 eval0 =
+ create_evaluator_tensor_product<Eval0>(fe_eval.get_shape_info());
+ Eval1 eval1 =
+ create_evaluator_tensor_product<Eval1>(fe_eval.get_shape_info());
+ Eval2 eval2 =
+ create_evaluator_tensor_product<Eval2>(fe_eval.get_shape_info());
+
+ Number *temp1 = fe_eval.get_scratch_data().begin();
+ Number *temp2;
+
+ temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
+ fe_eval.get_shape_info().data[0].fe_degree + 1),
+ Utilities::fixed_power<dim>(
+ fe_eval.get_shape_info().data[0].n_q_points_1d));
+
+ const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
+ const std::size_t dofs_per_comp =
+ fe_eval.get_shape_info().dofs_per_component_on_cell;
+
+ // Initial shift depending on component (normal_dir)
+ Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
+ Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
+ Number *gradients_quad =
+ fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
+ Number *hessians_quad =
+ (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
+ fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
+
+ switch (dim)
+ {
+ case 2:
+ if (evaluation_flag & EvaluationFlags::gradients)
+ {
+ eval0.template gradients<0, true, false, normal_dir>(values_dofs,
+ temp1);
+ eval1.template values<1, true, false, normal_dir>(temp1,
+ gradients_quad);
+ }
+ if (evaluation_flag & EvaluationFlags::hessians)
+ {
+ // The evaluation/integration here *should* work, however
+ // the piola transform is not implemented.
+ AssertThrow(false, ExcNotImplemented());
+ // grad xy
+ if (!(evaluation_flag & EvaluationFlags::gradients))
+ eval0.template gradients<0, true, false, normal_dir>(
+ values_dofs, temp1);
+ eval1.template gradients<1, true, false, normal_dir>(
+ temp1, hessians_quad + 2 * n_q_points);
+
+ // grad xx
+ eval0.template hessians<0, true, false, normal_dir>(values_dofs,
+ temp1);
+ eval1.template values<1, true, false, normal_dir>(temp1,
+ hessians_quad);
+ }
+
+ // grad y
+ eval0.template values<0, true, false, normal_dir>(values_dofs, temp1);
+ if (evaluation_flag & EvaluationFlags::gradients)
+ eval1.template gradients<1, true, false, normal_dir>(
+ temp1, gradients_quad + n_q_points);
+
+ // grad yy
+ if (evaluation_flag & EvaluationFlags::hessians)
+ eval1.template hessians<1, true, false, normal_dir>(temp1,
+ hessians_quad +
+ n_q_points);
+
+ // val: can use values applied in x
+ if (evaluation_flag & EvaluationFlags::values)
+ eval1.template values<1, true, false, normal_dir>(temp1,
+ values_quad);
+ break;
+ case 3:
+ if (evaluation_flag & EvaluationFlags::gradients)
+ {
+ // grad x
+ eval0.template gradients<0, true, false, normal_dir>(values_dofs,
+ temp1);
+ eval1.template values<1, true, false, normal_dir>(temp1, temp2);
+ eval2.template values<2, true, false, normal_dir>(temp2,
+ gradients_quad);
+ }
+
+ if (evaluation_flag & EvaluationFlags::hessians)
+ {
+ // The evaluation/integration here *should* work, however
+ // the piola transform is not implemented.
+ AssertThrow(false, ExcNotImplemented());
+ // grad xz
+ if (!(evaluation_flag & EvaluationFlags::gradients))
+ {
+ eval0.template gradients<0, true, false, normal_dir>(
+ values_dofs, temp1);
+ eval1.template values<1, true, false, normal_dir>(temp1,
+ temp2);
+ }
+ eval2.template gradients<2, true, false, normal_dir>(
+ temp2, hessians_quad + 4 * n_q_points);
+
+ // grad xy
+ eval1.template gradients<1, true, false, normal_dir>(temp1,
+ temp2);
+ eval2.template values<2, true, false, normal_dir>(
+ temp2, hessians_quad + 3 * n_q_points);
+
+ // grad xx
+ eval0.template hessians<0, true, false, normal_dir>(values_dofs,
+ temp1);
+ eval1.template values<1, true, false, normal_dir>(temp1, temp2);
+ eval2.template values<2, true, false, normal_dir>(temp2,
+ hessians_quad);
+ }
+
+ // grad y
+ eval0.template values<0, true, false, normal_dir>(values_dofs, temp1);
+ if (evaluation_flag & EvaluationFlags::gradients)
+ {
+ eval1.template gradients<1, true, false, normal_dir>(temp1,
+ temp2);
+ eval2.template values<2, true, false, normal_dir>(temp2,
+ gradients_quad +
+ n_q_points);
+ }
+
+ if (evaluation_flag & EvaluationFlags::hessians)
+ {
+ // grad yz
+ if (!(evaluation_flag & EvaluationFlags::gradients))
+ eval1.template gradients<1, true, false, normal_dir>(temp1,
+ temp2);
+ eval2.template gradients<2, true, false, normal_dir>(
+ temp2, hessians_quad + 5 * n_q_points);
+
+ // grad yy
+ eval1.template hessians<1, true, false, normal_dir>(temp1, temp2);
+ eval2.template values<2, true, false, normal_dir>(temp2,
+ hessians_quad +
+ n_q_points);
+ }
+
+ // grad z: can use the values applied in x direction stored in
+ // temp1
+ eval1.template values<1, true, false, normal_dir>(temp1, temp2);
+ if (evaluation_flag & EvaluationFlags::gradients)
+ eval2.template gradients<2, true, false, normal_dir>(
+ temp2, gradients_quad + 2 * n_q_points);
+
+ // grad zz: can use the values applied in x and y direction stored
+ // in temp2
+ if (evaluation_flag & EvaluationFlags::hessians)
+ eval2.template hessians<2, true, false, normal_dir>(
+ temp2, hessians_quad + 2 * n_q_points);
+
+ // val: can use the values applied in x & y direction stored in
+ // temp2
+ if (evaluation_flag & EvaluationFlags::values)
+ eval2.template values<2, true, false, normal_dir>(temp2,
+ values_quad);
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+ }
+
+ template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ template <int normal_dir>
+ inline void
+ FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
+ dim,
+ fe_degree,
+ n_q_points_1d,
+ Number>::
+ evaluate_tensor_product_per_component(
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ Number * values_dofs_actual,
+ FEEvaluationData<dim, Number, false> & fe_eval,
+ const bool add_into_values_array,
+ std::integral_constant<bool, true>)
+ {
+ using Eval0 =
+ typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
+ using Eval1 =
+ typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
+ using Eval2 =
+ typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
+
+ Eval0 eval0 =
+ create_evaluator_tensor_product<Eval0>(fe_eval.get_shape_info());
+ Eval1 eval1 =
+ create_evaluator_tensor_product<Eval1>(fe_eval.get_shape_info());
+ Eval2 eval2 =
+ create_evaluator_tensor_product<Eval2>(fe_eval.get_shape_info());
+
+ Number *temp1 = fe_eval.get_scratch_data().begin();
+ Number *temp2;
+
+ temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
+ fe_eval.get_shape_info().data[0].fe_degree + 1),
+ Utilities::fixed_power<dim>(
+ fe_eval.get_shape_info().data[0].n_q_points_1d));
+
+ const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
+ const std::size_t dofs_per_comp =
+ fe_eval.get_shape_info().dofs_per_component_on_cell;
+
+ // Initial shift depending on component (normal_dir)
+ Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
+ Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
+ Number *gradients_quad =
+ fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
+ Number *hessians_quad =
+ (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
+ fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
+
+ // Integrate path
+ switch (dim)
+ {
+ case 2:
+ if ((evaluation_flag & EvaluationFlags::values) &&
+ !(evaluation_flag & EvaluationFlags::gradients))
+ {
+ eval1.template values<1, false, false, normal_dir>(values_quad,
+ temp1);
+ if (add_into_values_array == false)
+ eval0.template values<0, false, false, normal_dir>(temp1,
+ values_dofs);
+ else
+ eval0.template values<0, false, true, normal_dir>(temp1,
+ values_dofs);
+ }
+ if (evaluation_flag & EvaluationFlags::gradients)
+ {
+ eval1.template gradients<1, false, false, normal_dir>(
+ gradients_quad + n_q_points, temp1);
+ if ((evaluation_flag & EvaluationFlags::values))
+ eval1.template values<1, false, true, normal_dir>(values_quad,
+ temp1);
+ if (add_into_values_array == false)
+ eval0.template values<0, false, false, normal_dir>(temp1,
+ values_dofs);
+ else
+ eval0.template values<0, false, true, normal_dir>(temp1,
+ values_dofs);
+ eval1.template values<1, false, false, normal_dir>(gradients_quad,
+ temp1);
+ eval0.template gradients<0, false, true, normal_dir>(temp1,
+ values_dofs);
+ }
+ if (evaluation_flag & EvaluationFlags::hessians)
+ {
+ // grad xx
+ eval1.template values<1, false, false, normal_dir>(hessians_quad,
+ temp1);
+
+ if ((evaluation_flag & EvaluationFlags::values) ||
+ (evaluation_flag & EvaluationFlags::gradients) ||
+ add_into_values_array == true)
+ eval0.template hessians<0, false, true, normal_dir>(
+ temp1, values_dofs);
+ else
+ eval0.template hessians<0, false, false, normal_dir>(
+ temp1, values_dofs);
+
+ // grad yy
+ eval1.template hessians<1, false, false, normal_dir>(
+ hessians_quad + n_q_points, temp1);
+ eval0.template values<0, false, true, normal_dir>(temp1,
+ values_dofs);
+
+ // grad xy
+ eval1.template gradients<1, false, false, normal_dir>(
+ hessians_quad + 2 * n_q_points, temp1);
+ eval0.template gradients<0, false, true, normal_dir>(temp1,
+ values_dofs);
+ }
+ break;
+
+ case 3:
+ if ((evaluation_flag & EvaluationFlags::values) &&
+ !(evaluation_flag & EvaluationFlags::gradients))
+ {
+ eval2.template values<2, false, false, normal_dir>(values_quad,
+ temp1);
+ eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+ if (add_into_values_array == false)
+ eval0.template values<0, false, false, normal_dir>(temp2,
+ values_dofs);
+ else
+ eval0.template values<0, false, true, normal_dir>(temp2,
+ values_dofs);
+ }
+ if (evaluation_flag & EvaluationFlags::gradients)
+ {
+ eval2.template gradients<2, false, false, normal_dir>(
+ gradients_quad + 2 * n_q_points, temp1);
+ if ((evaluation_flag & EvaluationFlags::values))
+ eval2.template values<2, false, true, normal_dir>(values_quad,
+ temp1);
+ eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+ eval2.template values<2, false, false, normal_dir>(
+ gradients_quad + n_q_points, temp1);
+ eval1.template gradients<1, false, true, normal_dir>(temp1,
+ temp2);
+ if (add_into_values_array == false)
+ eval0.template values<0, false, false, normal_dir>(temp2,
+ values_dofs);
+ else
+ eval0.template values<0, false, true, normal_dir>(temp2,
+ values_dofs);
+ eval2.template values<2, false, false, normal_dir>(gradients_quad,
+ temp1);
+ eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+ eval0.template gradients<0, false, true, normal_dir>(temp2,
+ values_dofs);
+ }
+ if (evaluation_flag & EvaluationFlags::hessians)
+ {
+ // grad xx
+ eval2.template values<2, false, false, normal_dir>(hessians_quad,
+ temp1);
+ eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+
+ if ((evaluation_flag & EvaluationFlags::values) ||
+ (evaluation_flag & EvaluationFlags::gradients) ||
+ add_into_values_array == true)
+ eval0.template hessians<0, false, true, normal_dir>(
+ temp2, values_dofs);
+ else
+ eval0.template hessians<0, false, false, normal_dir>(
+ temp2, values_dofs);
+
+ // grad yy
+ eval2.template values<2, false, false, normal_dir>(hessians_quad +
+ n_q_points,
+ temp1);
+ eval1.template hessians<1, false, false, normal_dir>(temp1,
+ temp2);
+ eval0.template values<0, false, true, normal_dir>(temp2,
+ values_dofs);
+
+ // grad zz
+ eval2.template hessians<2, false, false, normal_dir>(
+ hessians_quad + 2 * n_q_points, temp1);
+ eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+ eval0.template values<0, false, true, normal_dir>(temp2,
+ values_dofs);
+
+ // grad xy
+ eval2.template values<2, false, false, normal_dir>(
+ hessians_quad + 3 * n_q_points, temp1);
+ eval1.template gradients<1, false, false, normal_dir>(temp1,
+ temp2);
+ eval0.template gradients<0, false, true, normal_dir>(temp2,
+ values_dofs);
+
+ // grad xz
+ eval2.template gradients<2, false, false, normal_dir>(
+ hessians_quad + 4 * n_q_points, temp1);
+ eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+ eval0.template gradients<0, false, true, normal_dir>(temp2,
+ values_dofs);
+
+ // grad yz
+ eval2.template gradients<2, false, false, normal_dir>(
+ hessians_quad + 5 * n_q_points, temp1);
+ eval1.template gradients<1, false, false, normal_dir>(temp1,
+ temp2);
+ eval0.template values<0, false, true, normal_dir>(temp2,
+ values_dofs);
+ }
+
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+ }
/**
* This struct implements the change between two different bases. This is an
Assert(fe_eval.get_shape_info().data.size() == 1 ||
(fe_eval.get_shape_info().data.size() == dim &&
- element_type == ElementType::tensor_general),
+ element_type == ElementType::tensor_general) ||
+ element_type == ElementType::tensor_raviart_thomas,
ExcNotImplemented());
if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
values_dofs,
fe_eval);
}
+ else if (element_type == ElementType::tensor_raviart_thomas)
+ {
+ FEEvaluationImpl<
+ ElementType::tensor_raviart_thomas,
+ dim,
+ (fe_degree == -1) ? 1 : fe_degree,
+ (n_q_points_1d < 1) ? 1 : n_q_points_1d,
+ Number>::template evaluate_or_integrate<false>(evaluation_flag,
+ const_cast<Number *>(
+ values_dofs),
+ fe_eval);
+ }
else
{
FEEvaluationImpl<ElementType::tensor_general,
Assert(fe_eval.get_shape_info().data.size() == 1 ||
(fe_eval.get_shape_info().data.size() == dim &&
- element_type == ElementType::tensor_general),
+ element_type == ElementType::tensor_general) ||
+ element_type == ElementType::tensor_raviart_thomas,
ExcNotImplemented());
if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
fe_eval,
sum_into_values_array);
}
+ else if (element_type == ElementType::tensor_raviart_thomas)
+ {
+ FEEvaluationImpl<ElementType::tensor_raviart_thomas,
+ dim,
+ (fe_degree == -1) ? 1 : fe_degree,
+ (n_q_points_1d < 1) ? 1 : n_q_points_1d,
+ Number>::
+ template evaluate_or_integrate<true>(integration_flag,
+ const_cast<Number *>(
+ values_dofs),
+ fe_eval,
+ sum_into_values_array);
+ }
else
{
FEEvaluationImpl<ElementType::tensor_general,
}
};
+ template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ struct FEFaceEvaluationImplRaviartThomas
+ {
+ using EvalNormal =
+ EvaluatorTensorProduct<evaluate_raviart_thomas,
+ dim - 1,
+ (fe_degree == -1) ? 1 : fe_degree + 1,
+ n_q_points_1d,
+ Number>;
+ using EvalTangent =
+ EvaluatorTensorProduct<evaluate_raviart_thomas,
+ dim - 1,
+ (fe_degree == -1) ? 1 : fe_degree,
+ n_q_points_1d,
+ Number>;
+
+ using EvalGeneral = EvaluatorTensorProduct<evaluate_general,
+ dim - 1,
+ fe_degree,
+ n_q_points_1d,
+ Number>;
+ template <typename EvalType>
+ static EvalType
+ create_evaluator_tensor_product(
+ const MatrixFreeFunctions::UnivariateShapeData<Number> &data,
+ const unsigned int subface_index,
+ const unsigned int direction)
+ {
+ if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
+ return EvalType(data.shape_values,
+ data.shape_gradients,
+ data.shape_hessians);
+ else
+ {
+ Assert(false, ExcNotImplemented());
+
+ const unsigned int index =
+ direction == 0 ? subface_index % 2 : subface_index / 2;
+ return EvalType(data.values_within_subface[index],
+ data.gradients_within_subface[index],
+ data.hessians_within_subface[index]);
+ }
+ }
+
+ 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)
+ {
+ // 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<-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, 2>(
+ eval_general,
+ eval_general,
+ values_dofs,
+ fe_eval,
+ scratch_data,
+ evaluation_flag,
+ n_dofs_normal,
+ 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.
+ */
+ template <int normal_dir, int component, 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,
+ std::integral_constant<bool, false>)
+ {
+ constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
+
+ Number *values_quad = fe_eval.begin_values() + n_q_points * component;
+ Number *gradients_quad =
+ fe_eval.begin_gradients() + dim * n_q_points * component;
+ Number *hessians_quad =
+ fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * component;
+
+ // Evaluation path
+
+ if ((evaluation_flag & EvaluationFlags::values) &&
+ !(evaluation_flag & EvaluationFlags::gradients))
+ {
+ switch (dim)
+ {
+ case 3:
+ eval0.template values<0, true, false, normal_dir>(values_dofs,
+ values_quad);
+ eval1.template values<1, true, false, normal_dir>(values_quad,
+ values_quad);
+ break;
+ case 2:
+ eval0.template values<0, true, false, normal_dir>(values_dofs,
+ values_quad);
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ else if (evaluation_flag & EvaluationFlags::gradients)
+ {
+ switch (dim)
+ {
+ case 3:
+ // grad x
+ eval0.template gradients<0, true, false, normal_dir>(
+ values_dofs, scratch_data);
+ eval1.template values<1, true, false, normal_dir>(
+ scratch_data, gradients_quad);
+
+ // grad y
+ eval0.template values<0, true, false, normal_dir>(values_dofs,
+ scratch_data);
+ eval1.template gradients<1, true, false, normal_dir>(
+ scratch_data, gradients_quad + n_q_points);
+
+ if (evaluation_flag & EvaluationFlags::values)
+ eval1.template values<1, true, false, normal_dir>(
+ scratch_data, values_quad);
+
+ // grad z
+ eval0.template values<0, true, false, normal_dir>(values_dofs +
+ dofs_stride,
+ scratch_data);
+ eval1.template values<1, true, false, normal_dir>(
+ scratch_data, gradients_quad + 2 * n_q_points);
+
+ break;
+ case 2:
+ eval0.template values<0, true, false, normal_dir>(
+ values_dofs + dofs_stride, gradients_quad + n_q_points);
+ eval0.template gradients<0, true, false, normal_dir>(
+ values_dofs, gradients_quad);
+ if ((evaluation_flag & EvaluationFlags::values))
+ eval0.template values<0, true, false, normal_dir>(
+ values_dofs, values_quad);
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+ }
+
+ if (evaluation_flag & EvaluationFlags::hessians)
+ {
+ switch (dim)
+ {
+ case 3:
+ // grad xx
+ eval0.template hessians<0, true, false, normal_dir>(
+ values_dofs, scratch_data);
+ eval1.template values<1, true, false, normal_dir>(
+ scratch_data, hessians_quad);
+
+ // grad yy
+ eval0.template values<0, true, false, normal_dir>(values_dofs,
+ scratch_data);
+ eval1.template hessians<1, true, false, normal_dir>(
+ scratch_data, hessians_quad + n_q_points);
+
+ // grad zz
+ eval0.template values<0, true, false, normal_dir>(
+ values_dofs + 2 * dofs_stride, scratch_data);
+ eval1.template values<1, true, false, normal_dir>(
+ scratch_data, hessians_quad + 2 * n_q_points);
+
+ // grad xy
+ eval0.template gradients<0, true, false, normal_dir>(
+ values_dofs, scratch_data);
+ eval1.template gradients<1, true, false, normal_dir>(
+ scratch_data, hessians_quad + 3 * n_q_points);
+
+ // grad xz
+ eval0.template gradients<0, true, false, normal_dir>(
+ values_dofs + dofs_stride, scratch_data);
+ eval1.template values<1, true, false, normal_dir>(
+ scratch_data, hessians_quad + 4 * n_q_points);
+
+ // grad yz
+ eval0.template values<0, true, false, normal_dir>(values_dofs +
+ dofs_stride,
+ scratch_data);
+ eval1.template gradients<1, true, false, normal_dir>(
+ scratch_data, hessians_quad + 5 * n_q_points);
+
+ break;
+ case 2:
+ // grad xx
+ eval0.template hessians<0, true, false, normal_dir>(
+ values_dofs, hessians_quad);
+ // grad yy
+ eval0.template values<0, true, false, normal_dir>(
+ values_dofs + 2 * dofs_stride, hessians_quad + n_q_points);
+ // grad xy
+ eval0.template gradients<0, true, false, normal_dir>(
+ values_dofs + dofs_stride, hessians_quad + 2 * n_q_points);
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+ }
+ }
+
+ template <int normal_dir, int component, 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,
+ std::integral_constant<bool, true>)
+ {
+ constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
+
+ Number *values_quad = fe_eval.begin_values() + n_q_points * component;
+ Number *gradients_quad =
+ fe_eval.begin_gradients() + dim * n_q_points * component;
+ Number *hessians_quad =
+ fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * component;
+
+ // Integration path
+ if ((evaluation_flag & EvaluationFlags::values) &&
+ !(evaluation_flag & EvaluationFlags::gradients))
+ {
+ switch (dim)
+ {
+ case 3:
+ eval1.template values<1, false, false, normal_dir>(values_quad,
+ values_quad);
+ eval0.template values<0, false, false, normal_dir>(values_quad,
+ values_dofs);
+ break;
+ case 2:
+ eval0.template values<0, false, false, normal_dir>(values_quad,
+ values_dofs);
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ else if (evaluation_flag & EvaluationFlags::gradients)
+ {
+ switch (dim)
+ {
+ case 3:
+ // grad z
+ eval1.template values<1, false, false, normal_dir>(
+ gradients_quad + 2 * n_q_points,
+ gradients_quad + 2 * n_q_points);
+ eval0.template values<0, false, false, normal_dir>(
+ gradients_quad + 2 * n_q_points, values_dofs + dofs_stride);
+
+ if (evaluation_flag & EvaluationFlags::values)
+ {
+ eval1.template values<1, false, false, normal_dir>(
+ values_quad, scratch_data);
+ eval1.template gradients<1, false, true, normal_dir>(
+ gradients_quad + n_q_points, scratch_data);
+ }
+ else
+ eval1.template gradients<1, false, false, normal_dir>(
+ gradients_quad + n_q_points, scratch_data);
+
+ // grad y
+ eval0.template values<0, false, false, normal_dir>(scratch_data,
+ values_dofs);
+
+ // grad x
+ eval1.template values<1, false, false, normal_dir>(
+ gradients_quad, scratch_data);
+ eval0.template gradients<0, false, true, normal_dir>(
+ scratch_data, values_dofs);
+
+ break;
+ case 2:
+ eval0.template values<0, false, false, normal_dir>(
+ gradients_quad + n_q_points, values_dofs + dofs_stride);
+ eval0.template gradients<0, false, false, normal_dir>(
+ gradients_quad, values_dofs);
+ if (evaluation_flag & EvaluationFlags::values)
+ eval0.template values<0, false, true, normal_dir>(
+ values_quad, values_dofs);
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+ }
+
+ if (evaluation_flag & EvaluationFlags::hessians)
+ {
+ switch (dim)
+ {
+ case 3:
+ // grad xx
+ eval1.template values<1, false, false, normal_dir>(
+ hessians_quad, scratch_data);
+ if ((evaluation_flag &
+ (EvaluationFlags::values | EvaluationFlags::gradients)))
+ eval0.template hessians<0, false, true, normal_dir>(
+ scratch_data, values_dofs);
+ else
+ eval0.template hessians<0, false, false, normal_dir>(
+ scratch_data, values_dofs);
+
+ // grad yy
+ eval1.template hessians<1, false, false, normal_dir>(
+ hessians_quad + n_q_points, scratch_data);
+ eval0.template values<0, false, true, normal_dir>(scratch_data,
+ values_dofs);
+
+ // grad zz
+ eval1.template values<1, false, false, normal_dir>(
+ hessians_quad + 2 * n_q_points, scratch_data);
+ eval0.template values<0, false, false, normal_dir>(
+ scratch_data, values_dofs + 2 * dofs_stride);
+
+ // grad xy
+ eval1.template gradients<1, false, false, normal_dir>(
+ hessians_quad + 3 * n_q_points, scratch_data);
+ eval0.template gradients<0, false, true, normal_dir>(
+ scratch_data, values_dofs);
+
+ // grad xz
+ eval1.template values<1, false, false, normal_dir>(
+ hessians_quad + 4 * n_q_points, scratch_data);
+ if ((evaluation_flag & EvaluationFlags::gradients))
+ eval0.template gradients<0, false, true, normal_dir>(
+ scratch_data, values_dofs + dofs_stride);
+ else
+ eval0.template gradients<0, false, false, normal_dir>(
+ scratch_data, values_dofs + dofs_stride);
+
+ // grad yz
+ eval1.template gradients<1, false, false, normal_dir>(
+ hessians_quad + 5 * n_q_points, scratch_data);
+ eval0.template values<0, false, true, normal_dir>(
+ scratch_data, values_dofs + dofs_stride);
+
+ break;
+ case 2:
+ // grad xx
+ if (evaluation_flag &
+ (EvaluationFlags::values | EvaluationFlags::gradients))
+ eval0.template hessians<0, false, true, normal_dir>(
+ hessians_quad, values_dofs);
+ else
+ eval0.template hessians<0, false, false, normal_dir>(
+ hessians_quad, values_dofs);
+
+ // grad yy
+ eval0.template values<0, false, false, normal_dir>(
+ hessians_quad + n_q_points, values_dofs + 2 * dofs_stride);
+ // grad xy
+ if ((evaluation_flag & EvaluationFlags::gradients))
+ eval0.template gradients<0, false, true, normal_dir>(
+ hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
+ else
+ eval0.template gradients<0, false, false, normal_dir>(
+ hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+ }
+ }
+ };
template <int dim, int fe_degree, typename Number, bool lex_faces = false>
struct FEFaceNormalEvaluationImpl
{
+ using EvalNormal =
+ EvaluatorTensorProduct<evaluate_raviart_thomas,
+ dim,
+ (fe_degree == -1) ? 1 : fe_degree + 1,
+ 0,
+ Number>;
+ using EvalTangent =
+ EvaluatorTensorProduct<evaluate_raviart_thomas,
+ dim,
+ (fe_degree == -1) ? 1 : fe_degree,
+ 0,
+ Number>;
+
template <bool do_evaluate, bool add_into_output>
static void
interpolate(const unsigned int n_components,
shape_info.data.front().fe_degree ||
fe_degree == -1,
ExcInternalError());
-
- interpolate_generic<do_evaluate, add_into_output>(
- n_components,
- input,
- output,
- flags,
- face_no,
- shape_info.data.front().fe_degree + 1,
- shape_info.data.front().shape_data_on_face,
- shape_info.dofs_per_component_on_cell,
- 3 * shape_info.dofs_per_component_on_face);
+ if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
+ interpolate_generic_raviart_thomas<do_evaluate, add_into_output>(
+ n_components, input, output, flags, face_no, shape_info);
+ else
+ interpolate_generic<do_evaluate, add_into_output>(
+ n_components,
+ input,
+ output,
+ flags,
+ face_no,
+ shape_info.data.front().fe_degree + 1,
+ shape_info.data.front().shape_data_on_face,
+ shape_info.dofs_per_component_on_cell,
+ 3 * shape_info.dofs_per_component_on_face);
}
/**
dofs_per_component_on_face);
}
}
+
+ private:
+ template <typename EvalType>
+ static EvalType
+ create_evaluator_tensor_product(
+ const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const unsigned int face_no)
+ {
+ // Shape_data is of size 2 with normal direction in [0] and tangential
+ // direction in [1]. See shape_info.template.h
+ const MatrixFreeFunctions::UnivariateShapeData<Number> *shape_data;
+ if (std::is_same<EvalType, EvalNormal>::value)
+ {
+ shape_data = &shape_info.data[0];
+ }
+ else
+ {
+ shape_data = &shape_info.data[1];
+ }
+ return EvalType(shape_data->shape_data_on_face[face_no % 2],
+ AlignedVector<Number>(),
+ AlignedVector<Number>());
+ }
+
+ template <bool do_evaluate,
+ bool add_into_output,
+ int face_direction = 0,
+ int max_derivative = 0>
+ static void
+ interpolate_generic_raviart_thomas(
+ const unsigned int n_components,
+ const Number * input,
+ Number * output,
+ const EvaluationFlags::EvaluationFlags flag,
+ const unsigned int face_no,
+ const MatrixFreeFunctions::ShapeInfo<Number> &shape_info)
+ {
+ bool increase_max_der = false;
+ if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
+ (flag & EvaluationFlags::gradients && max_derivative < 1))
+ increase_max_der = true;
+
+ if (face_direction == face_no / 2 && !increase_max_der)
+ {
+ interpolate_generic_raviart_thomas_apply_face<do_evaluate,
+ add_into_output,
+ face_direction,
+ max_derivative>(
+ shape_info, face_no, input, output);
+ }
+ else if (face_direction == face_no / 2)
+ {
+ // Only increase max_derivative
+ interpolate_generic_raviart_thomas<do_evaluate,
+ add_into_output,
+ face_direction,
+ std::min(max_derivative + 1, 2)>(
+ n_components, input, output, flag, face_no, shape_info);
+ }
+ else if (face_direction < dim)
+ {
+ if (increase_max_der)
+ {
+ interpolate_generic_raviart_thomas<
+ do_evaluate,
+ add_into_output,
+ std::min(face_direction + 1, dim - 1),
+ std::min(max_derivative + 1, 2)>(
+ n_components, input, output, flag, face_no, shape_info);
+ }
+ else
+ {
+ interpolate_generic_raviart_thomas<do_evaluate,
+ add_into_output,
+ std::min(face_direction + 1,
+ dim - 1),
+ max_derivative>(
+ n_components, input, output, flag, face_no, shape_info);
+ }
+ }
+ }
+
+ /* Help function for interpolate_generic_raviart_thomas */
+ template <bool do_evaluate,
+ bool add_into_output,
+ int face_direction,
+ int max_derivative>
+ static inline void
+ interpolate_generic_raviart_thomas_apply_face(
+ const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const unsigned int face_no,
+ const Number * input,
+ Number * output)
+ {
+ using Evalf0 = typename std::
+ conditional<face_direction == 0, EvalNormal, EvalTangent>::type;
+ using Evalf1 = typename std::
+ conditional<face_direction == 1, EvalNormal, EvalTangent>::type;
+ using Evalf2 = typename std::
+ conditional<face_direction == 2, EvalNormal, EvalTangent>::type;
+
+ Evalf0 evalf0 =
+ create_evaluator_tensor_product<Evalf0>(shape_info, face_no);
+ Evalf1 evalf1 =
+ create_evaluator_tensor_product<Evalf1>(shape_info, face_no);
+ Evalf2 evalf2 =
+ create_evaluator_tensor_product<Evalf2>(shape_info, face_no);
+
+ const unsigned int dofs_per_component_on_cell =
+ shape_info.dofs_per_component_on_cell;
+ const unsigned int dofs_per_component_on_face =
+ 3 * shape_info.dofs_per_component_on_face;
+
+ // NOTE! dofs_per_component_on_face is in the tangent direction,
+ // i.e (fe.degree+1)*fe.degree. Normal faces are only
+ // fe.degree*fe.degree
+ const unsigned int in_stride =
+ do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face;
+ const unsigned int out_stride =
+ do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell;
+
+ const unsigned int in_stride_after_normal =
+ do_evaluate ?
+ dofs_per_component_on_cell :
+ dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2);
+ const unsigned int out_stride_after_normal =
+ do_evaluate ?
+ dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2) :
+ dofs_per_component_on_cell;
+
+ evalf0.template apply_face<face_direction,
+ do_evaluate,
+ add_into_output,
+ max_derivative,
+ lex_faces,
+ 0>(input, output);
+ // stride to next component
+ input += (face_direction == 0) ? in_stride_after_normal : in_stride;
+ output += (face_direction == 0) ? out_stride_after_normal : out_stride;
+
+ evalf1.template apply_face<face_direction,
+ do_evaluate,
+ add_into_output,
+ max_derivative,
+ lex_faces,
+ 1>(input, output);
+
+ if (dim == 3)
+ {
+ // stride to next component
+ input += (face_direction == 1) ? in_stride_after_normal : in_stride;
+ output +=
+ (face_direction == 1) ? out_stride_after_normal : out_stride;
+
+ evalf2.template apply_face<face_direction,
+ do_evaluate,
+ add_into_output,
+ max_derivative,
+ lex_faces,
+ 2>(input, output);
+ }
+ }
};
constexpr unsigned int n_q_points_1d_actual =
fe_degree > -1 ? n_q_points_1d : 0;
- if (fe_degree > -1 &&
- subface_index >= GeometryInfo<dim>::max_children_per_cell &&
- shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
+ if (fe_degree >= 1 &&
+ shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
+ {
+ FEFaceEvaluationImplRaviartThomas<dim,
+ (fe_degree == -1) ? 1 : fe_degree,
+ (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,
+ face_no);
+ }
+ else if (fe_degree > -1 &&
+ subface_index >= GeometryInfo<dim>::max_children_per_cell &&
+ shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
FEFaceEvaluationImpl<true,
dim,
fe_degree,
fe_degree > -1 ? n_q_points_1d : 0;
const unsigned int subface_index = fe_eval.get_subface_index();
- if (fe_degree > -1 &&
- fe_eval.get_subface_index() >=
- GeometryInfo<dim - 1>::max_children_per_cell &&
- shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
+ if (fe_degree >= 1 &&
+ shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
+ {
+ FEFaceEvaluationImplRaviartThomas<dim,
+ (fe_degree == -1) ? 1 : fe_degree,
+ (n_q_points_1d < 1) ? 1 :
+ n_q_points_1d,
+ Number>::
+ template evaluate_or_integrate_in_face<true>(integration_flag,
+ shape_info,
+ temp,
+ fe_eval,
+ scratch_data,
+ subface_index,
+ face_no);
+ }
+ else if (fe_degree > -1 &&
+ fe_eval.get_subface_index() >=
+ GeometryInfo<dim - 1>::max_children_per_cell &&
+ shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
FEFaceEvaluationImpl<
true,
dim,
--- /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.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_evaluation_template_factory_templates_h
+#define dealii_matrix_free_evaluation_template_factory_templates_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/matrix_free/evaluation_kernels.h>
+#include <deal.II/matrix_free/evaluation_selector.h>
+#include <deal.II/matrix_free/evaluation_template_factory.h>
+#include <deal.II/matrix_free/evaluation_template_factory_internal.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace internal
+{
+ template <int dim, typename Number>
+ void
+ FEFaceEvaluationFactory<dim, Number>::evaluate(
+ const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ const Number * values_dofs,
+ FEEvaluationData<dim, Number, true> & fe_eval)
+ {
+ instantiation_helper_run<1,
+ FEFaceEvaluationImplEvaluateSelector<dim, Number>>(
+ fe_eval.get_shape_info().data[0].fe_degree,
+ fe_eval.get_shape_info().data[0].n_q_points_1d,
+ n_components,
+ evaluation_flag,
+ values_dofs,
+ fe_eval);
+ }
+
+
+
+ template <int dim, typename Number>
+ void
+ FEFaceEvaluationFactory<dim, Number>::integrate(
+ const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags integration_flag,
+ Number * values_dofs,
+ FEEvaluationData<dim, Number, true> & fe_eval)
+ {
+ instantiation_helper_run<
+ 1,
+ FEFaceEvaluationImplIntegrateSelector<dim, Number>>(
+ fe_eval.get_shape_info().data[0].fe_degree,
+ fe_eval.get_shape_info().data[0].n_q_points_1d,
+ n_components,
+ integration_flag,
+ values_dofs,
+ fe_eval);
+ }
+
+
+
+ template <int dim, typename Number, typename VectorizedArrayType>
+ void
+ FEFaceEvaluationGatherFactory<dim, Number, VectorizedArrayType>::evaluate(
+ const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ const Number * src_ptr,
+ const std::vector<ArrayView<const Number>> * sm_ptr,
+ FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
+ {
+ instantiation_helper_run<
+ 1,
+ FEFaceEvaluationImplGatherEvaluateSelector<dim,
+ Number,
+ VectorizedArrayType>>(
+ fe_eval.get_shape_info().data[0].fe_degree,
+ fe_eval.get_shape_info().data[0].n_q_points_1d,
+ n_components,
+ evaluation_flag,
+ src_ptr,
+ sm_ptr,
+ fe_eval);
+ }
+
+
+
+ template <int dim, typename Number, typename VectorizedArrayType>
+ void
+ FEFaceEvaluationGatherFactory<dim, Number, VectorizedArrayType>::integrate(
+ const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags integration_flag,
+ Number * dst_ptr,
+ const std::vector<ArrayView<const Number>> * sm_ptr,
+ FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
+ {
+ instantiation_helper_run<
+ 1,
+ FEFaceEvaluationImplIntegrateScatterSelector<dim,
+ Number,
+ VectorizedArrayType>>(
+ fe_eval.get_shape_info().data[0].fe_degree,
+ fe_eval.get_shape_info().data[0].n_q_points_1d,
+ n_components,
+ integration_flag,
+ dst_ptr,
+ sm_ptr,
+ fe_eval);
+ }
+
+} // end of namespace internal
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
- template <int dim, typename Number>
- void
- FEFaceEvaluationFactory<dim, Number>::evaluate(
- const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- const Number * values_dofs,
- FEEvaluationData<dim, Number, true> & fe_eval)
- {
- instantiation_helper_run<1,
- FEFaceEvaluationImplEvaluateSelector<dim, Number>>(
- fe_eval.get_shape_info().data[0].fe_degree,
- fe_eval.get_shape_info().data[0].n_q_points_1d,
- n_components,
- evaluation_flag,
- values_dofs,
- fe_eval);
- }
-
-
-
- template <int dim, typename Number>
- void
- FEFaceEvaluationFactory<dim, Number>::integrate(
- const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags integration_flag,
- Number * values_dofs,
- FEEvaluationData<dim, Number, true> & fe_eval)
- {
- instantiation_helper_run<
- 1,
- FEFaceEvaluationImplIntegrateSelector<dim, Number>>(
- fe_eval.get_shape_info().data[0].fe_degree,
- fe_eval.get_shape_info().data[0].n_q_points_1d,
- n_components,
- integration_flag,
- values_dofs,
- fe_eval);
- }
-
-
-
- template <int dim, typename Number, typename VectorizedArrayType>
- void
- FEFaceEvaluationGatherFactory<dim, Number, VectorizedArrayType>::evaluate(
- const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- const Number * src_ptr,
- const std::vector<ArrayView<const Number>> * sm_ptr,
- FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
- {
- instantiation_helper_run<
- 1,
- FEFaceEvaluationImplGatherEvaluateSelector<dim,
- Number,
- VectorizedArrayType>>(
- fe_eval.get_shape_info().data[0].fe_degree,
- fe_eval.get_shape_info().data[0].n_q_points_1d,
- n_components,
- evaluation_flag,
- src_ptr,
- sm_ptr,
- fe_eval);
- }
-
-
-
- template <int dim, typename Number, typename VectorizedArrayType>
- void
- FEFaceEvaluationGatherFactory<dim, Number, VectorizedArrayType>::integrate(
- const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags integration_flag,
- Number * dst_ptr,
- const std::vector<ArrayView<const Number>> * sm_ptr,
- FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
- {
- instantiation_helper_run<
- 1,
- FEFaceEvaluationImplIntegrateScatterSelector<dim,
- Number,
- VectorizedArrayType>>(
- fe_eval.get_shape_info().data[0].fe_degree,
- fe_eval.get_shape_info().data[0].n_q_points_1d,
- n_components,
- integration_flag,
- dst_ptr,
- sm_ptr,
- fe_eval);
- }
-
-
-
template <int dim, typename Number>
void
CellwiseInverseMassFactory<dim, Number>::apply(
dof_info[no].fe_index_conversion[fe_index].clear();
for (unsigned int c = 0; c < dof_info[no].n_base_elements; ++c)
{
- dof_info[no].n_components[c] = fe.element_multiplicity(c);
+ dof_info[no].n_components[c] =
+ fe.element_multiplicity(c) *
+ fe.base_element(c).n_components();
for (unsigned int l = 0; l < dof_info[no].n_components[c]; ++l)
{
dof_info[no].component_to_base_index.push_back(c);
dof_info[no]
.component_dof_indices_offset[fe_index]
- .push_back(dof_info[no]
- .component_dof_indices_offset[fe_index]
- .back() +
- fe.base_element(c).n_dofs_per_cell());
+ .push_back(
+ dof_info[no]
+ .component_dof_indices_offset[fe_index]
+ .back() +
+ shape_infos(dof_info[no].global_base_element_offset + c,
+ fe_index)
+ .dofs_per_component_on_cell);
dof_info[no].fe_index_conversion[fe_index].push_back(
fe.base_element(c).degree);
}
tensor_symmetric_plus_dg0 = 5,
/**
- * Shape functions without an tensor product properties.
+ * Special case of the FE_RaviartThomasNodal element with anisotropic
+ * tensor product shape functions, i.e degree (k + 1) in normal direction,
+ * and k in tangential direction.
*/
- tensor_none = 6
+ tensor_raviart_thomas = 6,
+
+ /**
+ * Shape functions without a tensor product properties.
+ */
+ tensor_none = 7
+
+
};
#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_tools.h>
}
}
-
-
template <int dim_to, int dim, int spacedim>
std::unique_ptr<FiniteElement<dim_to, dim_to>>
create_fe(const FiniteElement<dim, spacedim> &fe)
}
-
// ----------------- actual ShapeInfo implementation --------------------
template <typename Number>
{
static_assert(dim == spacedim,
"Currently, only the case dim=spacedim is implemented");
- if (quad_in.is_tensor_product() == false ||
- dynamic_cast<const FE_SimplexP<dim> *>(
- &fe_in.base_element(base_element_number)) ||
- dynamic_cast<const FE_SimplexDGP<dim> *>(
- &fe_in.base_element(base_element_number)) ||
- dynamic_cast<const FE_WedgeP<dim> *>(
- &fe_in.base_element(base_element_number)) ||
- dynamic_cast<const FE_PyramidP<dim> *>(
+
+ // ShapeInfo for RT elements. Here, data is of size 2 instead of 1.
+ // data[0] is univariate_shape_data in normal direction and
+ // data[1] is univariate_shape_data in tangential direction
+ //
+ if (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
&fe_in.base_element(base_element_number)))
+ {
+ element_type = tensor_raviart_thomas;
+
+ const auto quad = quad_in.get_tensor_basis()[0];
+
+ const FiniteElement<dim> &fe =
+ fe_in.base_element(base_element_number);
+ n_dimensions = dim;
+ n_components = fe_in.n_components();
+
+ data.resize(2);
+ const unsigned int n_q_points_1d = quad.size();
+
+ n_q_points = Utilities::fixed_power<dim>(n_q_points_1d);
+ n_q_points_face = Utilities::fixed_power<dim - 1>(n_q_points_1d);
+
+ dofs_per_component_on_cell = fe_in.n_dofs_per_cell() / n_components;
+
+ // NOTE dofs_per_component_on_face is in tangential direction!
+ dofs_per_component_on_face =
+ fe_in.n_dofs_per_face() + Utilities::pow(fe_in.degree, dim - 2);
+ const unsigned int dofs_per_face_normal = fe_in.n_dofs_per_face();
+
+ lexicographic_numbering =
+ get_lexicographic_numbering_rt_nodal<dim>(fe_in.degree);
+
+ // To get the right shape_values of the RT element
+ std::vector<unsigned int> lex_normal, lex_tangent;
+ for (unsigned int i = 0; i < fe.degree; ++i)
+ lex_tangent.push_back(i);
+
+ lex_normal.push_back(0);
+ for (unsigned int i = dofs_per_face_normal * 2 * dim;
+ i < dofs_per_face_normal * 2 * dim + fe.degree - 1;
+ ++i)
+ lex_normal.push_back(i);
+ lex_normal.push_back(dofs_per_face_normal);
+
+ // 'direction' distingusishes between normal and tangential direction
+ for (unsigned int direction = 0; direction < 2; ++direction)
+ {
+ UnivariateShapeData<Number> &univariate_shape_data =
+ (direction == 0) ? data.front() : data.back();
+
+ univariate_shape_data.element_type = tensor_raviart_thomas;
+ univariate_shape_data.quadrature = quad;
+ univariate_shape_data.n_q_points_1d = n_q_points_1d;
+ univariate_shape_data.fe_degree = fe.degree - direction;
+
+ // grant write access to common univariate shape data
+ auto &shape_values = univariate_shape_data.shape_values;
+ auto &shape_gradients = univariate_shape_data.shape_gradients;
+ auto &shape_hessians = univariate_shape_data.shape_hessians;
+
+ auto &values_within_subface =
+ univariate_shape_data.values_within_subface;
+ auto &gradients_within_subface =
+ univariate_shape_data.gradients_within_subface;
+ auto &hessians_within_subface =
+ univariate_shape_data.hessians_within_subface;
+
+ auto &shape_data_on_face =
+ univariate_shape_data.shape_data_on_face;
+
+ const unsigned int n_dofs_1d = fe.degree + 1 - direction;
+ const unsigned int array_size = n_dofs_1d * n_q_points_1d;
+
+ shape_gradients.resize_fast(array_size);
+ shape_values.resize_fast(array_size);
+ shape_hessians.resize_fast(array_size);
+
+ values_within_subface[0].resize(array_size);
+ values_within_subface[1].resize(array_size);
+ gradients_within_subface[0].resize(array_size);
+ gradients_within_subface[1].resize(array_size);
+ hessians_within_subface[0].resize(array_size);
+ hessians_within_subface[1].resize(array_size);
+
+ shape_data_on_face[0].resize(3 * n_dofs_1d);
+ shape_data_on_face[1].resize(3 * n_dofs_1d);
+
+ Point<dim> unit_point;
+ for (unsigned int i = 0; i < n_dofs_1d; ++i)
+ {
+ // need to reorder from hierarchical to lexicographic to get
+ // the DoFs correct
+ const unsigned int my_i =
+ (direction == 0) ? lex_normal[i] : lex_tangent[i];
+ for (unsigned int q = 0; q < n_q_points_1d; ++q)
+ {
+ Point<dim> q_point = unit_point;
+ q_point[direction] = quad.get_points()[q][0];
+
+ shape_values[i * n_q_points_1d + q] =
+ fe.shape_value_component(my_i, q_point, 0);
+ shape_gradients[i * n_q_points_1d + q] =
+ fe.shape_grad_component(my_i, q_point, 0)[direction];
+ shape_hessians[i * n_q_points_1d + q] =
+ fe.shape_grad_grad_component(my_i,
+ q_point,
+ 0)[direction][direction];
+
+ // evaluate basis functions on the two 1D subfaces (i.e.,
+ // at the positions divided by one half and shifted by one
+ // half, respectively) for hanging nodes
+ q_point[direction] *= 0.5;
+ values_within_subface[0][i * n_q_points_1d + q] =
+ fe.shape_value_component(my_i, q_point, 0);
+ gradients_within_subface[0][i * n_q_points_1d + q] =
+ fe.shape_grad_component(my_i, q_point, 0)[direction];
+ hessians_within_subface[0][i * n_q_points_1d + q] =
+ fe.shape_grad_grad_component(my_i,
+ q_point,
+ 0)[direction][direction];
+ q_point[direction] += 0.5;
+ values_within_subface[1][i * n_q_points_1d + q] =
+ fe.shape_value_component(my_i, q_point, 0);
+ gradients_within_subface[1][i * n_q_points_1d + q] =
+ fe.shape_grad_component(my_i, q_point, 0)[direction];
+ hessians_within_subface[1][i * n_q_points_1d + q] =
+ fe.shape_grad_grad_component(my_i,
+ q_point,
+ 0)[direction][direction];
+ }
+ // evaluate basis functions on the 1D faces, i.e., in zero and
+ // one
+ Point<dim> q_point = unit_point;
+ q_point[direction] = 0;
+ shape_data_on_face[0][i] =
+ fe.shape_value_component(my_i, q_point, 0);
+ shape_data_on_face[0][i + n_dofs_1d] =
+ fe.shape_grad_component(my_i, q_point, 0)[direction];
+ shape_data_on_face[0][i + 2 * n_dofs_1d] =
+ fe.shape_grad_grad_component(my_i,
+ q_point,
+ 0)[direction][direction];
+ q_point[direction] = 1;
+ shape_data_on_face[1][i] =
+ fe.shape_value_component(my_i, q_point, 0);
+ shape_data_on_face[1][i + n_dofs_1d] =
+ fe.shape_grad_component(my_i, q_point, 0)[direction];
+ shape_data_on_face[1][i + 2 * n_dofs_1d] =
+ fe.shape_grad_grad_component(my_i,
+ q_point,
+ 0)[direction][direction];
+ }
+ }
+ return;
+ }
+ else if (quad_in.is_tensor_product() == false ||
+ dynamic_cast<const FE_SimplexP<dim> *>(
+ &fe_in.base_element(base_element_number)) ||
+ dynamic_cast<const FE_SimplexDGP<dim> *>(
+ &fe_in.base_element(base_element_number)) ||
+ dynamic_cast<const FE_WedgeP<dim> *>(
+ &fe_in.base_element(base_element_number)) ||
+ dynamic_cast<const FE_PyramidP<dim> *>(
+ &fe_in.base_element(base_element_number)))
{
// specialization for arbitrary finite elements and quadrature rules
// as needed in the context, e.g., of simplices
* coefficient arrays. See the documentation of the EvaluatorTensorProduct
* specialization for more information.
*/
- evaluate_symmetric_hierarchical
+ evaluate_symmetric_hierarchical,
+ /**
+ * Raviart-Thomas elements with anisotropic polynomials.
+ */
+ evaluate_raviart_thomas
};
(void)dummy2;
}
- template <int direction, bool contract_over_rows, bool add>
+ template <int direction,
+ bool contract_over_rows,
+ bool add,
+ int normal_dir = 0>
void
values(const Number in[], Number out[]) const
{
apply<direction, contract_over_rows, add>(shape_values, in, out);
}
- template <int direction, bool contract_over_rows, bool add>
+ template <int direction,
+ bool contract_over_rows,
+ bool add,
+ int normal_dir = 0>
void
gradients(const Number in[], Number out[]) const
{
apply<direction, contract_over_rows, add>(shape_gradients, in, out);
}
- template <int direction, bool contract_over_rows, bool add>
+ template <int direction,
+ bool contract_over_rows,
+ bool add,
+ int normal_dir = 0>
void
hessians(const Number in[], Number out[]) const
{
+ /**
+ * Internal evaluator for shape function in 2D and 3D using the
+ * tensor product form of the anisotropic basis functions of the
+ * raviart-thomas element, with degree k+1 in normal direction and
+ * k in tangential direction.
+ *
+ * @tparam dim Space dimension in which this class is applied
+ * @tparam n_rows Number of rows in the transformation matrix, which corresponds
+ * to the number of 1d shape functions in the usual tensor
+ * contraction setting
+ * @tparam n_columns Number of columns in the transformation matrix, which
+ * corresponds to the number of 1d shape functions in the
+ * usual tensor contraction setting
+ * @tparam Number Abstract number type for input and output arrays
+ * @tparam Number2 Abstract number type for coefficient arrays (defaults to
+ * same type as the input/output arrays); must implement
+ * operator* with Number and produce Number as an output to
+ * be a valid type
+ */
+ template <int dim,
+ int n_rows,
+ int n_columns,
+ typename Number,
+ typename Number2>
+ struct EvaluatorTensorProduct<evaluate_raviart_thomas,
+ dim,
+ n_rows,
+ n_columns,
+ Number,
+ Number2>
+ {
+ static constexpr unsigned int n_rows_of_product =
+ numbers::invalid_unsigned_int;
+ static constexpr unsigned int n_columns_of_product =
+ numbers::invalid_unsigned_int;
+
+ /**
+ * Empty constructor. Does nothing. Be careful when using 'values' and
+ * related methods because they need to be filled with the other pointer
+ */
+ EvaluatorTensorProduct()
+ : shape_values(nullptr)
+ , shape_gradients(nullptr)
+ , shape_hessians(nullptr)
+ {}
+
+ /**
+ * Constructor, taking the data from ShapeInfo
+ */
+ EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
+ const AlignedVector<Number2> &shape_gradients,
+ const AlignedVector<Number2> &shape_hessians,
+ const unsigned int dummy1 = 0,
+ const unsigned int dummy2 = 0)
+ : shape_values(shape_values.begin())
+ , shape_gradients(shape_gradients.begin())
+ , shape_hessians(shape_hessians.begin())
+ {
+ // We can enter this function either for the apply() path that has
+ // n_rows * n_columns entries or for the apply_face() path that only has
+ // n_rows * 3 entries in the array. Since we cannot decide about the use
+ // we must allow for both here.
+ Assert(shape_values.size() == 0 ||
+ shape_values.size() == n_rows * n_columns ||
+ shape_values.size() == 3 * n_rows,
+ ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
+ Assert(shape_gradients.size() == 0 ||
+ shape_gradients.size() == n_rows * n_columns,
+ ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
+ Assert(shape_hessians.size() == 0 ||
+ shape_hessians.size() == n_rows * n_columns,
+ ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
+ (void)dummy1;
+ (void)dummy2;
+ }
+
+ template <int direction, bool contract_over_rows, bool add, int normal_dir>
+ void
+ values(const Number in[], Number out[]) const
+ {
+ apply<direction, contract_over_rows, add, normal_dir>(shape_values,
+ in,
+ out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add, int normal_dir>
+ void
+ gradients(const Number in[], Number out[]) const
+ {
+ apply<direction, contract_over_rows, add, normal_dir>(shape_gradients,
+ in,
+ out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add, int normal_dir>
+ void
+ hessians(const Number in[], Number out[]) const
+ {
+ apply<direction, contract_over_rows, add, normal_dir>(shape_hessians,
+ in,
+ out);
+ }
+
+ /**
+ * This function applies the tensor product kernel, corresponding to a
+ * multiplication of 1D stripes, along the given @p direction of the tensor
+ * data in the input array. This function allows the @p in and @p out
+ * arrays to alias for the case n_rows == n_columns, i.e., it is safe to
+ * perform the contraction in place where @p in and @p out point to the
+ * same address. For the case n_rows != n_columns, the output is only
+ * correct if @p one_line is set to true.
+ *
+ * @tparam direction Direction that is evaluated
+ * @tparam contract_over_rows If true, the tensor contraction sums
+ * over the rows in the given @p shape_data
+ * array, otherwise it sums over the columns
+ * @tparam add If true, the result is added to the output vector, else
+ * the computed values overwrite the content in the output
+ * @tparam normal_dir Indicates the direction of the continuous component of the
+ * RT space in terms of the normal onto the face, e.g
+ * 0 if the is in x-direction, 1 if in y-direction
+ * etc.
+ * @tparam one_line If true, the kernel is only applied along a single 1D
+ * stripe within a dim-dimensional tensor, not the full
+ * n_rows^dim points as in the @p false case.
+ *
+ * @param shape_data Transformation matrix with @p n_rows rows and
+ * @p n_columns columns, stored in row-major format
+ * @param in Pointer to the start of the input data vector
+ * @param out Pointer to the start of the output data vector
+ */
+ template <int direction,
+ bool contract_over_rows,
+ bool add,
+ int normal_dir,
+ bool one_line = false>
+ static void
+ apply(const Number2 *DEAL_II_RESTRICT shape_data,
+ const Number * in,
+ Number * out);
+
+ template <int face_direction,
+ bool contract_onto_face,
+ bool add,
+ int max_derivative,
+ bool lex_faces = false,
+ int normal_direction>
+ void
+ apply_face(const Number *DEAL_II_RESTRICT in,
+ Number *DEAL_II_RESTRICT out) const;
+
+ const Number2 *shape_values;
+ const Number2 *shape_gradients;
+ const Number2 *shape_hessians;
+ };
+
+ template <int dim,
+ int n_rows,
+ int n_columns,
+ typename Number,
+ typename Number2>
+ template <int direction,
+ bool contract_over_rows,
+ bool add,
+ int normal_dir,
+ bool one_line>
+ inline void
+ EvaluatorTensorProduct<evaluate_raviart_thomas,
+ dim,
+ n_rows,
+ n_columns,
+ Number,
+ Number2>::apply(const Number2 *DEAL_II_RESTRICT
+ shape_data,
+ const Number *in,
+ Number * out)
+ {
+ static_assert(one_line == false || direction == dim - 1,
+ "Single-line evaluation only works for direction=dim-1.");
+ Assert(shape_data != nullptr,
+ ExcMessage(
+ "The given array shape_data must not be the null pointer!"));
+ Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
+ in != out,
+ ExcMessage("In-place operation only supported for "
+ "n_rows==n_columns or single-line interpolation"));
+ AssertIndexRange(direction, dim);
+ constexpr int mm = contract_over_rows ? n_rows : n_columns,
+ nn = contract_over_rows ? n_columns : n_rows;
+
+ constexpr int stride = Utilities::pow(n_columns, direction);
+ constexpr int n_blocks1 = one_line ? 1 : stride;
+
+ // The number of blocks depend on both direction and dimension.
+ constexpr int n_blocks2 =
+ (dim - direction - 1 == 0) ?
+ 1 :
+ ((direction == normal_dir) ?
+ Utilities::pow((n_rows - 1),
+ (direction >= dim) ? 0 : dim - direction - 1) :
+ (((direction < normal_dir) ? (n_rows + 1) : n_rows) *
+ ((dim - direction == 3) ? n_rows : 1)));
+
+ for (int i2 = 0; i2 < n_blocks2; ++i2)
+ {
+ for (int i1 = 0; i1 < n_blocks1; ++i1)
+ {
+ Number x[mm];
+ for (int i = 0; i < mm; ++i)
+ x[i] = in[stride * i];
+
+ for (int col = 0; col < nn; ++col)
+ {
+ Number2 val0;
+
+ if (contract_over_rows)
+ val0 = shape_data[col];
+ else
+ val0 = shape_data[col * n_columns];
+
+ Number res0 = val0 * x[0];
+ for (int i = 1; i < mm; ++i)
+ {
+ if (contract_over_rows)
+ val0 = shape_data[i * n_columns + col];
+ else
+ val0 = shape_data[col * n_columns + i];
+
+ res0 += val0 * x[i];
+ }
+ if (add)
+ out[stride * col] += res0;
+
+ else
+ out[stride * col] = res0;
+ }
+
+ if (one_line == false)
+ {
+ ++in;
+ ++out;
+ }
+ }
+ if (one_line == false)
+ {
+ in += stride * (mm - 1);
+ out += stride * (nn - 1);
+ }
+ }
+ }
+
+ template <int dim,
+ int n_rows,
+ int n_columns,
+ typename Number,
+ typename Number2>
+ template <int face_direction,
+ bool contract_onto_face,
+ bool add,
+ int max_derivative,
+ bool lex_faces,
+ int normal_direction>
+ inline void
+ EvaluatorTensorProduct<evaluate_raviart_thomas,
+ dim,
+ n_rows,
+ n_columns,
+ Number,
+ Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
+ Number *DEAL_II_RESTRICT
+ out) const
+ {
+ Assert(dim > 1 && (lex_faces || dim < 4),
+ ExcMessage("Only dim=2,3 supported"));
+ static_assert(max_derivative >= 0 && max_derivative < 3,
+ "Only derivative orders 0-2 implemented");
+ static_assert(!lex_faces,
+ "lex_faces = True is not implemented for Raviart-Thomas");
+ Assert(shape_values != nullptr,
+ ExcMessage(
+ "The given array shape_values must not be the null pointer."));
+
+ // Determine the number of blocks depending on the face and normaldirection,
+ // as well as dimension.
+ constexpr int n_blocks1 = (face_direction == normal_direction) ?
+ (n_rows - 1) :
+ ((face_direction == 0 && normal_direction == 2) ||
+ (face_direction == 1 && normal_direction == 2) ||
+ (face_direction == 2 && normal_direction == 1)) ?
+ n_rows :
+ (n_rows + 1);
+ constexpr int n_blocks2 =
+ (dim == 2) ? 1 :
+ ((face_direction == normal_direction) ?
+ (n_rows - 1) :
+ (((face_direction == 0 && normal_direction == 1) ||
+ (face_direction == 1 && normal_direction == 0) ||
+ (face_direction == 2 && normal_direction == 0)) ?
+ n_rows :
+ (n_rows + 1)));
+
+ AssertIndexRange(face_direction, dim);
+
+ constexpr int stride =
+ (face_direction == normal_direction) ?
+ Utilities::pow(n_rows - 1, face_direction) :
+ ((face_direction == 0) ?
+ 1 :
+ ((face_direction == 2) ?
+ n_rows * (n_rows + 1) :
+ ((face_direction == 1 && normal_direction == 0) ? (n_rows + 1) :
+ n_rows)));
+ constexpr int out_stride = n_blocks1 * n_blocks2;
+
+ const Number *DEAL_II_RESTRICT shape_values = this->shape_values;
+
+ for (int i2 = 0; i2 < n_blocks2; ++i2)
+ {
+ for (int i1 = 0; i1 < n_blocks1; ++i1)
+ {
+ if (contract_onto_face == true)
+ {
+ Number res0 = shape_values[0] * in[0];
+ Number res1, res2;
+
+ if (max_derivative > 0)
+ res1 = shape_values[n_rows] * in[0];
+
+ if (max_derivative > 1)
+ res2 = shape_values[2 * n_rows] * in[0];
+
+ for (int ind = 1; ind < n_rows; ++ind)
+ {
+ res0 += shape_values[ind] * in[stride * ind];
+ if (max_derivative > 0)
+ res1 += shape_values[ind + n_rows] * in[stride * ind];
+
+ if (max_derivative > 1)
+ res2 += shape_values[ind + 2 * n_rows] * in[stride * ind];
+ }
+ if (add)
+ {
+ out[0] += res0;
+
+ if (max_derivative > 0)
+ out[out_stride] += res1;
+
+ if (max_derivative > 1)
+ out[2 * out_stride] += res2;
+ }
+ else
+ {
+ out[0] = res0;
+
+ if (max_derivative > 0)
+ out[out_stride] = res1;
+
+ if (max_derivative > 1)
+ out[2 * out_stride] = res2;
+ }
+ }
+ else
+ {
+ for (int col = 0; col < n_rows; ++col)
+ {
+ if (add)
+ out[col * stride] += shape_values[col] * in[0];
+ else
+ out[col * stride] = shape_values[col] * in[0];
+
+ if (max_derivative > 0)
+ out[col * stride] +=
+ shape_values[col + n_rows] * in[out_stride];
+
+ if (max_derivative > 1)
+ out[col * stride] +=
+ shape_values[col + 2 * n_rows] * in[2 * out_stride];
+ }
+ }
+
+ // 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 (face_direction)
+ {
+ case 0:
+ in += contract_onto_face ? n_rows : 1;
+ out += contract_onto_face ? 1 : n_rows;
+ break;
+
+ case 1:
+ ++in;
+ ++out;
+ // faces 2 and 3 in 3D use local coordinate system zx, which
+ // is the other way around compared to the tensor
+ // product. Need to take that into account.
+ if (dim == 3)
+ {
+ if (normal_direction == 0)
+ {
+ if (contract_onto_face)
+ out += n_rows - 1;
+ else
+ in += n_rows - 1;
+ }
+ if (normal_direction == 1)
+ {
+ if (contract_onto_face)
+ out += n_rows - 2;
+ else
+ in += n_rows - 2;
+ }
+ if (normal_direction == 2)
+ {
+ if (contract_onto_face)
+ out += n_rows;
+ else
+ in += n_rows;
+ }
+ }
+ break;
+
+ case 2:
+ ++in;
+ ++out;
+ break;
+
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ if (face_direction == 1 && dim == 3)
+ {
+ // adjust for local coordinate system zx
+ if (contract_onto_face)
+ {
+ if (normal_direction == 0)
+ {
+ in += (n_rows + 1) * (n_rows - 1);
+ out -= n_rows * (n_rows + 1) - 1;
+ }
+ if (normal_direction == 1)
+ {
+ in += (n_rows - 1) * (n_rows - 1);
+ out -= (n_rows - 1) * (n_rows - 1) - 1;
+ }
+ if (normal_direction == 2)
+ {
+ in += (n_rows - 1) * (n_rows);
+ out -= (n_rows) * (n_rows + 1) - 1;
+ }
+ }
+ else
+ {
+ if (normal_direction == 0)
+ {
+ out += (n_rows + 1) * (n_rows - 1);
+ in -= n_rows * (n_rows + 1) - 1;
+ }
+ if (normal_direction == 1)
+ {
+ out += (n_rows - 1) * (n_rows - 1);
+ in -= (n_rows - 1) * (n_rows - 1) - 1;
+ }
+ if (normal_direction == 2)
+ {
+ out += (n_rows - 1) * (n_rows);
+ in -= (n_rows) * (n_rows + 1) - 1;
+ }
+ }
+ }
+ }
+ }
+
+
+
/**
* Struct to avoid using Tensor<1, dim, Point<dim2>> in
* evaluate_tensor_product_value_and_gradient because a Point cannot be used
}
-
- // set up the numbering of the rt polynomials
- std::vector<unsigned int>
- compute_rt_hierarchic_to_lexicographic(const unsigned int dim,
- const unsigned int degree)
- {
- const unsigned int n_pols =
- (degree + 2) * Utilities::pow(degree + 1, dim - 1);
-
- std::vector<unsigned int> hierarchic_to_lexicographic;
-
- // dofs on faces
- for (unsigned int face_no = 0; face_no < 2 * dim; ++face_no)
- {
- const unsigned int stride_x = face_no < 2 ? degree + 2 : 1;
- const unsigned int stride_y =
- face_no < 4 ? (degree + 2) * (degree + 1) : degree + 1;
- const unsigned int offset =
- (face_no % 2) * Utilities::pow(degree + 1, 1 + face_no / 2);
- for (unsigned int j = 0; j < (dim > 2 ? degree + 1 : 1); ++j)
- for (unsigned int i = 0; i < degree + 1; ++i)
- hierarchic_to_lexicographic.push_back(
- (face_no / 2) * n_pols + offset + i * stride_x + j * stride_y);
- }
- // dofs on cells, starting with x component...
- for (unsigned int k = 0; k < (dim > 2 ? degree + 1 : 1); ++k)
- for (unsigned int j = 0; j < (dim > 1 ? degree + 1 : 1); ++j)
- for (unsigned int i = 1; i < degree + 1; ++i)
- hierarchic_to_lexicographic.push_back(
- k * (degree + 1) * (degree + 2) + j * (degree + 2) + i);
- // ... then y component ...
- if (dim > 1)
- for (unsigned int k = 0; k < (dim > 2 ? degree + 1 : 1); ++k)
- for (unsigned int j = 1; j < degree + 1; ++j)
- for (unsigned int i = 0; i < degree + 1; ++i)
- hierarchic_to_lexicographic.push_back(
- n_pols + k * (degree + 1) * (degree + 2) + j * (degree + 1) + i);
- // ... and finally z component
- if (dim > 2)
- for (unsigned int k = 1; k < degree + 1; ++k)
- for (unsigned int j = 0; j < degree + 1; ++j)
- for (unsigned int i = 0; i < degree + 1; ++i)
- hierarchic_to_lexicographic.push_back(
- 2 * n_pols + k * (degree + 1) * (degree + 1) + j * (degree + 1) +
- i);
-
- AssertDimension(hierarchic_to_lexicographic.size(), n_pols * dim);
-
-#ifdef DEBUG
- // assert that we have a valid permutation
- std::vector<unsigned int> copy(hierarchic_to_lexicographic);
- std::sort(copy.begin(), copy.end());
- for (unsigned int i = 0; i < copy.size(); ++i)
- AssertDimension(i, copy[i]);
-#endif
-
- return hierarchic_to_lexicographic;
- }
-
-
-
template <int dim>
PolynomialsRaviartThomasNodal<dim>::PolynomialsRaviartThomasNodal(
const unsigned int degree)
// actual order required by the finite element class with unknowns on
// faces placed first
const unsigned int n_pols = polynomial_space.n();
- hierarchic_to_lexicographic =
- compute_rt_hierarchic_to_lexicographic(dim, degree);
-
lexicographic_to_hierarchic =
- Utilities::invert_permutation(hierarchic_to_lexicographic);
+ internal::get_lexicographic_numbering_rt_nodal<dim>(degree + 1);
+
+ hierarchic_to_lexicographic =
+ Utilities::invert_permutation(lexicographic_to_hierarchic);
// since we only store an anisotropic polynomial for the first component,
// we set up a second numbering to switch out the actual coordinate
}
} // namespace
+namespace internal
+{
+ template <int dim>
+ std::vector<unsigned int>
+ get_lexicographic_numbering_rt_nodal(const unsigned int degree)
+ {
+ const unsigned int n_dofs_face = Utilities::pow(degree, dim - 1);
+ std::vector<unsigned int> lexicographic_numbering;
+ // component 1
+ for (unsigned int j = 0; j < n_dofs_face; j++)
+ {
+ lexicographic_numbering.push_back(j);
+ for (unsigned int i = n_dofs_face * 2 * dim;
+ i < n_dofs_face * 2 * dim + degree - 1;
+ i++)
+ lexicographic_numbering.push_back(i + j * (degree - 1));
+ lexicographic_numbering.push_back(n_dofs_face + j);
+ }
+
+ // component 2
+ unsigned int layers = (dim == 3) ? degree : 1;
+ for (unsigned int k = 0; k < layers; k++)
+ {
+ unsigned int k_add = k * degree;
+ for (unsigned int j = n_dofs_face * 2; j < n_dofs_face * 2 + degree;
+ j++)
+ lexicographic_numbering.push_back(j + k_add);
+
+ for (unsigned int i = n_dofs_face * (2 * dim + (degree - 1));
+ i < n_dofs_face * (2 * dim + (degree - 1)) + (degree - 1) * degree;
+ i++)
+ {
+ lexicographic_numbering.push_back(i + k_add * (degree - 1));
+ }
+ for (unsigned int j = n_dofs_face * 3; j < n_dofs_face * 3 + degree;
+ j++)
+ lexicographic_numbering.push_back(j + k_add);
+ }
+ // component 3
+ if (dim == 3)
+ {
+ for (unsigned int i = 4 * n_dofs_face; i < 5 * n_dofs_face; i++)
+ lexicographic_numbering.push_back(i);
+ for (unsigned int i = 6 * n_dofs_face + n_dofs_face * 2 * (degree - 1);
+ i < 6 * n_dofs_face + n_dofs_face * 3 * (degree - 1);
+ i++)
+ lexicographic_numbering.push_back(i);
+ for (unsigned int i = 5 * n_dofs_face; i < 6 * n_dofs_face; i++)
+ lexicographic_numbering.push_back(i);
+ }
+
+ return lexicographic_numbering;
+ }
+} // namespace internal
// --------------------- actual implementation of element --------------------
evaluation_template_factory_inst4.cc
evaluation_template_factory_inst5.cc
evaluation_template_factory_inst6.cc
+ evaluation_template_factory_inst7.cc
+ evaluation_template_factory_inst8.cc
+ evaluation_template_factory_inst9.cc
+ evaluation_template_factory_inst10.cc
evaluation_template_factory_hanging_nodes.cc
+ evaluation_template_face_factory.cc
+ evaluation_template_face_factory_inst2.cc
+ evaluation_template_face_factory_inst3.cc
+ evaluation_template_face_factory_inst4.cc
+ evaluation_template_face_factory_inst5.cc
+ evaluation_template_face_factory_inst6.cc
fe_point_evaluation.cc
mapping_info.cc
mapping_info_inst2.cc
SET(_inst
evaluation_template_factory.inst.in
+ evaluation_template_face_factory.inst.in
evaluation_template_factory_hanging_nodes.inst.in
fe_point_evaluation.inst.in
mapping_info.inst.in
--- /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.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/matrix_free/evaluation_template_face_factory.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+#define SPLIT_INSTANTIATIONS_COUNT 6
+#ifndef SPLIT_INSTANTIATIONS_INDEX
+# define SPLIT_INSTANTIATIONS_INDEX 0
+#endif
+
+#include "evaluation_template_face_factory.inst"
+
+DEAL_II_NAMESPACE_CLOSE
--- /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.
+//
+// ---------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS;
+ deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED)
+ {
+ template struct dealii::internal::
+ FEFaceEvaluationFactory<deal_II_dimension, deal_II_scalar_vectorized>;
+
+ template struct dealii::internal::FEFaceEvaluationGatherFactory<
+ deal_II_dimension,
+ double,
+ deal_II_scalar_vectorized>;
+
+ template struct dealii::internal::FEFaceEvaluationGatherFactory<
+ deal_II_dimension,
+ float,
+ deal_II_scalar_vectorized>;
+ }
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 1
+#include "evaluation_template_face_factory.cc"
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 2
+#include "evaluation_template_face_factory.cc"
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 3
+#include "evaluation_template_face_factory.cc"
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 4
+#include "evaluation_template_face_factory.cc"
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 5
+#include "evaluation_template_face_factory.cc"
DEAL_II_NAMESPACE_OPEN
-#define SPLIT_INSTANTIATIONS_COUNT 6
+#define SPLIT_INSTANTIATIONS_COUNT 10
#ifndef SPLIT_INSTANTIATIONS_INDEX
# define SPLIT_INSTANTIATIONS_INDEX 0
#endif
template struct dealii::internal::
FEEvaluationFactory<deal_II_dimension, deal_II_scalar_vectorized>;
- template struct dealii::internal::
- FEFaceEvaluationFactory<deal_II_dimension, deal_II_scalar_vectorized>;
-
- template struct dealii::internal::FEFaceEvaluationGatherFactory<
- deal_II_dimension,
- double,
- deal_II_scalar_vectorized>;
-
- template struct dealii::internal::FEFaceEvaluationGatherFactory<
- deal_II_dimension,
- float,
- deal_II_scalar_vectorized>;
-
// inverse mass
template struct dealii::internal::
CellwiseInverseMassFactory<deal_II_dimension, deal_II_scalar_vectorized>;
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 9
+#include "evaluation_template_factory.cc"
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 6
+#include "evaluation_template_factory.cc"
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 7
+#include "evaluation_template_factory.cc"
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 8
+#include "evaluation_template_factory.cc"
--- /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 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
+
+#include "../tests.h"
+
+#include "matrix_vector_rt_common.h"
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(2);
+
+ FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+
+ DoFHandler<dim> dof(tria);
+
+ dof.distribute_dofs(fe);
+
+ AffineConstraints<double> constraints;
+ constraints.close();
+ do_test<dim, fe_degree, double>(dof, constraints);
+}
--- /dev/null
+
+DEAL:2d::Testing 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::
+DEAL:2d::Testing 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::
+DEAL:3d::Testing 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::
+DEAL:3d::Testing 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::
--- /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 template evaluates a simple operator on FE_PolyTensor
+// elements using FEEvaluation and compares the result with the output
+// of FEValues (which is considered to be the reference) on cartesian
+// meshes without hanging nodes. (It will be extended to also handle general
+// meshes and hanging nodes in the future.) The test do not include
+// multithreading because FEValues is not thread-safe.
+// See matrix_vector_rt_01.cc for an example that uses this template.
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q1_eulerian.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/matrix_free.templates.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+// forward declare this function. will be implemented in .cc files
+template <int dim, int fe_degree>
+void
+test();
+
+
+
+template <int dim,
+ int fe_degree,
+ int n_q_points_1d = fe_degree + 1,
+ typename Number = double>
+class MatrixFreeTest
+{
+public:
+ MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
+ : data(data_in){};
+
+ virtual ~MatrixFreeTest(){};
+
+ void
+ operator()(const MatrixFree<dim, Number> & data,
+ Vector<Number> & dst,
+ const Vector<Number> & src,
+ const std::pair<unsigned int, unsigned int> &cell_range) const
+ {
+ 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);
+
+ 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);
+ }
+
+ fe_eval.integrate_scatter(EvaluationFlags::values |
+ EvaluationFlags::gradients,
+ dst);
+ }
+ };
+
+
+ void
+ test_functions(Vector<Number> &dst, const Vector<Number> &src) const
+ {
+ data.cell_loop(&MatrixFreeTest::operator(), this, dst, src);
+ };
+
+protected:
+ const MatrixFree<dim, Number> &data;
+};
+
+
+template <int dim, int fe_degree, typename Number>
+void
+do_test(const DoFHandler<dim> & dof,
+ const AffineConstraints<double> &constraints)
+{
+ 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;
+
+
+ // constraints.distribute(solution);
+ MatrixFree<dim, Number> mf_data;
+ {
+ const QGaussLobatto<1> quad(fe_degree + 2);
+ typename MatrixFree<dim, Number>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim, Number>::AdditionalData::none;
+ data.mapping_update_flags = update_gradients | update_hessians;
+ mf_data.reinit(dof, constraints, quad, data);
+ }
+
+ // create vector with random entries
+ Vector<Number> solution, initial_condition;
+
+ mf_data.initialize_dof_vector(solution);
+ mf_data.initialize_dof_vector(initial_condition);
+
+ for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+ {
+ if (constraints.is_constrained(i))
+ continue;
+ initial_condition[i] = random_value<Number>();
+ }
+
+ // MatrixFree solution
+ MatrixFreeTest<dim, fe_degree, fe_degree + 2, Number> mf(mf_data);
+ mf.test_functions(solution, initial_condition);
+
+
+ SparsityPattern sp;
+ SparseMatrix<double> system_matrix;
+ DynamicSparsityPattern dsp(dof.n_dofs(), dof.n_dofs());
+ DoFTools::make_sparsity_pattern(dof, dsp);
+ sp.copy_from(dsp);
+ system_matrix.reinit(sp);
+
+ FEValues<dim> fe_val(dof.get_fe(),
+ Quadrature<dim>(mf_data.get_quadrature(0)),
+ update_values | update_gradients | update_JxW_values);
+
+
+ const unsigned int dofs_per_cell = fe_val.get_fe().dofs_per_cell;
+ std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+ FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
+
+
+ const FEValuesExtractors::Vector velocities(0);
+
+ // Assemble matrix
+ for (const auto &cell : dof.active_cell_iterators())
+ {
+ fe_val.reinit(cell);
+ local_matrix = 0;
+
+ for (const auto q : fe_val.quadrature_point_indices())
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ 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);
+
+ 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);
+ }
+ }
+ cell->get_dof_indices(local_dof_indices);
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ system_matrix.add(local_dof_indices[i],
+ local_dof_indices[j],
+ local_matrix(i, j));
+ }
+
+ Vector<Number> ref(solution.size());
+
+ // Compute reference
+ system_matrix.vmult(ref, initial_condition);
+
+ ref -= solution;
+
+ const double diff_norm = ref.linfty_norm() / solution.linfty_norm();
+ deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+int
+main()
+{
+ initlog();
+
+ {
+ deallog.push("2d");
+ test<2, 2>();
+ test<2, 3>();
+ deallog.pop();
+ deallog.push("3d");
+ test<3, 2>();
+ test<3, 3>();
+ deallog.pop();
+ }
+}
--- /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 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
+
+#include "../tests.h"
+
+#include "matrix_vector_rt_face_common.h"
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(1);
+
+ FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+
+ DoFHandler<dim> dof(tria);
+
+ dof.distribute_dofs(fe);
+
+ AffineConstraints<double> constraints;
+ constraints.close();
+ do_test<dim, fe_degree, double>(dof, constraints);
+}
--- /dev/null
+
+DEAL:2d::Testing 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::
+DEAL:2d::Testing 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::
+DEAL:3d::Testing 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::
--- /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 template evaluates a simple operator on FE_PolyTensor
+// elements using FEFaceEvaluation and compares the result with the output
+// of FEFaceValues (which is considered to be the reference) on cartesian
+// meshes without hanging nodes. (It will be extended to also handle general
+// meshes and hanging nodes in the future.) The test do not include
+// multithreading because FEFaceValues is not thread-safe.
+// See matrix_vector_rt_face_01.cc for an example that uses this template.
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q1_eulerian.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/matrix_free.templates.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+// forward declare this function. will be implemented in .cc files
+template <int dim, int fe_degree>
+void
+test();
+
+
+
+template <int dim,
+ int fe_degree,
+ int n_q_points_1d = fe_degree + 1,
+ typename Number = double>
+class MatrixFreeTest
+{
+public:
+ MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
+ : data(data_in){};
+
+ virtual ~MatrixFreeTest(){};
+
+ void
+ dummy(const MatrixFree<dim, Number> &,
+ Vector<Number> &,
+ const Vector<Number> &,
+ const std::pair<unsigned int, unsigned int> &) const
+ {}
+
+ void
+ operator_face(const MatrixFree<dim, Number> & data,
+ Vector<Number> & dst,
+ const Vector<Number> & src,
+ const std::pair<unsigned int, unsigned int> &face_range) const
+ {
+ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_eval(data,
+ true);
+ 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);
+
+ 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);
+ }
+
+ fe_eval.integrate_scatter(EvaluationFlags::values |
+ EvaluationFlags::gradients,
+ dst);
+ fe_eval_n.integrate_scatter(EvaluationFlags::values |
+ EvaluationFlags::gradients,
+ dst);
+ }
+ };
+
+ void
+ operator_boundary(
+ const MatrixFree<dim, Number> & data,
+ Vector<Number> & dst,
+ const Vector<Number> & src,
+ const std::pair<unsigned int, unsigned int> &face_range) const
+ {
+ 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);
+
+ 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);
+ }
+
+ fe_eval.integrate_scatter(EvaluationFlags::values |
+ EvaluationFlags::gradients,
+ dst);
+ }
+ };
+
+
+ void
+ test_functions(Vector<Number> &dst, const Vector<Number> &src) const
+ {
+ data.loop(&MatrixFreeTest::dummy,
+ &MatrixFreeTest::operator_face,
+ &MatrixFreeTest::operator_boundary,
+ this,
+ dst,
+ src);
+ };
+
+protected:
+ const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename Number>
+void
+do_test(const DoFHandler<dim> & dof,
+ const AffineConstraints<double> &constraints)
+{
+ 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;
+
+
+ // constraints.distribute(solution);
+ MatrixFree<dim, Number> mf_data;
+ {
+ const QGaussLobatto<1> quad(fe_degree + 2);
+ const MappingQ<dim> mapping(fe_degree);
+ typename MatrixFree<dim, Number>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim, Number>::AdditionalData::none;
+ data.mapping_update_flags = update_gradients | update_JxW_values;
+ data.mapping_update_flags_inner_faces =
+ (update_gradients | update_JxW_values);
+ data.mapping_update_flags_boundary_faces =
+ (update_gradients | update_JxW_values);
+ mf_data.reinit(mapping, dof, constraints, quad, data);
+ }
+
+ Vector<Number> solution, initial_condition;
+
+ // create vector with random entries
+ mf_data.initialize_dof_vector(solution);
+ mf_data.initialize_dof_vector(initial_condition);
+
+ for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+ {
+ if (constraints.is_constrained(i))
+ continue;
+ initial_condition[i] = random_value<Number>();
+ }
+
+ MatrixFreeTest<dim, fe_degree, fe_degree + 2, Number> mf(mf_data);
+ mf.test_functions(solution, initial_condition);
+
+
+ SparsityPattern sp;
+ SparseMatrix<double> system_matrix;
+ DynamicSparsityPattern dsp(dof.n_dofs(), dof.n_dofs());
+ DoFTools::make_sparsity_pattern(dof, dsp);
+ sp.copy_from(dsp);
+ system_matrix.reinit(sp);
+
+ FEFaceValues<dim> fe_val(
+ dof.get_fe(),
+ QGaussLobatto<dim - 1>(mf_data.get_quadrature(0).size()),
+ update_values | update_gradients | update_JxW_values);
+
+
+ const unsigned int dofs_per_cell = fe_val.get_fe().dofs_per_cell;
+ std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+ FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
+
+
+ const FEValuesExtractors::Vector velocities(0);
+ // Assemble matrix
+ for (const auto &cell : dof.active_cell_iterators())
+ for (const auto &face : cell->face_iterators())
+ {
+ fe_val.reinit(cell, face);
+ local_matrix = 0;
+
+ for (const auto q : fe_val.quadrature_point_indices())
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ 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);
+
+ 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);
+ }
+ }
+ cell->get_dof_indices(local_dof_indices);
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ system_matrix.add(local_dof_indices[i],
+ local_dof_indices[j],
+ local_matrix(i, j));
+ }
+
+ Vector<Number> ref(solution.size());
+
+ // Compute reference
+ system_matrix.vmult(ref, initial_condition);
+
+ ref -= solution;
+
+ const double diff_norm = ref.linfty_norm() / solution.linfty_norm();
+ deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+int
+main()
+{
+ initlog();
+
+ {
+ deallog.push("2d");
+ test<2, 2>();
+ test<2, 3>();
+ deallog.pop();
+ deallog.push("3d");
+ test<3, 2>();
+ // test<3, 3>(); //Takes way too long
+ deallog.pop();
+ }
+}
projects / dealii.git / commitdiff