Reference documentation for deal.II version GIT f3cf314fd8 2023-02-08 12:00:02+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 
484  void
485  evaluate(const ArrayView<const Number> & solution_values,
486  const EvaluationFlags::EvaluationFlags &evaluation_flags);
487 
511  void
512  integrate(const ArrayView<Number> & solution_values,
513  const EvaluationFlags::EvaluationFlags &integration_flags);
514 
522  const value_type &
523  get_value(const unsigned int point_index) const;
524 
533  void
534  submit_value(const value_type &value, const unsigned int point_index);
535 
545  const gradient_type &
546  get_gradient(const unsigned int point_index) const;
547 
557  const gradient_type &
558  get_unit_gradient(const unsigned int point_index) const;
559 
568  void
569  submit_gradient(const gradient_type &, const unsigned int point_index);
570 
577  jacobian(const unsigned int point_index) const;
578 
586  inverse_jacobian(const unsigned int point_index) const;
587 
593  real_point(const unsigned int point_index) const;
594 
599  Point<dim>
600  unit_point(const unsigned int point_index) const;
601 
602 private:
611  void
612  setup(const unsigned int first_selected_component);
613 
618 
623 
628  std::vector<Polynomials::Polynomial<double>> poly;
629 
634 
639  std::vector<unsigned int> renumber;
640 
647  std::vector<value_type> solution_renumbered;
648 
656  dim,
657  n_components,
660 
664  std::vector<value_type> values;
665 
669  std::vector<gradient_type> unit_gradients;
670 
674  std::vector<gradient_type> gradients;
675 
680  unsigned int dofs_per_component;
681 
686 
692  std::vector<std::array<bool, n_components>> nonzero_shape_function_component;
693 
698 
702  std::shared_ptr<FEValues<dim, spacedim>> fe_values;
703 
707  std::unique_ptr<NonMatching::MappingInfo<dim, spacedim>>
709 
715 
719  std::vector<Point<dim>> unit_points;
720 
724  bool fast_path;
725 };
726 
727 // ----------------------- template and inline function ----------------------
728 
729 
730 template <int n_components, int dim, int spacedim, typename Number>
732  const Mapping<dim> & mapping,
733  const FiniteElement<dim> &fe,
734  const UpdateFlags update_flags,
735  const unsigned int first_selected_component)
736  : mapping(&mapping)
737  , fe(&fe)
738  , update_flags(update_flags)
739  , mapping_info_on_the_fly(
740  std::make_unique<NonMatching::MappingInfo<dim, spacedim>>(mapping,
741  update_flags))
742  , mapping_info(mapping_info_on_the_fly.get())
743 {
744  setup(first_selected_component);
745 }
746 
747 
748 
749 template <int n_components, int dim, int spacedim, typename Number>
752  const FiniteElement<dim> & fe,
753  const unsigned int first_selected_component)
754  : mapping(&mapping_info.get_mapping())
755  , fe(&fe)
756  , update_flags(mapping_info.get_update_flags())
757  , mapping_info(&mapping_info)
758 {
759  setup(first_selected_component);
760 }
761 
762 
763 
764 template <int n_components, int dim, int spacedim, typename Number>
765 void
767  const unsigned int first_selected_component)
768 {
769  AssertIndexRange(first_selected_component + n_components,
770  fe->n_components() + 1);
771 
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)
779  {
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;
784  break;
785  }
786  else
787  component += fe->element_multiplicity(base_element_number);
788 
791  *fe, base_element_number) &&
792  same_base_element)
793  {
795 
796  shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
797  renumber = shape_info.lexicographic_numbering;
798  dofs_per_component = shape_info.dofs_per_component_on_cell;
800  fe->base_element(base_element_number));
801 
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.);
806 
807  fast_path = true;
808  }
809  else
810  {
811  nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
812  for (unsigned int d = 0; d < n_components; ++d)
813  {
814  const unsigned int component = first_selected_component + d;
815  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
816  {
817  const bool is_primitive =
818  fe->is_primitive() || fe->is_primitive(i);
819  if (is_primitive)
820  nonzero_shape_function_component[i][d] =
821  (component == fe->system_to_component_index(i).first);
822  else
823  nonzero_shape_function_component[i][d] =
824  (fe->get_nonzero_components(i)[component] == true);
825  }
826  }
827 
828  fast_path = false;
829  }
830 }
831 
832 
833 
834 template <int n_components, int dim, int spacedim, typename Number>
835 void
837  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
838  const ArrayView<const Point<dim>> & unit_points)
839 {
840  // reinit is only allowed for mapping computation on the fly
841  AssertThrow(mapping_info_on_the_fly.get() != nullptr, ExcNotImplemented());
842 
843  mapping_info->reinit(cell, unit_points);
844 
845  if (!fast_path)
846  {
847  fe_values = std::make_shared<FEValues<dim, spacedim>>(
848  *mapping,
849  *fe,
851  std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
852  update_flags);
853  fe_values->reinit(cell);
854  }
855 
856  this->unit_points =
857  std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
858 
859  if (update_flags & update_values)
860  values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
861  if (update_flags & update_gradients)
862  gradients.resize(unit_points.size(),
863  numbers::signaling_nan<gradient_type>());
864 }
865 
866 
867 
868 template <int n_components, int dim, int spacedim, typename Number>
869 void
871  const ArrayView<const Number> & solution_values,
872  const EvaluationFlags::EvaluationFlags &evaluation_flag)
873 {
874  const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr;
875  if (precomputed_mapping)
876  {
877  unit_points = mapping_info->get_unit_points();
878 
879  if (update_flags & update_values)
880  values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
881  if (update_flags & update_gradients)
882  gradients.resize(unit_points.size(),
883  numbers::signaling_nan<gradient_type>());
884  }
885 
886  if (unit_points.empty())
887  return;
888 
889  AssertDimension(solution_values.size(), fe->dofs_per_cell);
890  if (((evaluation_flag & EvaluationFlags::values) ||
891  (evaluation_flag & EvaluationFlags::gradients)) &&
892  fast_path)
893  {
894  // fast path with tensor product evaluation
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) *
902  dofs_per_component +
903  i]],
904  comp,
905  solution_renumbered[i]);
906 
907  // unit gradients are currently only implemented with the fast tensor
908  // path
909  unit_gradients.resize(unit_points.size(),
910  numbers::signaling_nan<gradient_type>());
911 
912  const std::size_t n_points = unit_points.size();
913  const std::size_t n_lanes = VectorizedArray<Number>::size();
914  for (unsigned int i = 0; i < n_points; i += n_lanes)
915  {
916  // convert to vectorized format
917  Point<dim, VectorizedArray<Number>> vectorized_points;
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];
921 
922  // compute
923  const auto val_and_grad =
925  poly,
926  solution_renumbered,
927  vectorized_points,
928  polynomials_are_hat_functions);
929 
930  // convert back to standard format
931  if (evaluation_flag & EvaluationFlags::values)
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]);
936  if (evaluation_flag & EvaluationFlags::gradients)
937  {
938  Assert(update_flags & update_gradients ||
939  update_flags & update_inverse_jacobians,
941  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
942  {
944  dim,
945  n_components,
946  Number>::set_gradient(val_and_grad.second,
947  j,
948  unit_gradients[i + j]);
949  gradients[i + j] =
950  apply_transformation(mapping_info->get_mapping_data()
951  .inverse_jacobians[i + j]
952  .transpose(),
953  unit_gradients[i + j]);
954  }
955  }
956  }
957  }
958  else if ((evaluation_flag & EvaluationFlags::values) ||
959  (evaluation_flag & EvaluationFlags::gradients))
960  {
961  // slow path with FEValues
962  Assert(fe_values.get() != nullptr,
963  ExcMessage(
964  "Not initialized. Please call FEPointEvaluation::reinit()!"));
965 
966  if (evaluation_flag & EvaluationFlags::values)
967  {
968  values.resize(unit_points.size());
969  std::fill(values.begin(), values.end(), value_type());
970  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
971  {
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(
984  values[q], d) +=
985  fe_values->shape_value_component(i, q, d) * value;
986  }
987  }
988 
989  if (evaluation_flag & EvaluationFlags::gradients)
990  {
991  gradients.resize(unit_points.size());
992  std::fill(gradients.begin(), gradients.end(), gradient_type());
993  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
994  {
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(
1007  gradients[q], d) +=
1008  fe_values->shape_grad_component(i, q, d) * value;
1009  }
1010  }
1011  }
1012 }
1013 
1014 
1015 
1016 template <int n_components, int dim, int spacedim, typename Number>
1017 void
1019  const ArrayView<Number> & solution_values,
1020  const EvaluationFlags::EvaluationFlags &integration_flags)
1021 {
1022  const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr;
1023  if (precomputed_mapping)
1024  {
1025  unit_points = mapping_info->get_unit_points();
1026 
1027  if (update_flags & update_values)
1028  values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
1029  if (update_flags & update_gradients)
1030  gradients.resize(unit_points.size(),
1031  numbers::signaling_nan<gradient_type>());
1032  }
1033 
1034  if (unit_points.size() == 0) // no evaluation points provided
1035  {
1036  std::fill(solution_values.begin(), solution_values.end(), 0.0);
1037  return;
1038  }
1039 
1040  AssertDimension(solution_values.size(), fe->dofs_per_cell);
1041  if (((integration_flags & EvaluationFlags::values) ||
1042  (integration_flags & EvaluationFlags::gradients)) &&
1043  fast_path)
1044  {
1045  // fast path with tensor product integration
1046 
1047  if (integration_flags & EvaluationFlags::values)
1048  AssertIndexRange(unit_points.size(), values.size() + 1);
1049  if (integration_flags & EvaluationFlags::gradients)
1050  AssertIndexRange(unit_points.size(), gradients.size() + 1);
1051 
1052  if (solution_renumbered_vectorized.size() != dofs_per_component)
1053  solution_renumbered_vectorized.resize(dofs_per_component);
1054  // zero content
1055  solution_renumbered_vectorized.fill(
1057  dim,
1058  n_components,
1060 
1061  const std::size_t n_points = unit_points.size();
1062  const std::size_t n_lanes = VectorizedArray<Number>::size();
1063  for (unsigned int i = 0; i < n_points; i += n_lanes)
1064  {
1065  // convert to vectorized format
1066  Point<dim, VectorizedArray<Number>> vectorized_points;
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];
1070 
1073  value = {};
1074  Tensor<1,
1075  dim,
1077  value_type,
1078  VectorizedArray<Number>>::type>
1079  gradient;
1080 
1081  if (integration_flags & EvaluationFlags::values)
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]);
1086  if (integration_flags & EvaluationFlags::gradients)
1087  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1088  {
1090  mapping_info->get_mapping_data().inverse_jacobians[i + j],
1091  gradients[i + j]);
1094  gradient, j, gradients[i + j]);
1095  }
1096 
1097  // compute
1099  poly,
1100  value,
1101  gradient,
1102  vectorized_points,
1103  solution_renumbered_vectorized);
1104  }
1105 
1106  // add between the lanes and write into the result
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)
1110  {
1111  VectorizedArray<Number> result;
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]] =
1119  result[0];
1120  }
1121  }
1122  else if ((integration_flags & EvaluationFlags::values) ||
1123  (integration_flags & EvaluationFlags::gradients))
1124  {
1125  // slow path with FEValues
1126 
1127  Assert(fe_values.get() != nullptr,
1128  ExcMessage(
1129  "Not initialized. Please call FEPointEvaluation::reinit()!"));
1130  std::fill(solution_values.begin(), solution_values.end(), 0.0);
1131 
1132  if (integration_flags & EvaluationFlags::values)
1133  {
1134  AssertIndexRange(unit_points.size(), values.size() + 1);
1135  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1136  {
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) *
1145  values[q], d);
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) *
1152  values[q], d);
1153  }
1154  }
1155 
1156  if (integration_flags & EvaluationFlags::gradients)
1157  {
1158  AssertIndexRange(unit_points.size(), gradients.size() + 1);
1159  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1160  {
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) *
1169  gradients[q], d);
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) *
1176  gradients[q], d);
1177  }
1178  }
1179  }
1180 }
1181 
1182 
1183 
1184 template <int n_components, int dim, int spacedim, typename Number>
1186  value_type &
1188  const unsigned int point_index) const
1189 {
1190  AssertIndexRange(point_index, values.size());
1191  return values[point_index];
1192 }
1193 
1194 
1195 
1196 template <int n_components, int dim, int spacedim, typename Number>
1198  gradient_type &
1200  const unsigned int point_index) const
1201 {
1202  AssertIndexRange(point_index, gradients.size());
1203  return gradients[point_index];
1204 }
1205 
1206 
1207 
1208 template <int n_components, int dim, int spacedim, typename Number>
1210  gradient_type &
1212  const unsigned int point_index) const
1213 {
1214  Assert(fast_path,
1215  ExcMessage("Unit gradients are currently only implemented for tensor "
1216  "product finite elements combined with MappingQ "
1217  "mappings"));
1218  AssertIndexRange(point_index, unit_gradients.size());
1219  return unit_gradients[point_index];
1220 }
1221 
1222 
1223 
1224 template <int n_components, int dim, int spacedim, typename Number>
1225 inline void
1227  const value_type & value,
1228  const unsigned int point_index)
1229 {
1230  AssertIndexRange(point_index, unit_points.size());
1231  values[point_index] = value;
1232 }
1233 
1234 
1235 
1236 template <int n_components, int dim, int spacedim, typename Number>
1237 inline void
1239  const gradient_type &gradient,
1240  const unsigned int point_index)
1241 {
1242  AssertIndexRange(point_index, unit_points.size());
1243  gradients[point_index] = gradient;
1244 }
1245 
1246 
1247 
1248 template <int n_components, int dim, int spacedim, typename Number>
1251  const unsigned int point_index) const
1252 {
1253  AssertIndexRange(point_index,
1254  mapping_info->get_mapping_data().jacobians.size());
1255  return mapping_info->get_mapping_data().jacobians[point_index];
1256 }
1257 
1258 
1259 
1260 template <int n_components, int dim, int spacedim, typename Number>
1263  const unsigned int point_index) const
1264 {
1265  AssertIndexRange(point_index,
1266  mapping_info->get_mapping_data().inverse_jacobians.size());
1267  return mapping_info->get_mapping_data().inverse_jacobians[point_index];
1268 }
1269 
1270 
1271 
1272 template <int n_components, int dim, int spacedim, typename Number>
1273 inline Point<spacedim>
1275  const unsigned int point_index) const
1276 {
1277  AssertIndexRange(point_index,
1278  mapping_info->get_mapping_data().quadrature_points.size());
1279  return mapping_info->get_mapping_data().quadrature_points[point_index];
1280 }
1281 
1282 
1283 
1284 template <int n_components, int dim, int spacedim, typename Number>
1285 inline Point<dim>
1287  const unsigned int point_index) const
1288 {
1289  AssertIndexRange(point_index, unit_points.size());
1290  return unit_points[point_index];
1291 }
1292 
1294 
1295 #endif
iterator begin() const
Definition: array_view.h:575
iterator end() const
Definition: array_view.h:584
std::size_t size() const
Definition: array_view.h:566
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:504
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1583
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1756
#define AssertIndexRange(index, range)
Definition: exceptions.h:1821
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1672
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:422