Reference documentation for deal.II version GIT bdf8bf8f35 2023-03-27 16:55:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_point_evaluation.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2020 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_fe_point_evaluation_h
17 #define dealii_fe_point_evaluation_h
18 
19 #include <deal.II/base/config.h>
20 
25 #include <deal.II/base/tensor.h>
27 
28 #include <deal.II/fe/fe_values.h>
29 #include <deal.II/fe/mapping.h>
30 
34 
36 
38 
39 namespace internal
40 {
42  {
47  template <int dim, int n_components, typename Number>
49  {
52 
53  static void
54  read_value(const Number vector_entry,
55  const unsigned int component,
56  value_type & result)
57  {
58  AssertIndexRange(component, n_components);
59  result[component] = vector_entry;
60  }
61 
62  static void
63  write_value(Number & vector_entry,
64  const unsigned int component,
65  const value_type & result)
66  {
67  AssertIndexRange(component, n_components);
68  vector_entry = result[component];
69  }
70 
71  static void
73  const Tensor<1, dim, Tensor<1, n_components, VectorizedArray<Number>>>
74  & value,
75  const unsigned int vector_lane,
76  gradient_type & result)
77  {
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];
81  }
82 
83  static void
85  Tensor<1, dim, Tensor<1, n_components, VectorizedArray<Number>>> &value,
86  const unsigned int vector_lane,
87  const gradient_type &result)
88  {
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];
92  }
93 
94  static void
95  set_value(const Tensor<1, n_components, VectorizedArray<Number>> &value,
96  const unsigned int vector_lane,
97  value_type & result)
98  {
99  for (unsigned int i = 0; i < n_components; ++i)
100  result[i] = value[i][vector_lane];
101  }
102 
103  static void
104  get_value(Tensor<1, n_components, VectorizedArray<Number>> &value,
105  const unsigned int vector_lane,
106  const value_type & result)
107  {
108  for (unsigned int i = 0; i < n_components; ++i)
109  value[i][vector_lane] = result[i];
110  }
111 
112  template <typename Number2>
113  static Number2 &
115  const unsigned int component)
116  {
117  return value[component];
118  }
119 
120  template <typename Number2>
121  static const Number2 &
123  const unsigned int component)
124  {
125  return value[component];
126  }
127  };
128 
129  template <int dim, typename Number>
130  struct EvaluatorTypeTraits<dim, 1, Number>
131  {
132  using value_type = Number;
134 
135  static void
136  read_value(const Number vector_entry,
137  const unsigned int,
138  value_type &result)
139  {
140  result = vector_entry;
141  }
142 
143  static void
144  write_value(Number &vector_entry,
145  const unsigned int,
146  const value_type &result)
147  {
148  vector_entry = result;
149  }
150 
151  static void
153  const unsigned int vector_lane,
154  gradient_type & result)
155  {
156  for (unsigned int d = 0; d < dim; ++d)
157  result[d] = value[d][vector_lane];
158  }
159 
160  static void
162  const unsigned int vector_lane,
163  const gradient_type & result)
164  {
165  for (unsigned int d = 0; d < dim; ++d)
166  value[d][vector_lane] = result[d];
167  }
168 
169  static void
171  const unsigned int vector_lane,
172  value_type & result)
173  {
174  result = value[vector_lane];
175  }
176 
177  static void
179  const unsigned int vector_lane,
180  const value_type & result)
181  {
182  value[vector_lane] = result;
183  }
184 
185  template <typename Number2>
186  static Number2 &
187  access(Number2 &value, const unsigned int)
188  {
189  return value;
190  }
191 
192  template <typename Number2>
193  static const Number2 &
194  access(const Number2 &value, const unsigned int)
195  {
196  return value;
197  }
198  };
199 
200  template <int dim, typename Number>
201  struct EvaluatorTypeTraits<dim, dim, Number>
202  {
205 
206  static void
207  read_value(const Number vector_entry,
208  const unsigned int component,
209  value_type & result)
210  {
211  result[component] = vector_entry;
212  }
213 
214  static void
215  write_value(Number & vector_entry,
216  const unsigned int component,
217  const value_type & result)
218  {
219  vector_entry = result[component];
220  }
221 
222  static void
224  const Tensor<1, dim, Tensor<1, dim, VectorizedArray<Number>>> &value,
225  const unsigned int vector_lane,
226  gradient_type & result)
227  {
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];
231  }
232 
233  static void
235  Tensor<1, dim, Tensor<1, dim, VectorizedArray<Number>>> &value,
236  const unsigned int vector_lane,
237  const gradient_type & result)
238  {
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];
242  }
243 
244  static void
245  set_value(const Tensor<1, dim, VectorizedArray<Number>> &value,
246  const unsigned int vector_lane,
247  value_type & result)
248  {
249  for (unsigned int i = 0; i < dim; ++i)
250  result[i] = value[i][vector_lane];
251  }
252 
253  static void
255  const unsigned int vector_lane,
256  const value_type & result)
257  {
258  for (unsigned int i = 0; i < dim; ++i)
259  value[i][vector_lane] = result[i];
260  }
261 
262  static Number &
263  access(value_type &value, const unsigned int component)
264  {
265  return value[component];
266  }
267 
268  static const Number &
269  access(const value_type &value, const unsigned int component)
270  {
271  return value[component];
272  }
273 
274  static Tensor<1, dim, Number> &
275  access(gradient_type &value, const unsigned int component)
276  {
277  return value[component];
278  }
279 
280  static const Tensor<1, dim, Number> &
281  access(const gradient_type &value, const unsigned int component)
282  {
283  return value[component];
284  }
285  };
286 
287  template <typename Number>
288  struct EvaluatorTypeTraits<1, 1, Number>
289  {
290  using value_type = Number;
292 
293  static void
294  read_value(const Number vector_entry,
295  const unsigned int,
296  value_type &result)
297  {
298  result = vector_entry;
299  }
300 
301  static void
302  write_value(Number &vector_entry,
303  const unsigned int,
304  const value_type &result)
305  {
306  vector_entry = result;
307  }
308 
309  static void
311  const unsigned int vector_lane,
312  gradient_type & result)
313  {
314  result[0] = value[0][vector_lane];
315  }
316 
317  static void
319  const unsigned int vector_lane,
320  const gradient_type & result)
321  {
322  value[0][vector_lane] = result[0];
323  }
324 
325  static void
327  const unsigned int vector_lane,
328  value_type & result)
329  {
330  result = value[vector_lane];
331  }
332 
333  static void
335  const unsigned int vector_lane,
336  const value_type & result)
337  {
338  value[vector_lane] = result;
339  }
340 
341  template <typename Number2>
342  static Number2 &
343  access(Number2 &value, const unsigned int)
344  {
345  return value;
346  }
347 
348  template <typename Number2>
349  static const Number2 &
350  access(const Number2 &value, const unsigned int)
351  {
352  return value;
353  }
354  };
355 
356  template <int dim, int spacedim>
357  bool
359  const unsigned int base_element_number);
360 
361  template <int dim, int spacedim>
362  bool
364 
365  template <int dim, int spacedim>
366  std::vector<Polynomials::Polynomial<double>>
368  } // namespace FEPointEvaluation
369 } // namespace internal
370 
371 
372 
403 template <int n_components,
404  int dim,
405  int spacedim = dim,
406  typename Number = double>
408 {
409 public:
414 
434  const FiniteElement<dim> &fe,
436  const unsigned int first_selected_component = 0);
437 
455  const FiniteElement<dim> & fe,
456  const unsigned int first_selected_component = 0);
457 
469  void
471  const ArrayView<const Point<dim>> &unit_points);
472 
477  void
478  reinit(const unsigned int cell_index);
479 
484  void
485  reinit(const unsigned int cell_index, const unsigned int face_number);
486 
498  void
499  evaluate(const ArrayView<const Number> & solution_values,
500  const EvaluationFlags::EvaluationFlags &evaluation_flags);
501 
525  void
526  integrate(const ArrayView<Number> & solution_values,
527  const EvaluationFlags::EvaluationFlags &integration_flags);
528 
536  const value_type &
537  get_value(const unsigned int point_index) const;
538 
547  void
548  submit_value(const value_type &value, const unsigned int point_index);
549 
559  const gradient_type &
560  get_gradient(const unsigned int point_index) const;
561 
571  const gradient_type &
572  get_unit_gradient(const unsigned int point_index) const;
573 
582  void
583  submit_gradient(const gradient_type &, const unsigned int point_index);
584 
591  jacobian(const unsigned int point_index) const;
592 
600  inverse_jacobian(const unsigned int point_index) const;
601 
607  Number
608  JxW(const unsigned int point_index) const;
609 
616  normal_vector(const unsigned int point_index) const;
617 
623  real_point(const unsigned int point_index) const;
624 
629  Point<dim>
630  unit_point(const unsigned int point_index) const;
631 
638  quadrature_point_indices() const;
639 
640 private:
649  void
650  setup(const unsigned int first_selected_component);
651 
655  const unsigned int n_q_points;
656 
661 
666 
671  std::vector<Polynomials::Polynomial<double>> poly;
672 
677 
682  std::vector<unsigned int> renumber;
683 
690  std::vector<value_type> solution_renumbered;
691 
699  dim,
700  n_components,
703 
707  std::vector<value_type> values;
708 
712  std::vector<gradient_type> unit_gradients;
713 
717  std::vector<gradient_type> gradients;
718 
723  unsigned int dofs_per_component;
724 
729 
735  std::vector<std::array<bool, n_components>> nonzero_shape_function_component;
736 
741 
745  std::shared_ptr<FEValues<dim, spacedim>> fe_values;
746 
750  std::unique_ptr<NonMatching::MappingInfo<dim, spacedim>>
752 
758 
762  unsigned int current_cell_index;
763 
767  unsigned int current_face_number;
768 
772  bool fast_path;
773 
778 };
779 
780 // ----------------------- template and inline function ----------------------
781 
782 
783 template <int n_components, int dim, int spacedim, typename Number>
785  const Mapping<dim> & mapping,
786  const FiniteElement<dim> &fe,
787  const UpdateFlags update_flags,
788  const unsigned int first_selected_component)
789  : n_q_points(numbers::invalid_unsigned_int)
790  , mapping(&mapping)
791  , fe(&fe)
792  , update_flags(update_flags)
793  , mapping_info_on_the_fly(
794  std::make_unique<NonMatching::MappingInfo<dim, spacedim>>(mapping,
795  update_flags))
796  , mapping_info(mapping_info_on_the_fly.get())
797  , current_cell_index(numbers::invalid_unsigned_int)
798  , current_face_number(numbers::invalid_unsigned_int)
799  , is_reinitialized(false)
800 {
801  setup(first_selected_component);
802 }
803 
804 
805 
806 template <int n_components, int dim, int spacedim, typename Number>
809  const FiniteElement<dim> & fe,
810  const unsigned int first_selected_component)
811  : n_q_points(numbers::invalid_unsigned_int)
812  , mapping(&mapping_info.get_mapping())
813  , fe(&fe)
814  , update_flags(mapping_info.get_update_flags())
815  , mapping_info(&mapping_info)
816  , current_cell_index(numbers::invalid_unsigned_int)
817  , current_face_number(numbers::invalid_unsigned_int)
818  , is_reinitialized(false)
819 {
820  setup(first_selected_component);
821 }
822 
823 
824 
825 template <int n_components, int dim, int spacedim, typename Number>
826 void
828  const unsigned int first_selected_component)
829 {
830  AssertIndexRange(first_selected_component + n_components,
831  fe->n_components() + 1);
832 
833  bool same_base_element = true;
834  unsigned int base_element_number = 0;
835  component_in_base_element = 0;
836  unsigned int component = 0;
837  for (; base_element_number < fe->n_base_elements(); ++base_element_number)
838  if (component + fe->element_multiplicity(base_element_number) >
839  first_selected_component)
840  {
841  if (first_selected_component + n_components >
842  component + fe->element_multiplicity(base_element_number))
843  same_base_element = false;
844  component_in_base_element = first_selected_component - component;
845  break;
846  }
847  else
848  component += fe->element_multiplicity(base_element_number);
849 
852  *fe, base_element_number) &&
853  same_base_element)
854  {
856 
857  shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
858  renumber = shape_info.lexicographic_numbering;
859  dofs_per_component = shape_info.dofs_per_component_on_cell;
861  fe->base_element(base_element_number));
862 
863  polynomials_are_hat_functions =
864  (poly.size() == 2 && poly[0].value(0.) == 1. &&
865  poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
866  poly[1].value(1.) == 1.);
867 
868  fast_path = true;
869  }
870  else
871  {
872  nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
873  for (unsigned int d = 0; d < n_components; ++d)
874  {
875  const unsigned int component = first_selected_component + d;
876  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
877  {
878  const bool is_primitive =
879  fe->is_primitive() || fe->is_primitive(i);
880  if (is_primitive)
881  nonzero_shape_function_component[i][d] =
882  (component == fe->system_to_component_index(i).first);
883  else
884  nonzero_shape_function_component[i][d] =
885  (fe->get_nonzero_components(i)[component] == true);
886  }
887  }
888 
889  fast_path = false;
890  }
891 }
892 
893 
894 
895 template <int n_components, int dim, int spacedim, typename Number>
896 void
898  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
899  const ArrayView<const Point<dim>> & unit_points)
900 {
901  // reinit is only allowed for mapping computation on the fly
902  AssertThrow(mapping_info_on_the_fly.get() != nullptr, ExcNotImplemented());
903 
904  mapping_info->reinit(cell, unit_points);
905 
906  if (!fast_path)
907  {
908  fe_values = std::make_shared<FEValues<dim, spacedim>>(
909  *mapping,
910  *fe,
912  std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
913  update_flags);
914  fe_values->reinit(cell);
915  }
916 
917  const_cast<unsigned int &>(n_q_points) = unit_points.size();
918 
919  if (update_flags & update_values)
920  values.resize(n_q_points, numbers::signaling_nan<value_type>());
921  if (update_flags & update_gradients)
922  gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
923 
924  is_reinitialized = true;
925 }
926 
927 
928 
929 template <int n_components, int dim, int spacedim, typename Number>
930 void
932  const unsigned int cell_index)
933 {
934  current_cell_index = cell_index;
935  current_face_number = numbers::invalid_unsigned_int;
936 
937  const_cast<unsigned int &>(n_q_points) =
938  mapping_info->get_unit_points(current_cell_index, current_face_number)
939  .size();
940 
941  if (update_flags & update_values)
942  values.resize(n_q_points, numbers::signaling_nan<value_type>());
943  if (update_flags & update_gradients)
944  gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
945 
946  is_reinitialized = true;
947 }
948 
949 
950 
951 template <int n_components, int dim, int spacedim, typename Number>
952 void
954  const unsigned int cell_index,
955  const unsigned int face_number)
956 {
957  current_cell_index = cell_index;
958  current_face_number = face_number;
959 
960  const_cast<unsigned int &>(n_q_points) =
961  mapping_info->get_unit_points(current_cell_index, current_face_number)
962  .size();
963 
964  if (update_flags & update_values)
965  values.resize(n_q_points, numbers::signaling_nan<value_type>());
966  if (update_flags & update_gradients)
967  gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
968 
969  is_reinitialized = true;
970 }
971 
972 
973 
974 template <int n_components, int dim, int spacedim, typename Number>
975 void
977  const ArrayView<const Number> & solution_values,
978  const EvaluationFlags::EvaluationFlags &evaluation_flag)
979 {
980  if (!is_reinitialized)
982 
983  if (n_q_points == 0)
984  return;
985 
986  AssertDimension(solution_values.size(), fe->dofs_per_cell);
987  if (((evaluation_flag & EvaluationFlags::values) ||
988  (evaluation_flag & EvaluationFlags::gradients)) &&
989  fast_path)
990  {
991  // fast path with tensor product evaluation
992  if (solution_renumbered.size() != dofs_per_component)
993  solution_renumbered.resize(dofs_per_component);
994  for (unsigned int comp = 0; comp < n_components; ++comp)
995  for (unsigned int i = 0; i < dofs_per_component; ++i)
997  EvaluatorTypeTraits<dim, n_components, Number>::read_value(
998  solution_values[renumber[(component_in_base_element + comp) *
999  dofs_per_component +
1000  i]],
1001  comp,
1002  solution_renumbered[i]);
1003 
1004  // unit gradients are currently only implemented with the fast tensor
1005  // path
1006  unit_gradients.resize(n_q_points,
1007  numbers::signaling_nan<gradient_type>());
1008 
1009  const auto unit_points =
1010  mapping_info->get_unit_points(current_cell_index, current_face_number);
1011 
1012  const std::size_t n_points = unit_points.size();
1013  const std::size_t n_lanes = VectorizedArray<Number>::size();
1014  for (unsigned int i = 0; i < n_points; i += n_lanes)
1015  {
1016  // convert to vectorized format
1017  Point<dim, VectorizedArray<Number>> vectorized_points;
1018  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1019  for (unsigned int d = 0; d < dim; ++d)
1020  vectorized_points[d][j] = unit_points[i + j][d];
1021 
1022  // compute
1023  const auto val_and_grad =
1025  poly,
1026  solution_renumbered,
1027  vectorized_points,
1028  polynomials_are_hat_functions);
1029 
1030  // convert back to standard format
1031  if (evaluation_flag & EvaluationFlags::values)
1032  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1034  EvaluatorTypeTraits<dim, n_components, Number>::set_value(
1035  val_and_grad.first, j, values[i + j]);
1036  if (evaluation_flag & EvaluationFlags::gradients)
1037  {
1038  Assert(update_flags & update_gradients ||
1039  update_flags & update_inverse_jacobians,
1040  ExcNotInitialized());
1041  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1042  {
1044  dim,
1045  n_components,
1046  Number>::set_gradient(val_and_grad.second,
1047  j,
1048  unit_gradients[i + j]);
1049  const auto &mapping_data =
1050  mapping_info->get_mapping_data(current_cell_index,
1051  current_face_number);
1053  mapping_data.inverse_jacobians[i + j].transpose(),
1054  unit_gradients[i + j]);
1055  }
1056  }
1057  }
1058  }
1059  else if ((evaluation_flag & EvaluationFlags::values) ||
1060  (evaluation_flag & EvaluationFlags::gradients))
1061  {
1062  // slow path with FEValues
1063  Assert(fe_values.get() != nullptr,
1064  ExcMessage(
1065  "Not initialized. Please call FEPointEvaluation::reinit()!"));
1066 
1067  if (evaluation_flag & EvaluationFlags::values)
1068  {
1069  values.resize(n_q_points);
1070  std::fill(values.begin(), values.end(), value_type());
1071  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1072  {
1073  const Number value = solution_values[i];
1074  for (unsigned int d = 0; d < n_components; ++d)
1075  if (nonzero_shape_function_component[i][d] &&
1076  (fe->is_primitive(i) || fe->is_primitive()))
1077  for (unsigned int q = 0; q < n_q_points; ++q)
1079  EvaluatorTypeTraits<dim, n_components, Number>::access(
1080  values[q], d) += fe_values->shape_value(i, q) * value;
1081  else if (nonzero_shape_function_component[i][d])
1082  for (unsigned int q = 0; q < n_q_points; ++q)
1084  EvaluatorTypeTraits<dim, n_components, Number>::access(
1085  values[q], d) +=
1086  fe_values->shape_value_component(i, q, d) * value;
1087  }
1088  }
1089 
1090  if (evaluation_flag & EvaluationFlags::gradients)
1091  {
1092  gradients.resize(n_q_points);
1093  std::fill(gradients.begin(), gradients.end(), gradient_type());
1094  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1095  {
1096  const Number value = solution_values[i];
1097  for (unsigned int d = 0; d < n_components; ++d)
1098  if (nonzero_shape_function_component[i][d] &&
1099  (fe->is_primitive(i) || fe->is_primitive()))
1100  for (unsigned int q = 0; q < n_q_points; ++q)
1102  EvaluatorTypeTraits<dim, n_components, Number>::access(
1103  gradients[q], d) += fe_values->shape_grad(i, q) * value;
1104  else if (nonzero_shape_function_component[i][d])
1105  for (unsigned int q = 0; q < n_q_points; ++q)
1107  EvaluatorTypeTraits<dim, n_components, Number>::access(
1108  gradients[q], d) +=
1109  fe_values->shape_grad_component(i, q, d) * value;
1110  }
1111  }
1112  }
1113 }
1114 
1115 
1116 
1117 template <int n_components, int dim, int spacedim, typename Number>
1118 void
1120  const ArrayView<Number> & solution_values,
1121  const EvaluationFlags::EvaluationFlags &integration_flags)
1122 {
1123  if (!is_reinitialized)
1125 
1126  if (n_q_points == 0) // no evaluation points provided
1127  {
1128  std::fill(solution_values.begin(), solution_values.end(), 0.0);
1129  return;
1130  }
1131 
1132  AssertDimension(solution_values.size(), fe->dofs_per_cell);
1133  if (((integration_flags & EvaluationFlags::values) ||
1134  (integration_flags & EvaluationFlags::gradients)) &&
1135  fast_path)
1136  {
1137  // fast path with tensor product integration
1138 
1139  if (integration_flags & EvaluationFlags::values)
1140  AssertIndexRange(n_q_points, values.size() + 1);
1141  if (integration_flags & EvaluationFlags::gradients)
1142  AssertIndexRange(n_q_points, gradients.size() + 1);
1143 
1144  if (solution_renumbered_vectorized.size() != dofs_per_component)
1145  solution_renumbered_vectorized.resize(dofs_per_component);
1146  // zero content
1147  solution_renumbered_vectorized.fill(
1149  dim,
1150  n_components,
1152 
1153  const auto unit_points =
1154  mapping_info->get_unit_points(current_cell_index, current_face_number);
1155 
1156  const std::size_t n_points = unit_points.size();
1157  const std::size_t n_lanes = VectorizedArray<Number>::size();
1158  for (unsigned int i = 0; i < n_points; i += n_lanes)
1159  {
1160  // convert to vectorized format
1161  Point<dim, VectorizedArray<Number>> vectorized_points;
1162  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1163  for (unsigned int d = 0; d < dim; ++d)
1164  vectorized_points[d][j] = unit_points[i + j][d];
1165 
1168  value = {};
1169  Tensor<1,
1170  dim,
1172  value_type,
1173  VectorizedArray<Number>>::type>
1174  gradient;
1175 
1176  if (integration_flags & EvaluationFlags::values)
1177  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1179  EvaluatorTypeTraits<dim, n_components, Number>::get_value(
1180  value, j, values[i + j]);
1181  if (integration_flags & EvaluationFlags::gradients)
1182  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1183  {
1184  const auto &mapping_data =
1185  mapping_info->get_mapping_data(current_cell_index,
1186  current_face_number);
1187  gradients[i + j] =
1188  apply_transformation(mapping_data.inverse_jacobians[i + j],
1189  gradients[i + j]);
1192  gradient, j, gradients[i + j]);
1193  }
1194 
1195  // compute
1197  poly,
1198  value,
1199  gradient,
1200  vectorized_points,
1201  solution_renumbered_vectorized);
1202  }
1203 
1204  // add between the lanes and write into the result
1205  std::fill(solution_values.begin(), solution_values.end(), Number());
1206  for (unsigned int comp = 0; comp < n_components; ++comp)
1207  for (unsigned int i = 0; i < dofs_per_component; ++i)
1208  {
1209  VectorizedArray<Number> result;
1210  internal::FEPointEvaluation::
1211  EvaluatorTypeTraits<dim, n_components, VectorizedArray<Number>>::
1212  write_value(result, comp, solution_renumbered_vectorized[i]);
1213  for (unsigned int lane = n_lanes / 2; lane > 0; lane /= 2)
1214  for (unsigned int j = 0; j < lane; ++j)
1215  result[j] += result[lane + j];
1216  solution_values[renumber[comp * dofs_per_component + i]] =
1217  result[0];
1218  }
1219  }
1220  else if ((integration_flags & EvaluationFlags::values) ||
1221  (integration_flags & EvaluationFlags::gradients))
1222  {
1223  // slow path with FEValues
1224 
1225  Assert(fe_values.get() != nullptr,
1226  ExcMessage(
1227  "Not initialized. Please call FEPointEvaluation::reinit()!"));
1228  std::fill(solution_values.begin(), solution_values.end(), 0.0);
1229 
1230  if (integration_flags & EvaluationFlags::values)
1231  {
1232  AssertIndexRange(n_q_points, values.size() + 1);
1233  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1234  {
1235  for (unsigned int d = 0; d < n_components; ++d)
1236  if (nonzero_shape_function_component[i][d] &&
1237  (fe->is_primitive(i) || fe->is_primitive()))
1238  for (unsigned int q = 0; q < n_q_points; ++q)
1239  solution_values[i] +=
1240  fe_values->shape_value(i, q) *
1243  values[q], d);
1244  else if (nonzero_shape_function_component[i][d])
1245  for (unsigned int q = 0; q < n_q_points; ++q)
1246  solution_values[i] +=
1247  fe_values->shape_value_component(i, q, d) *
1250  values[q], d);
1251  }
1252  }
1253 
1254  if (integration_flags & EvaluationFlags::gradients)
1255  {
1256  AssertIndexRange(n_q_points, gradients.size() + 1);
1257  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1258  {
1259  for (unsigned int d = 0; d < n_components; ++d)
1260  if (nonzero_shape_function_component[i][d] &&
1261  (fe->is_primitive(i) || fe->is_primitive()))
1262  for (unsigned int q = 0; q < n_q_points; ++q)
1263  solution_values[i] +=
1264  fe_values->shape_grad(i, q) *
1267  gradients[q], d);
1268  else if (nonzero_shape_function_component[i][d])
1269  for (unsigned int q = 0; q < n_q_points; ++q)
1270  solution_values[i] +=
1271  fe_values->shape_grad_component(i, q, d) *
1274  gradients[q], d);
1275  }
1276  }
1277  }
1278 }
1279 
1280 
1281 
1282 template <int n_components, int dim, int spacedim, typename Number>
1284  value_type &
1286  const unsigned int point_index) const
1287 {
1288  AssertIndexRange(point_index, values.size());
1289  return values[point_index];
1290 }
1291 
1292 
1293 
1294 template <int n_components, int dim, int spacedim, typename Number>
1296  gradient_type &
1298  const unsigned int point_index) const
1299 {
1300  AssertIndexRange(point_index, gradients.size());
1301  return gradients[point_index];
1302 }
1303 
1304 
1305 
1306 template <int n_components, int dim, int spacedim, typename Number>
1308  gradient_type &
1310  const unsigned int point_index) const
1311 {
1312  Assert(fast_path,
1313  ExcMessage("Unit gradients are currently only implemented for tensor "
1314  "product finite elements combined with MappingQ "
1315  "mappings"));
1316  AssertIndexRange(point_index, unit_gradients.size());
1317  return unit_gradients[point_index];
1318 }
1319 
1320 
1321 
1322 template <int n_components, int dim, int spacedim, typename Number>
1323 inline void
1325  const value_type & value,
1326  const unsigned int point_index)
1327 {
1328  AssertIndexRange(point_index, n_q_points);
1329  values[point_index] = value;
1330 }
1331 
1332 
1333 
1334 template <int n_components, int dim, int spacedim, typename Number>
1335 inline void
1337  const gradient_type &gradient,
1338  const unsigned int point_index)
1339 {
1340  AssertIndexRange(point_index, n_q_points);
1341  gradients[point_index] = gradient;
1342 }
1343 
1344 
1345 
1346 template <int n_components, int dim, int spacedim, typename Number>
1349  const unsigned int point_index) const
1350 {
1351  const auto &mapping_data =
1352  mapping_info->get_mapping_data(current_cell_index, current_face_number);
1353  AssertIndexRange(point_index, mapping_data.jacobians.size());
1354  return mapping_data.jacobians[point_index];
1355 }
1356 
1357 
1358 
1359 template <int n_components, int dim, int spacedim, typename Number>
1362  const unsigned int point_index) const
1363 {
1364  const auto &mapping_data =
1365  mapping_info->get_mapping_data(current_cell_index, current_face_number);
1366  AssertIndexRange(point_index, mapping_data.inverse_jacobians.size());
1367  return mapping_data.inverse_jacobians[point_index];
1368 }
1369 
1370 
1371 
1372 template <int n_components, int dim, int spacedim, typename Number>
1373 inline Number
1375  const unsigned int point_index) const
1376 {
1377  const auto &mapping_data =
1378  mapping_info->get_mapping_data(current_cell_index, current_face_number);
1379  AssertIndexRange(point_index, mapping_data.JxW_values.size());
1380  return mapping_data.JxW_values[point_index];
1381 }
1382 
1383 
1384 template <int n_components, int dim, int spacedim, typename Number>
1385 inline Tensor<1, spacedim>
1387  const unsigned int point_index) const
1388 {
1389  const auto &mapping_data =
1390  mapping_info->get_mapping_data(current_cell_index, current_face_number);
1391  AssertIndexRange(point_index, mapping_data.normal_vectors.size());
1392  return mapping_data.normal_vectors[point_index];
1393 }
1394 
1395 
1396 
1397 template <int n_components, int dim, int spacedim, typename Number>
1398 inline Point<spacedim>
1400  const unsigned int point_index) const
1401 {
1402  const auto &mapping_data =
1403  mapping_info->get_mapping_data(current_cell_index, current_face_number);
1404  AssertIndexRange(point_index, mapping_data.quadrature_points.size());
1405  return mapping_data.quadrature_points[point_index];
1406 }
1407 
1408 
1409 
1410 template <int n_components, int dim, int spacedim, typename Number>
1411 inline Point<dim>
1413  const unsigned int point_index) const
1414 {
1415  AssertIndexRange(point_index, n_q_points);
1416  const auto unit_points =
1417  mapping_info->get_unit_points(current_cell_index, current_face_number);
1418  return unit_points[point_index];
1419 }
1420 
1421 
1422 
1423 template <int n_components, int dim, int spacedim, typename Number>
1427 {
1428  return {0U, n_q_points};
1429 }
1430 
1432 
1433 #endif
iterator begin() const
Definition: array_view.h:594
iterator end() const
Definition: array_view.h:603
std::size_t size() const
Definition: array_view.h:576
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
std::vector< Polynomials::Polynomial< double > > poly
Number JxW(const unsigned int point_index) const
std::vector< value_type > solution_renumbered
const UpdateFlags update_flags
void submit_gradient(const gradient_type &, const unsigned int point_index)
std_cxx20::ranges::iota_view< unsigned int, unsigned int > quadrature_point_indices() const
std::vector< gradient_type > gradients
void setup(const unsigned int first_selected_component)
const value_type & get_value(const unsigned int point_index) const
Tensor< 1, spacedim > normal_vector(const unsigned int point_index) const
Point< dim > unit_point(const unsigned int point_index) const
unsigned int component_in_base_element
unsigned int current_face_number
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
const unsigned int n_q_points
std::vector< value_type > values
unsigned int current_cell_index
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
Definition: point.h:110
Definition: tensor.h:516
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
unsigned int cell_index
Definition: grid_tools.cc:1191
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
UpdateFlags
@ update_values
Shape function values.
@ update_inverse_jacobians
Volume element.
@ update_gradients
Shape function gradients.
EvaluationFlags
The EvaluationFlags enum.
static const char U
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={})
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
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 const unsigned int invalid_unsigned_int
Definition: types.h:213
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:46
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)
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
Definition: shape_info.h:435