Reference documentation for deal.II version GIT 7b2de2f2f9 2023-09-24 11: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 - 2023 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 
36 
38 
40 
41 namespace internal
42 {
44  {
47  std::string,
48  << "You are requesting information from an FEPointEvaluation "
49  << "object for which this kind of information has not been computed. "
50  << "What information these objects compute is determined by the update_* "
51  << "flags you pass to MappingInfo() in the Constructor. Here, "
52  << "the operation you are attempting requires the <" << arg1
53  << "> flag to be set, but it was apparently not specified "
54  << "upon initialization.");
55 
60  template <int dim, int n_components, typename Number>
62  {
63  using ScalarNumber =
66  typename ::internal::VectorizedArrayTrait<
79 
80  static void
81  read_value(const ScalarNumber vector_entry,
82  const unsigned int component,
83  scalar_value_type &result)
84  {
85  AssertIndexRange(component, n_components);
86  result[component] = vector_entry;
87  }
88 
89  static ScalarNumber
90  sum_value(const unsigned int component,
91  const vectorized_value_type &result)
92  {
93  AssertIndexRange(component, n_components);
94  return result[component].sum();
95  }
96 
97  static void
99  const unsigned int vector_lane,
100  gradient_type &result)
101  {
102  for (unsigned int i = 0; i < n_components; ++i)
103  for (unsigned int d = 0; d < dim; ++d)
104  result[i][d] =
106  value[d][i], vector_lane);
107  }
108 
109  static void
111  const unsigned int vector_lane,
112  const gradient_type &result)
113  {
114  for (unsigned int i = 0; i < n_components; ++i)
115  for (unsigned int d = 0; d < dim; ++d)
117  value[d][i], vector_lane) = result[i][d];
118  }
119 
120  static void
121  set_zero_gradient(gradient_type &value, const unsigned int vector_lane)
122  {
123  for (unsigned int i = 0; i < n_components; ++i)
124  for (unsigned int d = 0; d < dim; ++d)
126  vector_lane) = 0.;
127  }
128 
129  static void
131  const unsigned int vector_lane,
132  scalar_value_type &result)
133  {
134  for (unsigned int i = 0; i < n_components; ++i)
135  result[i] = value[i][vector_lane];
136  }
137 
138  static void
140  const unsigned int,
141  vectorized_value_type &result)
142  {
143  result = value;
144  }
145 
146  static void
148  const unsigned int vector_lane,
149  const scalar_value_type &result)
150  {
151  for (unsigned int i = 0; i < n_components; ++i)
152  value[i][vector_lane] = result[i];
153  }
154 
155  static void
157  const unsigned int,
158  const vectorized_value_type &result)
159  {
160  value = result;
161  }
162 
163  static void
164  set_zero_value(value_type &value, const unsigned int vector_lane)
165  {
166  for (unsigned int i = 0; i < n_components; ++i)
168  0.;
169  }
170 
171  static void
173  const unsigned int vector_lane,
174  const unsigned int component,
175  const ScalarNumber &shape_value)
176  {
178  vector_lane) += shape_value;
179  }
180 
181  static ScalarNumber
182  access(const value_type &value,
183  const unsigned int vector_lane,
184  const unsigned int component)
185  {
187  vector_lane);
188  }
189 
190  static void
192  const unsigned int vector_lane,
193  const unsigned int component,
194  const Tensor<1, dim, ScalarNumber> &shape_gradient)
195  {
196  for (unsigned int d = 0; d < dim; ++d)
198  vector_lane) +=
199  shape_gradient[d];
200  }
201 
203  access(const gradient_type &value,
204  const unsigned int vector_lane,
205  const unsigned int component)
206  {
208  for (unsigned int d = 0; d < dim; ++d)
209  result[d] =
211  vector_lane);
212  return result;
213  }
214  };
215 
216  template <int dim, typename Number>
217  struct EvaluatorTypeTraits<dim, 1, Number>
218  {
219  using ScalarNumber =
222  typename ::internal::VectorizedArrayTrait<
224  using value_type = Number;
231 
232  static void
233  read_value(const ScalarNumber vector_entry,
234  const unsigned int,
235  scalar_value_type &result)
236  {
237  result = vector_entry;
238  }
239 
240  static ScalarNumber
241  sum_value(const unsigned int, const vectorized_value_type &result)
242  {
243  return result.sum();
244  }
245 
246  static void
248  const unsigned int vector_lane,
249  scalar_gradient_type &result)
250  {
251  for (unsigned int d = 0; d < dim; ++d)
252  result[d] = value[d][vector_lane];
253  }
254 
255  static void
257  const unsigned int,
258  vectorized_gradient_type &result)
259  {
260  result = value;
261  }
262 
263  static void
265  const unsigned int vector_lane,
266  const scalar_gradient_type &result)
267  {
268  for (unsigned int d = 0; d < dim; ++d)
269  value[d][vector_lane] = result[d];
270  }
271 
272  static void
274  const unsigned int,
275  const vectorized_gradient_type &result)
276  {
277  value = result;
278  }
279 
280  static void
281  set_zero_gradient(gradient_type &value, const unsigned int vector_lane)
282  {
283  for (unsigned int d = 0; d < dim; ++d)
285  0.;
286  }
287 
288  static void
290  const unsigned int vector_lane,
291  scalar_value_type &result)
292  {
293  result = value[vector_lane];
294  }
295 
296  static void
298  const unsigned int,
299  vectorized_value_type &result)
300  {
301  result = value;
302  }
303 
304  static void
306  const unsigned int vector_lane,
307  const scalar_value_type &result)
308  {
309  value[vector_lane] = result;
310  }
311 
312  static void
314  const unsigned int,
315  const vectorized_value_type &result)
316  {
317  value = result;
318  }
319 
320  static void
321  set_zero_value(value_type &value, const unsigned int vector_lane)
322  {
324  }
325 
326  static void
328  const unsigned int vector_lane,
329  const unsigned int,
330  const ScalarNumber &shape_value)
331  {
333  shape_value;
334  }
335 
336  static ScalarNumber
337  access(const value_type &value,
338  const unsigned int vector_lane,
339  const unsigned int)
340  {
342  }
343 
344  static void
346  const unsigned int vector_lane,
347  const unsigned int,
348  const scalar_gradient_type &shape_gradient)
349  {
350  for (unsigned int d = 0; d < dim; ++d)
352  shape_gradient[d];
353  }
354 
355  static scalar_gradient_type
356  access(const gradient_type &value,
357  const unsigned int vector_lane,
358  const unsigned int)
359  {
360  scalar_gradient_type result;
361  for (unsigned int d = 0; d < dim; ++d)
362  result[d] =
364  return result;
365  }
366  };
367 
368  template <int dim, typename Number>
369  struct EvaluatorTypeTraits<dim, dim, Number>
370  {
371  using ScalarNumber =
374  typename ::internal::VectorizedArrayTrait<
384 
385  static void
386  read_value(const ScalarNumber vector_entry,
387  const unsigned int component,
388  scalar_value_type &result)
389  {
390  AssertIndexRange(component, dim);
391  result[component] = vector_entry;
392  }
393 
394  static ScalarNumber
395  sum_value(const unsigned int component,
396  const vectorized_value_type &result)
397  {
398  AssertIndexRange(component, dim);
399  return result[component].sum();
400  }
401 
402  static void
404  const unsigned int vector_lane,
405  gradient_type &result)
406  {
407  for (unsigned int i = 0; i < dim; ++i)
408  for (unsigned int d = 0; d < dim; ++d)
409  result[i][d] =
411  value[d][i], vector_lane);
412  }
413 
414  static void
416  const unsigned int vector_lane,
417  const gradient_type &result)
418  {
419  for (unsigned int i = 0; i < dim; ++i)
420  for (unsigned int d = 0; d < dim; ++d)
422  value[d][i], vector_lane) = result[i][d];
423  }
424 
425  static void
426  set_zero_gradient(gradient_type &value, const unsigned int vector_lane)
427  {
428  for (unsigned int i = 0; i < dim; ++i)
429  for (unsigned int d = 0; d < dim; ++d)
431  vector_lane) = 0.;
432  }
433 
434  static void
436  const unsigned int vector_lane,
437  scalar_value_type &result)
438  {
439  for (unsigned int i = 0; i < dim; ++i)
440  result[i] = value[i][vector_lane];
441  }
442 
443  static void
445  const unsigned int,
446  vectorized_value_type &result)
447  {
448  result = value;
449  }
450 
451  static void
453  const unsigned int vector_lane,
454  const scalar_value_type &result)
455  {
456  for (unsigned int i = 0; i < dim; ++i)
457  value[i][vector_lane] = result[i];
458  }
459 
460  static void
462  const unsigned int,
463  const vectorized_value_type &result)
464  {
465  value = result;
466  }
467 
468  static void
469  set_zero_value(value_type &value, const unsigned int vector_lane)
470  {
471  for (unsigned int i = 0; i < dim; ++i)
473  0.;
474  }
475 
476  static void
478  const unsigned int vector_lane,
479  const unsigned int component,
480  const ScalarNumber &shape_value)
481  {
483  vector_lane) += shape_value;
484  }
485 
486  static ScalarNumber
487  access(const value_type &value,
488  const unsigned int vector_lane,
489  const unsigned int component)
490  {
492  vector_lane);
493  }
494 
495  static void
497  const unsigned int vector_lane,
498  const unsigned int component,
499  const Tensor<1, dim, ScalarNumber> &shape_gradient)
500  {
501  for (unsigned int d = 0; d < dim; ++d)
503  vector_lane) +=
504  shape_gradient[d];
505  }
506 
508  access(const gradient_type &value,
509  const unsigned int vector_lane,
510  const unsigned int component)
511  {
513  for (unsigned int d = 0; d < dim; ++d)
514  result[d] =
516  vector_lane);
517  return result;
518  }
519  };
520 
521  template <typename Number>
522  struct EvaluatorTypeTraits<1, 1, Number>
523  {
524  using ScalarNumber =
527  typename ::internal::VectorizedArrayTrait<
529  using value_type = Number;
536 
537  static void
538  read_value(const ScalarNumber vector_entry,
539  const unsigned int,
540  scalar_value_type &result)
541  {
542  result = vector_entry;
543  }
544 
545  static ScalarNumber
546  sum_value(const unsigned int, const vectorized_value_type &result)
547  {
548  return result.sum();
549  }
550 
551  static void
553  const unsigned int vector_lane,
554  scalar_gradient_type &result)
555  {
556  result[0] = value[0][vector_lane];
557  }
558 
559  static void
561  const unsigned int,
562  vectorized_gradient_type &result)
563  {
564  result = value;
565  }
566 
567  static void
569  const unsigned int vector_lane,
570  const scalar_gradient_type &result)
571  {
572  value[0][vector_lane] = result[0];
573  }
574 
575  static void
577  const unsigned int,
578  const vectorized_gradient_type &result)
579  {
580  value = result;
581  }
582 
583  static void
584  set_zero_gradient(gradient_type &value, const unsigned int vector_lane)
585  {
587  }
588 
589  static void
591  const unsigned int vector_lane,
592  scalar_value_type &result)
593  {
594  result = value[vector_lane];
595  }
596 
597  static void
599  const unsigned int,
600  vectorized_value_type &result)
601  {
602  result = value;
603  }
604 
605  static void
607  const unsigned int vector_lane,
608  const scalar_value_type &result)
609  {
610  value[vector_lane] = result;
611  }
612 
613  static void
615  const unsigned int,
616  const vectorized_value_type &result)
617  {
618  value = result;
619  }
620 
621  static void
622  set_zero_value(value_type &value, const unsigned int vector_lane)
623  {
625  }
626 
627  static void
629  const unsigned int vector_lane,
630  const unsigned int,
631  const ScalarNumber &shape_value)
632  {
634  shape_value;
635  }
636 
637  static ScalarNumber
638  access(const value_type &value,
639  const unsigned int vector_lane,
640  const unsigned int)
641  {
643  }
644 
645  static void
647  const unsigned int vector_lane,
648  const unsigned int,
649  const scalar_gradient_type &shape_gradient)
650  {
652  shape_gradient[0];
653  }
654 
655  static scalar_gradient_type
656  access(const gradient_type &value,
657  const unsigned int vector_lane,
658  const unsigned int)
659  {
660  scalar_gradient_type result;
661  result[0] =
663  return result;
664  }
665  };
666 
667  template <int dim, int spacedim>
668  bool
670  const unsigned int base_element_number);
671 
672  template <int dim, int spacedim>
673  bool
675 
676  template <int dim, int spacedim>
677  std::vector<Polynomials::Polynomial<double>>
679  } // namespace FEPointEvaluation
680 } // namespace internal
681 
682 
683 
714 template <int n_components_,
715  int dim,
716  int spacedim = dim,
717  typename Number = double>
719 {
720 public:
721  static constexpr unsigned int dimension = dim;
722  static constexpr unsigned int n_components = n_components_;
723 
724  using number_type = Number;
725 
726  using ScalarNumber =
728  using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
730  using ETT = typename internal::FEPointEvaluation::
731  EvaluatorTypeTraits<dim, n_components, Number>;
732  using value_type = typename ETT::value_type;
733  using scalar_value_type = typename ETT::scalar_value_type;
734  using vectorized_value_type = typename ETT::vectorized_value_type;
735  using gradient_type = typename ETT::gradient_type;
737  typename ETT::interface_vectorized_gradient_type;
738 
758  const FiniteElement<dim> &fe,
760  const unsigned int first_selected_component = 0);
761 
780  const FiniteElement<dim> &fe,
781  const unsigned int first_selected_component = 0);
782 
786  FEPointEvaluation(FEPointEvaluation &other) noexcept;
787 
791  FEPointEvaluation(FEPointEvaluation &&other) noexcept;
792 
797 
809  void
811  const ArrayView<const Point<dim>> &unit_points);
812 
817  void
818  reinit();
819 
824  void
825  reinit(const unsigned int cell_index);
826 
831  void
832  reinit(const unsigned int cell_index, const unsigned int face_number);
833 
845  template <std::size_t stride_view>
846  void
847  evaluate(
849  const EvaluationFlags::EvaluationFlags &evaluation_flags);
850 
862  void
863  evaluate(const ArrayView<const ScalarNumber> &solution_values,
864  const EvaluationFlags::EvaluationFlags &evaluation_flags)
865  {
867  solution_values.size()),
868  evaluation_flags);
869  }
870 
889  template <std::size_t stride_view>
890  void
892  const EvaluationFlags::EvaluationFlags &integration_flags);
893 
894 
913  void
914  integrate(const ArrayView<ScalarNumber> &solution_values,
915  const EvaluationFlags::EvaluationFlags &integration_flags)
916  {
918  solution_values.size()),
919  integration_flags);
920  }
921 
944  template <std::size_t stride_view>
945  void
946  test_and_sum(
947  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
948  const EvaluationFlags::EvaluationFlags &integration_flags);
949 
972  void
973  test_and_sum(const ArrayView<ScalarNumber> &solution_values,
974  const EvaluationFlags::EvaluationFlags &integration_flags)
975  {
977  solution_values.size()),
978  integration_flags);
979  }
980 
988  const value_type &
989  get_value(const unsigned int point_index) const;
990 
999  void
1000  submit_value(const value_type &value, const unsigned int point_index);
1001 
1011  const gradient_type &
1012  get_gradient(const unsigned int point_index) const;
1013 
1022  void
1023  submit_gradient(const gradient_type &, const unsigned int point_index);
1024 
1031  jacobian(const unsigned int point_index) const;
1032 
1040  inverse_jacobian(const unsigned int point_index) const;
1041 
1047  Number
1048  JxW(const unsigned int point_index) const;
1049 
1056  normal_vector(const unsigned int point_index) const;
1057 
1063  real_point(const unsigned int point_index) const;
1064 
1070  unit_point(const unsigned int point_index) const;
1071 
1078  quadrature_point_indices() const;
1079 
1080 private:
1081  static constexpr std::size_t n_lanes_user_interface =
1083  static constexpr std::size_t n_lanes_internal =
1085  static constexpr std::size_t stride =
1087 
1096  void
1097  setup(const unsigned int first_selected_component);
1098 
1104  template <bool is_face, bool is_linear>
1105  void
1106  do_reinit();
1107 
1112  template <bool is_face_path, bool is_linear, std::size_t stride_view>
1113  void
1116  const EvaluationFlags::EvaluationFlags &evaluation_flags);
1117 
1122  template <bool is_face_path, bool is_linear, std::size_t stride_view>
1123  void
1126  const EvaluationFlags::EvaluationFlags &evaluation_flags,
1127  const unsigned int n_shapes,
1128  const unsigned int qb,
1129  vectorized_value_type &value,
1131 
1135  template <bool is_face_path, bool is_linear, std::size_t stride_view>
1136  void
1137  evaluate_fast(
1139  const EvaluationFlags::EvaluationFlags &evaluation_flags);
1140 
1144  template <std::size_t stride_view>
1145  void
1146  evaluate_slow(
1148  const EvaluationFlags::EvaluationFlags &evaluation_flags);
1149 
1155  template <bool is_face_path, bool is_linear>
1156  void
1158  const EvaluationFlags::EvaluationFlags &integration_flags,
1159  const unsigned int n_shapes,
1160  const unsigned int qb,
1161  const vectorized_value_type value,
1162  const interface_vectorized_gradient_type gradient,
1163  vectorized_value_type *solution_values_vectorized_linear);
1164 
1170  template <bool is_face_path, bool is_linear, std::size_t stride_view>
1171  void
1173  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1174  const EvaluationFlags::EvaluationFlags &integration_flags,
1175  vectorized_value_type *solution_values_vectorized_linear);
1176 
1180  template <bool do_JxW,
1181  bool is_face_path,
1182  bool is_linear,
1183  std::size_t stride_view>
1184  void
1186  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1187  const EvaluationFlags::EvaluationFlags &integration_flags);
1188 
1192  template <bool do_JxW, std::size_t stride_view>
1193  void
1195  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1196  const EvaluationFlags::EvaluationFlags &integration_flags);
1197 
1201  template <bool do_JxW, std::size_t stride_view>
1202  void
1203  do_integrate(
1204  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1205  const EvaluationFlags::EvaluationFlags &integration_flags);
1206 
1210  const unsigned int n_q_batches;
1211 
1215  const unsigned int n_q_points;
1216 
1220  const unsigned int n_q_points_scalar;
1221 
1226 
1231 
1236  std::vector<Polynomials::Polynomial<double>> poly;
1237 
1242 
1247  std::vector<unsigned int> renumber;
1248 
1256  std::vector<scalar_value_type> solution_renumbered;
1257 
1265 
1270 
1274  std::vector<value_type> values;
1275 
1279  std::vector<gradient_type> gradients;
1280 
1286 
1292 
1298 
1304 
1310 
1316 
1322  const Number *JxW_ptr;
1323 
1328 
1333  unsigned int dofs_per_component;
1334 
1341 
1347 
1352 
1357 
1363  std::vector<std::array<bool, n_components>> nonzero_shape_function_component;
1364 
1369 
1373  std::shared_ptr<FEValues<dim, spacedim>> fe_values;
1374 
1378  std::unique_ptr<NonMatching::MappingInfo<dim, spacedim, Number>>
1380 
1386 
1390  unsigned int current_cell_index;
1391 
1395  unsigned int current_face_number;
1396 
1401 
1406  boost::signals2::connection connection_is_reinitialized;
1407 
1412 
1418 
1424 };
1425 
1426 // ----------------------- template and inline function ----------------------
1427 
1428 
1429 template <int n_components_, int dim, int spacedim, typename Number>
1431  const Mapping<dim> &mapping,
1432  const FiniteElement<dim> &fe,
1433  const UpdateFlags update_flags,
1434  const unsigned int first_selected_component)
1435  : n_q_batches(numbers::invalid_unsigned_int)
1436  , n_q_points(numbers::invalid_unsigned_int)
1437  , n_q_points_scalar(numbers::invalid_unsigned_int)
1438  , mapping(&mapping)
1439  , fe(&fe)
1440  , JxW_ptr(nullptr)
1441  , use_face_path(false)
1442  , update_flags(update_flags)
1443  , mapping_info_on_the_fly(
1444  std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
1445  mapping,
1446  update_flags))
1447  , mapping_info(mapping_info_on_the_fly.get())
1448  , current_cell_index(numbers::invalid_unsigned_int)
1449  , current_face_number(numbers::invalid_unsigned_int)
1450  , is_reinitialized(false)
1451 {
1452  setup(first_selected_component);
1453 }
1454 
1455 
1456 
1457 template <int n_components_, int dim, int spacedim, typename Number>
1460  const FiniteElement<dim> &fe,
1461  const unsigned int first_selected_component)
1462  : n_q_batches(numbers::invalid_unsigned_int)
1463  , n_q_points(numbers::invalid_unsigned_int)
1464  , n_q_points_scalar(numbers::invalid_unsigned_int)
1465  , mapping(&mapping_info.get_mapping())
1466  , fe(&fe)
1467  , JxW_ptr(nullptr)
1468  , use_face_path(false)
1469  , update_flags(mapping_info.get_update_flags())
1470  , mapping_info(&mapping_info)
1471  , current_cell_index(numbers::invalid_unsigned_int)
1472  , current_face_number(numbers::invalid_unsigned_int)
1473  , is_reinitialized(false)
1474 {
1475  setup(first_selected_component);
1476  connection_is_reinitialized = mapping_info.connect_is_reinitialized(
1477  [this]() { this->is_reinitialized = false; });
1478 }
1479 
1480 
1481 
1482 template <int n_components_, int dim, int spacedim, typename Number>
1485  : n_q_batches(other.n_q_batches)
1486  , n_q_points(other.n_q_points)
1487  , n_q_points_scalar(other.n_q_points_scalar)
1488  , mapping(other.mapping)
1489  , fe(other.fe)
1490  , poly(other.poly)
1491  , use_linear_path(other.use_linear_path)
1492  , renumber(other.renumber)
1493  , solution_renumbered(other.solution_renumbered)
1494  , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1495  , values(other.values)
1496  , gradients(other.gradients)
1497  , dofs_per_component(other.dofs_per_component)
1498  , dofs_per_component_face(other.dofs_per_component_face)
1499  , use_face_path(false)
1500  , component_in_base_element(other.component_in_base_element)
1501  , nonzero_shape_function_component(other.nonzero_shape_function_component)
1502  , update_flags(other.update_flags)
1503  , fe_values(other.fe_values)
1504  , mapping_info_on_the_fly(
1505  other.mapping_info_on_the_fly ?
1507  *mapping,
1508  update_flags) :
1509  nullptr)
1510  , mapping_info(other.mapping_info)
1511  , current_cell_index(other.current_cell_index)
1512  , current_face_number(other.current_face_number)
1513  , fast_path(other.fast_path)
1514  , is_reinitialized(false)
1515  , shapes(other.shapes)
1516  , shapes_faces(other.shapes_faces)
1517 {
1518  connection_is_reinitialized = mapping_info->connect_is_reinitialized(
1519  [this]() { this->is_reinitialized = false; });
1520 }
1521 
1522 
1523 
1524 template <int n_components_, int dim, int spacedim, typename Number>
1527  : n_q_batches(other.n_q_batches)
1528  , n_q_points(other.n_q_points)
1529  , n_q_points_scalar(other.n_q_points_scalar)
1530  , mapping(other.mapping)
1531  , fe(other.fe)
1532  , poly(other.poly)
1533  , use_linear_path(other.use_linear_path)
1534  , renumber(other.renumber)
1535  , solution_renumbered(other.solution_renumbered)
1536  , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1537  , values(other.values)
1538  , gradients(other.gradients)
1539  , dofs_per_component(other.dofs_per_component)
1540  , dofs_per_component_face(other.dofs_per_component_face)
1541  , use_face_path(false)
1542  , component_in_base_element(other.component_in_base_element)
1543  , nonzero_shape_function_component(other.nonzero_shape_function_component)
1544  , update_flags(other.update_flags)
1545  , fe_values(other.fe_values)
1546  , mapping_info_on_the_fly(std::move(other.mapping_info_on_the_fly))
1547  , mapping_info(other.mapping_info)
1548  , current_cell_index(other.current_cell_index)
1549  , current_face_number(other.current_face_number)
1550  , fast_path(other.fast_path)
1551  , is_reinitialized(false)
1552  , shapes(other.shapes)
1553  , shapes_faces(other.shapes_faces)
1554 {
1555  connection_is_reinitialized = mapping_info->connect_is_reinitialized(
1556  [this]() { this->is_reinitialized = false; });
1557 }
1558 
1559 
1560 
1561 template <int n_components_, int dim, int spacedim, typename Number>
1563 {
1564  connection_is_reinitialized.disconnect();
1565 }
1566 
1567 
1568 
1569 template <int n_components_, int dim, int spacedim, typename Number>
1570 void
1572  const unsigned int first_selected_component)
1573 {
1574  AssertIndexRange(first_selected_component + n_components,
1575  fe->n_components() + 1);
1576 
1577  shapes.reserve(100);
1578 
1579  bool same_base_element = true;
1580  unsigned int base_element_number = 0;
1581  component_in_base_element = 0;
1582  unsigned int component = 0;
1583  for (; base_element_number < fe->n_base_elements(); ++base_element_number)
1584  if (component + fe->element_multiplicity(base_element_number) >
1585  first_selected_component)
1586  {
1587  if (first_selected_component + n_components >
1588  component + fe->element_multiplicity(base_element_number))
1589  same_base_element = false;
1590  component_in_base_element = first_selected_component - component;
1591  break;
1592  }
1593  else
1594  component += fe->element_multiplicity(base_element_number);
1595 
1598  *fe, base_element_number) &&
1599  same_base_element)
1600  {
1601  shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
1602  renumber = shape_info.lexicographic_numbering;
1603  dofs_per_component = shape_info.dofs_per_component_on_cell;
1604  dofs_per_component_face = shape_info.dofs_per_component_on_face;
1606  fe->base_element(base_element_number));
1607 
1608  bool is_lexicographic = true;
1609  for (unsigned int i = 0; i < renumber.size(); ++i)
1610  if (i != renumber[i])
1611  is_lexicographic = false;
1612 
1613  if (is_lexicographic)
1614  renumber.clear();
1615 
1616  use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. &&
1617  poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1618  poly[1].value(1.) == 1.) &&
1619  (fe->n_components() == n_components);
1620 
1621  const unsigned int size_face = 2 * dofs_per_component_face;
1622  const unsigned int size_cell = dofs_per_component;
1623  scratch_data_scalar.resize(size_face + size_cell);
1624 
1625  solution_renumbered.resize(dofs_per_component);
1626  solution_renumbered_vectorized.resize(dofs_per_component);
1627 
1628  fast_path = true;
1629  }
1630  else
1631  {
1632  nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
1633  for (unsigned int d = 0; d < n_components; ++d)
1634  {
1635  const unsigned int component = first_selected_component + d;
1636  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1637  {
1638  const bool is_primitive =
1639  fe->is_primitive() || fe->is_primitive(i);
1640  if (is_primitive)
1641  nonzero_shape_function_component[i][d] =
1642  (component == fe->system_to_component_index(i).first);
1643  else
1644  nonzero_shape_function_component[i][d] =
1645  (fe->get_nonzero_components(i)[component] == true);
1646  }
1647  }
1648 
1649  fast_path = false;
1650  }
1651 }
1652 
1653 
1654 
1655 template <int n_components_, int dim, int spacedim, typename Number>
1656 inline void
1658  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1659  const ArrayView<const Point<dim>> &unit_points)
1660 {
1661  // reinit is only allowed for mapping computation on the fly
1662  AssertThrow(mapping_info_on_the_fly.get() != nullptr, ExcNotImplemented());
1663 
1664  mapping_info->reinit(cell, unit_points);
1665 
1666  if (!fast_path)
1667  {
1668  fe_values = std::make_shared<FEValues<dim, spacedim>>(
1669  *mapping,
1670  *fe,
1672  std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
1673  update_flags);
1674  fe_values->reinit(cell);
1675  }
1676 
1677  if (use_linear_path)
1678  do_reinit<false, true>();
1679  else
1680  do_reinit<false, false>();
1681 }
1682 
1683 
1684 
1685 template <int n_components_, int dim, int spacedim, typename Number>
1686 inline void
1688 {
1689  current_cell_index = numbers::invalid_unsigned_int;
1690  current_face_number = numbers::invalid_unsigned_int;
1691 
1692  if (use_linear_path)
1693  do_reinit<false, true>();
1694  else
1695  do_reinit<false, false>();
1696 }
1697 
1698 
1699 
1700 template <int n_components_, int dim, int spacedim, typename Number>
1701 inline void
1703  const unsigned int cell_index)
1704 {
1705  current_cell_index = cell_index;
1706  current_face_number = numbers::invalid_unsigned_int;
1707 
1708  if (use_linear_path)
1709  do_reinit<false, true>();
1710  else
1711  do_reinit<false, false>();
1712 
1713  if (!fast_path)
1714  {
1715  std::vector<Point<dim>> unit_points(n_q_points_scalar);
1716 
1717  for (unsigned int v = 0; v < n_q_points_scalar; ++v)
1718  for (unsigned int d = 0; d < dim; ++d)
1719  unit_points[v][d] =
1720  unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
1721 
1722  fe_values = std::make_shared<FEValues<dim, spacedim>>(
1723  *mapping,
1724  *fe,
1726  std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
1727  update_flags);
1728 
1729  fe_values->reinit(mapping_info->get_cell_iterator(current_cell_index));
1730  }
1731 }
1732 
1733 
1734 
1735 template <int n_components_, int dim, int spacedim, typename Number>
1736 inline void
1738  const unsigned int cell_index,
1739  const unsigned int face_number)
1740 {
1741  current_cell_index = cell_index;
1742  current_face_number = face_number;
1743 
1744  if (use_linear_path)
1745  do_reinit<true, true>();
1746  else
1747  do_reinit<true, false>();
1748 }
1749 
1750 
1751 
1752 template <int n_components_, int dim, int spacedim, typename Number>
1753 template <bool is_face, bool is_linear>
1754 inline void
1756 {
1757  const unsigned int geometry_index =
1758  mapping_info->template compute_geometry_index_offset<is_face>(
1759  current_cell_index, current_face_number);
1760 
1761  cell_type = mapping_info->get_cell_type(geometry_index);
1762 
1763  const_cast<unsigned int &>(n_q_points_scalar) =
1764  mapping_info->get_n_q_points_unvectorized(geometry_index);
1765 
1766  // round up n_q_points_scalar / n_lanes_internal
1767  const_cast<unsigned int &>(n_q_batches) =
1768  (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
1769 
1770  const unsigned int n_q_points_before = n_q_points;
1771 
1772  const_cast<unsigned int &>(n_q_points) =
1773  (stride == 1) ? n_q_batches : n_q_points_scalar;
1774 
1775  if (n_q_points != n_q_points_before)
1776  {
1777  if (update_flags & update_values)
1778  values.resize(n_q_points);
1779  if (update_flags & update_gradients)
1780  gradients.resize(n_q_points);
1781  }
1782 
1783  if (n_q_points == 0)
1784  {
1785  is_reinitialized = true;
1786  return;
1787  }
1788 
1789  // use face path if mapping_info in face state and number of quadrature points
1790  // is large enough
1791  use_face_path =
1792  is_face && (mapping_info->is_face_state() && n_q_points_scalar >= 6);
1793 
1794  // set unit point pointer
1795  const unsigned int unit_point_offset =
1796  mapping_info->compute_unit_point_index_offset(geometry_index);
1797 
1798  if (use_face_path)
1799  unit_point_faces_ptr =
1800  mapping_info->get_unit_point_faces(unit_point_offset);
1801  else
1802  unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
1803 
1804  // set data pointers
1805  const unsigned int data_offset =
1806  mapping_info->compute_data_index_offset(geometry_index);
1807  const unsigned int compressed_data_offset =
1808  mapping_info->compute_compressed_data_index_offset(geometry_index);
1809 #ifdef DEBUG
1810  const UpdateFlags update_flags_mapping =
1811  mapping_info->get_update_flags_mapping();
1812  if (update_flags_mapping & UpdateFlags::update_quadrature_points)
1813  real_point_ptr = mapping_info->get_real_point(data_offset);
1814  if (update_flags_mapping & UpdateFlags::update_jacobians)
1815  jacobian_ptr = mapping_info->get_jacobian(compressed_data_offset);
1816  if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
1817  inverse_jacobian_ptr =
1818  mapping_info->get_inverse_jacobian(compressed_data_offset);
1819  if (is_face && update_flags_mapping & UpdateFlags::update_normal_vectors)
1820  normal_ptr = mapping_info->get_normal_vector(data_offset);
1821  if (update_flags_mapping & UpdateFlags::update_JxW_values)
1822  JxW_ptr = mapping_info->get_JxW(data_offset);
1823 #else
1824  real_point_ptr = mapping_info->get_real_point(data_offset);
1825  jacobian_ptr = mapping_info->get_jacobian(compressed_data_offset);
1826  inverse_jacobian_ptr =
1827  mapping_info->get_inverse_jacobian(compressed_data_offset);
1828  normal_ptr = mapping_info->get_normal_vector(data_offset);
1829  JxW_ptr = mapping_info->get_JxW(data_offset);
1830 #endif
1831 
1832  if (!is_linear && fast_path)
1833  {
1834  const std::size_t n_shapes = poly.size();
1835 
1836  for (unsigned int qb = 0; qb < n_q_batches; ++qb)
1837  if (use_face_path)
1838  {
1839  if (dim > 1)
1840  {
1841  shapes_faces.resize_fast(n_q_batches * n_shapes);
1843  shapes_faces.data() + qb * n_shapes,
1844  poly,
1845  unit_point_faces_ptr[qb],
1846  update_flags & UpdateFlags::update_gradients ? 1 : 0);
1847  }
1848  }
1849  else
1850  {
1851  shapes.resize_fast(n_q_batches * n_shapes);
1853  shapes.data() + qb * n_shapes,
1854  poly,
1855  unit_point_ptr[qb],
1856  update_flags & UpdateFlags::update_gradients ? 1 : 0);
1857  }
1858  }
1859 
1860  is_reinitialized = true;
1861 }
1862 
1863 
1864 
1865 template <int n_components_, int dim, int spacedim, typename Number>
1866 template <bool is_face_path, bool is_linear, std::size_t stride_view>
1867 inline void
1870  const EvaluationFlags::EvaluationFlags &evaluation_flags)
1871 {
1872  const unsigned int dofs_per_comp =
1873  is_linear ? Utilities::pow(2, dim) : dofs_per_component;
1874 
1875  for (unsigned int comp = 0; comp < n_components; ++comp)
1876  {
1877  const std::size_t offset =
1878  (component_in_base_element + comp) * dofs_per_comp;
1879 
1880  if (is_face_path)
1881  {
1882  const ScalarNumber *input;
1883  if (is_linear || renumber.empty())
1884  {
1885  for (unsigned int i = 0; i < dofs_per_comp; ++i)
1886  scratch_data_scalar[i] = solution_values[i + offset];
1887  input = scratch_data_scalar.data();
1888  }
1889  else
1890  {
1891  const unsigned int *renumber_ptr = renumber.data() + offset;
1892  for (unsigned int i = 0; i < dofs_per_comp; ++i)
1893  scratch_data_scalar[i] = solution_values[renumber_ptr[i]];
1894  input = scratch_data_scalar.data();
1895  }
1896 
1897  ScalarNumber *output = scratch_data_scalar.begin() + dofs_per_comp;
1898 
1900  template interpolate<true, false>(1,
1901  evaluation_flags,
1902  shape_info,
1903  input,
1904  output,
1905  current_face_number);
1906 
1907  const unsigned int dofs_per_comp_face =
1908  is_linear ? Utilities::pow(2, dim - 1) : dofs_per_component_face;
1909  for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
1910  ETT::read_value(output[i], comp, solution_renumbered[i]);
1911  }
1912  else
1913  {
1914  if ((is_linear && n_components == 1) || renumber.empty())
1915  {
1916  for (unsigned int i = 0; i < dofs_per_comp; ++i)
1917  ETT::read_value(solution_values[i + offset],
1918  comp,
1919  solution_renumbered[i]);
1920  }
1921  else
1922  {
1923  const unsigned int *renumber_ptr = renumber.data() + offset;
1924  for (unsigned int i = 0; i < dofs_per_comp; ++i)
1925  ETT::read_value(solution_values[renumber_ptr[i]],
1926  comp,
1927  solution_renumbered[i]);
1928  }
1929  }
1930  }
1931 }
1932 
1933 
1934 
1935 template <int n_components_, int dim, int spacedim, typename Number>
1936 template <bool is_face_path, bool is_linear, std::size_t stride_view>
1937 inline void
1940  const EvaluationFlags::EvaluationFlags &evaluation_flags,
1941  const unsigned int n_shapes,
1942  const unsigned int qb,
1943  vectorized_value_type &value,
1945 {
1946  if (is_face_path)
1947  {
1948  if (evaluation_flags & EvaluationFlags::gradients)
1949  {
1950  const std::array<vectorized_value_type, dim + 1> interpolated_value =
1951  is_linear ?
1953  dim - 1,
1956  2>(solution_renumbered.data(), unit_point_faces_ptr[qb]) :
1958  dim - 1,
1961  2,
1962  false>(shapes_faces.data() + qb * n_shapes,
1963  n_shapes,
1964  solution_renumbered.data());
1965 
1966  value = interpolated_value[dim - 1];
1967  // reorder derivative from tangential/normal derivatives into tensor
1968  // in physical coordinates
1969  if (current_face_number / 2 == 0)
1970  {
1971  gradient[0] = interpolated_value[dim];
1972  if (dim > 1)
1973  gradient[1] = interpolated_value[0];
1974  if (dim > 2)
1975  gradient[2] = interpolated_value[1];
1976  }
1977  else if (current_face_number / 2 == 1)
1978  {
1979  if (dim > 1)
1980  gradient[1] = interpolated_value[dim];
1981  if (dim == 3)
1982  {
1983  gradient[0] = interpolated_value[1];
1984  gradient[2] = interpolated_value[0];
1985  }
1986  else if (dim == 2)
1987  gradient[0] = interpolated_value[0];
1988  else
1989  Assert(false, ExcInternalError());
1990  }
1991  else if (current_face_number / 2 == 2)
1992  {
1993  if (dim > 2)
1994  {
1995  gradient[0] = interpolated_value[0];
1996  gradient[1] = interpolated_value[1];
1997  gradient[2] = interpolated_value[dim];
1998  }
1999  else
2000  Assert(false, ExcInternalError());
2001  }
2002  else
2003  Assert(false, ExcInternalError());
2004  }
2005  else
2006  {
2008  dim - 1,
2010  VectorizedArrayType>(solution_renumbered.data(),
2011  unit_point_faces_ptr[qb]) :
2013  dim - 1,
2016  false>(shapes_faces.data() + qb * n_shapes,
2017  n_shapes,
2018  solution_renumbered.data());
2019  }
2020  }
2021  else
2022  {
2023  if (evaluation_flags & EvaluationFlags::gradients)
2024  {
2025  std::array<vectorized_value_type, dim + 1> result;
2026  if constexpr (is_linear)
2027  {
2028  if constexpr (n_components == 1)
2029  result =
2031  dim,
2034  1,
2035  stride_view>(solution_values.data(), unit_point_ptr[qb]);
2036  else
2037  result =
2039  solution_renumbered.data(), unit_point_ptr[qb]);
2040  }
2041  else
2042  result =
2044  dim,
2047  1,
2048  false>(shapes.data() + qb * n_shapes,
2049  n_shapes,
2050  solution_renumbered.data());
2051  gradient[0] = result[0];
2052  if (dim > 1)
2053  gradient[1] = result[1];
2054  if (dim > 2)
2055  gradient[2] = result[2];
2056  value = result[dim];
2057  }
2058  else
2059  {
2060  if constexpr (is_linear)
2061  {
2062  if constexpr (n_components == 1)
2064  dim,
2067  stride_view>(solution_values.data(), unit_point_ptr[qb]);
2068  else
2070  solution_renumbered.data(), unit_point_ptr[qb]);
2071  }
2072  else
2074  dim,
2077  false>(shapes.data() + qb * n_shapes,
2078  n_shapes,
2079  solution_renumbered.data());
2080  }
2081  }
2082 }
2083 
2084 
2085 
2086 template <int n_components_, int dim, int spacedim, typename Number>
2087 template <bool is_face_path, bool is_linear, std::size_t stride_view>
2088 inline void
2091  const EvaluationFlags::EvaluationFlags &evaluation_flags)
2092 {
2093  if (!(is_linear && n_components == 1) || is_face_path)
2094  prepare_evaluate_fast<is_face_path, is_linear>(solution_values,
2095  evaluation_flags);
2096 
2097  // loop over quadrature batches qb
2098  const unsigned int n_shapes = is_linear ? 2 : poly.size();
2099 
2100  for (unsigned int qb = 0; qb < n_q_batches; ++qb)
2101  {
2102  vectorized_value_type value;
2104 
2105  compute_evaluate_fast<is_face_path, is_linear>(
2106  solution_values, evaluation_flags, n_shapes, qb, value, gradient);
2107 
2108  if (evaluation_flags & EvaluationFlags::values)
2109  {
2110  for (unsigned int v = 0, offset = qb * stride;
2111  v < stride && (stride == 1 || offset < n_q_points_scalar);
2112  ++v, ++offset)
2113  ETT::set_value(value, v, values[offset]);
2114  }
2115  if (evaluation_flags & EvaluationFlags::gradients)
2116  {
2117  Assert(update_flags & update_gradients ||
2118  update_flags & update_inverse_jacobians,
2119  ExcNotInitialized());
2120 
2121  for (unsigned int v = 0, offset = qb * stride;
2122  v < stride && (stride == 1 || offset < n_q_points_scalar);
2123  ++v, ++offset)
2124  {
2125  gradient_type unit_gradient;
2126  ETT::set_gradient(gradient, v, unit_gradient);
2127  gradients[offset] =
2128  cell_type <=
2130  apply_diagonal_transformation(inverse_jacobian_ptr[0],
2131  unit_gradient) :
2133  inverse_jacobian_ptr
2134  [cell_type <=
2136  0 :
2137  offset]
2138  .transpose(),
2139  unit_gradient);
2140  }
2141  }
2142  }
2143 }
2144 
2145 
2146 
2147 template <int n_components_, int dim, int spacedim, typename Number>
2148 template <std::size_t stride_view>
2149 inline void
2152  const EvaluationFlags::EvaluationFlags &evaluation_flags)
2153 {
2154  // slow path with FEValues
2155  Assert(fe_values.get() != nullptr,
2156  ExcMessage(
2157  "Not initialized. Please call FEPointEvaluation::reinit()!"));
2158 
2159  const std::size_t n_points = fe_values->get_quadrature().size();
2160 
2161  if (evaluation_flags & EvaluationFlags::values)
2162  {
2163  values.resize(n_q_points);
2164  std::fill(values.begin(), values.end(), value_type());
2165  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
2166  {
2167  const ScalarNumber value = solution_values[i];
2168  for (unsigned int d = 0; d < n_components; ++d)
2169  if (nonzero_shape_function_component[i][d] &&
2170  (fe->is_primitive(i) || fe->is_primitive()))
2171  for (unsigned int qb = 0, q = 0; q < n_points;
2172  ++qb, q += n_lanes_user_interface)
2173  for (unsigned int v = 0;
2174  v < n_lanes_user_interface && q + v < n_points;
2175  ++v)
2176  ETT::access(values[qb],
2177  v,
2178  d,
2179  fe_values->shape_value(i, q + v) * value);
2180  else if (nonzero_shape_function_component[i][d])
2181  for (unsigned int qb = 0, q = 0; q < n_points;
2182  ++qb, q += n_lanes_user_interface)
2183  for (unsigned int v = 0;
2184  v < n_lanes_user_interface && q + v < n_points;
2185  ++v)
2186  ETT::access(values[qb],
2187  v,
2188  d,
2189  fe_values->shape_value_component(i, q + v, d) *
2190  value);
2191  }
2192  }
2193 
2194  if (evaluation_flags & EvaluationFlags::gradients)
2195  {
2196  gradients.resize(n_q_points);
2197  std::fill(gradients.begin(), gradients.end(), gradient_type());
2198  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
2199  {
2200  const ScalarNumber value = solution_values[i];
2201  for (unsigned int d = 0; d < n_components; ++d)
2202  if (nonzero_shape_function_component[i][d] &&
2203  (fe->is_primitive(i) || fe->is_primitive()))
2204  for (unsigned int qb = 0, q = 0; q < n_points;
2205  ++qb, q += n_lanes_user_interface)
2206  for (unsigned int v = 0;
2207  v < n_lanes_user_interface && q + v < n_points;
2208  ++v)
2209  ETT::access(gradients[qb],
2210  v,
2211  d,
2212  fe_values->shape_grad(i, q + v) * value);
2213  else if (nonzero_shape_function_component[i][d])
2214  for (unsigned int qb = 0, q = 0; q < n_points;
2215  ++qb, q += n_lanes_user_interface)
2216  for (unsigned int v = 0;
2217  v < n_lanes_user_interface && q + v < n_points;
2218  ++v)
2219  ETT::access(gradients[qb],
2220  v,
2221  d,
2222  fe_values->shape_grad_component(i, q + v, d) *
2223  value);
2224  }
2225  }
2226 }
2227 
2228 
2229 
2230 template <int n_components_, int dim, int spacedim, typename Number>
2231 template <std::size_t stride_view>
2232 void
2235  const EvaluationFlags::EvaluationFlags &evaluation_flags)
2236 {
2237  if (!is_reinitialized)
2238  reinit();
2239 
2240  if (n_q_points == 0)
2241  return;
2242 
2243  Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
2244 
2245  if (!((evaluation_flags & EvaluationFlags::values) ||
2246  (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
2247  return;
2248 
2249  AssertDimension(solution_values.size(), fe->dofs_per_cell);
2250  if (fast_path)
2251  {
2252  if (use_face_path)
2253  {
2254  if (use_linear_path)
2255  evaluate_fast<true, true>(solution_values, evaluation_flags);
2256  else
2257  evaluate_fast<true, false>(solution_values, evaluation_flags);
2258  }
2259  else
2260  {
2261  if (use_linear_path)
2262  evaluate_fast<false, true>(solution_values, evaluation_flags);
2263  else
2264  evaluate_fast<false, false>(solution_values, evaluation_flags);
2265  }
2266  }
2267  else
2268  evaluate_slow(solution_values, evaluation_flags);
2269 }
2270 
2271 
2272 
2273 template <int n_components_, int dim, int spacedim, typename Number>
2274 template <bool is_face_path, bool is_linear>
2275 inline void
2277  const EvaluationFlags::EvaluationFlags &integration_flags,
2278  const unsigned int n_shapes,
2279  const unsigned int qb,
2280  const vectorized_value_type value,
2281  const interface_vectorized_gradient_type gradient,
2282  vectorized_value_type *solution_values_vectorized_linear)
2283 {
2284  if (is_face_path)
2285  {
2286  if (integration_flags & EvaluationFlags::gradients)
2287  {
2288  std::array<vectorized_value_type, 2> value_face = {};
2289  Tensor<1, dim - 1, vectorized_value_type> gradient_in_face;
2290 
2291  value_face[0] = value;
2292  // fill derivative in physical coordinates into tangential/normal
2293  // derivatives
2294  if (current_face_number / 2 == 0)
2295  {
2296  value_face[1] = gradient[0];
2297  if (dim > 1)
2298  gradient_in_face[0] = gradient[1];
2299  if (dim > 2)
2300  gradient_in_face[1] = gradient[2];
2301  }
2302  else if (current_face_number / 2 == 1)
2303  {
2304  if (dim > 1)
2305  value_face[1] = gradient[1];
2306  if (dim == 3)
2307  {
2308  gradient_in_face[0] = gradient[2];
2309  gradient_in_face[1] = gradient[0];
2310  }
2311  else if (dim == 2)
2312  gradient_in_face[0] = gradient[0];
2313  else
2314  Assert(false, ExcInternalError());
2315  }
2316  else if (current_face_number / 2 == 2)
2317  {
2318  if (dim > 2)
2319  {
2320  value_face[1] = gradient[2];
2321  gradient_in_face[0] = gradient[0];
2322  gradient_in_face[1] = gradient[1];
2323  }
2324  else
2325  Assert(false, ExcInternalError());
2326  }
2327  else
2328  Assert(false, ExcInternalError());
2329 
2331  is_linear,
2332  dim - 1,
2335  2>(shapes_faces.data() + qb * n_shapes,
2336  n_shapes,
2337  value_face.data(),
2338  gradient_in_face,
2339  is_linear ? solution_values_vectorized_linear :
2340  solution_renumbered_vectorized.data(),
2341  unit_point_faces_ptr[qb],
2342  qb != 0);
2343  }
2344  else
2346  dim - 1,
2349  shapes_faces.data() + qb * n_shapes,
2350  n_shapes,
2351  value,
2352  is_linear ? solution_values_vectorized_linear :
2353  solution_renumbered_vectorized.data(),
2354  unit_point_faces_ptr[qb],
2355  qb != 0);
2356  }
2357  else
2358  {
2359  if (integration_flags & EvaluationFlags::gradients)
2361  is_linear,
2362  dim,
2364  vectorized_value_type>(shapes.data() + qb * n_shapes,
2365  n_shapes,
2366  &value,
2367  gradient,
2368  is_linear ?
2369  solution_values_vectorized_linear :
2370  solution_renumbered_vectorized.data(),
2371  unit_point_ptr[qb],
2372  qb != 0);
2373  else
2375  dim,
2378  shapes.data() + qb * n_shapes,
2379  n_shapes,
2380  value,
2381  is_linear ? solution_values_vectorized_linear :
2382  solution_renumbered_vectorized.data(),
2383  unit_point_ptr[qb],
2384  qb != 0);
2385  }
2386 }
2387 
2388 
2389 
2390 template <int n_components_, int dim, int spacedim, typename Number>
2391 template <bool is_face_path, bool is_linear, std::size_t stride_view>
2392 inline void
2394  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2395  const EvaluationFlags::EvaluationFlags &integration_flags,
2396  vectorized_value_type *solution_values_vectorized_linear)
2397 {
2398  if (fe->n_components() > n_components)
2399  for (unsigned int i = 0; i < solution_values.size(); ++i)
2400  solution_values[i] = 0;
2401 
2402  const unsigned int dofs_per_comp =
2403  is_linear ? Utilities::pow(2, dim) : dofs_per_component;
2404 
2405  for (unsigned int comp = 0; comp < n_components; ++comp)
2406  {
2407  const std::size_t offset =
2408  (component_in_base_element + comp) * dofs_per_comp;
2409 
2410  if (is_face_path)
2411  {
2412  const unsigned int dofs_per_comp_face =
2413  is_linear ? Utilities::pow(2, dim - 1) : dofs_per_component_face;
2414 
2415  const unsigned int size_input = 2 * dofs_per_comp_face;
2416  ScalarNumber *input = scratch_data_scalar.begin();
2417  ScalarNumber *output = input + size_input;
2418 
2419  for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
2420  input[i] =
2421  ETT::sum_value(comp,
2422  is_linear ?
2423  *(solution_values_vectorized_linear + i) :
2424  solution_renumbered_vectorized[i]);
2425 
2427  template interpolate<false, false>(1,
2428  integration_flags,
2429  shape_info,
2430  input,
2431  output,
2432  current_face_number);
2433 
2434  if (is_linear || renumber.empty())
2435  {
2436  for (unsigned int i = 0; i < dofs_per_comp; ++i)
2437  solution_values[i + offset] = output[i];
2438  }
2439  else
2440  {
2441  const unsigned int *renumber_ptr = renumber.data() + offset;
2442  for (unsigned int i = 0; i < dofs_per_comp; ++i)
2443  solution_values[renumber_ptr[i]] = output[i];
2444  }
2445  }
2446  else
2447  {
2448  if (is_linear || renumber.empty())
2449  {
2450  for (unsigned int i = 0; i < dofs_per_comp; ++i)
2451  solution_values[i + offset] =
2452  ETT::sum_value(comp,
2453  is_linear ?
2454  *(solution_values_vectorized_linear + i) :
2455  solution_renumbered_vectorized[i]);
2456  }
2457  else
2458  {
2459  const unsigned int *renumber_ptr = renumber.data() + offset;
2460  for (unsigned int i = 0; i < dofs_per_comp; ++i)
2461  solution_values[renumber_ptr[i]] =
2462  ETT::sum_value(comp, solution_renumbered_vectorized[i]);
2463  }
2464  }
2465  }
2466 }
2467 
2468 
2469 
2470 template <int n_components_, int dim, int spacedim, typename Number>
2471 template <bool do_JxW,
2472  bool is_face_path,
2473  bool is_linear,
2474  std::size_t stride_view>
2475 inline void
2477  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2478  const EvaluationFlags::EvaluationFlags &integration_flags)
2479 {
2480  // zero out lanes of incomplete last quadrature point batch
2481  if constexpr (stride == 1)
2482  if (const unsigned int n_filled_lanes =
2483  n_q_points_scalar & (n_lanes_internal - 1);
2484  n_filled_lanes > 0)
2485  {
2486  if (integration_flags & EvaluationFlags::values)
2487  for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2488  ETT::set_zero_value(values.back(), v);
2489  if (integration_flags & EvaluationFlags::gradients)
2490  for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2491  ETT::set_zero_gradient(gradients.back(), v);
2492  }
2493 
2494  std::array<vectorized_value_type,
2495  is_linear ? (is_face_path ? 2 * Utilities::pow(2, dim - 1) :
2496  Utilities::pow(2, dim)) :
2497  0>
2498  solution_values_vectorized_linear = {};
2499 
2500  // loop over quadrature batches qb
2501  const unsigned int n_shapes = is_linear ? 2 : poly.size();
2502 
2503  const bool cartesian_cell =
2505  const bool affine_cell =
2507  for (unsigned int qb = 0; qb < n_q_batches; ++qb)
2508  {
2509  vectorized_value_type value = {};
2511 
2512  if (integration_flags & EvaluationFlags::values)
2513  for (unsigned int v = 0, offset = qb * stride;
2514  v < stride && (stride == 1 || offset < n_q_points_scalar);
2515  ++v, ++offset)
2516  ETT::get_value(value,
2517  v,
2518  do_JxW ? values[offset] * JxW_ptr[offset] :
2519  values[offset]);
2520 
2521  if (integration_flags & EvaluationFlags::gradients)
2522  for (unsigned int v = 0, offset = qb * stride;
2523  v < stride && (stride == 1 || offset < n_q_points_scalar);
2524  ++v, ++offset)
2525  {
2526  const auto grad_w =
2527  do_JxW ? gradients[offset] * JxW_ptr[offset] : gradients[offset];
2528  ETT::get_gradient(
2529  gradient,
2530  v,
2531  cartesian_cell ?
2532  apply_diagonal_transformation(inverse_jacobian_ptr[0], grad_w) :
2534  inverse_jacobian_ptr[affine_cell ? 0 : offset], grad_w));
2535  }
2536 
2537  compute_integrate_fast<is_face_path, is_linear>(
2538  integration_flags,
2539  n_shapes,
2540  qb,
2541  value,
2542  gradient,
2543  solution_values_vectorized_linear.data());
2544  }
2545 
2546  // add between the lanes and write into the result
2547  finish_integrate_fast<is_face_path, is_linear>(
2548  solution_values,
2549  integration_flags,
2550  solution_values_vectorized_linear.data());
2551 }
2552 
2553 
2554 
2555 template <int n_components_, int dim, int spacedim, typename Number>
2556 template <bool do_JxW, std::size_t stride_view>
2557 inline void
2559  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2560  const EvaluationFlags::EvaluationFlags &integration_flags)
2561 {
2562  // slow path with FEValues
2563  Assert(fe_values.get() != nullptr,
2564  ExcMessage(
2565  "Not initialized. Please call FEPointEvaluation::reinit()!"));
2566  for (unsigned int i = 0; i < solution_values.size(); ++i)
2567  solution_values[i] = 0;
2568 
2569  const std::size_t n_points = fe_values->get_quadrature().size();
2570 
2571  if (integration_flags & EvaluationFlags::values)
2572  {
2573  AssertIndexRange(n_q_points, values.size() + 1);
2574  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
2575  {
2576  for (unsigned int d = 0; d < n_components; ++d)
2577  if (nonzero_shape_function_component[i][d] &&
2578  (fe->is_primitive(i) || fe->is_primitive()))
2579  for (unsigned int qb = 0, q = 0; q < n_points;
2580  ++qb, q += n_lanes_user_interface)
2581  for (unsigned int v = 0;
2582  v < n_lanes_user_interface && q + v < n_points;
2583  ++v)
2584  solution_values[i] += fe_values->shape_value(i, q + v) *
2585  ETT::access(values[qb], v, d) *
2586  (do_JxW ? fe_values->JxW(q + v) : 1.);
2587  else if (nonzero_shape_function_component[i][d])
2588  for (unsigned int qb = 0, q = 0; q < n_points;
2589  ++qb, q += n_lanes_user_interface)
2590  for (unsigned int v = 0;
2591  v < n_lanes_user_interface && q + v < n_points;
2592  ++v)
2593  solution_values[i] +=
2594  fe_values->shape_value_component(i, q + v, d) *
2595  ETT::access(values[qb], v, d) *
2596  (do_JxW ? fe_values->JxW(q + v) : 1.);
2597  }
2598  }
2599 
2600  if (integration_flags & EvaluationFlags::gradients)
2601  {
2602  AssertIndexRange(n_q_points, gradients.size() + 1);
2603  for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
2604  {
2605  for (unsigned int d = 0; d < n_components; ++d)
2606  if (nonzero_shape_function_component[i][d] &&
2607  (fe->is_primitive(i) || fe->is_primitive()))
2608  for (unsigned int qb = 0, q = 0; q < n_points;
2609  ++qb, q += n_lanes_user_interface)
2610  for (unsigned int v = 0;
2611  v < n_lanes_user_interface && q + v < n_points;
2612  ++v)
2613  solution_values[i] += fe_values->shape_grad(i, q + v) *
2614  ETT::access(gradients[qb], v, d) *
2615  (do_JxW ? fe_values->JxW(q + v) : 1.);
2616  else if (nonzero_shape_function_component[i][d])
2617  for (unsigned int qb = 0, q = 0; q < n_points;
2618  ++qb, q += n_lanes_user_interface)
2619  for (unsigned int v = 0;
2620  v < n_lanes_user_interface && q + v < n_points;
2621  ++v)
2622  solution_values[i] +=
2623  fe_values->shape_grad_component(i, q + v, d) *
2624  ETT::access(gradients[qb], v, d) *
2625  (do_JxW ? fe_values->JxW(q + v) : 1.);
2626  }
2627  }
2628 }
2629 
2630 
2631 
2632 template <int n_components_, int dim, int spacedim, typename Number>
2633 template <bool do_JxW, std::size_t stride_view>
2634 void
2636  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2637  const EvaluationFlags::EvaluationFlags &integration_flags)
2638 {
2639  if (!is_reinitialized)
2640  reinit();
2641 
2642  Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
2643 
2644  if (n_q_points == 0 || // no evaluation points provided
2645  !((integration_flags & EvaluationFlags::values) ||
2646  (integration_flags &
2647  EvaluationFlags::gradients))) // no integration flags
2648  {
2649  for (unsigned int i = 0; i < solution_values.size(); ++i)
2650  solution_values[i] = 0;
2651  return;
2652  }
2653 
2654  Assert(
2655  !do_JxW || JxW_ptr != nullptr,
2656  ExcMessage(
2657  "JxW pointer is not set! If you do not want to integrate() use test_and_sum()"));
2658 
2659  AssertDimension(solution_values.size(), fe->dofs_per_cell);
2660  if (fast_path)
2661  {
2662  if (use_face_path)
2663  {
2664  if (use_linear_path)
2665  integrate_fast<do_JxW, true, true>(solution_values,
2666  integration_flags);
2667  else
2668  integrate_fast<do_JxW, true, false>(solution_values,
2669  integration_flags);
2670  }
2671  else
2672  {
2673  if (use_linear_path)
2674  integrate_fast<do_JxW, false, true>(solution_values,
2675  integration_flags);
2676  else
2677  integrate_fast<do_JxW, false, false>(solution_values,
2678  integration_flags);
2679  }
2680  }
2681  else
2682  integrate_slow<do_JxW>(solution_values, integration_flags);
2683 }
2684 
2685 
2686 
2687 template <int n_components_, int dim, int spacedim, typename Number>
2688 template <std::size_t stride_view>
2689 void
2691  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2692  const EvaluationFlags::EvaluationFlags &integration_flags)
2693 {
2694  do_integrate<true>(solution_values, integration_flags);
2695 }
2696 
2697 
2698 
2699 template <int n_components_, int dim, int spacedim, typename Number>
2700 template <std::size_t stride_view>
2701 void
2703  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2704  const EvaluationFlags::EvaluationFlags &integration_flags)
2705 {
2706  do_integrate<false>(solution_values, integration_flags);
2707 }
2708 
2709 
2710 
2711 template <int n_components_, int dim, int spacedim, typename Number>
2713  value_type &
2715  const unsigned int point_index) const
2716 {
2717  AssertIndexRange(point_index, values.size());
2718  return values[point_index];
2719 }
2720 
2721 
2722 
2723 template <int n_components_, int dim, int spacedim, typename Number>
2725  gradient_type &
2727  const unsigned int point_index) const
2728 {
2729  AssertIndexRange(point_index, gradients.size());
2730  return gradients[point_index];
2731 }
2732 
2733 
2734 
2735 template <int n_components_, int dim, int spacedim, typename Number>
2736 inline void
2738  const value_type &value,
2739  const unsigned int point_index)
2740 {
2741  AssertIndexRange(point_index, n_q_points);
2742  values[point_index] = value;
2743 }
2744 
2745 
2746 
2747 template <int n_components_, int dim, int spacedim, typename Number>
2748 inline void
2750  const gradient_type &gradient,
2751  const unsigned int point_index)
2752 {
2753  AssertIndexRange(point_index, n_q_points);
2754  gradients[point_index] = gradient;
2755 }
2756 
2757 
2758 
2759 template <int n_components_, int dim, int spacedim, typename Number>
2762  const unsigned int point_index) const
2763 {
2764  AssertIndexRange(point_index, n_q_points);
2765  Assert(jacobian_ptr != nullptr,
2768  "update_jacobians"));
2769  return jacobian_ptr[cell_type <= ::internal::MatrixFreeFunctions::
2771  0 :
2772  point_index];
2773 }
2774 
2775 
2776 
2777 template <int n_components_, int dim, int spacedim, typename Number>
2780  const unsigned int point_index) const
2781 {
2782  AssertIndexRange(point_index, n_q_points);
2783  Assert(inverse_jacobian_ptr != nullptr,
2786  "update_inverse_jacobians"));
2787  return inverse_jacobian_ptr
2788  [cell_type <=
2790  0 :
2791  point_index];
2792 }
2793 
2794 
2795 
2796 template <int n_components_, int dim, int spacedim, typename Number>
2797 inline Number
2799  const unsigned int point_index) const
2800 {
2801  AssertIndexRange(point_index, n_q_points);
2802  Assert(JxW_ptr != nullptr,
2805  "update_JxW_values"));
2806  return JxW_ptr[point_index];
2807 }
2808 
2809 
2810 
2811 template <int n_components_, int dim, int spacedim, typename Number>
2814  const unsigned int point_index) const
2815 {
2816  AssertIndexRange(point_index, n_q_points);
2817  Assert(normal_ptr != nullptr,
2820  "update_normal_vectors"));
2821  return normal_ptr[point_index];
2822 }
2823 
2824 
2825 
2826 template <int n_components_, int dim, int spacedim, typename Number>
2829  const unsigned int point_index) const
2830 {
2831  AssertIndexRange(point_index, n_q_points);
2832  Assert(real_point_ptr != nullptr,
2835  "update_quadrature_points"));
2836  return real_point_ptr[point_index];
2837 }
2838 
2839 
2840 
2841 template <int n_components_, int dim, int spacedim, typename Number>
2842 inline Point<dim, Number>
2844  const unsigned int point_index) const
2845 {
2846  AssertIndexRange(point_index, n_q_points);
2847  Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
2848  Point<dim, Number> unit_point;
2849  for (unsigned int d = 0; d < dim; ++d)
2851  unit_point_ptr[point_index / stride][d], point_index % stride);
2852  return unit_point;
2853 }
2854 
2855 
2856 
2857 template <int n_components_, int dim, int spacedim, typename Number>
2861 {
2862  return {0U, n_q_points};
2863 }
2864 
2866 
2867 #endif
value_type * data() const noexcept
Definition: array_view.h:550
std::size_t size() const
Definition: array_view.h:573
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
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)
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_diagonal_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
static constexpr std::size_t stride
SmartPointer< const FiniteElement< dim > > fe
void integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags)
void integrate(const ArrayView< ScalarNumber > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags)
AlignedVector< vectorized_value_type > solution_renumbered_vectorized
void do_integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags)
unsigned int current_cell_index
void submit_gradient(const gradient_type &, const unsigned int point_index)
AlignedVector<::ndarray< VectorizedArrayType, 2, dim - 1 > > shapes_faces
std::vector< gradient_type > gradients
typename ETT::gradient_type gradient_type
std::vector< std::array< bool, n_components > > nonzero_shape_function_component
unsigned int current_face_number
static constexpr unsigned int dimension
const Point< dim - 1, VectorizedArrayType > * unit_point_faces_ptr
typename ETT::value_type value_type
std::unique_ptr< NonMatching::MappingInfo< dim, spacedim, Number > > mapping_info_on_the_fly
void test_and_sum(const ArrayView< ScalarNumber > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags)
typename ETT::scalar_value_type scalar_value_type
void compute_integrate_fast(const EvaluationFlags::EvaluationFlags &integration_flags, const unsigned int n_shapes, const unsigned int qb, const vectorized_value_type value, const interface_vectorized_gradient_type gradient, vectorized_value_type *solution_values_vectorized_linear)
const Point< dim, VectorizedArrayType > * unit_point_ptr
void evaluate(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
std::vector< Polynomials::Polynomial< double > > poly
void test_and_sum(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags)
std::vector< value_type > values
AlignedVector<::ndarray< VectorizedArrayType, 2, dim > > shapes
const DerivativeForm< 1, spacedim, dim, Number > * inverse_jacobian_ptr
const value_type & get_value(const unsigned int point_index) const
void evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
const unsigned int n_q_batches
typename ETT::vectorized_value_type vectorized_value_type
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
unsigned int dofs_per_component
std_cxx20::ranges::iota_view< unsigned int, unsigned int > quadrature_point_indices() const
DerivativeForm< 1, dim, spacedim, Number > jacobian(const unsigned int point_index) const
unsigned int dofs_per_component_face
void finish_integrate_fast(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, vectorized_value_type *solution_values_vectorized_linear)
const DerivativeForm< 1, dim, spacedim, Number > * jacobian_ptr
const Point< spacedim, Number > * real_point_ptr
void prepare_evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
DerivativeForm< 1, spacedim, dim, Number > inverse_jacobian(const unsigned int point_index) const
Point< dim, Number > unit_point(const unsigned int point_index) const
const UpdateFlags update_flags
std::vector< unsigned int > renumber
AlignedVector< ScalarNumber > scratch_data_scalar
void evaluate_slow(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
void evaluate(const ArrayView< const ScalarNumber > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
Number JxW(const unsigned int point_index) const
void setup(const unsigned int first_selected_component)
void submit_value(const value_type &value, const unsigned int point_index)
typename ETT::interface_vectorized_gradient_type interface_vectorized_gradient_type
static constexpr std::size_t n_lanes_internal
const unsigned int n_q_points
static constexpr unsigned int n_components
internal::MatrixFreeFunctions::GeometryType cell_type
internal::MatrixFreeFunctions::ShapeInfo< ScalarNumber > shape_info
void integrate_fast(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags)
const unsigned int n_q_points_scalar
unsigned int component_in_base_element
std::vector< scalar_value_type > solution_renumbered
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, n_components, Number > ETT
void integrate_slow(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags)
std::shared_ptr< FEValues< dim, spacedim > > fe_values
static constexpr std::size_t n_lanes_user_interface
const Tensor< 1, spacedim, Number > * normal_ptr
FEPointEvaluation(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
boost::signals2::connection connection_is_reinitialized
void compute_evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags, const unsigned int n_shapes, const unsigned int qb, vectorized_value_type &value, interface_vectorized_gradient_type &gradient)
Point< spacedim, Number > real_point(const unsigned int point_index) const
SmartPointer< const Mapping< dim, spacedim > > mapping
SmartPointer< NonMatching::MappingInfo< dim, spacedim, Number > > mapping_info
Tensor< 1, spacedim, Number > normal_vector(const unsigned int point_index) const
const gradient_type & get_gradient(const unsigned int point_index) const
Definition: point.h:112
value_type * data() const noexcept
Definition: array_view.h:735
std::size_t size() const
Definition: array_view.h:756
Definition: tensor.h:516
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
unsigned int cell_index
Definition: grid_tools.cc:1192
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcFEPointEvaluationAccessToUninitializedMappingField(std::string arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:512
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
UpdateFlags
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
EvaluationFlags
The EvaluationFlags enum.
static const char U
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:447
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)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
void integrate_tensor_product_value_and_gradient(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values, const Point< dim, Number > &p, const bool do_add)
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_linear(const Number *values, const Point< dim, Number2 > &p)
void compute_values_of_array(::ndarray< Number, 2, dim > *shapes, const std::vector< Polynomials::Polynomial< double >> &poly, const Point< dim, Number > &p, const unsigned int derivative=1)
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_linear(const Number *values, const Point< dim, Number2 > &p)
void integrate_tensor_product_value(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 &value, Number2 *values, const Point< dim, Number > &p, const bool do_add)
static const unsigned int invalid_unsigned_int
Definition: types.h:213
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:46
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static void access(gradient_type &value, const unsigned int vector_lane, const unsigned int, const scalar_gradient_type &shape_gradient)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static scalar_gradient_type access(const gradient_type &value, const unsigned int vector_lane, const unsigned int)
static void get_gradient(vectorized_gradient_type &value, const unsigned int, const vectorized_gradient_type &result)
static void set_gradient(const vectorized_gradient_type &value, const unsigned int vector_lane, scalar_gradient_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int, scalar_value_type &result)
static ScalarNumber sum_value(const unsigned int, const vectorized_value_type &result)
static void set_gradient(const vectorized_gradient_type &value, const unsigned int, vectorized_gradient_type &result)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
static void set_zero_value(value_type &value, const unsigned int vector_lane)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int)
static void set_zero_gradient(gradient_type &value, const unsigned int vector_lane)
static void access(value_type &value, const unsigned int vector_lane, const unsigned int, const ScalarNumber &shape_value)
static void get_gradient(vectorized_gradient_type &value, const unsigned int vector_lane, const scalar_gradient_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int, scalar_value_type &result)
static void set_gradient(const vectorized_gradient_type &value, const unsigned int, vectorized_gradient_type &result)
static void get_gradient(vectorized_gradient_type &value, const unsigned int vector_lane, const scalar_gradient_type &result)
static scalar_gradient_type access(const gradient_type &value, const unsigned int vector_lane, const unsigned int)
static ScalarNumber sum_value(const unsigned int, const vectorized_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static void set_zero_value(value_type &value, const unsigned int vector_lane)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
static void set_zero_gradient(gradient_type &value, const unsigned int vector_lane)
static void set_gradient(const vectorized_gradient_type &value, const unsigned int vector_lane, scalar_gradient_type &result)
static void access(gradient_type &value, const unsigned int vector_lane, const unsigned int, const scalar_gradient_type &shape_gradient)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static void get_gradient(vectorized_gradient_type &value, const unsigned int, const vectorized_gradient_type &result)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
static void access(value_type &value, const unsigned int vector_lane, const unsigned int, const ScalarNumber &shape_value)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
static void get_gradient(interface_vectorized_gradient_type &value, const unsigned int vector_lane, const gradient_type &result)
static void set_zero_gradient(gradient_type &value, const unsigned int vector_lane)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static Tensor< 1, dim, ScalarNumber > access(const gradient_type &value, const unsigned int vector_lane, const unsigned int component)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static void set_gradient(const interface_vectorized_gradient_type &value, const unsigned int vector_lane, gradient_type &result)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static void access(gradient_type &value, const unsigned int vector_lane, const unsigned int component, const Tensor< 1, dim, ScalarNumber > &shape_gradient)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static void access(value_type &value, const unsigned int vector_lane, const unsigned int component, const ScalarNumber &shape_value)
static void read_value(const ScalarNumber vector_entry, const unsigned int component, scalar_value_type &result)
static void set_zero_value(value_type &value, const unsigned int vector_lane)
static ScalarNumber sum_value(const unsigned int component, const vectorized_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static ScalarNumber sum_value(const unsigned int component, const vectorized_value_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static void access(value_type &value, const unsigned int vector_lane, const unsigned int component, const ScalarNumber &shape_value)
static Tensor< 1, dim, ScalarNumber > access(const gradient_type &value, const unsigned int vector_lane, const unsigned int component)
static void set_zero_value(value_type &value, const unsigned int vector_lane)
Tensor< 1, n_components, VectorizedArrayType > vectorized_value_type
Tensor< 1, n_components, Tensor< 1, dim, VectorizedArrayType > > vectorized_gradient_type
static void access(gradient_type &value, const unsigned int vector_lane, const unsigned int component, const Tensor< 1, dim, ScalarNumber > &shape_gradient)
static void read_value(const ScalarNumber vector_entry, const unsigned int component, scalar_value_type &result)
Tensor< 1, n_components, Tensor< 1, dim, ScalarNumber > > scalar_gradient_type
static void set_zero_gradient(gradient_type &value, const unsigned int vector_lane)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
static void get_gradient(interface_vectorized_gradient_type &value, const unsigned int vector_lane, const gradient_type &result)
static void set_gradient(const interface_vectorized_gradient_type &value, const unsigned int vector_lane, gradient_type &result)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
static constexpr std::size_t width()
static constexpr std::size_t stride()
static value_type & get(value_type &value, unsigned int c)