Reference documentation for deal.II version 9.3.0
\(\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 - 2021 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 
34 
36 
37 namespace internal
38 {
40  {
45  template <int dim, int n_components, typename Number>
47  {
50 
51  static void
52  read_value(const Number vector_entry,
53  const unsigned int component,
54  value_type & result)
55  {
56  AssertIndexRange(component, n_components);
57  result[component] = vector_entry;
58  }
59 
60  static void
61  write_value(Number & vector_entry,
62  const unsigned int component,
63  const value_type & result)
64  {
65  AssertIndexRange(component, n_components);
66  vector_entry = result[component];
67  }
68 
69  static void
71  const Tensor<1, dim, Tensor<1, n_components, VectorizedArray<Number>>>
72  & value,
73  const unsigned int vector_lane,
74  gradient_type & result)
75  {
76  for (unsigned int i = 0; i < n_components; ++i)
77  for (unsigned int d = 0; d < dim; ++d)
78  result[i][d] = value[d][i][vector_lane];
79  }
80 
81  static void get_gradient(
82  Tensor<1, dim, Tensor<1, n_components, VectorizedArray<Number>>> &value,
83  const unsigned int vector_lane,
84  const gradient_type &result)
85  {
86  for (unsigned int i = 0; i < n_components; ++i)
87  for (unsigned int d = 0; d < dim; ++d)
88  value[d][i][vector_lane] = result[i][d];
89  }
90 
91  static void
92  set_value(const Tensor<1, n_components, VectorizedArray<Number>> &value,
93  const unsigned int vector_lane,
94  value_type & result)
95  {
96  for (unsigned int i = 0; i < n_components; ++i)
97  result[i] = value[i][vector_lane];
98  }
99 
100  static void
101  get_value(Tensor<1, n_components, VectorizedArray<Number>> &value,
102  const unsigned int vector_lane,
103  const value_type & result)
104  {
105  for (unsigned int i = 0; i < n_components; ++i)
106  value[i][vector_lane] = result[i];
107  }
108 
109  template <typename Number2>
110  static Number2 &access(Tensor<1, n_components, Number2> &value,
111  const unsigned int component)
112  {
113  return value[component];
114  }
115 
116  template <typename Number2>
117  static const Number2 &
119  const unsigned int component)
120  {
121  return value[component];
122  }
123  };
124 
125  template <int dim, typename Number>
126  struct EvaluatorTypeTraits<dim, 1, Number>
127  {
128  using value_type = Number;
130 
131  static void
132  read_value(const Number vector_entry,
133  const unsigned int,
134  value_type &result)
135  {
136  result = vector_entry;
137  }
138 
139  static void
140  write_value(Number &vector_entry,
141  const unsigned int,
142  const value_type &result)
143  {
144  vector_entry = result;
145  }
146 
147  static void
149  const unsigned int vector_lane,
150  gradient_type & result)
151  {
152  for (unsigned int d = 0; d < dim; ++d)
153  result[d] = value[d][vector_lane];
154  }
155 
156  static void get_gradient(Tensor<1, dim, VectorizedArray<Number>> &value,
157  const unsigned int vector_lane,
158  const gradient_type &result)
159  {
160  for (unsigned int d = 0; d < dim; ++d)
161  value[d][vector_lane] = result[d];
162  }
163 
164  static void
166  const unsigned int vector_lane,
167  value_type & result)
168  {
169  result = value[vector_lane];
170  }
171 
172  static void
174  const unsigned int vector_lane,
175  const value_type & result)
176  {
177  value[vector_lane] = result;
178  }
179 
180  template <typename Number2>
181  static Number2 &
182  access(Number2 &value, const unsigned int)
183  {
184  return value;
185  }
186 
187  template <typename Number2>
188  static const Number2 &
189  access(const Number2 &value, const unsigned int)
190  {
191  return value;
192  }
193  };
194 
195  template <int dim, typename Number>
196  struct EvaluatorTypeTraits<dim, dim, Number>
197  {
200 
201  static void
202  read_value(const Number vector_entry,
203  const unsigned int component,
204  value_type & result)
205  {
206  result[component] = vector_entry;
207  }
208 
209  static void
210  write_value(Number & vector_entry,
211  const unsigned int component,
212  const value_type & result)
213  {
214  vector_entry = result[component];
215  }
216 
217  static void
219  const Tensor<1, dim, Tensor<1, dim, VectorizedArray<Number>>> &value,
220  const unsigned int vector_lane,
221  gradient_type & result)
222  {
223  for (unsigned int i = 0; i < dim; ++i)
224  for (unsigned int d = 0; d < dim; ++d)
225  result[i][d] = value[d][i][vector_lane];
226  }
227 
228  static void get_gradient(
229  Tensor<1, dim, Tensor<1, dim, VectorizedArray<Number>>> &value,
230  const unsigned int vector_lane,
231  const gradient_type & result)
232  {
233  for (unsigned int i = 0; i < dim; ++i)
234  for (unsigned int d = 0; d < dim; ++d)
235  value[d][i][vector_lane] = result[i][d];
236  }
237 
238  static void
239  set_value(const Tensor<1, dim, VectorizedArray<Number>> &value,
240  const unsigned int vector_lane,
241  value_type & result)
242  {
243  for (unsigned int i = 0; i < dim; ++i)
244  result[i] = value[i][vector_lane];
245  }
246 
247  static void get_value(Tensor<1, dim, VectorizedArray<Number>> &value,
248  const unsigned int vector_lane,
249  const value_type & result)
250  {
251  for (unsigned int i = 0; i < dim; ++i)
252  value[i][vector_lane] = result[i];
253  }
254 
255  static Number &
256  access(value_type &value, const unsigned int component)
257  {
258  return value[component];
259  }
260 
261  static const Number &
262  access(const value_type &value, const unsigned int component)
263  {
264  return value[component];
265  }
266 
267  static Tensor<1, dim> &
268  access(gradient_type &value, const unsigned int component)
269  {
270  return value[component];
271  }
272 
273  static const Tensor<1, dim> &
274  access(const gradient_type &value, const unsigned int component)
275  {
276  return value[component];
277  }
278  };
279 
280  template <typename Number>
281  struct EvaluatorTypeTraits<1, 1, Number>
282  {
283  using value_type = Number;
285 
286  static void
287  read_value(const Number vector_entry,
288  const unsigned int,
289  value_type &result)
290  {
291  result = vector_entry;
292  }
293 
294  static void
295  write_value(Number &vector_entry,
296  const unsigned int,
297  const value_type &result)
298  {
299  vector_entry = result;
300  }
301 
302  static void
304  const unsigned int vector_lane,
305  gradient_type & result)
306  {
307  result[0] = value[0][vector_lane];
308  }
309 
310  static void get_gradient(Tensor<1, 1, VectorizedArray<Number>> &value,
311  const unsigned int vector_lane,
312  const gradient_type &result)
313  {
314  value[0][vector_lane] = result[0];
315  }
316 
317  static void
319  const unsigned int vector_lane,
320  value_type & result)
321  {
322  result = value[vector_lane];
323  }
324 
325  static void
327  const unsigned int vector_lane,
328  const value_type & result)
329  {
330  value[vector_lane] = result;
331  }
332 
333  template <typename Number2>
334  static Number2 &
335  access(Number2 &value, const unsigned int)
336  {
337  return value;
338  }
339 
340  template <typename Number2>
341  static const Number2 &
342  access(const Number2 &value, const unsigned int)
343  {
344  return value;
345  }
346  };
347 
348  template <int dim, int spacedim>
349  bool
351  const unsigned int base_element_number);
352 
353  template <int dim, int spacedim>
354  std::vector<Polynomials::Polynomial<double>>
356  } // namespace FEPointEvaluation
357 } // namespace internal
358 
359 
360 
389 template <int n_components,
390  int dim,
391  int spacedim = dim,
392  typename Number = double>
394 {
395 public:
400 
419  FEPointEvaluation(const Mapping<dim> & mapping,
420  const FiniteElement<dim> &fe,
421  const UpdateFlags update_flags,
422  const unsigned int first_selected_component = 0);
423 
435  void
437  const ArrayView<const Point<dim>> &unit_points);
438 
450  void
451  evaluate(const ArrayView<const Number> & solution_values,
452  const EvaluationFlags::EvaluationFlags &evaluation_flags);
453 
477  void
478  integrate(const ArrayView<Number> & solution_values,
479  const EvaluationFlags::EvaluationFlags &integration_flags);
480 
488  const value_type &
489  get_value(const unsigned int point_index) const;
490 
499  void
500  submit_value(const value_type &value, const unsigned int point_index);
501 
511  const gradient_type &
512  get_gradient(const unsigned int point_index) const;
513 
523  const gradient_type &
524  get_unit_gradient(const unsigned int point_index) const;
525 
534  void
535  submit_gradient(const gradient_type &, const unsigned int point_index);
536 
543  jacobian(const unsigned int point_index) const;
544 
552  inverse_jacobian(const unsigned int point_index) const;
553 
559  real_point(const unsigned int point_index) const;
560 
565  Point<dim>
566  unit_point(const unsigned int point_index) const;
567 
568 private:
573 
579 
584 
589  std::vector<Polynomials::Polynomial<double>> poly;
590 
595 
600  std::vector<unsigned int> renumber;
601 
608  std::vector<value_type> solution_renumbered;
609 
617  dim,
618  n_components,
621 
625  std::vector<value_type> values;
626 
630  std::vector<gradient_type> unit_gradients;
631 
635  std::vector<gradient_type> gradients;
636 
641  unsigned int dofs_per_component;
642 
648  std::vector<std::array<bool, n_components>> nonzero_shape_function_component;
649 
654 
659 
663  std::shared_ptr<FEValues<dim, spacedim>> fe_values;
664 
671 
675  std::vector<Point<dim>> unit_points;
676 };
677 
678 // ----------------------- template and inline function ----------------------
679 
680 
681 template <int n_components, int dim, int spacedim, typename Number>
683  const Mapping<dim> & mapping,
684  const FiniteElement<dim> &fe,
685  const UpdateFlags update_flags,
686  const unsigned int first_selected_component)
687  : mapping(&mapping)
688  , mapping_q_generic(
689  dynamic_cast<const MappingQGeneric<dim, spacedim> *>(&mapping))
690  , fe(&fe)
691  , update_flags(update_flags)
692  , update_flags_mapping(update_default)
693 {
694  bool same_base_element = true;
695  unsigned int base_element_number = 0;
696  unsigned int component = 0;
697  for (; base_element_number < fe.n_base_elements(); ++base_element_number)
698  if (component + fe.element_multiplicity(base_element_number) >
699  first_selected_component)
700  {
701  if (first_selected_component + n_components >
702  component + fe.element_multiplicity(base_element_number))
703  same_base_element = false;
704  break;
705  }
706  else
707  component += fe.element_multiplicity(base_element_number);
708  if (mapping_q_generic != nullptr &&
710  fe, base_element_number) &&
711  same_base_element)
712  {
714 
715  shape_info.reinit(QMidpoint<1>(), fe, base_element_number);
716  renumber = shape_info.lexicographic_numbering;
719  fe.base_element(base_element_number));
720 
722  (poly.size() == 2 && poly[0].value(0.) == 1. &&
723  poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
724  poly[1].value(1.) == 1.);
725  }
726  else
727  {
729  for (unsigned int d = 0; d < n_components; ++d)
730  {
731  const unsigned int component = first_selected_component + d;
732  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
733  {
734  const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
735  if (is_primitive)
737  (component == fe.system_to_component_index(i).first);
738  else
740  (fe.get_nonzero_components(i)[component] == true);
741  }
742  }
743  }
744 
745  // translate update flags
746  if (update_flags & update_jacobians)
748  if (update_flags & update_gradients ||
749  update_flags & update_inverse_jacobians)
751  if (update_flags & update_quadrature_points)
753 }
754 
755 
756 
757 template <int n_components, int dim, int spacedim, typename Number>
758 void
760  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
761  const ArrayView<const Point<dim>> & unit_points)
762 {
763  this->unit_points.resize(unit_points.size());
764  std::copy(unit_points.begin(), unit_points.end(), this->unit_points.begin());
765 
766  if (!poly.empty())
767  mapping_q_generic->fill_mapping_data_for_generic_points(
769  else
770  {
771  fe_values = std::make_shared<FEValues<dim, spacedim>>(
772  *mapping,
773  *fe,
775  std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
777  fe_values->reinit(cell);
778  mapping_data.initialize(unit_points.size(), update_flags_mapping);
780  for (unsigned int q = 0; q < unit_points.size(); ++q)
781  mapping_data.jacobians[q] = fe_values->jacobian(q);
783  for (unsigned int q = 0; q < unit_points.size(); ++q)
784  mapping_data.inverse_jacobians[q] = fe_values->inverse_jacobian(q);
786  for (unsigned int q = 0; q < unit_points.size(); ++q)
787  mapping_data.quadrature_points[q] = fe_values->quadrature_point(q);
788  }
789 
791  values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
793  gradients.resize(unit_points.size(),
794  numbers::signaling_nan<gradient_type>());
795 }
796 
797 
798 
799 template <int n_components, int dim, int spacedim, typename Number>
800 void
802  const ArrayView<const Number> & solution_values,
803  const EvaluationFlags::EvaluationFlags &evaluation_flag)
804 {
805  if (unit_points.empty())
806  return;
807 
808  AssertDimension(solution_values.size(), fe->dofs_per_cell);
809  if (((evaluation_flag & EvaluationFlags::values) ||
810  (evaluation_flag & EvaluationFlags::gradients)) &&
811  !poly.empty())
812  {
813  // fast path with tensor product evaluation
816  for (unsigned int comp = 0; comp < n_components; ++comp)
817  for (unsigned int i = 0; i < dofs_per_component; ++i)
819  EvaluatorTypeTraits<dim, n_components, Number>::read_value(
820  solution_values[renumber[comp * dofs_per_component + i]],
821  comp,
823 
824  // unit gradients are currently only implemented with the fast tensor
825  // path
826  unit_gradients.resize(unit_points.size(),
827  numbers::signaling_nan<gradient_type>());
828 
829  const std::size_t n_points = unit_points.size();
830  const std::size_t n_lanes = VectorizedArray<Number>::size();
831  for (unsigned int i = 0; i < n_points; i += n_lanes)
832  {
833  // convert to vectorized format
834  Point<dim, VectorizedArray<Number>> vectorized_points;
835  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
836  for (unsigned int d = 0; d < dim; ++d)
837  vectorized_points[d][j] = unit_points[i + j][d];
838 
839  // compute
840  const auto val_and_grad =
842  poly,
844  vectorized_points,
846 
847  // convert back to standard format
848  if (evaluation_flag & EvaluationFlags::values)
849  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
851  EvaluatorTypeTraits<dim, n_components, Number>::set_value(
852  val_and_grad.first, j, values[i + j]);
853  if (evaluation_flag & EvaluationFlags::gradients)
854  {
858  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
859  {
860  Assert(update_flags_mapping & update_inverse_jacobians,
863  dim,
864  n_components,
865  Number>::set_gradient(val_and_grad.second,
866  j,
867  unit_gradients[i + j]);
869  mapping_data.inverse_jacobians[i + j].transpose(),
870  unit_gradients[i + j]);
871  }
872  }
873  }
874  }
875  else if ((evaluation_flag & EvaluationFlags::values) ||
876  (evaluation_flag & EvaluationFlags::gradients))
877  {
878  // slow path with FEValues
879  Assert(fe_values.get() != nullptr,
880  ExcMessage(
881  "Not initialized. Please call FEPointEvaluation::reinit()!"));
882 
883  if (evaluation_flag & EvaluationFlags::values)
884  {
885  values.resize(unit_points.size());
886  std::fill(values.begin(), values.end(), value_type());
887  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
888  {
889  const Number value = solution_values[i];
890  for (unsigned int d = 0; d < n_components; ++d)
892  (fe->is_primitive(i) || fe->is_primitive()))
893  for (unsigned int q = 0; q < unit_points.size(); ++q)
895  EvaluatorTypeTraits<dim, n_components, Number>::access(
896  values[q], d) += fe_values->shape_value(i, q) * value;
897  else if (nonzero_shape_function_component[i][d])
898  for (unsigned int q = 0; q < unit_points.size(); ++q)
900  EvaluatorTypeTraits<dim, n_components, Number>::access(
901  values[q], d) +=
902  fe_values->shape_value_component(i, q, d) * value;
903  }
904  }
905 
906  if (evaluation_flag & EvaluationFlags::gradients)
907  {
908  gradients.resize(unit_points.size());
909  std::fill(gradients.begin(), gradients.end(), gradient_type());
910  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
911  {
912  const Number value = solution_values[i];
913  for (unsigned int d = 0; d < n_components; ++d)
915  (fe->is_primitive(i) || fe->is_primitive()))
916  for (unsigned int q = 0; q < unit_points.size(); ++q)
918  EvaluatorTypeTraits<dim, n_components, Number>::access(
919  gradients[q], d) += fe_values->shape_grad(i, q) * value;
920  else if (nonzero_shape_function_component[i][d])
921  for (unsigned int q = 0; q < unit_points.size(); ++q)
923  EvaluatorTypeTraits<dim, n_components, Number>::access(
924  gradients[q], d) +=
925  fe_values->shape_grad_component(i, q, d) * value;
926  }
927  }
928  }
929 }
930 
931 
932 
933 template <int n_components, int dim, int spacedim, typename Number>
934 void
936  const ArrayView<Number> & solution_values,
937  const EvaluationFlags::EvaluationFlags &integration_flags)
938 {
939  if (unit_points.size() == 0) // no evaluation points provided
940  {
941  std::fill(solution_values.begin(), solution_values.end(), 0.0);
942  return;
943  }
944 
945  AssertDimension(solution_values.size(), fe->dofs_per_cell);
946  if (((integration_flags & EvaluationFlags::values) ||
947  (integration_flags & EvaluationFlags::gradients)) &&
948  !poly.empty())
949  {
950  // fast path with tensor product integration
951 
952  if (integration_flags & EvaluationFlags::values)
953  AssertIndexRange(unit_points.size(), values.size() + 1);
954  if (integration_flags & EvaluationFlags::gradients)
955  AssertIndexRange(unit_points.size(), gradients.size() + 1);
956 
959  // zero content
962  dim,
963  n_components,
965 
966  const std::size_t n_points = unit_points.size();
967  const std::size_t n_lanes = VectorizedArray<Number>::size();
968  for (unsigned int i = 0; i < n_points; i += n_lanes)
969  {
970  // convert to vectorized format
971  Point<dim, VectorizedArray<Number>> vectorized_points;
972  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
973  for (unsigned int d = 0; d < dim; ++d)
974  vectorized_points[d][j] = unit_points[i + j][d];
975 
978  value = {};
979  Tensor<1,
980  dim,
982  value_type,
984  gradient;
985 
986  if (integration_flags & EvaluationFlags::values)
987  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
989  EvaluatorTypeTraits<dim, n_components, Number>::get_value(
990  value, j, values[i + j]);
991  if (integration_flags & EvaluationFlags::gradients)
992  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
993  {
996  gradients[i + j] =
997  apply_transformation(mapping_data.inverse_jacobians[i + j],
998  gradients[i + j]);
1001  gradient, j, gradients[i + j]);
1002  }
1003 
1004  // compute
1006  poly,
1007  value,
1008  gradient,
1009  vectorized_points,
1011  }
1012 
1013  // add between the lanes and write into the result
1014  std::fill(solution_values.begin(), solution_values.end(), Number());
1015  for (unsigned int comp = 0; comp < n_components; ++comp)
1016  for (unsigned int i = 0; i < dofs_per_component; ++i)
1017  {
1018  VectorizedArray<Number> result;
1019  internal::FEPointEvaluation::
1020  EvaluatorTypeTraits<dim, n_components, VectorizedArray<Number>>::
1021  write_value(result, comp, solution_renumbered_vectorized[i]);
1022  for (unsigned int lane = n_lanes / 2; lane > 0; lane /= 2)
1023  for (unsigned int j = 0; j < lane; ++j)
1024  result[j] += result[lane + j];
1025  solution_values[renumber[comp * dofs_per_component + i]] =
1026  result[0];
1027  }
1028  }
1029  else if ((integration_flags & EvaluationFlags::values) ||
1030  (integration_flags & EvaluationFlags::gradients))
1031  {
1032  // slow path with FEValues
1033 
1034  Assert(fe_values.get() != nullptr,
1035  ExcMessage(
1036  "Not initialized. Please call FEPointEvaluation::reinit()!"));
1037  std::fill(solution_values.begin(), solution_values.end(), 0.0);
1038 
1039  if (integration_flags & EvaluationFlags::values)
1040  {
1041  AssertIndexRange(unit_points.size(), values.size() + 1);
1042  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1043  {
1044  for (unsigned int d = 0; d < n_components; ++d)
1046  (fe->is_primitive(i) || fe->is_primitive()))
1047  for (unsigned int q = 0; q < unit_points.size(); ++q)
1048  solution_values[i] +=
1049  fe_values->shape_value(i, q) *
1052  values[q], d);
1053  else if (nonzero_shape_function_component[i][d])
1054  for (unsigned int q = 0; q < unit_points.size(); ++q)
1055  solution_values[i] +=
1056  fe_values->shape_value_component(i, q, d) *
1059  values[q], d);
1060  }
1061  }
1062 
1063  if (integration_flags & EvaluationFlags::gradients)
1064  {
1065  AssertIndexRange(unit_points.size(), gradients.size() + 1);
1066  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1067  {
1068  for (unsigned int d = 0; d < n_components; ++d)
1070  (fe->is_primitive(i) || fe->is_primitive()))
1071  for (unsigned int q = 0; q < unit_points.size(); ++q)
1072  solution_values[i] +=
1073  fe_values->shape_grad(i, q) *
1076  gradients[q], d);
1077  else if (nonzero_shape_function_component[i][d])
1078  for (unsigned int q = 0; q < unit_points.size(); ++q)
1079  solution_values[i] +=
1080  fe_values->shape_grad_component(i, q, d) *
1083  gradients[q], d);
1084  }
1085  }
1086  }
1087 }
1088 
1089 
1090 
1091 template <int n_components, int dim, int spacedim, typename Number>
1093  value_type &
1095  const unsigned int point_index) const
1096 {
1097  AssertIndexRange(point_index, values.size());
1098  return values[point_index];
1099 }
1100 
1101 
1102 
1103 template <int n_components, int dim, int spacedim, typename Number>
1105  gradient_type &
1107  const unsigned int point_index) const
1108 {
1109  AssertIndexRange(point_index, gradients.size());
1110  return gradients[point_index];
1111 }
1112 
1113 
1114 
1115 template <int n_components, int dim, int spacedim, typename Number>
1117  gradient_type &
1119  const unsigned int point_index) const
1120 {
1121  Assert(!poly.empty(),
1122  ExcMessage("Unit gradients are currently only implemented for tensor "
1123  "product finite elements combined with MappingQGeneric "
1124  "mappings"));
1125  AssertIndexRange(point_index, unit_gradients.size());
1126  return unit_gradients[point_index];
1127 }
1128 
1129 
1130 
1131 template <int n_components, int dim, int spacedim, typename Number>
1132 inline void
1134  const value_type & value,
1135  const unsigned int point_index)
1136 {
1137  AssertIndexRange(point_index, unit_points.size());
1138  values[point_index] = value;
1139 }
1140 
1141 
1142 
1143 template <int n_components, int dim, int spacedim, typename Number>
1144 inline void
1146  const gradient_type &gradient,
1147  const unsigned int point_index)
1148 {
1149  AssertIndexRange(point_index, unit_points.size());
1150  gradients[point_index] = gradient;
1151 }
1152 
1153 
1154 
1155 template <int n_components, int dim, int spacedim, typename Number>
1158  const unsigned int point_index) const
1159 {
1161  AssertIndexRange(point_index, mapping_data.jacobians.size());
1162  return mapping_data.jacobians[point_index];
1163 }
1164 
1165 
1166 
1167 template <int n_components, int dim, int spacedim, typename Number>
1170  const unsigned int point_index) const
1171 {
1174  ExcNotInitialized());
1175  AssertIndexRange(point_index, mapping_data.inverse_jacobians.size());
1176  return mapping_data.inverse_jacobians[point_index];
1177 }
1178 
1179 
1180 
1181 template <int n_components, int dim, int spacedim, typename Number>
1182 inline Point<spacedim>
1184  const unsigned int point_index) const
1185 {
1187  AssertIndexRange(point_index, mapping_data.quadrature_points.size());
1188  return mapping_data.quadrature_points[point_index];
1189 }
1190 
1191 
1192 
1193 template <int n_components, int dim, int spacedim, typename Number>
1194 inline Point<dim>
1196  const unsigned int point_index) const
1197 {
1198  AssertIndexRange(point_index, unit_points.size());
1199  return unit_points[point_index];
1200 }
1201 
1203 
1204 #endif
AlignedVector< typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, n_components, VectorizedArray< Number > >::value_type > solution_renumbered_vectorized
static void set_gradient(const Tensor< 1, dim, VectorizedArray< Number >> &value, const unsigned int vector_lane, gradient_type &result)
void resize(const size_type new_size)
Shape function values.
static void get_gradient(Tensor< 1, dim, VectorizedArray< Number >> &value, const unsigned int vector_lane, const gradient_type &result)
void submit_gradient(const gradient_type &, const unsigned int point_index)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
static Number2 & access(Number2 &value, const unsigned int)
DerivativeForm< 1, spacedim, dim > inverse_jacobian(const unsigned int point_index) const
std::vector< Polynomials::Polynomial< double > > get_polynomial_space(const FiniteElement< dim, spacedim > &fe)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
std::vector< gradient_type > unit_gradients
void evaluate(const ArrayView< const Number > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
UpdateFlags update_flags_mapping
std::vector< Polynomials::Polynomial< double > > poly
SmartPointer< const Mapping< dim, spacedim > > mapping
Volume element.
const gradient_type & get_unit_gradient(const unsigned int point_index) const
iterator end() const
Definition: array_view.h:592
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
SmartPointer< const FiniteElement< dim > > fe
void submit_value(const value_type &value, const unsigned int point_index)
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< gradient_type > gradients
static Number2 & access(Number2 &value, const unsigned int)
static void read_value(const Number vector_entry, const unsigned int component, value_type &result)
static void get_gradient(Tensor< 1, dim, Tensor< 1, dim, VectorizedArray< Number >>> &value, const unsigned int vector_lane, const gradient_type &result)
void reinit(const Quadrature< dim_q > &quad, const FiniteElement< dim > &fe_dim, const unsigned int base_element=0)
const MappingQGeneric< dim, spacedim > * mapping_q_generic
static const Number2 & access(const Number2 &value, const unsigned int)
Transformed quadrature points.
static ::ExceptionBase & ExcNotInitialized()
static void get_gradient(Tensor< 1, dim, Tensor< 1, n_components, VectorizedArray< Number >>> &value, const unsigned int vector_lane, const gradient_type &result)
bool is_primitive() const
Definition: fe.h:3302
internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > mapping_data
static void set_gradient(const Tensor< 1, 1, VectorizedArray< Number >> &value, const unsigned int vector_lane, gradient_type &result)
std::vector< Point< dim > > unit_points
std::size_t size() const
Definition: array_view.h:574
static ::ExceptionBase & ExcMessage(std::string arg1)
static Number & access(value_type &value, const unsigned int component)
No update.
static const Tensor< 1, dim > & access(const gradient_type &value, const unsigned int component)
Point< dim > unit_point(const unsigned int point_index) const
bool is_fast_path_supported(const FiniteElement< dim, spacedim > &fe, const unsigned int base_element_number)
#define Assert(cond, exc)
Definition: exceptions.h:1465
unsigned int element_multiplicity(const unsigned int index) const
Definition: fe.h:3123
UpdateFlags
static const Number2 & access(const Number2 &value, const unsigned int)
static Tensor< 1, dim > & access(gradient_type &value, const unsigned int component)
static void set_value(const Tensor< 1, dim, VectorizedArray< Number >> &value, const unsigned int vector_lane, value_type &result)
const ComponentMask & get_nonzero_components(const unsigned int i) const
Definition: fe.h:3282
Tensor< 1, n_components, Number > value_type
static Number2 & access(Tensor< 1, n_components, Number2 > &value, const unsigned int component)
static void set_value(const Tensor< 1, n_components, VectorizedArray< Number >> &value, const unsigned int vector_lane, value_type &result)
static void write_value(Number &vector_entry, const unsigned int component, const value_type &result)
static void write_value(Number &vector_entry, const unsigned int component, const value_type &result)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
const value_type & get_value(const unsigned int point_index) const
static const Number & access(const value_type &value, const unsigned int component)
std::vector< std::array< bool, n_components > > nonzero_shape_function_component
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
Definition: fe.cc:1402
Tensor< 1, n_components, Tensor< 1, dim, Number > > gradient_type
static void get_value(Tensor< 1, dim, VectorizedArray< Number >> &value, const unsigned int vector_lane, const 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)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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={})
const gradient_type & get_gradient(const unsigned int point_index) const
static void read_value(const Number vector_entry, const unsigned int, value_type &result)
const unsigned int dofs_per_cell
Definition: fe_base.h:419
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, n_components, Number >::value_type value_type
std::vector< value_type > solution_renumbered
static const Number2 & access(const Tensor< 1, n_components, Number2 > &value, const unsigned int component)
static void write_value(Number &vector_entry, const unsigned int, const value_type &result)
void integrate(const ArrayView< Number > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags)
static void set_gradient(const Tensor< 1, dim, Tensor< 1, n_components, VectorizedArray< Number >>> &value, const unsigned int vector_lane, gradient_type &result)
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:3100
unsigned int n_dofs_per_cell() const
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={})
Definition: tensor.h:462
DerivativeForm< 1, dim, spacedim > jacobian(const unsigned int point_index) const
static void get_gradient(Tensor< 1, 1, VectorizedArray< Number >> &value, const unsigned int vector_lane, const gradient_type &result)
void reinit(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< dim >> &unit_points)
static void read_value(const Number vector_entry, const unsigned int component, value_type &result)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
Shape function gradients.
Point< spacedim > real_point(const unsigned int point_index) const
static void read_value(const Number vector_entry, const unsigned int, value_type &result)
size_type size() const
std::vector< value_type > values
static void set_value(const VectorizedArray< Number > &value, const unsigned int vector_lane, value_type &result)
std::shared_ptr< FEValues< dim, spacedim > > fe_values
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, n_components, Number >::gradient_type gradient_type
static void get_value(VectorizedArray< Number > &value, const unsigned int vector_lane, const value_type &result)
static void write_value(Number &vector_entry, const unsigned int, const value_type &result)
iterator begin() const
Definition: array_view.h:583
std::vector< unsigned int > lexicographic_numbering
Definition: shape_info.h:380
static void get_value(VectorizedArray< Number > &value, const unsigned int vector_lane, const value_type &result)
static void get_value(Tensor< 1, n_components, VectorizedArray< Number >> &value, const unsigned int vector_lane, const value_type &result)
unsigned int dofs_per_component
unsigned int n_base_elements() const
Definition: fe.h:3114
void copy(const T *begin, const T *end, U *dest)
EvaluationFlags
The EvaluationFlags enum.
std::vector< unsigned int > renumber
static void set_value(const VectorizedArray< Number > &value, const unsigned int vector_lane, value_type &result)
FEPointEvaluation(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const UpdateFlags update_flags, const unsigned int first_selected_component=0)