16 #ifndef dealii_fe_point_evaluation_h
17 #define dealii_fe_point_evaluation_h
47 template <
int dim,
int n_components,
typename Number>
55 const unsigned int component,
59 result[component] = vector_entry;
64 const unsigned int component,
68 vector_entry = result[component];
75 const unsigned int vector_lane,
78 for (
unsigned int i = 0; i < n_components; ++i)
79 for (
unsigned int d = 0;
d < dim; ++
d)
80 result[i][
d] =
value[
d][i][vector_lane];
86 const unsigned int vector_lane,
89 for (
unsigned int i = 0; i < n_components; ++i)
90 for (
unsigned int d = 0;
d < dim; ++
d)
91 value[
d][i][vector_lane] = result[i][
d];
96 const unsigned int vector_lane,
99 for (
unsigned int i = 0; i < n_components; ++i)
100 result[i] =
value[i][vector_lane];
105 const unsigned int vector_lane,
108 for (
unsigned int i = 0; i < n_components; ++i)
109 value[i][vector_lane] = result[i];
112 template <
typename Number2>
115 const unsigned int component)
117 return value[component];
120 template <
typename Number2>
121 static const Number2 &
123 const unsigned int component)
125 return value[component];
129 template <
int dim,
typename Number>
140 result = vector_entry;
148 vector_entry = result;
153 const unsigned int vector_lane,
156 for (
unsigned int d = 0;
d < dim; ++
d)
157 result[
d] =
value[
d][vector_lane];
162 const unsigned int vector_lane,
165 for (
unsigned int d = 0;
d < dim; ++
d)
166 value[
d][vector_lane] = result[
d];
171 const unsigned int vector_lane,
174 result =
value[vector_lane];
179 const unsigned int vector_lane,
182 value[vector_lane] = result;
185 template <
typename Number2>
187 access(Number2 &value,
const unsigned int)
192 template <
typename Number2>
193 static const Number2 &
194 access(
const Number2 &value,
const unsigned int)
200 template <
int dim,
typename Number>
208 const unsigned int component,
211 result[component] = vector_entry;
216 const unsigned int component,
219 vector_entry = result[component];
225 const unsigned int vector_lane,
228 for (
unsigned int i = 0; i < dim; ++i)
229 for (
unsigned int d = 0;
d < dim; ++
d)
230 result[i][
d] =
value[
d][i][vector_lane];
236 const unsigned int vector_lane,
239 for (
unsigned int i = 0; i < dim; ++i)
240 for (
unsigned int d = 0;
d < dim; ++
d)
241 value[
d][i][vector_lane] = result[i][
d];
246 const unsigned int vector_lane,
249 for (
unsigned int i = 0; i < dim; ++i)
250 result[i] =
value[i][vector_lane];
255 const unsigned int vector_lane,
258 for (
unsigned int i = 0; i < dim; ++i)
259 value[i][vector_lane] = result[i];
265 return value[component];
268 static const Number &
271 return value[component];
277 return value[component];
283 return value[component];
287 template <
typename Number>
298 result = vector_entry;
306 vector_entry = result;
311 const unsigned int vector_lane,
314 result[0] =
value[0][vector_lane];
319 const unsigned int vector_lane,
322 value[0][vector_lane] = result[0];
327 const unsigned int vector_lane,
330 result =
value[vector_lane];
335 const unsigned int vector_lane,
338 value[vector_lane] = result;
341 template <
typename Number2>
343 access(Number2 &value,
const unsigned int)
348 template <
typename Number2>
349 static const Number2 &
350 access(
const Number2 &value,
const unsigned int)
356 template <
int dim,
int spacedim>
359 const unsigned int base_element_number);
361 template <
int dim,
int spacedim>
365 template <
int dim,
int spacedim>
366 std::vector<Polynomials::Polynomial<double>>
403 template <
int n_components,
406 typename Number =
double>
436 const unsigned int first_selected_component = 0);
456 const unsigned int first_selected_component = 0);
523 get_value(
const unsigned int point_index)
const;
577 jacobian(
const unsigned int point_index)
const;
593 real_point(
const unsigned int point_index)
const;
600 unit_point(
const unsigned int point_index)
const;
612 setup(
const unsigned int first_selected_component);
628 std::vector<Polynomials::Polynomial<double>>
poly;
707 std::unique_ptr<NonMatching::MappingInfo<dim, spacedim>>
730 template <
int n_components,
int dim,
int spacedim,
typename Number>
735 const unsigned int first_selected_component)
738 , update_flags(update_flags)
739 , mapping_info_on_the_fly(
740 std::make_unique<
NonMatching::MappingInfo<dim, spacedim>>(mapping,
742 , mapping_info(mapping_info_on_the_fly.get())
744 setup(first_selected_component);
749 template <
int n_components,
int dim,
int spacedim,
typename Number>
753 const unsigned int first_selected_component)
754 : mapping(&mapping_info.get_mapping())
756 , update_flags(mapping_info.get_update_flags())
757 , mapping_info(&mapping_info)
759 setup(first_selected_component);
764 template <
int n_components,
int dim,
int spacedim,
typename Number>
767 const unsigned int first_selected_component)
770 fe->n_components() + 1);
772 bool same_base_element =
true;
773 unsigned int base_element_number = 0;
774 component_in_base_element = 0;
775 unsigned int component = 0;
776 for (; base_element_number < fe->n_base_elements(); ++base_element_number)
777 if (component + fe->element_multiplicity(base_element_number) >
778 first_selected_component)
780 if (first_selected_component + n_components >
781 component + fe->element_multiplicity(base_element_number))
782 same_base_element =
false;
783 component_in_base_element = first_selected_component - component;
787 component += fe->element_multiplicity(base_element_number);
791 *fe, base_element_number) &&
800 fe->base_element(base_element_number));
802 polynomials_are_hat_functions =
803 (poly.size() == 2 && poly[0].value(0.) == 1. &&
804 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
805 poly[1].value(1.) == 1.);
811 nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
812 for (
unsigned int d = 0;
d < n_components; ++
d)
814 const unsigned int component = first_selected_component +
d;
815 for (
unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
817 const bool is_primitive =
818 fe->is_primitive() || fe->is_primitive(i);
820 nonzero_shape_function_component[i][
d] =
821 (component == fe->system_to_component_index(i).first);
823 nonzero_shape_function_component[i][
d] =
824 (fe->get_nonzero_components(i)[component] ==
true);
834 template <
int n_components,
int dim,
int spacedim,
typename Number>
843 mapping_info->reinit(cell, unit_points);
847 fe_values = std::make_shared<FEValues<dim, spacedim>>(
851 std::vector<
Point<dim>>(unit_points.begin(), unit_points.end())),
853 fe_values->reinit(cell);
857 std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
860 values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
863 numbers::signaling_nan<gradient_type>());
868 template <
int n_components,
int dim,
int spacedim,
typename Number>
874 const bool precomputed_mapping = mapping_info_on_the_fly.get() ==
nullptr;
875 if (precomputed_mapping)
877 unit_points = mapping_info->get_unit_points();
880 values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
883 numbers::signaling_nan<gradient_type>());
886 if (unit_points.empty())
895 if (solution_renumbered.size() != dofs_per_component)
896 solution_renumbered.resize(dofs_per_component);
897 for (
unsigned int comp = 0; comp < n_components; ++comp)
898 for (
unsigned int i = 0; i < dofs_per_component; ++i)
900 EvaluatorTypeTraits<dim, n_components, Number>::read_value(
901 solution_values[renumber[(component_in_base_element + comp) *
905 solution_renumbered[i]);
909 unit_gradients.resize(unit_points.size(),
910 numbers::signaling_nan<gradient_type>());
912 const std::size_t n_points = unit_points.size();
914 for (
unsigned int i = 0; i < n_points; i += n_lanes)
918 for (
unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
919 for (
unsigned int d = 0;
d < dim; ++
d)
920 vectorized_points[
d][j] = unit_points[i + j][
d];
923 const auto val_and_grad =
928 polynomials_are_hat_functions);
932 for (
unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
934 EvaluatorTypeTraits<dim, n_components, Number>::set_value(
935 val_and_grad.first, j,
values[i + j]);
941 for (
unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
946 Number>::set_gradient(val_and_grad.second,
948 unit_gradients[i + j]);
951 .inverse_jacobians[i + j]
953 unit_gradients[i + j]);
962 Assert(fe_values.get() !=
nullptr,
964 "Not initialized. Please call FEPointEvaluation::reinit()!"));
968 values.resize(unit_points.size());
970 for (
unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
972 const Number value = solution_values[i];
973 for (
unsigned int d = 0;
d < n_components; ++
d)
974 if (nonzero_shape_function_component[i][
d] &&
975 (fe->is_primitive(i) || fe->is_primitive()))
976 for (
unsigned int q = 0; q < unit_points.size(); ++q)
978 EvaluatorTypeTraits<dim, n_components, Number>::access(
979 values[q],
d) += fe_values->shape_value(i, q) * value;
980 else if (nonzero_shape_function_component[i][
d])
981 for (
unsigned int q = 0; q < unit_points.size(); ++q)
983 EvaluatorTypeTraits<dim, n_components, Number>::access(
985 fe_values->shape_value_component(i, q,
d) * value;
993 for (
unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
995 const Number value = solution_values[i];
996 for (
unsigned int d = 0;
d < n_components; ++
d)
997 if (nonzero_shape_function_component[i][
d] &&
998 (fe->is_primitive(i) || fe->is_primitive()))
999 for (
unsigned int q = 0; q < unit_points.size(); ++q)
1001 EvaluatorTypeTraits<dim, n_components, Number>::access(
1002 gradients[q],
d) += fe_values->shape_grad(i, q) * value;
1003 else if (nonzero_shape_function_component[i][
d])
1004 for (
unsigned int q = 0; q < unit_points.size(); ++q)
1006 EvaluatorTypeTraits<dim, n_components, Number>::access(
1008 fe_values->shape_grad_component(i, q,
d) * value;
1016 template <
int n_components,
int dim,
int spacedim,
typename Number>
1022 const bool precomputed_mapping = mapping_info_on_the_fly.get() ==
nullptr;
1023 if (precomputed_mapping)
1025 unit_points = mapping_info->get_unit_points();
1028 values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
1031 numbers::signaling_nan<gradient_type>());
1034 if (unit_points.size() == 0)
1036 std::fill(solution_values.
begin(), solution_values.
end(), 0.0);
1052 if (solution_renumbered_vectorized.size() != dofs_per_component)
1053 solution_renumbered_vectorized.resize(dofs_per_component);
1055 solution_renumbered_vectorized.fill(
1061 const std::size_t n_points = unit_points.size();
1063 for (
unsigned int i = 0; i < n_points; i += n_lanes)
1067 for (
unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1068 for (
unsigned int d = 0;
d < dim; ++
d)
1069 vectorized_points[
d][j] = unit_points[i + j][
d];
1082 for (
unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1084 EvaluatorTypeTraits<dim, n_components, Number>::get_value(
1085 value, j,
values[i + j]);
1087 for (
unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1090 mapping_info->get_mapping_data().inverse_jacobians[i + j],
1103 solution_renumbered_vectorized);
1107 std::fill(solution_values.
begin(), solution_values.
end(), Number());
1108 for (
unsigned int comp = 0; comp < n_components; ++comp)
1109 for (
unsigned int i = 0; i < dofs_per_component; ++i)
1112 internal::FEPointEvaluation::
1113 EvaluatorTypeTraits<dim, n_components, VectorizedArray<Number>>::
1114 write_value(result, comp, solution_renumbered_vectorized[i]);
1115 for (
unsigned int lane = n_lanes / 2; lane > 0; lane /= 2)
1116 for (
unsigned int j = 0; j < lane; ++j)
1117 result[j] += result[lane + j];
1118 solution_values[renumber[comp * dofs_per_component + i]] =
1127 Assert(fe_values.get() !=
nullptr,
1129 "Not initialized. Please call FEPointEvaluation::reinit()!"));
1130 std::fill(solution_values.
begin(), solution_values.
end(), 0.0);
1135 for (
unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1137 for (
unsigned int d = 0;
d < n_components; ++
d)
1138 if (nonzero_shape_function_component[i][
d] &&
1139 (fe->is_primitive(i) || fe->is_primitive()))
1140 for (
unsigned int q = 0; q < unit_points.size(); ++q)
1141 solution_values[i] +=
1142 fe_values->shape_value(i, q) *
1146 else if (nonzero_shape_function_component[i][
d])
1147 for (
unsigned int q = 0; q < unit_points.size(); ++q)
1148 solution_values[i] +=
1149 fe_values->shape_value_component(i, q,
d) *
1159 for (
unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1161 for (
unsigned int d = 0;
d < n_components; ++
d)
1162 if (nonzero_shape_function_component[i][
d] &&
1163 (fe->is_primitive(i) || fe->is_primitive()))
1164 for (
unsigned int q = 0; q < unit_points.size(); ++q)
1165 solution_values[i] +=
1166 fe_values->shape_grad(i, q) *
1170 else if (nonzero_shape_function_component[i][
d])
1171 for (
unsigned int q = 0; q < unit_points.size(); ++q)
1172 solution_values[i] +=
1173 fe_values->shape_grad_component(i, q,
d) *
1184 template <
int n_components,
int dim,
int spacedim,
typename Number>
1188 const unsigned int point_index)
const
1191 return values[point_index];
1196 template <
int n_components,
int dim,
int spacedim,
typename Number>
1200 const unsigned int point_index)
const
1208 template <
int n_components,
int dim,
int spacedim,
typename Number>
1212 const unsigned int point_index)
const
1215 ExcMessage(
"Unit gradients are currently only implemented for tensor "
1216 "product finite elements combined with MappingQ "
1219 return unit_gradients[point_index];
1224 template <
int n_components,
int dim,
int spacedim,
typename Number>
1228 const unsigned int point_index)
1231 values[point_index] = value;
1236 template <
int n_components,
int dim,
int spacedim,
typename Number>
1240 const unsigned int point_index)
1248 template <
int n_components,
int dim,
int spacedim,
typename Number>
1251 const unsigned int point_index)
const
1254 mapping_info->get_mapping_data().jacobians.size());
1255 return mapping_info->get_mapping_data().jacobians[point_index];
1260 template <
int n_components,
int dim,
int spacedim,
typename Number>
1263 const unsigned int point_index)
const
1266 mapping_info->get_mapping_data().inverse_jacobians.size());
1267 return mapping_info->get_mapping_data().inverse_jacobians[point_index];
1272 template <
int n_components,
int dim,
int spacedim,
typename Number>
1275 const unsigned int point_index)
const
1278 mapping_info->get_mapping_data().quadrature_points.size());
1279 return mapping_info->get_mapping_data().quadrature_points[point_index];
1284 template <
int n_components,
int dim,
int spacedim,
typename Number>
1287 const unsigned int point_index)
const
1290 return unit_points[point_index];
std::vector< Polynomials::Polynomial< double > > poly
std::vector< value_type > solution_renumbered
const UpdateFlags update_flags
void submit_gradient(const gradient_type &, const unsigned int point_index)
bool polynomials_are_hat_functions
std::vector< gradient_type > gradients
void setup(const unsigned int first_selected_component)
const value_type & get_value(const unsigned int point_index) const
std::vector< Point< dim > > unit_points
Point< dim > unit_point(const unsigned int point_index) const
unsigned int component_in_base_element
SmartPointer< NonMatching::MappingInfo< dim, spacedim > > mapping_info
SmartPointer< const Mapping< dim, spacedim > > mapping
void evaluate(const ArrayView< const Number > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
std::vector< std::array< bool, n_components > > nonzero_shape_function_component
const gradient_type & get_unit_gradient(const unsigned int point_index) const
DerivativeForm< 1, spacedim, dim > inverse_jacobian(const unsigned int point_index) const
Point< spacedim > real_point(const unsigned int point_index) const
std::vector< gradient_type > unit_gradients
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, n_components, Number >::value_type value_type
const gradient_type & get_gradient(const unsigned int point_index) const
void reinit(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< dim >> &unit_points)
FEPointEvaluation(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
std::unique_ptr< NonMatching::MappingInfo< dim, spacedim > > mapping_info_on_the_fly
DerivativeForm< 1, dim, spacedim > jacobian(const unsigned int point_index) const
AlignedVector< typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, n_components, VectorizedArray< Number > >::value_type > solution_renumbered_vectorized
std::vector< value_type > values
void submit_value(const value_type &value, const unsigned int point_index)
void integrate(const ArrayView< Number > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags)
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, n_components, Number >::gradient_type gradient_type
unsigned int dofs_per_component
std::vector< unsigned int > renumber
std::shared_ptr< FEValues< dim, spacedim > > fe_values
SmartPointer< const FiniteElement< dim > > fe
static constexpr std::size_t size()
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_values
Shape function values.
@ update_inverse_jacobians
Volume element.
@ update_gradients
Shape function gradients.
EvaluationFlags
The EvaluationFlags enum.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< Polynomials::Polynomial< double > > get_polynomial_space(const FiniteElement< dim, spacedim > &fe)
bool is_fast_path_supported(const FiniteElement< dim, spacedim > &fe, const unsigned int base_element_number)
void integrate_add_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double >> &poly, const Number2 &value, const Tensor< 1, dim, Number2 > &gradient, const Point< dim, Number > &p, AlignedVector< Number2 > &values, const std::vector< unsigned int > &renumber={})
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
static Number2 & access(Number2 &value, const unsigned int)
static void get_value(VectorizedArray< Number > &value, const unsigned int vector_lane, const value_type &result)
static void get_gradient(Tensor< 1, 1, VectorizedArray< Number >> &value, const unsigned int vector_lane, const gradient_type &result)
static void set_value(const VectorizedArray< Number > &value, const unsigned int vector_lane, value_type &result)
static void write_value(Number &vector_entry, const unsigned int, const value_type &result)
static void set_gradient(const Tensor< 1, 1, VectorizedArray< Number >> &value, const unsigned int vector_lane, gradient_type &result)
static const Number2 & access(const Number2 &value, const unsigned int)
static void read_value(const Number vector_entry, const unsigned int, value_type &result)
static Number2 & access(Number2 &value, const unsigned int)
static void get_value(VectorizedArray< Number > &value, const unsigned int vector_lane, const value_type &result)
static void set_gradient(const Tensor< 1, dim, VectorizedArray< Number >> &value, const unsigned int vector_lane, gradient_type &result)
static void read_value(const Number vector_entry, const unsigned int, value_type &result)
static const Number2 & access(const Number2 &value, const unsigned int)
static void get_gradient(Tensor< 1, dim, VectorizedArray< Number >> &value, const unsigned int vector_lane, const gradient_type &result)
static void write_value(Number &vector_entry, const unsigned int, const value_type &result)
static void set_value(const VectorizedArray< Number > &value, const unsigned int vector_lane, value_type &result)
static void get_value(Tensor< 1, dim, VectorizedArray< Number >> &value, const unsigned int vector_lane, const value_type &result)
static void write_value(Number &vector_entry, const unsigned int component, const value_type &result)
static const Number & access(const value_type &value, const unsigned int component)
static const Tensor< 1, dim, Number > & access(const gradient_type &value, const unsigned int component)
static void get_gradient(Tensor< 1, dim, Tensor< 1, dim, VectorizedArray< Number >>> &value, const unsigned int vector_lane, const gradient_type &result)
static void set_value(const Tensor< 1, dim, VectorizedArray< Number >> &value, const unsigned int vector_lane, value_type &result)
static void set_gradient(const Tensor< 1, dim, Tensor< 1, dim, VectorizedArray< Number >>> &value, const unsigned int vector_lane, gradient_type &result)
static Tensor< 1, dim, Number > & access(gradient_type &value, const unsigned int component)
static Number & access(value_type &value, const unsigned int component)
static void read_value(const Number vector_entry, const unsigned int component, value_type &result)
static void get_value(Tensor< 1, n_components, VectorizedArray< Number >> &value, const unsigned int vector_lane, const value_type &result)
static void set_value(const Tensor< 1, n_components, VectorizedArray< Number >> &value, const unsigned int vector_lane, value_type &result)
static const Number2 & access(const Tensor< 1, n_components, Number2 > &value, const unsigned int component)
static Number2 & access(Tensor< 1, n_components, Number2 > &value, const unsigned int component)
static void set_gradient(const Tensor< 1, dim, Tensor< 1, n_components, VectorizedArray< Number >>> &value, const unsigned int vector_lane, gradient_type &result)
Tensor< 1, n_components, Number > value_type
static void write_value(Number &vector_entry, const unsigned int component, const value_type &result)
Tensor< 1, n_components, Tensor< 1, dim, Number > > gradient_type
static void read_value(const Number vector_entry, const unsigned int component, value_type &result)
static void get_gradient(Tensor< 1, dim, Tensor< 1, n_components, VectorizedArray< Number >>> &value, const unsigned int vector_lane, const gradient_type &result)
unsigned int dofs_per_component_on_cell
void reinit(const Quadrature< dim_q > &quad, const FiniteElement< dim, spacedim > &fe_dim, const unsigned int base_element=0)
std::vector< unsigned int > lexicographic_numbering