Reference documentation for deal.II version GIT d8dacc551e 2022-08-19 06:50:03+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>
30 #include <deal.II/fe/mapping_q.h>
31 
35 
37 
39 
40 namespace internal
41 {
43  {
48  template <int dim, int n_components, typename Number>
50  {
53 
54  static void
55  read_value(const Number vector_entry,
56  const unsigned int component,
57  value_type & result)
58  {
59  AssertIndexRange(component, n_components);
60  result[component] = vector_entry;
61  }
62 
63  static void
64  write_value(Number & vector_entry,
65  const unsigned int component,
66  const value_type & result)
67  {
68  AssertIndexRange(component, n_components);
69  vector_entry = result[component];
70  }
71 
72  static void
74  const Tensor<1, dim, Tensor<1, n_components, VectorizedArray<Number>>>
75  & value,
76  const unsigned int vector_lane,
77  gradient_type & result)
78  {
79  for (unsigned int i = 0; i < n_components; ++i)
80  for (unsigned int d = 0; d < dim; ++d)
81  result[i][d] = value[d][i][vector_lane];
82  }
83 
84  static void
86  Tensor<1, dim, Tensor<1, n_components, VectorizedArray<Number>>> &value,
87  const unsigned int vector_lane,
88  const gradient_type &result)
89  {
90  for (unsigned int i = 0; i < n_components; ++i)
91  for (unsigned int d = 0; d < dim; ++d)
92  value[d][i][vector_lane] = result[i][d];
93  }
94 
95  static void
96  set_value(const Tensor<1, n_components, VectorizedArray<Number>> &value,
97  const unsigned int vector_lane,
98  value_type & result)
99  {
100  for (unsigned int i = 0; i < n_components; ++i)
101  result[i] = value[i][vector_lane];
102  }
103 
104  static void
105  get_value(Tensor<1, n_components, VectorizedArray<Number>> &value,
106  const unsigned int vector_lane,
107  const value_type & result)
108  {
109  for (unsigned int i = 0; i < n_components; ++i)
110  value[i][vector_lane] = result[i];
111  }
112 
113  template <typename Number2>
114  static Number2 &
116  const unsigned int component)
117  {
118  return value[component];
119  }
120 
121  template <typename Number2>
122  static const Number2 &
124  const unsigned int component)
125  {
126  return value[component];
127  }
128  };
129 
130  template <int dim, typename Number>
131  struct EvaluatorTypeTraits<dim, 1, Number>
132  {
133  using value_type = Number;
135 
136  static void
137  read_value(const Number vector_entry,
138  const unsigned int,
139  value_type &result)
140  {
141  result = vector_entry;
142  }
143 
144  static void
145  write_value(Number &vector_entry,
146  const unsigned int,
147  const value_type &result)
148  {
149  vector_entry = result;
150  }
151 
152  static void
154  const unsigned int vector_lane,
155  gradient_type & result)
156  {
157  for (unsigned int d = 0; d < dim; ++d)
158  result[d] = value[d][vector_lane];
159  }
160 
161  static void
163  const unsigned int vector_lane,
164  const gradient_type & result)
165  {
166  for (unsigned int d = 0; d < dim; ++d)
167  value[d][vector_lane] = result[d];
168  }
169 
170  static void
172  const unsigned int vector_lane,
173  value_type & result)
174  {
175  result = value[vector_lane];
176  }
177 
178  static void
180  const unsigned int vector_lane,
181  const value_type & result)
182  {
183  value[vector_lane] = result;
184  }
185 
186  template <typename Number2>
187  static Number2 &
188  access(Number2 &value, const unsigned int)
189  {
190  return value;
191  }
192 
193  template <typename Number2>
194  static const Number2 &
195  access(const Number2 &value, const unsigned int)
196  {
197  return value;
198  }
199  };
200 
201  template <int dim, typename Number>
202  struct EvaluatorTypeTraits<dim, dim, Number>
203  {
206 
207  static void
208  read_value(const Number vector_entry,
209  const unsigned int component,
210  value_type & result)
211  {
212  result[component] = vector_entry;
213  }
214 
215  static void
216  write_value(Number & vector_entry,
217  const unsigned int component,
218  const value_type & result)
219  {
220  vector_entry = result[component];
221  }
222 
223  static void
225  const Tensor<1, dim, Tensor<1, dim, VectorizedArray<Number>>> &value,
226  const unsigned int vector_lane,
227  gradient_type & result)
228  {
229  for (unsigned int i = 0; i < dim; ++i)
230  for (unsigned int d = 0; d < dim; ++d)
231  result[i][d] = value[d][i][vector_lane];
232  }
233 
234  static void
236  Tensor<1, dim, Tensor<1, dim, VectorizedArray<Number>>> &value,
237  const unsigned int vector_lane,
238  const gradient_type & result)
239  {
240  for (unsigned int i = 0; i < dim; ++i)
241  for (unsigned int d = 0; d < dim; ++d)
242  value[d][i][vector_lane] = result[i][d];
243  }
244 
245  static void
246  set_value(const Tensor<1, dim, VectorizedArray<Number>> &value,
247  const unsigned int vector_lane,
248  value_type & result)
249  {
250  for (unsigned int i = 0; i < dim; ++i)
251  result[i] = value[i][vector_lane];
252  }
253 
254  static void
256  const unsigned int vector_lane,
257  const value_type & result)
258  {
259  for (unsigned int i = 0; i < dim; ++i)
260  value[i][vector_lane] = result[i];
261  }
262 
263  static Number &
264  access(value_type &value, const unsigned int component)
265  {
266  return value[component];
267  }
268 
269  static const Number &
270  access(const value_type &value, const unsigned int component)
271  {
272  return value[component];
273  }
274 
275  static Tensor<1, dim, Number> &
276  access(gradient_type &value, const unsigned int component)
277  {
278  return value[component];
279  }
280 
281  static const Tensor<1, dim, Number> &
282  access(const gradient_type &value, const unsigned int component)
283  {
284  return value[component];
285  }
286  };
287 
288  template <typename Number>
289  struct EvaluatorTypeTraits<1, 1, Number>
290  {
291  using value_type = Number;
293 
294  static void
295  read_value(const Number vector_entry,
296  const unsigned int,
297  value_type &result)
298  {
299  result = vector_entry;
300  }
301 
302  static void
303  write_value(Number &vector_entry,
304  const unsigned int,
305  const value_type &result)
306  {
307  vector_entry = result;
308  }
309 
310  static void
312  const unsigned int vector_lane,
313  gradient_type & result)
314  {
315  result[0] = value[0][vector_lane];
316  }
317 
318  static void
320  const unsigned int vector_lane,
321  const gradient_type & result)
322  {
323  value[0][vector_lane] = result[0];
324  }
325 
326  static void
328  const unsigned int vector_lane,
329  value_type & result)
330  {
331  result = value[vector_lane];
332  }
333 
334  static void
336  const unsigned int vector_lane,
337  const value_type & result)
338  {
339  value[vector_lane] = result;
340  }
341 
342  template <typename Number2>
343  static Number2 &
344  access(Number2 &value, const unsigned int)
345  {
346  return value;
347  }
348 
349  template <typename Number2>
350  static const Number2 &
351  access(const Number2 &value, const unsigned int)
352  {
353  return value;
354  }
355  };
356 
357  template <int dim, int spacedim>
358  bool
360  const unsigned int base_element_number);
361 
362  template <int dim, int spacedim>
363  bool
365 
366  template <int dim, int spacedim>
367  std::vector<Polynomials::Polynomial<double>>
369  } // namespace FEPointEvaluation
370 } // namespace internal
371 
372 
373 
404 template <int n_components,
405  int dim,
406  int spacedim = dim,
407  typename Number = double>
409 {
410 public:
415 
435  const FiniteElement<dim> &fe,
437  const unsigned int first_selected_component = 0);
438 
456  const FiniteElement<dim> & fe,
457  const unsigned int first_selected_component = 0);
458 
470  void
472  const ArrayView<const Point<dim>> &unit_points);
473 
485  void
486  evaluate(const ArrayView<const Number> & solution_values,
487  const EvaluationFlags::EvaluationFlags &evaluation_flags);
488 
512  void
513  integrate(const ArrayView<Number> & solution_values,
514  const EvaluationFlags::EvaluationFlags &integration_flags);
515 
523  const value_type &
524  get_value(const unsigned int point_index) const;
525 
534  void
535  submit_value(const value_type &value, const unsigned int point_index);
536 
546  const gradient_type &
547  get_gradient(const unsigned int point_index) const;
548 
558  const gradient_type &
559  get_unit_gradient(const unsigned int point_index) const;
560 
569  void
570  submit_gradient(const gradient_type &, const unsigned int point_index);
571 
578  jacobian(const unsigned int point_index) const;
579 
587  inverse_jacobian(const unsigned int point_index) const;
588 
594  real_point(const unsigned int point_index) const;
595 
600  Point<dim>
601  unit_point(const unsigned int point_index) const;
602 
603 private:
612  void
613  setup(const unsigned int first_selected_component);
614 
619 
624 
629  std::vector<Polynomials::Polynomial<double>> poly;
630 
635 
640  std::vector<unsigned int> renumber;
641 
648  std::vector<value_type> solution_renumbered;
649 
657  dim,
658  n_components,
661 
665  std::vector<value_type> values;
666 
670  std::vector<gradient_type> unit_gradients;
671 
675  std::vector<gradient_type> gradients;
676 
681  unsigned int dofs_per_component;
682 
687 
693  std::vector<std::array<bool, n_components>> nonzero_shape_function_component;
694 
699 
703  std::shared_ptr<FEValues<dim, spacedim>> fe_values;
704 
708  std::unique_ptr<NonMatching::MappingInfo<dim, spacedim>>
710 
716 
720  std::vector<Point<dim>> unit_points;
721 
725  bool fast_path;
726 };
727 
728 // ----------------------- template and inline function ----------------------
729 
730 
731 template <int n_components, int dim, int spacedim, typename Number>
733  const Mapping<dim> & mapping,
734  const FiniteElement<dim> &fe,
735  const UpdateFlags update_flags,
736  const unsigned int first_selected_component)
737  : mapping(&mapping)
738  , fe(&fe)
739  , update_flags(update_flags)
740  , mapping_info_on_the_fly(
741  std::make_unique<NonMatching::MappingInfo<dim, spacedim>>(mapping,
742  update_flags))
743  , mapping_info(mapping_info_on_the_fly.get())
744 {
745  setup(first_selected_component);
746 }
747 
748 
749 
750 template <int n_components, int dim, int spacedim, typename Number>
753  const FiniteElement<dim> & fe,
754  const unsigned int first_selected_component)
755  : mapping(&mapping_info.get_mapping())
756  , fe(&fe)
757  , update_flags(mapping_info.get_update_flags())
758  , mapping_info(&mapping_info)
759 {
760  setup(first_selected_component);
761 }
762 
763 
764 
765 template <int n_components, int dim, int spacedim, typename Number>
766 void
768  const unsigned int first_selected_component)
769 {
770  AssertIndexRange(first_selected_component + n_components,
771  fe->n_components() + 1);
772 
773  bool same_base_element = true;
774  unsigned int base_element_number = 0;
775  component_in_base_element = 0;
776  unsigned int component = 0;
777  for (; base_element_number < fe->n_base_elements(); ++base_element_number)
778  if (component + fe->element_multiplicity(base_element_number) >
779  first_selected_component)
780  {
781  if (first_selected_component + n_components >
782  component + fe->element_multiplicity(base_element_number))
783  same_base_element = false;
784  component_in_base_element = first_selected_component - component;
785  break;
786  }
787  else
788  component += fe->element_multiplicity(base_element_number);
789 
792  *fe, base_element_number) &&
793  same_base_element)
794  {
796 
797  shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
798  renumber = shape_info.lexicographic_numbering;
799  dofs_per_component = shape_info.dofs_per_component_on_cell;
801  fe->base_element(base_element_number));
802 
803  polynomials_are_hat_functions =
804  (poly.size() == 2 && poly[0].value(0.) == 1. &&
805  poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
806  poly[1].value(1.) == 1.);
807 
808  fast_path = true;
809  }
810  else
811  {
812  nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
813  for (unsigned int d = 0; d < n_components; ++d)
814  {
815  const unsigned int component = first_selected_component + d;
816  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
817  {
818  const bool is_primitive =
819  fe->is_primitive() || fe->is_primitive(i);
820  if (is_primitive)
821  nonzero_shape_function_component[i][d] =
822  (component == fe->system_to_component_index(i).first);
823  else
824  nonzero_shape_function_component[i][d] =
825  (fe->get_nonzero_components(i)[component] == true);
826  }
827  }
828 
829  fast_path = false;
830  }
831 }
832 
833 
834 
835 template <int n_components, int dim, int spacedim, typename Number>
836 void
838  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
839  const ArrayView<const Point<dim>> & unit_points)
840 {
841  // reinit is only allowed for mapping computation on the fly
842  AssertThrow(mapping_info_on_the_fly.get() != nullptr, ExcNotImplemented());
843 
844  mapping_info->reinit(cell, unit_points);
845 
846  if (!fast_path)
847  {
848  fe_values = std::make_shared<FEValues<dim, spacedim>>(
849  *mapping,
850  *fe,
852  std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
853  update_flags);
854  fe_values->reinit(cell);
855  }
856 
857  this->unit_points =
858  std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
859 
860  if (update_flags & update_values)
861  values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
862  if (update_flags & update_gradients)
863  gradients.resize(unit_points.size(),
864  numbers::signaling_nan<gradient_type>());
865 }
866 
867 
868 
869 template <int n_components, int dim, int spacedim, typename Number>
870 void
872  const ArrayView<const Number> & solution_values,
873  const EvaluationFlags::EvaluationFlags &evaluation_flag)
874 {
875  const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr;
876  if (precomputed_mapping)
877  {
878  unit_points = mapping_info->get_unit_points();
879 
880  if (update_flags & update_values)
881  values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
882  if (update_flags & update_gradients)
883  gradients.resize(unit_points.size(),
884  numbers::signaling_nan<gradient_type>());
885  }
886 
887  if (unit_points.empty())
888  return;
889 
890  AssertDimension(solution_values.size(), fe->dofs_per_cell);
891  if (((evaluation_flag & EvaluationFlags::values) ||
892  (evaluation_flag & EvaluationFlags::gradients)) &&
893  fast_path)
894  {
895  // fast path with tensor product evaluation
896  if (solution_renumbered.size() != dofs_per_component)
897  solution_renumbered.resize(dofs_per_component);
898  for (unsigned int comp = 0; comp < n_components; ++comp)
899  for (unsigned int i = 0; i < dofs_per_component; ++i)
901  EvaluatorTypeTraits<dim, n_components, Number>::read_value(
902  solution_values[renumber[(component_in_base_element + comp) *
903  dofs_per_component +
904  i]],
905  comp,
906  solution_renumbered[i]);
907 
908  // unit gradients are currently only implemented with the fast tensor
909  // path
910  unit_gradients.resize(unit_points.size(),
911  numbers::signaling_nan<gradient_type>());
912 
913  const std::size_t n_points = unit_points.size();
914  const std::size_t n_lanes = VectorizedArray<Number>::size();
915  for (unsigned int i = 0; i < n_points; i += n_lanes)
916  {
917  // convert to vectorized format
918  Point<dim, VectorizedArray<Number>> vectorized_points;
919  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
920  for (unsigned int d = 0; d < dim; ++d)
921  vectorized_points[d][j] = unit_points[i + j][d];
922 
923  // compute
924  const auto val_and_grad =
926  poly,
927  solution_renumbered,
928  vectorized_points,
929  polynomials_are_hat_functions);
930 
931  // convert back to standard format
932  if (evaluation_flag & EvaluationFlags::values)
933  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
935  EvaluatorTypeTraits<dim, n_components, Number>::set_value(
936  val_and_grad.first, j, values[i + j]);
937  if (evaluation_flag & EvaluationFlags::gradients)
938  {
939  Assert(update_flags & update_gradients ||
940  update_flags & update_inverse_jacobians,
942  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
943  {
945  dim,
946  n_components,
947  Number>::set_gradient(val_and_grad.second,
948  j,
949  unit_gradients[i + j]);
950  gradients[i + j] =
951  static_cast<typename internal::FEPointEvaluation::
952  EvaluatorTypeTraits<dim,
953  n_components,
954  Number>::gradient_type>(
955  apply_transformation(mapping_info->get_mapping_data()
956  .inverse_jacobians[i + j]
957  .transpose(),
958  unit_gradients[i + j]));
959  }
960  }
961  }
962  }
963  else if ((evaluation_flag & EvaluationFlags::values) ||
964  (evaluation_flag & EvaluationFlags::gradients))
965  {
966  // slow path with FEValues
967  Assert(fe_values.get() != nullptr,
968  ExcMessage(
969  "Not initialized. Please call FEPointEvaluation::reinit()!"));
970 
971  if (evaluation_flag & EvaluationFlags::values)
972  {
973  values.resize(unit_points.size());
974  std::fill(values.begin(), values.end(), value_type());
975  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
976  {
977  const Number value = solution_values[i];
978  for (unsigned int d = 0; d < n_components; ++d)
979  if (nonzero_shape_function_component[i][d] &&
980  (fe->is_primitive(i) || fe->is_primitive()))
981  for (unsigned int q = 0; q < unit_points.size(); ++q)
983  EvaluatorTypeTraits<dim, n_components, Number>::access(
984  values[q], d) += fe_values->shape_value(i, q) * value;
985  else if (nonzero_shape_function_component[i][d])
986  for (unsigned int q = 0; q < unit_points.size(); ++q)
988  EvaluatorTypeTraits<dim, n_components, Number>::access(
989  values[q], d) +=
990  fe_values->shape_value_component(i, q, d) * value;
991  }
992  }
993 
994  if (evaluation_flag & EvaluationFlags::gradients)
995  {
996  gradients.resize(unit_points.size());
997  std::fill(gradients.begin(), gradients.end(), gradient_type());
998  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
999  {
1000  const Number value = solution_values[i];
1001  for (unsigned int d = 0; d < n_components; ++d)
1002  if (nonzero_shape_function_component[i][d] &&
1003  (fe->is_primitive(i) || fe->is_primitive()))
1004  for (unsigned int q = 0; q < unit_points.size(); ++q)
1006  EvaluatorTypeTraits<dim, n_components, Number>::access(
1007  gradients[q], d) += fe_values->shape_grad(i, q) * value;
1008  else if (nonzero_shape_function_component[i][d])
1009  for (unsigned int q = 0; q < unit_points.size(); ++q)
1011  EvaluatorTypeTraits<dim, n_components, Number>::access(
1012  gradients[q], d) +=
1013  fe_values->shape_grad_component(i, q, d) * value;
1014  }
1015  }
1016  }
1017 }
1018 
1019 
1020 
1021 template <int n_components, int dim, int spacedim, typename Number>
1022 void
1024  const ArrayView<Number> & solution_values,
1025  const EvaluationFlags::EvaluationFlags &integration_flags)
1026 {
1027  const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr;
1028  if (precomputed_mapping)
1029  {
1030  unit_points = mapping_info->get_unit_points();
1031 
1032  if (update_flags & update_values)
1033  values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
1034  if (update_flags & update_gradients)
1035  gradients.resize(unit_points.size(),
1036  numbers::signaling_nan<gradient_type>());
1037  }
1038 
1039  if (unit_points.size() == 0) // no evaluation points provided
1040  {
1041  std::fill(solution_values.begin(), solution_values.end(), 0.0);
1042  return;
1043  }
1044 
1045  AssertDimension(solution_values.size(), fe->dofs_per_cell);
1046  if (((integration_flags & EvaluationFlags::values) ||
1047  (integration_flags & EvaluationFlags::gradients)) &&
1048  fast_path)
1049  {
1050  // fast path with tensor product integration
1051 
1052  if (integration_flags & EvaluationFlags::values)
1053  AssertIndexRange(unit_points.size(), values.size() + 1);
1054  if (integration_flags & EvaluationFlags::gradients)
1055  AssertIndexRange(unit_points.size(), gradients.size() + 1);
1056 
1057  if (solution_renumbered_vectorized.size() != dofs_per_component)
1058  solution_renumbered_vectorized.resize(dofs_per_component);
1059  // zero content
1060  solution_renumbered_vectorized.fill(
1062  dim,
1063  n_components,
1065 
1066  const std::size_t n_points = unit_points.size();
1067  const std::size_t n_lanes = VectorizedArray<Number>::size();
1068  for (unsigned int i = 0; i < n_points; i += n_lanes)
1069  {
1070  // convert to vectorized format
1071  Point<dim, VectorizedArray<Number>> vectorized_points;
1072  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1073  for (unsigned int d = 0; d < dim; ++d)
1074  vectorized_points[d][j] = unit_points[i + j][d];
1075 
1078  value = {};
1079  Tensor<1,
1080  dim,
1082  value_type,
1083  VectorizedArray<Number>>::type>
1084  gradient;
1085 
1086  if (integration_flags & EvaluationFlags::values)
1087  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1089  EvaluatorTypeTraits<dim, n_components, Number>::get_value(
1090  value, j, values[i + j]);
1091  if (integration_flags & EvaluationFlags::gradients)
1092  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1093  {
1094  gradients[i + j] =
1095  static_cast<typename internal::FEPointEvaluation::
1096  EvaluatorTypeTraits<dim, n_components, Number>::
1098  mapping_info->get_mapping_data().inverse_jacobians[i + j],
1099  gradients[i + j]));
1102  gradient, j, gradients[i + j]);
1103  }
1104 
1105  // compute
1107  poly,
1108  value,
1109  gradient,
1110  vectorized_points,
1111  solution_renumbered_vectorized);
1112  }
1113 
1114  // add between the lanes and write into the result
1115  std::fill(solution_values.begin(), solution_values.end(), Number());
1116  for (unsigned int comp = 0; comp < n_components; ++comp)
1117  for (unsigned int i = 0; i < dofs_per_component; ++i)
1118  {
1119  VectorizedArray<Number> result;
1120  internal::FEPointEvaluation::
1121  EvaluatorTypeTraits<dim, n_components, VectorizedArray<Number>>::
1122  write_value(result, comp, solution_renumbered_vectorized[i]);
1123  for (unsigned int lane = n_lanes / 2; lane > 0; lane /= 2)
1124  for (unsigned int j = 0; j < lane; ++j)
1125  result[j] += result[lane + j];
1126  solution_values[renumber[comp * dofs_per_component + i]] =
1127  result[0];
1128  }
1129  }
1130  else if ((integration_flags & EvaluationFlags::values) ||
1131  (integration_flags & EvaluationFlags::gradients))
1132  {
1133  // slow path with FEValues
1134 
1135  Assert(fe_values.get() != nullptr,
1136  ExcMessage(
1137  "Not initialized. Please call FEPointEvaluation::reinit()!"));
1138  std::fill(solution_values.begin(), solution_values.end(), 0.0);
1139 
1140  if (integration_flags & EvaluationFlags::values)
1141  {
1142  AssertIndexRange(unit_points.size(), values.size() + 1);
1143  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1144  {
1145  for (unsigned int d = 0; d < n_components; ++d)
1146  if (nonzero_shape_function_component[i][d] &&
1147  (fe->is_primitive(i) || fe->is_primitive()))
1148  for (unsigned int q = 0; q < unit_points.size(); ++q)
1149  solution_values[i] +=
1150  fe_values->shape_value(i, q) *
1153  values[q], d);
1154  else if (nonzero_shape_function_component[i][d])
1155  for (unsigned int q = 0; q < unit_points.size(); ++q)
1156  solution_values[i] +=
1157  fe_values->shape_value_component(i, q, d) *
1160  values[q], d);
1161  }
1162  }
1163 
1164  if (integration_flags & EvaluationFlags::gradients)
1165  {
1166  AssertIndexRange(unit_points.size(), gradients.size() + 1);
1167  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1168  {
1169  for (unsigned int d = 0; d < n_components; ++d)
1170  if (nonzero_shape_function_component[i][d] &&
1171  (fe->is_primitive(i) || fe->is_primitive()))
1172  for (unsigned int q = 0; q < unit_points.size(); ++q)
1173  solution_values[i] +=
1174  fe_values->shape_grad(i, q) *
1177  gradients[q], d);
1178  else if (nonzero_shape_function_component[i][d])
1179  for (unsigned int q = 0; q < unit_points.size(); ++q)
1180  solution_values[i] +=
1181  fe_values->shape_grad_component(i, q, d) *
1184  gradients[q], d);
1185  }
1186  }
1187  }
1188 }
1189 
1190 
1191 
1192 template <int n_components, int dim, int spacedim, typename Number>
1194  value_type &
1196  const unsigned int point_index) const
1197 {
1198  AssertIndexRange(point_index, values.size());
1199  return values[point_index];
1200 }
1201 
1202 
1203 
1204 template <int n_components, int dim, int spacedim, typename Number>
1206  gradient_type &
1208  const unsigned int point_index) const
1209 {
1210  AssertIndexRange(point_index, gradients.size());
1211  return gradients[point_index];
1212 }
1213 
1214 
1215 
1216 template <int n_components, int dim, int spacedim, typename Number>
1218  gradient_type &
1220  const unsigned int point_index) const
1221 {
1222  Assert(fast_path,
1223  ExcMessage("Unit gradients are currently only implemented for tensor "
1224  "product finite elements combined with MappingQ "
1225  "mappings"));
1226  AssertIndexRange(point_index, unit_gradients.size());
1227  return unit_gradients[point_index];
1228 }
1229 
1230 
1231 
1232 template <int n_components, int dim, int spacedim, typename Number>
1233 inline void
1235  const value_type & value,
1236  const unsigned int point_index)
1237 {
1238  AssertIndexRange(point_index, unit_points.size());
1239  values[point_index] = value;
1240 }
1241 
1242 
1243 
1244 template <int n_components, int dim, int spacedim, typename Number>
1245 inline void
1247  const gradient_type &gradient,
1248  const unsigned int point_index)
1249 {
1250  AssertIndexRange(point_index, unit_points.size());
1251  gradients[point_index] = gradient;
1252 }
1253 
1254 
1255 
1256 template <int n_components, int dim, int spacedim, typename Number>
1259  const unsigned int point_index) const
1260 {
1261  AssertIndexRange(point_index,
1262  mapping_info->get_mapping_data().jacobians.size());
1263  return mapping_info->get_mapping_data().jacobians[point_index];
1264 }
1265 
1266 
1267 
1268 template <int n_components, int dim, int spacedim, typename Number>
1271  const unsigned int point_index) const
1272 {
1273  AssertIndexRange(point_index,
1274  mapping_info->get_mapping_data().inverse_jacobians.size());
1275  return mapping_info->get_mapping_data().inverse_jacobians[point_index];
1276 }
1277 
1278 
1279 
1280 template <int n_components, int dim, int spacedim, typename Number>
1281 inline Point<spacedim>
1283  const unsigned int point_index) const
1284 {
1285  AssertIndexRange(point_index,
1286  mapping_info->get_mapping_data().quadrature_points.size());
1287  return mapping_info->get_mapping_data().quadrature_points[point_index];
1288 }
1289 
1290 
1291 
1292 template <int n_components, int dim, int spacedim, typename Number>
1293 inline Point<dim>
1295  const unsigned int point_index) const
1296 {
1297  AssertIndexRange(point_index, unit_points.size());
1298  return unit_points[point_index];
1299 }
1300 
1302 
1303 #endif
iterator begin() const
Definition: array_view.h:585
iterator end() const
Definition: array_view.h:594
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
std::vector< value_type > solution_renumbered
const UpdateFlags update_flags
void submit_gradient(const gradient_type &, const unsigned int point_index)
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
Definition: point.h:111
Definition: tensor.h:503
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
UpdateFlags
@ 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)
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:423