Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2968-g5f01c80b02 2025-03-29 13:10:00+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
fe_evaluation.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2012 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_matrix_free_fe_evaluation_h
17#define dealii_matrix_free_fe_evaluation_h
18
19
20#include <deal.II/base/config.h>
21
28
30
43
44#include <type_traits>
45
46
48
49
50
88template <int dim,
89 int n_components_,
90 typename Number,
91 bool is_face,
92 typename VectorizedArrayType>
94 : public FEEvaluationData<dim, VectorizedArrayType, is_face>
95{
96public:
97 using number_type = Number;
98 using value_type =
99 std::conditional_t<n_components_ == 1,
100 VectorizedArrayType,
102 using gradient_type = std::conditional_t<
103 n_components_ == 1,
105 std::conditional_t<
106 n_components_ == dim,
109 using hessian_type = std::conditional_t<
110 n_components_ == 1,
112 std::conditional_t<
113 n_components_ == dim,
116 static constexpr unsigned int dimension = dim;
117 static constexpr unsigned int n_components = n_components_;
118 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
119
156 template <typename VectorType>
157 void
159 const VectorType &src,
160 const unsigned int first_index = 0,
161 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
162
191 template <typename VectorType>
192 void
194 const VectorType &src,
195 const unsigned int first_index = 0,
196 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
197
229 template <typename VectorType>
230 void
232 VectorType &dst,
233 const unsigned int first_index = 0,
234 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
235
274 template <typename VectorType>
275 void
277 VectorType &dst,
278 const unsigned int first_index = 0,
279 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
280
284 template <typename VectorType>
285 void
287 VectorType &dst,
288 const unsigned int first_index = 0,
289 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
290
313 get_dof_value(const unsigned int dof) const;
314
323 void
324 submit_dof_value(const value_type val_in, const unsigned int dof);
325
337 get_value(const unsigned int q_point) const;
338
354 void
355 submit_value(const value_type val_in, const unsigned int q_point);
356
364 template <int n_components_local = n_components,
365 typename = std::enable_if_t<n_components == n_components_local>>
366 void
368 const unsigned int q_point);
369
380 get_gradient(const unsigned int q_point) const;
381
393 get_normal_derivative(const unsigned int q_point) const;
394
409 void
410 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
411
419 template <int dim_ = dim,
420 typename = std::enable_if_t<dim_ == 1 && n_components == dim_>>
421 void
423 const unsigned int q_point);
424
440 void
442 const unsigned int q_point);
443
452 get_hessian(const unsigned int q_point) const;
453
460 get_hessian_diagonal(const unsigned int q_point) const;
461
470 get_laplacian(const unsigned int q_point) const;
471
482 get_normal_hessian(const unsigned int q_point) const;
483
498 void
499 submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
500
516 void
517 submit_normal_hessian(const value_type normal_hessian_in,
518 const unsigned int q_point);
519
527 template <int dim_ = dim, typename = std::enable_if_t<n_components_ == dim_>>
528 VectorizedArrayType
529 get_divergence(const unsigned int q_point) const;
530
546 template <int dim_ = dim, typename = std::enable_if_t<n_components_ == dim_>>
547 void
548 submit_divergence(const VectorizedArrayType div_in,
549 const unsigned int q_point);
550
559 template <int dim_ = dim, typename = std::enable_if_t<n_components_ == dim_>>
561 get_symmetric_gradient(const unsigned int q_point) const;
562
578 template <int dim_ = dim, typename = std::enable_if_t<n_components_ == dim_>>
579 void
582 const unsigned int q_point);
583
592 template <int dim_ = dim,
593 typename = std::enable_if_t<n_components_ == dim_ && dim_ != 1>>
594 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
595 get_curl(const unsigned int q_point) const;
596
610 template <int dim_ = dim,
611 typename = std::enable_if_t<n_components_ == dim_ && dim != 1>>
612 void
614 const unsigned int q_point);
615
635
643
644protected:
655 const unsigned int dof_no,
656 const unsigned int first_selected_component,
657 const unsigned int quad_no,
658 const unsigned int fe_degree,
659 const unsigned int n_q_points,
660 const bool is_interior_face,
661 const unsigned int active_fe_index,
662 const unsigned int active_quad_index,
663 const unsigned int face_type);
664
702 const Mapping<dim> &mapping,
703 const FiniteElement<dim> &fe,
704 const Quadrature<1> &quadrature,
705 const UpdateFlags update_flags,
706 const unsigned int first_selected_component,
708
716
725
730
737 template <typename VectorType, typename VectorOperation>
738 void
740 const VectorOperation &operation,
741 const std::array<VectorType *, n_components_> &vectors,
742 const std::array<
744 n_components_> &vectors_sm,
745 const std::bitset<n_lanes> &mask,
746 const bool apply_constraints = true) const;
747
755 template <typename VectorType, typename VectorOperation>
756 void
758 const VectorOperation &operation,
759 const std::array<VectorType *, n_components_> &vectors,
760 const std::array<
762 n_components_> &vectors_sm,
763 const std::bitset<n_lanes> &mask) const;
764
772 template <typename VectorType, typename VectorOperation>
773 void
775 const VectorOperation &operation,
776 const std::array<VectorType *, n_components_> &vectors) const;
777
781 void
783
788
793
798 mutable std::vector<types::global_dof_index> local_dof_indices;
799};
800
801
802
803// backward compatibility
804template <int dim,
805 int n_components_,
806 typename Number,
807 bool is_face,
808 typename VectorizedArrayType = VectorizedArray<Number>>
811
1375template <int dim,
1376 int fe_degree,
1377 int n_q_points_1d,
1378 int n_components_,
1379 typename Number,
1380 typename VectorizedArrayType>
1382 n_components_,
1383 Number,
1384 false,
1385 VectorizedArrayType>
1386{
1387 static_assert(
1388 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1389 "Type of Number and of VectorizedArrayType do not match.");
1390
1391public:
1397
1401 using number_type = Number;
1402
1409
1416
1420 static constexpr unsigned int dimension = dim;
1421
1426 static constexpr unsigned int n_components = n_components_;
1427
1431 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1432
1441 static constexpr unsigned int static_n_q_points =
1442 Utilities::pow(n_q_points_1d, dim);
1443
1453 static constexpr unsigned int static_dofs_per_component =
1454 Utilities::pow(fe_degree + 1, dim);
1455
1465 static constexpr unsigned int tensor_dofs_per_cell =
1467
1477 static constexpr unsigned int static_dofs_per_cell =
1479
1516 const unsigned int dof_no = 0,
1517 const unsigned int quad_no = 0,
1518 const unsigned int first_selected_component = 0,
1521
1530 const std::pair<unsigned int, unsigned int> &range,
1531 const unsigned int dof_no = 0,
1532 const unsigned int quad_no = 0,
1533 const unsigned int first_selected_component = 0);
1534
1564 const FiniteElement<dim> &fe,
1565 const Quadrature<1> &quadrature,
1566 const UpdateFlags update_flags,
1567 const unsigned int first_selected_component = 0);
1568
1575 const Quadrature<1> &quadrature,
1576 const UpdateFlags update_flags,
1577 const unsigned int first_selected_component = 0);
1578
1591 const unsigned int first_selected_component = 0);
1592
1600
1607 FEEvaluation &
1609
1618 void
1619 reinit(const unsigned int cell_batch_index);
1620
1627 void
1628 reinit(const std::array<unsigned int, n_lanes> &cell_ids);
1629
1642 template <bool level_dof_access>
1643 void
1645
1656 void
1658
1662 static bool
1663 fast_evaluation_supported(const unsigned int given_degree,
1664 const unsigned int given_n_q_points_1d);
1665
1675 void
1677
1690 void
1691 evaluate(const VectorizedArrayType *values_array,
1692 const EvaluationFlags::EvaluationFlags evaluation_flag);
1693
1707 template <typename VectorType>
1708 void
1709 gather_evaluate(const VectorType &input_vector,
1710 const EvaluationFlags::EvaluationFlags evaluation_flag);
1711
1721 void
1723
1735 void
1737 VectorizedArrayType *values_array,
1738 const bool sum_into_values = false);
1739
1753 template <typename VectorType>
1754 void
1756 VectorType &output_vector);
1757
1765
1772 const unsigned int dofs_per_component;
1773
1780 const unsigned int dofs_per_cell;
1781
1789 const unsigned int n_q_points;
1790
1791private:
1796 void
1797 check_template_arguments(const unsigned int fe_no,
1798 const unsigned int first_selected_component);
1799};
1800
1801
1802
1838template <int dim,
1839 int fe_degree,
1840 int n_q_points_1d = fe_degree + 1,
1841 int n_components_ = 1,
1842 typename Number = double,
1843 typename VectorizedArrayType = VectorizedArray<Number>>
1845 n_components_,
1846 Number,
1847 true,
1848 VectorizedArrayType>
1849{
1850 static_assert(
1851 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1852 "Type of Number and of VectorizedArrayType do not match.");
1853
1854public:
1860
1864 using number_type = Number;
1865
1872
1879
1883 static constexpr unsigned int dimension = dim;
1884
1889 static constexpr unsigned int n_components = n_components_;
1890
1894 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1895
1905 static constexpr unsigned int static_n_q_points =
1906 Utilities::pow(n_q_points_1d, dim - 1);
1907
1916 static constexpr unsigned int static_n_q_points_cell =
1917 Utilities::pow(n_q_points_1d, dim);
1918
1927 static constexpr unsigned int static_dofs_per_component =
1928 Utilities::pow(fe_degree + 1, dim);
1929
1938 static constexpr unsigned int tensor_dofs_per_cell =
1940
1949 static constexpr unsigned int static_dofs_per_cell =
1951
1995 const bool is_interior_face = true,
1996 const unsigned int dof_no = 0,
1997 const unsigned int quad_no = 0,
1998 const unsigned int first_selected_component = 0,
2001 const unsigned int face_type = numbers::invalid_unsigned_int);
2002
2012 const std::pair<unsigned int, unsigned int> &range,
2013 const bool is_interior_face = true,
2014 const unsigned int dof_no = 0,
2015 const unsigned int quad_no = 0,
2016 const unsigned int first_selected_component = 0);
2017
2028 void
2029 reinit(const unsigned int face_batch_number);
2030
2038 void
2039 reinit(const unsigned int cell_batch_number, const unsigned int face_number);
2040
2044 static bool
2045 fast_evaluation_supported(const unsigned int given_degree,
2046 const unsigned int given_n_q_points_1d);
2047
2058 void
2060
2073 void
2074 evaluate(const VectorizedArrayType *values_array,
2075 const EvaluationFlags::EvaluationFlags evaluation_flag);
2076
2081 void
2083
2088 void
2089 project_to_face(const VectorizedArrayType *values_array,
2090 const EvaluationFlags::EvaluationFlags evaluation_flag);
2091
2096 void
2098
2110 template <typename VectorType>
2111 void
2112 gather_evaluate(const VectorType &input_vector,
2113 const EvaluationFlags::EvaluationFlags evaluation_flag);
2114
2124 void
2126 const bool sum_into_values = false);
2127
2137 void
2139 VectorizedArrayType *values_array,
2140 const bool sum_into_values = false);
2141
2148 void
2150
2155 void
2157 const bool sum_into_values = false);
2158
2163 void
2165 VectorizedArrayType *values_array,
2166 const bool sum_into_values = false);
2167
2179 template <typename VectorType>
2180 void
2182 VectorType &output_vector);
2183
2187 template <typename VectorType>
2188 void
2189 integrate_scatter(const bool integrate_values,
2190 const bool integrate_gradients,
2191 VectorType &output_vector);
2192
2200
2205 bool
2207
2222
2227 unsigned int
2229
2234 unsigned int
2236
2243 const unsigned int dofs_per_component;
2244
2251 const unsigned int dofs_per_cell;
2252
2260 const unsigned int n_q_points;
2261};
2262
2263
2264
2265namespace internal
2266{
2267 namespace MatrixFreeFunctions
2268 {
2269 // a helper function to compute the number of DoFs of a DGP element at
2270 // compile time, depending on the degree
2271 template <int dim, int degree>
2273 {
2274 // this division is always without remainder
2275 static constexpr unsigned int value =
2276 (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
2277 };
2278
2279 // base specialization: 1d elements have 'degree+1' degrees of freedom
2280 template <int degree>
2281 struct DGP_dofs_per_component<1, degree>
2282 {
2283 static constexpr unsigned int value = degree + 1;
2284 };
2285 } // namespace MatrixFreeFunctions
2286} // namespace internal
2287
2288
2289/*----------------------- Inline functions ----------------------------------*/
2290
2291#ifndef DOXYGEN
2292
2293
2294namespace internal
2295{
2296 // Extract all internal data pointers and indices in a single function that
2297 // get passed on to the constructor of FEEvaluationData, avoiding to look
2298 // things up multiple times
2299 template <bool is_face,
2300 int dim,
2301 typename Number,
2302 typename VectorizedArrayType>
2304 InitializationData
2305 extract_initialization_data(
2307 const unsigned int dof_no,
2308 const unsigned int first_selected_component,
2309 const unsigned int quad_no,
2310 const unsigned int fe_degree,
2311 const unsigned int n_q_points,
2312 const unsigned int active_fe_index_given,
2313 const unsigned int active_quad_index_given,
2314 const unsigned int face_type)
2315 {
2317 InitializationData init_data;
2318
2319 init_data.dof_info = &matrix_free.get_dof_info(dof_no);
2320 init_data.mapping_data =
2321 &internal::MatrixFreeFunctions::
2322 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2323 matrix_free.get_mapping_info(), quad_no);
2324
2325 init_data.active_fe_index =
2326 fe_degree != numbers::invalid_unsigned_int ?
2327 init_data.dof_info->fe_index_from_degree(first_selected_component,
2328 fe_degree) :
2329 (active_fe_index_given != numbers::invalid_unsigned_int ?
2330 active_fe_index_given :
2331 0);
2332
2333 init_data.active_quad_index =
2334 fe_degree == numbers::invalid_unsigned_int ?
2335 (active_quad_index_given != numbers::invalid_unsigned_int ?
2336 active_quad_index_given :
2337 std::min<unsigned int>(
2338 init_data.active_fe_index,
2339 init_data.mapping_data->descriptor.size() /
2340 (is_face ? std::max<unsigned int>(1, dim - 1) : 1) -
2341 1)) :
2342 init_data.mapping_data->quad_index_from_n_q_points(n_q_points);
2343
2344 init_data.shape_info = &matrix_free.get_shape_info(
2345 dof_no,
2346 quad_no,
2347 init_data.dof_info->component_to_base_index[first_selected_component],
2348 init_data.active_fe_index,
2349 init_data.active_quad_index);
2350 init_data.descriptor =
2351 &init_data.mapping_data->descriptor
2352 [is_face ?
2353 (init_data.active_quad_index * std::max<unsigned int>(1, dim - 1) +
2354 (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
2355 init_data.active_quad_index];
2356
2357 return init_data;
2358 }
2359} // namespace internal
2360
2361
2362
2363/*----------------------- FEEvaluationBase ----------------------------------*/
2364
2365template <int dim,
2366 int n_components_,
2367 typename Number,
2368 bool is_face,
2369 typename VectorizedArrayType>
2370inline FEEvaluationBase<dim,
2371 n_components_,
2372 Number,
2373 is_face,
2374 VectorizedArrayType>::
2375 FEEvaluationBase(
2377 const unsigned int dof_no,
2378 const unsigned int first_selected_component,
2379 const unsigned int quad_no,
2380 const unsigned int fe_degree,
2381 const unsigned int n_q_points,
2382 const bool is_interior_face,
2383 const unsigned int active_fe_index,
2384 const unsigned int active_quad_index,
2385 const unsigned int face_type)
2386 : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2387 internal::extract_initialization_data<is_face>(matrix_free,
2388 dof_no,
2389 first_selected_component,
2390 quad_no,
2391 fe_degree,
2392 n_q_points,
2393 active_fe_index,
2394 active_quad_index,
2395 face_type),
2396 is_interior_face,
2397 quad_no,
2398 first_selected_component)
2399 , scratch_data_array(matrix_free.acquire_scratch_data())
2400 , matrix_free(&matrix_free)
2401{
2402 this->set_data_pointers(scratch_data_array, n_components_);
2403 Assert(
2404 this->dof_info->start_components.back() == 1 ||
2405 static_cast<int>(n_components_) <=
2406 static_cast<int>(
2407 this->dof_info->start_components
2408 [this->dof_info->component_to_base_index[first_selected_component] +
2409 1]) -
2410 first_selected_component,
2411 ExcMessage(
2412 "You tried to construct a vector-valued evaluator with " +
2413 std::to_string(n_components) +
2414 " components. However, "
2415 "the current base element has only " +
2416 std::to_string(
2417 this->dof_info->start_components
2418 [this->dof_info->component_to_base_index[first_selected_component] +
2419 1] -
2420 first_selected_component) +
2421 " components left when starting from local element index " +
2422 std::to_string(
2423 first_selected_component -
2424 this->dof_info->start_components
2425 [this->dof_info->component_to_base_index[first_selected_component]]) +
2426 " (global index " + std::to_string(first_selected_component) + ")"));
2427
2428 // do not check for correct dimensions of data fields here, should be done
2429 // in derived classes
2430}
2431
2432
2433
2434template <int dim,
2435 int n_components_,
2436 typename Number,
2437 bool is_face,
2438 typename VectorizedArrayType>
2439inline FEEvaluationBase<dim,
2440 n_components_,
2441 Number,
2442 is_face,
2443 VectorizedArrayType>::
2444 FEEvaluationBase(
2445 const Mapping<dim> &mapping,
2446 const FiniteElement<dim> &fe,
2447 const Quadrature<1> &quadrature,
2448 const UpdateFlags update_flags,
2449 const unsigned int first_selected_component,
2451 : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2452 other != nullptr &&
2453 other->mapped_geometry->get_quadrature() == quadrature ?
2454 other->mapped_geometry :
2455 std::make_shared<internal::MatrixFreeFunctions::
2456 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2457 mapping,
2458 quadrature,
2459 update_flags),
2460 n_components_,
2461 first_selected_component)
2462 , scratch_data_array(new AlignedVector<VectorizedArrayType>())
2463 , matrix_free(nullptr)
2464{
2465 const unsigned int base_element_number =
2466 fe.component_to_base_index(first_selected_component).first;
2467 Assert(fe.element_multiplicity(base_element_number) == 1 ||
2468 fe.element_multiplicity(base_element_number) -
2469 first_selected_component >=
2470 n_components_,
2471 ExcMessage("The underlying element must at least contain as many "
2472 "components as requested by this class"));
2473 (void)base_element_number;
2474
2475 Assert(this->data == nullptr, ExcInternalError());
2477 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2478 fe,
2479 fe.component_to_base_index(first_selected_component).first);
2480
2481 this->set_data_pointers(scratch_data_array, n_components_);
2482}
2483
2484
2485
2486template <int dim,
2487 int n_components_,
2488 typename Number,
2489 bool is_face,
2490 typename VectorizedArrayType>
2491inline FEEvaluationBase<dim,
2492 n_components_,
2493 Number,
2494 is_face,
2495 VectorizedArrayType>::
2496 FEEvaluationBase(const FEEvaluationBase<dim,
2497 n_components_,
2498 Number,
2499 is_face,
2500 VectorizedArrayType> &other)
2501 : FEEvaluationData<dim, VectorizedArrayType, is_face>(other)
2502 , scratch_data_array(other.matrix_free == nullptr ?
2503 new AlignedVector<VectorizedArrayType>() :
2504 other.matrix_free->acquire_scratch_data())
2505 , matrix_free(other.matrix_free)
2506{
2507 if (other.matrix_free == nullptr)
2508 {
2509 Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
2510 this->data =
2512
2513 // Create deep copy of mapped geometry for use in parallel
2514 this->mapped_geometry =
2515 std::make_shared<internal::MatrixFreeFunctions::
2516 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2517 other.mapped_geometry->get_fe_values().get_mapping(),
2518 other.mapped_geometry->get_quadrature(),
2519 other.mapped_geometry->get_fe_values().get_update_flags());
2520
2521 if constexpr (is_face == false)
2522 this->mapping_data = &this->mapped_geometry->get_data_storage();
2523 else
2524 Assert(false,
2525 ExcNotImplemented("On-the-fly geometry evaluation with "
2526 "face evaluators is not currently "
2527 "implemented!"));
2528
2529 this->cell = 0;
2530
2531 this->jacobian =
2532 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2533 this->J_value =
2534 this->mapped_geometry->get_data_storage().JxW_values.begin();
2535 this->jacobian_gradients =
2536 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2537 this->jacobian_gradients_non_inverse =
2538 this->mapped_geometry->get_data_storage()
2539 .jacobian_gradients_non_inverse[0]
2540 .begin();
2541 this->quadrature_points =
2542 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2543 }
2544
2545 this->set_data_pointers(scratch_data_array, n_components_);
2546}
2547
2548
2549
2550template <int dim,
2551 int n_components_,
2552 typename Number,
2553 bool is_face,
2554 typename VectorizedArrayType>
2555inline FEEvaluationBase<dim,
2556 n_components_,
2557 Number,
2558 is_face,
2559 VectorizedArrayType> &
2561operator=(const FEEvaluationBase<dim,
2562 n_components_,
2563 Number,
2564 is_face,
2565 VectorizedArrayType> &other)
2566{
2567 // release old memory
2568 if (matrix_free == nullptr)
2569 {
2570 delete this->data;
2571 delete scratch_data_array;
2572 }
2573 else
2574 {
2575 matrix_free->release_scratch_data(scratch_data_array);
2576 }
2577
2579
2580 matrix_free = other.matrix_free;
2581
2582 if (other.matrix_free == nullptr)
2583 {
2584 Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
2585 this->data =
2587 scratch_data_array = new AlignedVector<VectorizedArrayType>();
2588
2589 // Create deep copy of mapped geometry for use in parallel
2590 this->mapped_geometry =
2591 std::make_shared<internal::MatrixFreeFunctions::
2592 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2593 other.mapped_geometry->get_fe_values().get_mapping(),
2594 other.mapped_geometry->get_quadrature(),
2595 other.mapped_geometry->get_fe_values().get_update_flags());
2596
2597 if constexpr (is_face == false)
2598 this->mapping_data = &this->mapped_geometry->get_data_storage();
2599 else
2600 Assert(false,
2601 ExcNotImplemented("On-the-fly geometry evaluation with "
2602 "face evaluators is not currently "
2603 "implemented!"));
2604 this->cell = 0;
2605
2606 this->jacobian =
2607 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2608 this->J_value =
2609 this->mapped_geometry->get_data_storage().JxW_values.begin();
2610 this->jacobian_gradients =
2611 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2612 this->jacobian_gradients_non_inverse =
2613 this->mapped_geometry->get_data_storage()
2614 .jacobian_gradients_non_inverse[0]
2615 .begin();
2616 this->quadrature_points =
2617 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2618 }
2619 else
2620 {
2621 scratch_data_array = matrix_free->acquire_scratch_data();
2622 }
2623
2624 this->set_data_pointers(scratch_data_array, n_components_);
2625
2626 return *this;
2627}
2628
2629
2630
2631template <int dim,
2632 int n_components_,
2633 typename Number,
2634 bool is_face,
2635 typename VectorizedArrayType>
2636inline FEEvaluationBase<dim,
2637 n_components_,
2638 Number,
2639 is_face,
2640 VectorizedArrayType>::~FEEvaluationBase()
2641{
2642 if (matrix_free != nullptr)
2643 {
2644 try
2645 {
2646 matrix_free->release_scratch_data(scratch_data_array);
2647 }
2648 catch (...)
2649 {}
2650 }
2651 else
2652 {
2653 delete scratch_data_array;
2654 delete this->data;
2655 }
2656}
2657
2658
2659
2660template <int dim,
2661 int n_components_,
2662 typename Number,
2663 bool is_face,
2664 typename VectorizedArrayType>
2667 get_matrix_free() const
2668{
2669 Assert(matrix_free != nullptr,
2670 ExcMessage(
2671 "FEEvaluation was not initialized with a MatrixFree object!"));
2672 return *matrix_free;
2673}
2674
2675
2676
2677namespace internal
2678{
2679 // given a block vector return the underlying vector type
2680 // including constness (specified by bool)
2681 template <typename VectorType, bool>
2682 struct ConstBlockVectorSelector;
2683
2684 template <typename VectorType>
2685 struct ConstBlockVectorSelector<VectorType, true>
2686 {
2687 using BaseVectorType = const typename VectorType::BlockType;
2688 };
2689
2690 template <typename VectorType>
2691 struct ConstBlockVectorSelector<VectorType, false>
2692 {
2693 using BaseVectorType = typename VectorType::BlockType;
2694 };
2695
2696 // allows to select between block vectors and non-block vectors, which
2697 // allows to use a unified interface for extracting blocks on block vectors
2698 // and doing nothing on usual vectors
2699 template <typename VectorType, bool>
2700 struct BlockVectorSelector;
2701
2702 template <typename VectorType>
2703 struct BlockVectorSelector<VectorType, true>
2704 {
2705 using BaseVectorType = typename ConstBlockVectorSelector<
2706 VectorType,
2707 std::is_const_v<VectorType>>::BaseVectorType;
2708
2709 static BaseVectorType *
2710 get_vector_component(VectorType &vec, const unsigned int component)
2711 {
2712 AssertIndexRange(component, vec.n_blocks());
2713 return &vec.block(component);
2714 }
2715 };
2716
2717 template <typename VectorType>
2718 struct BlockVectorSelector<VectorType, false>
2719 {
2720 using BaseVectorType = VectorType;
2721
2722 static BaseVectorType *
2723 get_vector_component(VectorType &vec, const unsigned int component)
2724 {
2725 // FEEvaluation allows to combine several vectors from a scalar
2726 // FiniteElement into a "vector-valued" FEEvaluation object with
2727 // multiple components. These components can be extracted with the other
2728 // get_vector_component functions. If we do not get a vector of vectors
2729 // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
2730 // must make sure that we do not duplicate the components in input
2731 // and/or duplicate the resulting integrals. In such a case, we should
2732 // only get the zeroth component in the vector contained set nullptr for
2733 // the others which allows us to catch unintended use in
2734 // read_write_operation.
2735 if (component == 0)
2736 return &vec;
2737 else
2738 return nullptr;
2739 }
2740 };
2741
2742 template <typename VectorType>
2743 struct BlockVectorSelector<std::vector<VectorType>, false>
2744 {
2745 using BaseVectorType = VectorType;
2746
2747 static BaseVectorType *
2748 get_vector_component(std::vector<VectorType> &vec,
2749 const unsigned int component)
2750 {
2751 AssertIndexRange(component, vec.size());
2752 return &vec[component];
2753 }
2754 };
2755
2756 template <typename VectorType>
2757 struct BlockVectorSelector<const std::vector<VectorType>, false>
2758 {
2759 using BaseVectorType = const VectorType;
2760
2761 static const BaseVectorType *
2762 get_vector_component(const std::vector<VectorType> &vec,
2763 const unsigned int component)
2764 {
2765 AssertIndexRange(component, vec.size());
2766 return &vec[component];
2767 }
2768 };
2769
2770 template <typename VectorType>
2771 struct BlockVectorSelector<std::vector<VectorType *>, false>
2772 {
2773 using BaseVectorType = VectorType;
2774
2775 static BaseVectorType *
2776 get_vector_component(std::vector<VectorType *> &vec,
2777 const unsigned int component)
2778 {
2779 AssertIndexRange(component, vec.size());
2780 return vec[component];
2781 }
2782 };
2783
2784 template <typename VectorType>
2785 struct BlockVectorSelector<const std::vector<VectorType *>, false>
2786 {
2787 using BaseVectorType = const VectorType;
2788
2789 static const BaseVectorType *
2790 get_vector_component(const std::vector<VectorType *> &vec,
2791 const unsigned int component)
2792 {
2793 AssertIndexRange(component, vec.size());
2794 return vec[component];
2795 }
2796 };
2797
2798 template <typename VectorType, std::size_t N>
2799 struct BlockVectorSelector<std::array<VectorType *, N>, false>
2800 {
2801 using BaseVectorType = VectorType;
2802
2803 static BaseVectorType *
2804 get_vector_component(std::array<VectorType *, N> &vec,
2805 const unsigned int component)
2806 {
2807 AssertIndexRange(component, vec.size());
2808 return vec[component];
2809 }
2810 };
2811} // namespace internal
2812
2813
2814
2815template <int dim,
2816 int n_components_,
2817 typename Number,
2818 bool is_face,
2819 typename VectorizedArrayType>
2820template <typename VectorType, typename VectorOperation>
2821inline void
2824 const VectorOperation &operation,
2825 const std::array<VectorType *, n_components_> &src,
2826 const std::array<
2828 n_components_> &src_sm,
2829 const std::bitset<n_lanes> &mask,
2830 const bool apply_constraints) const
2831{
2832 // Case 1: No MatrixFree object given, simple case because we do not need to
2833 // process constraints and need not care about vectorization -> go to
2834 // separate function
2835 if (this->matrix_free == nullptr)
2836 {
2837 read_write_operation_global(operation, src);
2838 return;
2839 }
2840
2841 Assert(this->dof_info != nullptr, ExcNotInitialized());
2842 const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
2843 Assert(this->matrix_free->indices_initialized() == true, ExcNotInitialized());
2844 if (this->n_fe_components == 1)
2845 for (unsigned int comp = 0; comp < n_components; ++comp)
2846 {
2847 Assert(src[comp] != nullptr,
2848 ExcMessage("The finite element underlying this FEEvaluation "
2849 "object is scalar, but you requested " +
2850 std::to_string(n_components) +
2851 " components via the template argument in "
2852 "FEEvaluation. In that case, you must pass an "
2853 "std::vector<VectorType> or a BlockVector to " +
2854 "read_dof_values and distribute_local_to_global."));
2856 *this->matrix_free,
2857 *this->dof_info);
2858 }
2859 else
2860 {
2862 *this->matrix_free,
2863 *this->dof_info);
2864 }
2865
2866 const bool accesses_exterior_dofs =
2867 this->dof_access_index ==
2869 this->is_interior_face() == false;
2870
2871 // Case 2: contiguous indices which use reduced storage of indices and can
2872 // use vectorized load/store operations -> go to separate function
2873 if (this->cell != numbers::invalid_unsigned_int)
2874 {
2876 this->cell,
2877 dof_info.index_storage_variants[this->dof_access_index].size());
2878
2879 bool is_contiguous = true;
2880 // check if exterior cells are not contiguous (ECL case)
2881 if (accesses_exterior_dofs)
2882 {
2883 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
2884 const unsigned int n_filled_lanes =
2887 [this->cell];
2888 // we have to check all filled lanes which are active in the mask
2889 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2890 if (mask[v] == true &&
2891 dof_info.index_storage_variants
2893 [cells[v] / n_lanes] <
2895 contiguous)
2896 is_contiguous = false;
2897 } // or if cell/face batch is not contiguous
2898 else if (dof_info.index_storage_variants
2899 [is_face ?
2900 this->dof_access_index :
2901 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2902 [this->cell] < internal::MatrixFreeFunctions::DoFInfo::
2903 IndexStorageVariants::contiguous)
2904 {
2905 is_contiguous = false;
2906 }
2907
2908 if (is_contiguous)
2909 {
2910 read_write_operation_contiguous(operation, src, src_sm, mask);
2911 return;
2912 }
2913 }
2914
2915 // Case 3: standard operation with one index per degree of freedom -> go on
2916 // here
2917 std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
2918
2919 const bool masking_is_active = mask.count() < n_lanes;
2920 if (masking_is_active)
2921 for (unsigned int v = 0; v < n_lanes; ++v)
2922 if (mask[v] == false)
2924
2925 bool has_hn_constraints = false;
2926
2927 if (is_face == false)
2928 {
2929 if (!dof_info.hanging_node_constraint_masks.empty() &&
2930 !dof_info.hanging_node_constraint_masks_comp.empty() &&
2931 dof_info
2932 .hanging_node_constraint_masks_comp[this->active_fe_index]
2933 [this->first_selected_component])
2934 for (unsigned int v = 0; v < n_lanes; ++v)
2935 if (cells[v] != numbers::invalid_unsigned_int &&
2936 dof_info.hanging_node_constraint_masks[cells[v]] !=
2939 has_hn_constraints = true;
2940 }
2941
2942 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
2943 vector_selector;
2944
2945 const bool use_vectorized_path =
2946 !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
2947
2948 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
2949 std::array<VectorizedArrayType *, n_components> values_dofs;
2950 for (unsigned int c = 0; c < n_components; ++c)
2951 values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
2952 c * dofs_per_component;
2953
2954 if (this->cell != numbers::invalid_unsigned_int &&
2955 dof_info.index_storage_variants
2956 [is_face ? this->dof_access_index :
2957 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2958 [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
2959 IndexStorageVariants::interleaved &&
2960 use_vectorized_path)
2961 {
2962 const unsigned int *dof_indices =
2963 dof_info.dof_indices_interleaved.data() +
2964 dof_info.row_starts[this->cell * this->n_fe_components * n_lanes]
2965 .first +
2966 this->dof_info
2967 ->component_dof_indices_offset[this->active_fe_index]
2968 [this->first_selected_component] *
2969 n_lanes;
2970
2971 std::array<typename VectorType::value_type *, n_components> src_ptrs;
2972 if (n_components == 1 || this->n_fe_components == 1)
2973 for (unsigned int comp = 0; comp < n_components; ++comp)
2974 src_ptrs[comp] =
2975 const_cast<typename VectorType::value_type *>(src[comp]->begin());
2976 else
2977 src_ptrs[0] =
2978 const_cast<typename VectorType::value_type *>(src[0]->begin());
2979
2980 if (n_components == 1 || this->n_fe_components == 1)
2981 for (unsigned int i = 0; i < dofs_per_component;
2982 ++i, dof_indices += n_lanes)
2983 for (unsigned int comp = 0; comp < n_components; ++comp)
2984 operation.process_dof_gather(dof_indices,
2985 *src[comp],
2986 0,
2987 src_ptrs[comp],
2988 values_dofs[comp][i],
2989 vector_selector);
2990 else
2991 for (unsigned int comp = 0; comp < n_components; ++comp)
2992 for (unsigned int i = 0; i < dofs_per_component;
2993 ++i, dof_indices += n_lanes)
2994 operation.process_dof_gather(dof_indices,
2995 *src[0],
2996 0,
2997 src_ptrs[0],
2998 values_dofs[comp][i],
2999 vector_selector);
3000 return;
3001 }
3002
3003 // Allocate pointers, then initialize all of them to nullptrs and
3004 // below overwrite the ones we actually use:
3005 std::array<const unsigned int *, n_lanes> dof_indices;
3006 dof_indices.fill(nullptr);
3007
3008 // Assign the appropriate cell ids for face/cell case and get the pointers
3009 // to the dof indices of the cells on all lanes
3010
3011 bool has_constraints = false;
3012 const unsigned int n_components_read =
3013 this->n_fe_components > 1 ? n_components : 1;
3014
3015 if (is_face)
3016 {
3017 for (unsigned int v = 0; v < n_lanes; ++v)
3018 {
3019 if (cells[v] == numbers::invalid_unsigned_int)
3020 continue;
3021
3022 Assert(cells[v] < dof_info.row_starts.size() - 1, ExcInternalError());
3023 const std::pair<unsigned int, unsigned int> *my_index_start =
3024 &dof_info.row_starts[cells[v] * this->n_fe_components +
3025 this->first_selected_component];
3026
3027 // check whether any of the SIMD lanes has constraints, i.e., the
3028 // constraint indicator which is the second entry of row_starts
3029 // increments on this cell
3030 if (my_index_start[n_components_read].second !=
3031 my_index_start[0].second)
3032 has_constraints = true;
3033
3034 dof_indices[v] =
3035 dof_info.dof_indices.data() + my_index_start[0].first;
3036 }
3037 }
3038 else
3039 {
3040 for (unsigned int v = 0; v < n_lanes; ++v)
3041 {
3042 if (cells[v] == numbers::invalid_unsigned_int)
3043 continue;
3044
3045 const std::pair<unsigned int, unsigned int> *my_index_start =
3046 &dof_info.row_starts[cells[v] * this->n_fe_components +
3047 this->first_selected_component];
3048 if (my_index_start[n_components_read].second !=
3049 my_index_start[0].second)
3050 has_constraints = true;
3051
3052 if (dof_info.hanging_node_constraint_masks.size() > 0 &&
3053 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
3054 dof_info.hanging_node_constraint_masks[cells[v]] !=
3057 dof_info.hanging_node_constraint_masks_comp
3058 [this->active_fe_index][this->first_selected_component])
3059 has_hn_constraints = true;
3060
3061 Assert(my_index_start[n_components_read].first ==
3062 my_index_start[0].first ||
3063 my_index_start[0].first < dof_info.dof_indices.size(),
3064 ExcIndexRange(0,
3065 my_index_start[0].first,
3066 dof_info.dof_indices.size()));
3067 dof_indices[v] =
3068 dof_info.dof_indices.data() + my_index_start[0].first;
3069 }
3070 }
3071
3072 if (std::count_if(cells.begin(), cells.end(), [](const auto i) {
3073 return i != numbers::invalid_unsigned_int;
3074 }) < n_lanes)
3075 for (unsigned int comp = 0; comp < n_components; ++comp)
3076 for (unsigned int i = 0; i < dofs_per_component; ++i)
3077 operation.process_empty(values_dofs[comp][i]);
3078
3079 // Case where we have no constraints throughout the whole cell: Can go
3080 // through the list of DoFs directly
3081 if (!has_constraints && apply_constraints)
3082 {
3083 if (n_components == 1 || this->n_fe_components == 1)
3084 {
3085 for (unsigned int v = 0; v < n_lanes; ++v)
3086 {
3087 if (cells[v] == numbers::invalid_unsigned_int)
3088 continue;
3089
3090 for (unsigned int i = 0; i < dofs_per_component; ++i)
3091 for (unsigned int comp = 0; comp < n_components; ++comp)
3092 operation.process_dof(dof_indices[v][i],
3093 *src[comp],
3094 values_dofs[comp][i][v]);
3095 }
3096 }
3097 else
3098 {
3099 for (unsigned int comp = 0; comp < n_components; ++comp)
3100 for (unsigned int v = 0; v < n_lanes; ++v)
3101 {
3102 if (cells[v] == numbers::invalid_unsigned_int)
3103 continue;
3104
3105 for (unsigned int i = 0; i < dofs_per_component; ++i)
3106 operation.process_dof(
3107 dof_indices[v][comp * dofs_per_component + i],
3108 *src[0],
3109 values_dofs[comp][i][v]);
3110 }
3111 }
3112 return;
3113 }
3114
3115 // In the case where there are some constraints to be resolved, loop over
3116 // all vector components that are filled and then over local dofs. ind_local
3117 // holds local number on cell, index iterates over the elements of
3118 // index_local_to_global and dof_indices points to the global indices stored
3119 // in index_local_to_global
3120
3121 for (unsigned int v = 0; v < n_lanes; ++v)
3122 {
3123 if (cells[v] == numbers::invalid_unsigned_int)
3124 continue;
3125
3126 const unsigned int cell_index = cells[v];
3127 const unsigned int cell_dof_index =
3128 cell_index * this->n_fe_components + this->first_selected_component;
3129 const unsigned int n_components_read =
3130 this->n_fe_components > 1 ? n_components : 1;
3131 unsigned int index_indicators =
3132 dof_info.row_starts[cell_dof_index].second;
3133 unsigned int next_index_indicators =
3134 dof_info.row_starts[cell_dof_index + 1].second;
3135
3136 // For read_dof_values_plain, redirect the dof_indices field to the
3137 // unconstrained indices
3138 if (apply_constraints == false &&
3139 (dof_info.row_starts[cell_dof_index].second !=
3140 dof_info.row_starts[cell_dof_index + n_components_read].second ||
3141 ((dof_info.hanging_node_constraint_masks.size() > 0 &&
3142 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
3146 dof_info.hanging_node_constraint_masks_comp
3147 [this->active_fe_index][this->first_selected_component])))
3148 {
3152 dof_indices[v] =
3153 dof_info.plain_dof_indices.data() +
3154 this->dof_info
3155 ->component_dof_indices_offset[this->active_fe_index]
3156 [this->first_selected_component] +
3158 next_index_indicators = index_indicators;
3159 }
3160
3161 if (n_components == 1 || this->n_fe_components == 1)
3162 {
3163 unsigned int ind_local = 0;
3164 for (; index_indicators != next_index_indicators; ++index_indicators)
3165 {
3166 const std::pair<unsigned short, unsigned short> indicator =
3167 dof_info.constraint_indicator[index_indicators];
3168 // run through values up to next constraint
3169 for (unsigned int j = 0; j < indicator.first; ++j)
3170 for (unsigned int comp = 0; comp < n_components; ++comp)
3171 operation.process_dof(dof_indices[v][j],
3172 *src[comp],
3173 values_dofs[comp][ind_local + j][v]);
3174
3175 ind_local += indicator.first;
3176 dof_indices[v] += indicator.first;
3177
3178 // constrained case: build the local value as a linear
3179 // combination of the global value according to constraints
3180 Number value[n_components];
3181 for (unsigned int comp = 0; comp < n_components; ++comp)
3182 operation.pre_constraints(values_dofs[comp][ind_local][v],
3183 value[comp]);
3184
3185 const Number *data_val =
3186 this->matrix_free->constraint_pool_begin(indicator.second);
3187 const Number *end_pool =
3188 this->matrix_free->constraint_pool_end(indicator.second);
3189 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3190 for (unsigned int comp = 0; comp < n_components; ++comp)
3191 operation.process_constraint(*dof_indices[v],
3192 *data_val,
3193 *src[comp],
3194 value[comp]);
3195
3196 for (unsigned int comp = 0; comp < n_components; ++comp)
3197 operation.post_constraints(value[comp],
3198 values_dofs[comp][ind_local][v]);
3199 ++ind_local;
3200 }
3201
3202 AssertIndexRange(ind_local, dofs_per_component + 1);
3203
3204 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3205 for (unsigned int comp = 0; comp < n_components; ++comp)
3206 operation.process_dof(*dof_indices[v],
3207 *src[comp],
3208 values_dofs[comp][ind_local][v]);
3209 }
3210 else
3211 {
3212 // case with vector-valued finite elements where all components are
3213 // included in one single vector. Assumption: first come all entries
3214 // to the first component, then all entries to the second one, and
3215 // so on. This is ensured by the way MatrixFree reads out the
3216 // indices.
3217 for (unsigned int comp = 0; comp < n_components; ++comp)
3218 {
3219 unsigned int ind_local = 0;
3220
3221 // check whether there is any constraint on the current cell
3222 for (; index_indicators != next_index_indicators;
3223 ++index_indicators)
3224 {
3225 const std::pair<unsigned short, unsigned short> indicator =
3226 dof_info.constraint_indicator[index_indicators];
3227
3228 // run through values up to next constraint
3229 for (unsigned int j = 0; j < indicator.first; ++j)
3230 operation.process_dof(dof_indices[v][j],
3231 *src[0],
3232 values_dofs[comp][ind_local + j][v]);
3233 ind_local += indicator.first;
3234 dof_indices[v] += indicator.first;
3235
3236 // constrained case: build the local value as a linear
3237 // combination of the global value according to constraints
3238 Number value;
3239 operation.pre_constraints(values_dofs[comp][ind_local][v],
3240 value);
3241
3242 const Number *data_val =
3243 this->matrix_free->constraint_pool_begin(indicator.second);
3244 const Number *end_pool =
3245 this->matrix_free->constraint_pool_end(indicator.second);
3246
3247 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3248 operation.process_constraint(*dof_indices[v],
3249 *data_val,
3250 *src[0],
3251 value);
3252
3253 operation.post_constraints(value,
3254 values_dofs[comp][ind_local][v]);
3255 ++ind_local;
3256 }
3257
3258 AssertIndexRange(ind_local, dofs_per_component + 1);
3259
3260 // get the dof values past the last constraint
3261 for (; ind_local < dofs_per_component;
3262 ++dof_indices[v], ++ind_local)
3263 {
3264 AssertIndexRange(*dof_indices[v], src[0]->size());
3265 operation.process_dof(*dof_indices[v],
3266 *src[0],
3267 values_dofs[comp][ind_local][v]);
3268 }
3269
3270 if (apply_constraints == true && comp + 1 < n_components)
3271 next_index_indicators =
3272 dof_info.row_starts[cell_dof_index + comp + 2].second;
3273 }
3274 }
3275 }
3276}
3277
3278
3279
3280template <int dim,
3281 int n_components_,
3282 typename Number,
3283 bool is_face,
3284 typename VectorizedArrayType>
3285template <typename VectorType, typename VectorOperation>
3286inline void
3289 const VectorOperation &operation,
3290 const std::array<VectorType *, n_components_> &src) const
3291{
3292 Assert(!local_dof_indices.empty(), ExcNotInitialized());
3293
3294 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3295 unsigned int index = this->first_selected_component * dofs_per_component;
3296 for (unsigned int comp = 0; comp < n_components; ++comp)
3297 {
3298 for (unsigned int i = 0; i < dofs_per_component; ++i, ++index)
3299 {
3300 operation.process_empty(
3301 this->values_dofs[comp * dofs_per_component + i]);
3302 operation.process_dof_global(
3303 local_dof_indices[this->data->lexicographic_numbering[index]],
3304 *src[0],
3305 this->values_dofs[comp * dofs_per_component + i][0]);
3306 }
3307 }
3308}
3309
3310
3311
3312template <int dim,
3313 int n_components_,
3314 typename Number,
3315 bool is_face,
3316 typename VectorizedArrayType>
3317template <typename VectorType, typename VectorOperation>
3318inline void
3321 const VectorOperation &operation,
3322 const std::array<VectorType *, n_components_> &src,
3323 const std::array<
3325 n_components_> &vectors_sm,
3326 const std::bitset<n_lanes> &mask) const
3327{
3328 // This functions processes the functions read_dof_values,
3329 // distribute_local_to_global, and set_dof_values with the same code for
3330 // contiguous cell indices (DG case). The distinction between these three
3331 // cases is made by the input VectorOperation that either reads values from
3332 // a vector and puts the data into the local data field or write local data
3333 // into the vector. Certain operations are no-ops for the given use case.
3334
3335 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
3336 vector_selector;
3338 is_face ? this->dof_access_index :
3340 const unsigned int n_active_lanes = mask.count();
3341
3342 const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
3343 const std::vector<unsigned int> &dof_indices_cont =
3344 dof_info.dof_indices_contiguous[ind];
3345
3346 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3347 std::array<VectorizedArrayType *, n_components> values_dofs{{nullptr}};
3348 for (unsigned int c = 0; c < n_components; ++c)
3349 values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3350 c * dofs_per_component;
3351
3353
3354 const bool accesses_exterior_dofs =
3355 this->dof_access_index ==
3357 this->is_interior_face() == false;
3358
3359 // Simple case: We have contiguous storage, so we can simply copy out the
3360 // data
3361 if (dof_info.index_storage_variants[ind][this->cell] ==
3363 interleaved_contiguous &&
3364 n_active_lanes == n_lanes && !accesses_exterior_dofs)
3365 {
3366 const unsigned int dof_index =
3367 dof_indices_cont[this->cell * n_lanes] +
3368 this->dof_info
3369 ->component_dof_indices_offset[this->active_fe_index]
3370 [this->first_selected_component] *
3371 n_lanes;
3372 if (n_components == 1 || this->n_fe_components == 1)
3373 for (unsigned int comp = 0; comp < n_components; ++comp)
3374 operation.process_dofs_vectorized(dofs_per_component,
3375 dof_index,
3376 *src[comp],
3377 values_dofs[comp],
3378 vector_selector);
3379 else
3380 operation.process_dofs_vectorized(dofs_per_component * n_components,
3381 dof_index,
3382 *src[0],
3383 values_dofs[0],
3384 vector_selector);
3385 return;
3386 }
3387
3388 const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3389
3390 // More general case: Must go through the components one by one and apply
3391 // some transformations
3392 const unsigned int n_filled_lanes =
3393 dof_info.n_vectorization_lanes_filled[ind][this->cell];
3394
3395 const bool use_vectorized_path = n_filled_lanes == n_lanes &&
3396 n_active_lanes == n_lanes &&
3397 !accesses_exterior_dofs;
3398
3399 if (vectors_sm[0] != nullptr)
3400 {
3401 const auto compute_vector_ptrs = [&](const unsigned int comp) {
3402 std::array<typename VectorType::value_type *, n_lanes> vector_ptrs{
3403 {nullptr}};
3404
3405 const auto upper_bound =
3406 std::min<unsigned int>(n_filled_lanes, n_lanes);
3407 for (unsigned int v = 0; v < upper_bound; ++v)
3408 {
3409 if (mask[v] == false)
3410 {
3411 vector_ptrs[v] = nullptr;
3412 continue;
3413 }
3414
3417 Assert(ind < dof_info.dof_indices_contiguous_sm.size(),
3418 ExcIndexRange(ind,
3419 0,
3420 dof_info.dof_indices_contiguous_sm.size()));
3421 Assert(
3422 cells[v] < dof_info.dof_indices_contiguous_sm[ind].size(),
3423 ExcIndexRange(cells[v],
3424 0,
3425 dof_info.dof_indices_contiguous_sm[ind].size()));
3426
3427 const auto &temp =
3428 dof_info.dof_indices_contiguous_sm[ind][cells[v]];
3429
3430 if (temp.first != numbers::invalid_unsigned_int)
3431 vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
3432 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3434 [this->active_fe_index][this->first_selected_component]);
3435 else
3436 vector_ptrs[v] = nullptr;
3437 }
3438 for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
3439 vector_ptrs[v] = nullptr;
3440
3441 return vector_ptrs;
3442 };
3443
3444 if (use_vectorized_path)
3445 {
3446 if (n_components == 1 || this->n_fe_components == 1)
3447 {
3448 for (unsigned int comp = 0; comp < n_components; ++comp)
3449 {
3450 auto vector_ptrs = compute_vector_ptrs(comp);
3451 operation.process_dofs_vectorized_transpose(
3452 dofs_per_component,
3453 vector_ptrs,
3454 values_dofs[comp],
3455 vector_selector);
3456 }
3457 }
3458 else
3459 {
3460 auto vector_ptrs = compute_vector_ptrs(0);
3461 operation.process_dofs_vectorized_transpose(dofs_per_component *
3462 n_components,
3463 vector_ptrs,
3464 &values_dofs[0][0],
3465 vector_selector);
3466 }
3467 }
3468 else
3469 for (unsigned int comp = 0; comp < n_components; ++comp)
3470 {
3471 auto vector_ptrs = compute_vector_ptrs(
3472 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3473
3474 for (unsigned int i = 0; i < dofs_per_component; ++i)
3475 operation.process_empty(values_dofs[comp][i]);
3476
3477 if (n_components == 1 || this->n_fe_components == 1)
3478 {
3479 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3480 if (mask[v] == true)
3481 for (unsigned int i = 0; i < dofs_per_component; ++i)
3482 operation.process_dof(vector_ptrs[v][i],
3483 values_dofs[comp][i][v]);
3484 }
3485 else
3486 {
3487 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3488 if (mask[v] == true)
3489 for (unsigned int i = 0; i < dofs_per_component; ++i)
3490 operation.process_dof(
3491 vector_ptrs[v][i + comp * dofs_per_component],
3492 values_dofs[comp][i][v]);
3493 }
3494 }
3495 return;
3496 }
3497
3498 std::array<unsigned int, n_lanes> dof_indices{
3500 Assert(n_filled_lanes <= n_lanes, ExcInternalError());
3501 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3502 {
3503 Assert(mask[v] == false || cells[v] != numbers::invalid_unsigned_int,
3505 if (mask[v] == true)
3506 dof_indices[v] =
3507 dof_indices_cont[cells[v]] +
3508 this->dof_info
3509 ->component_dof_indices_offset[this->active_fe_index]
3510 [this->first_selected_component] *
3511 dof_info.dof_indices_interleave_strides[ind][cells[v]];
3512 }
3513
3514 // In the case with contiguous cell indices, we know that there are no
3515 // constraints and that the indices within each element are contiguous
3516 if (use_vectorized_path)
3517 {
3518 if (dof_info.index_storage_variants[ind][this->cell] ==
3520 contiguous)
3521 {
3522 if (n_components == 1 || this->n_fe_components == 1)
3523 for (unsigned int comp = 0; comp < n_components; ++comp)
3524 operation.process_dofs_vectorized_transpose(dofs_per_component,
3525 dof_indices.data(),
3526 *src[comp],
3527 values_dofs[comp],
3528 vector_selector);
3529 else
3530 operation.process_dofs_vectorized_transpose(dofs_per_component *
3531 n_components,
3532 dof_indices.data(),
3533 *src[0],
3534 &values_dofs[0][0],
3535 vector_selector);
3536 }
3537 else if (dof_info.index_storage_variants[ind][this->cell] ==
3539 interleaved_contiguous_strided)
3540 {
3541 std::array<typename VectorType::value_type *, n_components> src_ptrs{
3542 {nullptr}};
3543 if (n_components == 1 || this->n_fe_components == 1)
3544 for (unsigned int comp = 0; comp < n_components; ++comp)
3545 src_ptrs[comp] = const_cast<typename VectorType::value_type *>(
3546 src[comp]->begin());
3547 else
3548 src_ptrs[0] =
3549 const_cast<typename VectorType::value_type *>(src[0]->begin());
3550
3551 if (n_components == 1 || this->n_fe_components == 1)
3552 for (unsigned int i = 0; i < dofs_per_component; ++i)
3553 {
3554 for (unsigned int comp = 0; comp < n_components; ++comp)
3555 operation.process_dof_gather(dof_indices.data(),
3556 *src[comp],
3557 i * n_lanes,
3558 src_ptrs[comp] + i * n_lanes,
3559 values_dofs[comp][i],
3560 vector_selector);
3561 }
3562 else
3563 for (unsigned int comp = 0; comp < n_components; ++comp)
3564 for (unsigned int i = 0; i < dofs_per_component; ++i)
3565 {
3566 operation.process_dof_gather(
3567 dof_indices.data(),
3568 *src[0],
3569 (comp * dofs_per_component + i) * n_lanes,
3570 src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
3571 values_dofs[comp][i],
3572 vector_selector);
3573 }
3574 }
3575 else
3576 {
3577 Assert(dof_info.index_storage_variants[ind][this->cell] ==
3579 IndexStorageVariants::interleaved_contiguous_mixed_strides,
3581 std::array<typename VectorType::value_type *, n_components> src_ptrs{
3582 {nullptr}};
3583 if (n_components == 1 || this->n_fe_components == 1)
3584 for (unsigned int comp = 0; comp < n_components; ++comp)
3585 src_ptrs[comp] = const_cast<typename VectorType::value_type *>(
3586 src[comp]->begin());
3587 else
3588 src_ptrs[0] =
3589 const_cast<typename VectorType::value_type *>(src[0]->begin());
3590
3591 const unsigned int *offsets =
3592 &dof_info.dof_indices_interleave_strides[ind][n_lanes * this->cell];
3593 if (n_components == 1 || this->n_fe_components == 1)
3594 for (unsigned int i = 0; i < dofs_per_component; ++i)
3595 {
3596 for (unsigned int comp = 0; comp < n_components; ++comp)
3597 operation.process_dof_gather(dof_indices.data(),
3598 *src[comp],
3599 0,
3600 src_ptrs[comp],
3601 values_dofs[comp][i],
3602 vector_selector);
3604 for (unsigned int v = 0; v < n_lanes; ++v)
3605 dof_indices[v] += offsets[v];
3606 }
3607 else
3608 for (unsigned int comp = 0; comp < n_components; ++comp)
3609 for (unsigned int i = 0; i < dofs_per_component; ++i)
3610 {
3611 operation.process_dof_gather(dof_indices.data(),
3612 *src[0],
3613 0,
3614 src_ptrs[0],
3615 values_dofs[comp][i],
3616 vector_selector);
3618 for (unsigned int v = 0; v < n_lanes; ++v)
3619 dof_indices[v] += offsets[v];
3620 }
3621 }
3622 }
3623 else
3624 for (unsigned int comp = 0; comp < n_components; ++comp)
3625 {
3626 for (unsigned int i = 0; i < dofs_per_component; ++i)
3627 operation.process_empty(values_dofs[comp][i]);
3628 if (accesses_exterior_dofs)
3629 {
3630 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3631 if (mask[v] == true)
3632 {
3633 if (dof_info.index_storage_variants
3634 [ind][cells[v] / VectorizedArrayType::size()] ==
3637 {
3638 if (n_components == 1 || this->n_fe_components == 1)
3639 {
3640 for (unsigned int i = 0; i < dofs_per_component; ++i)
3641 operation.process_dof(dof_indices[v] + i,
3642 *src[comp],
3643 values_dofs[comp][i][v]);
3644 }
3645 else
3646 {
3647 for (unsigned int i = 0; i < dofs_per_component; ++i)
3648 operation.process_dof(dof_indices[v] + i +
3649 comp * dofs_per_component,
3650 *src[0],
3651 values_dofs[comp][i][v]);
3652 }
3653 }
3654 else
3655 {
3656 const unsigned int offset =
3657 dof_info.dof_indices_interleave_strides[ind][cells[v]];
3658 AssertIndexRange(offset, VectorizedArrayType::size() + 1);
3659 if (n_components == 1 || this->n_fe_components == 1)
3660 {
3661 for (unsigned int i = 0; i < dofs_per_component; ++i)
3662 operation.process_dof(dof_indices[v] + i * offset,
3663 *src[comp],
3664 values_dofs[comp][i][v]);
3665 }
3666 else
3667 {
3668 for (unsigned int i = 0; i < dofs_per_component; ++i)
3669 operation.process_dof(
3670 dof_indices[v] +
3671 (i + comp * dofs_per_component) * offset,
3672 *src[0],
3673 values_dofs[comp][i][v]);
3674 }
3675 }
3676 }
3677 }
3678 else
3679 {
3680 if (dof_info.index_storage_variants[ind][this->cell] ==
3682 contiguous)
3683 {
3684 if (n_components == 1 || this->n_fe_components == 1)
3685 {
3686 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3687 if (mask[v] == true)
3688 for (unsigned int i = 0; i < dofs_per_component; ++i)
3689 operation.process_dof(dof_indices[v] + i,
3690 *src[comp],
3691 values_dofs[comp][i][v]);
3692 }
3693 else
3694 {
3695 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3696 if (mask[v] == true)
3697 for (unsigned int i = 0; i < dofs_per_component; ++i)
3698 operation.process_dof(dof_indices[v] + i +
3699 comp * dofs_per_component,
3700 *src[0],
3701 values_dofs[comp][i][v]);
3702 }
3703 }
3704 else
3705 {
3706 const unsigned int *offsets =
3708 [ind][VectorizedArrayType::size() * this->cell];
3709 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3710 AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
3711 if (n_components == 1 || this->n_fe_components == 1)
3712 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3713 {
3714 if (mask[v] == true)
3715 for (unsigned int i = 0; i < dofs_per_component; ++i)
3716 operation.process_dof(dof_indices[v] + i * offsets[v],
3717 *src[comp],
3718 values_dofs[comp][i][v]);
3719 }
3720 else
3721 {
3722 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3723 if (mask[v] == true)
3724 for (unsigned int i = 0; i < dofs_per_component; ++i)
3725 operation.process_dof(
3726 dof_indices[v] +
3727 (i + comp * dofs_per_component) * offsets[v],
3728 *src[0],
3729 values_dofs[comp][i][v]);
3730 }
3731 }
3732 }
3733 }
3734}
3735
3736namespace internal
3737{
3738 template <
3739 typename Number,
3740 typename VectorType,
3741 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType> * = nullptr>
3742 decltype(std::declval<VectorType>().begin())
3743 get_beginning(VectorType &vec)
3744 {
3745 return vec.begin();
3746 }
3747
3748 template <
3749 typename Number,
3750 typename VectorType,
3751 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * = nullptr>
3752 typename VectorType::value_type *
3753 get_beginning(VectorType &)
3754 {
3755 return nullptr;
3756 }
3757
3758 template <typename VectorType,
3759 std::enable_if_t<has_shared_vector_data<VectorType>, VectorType> * =
3760 nullptr>
3761 const std::vector<ArrayView<const typename VectorType::value_type>> *
3762 get_shared_vector_data(VectorType *vec,
3763 const bool is_valid_mode_for_sm,
3764 const unsigned int active_fe_index,
3766 {
3767 // note: no hp is supported
3768 if (is_valid_mode_for_sm &&
3769 dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
3770 .size() > 0 &&
3771 active_fe_index == 0)
3772 return &vec->shared_vector_data();
3773 else
3774 return nullptr;
3775 }
3776
3777 template <typename VectorType,
3778 std::enable_if_t<!has_shared_vector_data<VectorType>, VectorType>
3779 * = nullptr>
3780 const std::vector<ArrayView<const typename VectorType::value_type>> *
3781 get_shared_vector_data(VectorType *,
3782 const bool,
3783 const unsigned int,
3785 {
3786 return nullptr;
3787 }
3788
3789 template <int n_components, typename VectorType>
3790 std::pair<
3791 std::array<typename internal::BlockVectorSelector<
3792 VectorType,
3793 IsBlockVector<VectorType>::value>::BaseVectorType *,
3794 n_components>,
3795 std::array<
3796 const std::vector<ArrayView<const typename internal::BlockVectorSelector<
3797 VectorType,
3798 IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
3799 n_components>>
3800 get_vector_data(VectorType &src,
3801 const unsigned int first_index,
3802 const bool is_valid_mode_for_sm,
3803 const unsigned int active_fe_index,
3805 {
3806 // select between block vectors and non-block vectors. Note that the number
3807 // of components is checked in the internal data
3808 std::pair<
3809 std::array<typename internal::BlockVectorSelector<
3810 VectorType,
3811 IsBlockVector<VectorType>::value>::BaseVectorType *,
3812 n_components>,
3813 std::array<
3814 const std::vector<
3815 ArrayView<const typename internal::BlockVectorSelector<
3816 VectorType,
3817 IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
3818 n_components>>
3819 src_data;
3820
3821 for (unsigned int d = 0; d < n_components; ++d)
3822 src_data.first[d] = internal::BlockVectorSelector<
3823 VectorType,
3824 IsBlockVector<VectorType>::value>::get_vector_component(src,
3825 d +
3826 first_index);
3827
3828 for (unsigned int d = 0; d < n_components; ++d)
3829 src_data.second[d] = get_shared_vector_data(
3830 const_cast<typename internal::BlockVectorSelector<
3831 std::remove_const_t<VectorType>,
3833 *>(src_data.first[d]),
3834 is_valid_mode_for_sm,
3835 active_fe_index,
3836 dof_info);
3837
3838 return src_data;
3839 }
3840} // namespace internal
3841
3842
3843
3844template <int dim,
3845 int n_components_,
3846 typename Number,
3847 bool is_face,
3848 typename VectorizedArrayType>
3849inline void
3852{
3853 if (this->dof_info == nullptr ||
3854 this->dof_info->hanging_node_constraint_masks.empty() ||
3855 this->dof_info->hanging_node_constraint_masks_comp.empty() ||
3856 this->dof_info->hanging_node_constraint_masks_comp
3857 [this->active_fe_index][this->first_selected_component] == false)
3858 return; // nothing to do with faces
3859
3860 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
3861 constraint_mask{{internal::MatrixFreeFunctions::
3862 unconstrained_compressed_constraint_kind}};
3863
3864 bool hn_available = false;
3865
3866 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3867
3868 for (unsigned int v = 0; v < n_lanes; ++v)
3869 {
3870 if (cells[v] == numbers::invalid_unsigned_int)
3871 {
3872 constraint_mask[v] = internal::MatrixFreeFunctions::
3874 continue;
3875 }
3876
3877 const unsigned int cell_index = cells[v];
3878 const auto mask =
3880 constraint_mask[v] = mask;
3881
3882 hn_available |= (mask != internal::MatrixFreeFunctions::
3884 }
3885
3886 if (hn_available == false)
3887 return; // no hanging node on cell batch -> nothing to do
3888
3890 apply(n_components,
3891 this->data->data.front().fe_degree,
3892 this->get_shape_info(),
3893 transpose,
3894 constraint_mask,
3895 this->values_dofs);
3896}
3897
3898
3899
3900template <int dim,
3901 int n_components_,
3902 typename Number,
3903 bool is_face,
3904 typename VectorizedArrayType>
3905template <typename VectorType>
3906inline void
3908 read_dof_values(const VectorType &src,
3909 const unsigned int first_index,
3910 const std::bitset<n_lanes> &mask)
3911{
3912 const auto src_data = internal::get_vector_data<n_components_>(
3913 src,
3914 first_index,
3915 this->dof_info != nullptr &&
3916 this->dof_access_index ==
3918 this->active_fe_index,
3919 this->dof_info);
3920
3922 read_write_operation(reader, src_data.first, src_data.second, mask, true);
3923
3924 apply_hanging_node_constraints(false);
3925
3926 if constexpr (running_in_debug_mode())
3927 {
3928 this->dof_values_initialized = true;
3929 }
3930}
3931
3932
3933
3934template <int dim,
3935 int n_components_,
3936 typename Number,
3937 bool is_face,
3938 typename VectorizedArrayType>
3939template <typename VectorType>
3940inline void
3942 read_dof_values_plain(const VectorType &src,
3943 const unsigned int first_index,
3944 const std::bitset<n_lanes> &mask)
3945{
3946 const auto src_data = internal::get_vector_data<n_components_>(
3947 src,
3948 first_index,
3949 this->dof_access_index ==
3951 this->active_fe_index,
3952 this->dof_info);
3953
3955 read_write_operation(reader, src_data.first, src_data.second, mask, false);
3956
3957 if constexpr (running_in_debug_mode())
3958 {
3959 this->dof_values_initialized = true;
3960 }
3961}
3962
3963
3964
3965template <int dim,
3966 int n_components_,
3967 typename Number,
3968 bool is_face,
3969 typename VectorizedArrayType>
3970template <typename VectorType>
3971inline void
3973 distribute_local_to_global(VectorType &dst,
3974 const unsigned int first_index,
3975 const std::bitset<n_lanes> &mask) const
3976{
3977 if constexpr (running_in_debug_mode())
3978 {
3979 Assert(this->dof_values_initialized == true,
3981 }
3982
3983 apply_hanging_node_constraints(true);
3984
3985 const auto dst_data = internal::get_vector_data<n_components_>(
3986 dst,
3987 first_index,
3988 this->dof_access_index ==
3990 this->active_fe_index,
3991 this->dof_info);
3992
3994 distributor;
3995 read_write_operation(distributor, dst_data.first, dst_data.second, mask);
3996}
3997
3998
3999
4000template <int dim,
4001 int n_components_,
4002 typename Number,
4003 bool is_face,
4004 typename VectorizedArrayType>
4005template <typename VectorType>
4006inline void
4008 set_dof_values(VectorType &dst,
4009 const unsigned int first_index,
4010 const std::bitset<n_lanes> &mask) const
4011{
4012 if constexpr (running_in_debug_mode())
4013 {
4014 Assert(this->dof_values_initialized == true,
4016 }
4017
4018 const auto dst_data = internal::get_vector_data<n_components_>(
4019 dst,
4020 first_index,
4021 this->dof_access_index ==
4023 this->active_fe_index,
4024 this->dof_info);
4025
4027 read_write_operation(setter, dst_data.first, dst_data.second, mask);
4028}
4029
4030
4031
4032template <int dim,
4033 int n_components_,
4034 typename Number,
4035 bool is_face,
4036 typename VectorizedArrayType>
4037template <typename VectorType>
4038inline void
4040 set_dof_values_plain(VectorType &dst,
4041 const unsigned int first_index,
4042 const std::bitset<n_lanes> &mask) const
4043{
4044 if constexpr (running_in_debug_mode())
4045 {
4046 Assert(this->dof_values_initialized == true,
4048 }
4049
4050 const auto dst_data = internal::get_vector_data<n_components_>(
4051 dst,
4052 first_index,
4053 this->dof_access_index ==
4055 this->active_fe_index,
4056 this->dof_info);
4057
4059 read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
4060}
4061
4062
4063
4064/*------------------------------ access to data fields ----------------------*/
4065
4066
4067
4068template <int dim,
4069 int n_components_,
4070 typename Number,
4071 bool is_face,
4072 typename VectorizedArrayType>
4074 typename FEEvaluationBase<dim,
4075 n_components_,
4076 Number,
4077 is_face,
4078 VectorizedArrayType>::value_type
4080 get_dof_value(const unsigned int dof) const
4081{
4082 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4083 if constexpr (n_components == 1)
4084 return this->values_dofs[dof];
4085 else
4086 {
4087 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4088 Tensor<1, n_components_, VectorizedArrayType> return_value;
4089 for (unsigned int comp = 0; comp < n_components; ++comp)
4090 return_value[comp] = this->values_dofs[comp * dofs + dof];
4091 return return_value;
4092 }
4093}
4094
4095
4096
4097template <int dim,
4098 int n_components_,
4099 typename Number,
4100 bool is_face,
4101 typename VectorizedArrayType>
4103 typename FEEvaluationBase<dim,
4104 n_components_,
4105 Number,
4106 is_face,
4107 VectorizedArrayType>::value_type
4109 get_value(const unsigned int q_point) const
4110{
4111 if constexpr (running_in_debug_mode())
4112 {
4113 Assert(this->values_quad_initialized == true,
4115 }
4116
4117 AssertIndexRange(q_point, this->n_quadrature_points);
4118 if constexpr (n_components == 1)
4119 return this->values_quad[q_point];
4120 else
4121 {
4122 if (n_components == dim &&
4123 this->data->element_type ==
4125 {
4126 // Piola transform is required
4127 if constexpr (running_in_debug_mode())
4128 {
4129 Assert(this->values_quad_initialized == true,
4131 }
4132
4133 AssertIndexRange(q_point, this->n_quadrature_points);
4134 Assert(this->J_value != nullptr,
4136 "update_values"));
4137 const std::size_t nqp = this->n_quadrature_points;
4139
4140 if (!is_face &&
4142 {
4143 // Cartesian cell
4144 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
4145 const VectorizedArrayType inv_det =
4146 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4147 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4148 this->jacobian[0][2][2];
4149
4150 // J * u * det(J^-1)
4151 for (unsigned int comp = 0; comp < n_components; ++comp)
4152 value_out[comp] = this->values_quad[comp * nqp + q_point] *
4153 jac[comp][comp] * inv_det;
4154 }
4155 else
4156 {
4157 // Affine or general cell
4158 const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
4159 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
4160 this->jacobian[q_point] :
4161 this->jacobian[0];
4163 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
4164 transpose(invert(inv_t_jac)) :
4165 this->jacobian[1];
4166
4167 // Derivatives are reordered for faces. Need to take this into
4168 // account
4169 const VectorizedArrayType inv_det =
4170 (is_face && dim == 2 && this->get_face_no() < 2) ?
4171 -determinant(inv_t_jac) :
4172 determinant(inv_t_jac);
4173 // J * u * det(J^-1)
4174 for (unsigned int comp = 0; comp < n_components; ++comp)
4175 {
4176 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
4177 for (unsigned int e = 1; e < dim; ++e)
4178 value_out[comp] +=
4179 this->values_quad[e * nqp + q_point] * jac[comp][e];
4180 value_out[comp] *= inv_det;
4181 }
4182 }
4183 return value_out;
4184 }
4185 else
4186 {
4187 const std::size_t nqp = this->n_quadrature_points;
4189 for (unsigned int comp = 0; comp < n_components; ++comp)
4190 return_value[comp] = this->values_quad[comp * nqp + q_point];
4191 return return_value;
4192 }
4193 }
4194}
4195
4196
4197
4198template <int dim,
4199 int n_components_,
4200 typename Number,
4201 bool is_face,
4202 typename VectorizedArrayType>
4204 typename FEEvaluationBase<dim,
4205 n_components_,
4206 Number,
4207 is_face,
4208 VectorizedArrayType>::gradient_type
4210 get_gradient(const unsigned int q_point) const
4211{
4212 if constexpr (running_in_debug_mode())
4213 {
4214 Assert(this->gradients_quad_initialized == true,
4216 }
4217
4218 AssertIndexRange(q_point, this->n_quadrature_points);
4219 Assert(this->jacobian != nullptr,
4221 "update_gradients"));
4222 const std::size_t nqp = this->n_quadrature_points;
4223
4224 if constexpr (n_components == dim && dim > 1)
4225 {
4226 if (this->data->element_type ==
4228 {
4229 // Piola transform is required
4230 if constexpr (running_in_debug_mode())
4231 {
4232 Assert(this->gradients_quad_initialized == true,
4234 }
4235
4236 AssertIndexRange(q_point, this->n_quadrature_points);
4237 Assert(this->jacobian != nullptr,
4239 "update_gradients"));
4240 const std::size_t nqp = this->n_quadrature_points;
4241 const std::size_t nqp_d = nqp * dim;
4243 const VectorizedArrayType *gradients =
4244 this->gradients_quad + q_point * dim;
4245
4246
4247 if (!is_face &&
4249 {
4250 // Cartesian cell
4251 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
4252 this->jacobian[0];
4254 this->jacobian[1];
4255 const VectorizedArrayType inv_det =
4256 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4257 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4258 this->jacobian[0][2][2];
4259
4260 // J * grad_quad * J^-1 * det(J^-1)
4261 for (unsigned int d = 0; d < dim; ++d)
4262 for (unsigned int comp = 0; comp < n_components; ++comp)
4263 grad_out[comp][d] = gradients[comp * nqp_d + d] *
4264 inv_t_jac[d][d] *
4265 (jac[comp][comp] * inv_det);
4266 }
4267 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4268 {
4269 // Affine cell
4270 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
4271 this->jacobian[0];
4273 this->jacobian[1];
4274
4275 // Derivatives are reordered for faces. Need to take this into
4276 // account
4277 const VectorizedArrayType inv_det =
4278 (is_face && dim == 2 && this->get_face_no() < 2) ?
4279 -determinant(inv_t_jac) :
4280 determinant(inv_t_jac);
4281
4282 VectorizedArrayType tmp[dim][dim];
4283 // J * grad_quad * J^-1 * det(J^-1)
4284 for (unsigned int d = 0; d < dim; ++d)
4285 for (unsigned int e = 0; e < dim; ++e)
4286 {
4287 tmp[d][e] = inv_t_jac[d][0] * gradients[e * nqp_d + 0];
4288 for (unsigned int f = 1; f < dim; ++f)
4289 tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
4290 }
4291 for (unsigned int comp = 0; comp < n_components; ++comp)
4292 for (unsigned int d = 0; d < dim; ++d)
4293 {
4294 VectorizedArrayType res = jac[comp][0] * tmp[d][0];
4295 for (unsigned int f = 1; f < dim; ++f)
4296 res += jac[comp][f] * tmp[d][f];
4297
4298 grad_out[comp][d] = res * inv_det;
4299 }
4300 }
4301 else
4302 {
4303 // General cell
4304
4305 // This assert could be removed if we make sure that this is
4306 // updated even though update_hessians or update_jacobian_grads is
4307 // not passed, i.e make the necessary changes in
4308 // MatrixFreeFunctions::MappingInfoStorage::compute_update_flags
4309 Assert(this->jacobian_gradients_non_inverse != nullptr,
4311 "update_hessians"));
4312
4313 const auto jac_grad =
4314 this->jacobian_gradients_non_inverse[q_point];
4315 const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
4316 this->jacobian[q_point];
4317
4318 // Derivatives are reordered for faces. Need to take this into
4319 // account
4320 const VectorizedArrayType inv_det =
4321 (is_face && dim == 2 && this->get_face_no() < 2) ?
4322 -determinant(inv_t_jac) :
4323 determinant(inv_t_jac);
4325 invert(inv_t_jac);
4326
4327 // (J * grad_quad) * J^-1 * det(J^-1), part in braces
4328 VectorizedArrayType tmp[dim][dim];
4329 for (unsigned int d = 0; d < dim; ++d)
4330 for (unsigned int e = 0; e < dim; ++e)
4331 {
4332 tmp[e][d] = t_jac[0][d] * gradients[0 * nqp_d + e];
4333 for (unsigned int f = 1; f < dim; ++f)
4334 tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
4335 }
4336
4337 // Add (jac_grad * values) * J^{-1} * det(J^{-1}), combine terms
4338 // outside braces with gradient part from above
4339 for (unsigned int d = 0; d < dim; ++d)
4340 {
4341 for (unsigned int e = 0; e < dim; ++e)
4342 tmp[e][d] +=
4343 jac_grad[e][d] * this->values_quad[e * nqp + q_point];
4344 for (unsigned int f = 0, r = dim; f < dim; ++f)
4345 for (unsigned int k = f + 1; k < dim; ++k, ++r)
4346 {
4347 tmp[k][d] +=
4348 jac_grad[r][d] * this->values_quad[f * nqp + q_point];
4349 tmp[f][d] +=
4350 jac_grad[r][d] * this->values_quad[k * nqp + q_point];
4351 }
4352 }
4353
4354 // Apply J^{-1} appearing in both terms outside braces above
4355 for (unsigned int d = 0; d < dim; ++d)
4356 for (unsigned int e = 0; e < dim; ++e)
4357 {
4358 VectorizedArrayType res = tmp[0][d] * inv_t_jac[e][0];
4359 for (unsigned int f = 1; f < dim; ++f)
4360 res += tmp[f][d] * inv_t_jac[e][f];
4361 grad_out[d][e] = res;
4362 }
4363
4364 // Add -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1})),
4365 // which can be expressed as a rank-1 update tmp[d] * tmp4[e],
4366 // where tmp = J * values and tmp4 = (J^{-T} * jac_grad * J^{-1})
4367 VectorizedArrayType tmp3[dim], tmp4[dim];
4368 for (unsigned int d = 0; d < dim; ++d)
4369 {
4370 tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0];
4371 for (unsigned int e = 1; e < dim; ++e)
4372 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
4373 }
4374 for (unsigned int e = 0, k = dim; e < dim; ++e)
4375 for (unsigned int f = e + 1; f < dim; ++k, ++f)
4376 for (unsigned int d = 0; d < dim; ++d)
4377 {
4378 tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d];
4379 tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d];
4380 }
4381 for (unsigned int d = 0; d < dim; ++d)
4382 {
4383 tmp4[d] = tmp3[0] * inv_t_jac[d][0];
4384 for (unsigned int e = 1; e < dim; ++e)
4385 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
4386 }
4387
4388 VectorizedArrayType tmp2[dim];
4389 for (unsigned int d = 0; d < dim; ++d)
4390 {
4391 tmp2[d] = t_jac[0][d] * this->values_quad[q_point];
4392 for (unsigned e = 1; e < dim; ++e)
4393 tmp2[d] +=
4394 t_jac[e][d] * this->values_quad[e * nqp + q_point];
4395 }
4396
4397 for (unsigned int d = 0; d < dim; ++d)
4398 for (unsigned int e = 0; e < dim; ++e)
4399 {
4400 grad_out[d][e] -= tmp4[e] * tmp2[d];
4401
4402 // finally multiply by det(J^{-1}) necessary in all
4403 // contributions above
4404 grad_out[d][e] *= inv_det;
4405 }
4406 }
4407 return grad_out;
4408 }
4409 }
4411
4412 // Cartesian cell
4413 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4414 {
4415 for (unsigned int comp = 0; comp < n_components; ++comp)
4416 for (unsigned int d = 0; d < dim; ++d)
4417 grad_out[comp][d] =
4418 this->gradients_quad[(comp * nqp + q_point) * dim + d] *
4419 this->jacobian[0][d][d];
4420 }
4421 // cell with general/affine Jacobian
4422 else
4423 {
4425 this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
4426 q_point :
4427 0];
4428 for (unsigned int comp = 0; comp < n_components; ++comp)
4429 for (unsigned int d = 0; d < dim; ++d)
4430 {
4431 grad_out[comp][d] =
4432 jac[d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4433 for (unsigned int e = 1; e < dim; ++e)
4434 grad_out[comp][d] +=
4435 jac[d][e] *
4436 this->gradients_quad[(comp * nqp + q_point) * dim + e];
4437 }
4438 }
4439 if constexpr (n_components == 1)
4440 return grad_out[0];
4441 else
4442 return grad_out;
4443}
4444
4445
4446
4447template <int dim,
4448 int n_components_,
4449 typename Number,
4450 bool is_face,
4451 typename VectorizedArrayType>
4453 typename FEEvaluationBase<dim,
4454 n_components_,
4455 Number,
4456 is_face,
4457 VectorizedArrayType>::value_type
4459 get_normal_derivative(const unsigned int q_point) const
4460{
4461 AssertIndexRange(q_point, this->n_quadrature_points);
4462 if constexpr (running_in_debug_mode())
4463 {
4464 Assert(this->gradients_quad_initialized == true,
4466 }
4467
4468 Assert(this->normal_x_jacobian != nullptr,
4470 "update_gradients"));
4471
4472 const std::size_t nqp = this->n_quadrature_points;
4474
4475 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4476 for (unsigned int comp = 0; comp < n_components; ++comp)
4477 grad_out[comp] =
4478 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4479 (this->normal_x_jacobian[0][dim - 1]);
4480 else
4481 {
4482 const std::size_t index =
4483 this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
4484 for (unsigned int comp = 0; comp < n_components; ++comp)
4485 {
4486 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4487 this->normal_x_jacobian[index][0];
4488 for (unsigned int d = 1; d < dim; ++d)
4489 grad_out[comp] +=
4490 this->gradients_quad[(comp * nqp + q_point) * dim + d] *
4491 this->normal_x_jacobian[index][d];
4492 }
4493 }
4494 if constexpr (n_components == 1)
4495 return grad_out[0];
4496 else
4497 return grad_out;
4498}
4499
4500
4501
4502namespace internal
4503{
4504 // compute tmp = hess_unit(u) * J^T. do this manually because we do not
4505 // store the lower diagonal because of symmetry
4506 template <typename VectorizedArrayType>
4507 inline void
4508 hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
4509 const VectorizedArrayType *const hessians,
4510 const unsigned int,
4511 VectorizedArrayType (&tmp)[1][1])
4512 {
4513 tmp[0][0] = jac[0][0] * hessians[0];
4514 }
4515
4516 template <typename VectorizedArrayType>
4517 inline void
4518 hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
4519 const VectorizedArrayType *const hessians,
4520 const unsigned int nqp,
4521 VectorizedArrayType (&tmp)[2][2])
4522 {
4523 for (unsigned int d = 0; d < 2; ++d)
4524 {
4525 tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
4526 tmp[1][d] =
4527 (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
4528 }
4529 }
4530
4531 template <typename VectorizedArrayType>
4532 inline void
4533 hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
4534 const VectorizedArrayType *const hessians,
4535 const unsigned int nqp,
4536 VectorizedArrayType (&tmp)[3][3])
4537 {
4538 for (unsigned int d = 0; d < 3; ++d)
4539 {
4540 tmp[0][d] =
4541 (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
4542 jac[d][2] * hessians[4 * nqp]);
4543 tmp[1][d] =
4544 (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
4545 jac[d][2] * hessians[5 * nqp]);
4546 tmp[2][d] =
4547 (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
4548 jac[d][2] * hessians[2 * nqp]);
4549 }
4550 }
4551} // namespace internal
4552
4553
4554
4555template <int dim,
4556 int n_components_,
4557 typename Number,
4558 bool is_face,
4559 typename VectorizedArrayType>
4560inline typename FEEvaluationBase<dim,
4561 n_components_,
4562 Number,
4563 is_face,
4564 VectorizedArrayType>::hessian_type
4566 get_hessian(const unsigned int q_point) const
4567{
4568 if constexpr (running_in_debug_mode())
4569 {
4570 Assert(this->hessians_quad_initialized == true,
4572 }
4573 AssertIndexRange(q_point, this->n_quadrature_points);
4574
4575 Assert(this->jacobian != nullptr,
4577 "update_hessian"));
4579 this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4580 0 :
4581 q_point];
4582
4584
4585 const std::size_t nqp = this->n_quadrature_points;
4586 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4587
4588 // Cartesian cell
4589 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4590 {
4591 for (unsigned int comp = 0; comp < n_components; ++comp)
4592 {
4593 for (unsigned int d = 0; d < dim; ++d)
4594 hessian_out[comp][d][d] =
4595 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4596 (jac[d][d] * jac[d][d]);
4597 switch (dim)
4598 {
4599 case 1:
4600 break;
4601 case 2:
4602 hessian_out[comp][0][1] =
4603 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4604 (jac[0][0] * jac[1][1]);
4605 break;
4606 case 3:
4607 hessian_out[comp][0][1] =
4608 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4609 (jac[0][0] * jac[1][1]);
4610 hessian_out[comp][0][2] =
4611 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4612 (jac[0][0] * jac[2][2]);
4613 hessian_out[comp][1][2] =
4614 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4615 (jac[1][1] * jac[2][2]);
4616 break;
4617 default:
4619 }
4620 for (unsigned int d = 0; d < dim; ++d)
4621 for (unsigned int e = d + 1; e < dim; ++e)
4622 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4623 }
4624 }
4625 // cell with general Jacobian, but constant within the cell
4626 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4627 {
4628 for (unsigned int comp = 0; comp < n_components; ++comp)
4629 {
4630 VectorizedArrayType tmp[dim][dim];
4631 internal::hessian_unit_times_jac(
4632 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4633
4634 // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4635 for (unsigned int d = 0; d < dim; ++d)
4636 for (unsigned int e = d; e < dim; ++e)
4637 {
4638 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4639 for (unsigned int f = 1; f < dim; ++f)
4640 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4641 }
4642
4643 // no J' * grad(u) part here because the Jacobian is constant
4644 // throughout the cell and hence, its derivative is zero
4645
4646 // take symmetric part
4647 for (unsigned int d = 0; d < dim; ++d)
4648 for (unsigned int e = d + 1; e < dim; ++e)
4649 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4650 }
4651 }
4652 // cell with general Jacobian
4653 else
4654 {
4655 const auto &jac_grad = this->jacobian_gradients[q_point];
4656 for (unsigned int comp = 0; comp < n_components; ++comp)
4657 {
4658 VectorizedArrayType tmp[dim][dim];
4659 internal::hessian_unit_times_jac(
4660 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4661
4662 // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4663 for (unsigned int d = 0; d < dim; ++d)
4664 for (unsigned int e = d; e < dim; ++e)
4665 {
4666 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4667 for (unsigned int f = 1; f < dim; ++f)
4668 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4669 }
4670
4671 // add diagonal part of J' * grad(u)
4672 for (unsigned int d = 0; d < dim; ++d)
4673 for (unsigned int e = 0; e < dim; ++e)
4674 hessian_out[comp][d][d] +=
4675 jac_grad[d][e] *
4676 this->gradients_quad[(comp * nqp + q_point) * dim + e];
4677
4678 // add off-diagonal part of J' * grad(u)
4679 for (unsigned int d = 0, count = dim; d < dim; ++d)
4680 for (unsigned int e = d + 1; e < dim; ++e, ++count)
4681 for (unsigned int f = 0; f < dim; ++f)
4682 hessian_out[comp][d][e] +=
4683 jac_grad[count][f] *
4684 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4685
4686 // take symmetric part
4687 for (unsigned int d = 0; d < dim; ++d)
4688 for (unsigned int e = d + 1; e < dim; ++e)
4689 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4690 }
4691 }
4692 if constexpr (n_components == 1)
4693 return hessian_out[0];
4694 else
4695 return hessian_out;
4696}
4697
4698
4699
4700template <int dim,
4701 int n_components_,
4702 typename Number,
4703 bool is_face,
4704 typename VectorizedArrayType>
4705inline typename FEEvaluationBase<dim,
4706 n_components_,
4707 Number,
4708 is_face,
4709 VectorizedArrayType>::gradient_type
4711 get_hessian_diagonal(const unsigned int q_point) const
4712{
4713 Assert(!is_face, ExcNotImplemented());
4714 if constexpr (running_in_debug_mode())
4715 {
4716 Assert(this->hessians_quad_initialized == true,
4718 }
4719 AssertIndexRange(q_point, this->n_quadrature_points);
4720
4721 Assert(this->jacobian != nullptr, ExcNotImplemented());
4723 this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4724 0 :
4725 q_point];
4726
4727 const std::size_t nqp = this->n_quadrature_points;
4728 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4730
4731 // Cartesian cell
4732 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4733 {
4734 for (unsigned int comp = 0; comp < n_components; ++comp)
4735 for (unsigned int d = 0; d < dim; ++d)
4736 hessian_out[comp][d] =
4737 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4738 (jac[d][d] * jac[d][d]);
4739 }
4740 // cell with general Jacobian, but constant within the cell
4741 else if (this->cell_type == internal::MatrixFreeFunctions::affine)
4742 {
4743 for (unsigned int comp = 0; comp < n_components; ++comp)
4744 {
4745 // compute laplacian before the gradient because it needs to access
4746 // unscaled gradient data
4747 VectorizedArrayType tmp[dim][dim];
4748 internal::hessian_unit_times_jac(
4749 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4750
4751 // compute only the trace part of hessian, J * tmp = J *
4752 // hess_unit(u) * J^T
4753 for (unsigned int d = 0; d < dim; ++d)
4754 {
4755 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4756 for (unsigned int f = 1; f < dim; ++f)
4757 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4758 }
4759 }
4760 }
4761 // cell with general Jacobian
4762 else
4763 {
4764 const auto &jac_grad = this->jacobian_gradients[q_point];
4765 for (unsigned int comp = 0; comp < n_components; ++comp)
4766 {
4767 // compute laplacian before the gradient because it needs to access
4768 // unscaled gradient data
4769 VectorizedArrayType tmp[dim][dim];
4770 internal::hessian_unit_times_jac(
4771 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4772
4773 // compute only the trace part of hessian, J * tmp = J *
4774 // hess_unit(u) * J^T
4775 for (unsigned int d = 0; d < dim; ++d)
4776 {
4777 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4778 for (unsigned int f = 1; f < dim; ++f)
4779 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4780 }
4781
4782 for (unsigned int d = 0; d < dim; ++d)
4783 for (unsigned int e = 0; e < dim; ++e)
4784 hessian_out[comp][d] +=
4785 jac_grad[d][e] *
4786 this->gradients_quad[(comp * nqp + q_point) * dim + e];
4787 }
4788 }
4789
4790 if constexpr (n_components == 1)
4791 return hessian_out[0];
4792 else
4793 return hessian_out;
4794}
4795
4796
4797
4798template <int dim,
4799 int n_components_,
4800 typename Number,
4801 bool is_face,
4802 typename VectorizedArrayType>
4803inline typename FEEvaluationBase<dim,
4804 n_components_,
4805 Number,
4806 is_face,
4807 VectorizedArrayType>::value_type
4809 get_laplacian(const unsigned int q_point) const
4810{
4811 Assert(is_face == false, ExcNotImplemented());
4812 if constexpr (running_in_debug_mode())
4813 {
4814 Assert(this->hessians_quad_initialized == true,
4816 }
4817 AssertIndexRange(q_point, this->n_quadrature_points);
4818
4819 const gradient_type hess_diag = get_hessian_diagonal(q_point);
4820 if constexpr (n_components == 1)
4821 {
4822 VectorizedArrayType sum = hess_diag[0];
4823 for (unsigned int d = 1; d < dim; ++d)
4824 sum += hess_diag[d];
4825 return sum;
4826 }
4827 else
4828 {
4830 for (unsigned int comp = 0; comp < n_components; ++comp)
4831 {
4832 laplacian_out[comp] = hess_diag[comp][0];
4833 for (unsigned int d = 1; d < dim; ++d)
4834 laplacian_out[comp] += hess_diag[comp][d];
4835 }
4836 return laplacian_out;
4837 }
4838}
4839
4840
4841
4842template <int dim,
4843 int n_components_,
4844 typename Number,
4845 bool is_face,
4846 typename VectorizedArrayType>
4847inline typename FEEvaluationBase<dim,
4848 n_components_,
4849 Number,
4850 is_face,
4851 VectorizedArrayType>::value_type
4853 get_normal_hessian(const unsigned int q_point) const
4854{
4855 if constexpr (running_in_debug_mode())
4856 {
4857 Assert(this->hessians_quad_initialized == true,
4859 }
4860 AssertIndexRange(q_point, this->n_quadrature_points);
4861
4862 Assert(this->normal_x_jacobian != nullptr,
4864 "update_hessians"));
4865
4867
4868 const std::size_t nqp = this->n_quadrature_points;
4869 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4870
4871 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4872 {
4873 const auto nxj = this->normal_x_jacobian[0];
4874
4875 for (unsigned int comp = 0; comp < n_components; ++comp)
4876 {
4877 for (unsigned int d = 0; d < dim; ++d)
4878 hessian_out[comp] +=
4879 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4880 (nxj[d]) * (nxj[d]);
4881
4882 switch (dim)
4883 {
4884 case 1:
4885 break;
4886 case 2:
4887 hessian_out[comp] +=
4888 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4889 (nxj[0] * nxj[1]);
4890 break;
4891 case 3:
4892 hessian_out[comp] +=
4893 2. * this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4894 (nxj[0] * nxj[1]);
4895 hessian_out[comp] +=
4896 2. * this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4897 (nxj[0] * nxj[2]);
4898 hessian_out[comp] +=
4899 2. * this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4900 (nxj[1] * nxj[2]);
4901 break;
4902 default:
4904 }
4905 }
4906 }
4907 // cell with general Jacobian
4908 else
4909 {
4910 const auto normal = this->normal_vector(q_point);
4911 const auto hessian = get_hessian(q_point);
4912
4913 if constexpr (n_components == 1)
4914 hessian_out[0] = hessian * normal * normal;
4915 else
4916 for (unsigned int comp = 0; comp < n_components; ++comp)
4917 hessian_out[comp] = hessian[comp] * normal * normal;
4918 }
4919 if constexpr (n_components == 1)
4920 return hessian_out[0];
4921 else
4922 return hessian_out;
4923}
4924
4925
4926
4927template <int dim,
4928 int n_components_,
4929 typename Number,
4930 bool is_face,
4931 typename VectorizedArrayType>
4932inline DEAL_II_ALWAYS_INLINE void
4934 submit_dof_value(const value_type val_in, const unsigned int dof)
4935{
4936 if constexpr (running_in_debug_mode())
4937 {
4938 this->dof_values_initialized = true;
4939 }
4940 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4941 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4942 for (unsigned int comp = 0; comp < n_components; ++comp)
4943 if constexpr (n_components == 1)
4944 this->values_dofs[comp * dofs + dof] = val_in;
4945 else
4946 this->values_dofs[comp * dofs + dof] = val_in[comp];
4947}
4948
4949
4950
4951template <int dim,
4952 int n_components_,
4953 typename Number,
4954 bool is_face,
4955 typename VectorizedArrayType>
4956inline DEAL_II_ALWAYS_INLINE void
4958 submit_value(const value_type val_in, const unsigned int q_point)
4959{
4960 if constexpr (running_in_debug_mode())
4961 {
4962 Assert(this->is_reinitialized, ExcNotInitialized());
4963 }
4964 AssertIndexRange(q_point, this->n_quadrature_points);
4965 Assert(this->J_value != nullptr,
4967 "update_values"));
4968 if constexpr (running_in_debug_mode())
4969 {
4970 this->values_quad_submitted = true;
4971 }
4972
4973 const std::size_t nqp = this->n_quadrature_points;
4974 VectorizedArrayType *values = this->values_quad + q_point;
4975
4976 const VectorizedArrayType JxW =
4977 this->cell_type <= internal::MatrixFreeFunctions::affine ?
4978 this->J_value[0] * this->quadrature_weights[q_point] :
4979 this->J_value[q_point];
4980 if constexpr (n_components == 1)
4981 values[0] = val_in * JxW;
4982 else
4983 {
4984 if (n_components == dim &&
4985 this->data->element_type ==
4987 {
4988 // Piola transform is required
4989 AssertIndexRange(q_point, this->n_quadrature_points);
4990 Assert(this->J_value != nullptr,
4992 "update_value"));
4993 if constexpr (running_in_debug_mode())
4994 {
4995 Assert(this->is_reinitialized, ExcNotInitialized());
4996 this->values_quad_submitted = true;
4997 }
4998
4999 VectorizedArrayType *values = this->values_quad + q_point;
5000 const std::size_t nqp = this->n_quadrature_points;
5001
5002 if (!is_face &&
5004 {
5005 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
5006 const VectorizedArrayType weight =
5007 this->quadrature_weights[q_point];
5008
5009 for (unsigned int comp = 0; comp < n_components; ++comp)
5010 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
5011 }
5012 else
5013 {
5014 // Affine or general cell
5015 const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
5016 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5017 this->jacobian[q_point] :
5018 this->jacobian[0];
5019
5020 // Derivatives are reordered for faces. Need to take this into
5021 // account and 1/inv_det != J_value for faces
5022 const VectorizedArrayType fac =
5023 (!is_face) ?
5024 this->quadrature_weights[q_point] :
5025 (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
5026 this->J_value[q_point] :
5027 this->J_value[0] * this->quadrature_weights[q_point]) *
5028 ((dim == 2 && this->get_face_no() < 2) ?
5029 -determinant(inv_t_jac) :
5030 determinant(inv_t_jac)));
5032 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5033 transpose(invert(inv_t_jac)) :
5034 this->jacobian[1];
5035
5036 // J^T * u * factor
5037 for (unsigned int comp = 0; comp < n_components; ++comp)
5038 {
5039 values[comp * nqp] = val_in[0] * jac[0][comp];
5040 for (unsigned int e = 1; e < dim; ++e)
5041 values[comp * nqp] += val_in[e] * jac[e][comp];
5042 values[comp * nqp] *= fac;
5043 }
5044 }
5045 }
5046 else
5047 for (unsigned int comp = 0; comp < n_components; ++comp)
5048 values[comp * nqp] = val_in[comp] * JxW;
5049 }
5050}
5051
5052
5053
5054template <int dim,
5055 int n_components_,
5056 typename Number,
5057 bool is_face,
5058 typename VectorizedArrayType>
5059template <int, typename>
5060inline DEAL_II_ALWAYS_INLINE void
5063 const unsigned int q_point)
5064{
5065 static_assert(n_components == 1,
5066 "Do not try to modify the default template parameters used for"
5067 " selectively enabling this function via std::enable_if!");
5068 submit_value(val_in[0], q_point);
5069}
5070
5071
5072
5073template <int dim,
5074 int n_components_,
5075 typename Number,
5076 bool is_face,
5077 typename VectorizedArrayType>
5078inline DEAL_II_ALWAYS_INLINE void
5080 submit_gradient(const gradient_type grad_in, const unsigned int q_point)
5081{
5082 if constexpr (running_in_debug_mode())
5083 {
5084 Assert(this->is_reinitialized, ExcNotInitialized());
5085 }
5086 AssertIndexRange(q_point, this->n_quadrature_points);
5087 Assert(this->J_value != nullptr,
5089 "update_gradients"));
5090 Assert(this->jacobian != nullptr,
5092 "update_gradients"));
5093 if constexpr (running_in_debug_mode())
5094 {
5095 this->gradients_quad_submitted = true;
5096 }
5097
5098 if constexpr (dim > 1 && n_components == dim)
5099 {
5100 if (this->data->element_type ==
5102 {
5103 // Piola transform is required
5104
5105 if constexpr (running_in_debug_mode())
5106 {
5107 Assert(this->is_reinitialized, ExcNotInitialized());
5108 }
5109 AssertIndexRange(q_point, this->n_quadrature_points);
5110 Assert(this->J_value != nullptr,
5112 "update_gradients"));
5113 Assert(this->jacobian != nullptr,
5115 "update_gradients"));
5116 if constexpr (running_in_debug_mode())
5117 {
5118 this->gradients_quad_submitted = true;
5119 }
5120
5121 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5122 VectorizedArrayType *values =
5123 this->values_from_gradients_quad + q_point;
5124 const std::size_t nqp = this->n_quadrature_points;
5125 const std::size_t nqp_d = nqp * dim;
5126
5127 if (!is_face &&
5129 {
5130 // Cartesian cell
5131 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5132 this->jacobian[0];
5134 this->jacobian[1];
5135 const VectorizedArrayType weight =
5136 this->quadrature_weights[q_point];
5137 for (unsigned int d = 0; d < dim; ++d)
5138 for (unsigned int comp = 0; comp < n_components; ++comp)
5139 gradients[comp * nqp_d + d] = grad_in[comp][d] *
5140 inv_t_jac[d][d] *
5141 (jac[comp][comp] * weight);
5142 }
5143 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5144 {
5145 // Affine cell
5146 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5147 this->jacobian[0];
5149 this->jacobian[1];
5150
5151 // Derivatives are reordered for faces. Need to take this into
5152 // account and 1/inv_det != J_value for faces
5153 const VectorizedArrayType fac =
5154 (!is_face) ?
5155 this->quadrature_weights[q_point] :
5156 this->J_value[0] * this->quadrature_weights[q_point] *
5157 ((dim == 2 && this->get_face_no() < 2) ?
5158 -determinant(inv_t_jac) :
5159 determinant(inv_t_jac));
5160
5161 // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
5162 VectorizedArrayType tmp[dim][dim];
5163 for (unsigned int d = 0; d < dim; ++d)
5164 for (unsigned int e = 0; e < dim; ++e)
5165 {
5166 tmp[d][e] = inv_t_jac[0][d] * grad_in[e][0];
5167 for (unsigned int f = 1; f < dim; ++f)
5168 tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
5169 }
5170 for (unsigned int comp = 0; comp < n_components; ++comp)
5171 for (unsigned int d = 0; d < dim; ++d)
5172 {
5173 VectorizedArrayType res = jac[0][comp] * tmp[d][0];
5174 for (unsigned int f = 1; f < dim; ++f)
5175 res += jac[f][comp] * tmp[d][f];
5176
5177 gradients[comp * nqp_d + d] = res * fac;
5178 }
5179 }
5180 else
5181 {
5182 // General cell
5183
5184 const auto jac_grad =
5185 this->jacobian_gradients_non_inverse[q_point];
5186 const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
5187 this->jacobian[q_point];
5188
5189 // Derivatives are reordered for faces. Need to take this into
5190 // account and 1/inv_det != J_value for faces
5191 const VectorizedArrayType fac =
5192 (!is_face) ? this->quadrature_weights[q_point] :
5193 this->J_value[q_point] *
5194 ((dim == 2 && this->get_face_no() < 2) ?
5195 -determinant(inv_t_jac) :
5196 determinant(inv_t_jac));
5198 invert(inv_t_jac);
5199
5200 // Start evaluation for values part below to enable the compiler
5201 // to possibly re-use the same computation in get_gradient()
5202 // without interfering with stores to 'gradients'
5203 VectorizedArrayType tmp3[dim], tmp4[dim];
5204 for (unsigned int d = 0; d < dim; ++d)
5205 {
5206 tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0];
5207 for (unsigned int e = 1; e < dim; ++e)
5208 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
5209 }
5210 for (unsigned int e = 0, k = dim; e < dim; ++e)
5211 for (unsigned int f = e + 1; f < dim; ++k, ++f)
5212 for (unsigned int d = 0; d < dim; ++d)
5213 {
5214 tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d];
5215 tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d];
5216 }
5217 for (unsigned int d = 0; d < dim; ++d)
5218 {
5219 tmp4[d] = tmp3[0] * inv_t_jac[d][0];
5220 for (unsigned int e = 1; e < dim; ++e)
5221 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
5222 }
5223
5224 const Tensor<2, dim, VectorizedArrayType> grad_in_scaled =
5225 fac * grad_in;
5226
5227 VectorizedArrayType tmp[dim][dim];
5228
5229 // J * (J^{-1} * (grad_in * factor))
5230 for (unsigned int d = 0; d < dim; ++d)
5231 for (unsigned int e = 0; e < dim; ++e)
5232 {
5233 tmp[d][e] = inv_t_jac[0][d] * grad_in_scaled[e][0];
5234 for (unsigned int f = 1; f < dim; ++f)
5235 tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
5236 }
5237
5238 for (unsigned int d = 0; d < dim; ++d)
5239 for (unsigned int e = 0; e < dim; ++e)
5240 {
5241 VectorizedArrayType res = t_jac[d][0] * tmp[e][0];
5242 for (unsigned int f = 1; f < dim; ++f)
5243 res += t_jac[d][f] * tmp[e][f];
5244
5245 gradients[d * nqp_d + e] = res;
5246 }
5247
5248 // jac_grad * (J^{-1} * (grad_in * factor)), re-use part in braces
5249 // as 'tmp' from above
5250 VectorizedArrayType value[dim];
5251 for (unsigned int d = 0; d < dim; ++d)
5252 {
5253 value[d] = tmp[d][0] * jac_grad[d][0];
5254 for (unsigned int e = 1; e < dim; ++e)
5255 value[d] += tmp[d][e] * jac_grad[d][e];
5256 }
5257 for (unsigned int e = 0, k = dim; e < dim; ++e)
5258 for (unsigned int f = e + 1; f < dim; ++k, ++f)
5259 for (unsigned int d = 0; d < dim; ++d)
5260 {
5261 value[e] += tmp[f][d] * jac_grad[k][d];
5262 value[f] += tmp[e][d] * jac_grad[k][d];
5263 }
5264
5265 // -(grad_in * factor) * J * (J^{-T} * jac_grad * J^{-1})
5266 // = -(grad_in * factor) * J * ( \------- tmp4 ---------/ )
5267 for (unsigned int d = 0; d < dim; ++d)
5268 {
5269 VectorizedArrayType tmp2 = grad_in_scaled[d][0] * tmp4[0];
5270 for (unsigned int e = 1; e < dim; ++e)
5271 tmp2 += grad_in_scaled[d][e] * tmp4[e];
5272 for (unsigned int e = 0; e < dim; ++e)
5273 value[e] -= t_jac[e][d] * tmp2;
5274 }
5275
5276 for (unsigned int d = 0; d < dim; ++d)
5277 values[d * nqp] = value[d];
5278 }
5279 return;
5280 }
5281 }
5282
5283 const std::size_t nqp_d = this->n_quadrature_points * dim;
5284 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5285
5286 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5287 {
5288 const VectorizedArrayType JxW =
5289 this->J_value[0] * this->quadrature_weights[q_point];
5290
5291 // Load all entries before starting to write back to make sure the
5292 // compiler sees opportunity of loads in a possibly nearby
5293 // get_gradient() function (i.e., the compiler should not think that
5294 // 'jacobian' could alias with 'gradients').
5295 std::array<VectorizedArrayType, dim> jac;
5296 for (unsigned int d = 0; d < dim; ++d)
5297 jac[d] = this->jacobian[0][d][d];
5298
5299 for (unsigned int d = 0; d < dim; ++d)
5300 {
5301 const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW;
5302 if constexpr (n_components == 1)
5303 gradients[d] = grad_in[d] * factor;
5304 else
5305 for (unsigned int comp = 0; comp < n_components; ++comp)
5306 gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5307 }
5308 }
5309 else
5310 {
5312 this->cell_type > internal::MatrixFreeFunctions::affine ?
5313 this->jacobian[q_point] :
5314 this->jacobian[0];
5315 const VectorizedArrayType JxW =
5316 this->cell_type > internal::MatrixFreeFunctions::affine ?
5317 this->J_value[q_point] :
5318 this->J_value[0] * this->quadrature_weights[q_point];
5319 if constexpr (n_components == 1)
5320 for (unsigned int d = 0; d < dim; ++d)
5321 {
5322 VectorizedArrayType new_val = jac[0][d] * grad_in[0];
5323 for (unsigned int e = 1; e < dim; ++e)
5324 new_val += (jac[e][d] * grad_in[e]);
5325 gradients[d] = new_val * JxW;
5326 }
5327 else
5328 for (unsigned int comp = 0; comp < n_components; ++comp)
5329 for (unsigned int d = 0; d < dim; ++d)
5330 {
5331 VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
5332 for (unsigned int e = 1; e < dim; ++e)
5333 new_val += (jac[e][d] * grad_in[comp][e]);
5334 gradients[comp * nqp_d + d] = new_val * JxW;
5335 }
5336 }
5337}
5338
5339
5340
5341template <int dim,
5342 int n_components_,
5343 typename Number,
5344 bool is_face,
5345 typename VectorizedArrayType>
5346template <int, typename>
5347inline DEAL_II_ALWAYS_INLINE void
5350 const unsigned int q_point)
5351{
5352 static_assert(n_components == 1 && dim == 1,
5353 "Do not try to modify the default template parameters used for"
5354 " selectively enabling this function via std::enable_if!");
5355 submit_gradient(grad_in[0], q_point);
5356}
5357
5358
5359
5360template <int dim,
5361 int n_components_,
5362 typename Number,
5363 bool is_face,
5364 typename VectorizedArrayType>
5365inline DEAL_II_ALWAYS_INLINE void
5367 submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
5368{
5369 AssertIndexRange(q_point, this->n_quadrature_points);
5370 Assert(this->normal_x_jacobian != nullptr,
5372 "update_gradients"));
5373 if constexpr (running_in_debug_mode())
5374 {
5375 this->gradients_quad_submitted = true;
5376 }
5377
5378 const std::size_t nqp_d = this->n_quadrature_points * dim;
5379 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5380
5381 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5382 {
5383 const VectorizedArrayType JxW_jac = this->J_value[0] *
5384 this->quadrature_weights[q_point] *
5385 this->normal_x_jacobian[0][dim - 1];
5386 for (unsigned int comp = 0; comp < n_components; ++comp)
5387 {
5388 for (unsigned int d = 0; d < dim - 1; ++d)
5389 gradients[comp * nqp_d + d] = VectorizedArrayType();
5390 if constexpr (n_components == 1)
5391 gradients[dim - 1] = grad_in * JxW_jac;
5392 else
5393 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5394 }
5395 }
5396 else
5397 {
5398 const unsigned int index =
5399 this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5401 this->normal_x_jacobian[index];
5402 const VectorizedArrayType JxW =
5403 (this->cell_type <= internal::MatrixFreeFunctions::affine) ?
5404 this->J_value[index] * this->quadrature_weights[q_point] :
5405 this->J_value[index];
5406 for (unsigned int comp = 0; comp < n_components; ++comp)
5407 for (unsigned int d = 0; d < dim; ++d)
5408 if constexpr (n_components == 1)
5409 gradients[d] = (grad_in * JxW) * jac[d];
5410 else
5411 gradients[comp * nqp_d + d] = (grad_in[comp] * JxW) * jac[d];
5412 }
5413}
5414
5415
5416
5417template <int dim,
5418 int n_components_,
5419 typename Number,
5420 bool is_face,
5421 typename VectorizedArrayType>
5422inline DEAL_II_ALWAYS_INLINE void
5424 submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
5425{
5426 if constexpr (running_in_debug_mode())
5427 {
5428 Assert(this->is_reinitialized, ExcNotInitialized());
5429 }
5430 AssertIndexRange(q_point, this->n_quadrature_points);
5431 Assert(this->J_value != nullptr,
5433 "update_hessians"));
5434 Assert(this->jacobian != nullptr,
5436 "update_hessians"));
5437 if constexpr (running_in_debug_mode())
5438 {
5439 this->hessians_quad_submitted = true;
5440 }
5441
5442 // compute hessian_unit = J^T * hessian_in(u) * J
5443 const std::size_t nqp = this->n_quadrature_points;
5444 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5445 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5446 {
5447 const VectorizedArrayType JxW =
5448 this->J_value[0] * this->quadrature_weights[q_point];
5449
5450 // diagonal part
5451 for (unsigned int d = 0; d < dim; ++d)
5452 {
5453 const auto jac_d = this->jacobian[0][d][d];
5454 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5455 for (unsigned int comp = 0; comp < n_components; ++comp)
5456 if constexpr (n_components == 1)
5457 this->hessians_quad[d * nqp + q_point] =
5458 hessian_in[d][d] * factor;
5459 else
5460 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5461 hessian_in[comp][d][d] * factor;
5462 }
5463
5464 // off diagonal part
5465 for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5466 for (unsigned int e = 0; e < d; ++e, ++off_dia)
5467 {
5468 const auto jac_d = this->jacobian[0][d][d];
5469 const auto jac_e = this->jacobian[0][e][e];
5470 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5471 for (unsigned int comp = 0; comp < n_components; ++comp)
5472 if constexpr (n_components == 1)
5473 this->hessians_quad[off_dia * nqp + q_point] =
5474 (hessian_in[d][e] + hessian_in[e][d]) * factor;
5475 else
5476 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5477 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5478 }
5479 }
5480 // cell with general Jacobian, but constant within the cell
5481 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5482 {
5483 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
5484 const VectorizedArrayType JxW =
5485 this->J_value[0] * this->quadrature_weights[q_point];
5486 for (unsigned int comp = 0; comp < n_components; ++comp)
5487 {
5489 if constexpr (n_components == 1)
5490 hessian_c = hessian_in;
5491 else
5492 hessian_c = hessian_in[comp];
5493
5494 // 1. tmp = hessian(u) * J
5495 VectorizedArrayType tmp[dim][dim];
5496 for (unsigned int i = 0; i < dim; ++i)
5497 for (unsigned int j = 0; j < dim; ++j)
5498 {
5499 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5500 for (unsigned int k = 1; k < dim; ++k)
5501 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5502 }
5503
5504 // 2. hessian_unit = J^T * tmp
5505 VectorizedArrayType tmp2[dim][dim];
5506 for (unsigned int i = 0; i < dim; ++i)
5507 for (unsigned int j = 0; j < dim; ++j)
5508 {
5509 tmp2[i][j] = jac[0][i] * tmp[0][j];
5510 for (unsigned int k = 1; k < dim; ++k)
5511 tmp2[i][j] += jac[k][i] * tmp[k][j];
5512 }
5513
5514 // diagonal part
5515 for (unsigned int d = 0; d < dim; ++d)
5516 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5517 tmp2[d][d] * JxW;
5518
5519 // off diagonal part
5520 for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5521 for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5522 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5523 (tmp2[d][e] + tmp2[e][d]) * JxW;
5524 }
5525 }
5526 else
5527 {
5528 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
5529 const VectorizedArrayType JxW = this->J_value[q_point];
5530 const auto &jac_grad = this->jacobian_gradients[q_point];
5531 for (unsigned int comp = 0; comp < n_components; ++comp)
5532 {
5534 if constexpr (n_components == 1)
5535 hessian_c = hessian_in;
5536 else
5537 hessian_c = hessian_in[comp];
5538
5539 // 1. tmp = hessian(u) * J
5540 VectorizedArrayType tmp[dim][dim];
5541 for (unsigned int i = 0; i < dim; ++i)
5542 for (unsigned int j = 0; j < dim; ++j)
5543 {
5544 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5545 for (unsigned int k = 1; k < dim; ++k)
5546 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5547 }
5548
5549 // 2. hessian_unit = J^T * tmp
5550 VectorizedArrayType tmp2[dim][dim];
5551 for (unsigned int i = 0; i < dim; ++i)
5552 for (unsigned int j = 0; j < dim; ++j)
5553 {
5554 tmp2[i][j] = jac[0][i] * tmp[0][j];
5555 for (unsigned int k = 1; k < dim; ++k)
5556 tmp2[i][j] += jac[k][i] * tmp[k][j];
5557 }
5558
5559 // diagonal part
5560 for (unsigned int d = 0; d < dim; ++d)
5561 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5562 tmp2[d][d] * JxW;
5563
5564 // off diagonal part
5565 for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5566 for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5567 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5568 (tmp2[d][e] + tmp2[e][d]) * JxW;
5569
5570 // 3. gradient_unit = J' * hessian
5571 for (unsigned int d = 0; d < dim; ++d)
5572 {
5573 VectorizedArrayType sum = 0;
5574 for (unsigned int e = 0; e < dim; ++e)
5575 sum += hessian_c[e][e] * jac_grad[e][d];
5576 for (unsigned int e = 0, count = dim; e < dim; ++e)
5577 for (unsigned int f = e + 1; f < dim; ++f, ++count)
5578 sum +=
5579 (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][d];
5580 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5581 d] = sum * JxW;
5582 }
5583 }
5584 }
5585}
5586
5587
5588
5589template <int dim,
5590 int n_components_,
5591 typename Number,
5592 bool is_face,
5593 typename VectorizedArrayType>
5594inline DEAL_II_ALWAYS_INLINE void
5596 submit_normal_hessian(const value_type normal_hessian_in,
5597 const unsigned int q_point)
5598{
5599 if constexpr (running_in_debug_mode())
5600 {
5601 Assert(this->is_reinitialized, ExcNotInitialized());
5602 }
5603 AssertIndexRange(q_point, this->n_quadrature_points);
5604 Assert(this->J_value != nullptr,
5606 "update_hessians"));
5607 Assert(this->jacobian != nullptr,
5609 "update_hessians"));
5610 if constexpr (running_in_debug_mode())
5611 {
5612 this->hessians_quad_submitted = true;
5613 }
5614
5615 // compute hessian_unit = J^T * hessian_in(u) * J
5616 const std::size_t nqp = this->n_quadrature_points;
5617 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5618 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5619 {
5620 const VectorizedArrayType JxW =
5621 this->J_value[0] * this->quadrature_weights[q_point];
5622
5623 const auto nxj = this->normal_x_jacobian[0];
5624
5625 // diagonal part
5626 for (unsigned int d = 0; d < dim; ++d)
5627 {
5628 const auto nxj_d = nxj[d];
5629 const VectorizedArrayType factor = nxj_d * nxj_d * JxW;
5630 for (unsigned int comp = 0; comp < n_components; ++comp)
5631 if constexpr (n_components == 1)
5632 this->hessians_quad[d * nqp + q_point] =
5633 normal_hessian_in * factor;
5634 else
5635 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5636 normal_hessian_in[comp] * factor;
5637 }
5638
5639 // off diagonal part
5640 for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5641 for (unsigned int e = 0; e < d; ++e, ++off_dia)
5642 {
5643 const auto jac_d = nxj[d];
5644 const auto jac_e = nxj[e];
5645 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5646 for (unsigned int comp = 0; comp < n_components; ++comp)
5647 if constexpr (n_components == 1)
5648 this->hessians_quad[off_dia * nqp + q_point] =
5649 2. * normal_hessian_in * factor;
5650 else
5651 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5652 2. * normal_hessian_in[comp] * factor;
5653 }
5654 }
5655 else
5656 {
5657 const auto normal = this->normal_vector(q_point);
5658 const auto normal_projector = outer_product(normal, normal);
5659 if constexpr (n_components == 1)
5660 submit_hessian(normal_hessian_in * normal_projector, q_point);
5661 else
5662 {
5663 hessian_type tmp;
5664 for (unsigned int comp = 0; comp < n_components; ++comp)
5665 tmp[comp] = normal_hessian_in[comp] * normal_projector;
5666 submit_hessian(tmp, q_point);
5667 }
5668 }
5669}
5670
5671
5672
5673template <int dim,
5674 int n_components_,
5675 typename Number,
5676 bool is_face,
5677 typename VectorizedArrayType>
5678inline typename FEEvaluationBase<dim,
5679 n_components_,
5680 Number,
5681 is_face,
5682 VectorizedArrayType>::value_type
5684 integrate_value() const
5685{
5686 if constexpr (running_in_debug_mode())
5687 {
5688 Assert(this->is_reinitialized, ExcNotInitialized());
5689 Assert(this->values_quad_submitted == true,
5691 }
5692
5694 const std::size_t nqp = this->n_quadrature_points;
5695 for (unsigned int q = 0; q < nqp; ++q)
5696 for (unsigned int comp = 0; comp < n_components; ++comp)
5697 return_value[comp] += this->values_quad[comp * nqp + q];
5698 if constexpr (n_components == 1)
5699 return return_value[0];
5700 else
5701 return return_value;
5702}
5703
5704
5705
5706template <int dim,
5707 int n_components_,
5708 typename Number,
5709 bool is_face,
5710 typename VectorizedArrayType>
5711template <int, typename>
5712inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5714 get_divergence(const unsigned int q_point) const
5715{
5716 static_assert(n_components == dim,
5717 "Do not try to modify the default template parameters used for"
5718 " selectively enabling this function via std::enable_if!");
5719
5720 if constexpr (running_in_debug_mode())
5721 {
5722 Assert(this->gradients_quad_initialized == true,
5724 }
5725 AssertIndexRange(q_point, this->n_quadrature_points);
5726 Assert(this->jacobian != nullptr,
5728 "update_gradients"));
5729
5730 VectorizedArrayType divergence;
5731 const std::size_t nqp = this->n_quadrature_points;
5732
5733 if (dim > 1 &&
5734 this->data->element_type ==
5736 {
5737 VectorizedArrayType inv_det =
5738 (!is_face &&
5739 this->cell_type == internal::MatrixFreeFunctions::cartesian) ?
5740 this->jacobian[0][0][0] *
5741 ((dim == 2) ? this->jacobian[0][1][1] :
5742 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
5743 determinant(this->jacobian[this->cell_type >
5744 internal::MatrixFreeFunctions::affine ?
5745 q_point :
5746 0]);
5747
5748 // on faces in 2d, the determinant has the wrong sign due to ordering of
5749 // derivatives
5750 if (is_face && dim == 2 && this->get_face_no() < 2)
5751 inv_det = -inv_det;
5752
5753 // div * det(J^-1)
5754 divergence = this->gradients_quad[q_point * dim];
5755 for (unsigned int d = 1; d < dim; ++d)
5756 divergence += this->gradients_quad[(d * nqp + q_point) * dim + d];
5757 divergence *= inv_det;
5758 }
5759 else
5760 {
5761 if (!is_face &&
5763 {
5764 // Cartesian cell
5765 divergence =
5766 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
5767 for (unsigned int d = 1; d < dim; ++d)
5768 divergence += this->gradients_quad[(d * nqp + q_point) * dim + d] *
5769 this->jacobian[0][d][d];
5770 }
5771 else
5772 {
5773 // cell with general/constant Jacobian
5775 this->cell_type == internal::MatrixFreeFunctions::general ?
5776 this->jacobian[q_point] :
5777 this->jacobian[0];
5778 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
5779 for (unsigned int e = 1; e < dim; ++e)
5780 divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
5781 for (unsigned int d = 1; d < dim; ++d)
5782 for (unsigned int e = 0; e < dim; ++e)
5783 divergence +=
5784 jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim + e];
5785 }
5786 }
5787 return divergence;
5788}
5789
5790
5791
5792template <int dim,
5793 int n_components_,
5794 typename Number,
5795 bool is_face,
5796 typename VectorizedArrayType>
5797template <int, typename>
5800 get_symmetric_gradient(const unsigned int q_point) const
5801{
5802 static_assert(n_components == dim,
5803 "Do not try to modify the default template parameters used for"
5804 " selectively enabling this function via std::enable_if!");
5805
5806 // copy from generic function into dim-specialization function
5807 const auto grad = get_gradient(q_point);
5808 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5809 VectorizedArrayType half = Number(0.5);
5810 for (unsigned int d = 0; d < dim; ++d)
5811 symmetrized[d] = grad[d][d];
5812 switch (dim)
5813 {
5814 case 1:
5815 break;
5816 case 2:
5817 symmetrized[2] = grad[0][1] + grad[1][0];
5818 symmetrized[2] *= half;
5819 break;
5820 case 3:
5821 symmetrized[3] = grad[0][1] + grad[1][0];
5822 symmetrized[3] *= half;
5823 symmetrized[4] = grad[0][2] + grad[2][0];
5824 symmetrized[4] *= half;
5825 symmetrized[5] = grad[1][2] + grad[2][1];
5826 symmetrized[5] *= half;
5827 break;
5828 default:
5830 }
5832}
5833
5834
5835
5836template <int dim,
5837 int n_components_,
5838 typename Number,
5839 bool is_face,
5840 typename VectorizedArrayType>
5841template <int, typename>
5843 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
5845 get_curl(const unsigned int q_point) const
5846{
5847 static_assert(dim > 1 && n_components == dim,
5848 "Do not try to modify the default template parameters used for"
5849 " selectively enabling this function via std::enable_if!");
5850
5851 // copy from generic function into dim-specialization function
5852 const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
5853 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
5854 switch (dim)
5855 {
5856 case 2:
5857 curl[0] = grad[1][0] - grad[0][1];
5858 break;
5859 case 3:
5860 curl[0] = grad[2][1] - grad[1][2];
5861 curl[1] = grad[0][2] - grad[2][0];
5862 curl[2] = grad[1][0] - grad[0][1];
5863 break;
5864 default:
5866 }
5867 return curl;
5868}
5869
5870
5871
5872template <int dim,
5873 int n_components_,
5874 typename Number,
5875 bool is_face,
5876 typename VectorizedArrayType>
5877template <int, typename>
5878inline DEAL_II_ALWAYS_INLINE void
5880 submit_divergence(const VectorizedArrayType div_in,
5881 const unsigned int q_point)
5882{
5883 static_assert(n_components == dim,
5884 "Do not try to modify the default template parameters used for"
5885 " selectively enabling this function via std::enable_if!");
5886
5887 if constexpr (running_in_debug_mode())
5888 {
5889 Assert(this->is_reinitialized, ExcNotInitialized());
5890 }
5891 AssertIndexRange(q_point, this->n_quadrature_points);
5892 Assert(this->J_value != nullptr,
5894 "update_gradients"));
5895 Assert(this->jacobian != nullptr,
5897 "update_gradients"));
5898 if constexpr (running_in_debug_mode())
5899 {
5900 this->gradients_quad_submitted = true;
5901 }
5902
5903 const std::size_t nqp_d = this->n_quadrature_points * dim;
5904 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5905
5906 if (this->data->element_type ==
5908 {
5909 // General cell
5910
5911 // Derivatives are reordered for faces. Need to take this into account
5912 // and 1/inv_det != J_value for faces
5913 const VectorizedArrayType fac =
5914 (!is_face) ?
5915 this->quadrature_weights[q_point] * div_in :
5916 (this->cell_type > internal::MatrixFreeFunctions::affine ?
5917 this->J_value[q_point] :
5918 this->J_value[0] * this->quadrature_weights[q_point]) *
5919 div_in *
5921 this->jacobian[this->cell_type >
5922 internal::MatrixFreeFunctions::affine ?
5923 q_point :
5924 0]) *
5925 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
5926
5927 for (unsigned int d = 0; d < dim; ++d)
5928 {
5929 for (unsigned int e = 0; e < dim; ++e)
5930 gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
5931 }
5932 this->divergence_is_requested = true;
5933 }
5934 else
5935 {
5936 if (!is_face &&
5938 {
5939 const VectorizedArrayType fac =
5940 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
5941 for (unsigned int d = 0; d < dim; ++d)
5942 {
5943 const VectorizedArrayType jac_dd = this->jacobian[0][d][d];
5944 for (unsigned int e = 0; e < dim; ++e)
5945 gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
5946 }
5947 }
5948 else
5949 {
5951 this->cell_type == internal::MatrixFreeFunctions::general ?
5952 this->jacobian[q_point] :
5953 this->jacobian[0];
5954 const VectorizedArrayType fac =
5955 (this->cell_type == internal::MatrixFreeFunctions::general ?
5956 this->J_value[q_point] :
5957 this->J_value[0] * this->quadrature_weights[q_point]) *
5958 div_in;
5959 for (unsigned int d = 0; d < dim; ++d)
5960 {
5961 for (unsigned int e = 0; e < dim; ++e)
5962 gradients[d * nqp_d + e] = jac[d][e] * fac;
5963 }
5964 }
5965 }
5966}
5967
5968
5969
5970template <int dim,
5971 int n_components_,
5972 typename Number,
5973 bool is_face,
5974 typename VectorizedArrayType>
5975template <int, typename>
5976inline DEAL_II_ALWAYS_INLINE void
5980 const unsigned int q_point)
5981{
5982 static_assert(n_components == dim,
5983 "Do not try to modify the default template parameters used for"
5984 " selectively enabling this function via std::enable_if!");
5985
5987 this->data->element_type !=
5990
5991 // could have used base class operator, but that involves some overhead
5992 // which is inefficient. it is nice to have the symmetric tensor because
5993 // that saves some operations
5994 if constexpr (running_in_debug_mode())
5995 {
5996 Assert(this->is_reinitialized, ExcNotInitialized());
5997 }
5998 AssertIndexRange(q_point, this->n_quadrature_points);
5999 Assert(this->J_value != nullptr,
6001 "update_gradients"));
6002 Assert(this->jacobian != nullptr,
6004 "update_gradients"));
6005 if constexpr (running_in_debug_mode())
6006 {
6007 this->gradients_quad_submitted = true;
6008 }
6009
6010 const std::size_t nqp_d = this->n_quadrature_points * dim;
6011 VectorizedArrayType *gradients = this->gradients_quad + dim * q_point;
6012 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6013 {
6014 const VectorizedArrayType JxW =
6015 this->J_value[0] * this->quadrature_weights[q_point];
6016 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
6017 for (unsigned int d = 0; d < dim; ++d)
6018 gradients[d * nqp_d + d] =
6019 (sym_grad.access_raw_entry(d) * JxW * jac[d][d]);
6020 for (unsigned int e = 0, counter = dim; e < dim; ++e)
6021 for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6022 {
6023 const VectorizedArrayType value =
6024 sym_grad.access_raw_entry(counter) * JxW;
6025 gradients[e * nqp_d + d] = value * jac[d][d];
6026 gradients[d * nqp_d + e] = value * jac[e][e];
6027 }
6028 }
6029 // general/affine cell type
6030 else
6031 {
6032 const VectorizedArrayType JxW =
6033 this->cell_type == internal::MatrixFreeFunctions::general ?
6034 this->J_value[q_point] :
6035 this->J_value[0] * this->quadrature_weights[q_point];
6037 this->cell_type == internal::MatrixFreeFunctions::general ?
6038 this->jacobian[q_point] :
6039 this->jacobian[0];
6040 VectorizedArrayType weighted[dim][dim];
6041 for (unsigned int i = 0; i < dim; ++i)
6042 weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6043 for (unsigned int i = 0, counter = dim; i < dim; ++i)
6044 for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6045 {
6046 const VectorizedArrayType value =
6047 sym_grad.access_raw_entry(counter) * JxW;
6048 weighted[i][j] = value;
6049 weighted[j][i] = value;
6050 }
6051 for (unsigned int comp = 0; comp < dim; ++comp)
6052 for (unsigned int d = 0; d < dim; ++d)
6053 {
6054 VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6055 for (unsigned int e = 1; e < dim; ++e)
6056 new_val += jac[e][d] * weighted[comp][e];
6057 gradients[comp * nqp_d + d] = new_val;
6058 }
6059 }
6060}
6061
6062
6063
6064template <int dim,
6065 int n_components_,
6066 typename Number,
6067 bool is_face,
6068 typename VectorizedArrayType>
6069template <int, typename>
6070inline DEAL_II_ALWAYS_INLINE void
6073 const unsigned int q_point)
6074{
6075 static_assert(n_components == dim,
6076 "Do not try to modify the default template parameters used for"
6077 " selectively enabling this function via std::enable_if!");
6078
6080 switch (dim)
6081 {
6082 case 2:
6083 grad[1][0] = curl[0];
6084 grad[0][1] = -curl[0];
6085 break;
6086 case 3:
6087 grad[2][1] = curl[0];
6088 grad[1][2] = -curl[0];
6089 grad[0][2] = curl[1];
6090 grad[2][0] = -curl[1];
6091 grad[1][0] = curl[2];
6092 grad[0][1] = -curl[2];
6093 break;
6094 default:
6096 }
6097 submit_gradient(grad, q_point);
6098}
6099
6100
6101
6102/*-------------------------- FEEvaluation -----------------------------------*/
6103
6104
6105template <int dim,
6106 int fe_degree,
6107 int n_q_points_1d,
6108 int n_components_,
6109 typename Number,
6110 typename VectorizedArrayType>
6111inline FEEvaluation<dim,
6112 fe_degree,
6113 n_q_points_1d,
6114 n_components_,
6115 Number,
6116 VectorizedArrayType>::
6117 FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
6118 const unsigned int fe_no,
6119 const unsigned int quad_no,
6120 const unsigned int first_selected_component,
6121 const unsigned int active_fe_index,
6122 const unsigned int active_quad_index)
6123 : BaseClass(matrix_free,
6124 fe_no,
6125 first_selected_component,
6126 quad_no,
6127 fe_degree,
6128 static_n_q_points,
6129 true /*note: this is not a face*/,
6130 active_fe_index,
6131 active_quad_index,
6132 numbers::invalid_unsigned_int /*face_type*/)
6133 , dofs_per_component(this->data->dofs_per_component_on_cell)
6134 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6135 , n_q_points(this->data->n_q_points)
6136{
6137 check_template_arguments(fe_no, 0);
6138}
6139
6140
6141
6142template <int dim,
6143 int fe_degree,
6144 int n_q_points_1d,
6145 int n_components_,
6146 typename Number,
6147 typename VectorizedArrayType>
6148inline FEEvaluation<dim,
6149 fe_degree,
6150 n_q_points_1d,
6151 n_components_,
6152 Number,
6153 VectorizedArrayType>::
6154 FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
6155 const std::pair<unsigned int, unsigned int> &range,
6156 const unsigned int dof_no,
6157 const unsigned int quad_no,
6158 const unsigned int first_selected_component)
6159 : FEEvaluation(matrix_free,
6160 dof_no,
6161 quad_no,
6162 first_selected_component,
6163 matrix_free.get_cell_active_fe_index(range))
6164{}
6165
6166
6167
6168template <int dim,
6169 int fe_degree,
6170 int n_q_points_1d,
6171 int n_components_,
6172 typename Number,
6173 typename VectorizedArrayType>
6174inline FEEvaluation<dim,
6175 fe_degree,
6176 n_q_points_1d,
6177 n_components_,
6178 Number,
6179 VectorizedArrayType>::
6180 FEEvaluation(const Mapping<dim> &mapping,
6181 const FiniteElement<dim> &fe,
6182 const Quadrature<1> &quadrature,
6183 const UpdateFlags update_flags,
6184 const unsigned int first_selected_component)
6185 : BaseClass(mapping,
6186 fe,
6187 quadrature,
6188 update_flags,
6189 first_selected_component,
6190 nullptr)
6191 , dofs_per_component(this->data->dofs_per_component_on_cell)
6192 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6193 , n_q_points(this->data->n_q_points)
6194{
6195 check_template_arguments(numbers::invalid_unsigned_int, 0);
6196}
6197
6198
6199
6200template <int dim,
6201 int fe_degree,
6202 int n_q_points_1d,
6203 int n_components_,
6204 typename Number,
6205 typename VectorizedArrayType>
6206inline FEEvaluation<dim,
6207 fe_degree,
6208 n_q_points_1d,
6209 n_components_,
6210 Number,
6211 VectorizedArrayType>::
6212 FEEvaluation(const FiniteElement<dim> &fe,
6213 const Quadrature<1> &quadrature,
6214 const UpdateFlags update_flags,
6215 const unsigned int first_selected_component)
6216 : BaseClass(StaticMappingQ1<dim>::mapping,
6217 fe,
6218 quadrature,
6219 update_flags,
6220 first_selected_component,
6221 nullptr)
6222 , dofs_per_component(this->data->dofs_per_component_on_cell)
6223 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6224 , n_q_points(this->data->n_q_points)
6225{
6226 check_template_arguments(numbers::invalid_unsigned_int, 0);
6227}
6228
6229
6230
6231template <int dim,
6232 int fe_degree,
6233 int n_q_points_1d,
6234 int n_components_,
6235 typename Number,
6236 typename VectorizedArrayType>
6237inline FEEvaluation<dim,
6238 fe_degree,
6239 n_q_points_1d,
6240 n_components_,
6241 Number,
6242 VectorizedArrayType>::
6243 FEEvaluation(const FiniteElement<dim> &fe,
6245 const unsigned int first_selected_component)
6246 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6247 fe,
6248 other.mapped_geometry->get_quadrature(),
6249 other.mapped_geometry->get_fe_values().get_update_flags(),
6250 first_selected_component,
6251 &other)
6252 , dofs_per_component(this->data->dofs_per_component_on_cell)
6253 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6254 , n_q_points(this->data->n_q_points)
6255{
6256 check_template_arguments(numbers::invalid_unsigned_int, 0);
6257}
6258
6259
6260
6261template <int dim,
6262 int fe_degree,
6263 int n_q_points_1d,
6264 int n_components_,
6265 typename Number,
6266 typename VectorizedArrayType>
6267inline FEEvaluation<dim,
6268 fe_degree,
6269 n_q_points_1d,
6270 n_components_,
6271 Number,
6272 VectorizedArrayType>::FEEvaluation(const FEEvaluation
6273 &other)
6274 : BaseClass(other)
6275 , dofs_per_component(this->data->dofs_per_component_on_cell)
6276 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6277 , n_q_points(this->data->n_q_points)
6278{
6279 check_template_arguments(numbers::invalid_unsigned_int, 0);
6280}
6281
6282
6283
6284template <int dim,
6285 int fe_degree,
6286 int n_q_points_1d,
6287 int n_components_,
6288 typename Number,
6289 typename VectorizedArrayType>
6290inline FEEvaluation<dim,
6291 fe_degree,
6292 n_q_points_1d,
6293 n_components_,
6294 Number,
6295 VectorizedArrayType> &
6296FEEvaluation<dim,
6297 fe_degree,
6298 n_q_points_1d,
6299 n_components_,
6300 Number,
6301 VectorizedArrayType>::operator=(const FEEvaluation &other)
6302{
6303 BaseClass::operator=(other);
6304 check_template_arguments(numbers::invalid_unsigned_int, 0);
6305 return *this;
6306}
6307
6308
6309
6310template <int dim,
6311 int fe_degree,
6312 int n_q_points_1d,
6313 int n_components_,
6314 typename Number,
6315 typename VectorizedArrayType>
6316inline void
6317FEEvaluation<dim,
6318 fe_degree,
6319 n_q_points_1d,
6320 n_components_,
6321 Number,
6322 VectorizedArrayType>::
6323 check_template_arguments(const unsigned int dof_no,
6324 const unsigned int first_selected_component)
6325{
6326 (void)dof_no;
6327 (void)first_selected_component;
6328
6329 Assert(
6330 this->data->dofs_per_component_on_cell > 0,
6331 ExcMessage(
6332 "There is nothing useful you can do with an FEEvaluation object with "
6333 "FE_Nothing, i.e., without DoFs! If you have passed to "
6334 "MatrixFree::reinit() a collection of finite elements also containing "
6335 "FE_Nothing, please check - before creating FEEvaluation - the category "
6336 "of the current range by calling either "
6337 "MatrixFree::get_cell_range_category(range) or "
6338 "MatrixFree::get_face_range_category(range). The returned category "
6339 "is the index of the active FE, which you can use to exclude "
6340 "FE_Nothing."));
6341
6342 if constexpr (running_in_debug_mode())
6343 {
6344 // print error message when the dimensions do not match. Propose a
6345 // possible fix
6346 if ((static_cast<unsigned int>(fe_degree) !=
6348 static_cast<unsigned int>(fe_degree) !=
6349 this->data->data.front().fe_degree) ||
6350 n_q_points != this->n_quadrature_points)
6351 {
6352 std::string message =
6353 "-------------------------------------------------------\n";
6354 message +=
6355 "Illegal arguments in constructor/wrong template arguments!\n";
6356 message += " Called --> FEEvaluation<dim,";
6357 message += Utilities::int_to_string(fe_degree) + ",";
6358 message += Utilities::int_to_string(n_q_points_1d);
6359 message += "," + Utilities::int_to_string(n_components);
6360 message += ",Number>(data";
6361 if (first_selected_component != numbers::invalid_unsigned_int)
6362 {
6363 message += ", " + Utilities::int_to_string(dof_no) + ", ";
6364 message += Utilities::int_to_string(this->quad_no) + ", ";
6365 message += Utilities::int_to_string(first_selected_component);
6366 }
6367 message += ")\n";
6368
6369 // check whether some other vector component has the correct number of
6370 // points
6371 unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
6372 proposed_fe_comp = numbers::invalid_unsigned_int,
6373 proposed_quad_comp = numbers::invalid_unsigned_int;
6374 if (dof_no != numbers::invalid_unsigned_int)
6375 {
6376 if (static_cast<unsigned int>(fe_degree) ==
6377 this->data->data.front().fe_degree)
6378 {
6379 proposed_dof_comp = dof_no;
6380 proposed_fe_comp = first_selected_component;
6381 }
6382 else
6383 for (unsigned int no = 0;
6384 no < this->matrix_free->n_components();
6385 ++no)
6386 for (unsigned int nf = 0;
6387 nf < this->matrix_free->n_base_elements(no);
6388 ++nf)
6389 if (this->matrix_free
6390 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6391 .data.front()
6392 .fe_degree == static_cast<unsigned int>(fe_degree))
6393 {
6394 proposed_dof_comp = no;
6395 proposed_fe_comp = nf;
6396 break;
6397 }
6398 if (n_q_points ==
6399 this->mapping_data->descriptor[this->active_quad_index]
6400 .n_q_points)
6401 proposed_quad_comp = this->quad_no;
6402 else
6403 for (unsigned int no = 0;
6404 no <
6405 this->matrix_free->get_mapping_info().cell_data.size();
6406 ++no)
6407 if (this->matrix_free->get_mapping_info()
6408 .cell_data[no]
6409 .descriptor[this->active_quad_index]
6410 .n_q_points == n_q_points)
6411 {
6412 proposed_quad_comp = no;
6413 break;
6414 }
6415 }
6416 if (proposed_dof_comp != numbers::invalid_unsigned_int &&
6417 proposed_quad_comp != numbers::invalid_unsigned_int)
6418 {
6419 if (proposed_dof_comp != first_selected_component)
6420 message += "Wrong vector component selection:\n";
6421 else
6422 message += "Wrong quadrature formula selection:\n";
6423 message += " Did you mean FEEvaluation<dim,";
6424 message += Utilities::int_to_string(fe_degree) + ",";
6425 message += Utilities::int_to_string(n_q_points_1d);
6426 message += "," + Utilities::int_to_string(n_components);
6427 message += ",Number>(data";
6428 if (dof_no != numbers::invalid_unsigned_int)
6429 {
6430 message +=
6431 ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
6432 message +=
6433 Utilities::int_to_string(proposed_quad_comp) + ", ";
6434 message += Utilities::int_to_string(proposed_fe_comp);
6435 }
6436 message += ")?\n";
6437 std::string correct_pos;
6438 if (proposed_dof_comp != dof_no)
6439 correct_pos = " ^ ";
6440 else
6441 correct_pos = " ";
6442 if (proposed_quad_comp != this->quad_no)
6443 correct_pos += " ^ ";
6444 else
6445 correct_pos += " ";
6446 if (proposed_fe_comp != first_selected_component)
6447 correct_pos += " ^\n";
6448 else
6449 correct_pos += " \n";
6450 message +=
6451 " " +
6452 correct_pos;
6453 }
6454 // ok, did not find the numbers specified by the template arguments in
6455 // the given list. Suggest correct template arguments
6456 const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
6457 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6458 message += "Wrong template arguments:\n";
6459 message += " Did you mean FEEvaluation<dim,";
6460 message +=
6461 Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
6462 message += Utilities::int_to_string(proposed_n_q_points_1d);
6463 message += "," + Utilities::int_to_string(n_components);
6464 message += ",Number>(data";
6465 if (dof_no != numbers::invalid_unsigned_int)
6466 {
6467 message += ", " + Utilities::int_to_string(dof_no) + ", ";
6468 message += Utilities::int_to_string(this->quad_no);
6469 message +=
6470 ", " + Utilities::int_to_string(first_selected_component);
6471 }
6472 message += ")?\n";
6473 std::string correct_pos;
6474 if (this->data->data.front().fe_degree !=
6475 static_cast<unsigned int>(fe_degree))
6476 correct_pos = " ^";
6477 else
6478 correct_pos = " ";
6479 if (proposed_n_q_points_1d != n_q_points_1d)
6480 correct_pos += " ^\n";
6481 else
6482 correct_pos += " \n";
6483 message += " " + correct_pos;
6484
6485 Assert(static_cast<unsigned int>(fe_degree) ==
6486 this->data->data.front().fe_degree &&
6487 n_q_points == this->n_quadrature_points,
6488 ExcMessage(message));
6489 }
6490 if (dof_no != numbers::invalid_unsigned_int)
6492 n_q_points,
6493 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6494 }
6495}
6496
6497
6498
6499template <int dim,
6500 int fe_degree,
6501 int n_q_points_1d,
6502 int n_components_,
6503 typename Number,
6504 typename VectorizedArrayType>
6505inline void
6506FEEvaluation<dim,
6507 fe_degree,
6508 n_q_points_1d,
6509 n_components_,
6510 Number,
6511 VectorizedArrayType>::reinit(const unsigned int cell_index)
6512{
6513 Assert(this->matrix_free != nullptr,
6514 ExcMessage("FEEvaluation was initialized without a matrix-free object."
6515 " Integer indexing is not possible."));
6516
6517 Assert(this->dof_info != nullptr, ExcNotInitialized());
6518 Assert(this->mapping_data != nullptr, ExcNotInitialized());
6519 this->cell = cell_index;
6520 this->cell_type =
6521 this->matrix_free->get_mapping_info().get_cell_type(cell_index);
6522
6523 const unsigned int offsets =
6524 this->mapping_data->data_index_offsets[cell_index];
6525 this->jacobian = &this->mapping_data->jacobians[0][offsets];
6526 this->J_value = &this->mapping_data->JxW_values[offsets];
6527 if (!this->mapping_data->jacobian_gradients[0].empty())
6528 {
6529 this->jacobian_gradients =
6530 this->mapping_data->jacobian_gradients[0].data() + offsets;
6531 this->jacobian_gradients_non_inverse =
6532 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
6533 }
6534
6535 if (this->matrix_free->n_active_entries_per_cell_batch(this->cell) == n_lanes)
6536 {
6538 for (unsigned int i = 0; i < n_lanes; ++i)
6539 this->cell_ids[i] = cell_index * n_lanes + i;
6540 }
6541 else
6542 {
6543 unsigned int i = 0;
6544 for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell);
6545 ++i)
6546 this->cell_ids[i] = cell_index * n_lanes + i;
6547 for (; i < n_lanes; ++i)
6548 this->cell_ids[i] = numbers::invalid_unsigned_int;
6549 }
6550
6551 if (this->mapping_data->quadrature_points.empty() == false)
6552 this->quadrature_points =
6553 &this->mapping_data->quadrature_points
6554 [this->mapping_data->quadrature_point_offsets[this->cell]];
6555
6556 if constexpr (running_in_debug_mode())
6557 {
6558 this->is_reinitialized = true;
6559 this->dof_values_initialized = false;
6560 this->values_quad_initialized = false;
6561 this->gradients_quad_initialized = false;
6562 this->hessians_quad_initialized = false;
6563 }
6564}
6565
6566
6567
6568template <int dim,
6569 int fe_degree,
6570 int n_q_points_1d,
6571 int n_components_,
6572 typename Number,
6573 typename VectorizedArrayType>
6574inline void
6575FEEvaluation<dim,
6576 fe_degree,
6577 n_q_points_1d,
6578 n_components_,
6579 Number,
6580 VectorizedArrayType>::reinit(const std::array<unsigned int,
6581 n_lanes> &cell_ids)
6582{
6583 Assert(this->dof_info != nullptr, ExcNotInitialized());
6584 Assert(this->mapping_data != nullptr, ExcNotInitialized());
6585
6586 this->cell = numbers::invalid_unsigned_int;
6587 this->cell_ids = cell_ids;
6588
6589 // determine type of cell batch
6591
6592 for (unsigned int v = 0; v < n_lanes; ++v)
6593 {
6594 const unsigned int cell_index = cell_ids[v];
6595
6597 continue;
6598
6599 this->cell_type =
6600 std::max(this->cell_type,
6601 this->matrix_free->get_mapping_info().get_cell_type(
6602 cell_index / n_lanes));
6603 }
6604
6605 // allocate memory for internal data storage
6606 if (this->mapped_geometry == nullptr)
6607 this->mapped_geometry =
6608 std::make_shared<internal::MatrixFreeFunctions::
6609 MappingDataOnTheFly<dim, VectorizedArrayType>>();
6610
6611 auto &mapping_storage = this->mapped_geometry->get_data_storage();
6612
6613 auto &this_jacobian_data = mapping_storage.jacobians[0];
6614 auto &this_J_value_data = mapping_storage.JxW_values;
6615 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
6616 auto &this_jacobian_gradients_non_inverse_data =
6617 mapping_storage.jacobian_gradients_non_inverse[0];
6618 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
6619
6621 {
6622 if (this_jacobian_data.size() != 2)
6623 this_jacobian_data.resize_fast(2);
6624
6625 if (this_J_value_data.size() != 1)
6626 this_J_value_data.resize_fast(1);
6627
6628 const auto &update_flags_cells =
6629 this->matrix_free->get_mapping_info().update_flags_cells;
6630
6631 if (update_flags_cells & update_jacobian_grads &&
6632 this_jacobian_gradients_data.size() != 1)
6633 {
6634 this_jacobian_gradients_data.resize_fast(1);
6635 this_jacobian_gradients_non_inverse_data.resize_fast(1);
6636 }
6637
6638 if (update_flags_cells & update_quadrature_points &&
6639 this_quadrature_points_data.size() != 1)
6640 this_quadrature_points_data.resize_fast(1);
6641 }
6642 else
6643 {
6644 if (this_jacobian_data.size() != this->n_quadrature_points)
6645 this_jacobian_data.resize_fast(this->n_quadrature_points);
6646
6647 if (this_J_value_data.size() != this->n_quadrature_points)
6648 this_J_value_data.resize_fast(this->n_quadrature_points);
6649
6650 const auto &update_flags_cells =
6651 this->matrix_free->get_mapping_info().update_flags_cells;
6652
6653 if (update_flags_cells & update_jacobian_grads &&
6654 this_jacobian_gradients_data.size() != this->n_quadrature_points)
6655 {
6656 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
6657 this_jacobian_gradients_non_inverse_data.resize_fast(
6658 this->n_quadrature_points);
6659 }
6660
6661 if (update_flags_cells & update_quadrature_points &&
6662 this_quadrature_points_data.size() != this->n_quadrature_points)
6663 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
6664 }
6665
6666 // set pointers to internal data storage
6667 this->jacobian = this_jacobian_data.data();
6668 this->J_value = this_J_value_data.data();
6669 this->jacobian_gradients = this_jacobian_gradients_data.data();
6670 this->jacobian_gradients_non_inverse =
6671 this_jacobian_gradients_non_inverse_data.data();
6672 this->quadrature_points = this_quadrature_points_data.data();
6673
6674 // fill internal data storage lane by lane
6675 for (unsigned int v = 0; v < n_lanes; ++v)
6676 {
6677 const unsigned int cell_index = cell_ids[v];
6678
6680 continue;
6681
6682 const unsigned int cell_batch_index = cell_index / n_lanes;
6683 const unsigned int offsets =
6684 this->mapping_data->data_index_offsets[cell_batch_index];
6685 const unsigned int lane = cell_index % n_lanes;
6686
6687 if (this->cell_type <=
6689 {
6690 // case that all cells are Cartesian or affine
6691 for (unsigned int q = 0; q < 2; ++q)
6692 for (unsigned int i = 0; i < dim; ++i)
6693 for (unsigned int j = 0; j < dim; ++j)
6694 this_jacobian_data[q][i][j][v] =
6695 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
6696
6697 const unsigned int q = 0;
6698
6699 this_J_value_data[q][v] =
6700 this->mapping_data->JxW_values[offsets + q][lane];
6701
6702 const auto &update_flags_cells =
6703 this->matrix_free->get_mapping_info().update_flags_cells;
6704
6705 if (update_flags_cells & update_jacobian_grads)
6706 {
6707 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6708 for (unsigned int j = 0; j < dim; ++j)
6709 this_jacobian_gradients_data[q][i][j][v] =
6710 this->mapping_data
6711 ->jacobian_gradients[0][offsets + q][i][j][lane];
6712
6713 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6714 for (unsigned int j = 0; j < dim; ++j)
6715 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6716 this->mapping_data
6717 ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
6718 [lane];
6719 }
6720
6721 if (update_flags_cells & update_quadrature_points)
6722 for (unsigned int i = 0; i < dim; ++i)
6723 this_quadrature_points_data[q][i][v] =
6724 this->mapping_data->quadrature_points
6725 [this->mapping_data
6726 ->quadrature_point_offsets[cell_batch_index] +
6727 q][i][lane];
6728 }
6729 else
6730 {
6731 // general case that at least one cell is not Cartesian or affine
6732 const auto cell_type =
6733 this->matrix_free->get_mapping_info().get_cell_type(
6734 cell_batch_index);
6735
6736 for (unsigned int q = 0; q < this->n_quadrature_points; ++q)
6737 {
6738 const unsigned int q_src =
6739 (cell_type <=
6741 0 :
6742 q;
6743
6744 this_J_value_data[q][v] =
6745 this->mapping_data->JxW_values[offsets + q_src][lane];
6746
6747 for (unsigned int i = 0; i < dim; ++i)
6748 for (unsigned int j = 0; j < dim; ++j)
6749 this_jacobian_data[q][i][j][v] =
6750 this->mapping_data
6751 ->jacobians[0][offsets + q_src][i][j][lane];
6752
6753 const auto &update_flags_cells =
6754 this->matrix_free->get_mapping_info().update_flags_cells;
6755
6756 if (update_flags_cells & update_jacobian_grads)
6757 {
6758 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6759 for (unsigned int j = 0; j < dim; ++j)
6760 this_jacobian_gradients_data[q][i][j][v] =
6761 this->mapping_data
6762 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
6763
6764 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6765 for (unsigned int j = 0; j < dim; ++j)
6766 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6767 this->mapping_data
6768 ->jacobian_gradients_non_inverse[0][offsets + q_src]
6769 [i][j][lane];
6770 }
6771
6772 if (update_flags_cells & update_quadrature_points)
6773 {
6774 if (cell_type <=
6776 {
6777 // affine case: quadrature points are not available but
6778 // have to be computed from the corner point and the
6779 // Jacobian
6781 this->mapping_data->quadrature_points
6782 [this->mapping_data
6783 ->quadrature_point_offsets[cell_batch_index] +
6784 0];
6785
6787 this->mapping_data->jacobians[0][offsets + 1];
6789 for (unsigned int d = 0; d < dim; ++d)
6790 point[d] +=
6791 jac[d][d] *
6792 static_cast<Number>(
6793 this->descriptor->quadrature.point(q)[d]);
6794 else
6795 for (unsigned int d = 0; d < dim; ++d)
6796 for (unsigned int e = 0; e < dim; ++e)
6797 point[d] +=
6798 jac[d][e] *
6799 static_cast<Number>(
6800 this->descriptor->quadrature.point(q)[e]);
6801
6802 for (unsigned int i = 0; i < dim; ++i)
6803 this_quadrature_points_data[q][i][v] = point[i][lane];
6804 }
6805 else
6806 {
6807 // general case: quadrature points are available
6808 for (unsigned int i = 0; i < dim; ++i)
6809 this_quadrature_points_data[q][i][v] =
6810 this->mapping_data->quadrature_points
6811 [this->mapping_data
6812 ->quadrature_point_offsets[cell_batch_index] +
6813 q][i][lane];
6814 }
6815 }
6816 }
6817 }
6818 }
6819
6820 if constexpr (running_in_debug_mode())
6821 {
6822 this->is_reinitialized = true;
6823 this->dof_values_initialized = false;
6824 this->values_quad_initialized = false;
6825 this->gradients_quad_initialized = false;
6826 this->hessians_quad_initialized = false;
6827 }
6828}
6829
6830
6831
6832template <int dim,
6833 int fe_degree,
6834 int n_q_points_1d,
6835 int n_components_,
6836 typename Number,
6837 typename VectorizedArrayType>
6838template <bool level_dof_access>
6839inline void
6840FEEvaluation<dim,
6841 fe_degree,
6842 n_q_points_1d,
6843 n_components_,
6844 Number,
6845 VectorizedArrayType>::
6847{
6848 Assert(this->matrix_free == nullptr,
6849 ExcMessage("Cannot use initialization from cell iterator if "
6850 "initialized from MatrixFree object. Use variant for "
6851 "on the fly computation with arguments as for FEValues "
6852 "instead"));
6853 Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
6854 this->mapped_geometry->reinit(
6855 static_cast<typename Triangulation<dim>::cell_iterator>(cell));
6856 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6857 if (level_dof_access)
6858 cell->get_mg_dof_indices(this->local_dof_indices);
6859 else
6860 cell->get_dof_indices(this->local_dof_indices);
6861
6862 if constexpr (running_in_debug_mode())
6863 {
6864 this->is_reinitialized = true;
6865 }
6866}
6867
6868
6869
6870template <int dim,
6871 int fe_degree,
6872 int n_q_points_1d,
6873 int n_components_,
6874 typename Number,
6875 typename VectorizedArrayType>
6876inline void
6877FEEvaluation<dim,
6878 fe_degree,
6879 n_q_points_1d,
6880 n_components_,
6881 Number,
6882 VectorizedArrayType>::
6883 reinit(const typename Triangulation<dim>::cell_iterator &cell)
6884{
6885 Assert(this->matrix_free == 0,
6886 ExcMessage("Cannot use initialization from cell iterator if "
6887 "initialized from MatrixFree object. Use variant for "
6888 "on the fly computation with arguments as for FEValues "
6889 "instead"));
6890 Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
6891 this->mapped_geometry->reinit(cell);
6892
6893 if constexpr (running_in_debug_mode())
6894 {
6895 this->is_reinitialized = true;
6896 }
6897}
6898
6899
6900
6901template <int dim,
6902 int fe_degree,
6903 int n_q_points_1d,
6904 int n_components_,
6905 typename Number,
6906 typename VectorizedArrayType>
6907inline void
6908FEEvaluation<dim,
6909 fe_degree,
6910 n_q_points_1d,
6911 n_components_,
6912 Number,
6913 VectorizedArrayType>::
6914 evaluate(const EvaluationFlags::EvaluationFlags evaluation_flags)
6915{
6916 if constexpr (running_in_debug_mode())
6917 {
6918 Assert(this->dof_values_initialized == true,
6920 }
6921 evaluate(this->values_dofs, evaluation_flags);
6922}
6923
6924
6925
6926template <int dim,
6927 int fe_degree,
6928 int n_q_points_1d,
6929 int n_components_,
6930 typename Number,
6931 typename VectorizedArrayType>
6932inline void
6933FEEvaluation<dim,
6934 fe_degree,
6935 n_q_points_1d,
6936 n_components_,
6937 Number,
6938 VectorizedArrayType>::
6939 evaluate(const VectorizedArrayType *values_array,
6940 const EvaluationFlags::EvaluationFlags evaluation_flag)
6941{
6942 const bool hessians_on_general_cells =
6943 evaluation_flag & EvaluationFlags::hessians &&
6944 (this->cell_type > internal::MatrixFreeFunctions::affine);
6945 EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
6946 if (hessians_on_general_cells)
6947 evaluation_flag_actual |= EvaluationFlags::gradients;
6948
6949 if (this->data->element_type ==
6951 evaluation_flag & EvaluationFlags::gradients &&
6952 (this->cell_type > internal::MatrixFreeFunctions::affine))
6953 evaluation_flag_actual |= EvaluationFlags::values;
6954
6955 if constexpr (fe_degree > -1)
6956 {
6958 template run<fe_degree, n_q_points_1d>(n_components,
6959 evaluation_flag_actual,
6960 values_array,
6961 *this);
6962 }
6963 else
6964 {
6966 n_components,
6967 evaluation_flag_actual,
6968 const_cast<VectorizedArrayType *>(values_array),
6969 *this);
6970 }
6971
6972 if constexpr (running_in_debug_mode())
6973 {
6974 this->values_quad_initialized =
6975 evaluation_flag_actual & EvaluationFlags::values;
6976 this->gradients_quad_initialized =
6977 evaluation_flag_actual & EvaluationFlags::gradients;
6978 this->hessians_quad_initialized =
6979 evaluation_flag_actual & EvaluationFlags::hessians;
6980 }
6981}
6982
6983
6984namespace internal
6985{
6989 template <typename Number,
6990 typename VectorizedArrayType,
6991 typename VectorType,
6992 typename EvaluatorType,
6993 std::enable_if_t<internal::has_begin<VectorType> &&
6995 VectorType> * = nullptr>
6996 VectorizedArrayType *
6997 check_vector_access_inplace(const EvaluatorType &fe_eval, VectorType &vector)
6998 {
6999 // for user-defined cell batches this functionality is not supported
7000 if (fe_eval.get_current_cell_index() == numbers::invalid_unsigned_int)
7001 return nullptr;
7002
7003 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
7004 const auto &dof_info = fe_eval.get_dof_info();
7005
7006 // If the index storage is interleaved and contiguous and the vector
7007 // storage has the correct alignment, we can directly pass the pointer
7008 // into the vector to the evaluate() and integrate() calls, without
7009 // reading the vector entries into a separate data field. This saves some
7010 // operations.
7011 if (std::is_same_v<typename VectorType::value_type, Number> &&
7012 dof_info.index_storage_variants
7015 interleaved_contiguous &&
7016 reinterpret_cast<std::size_t>(
7017 vector.begin() +
7018 dof_info.dof_indices_contiguous
7019 [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
7020 [cell * VectorizedArrayType::size()]) %
7021 sizeof(VectorizedArrayType) ==
7022 0)
7023 {
7024 return reinterpret_cast<VectorizedArrayType *>(
7025 vector.begin() +
7026 dof_info.dof_indices_contiguous
7028 [cell * VectorizedArrayType::size()] +
7030 [fe_eval.get_active_fe_index()]
7031 [fe_eval.get_first_selected_component()] *
7032 VectorizedArrayType::size());
7033 }
7034 else
7035 return nullptr;
7036 }
7037
7041 template <typename Number,
7042 typename VectorizedArrayType,
7043 typename VectorType,
7044 typename EvaluatorType,
7045 std::enable_if_t<!internal::has_begin<VectorType> ||
7047 VectorType> * = nullptr>
7048 VectorizedArrayType *
7049 check_vector_access_inplace(const EvaluatorType &, VectorType &)
7050 {
7051 return nullptr;
7052 }
7053} // namespace internal
7054
7055
7056
7057template <int dim,
7058 int fe_degree,
7059 int n_q_points_1d,
7060 int n_components_,
7061 typename Number,
7062 typename VectorizedArrayType>
7063template <typename VectorType>
7064inline void
7065FEEvaluation<dim,
7066 fe_degree,
7067 n_q_points_1d,
7068 n_components_,
7069 Number,
7070 VectorizedArrayType>::
7071 gather_evaluate(const VectorType &input_vector,
7072 const EvaluationFlags::EvaluationFlags evaluation_flag)
7073{
7074 const VectorizedArrayType *src_ptr =
7075 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7076 *this, input_vector);
7077 if (src_ptr != nullptr)
7078 evaluate(src_ptr, evaluation_flag);
7079 else
7080 {
7081 this->read_dof_values(input_vector);
7082 evaluate(this->begin_dof_values(), evaluation_flag);
7083 }
7084}
7085
7086
7087
7088template <int dim,
7089 int fe_degree,
7090 int n_q_points_1d,
7091 int n_components_,
7092 typename Number,
7093 typename VectorizedArrayType>
7094inline void
7095FEEvaluation<dim,
7096 fe_degree,
7097 n_q_points_1d,
7098 n_components_,
7099 Number,
7100 VectorizedArrayType>::
7101 integrate(const EvaluationFlags::EvaluationFlags integration_flag)
7102{
7103 integrate(integration_flag, this->values_dofs);
7104
7105 if constexpr (running_in_debug_mode())
7106 {
7107 this->dof_values_initialized = true;
7108 }
7109}
7110
7111
7112
7113template <int dim,
7114 int fe_degree,
7115 int n_q_points_1d,
7116 int n_components_,
7117 typename Number,
7118 typename VectorizedArrayType>
7119inline void
7120FEEvaluation<dim,
7121 fe_degree,
7122 n_q_points_1d,
7123 n_components_,
7124 Number,
7125 VectorizedArrayType>::
7126 integrate(const EvaluationFlags::EvaluationFlags integration_flag,
7127 VectorizedArrayType *values_array,
7128 const bool sum_into_values_array)
7129{
7130 if constexpr (running_in_debug_mode())
7131 {
7132 if (integration_flag & EvaluationFlags::values)
7133 Assert(this->values_quad_submitted == true,
7135 if (integration_flag & EvaluationFlags::gradients)
7136 Assert(this->gradients_quad_submitted == true,
7138 if ((integration_flag & EvaluationFlags::hessians) != 0u)
7139 Assert(this->hessians_quad_submitted == true,
7141 }
7142 Assert(this->matrix_free != nullptr ||
7143 this->mapped_geometry->is_initialized(),
7145
7146 Assert(
7147 (integration_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients |
7149 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7150 "EvaluationFlags::hessians are supported."));
7151
7152 EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
7153 if (integration_flag & EvaluationFlags::hessians &&
7154 (this->cell_type > internal::MatrixFreeFunctions::affine))
7155 {
7156 unsigned int size = n_components * dim * n_q_points;
7157 if ((integration_flag & EvaluationFlags::gradients) != 0u)
7158 {
7159 for (unsigned int i = 0; i < size; ++i)
7160 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7161 }
7162 else
7163 {
7164 for (unsigned int i = 0; i < size; ++i)
7165 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7166 integration_flag_actual |= EvaluationFlags::gradients;
7167 }
7168 }
7169
7170 if (n_components == dim &&
7171 this->data->element_type ==
7173 integration_flag & EvaluationFlags::gradients &&
7174 this->cell_type > internal::MatrixFreeFunctions::affine &&
7175 this->divergence_is_requested == false)
7176 {
7177 unsigned int size = n_components * n_q_points;
7178 if ((integration_flag & EvaluationFlags::values) != 0u)
7179 {
7180 for (unsigned int i = 0; i < size; ++i)
7181 this->values_quad[i] += this->values_from_gradients_quad[i];
7182 }
7183 else
7184 {
7185 for (unsigned int i = 0; i < size; ++i)
7186 this->values_quad[i] = this->values_from_gradients_quad[i];
7187 integration_flag_actual |= EvaluationFlags::values;
7188 }
7189 }
7190
7191 if constexpr (fe_degree > -1)
7192 {
7194 template run<fe_degree, n_q_points_1d>(n_components,
7195 integration_flag_actual,
7196 values_array,
7197 *this,
7198 sum_into_values_array);
7199 }
7200 else
7201 {
7203 n_components,
7204 integration_flag_actual,
7205 values_array,
7206 *this,
7207 sum_into_values_array);
7208 }
7209
7210 if constexpr (running_in_debug_mode())
7211 {
7212 this->dof_values_initialized = true;
7213 }
7214}
7215
7216
7217
7218template <int dim,
7219 int fe_degree,
7220 int n_q_points_1d,
7221 int n_components_,
7222 typename Number,
7223 typename VectorizedArrayType>
7224template <typename VectorType>
7225inline void
7226FEEvaluation<dim,
7227 fe_degree,
7228 n_q_points_1d,
7229 n_components_,
7230 Number,
7231 VectorizedArrayType>::
7232 integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
7233 VectorType &destination)
7234{
7235 VectorizedArrayType *dst_ptr =
7236 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
7237 *this, destination);
7238 if (dst_ptr != nullptr)
7239 integrate(integration_flag, dst_ptr, true);
7240 else
7241 {
7242 integrate(integration_flag, this->begin_dof_values());
7243 this->distribute_local_to_global(destination);
7244 }
7245}
7246
7247
7248
7249template <int dim,
7250 int fe_degree,
7251 int n_q_points_1d,
7252 int n_components_,
7253 typename Number,
7254 typename VectorizedArrayType>
7256FEEvaluation<dim,
7257 fe_degree,
7258 n_q_points_1d,
7259 n_components_,
7260 Number,
7261 VectorizedArrayType>::dof_indices() const
7262{
7264 0U, dofs_per_cell);
7265}
7266
7267
7268
7269/*-------------------------- FEFaceEvaluation ---------------------------*/
7270
7271
7272
7273template <int dim,
7274 int fe_degree,
7275 int n_q_points_1d,
7276 int n_components_,
7277 typename Number,
7278 typename VectorizedArrayType>
7279inline FEFaceEvaluation<dim,
7280 fe_degree,
7281 n_q_points_1d,
7282 n_components_,
7283 Number,
7284 VectorizedArrayType>::
7285 FEFaceEvaluation(
7287 const bool is_interior_face,
7288 const unsigned int dof_no,
7289 const unsigned int quad_no,
7290 const unsigned int first_selected_component,
7291 const unsigned int active_fe_index,
7292 const unsigned int active_quad_index,
7293 const unsigned int face_type)
7294 : BaseClass(matrix_free,
7295 dof_no,
7296 first_selected_component,
7297 quad_no,
7298 fe_degree,
7299 static_n_q_points,
7300 is_interior_face,
7301 active_fe_index,
7302 active_quad_index,
7303 face_type)
7304 , dofs_per_component(this->data->dofs_per_component_on_cell)
7305 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7306 , n_q_points(this->n_quadrature_points)
7307{}
7308
7309
7310
7311template <int dim,
7312 int fe_degree,
7313 int n_q_points_1d,
7314 int n_components_,
7315 typename Number,
7316 typename VectorizedArrayType>
7317inline FEFaceEvaluation<dim,
7318 fe_degree,
7319 n_q_points_1d,
7320 n_components_,
7321 Number,
7322 VectorizedArrayType>::
7323 FEFaceEvaluation(
7325 const std::pair<unsigned int, unsigned int> &range,
7326 const bool is_interior_face,
7327 const unsigned int dof_no,
7328 const unsigned int quad_no,
7329 const unsigned int first_selected_component)
7330 : FEFaceEvaluation(matrix_free,
7331 is_interior_face,
7332 dof_no,
7333 quad_no,
7334 first_selected_component,
7335 matrix_free.get_face_active_fe_index(range,
7336 is_interior_face),
7338 matrix_free.get_face_info(range.first).face_type)
7339{}
7340
7341
7342
7343template <int dim,
7344 int fe_degree,
7345 int n_q_points_1d,
7346 int n_components_,
7347 typename Number,
7348 typename VectorizedArrayType>
7349inline void
7351 fe_degree,
7352 n_q_points_1d,
7353 n_components_,
7354 Number,
7355 VectorizedArrayType>::reinit(const unsigned int face_index)
7356{
7357 Assert(this->mapped_geometry == nullptr,
7358 ExcMessage("FEEvaluation was initialized without a matrix-free object."
7359 " Integer indexing is not possible"));
7360 if (this->mapped_geometry != nullptr)
7361 return;
7362
7363 this->cell = face_index;
7364 this->dof_access_index =
7365 this->is_interior_face() ?
7368 Assert(this->mapping_data != nullptr, ExcNotInitialized());
7369
7370 if (face_index >=
7371 this->matrix_free->get_task_info().face_partition_data.back() &&
7372 face_index <
7373 this->matrix_free->get_task_info().boundary_partition_data.back())
7374 Assert(this->is_interior_face(),
7375 ExcMessage(
7376 "Boundary faces do not have a neighbor. When looping over "
7377 "boundary faces use FEFaceEvaluation with the parameter "
7378 "is_interior_face set to true. "));
7379
7380 this->reinit_face(this->matrix_free->get_face_info(face_index));
7381
7382 unsigned int i = 0;
7383 for (; i < this->matrix_free->n_active_entries_per_face_batch(this->cell);
7384 ++i)
7385 this->face_ids[i] = face_index * n_lanes + i;
7386 for (; i < n_lanes; ++i)
7387 this->face_ids[i] = numbers::invalid_unsigned_int;
7388
7389 this->cell_type = this->matrix_free->get_mapping_info().face_type[face_index];
7390 const unsigned int offsets =
7391 this->mapping_data->data_index_offsets[face_index];
7392 this->J_value = &this->mapping_data->JxW_values[offsets];
7393 this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
7394 this->jacobian =
7395 &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
7396 this->normal_x_jacobian =
7397 &this->mapping_data
7398 ->normals_times_jacobians[!this->is_interior_face()][offsets];
7399 this->jacobian_gradients =
7400 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7401 offsets;
7402 this->jacobian_gradients_non_inverse =
7403 this->mapping_data
7404 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7405 .data() +
7406 offsets;
7407
7408 if (this->mapping_data->quadrature_point_offsets.empty() == false)
7409 {
7410 AssertIndexRange(this->cell,
7411 this->mapping_data->quadrature_point_offsets.size());
7412 this->quadrature_points =
7413 this->mapping_data->quadrature_points.data() +
7414 this->mapping_data->quadrature_point_offsets[this->cell];
7415 }
7416
7417 if constexpr (running_in_debug_mode())
7418 {
7419 this->is_reinitialized = true;
7420 this->dof_values_initialized = false;
7421 this->values_quad_initialized = false;
7422 this->gradients_quad_initialized = false;
7423 this->hessians_quad_initialized = false;
7424 }
7425}
7426
7427
7428
7429template <int dim,
7430 int fe_degree,
7431 int n_q_points_1d,
7432 int n_components_,
7433 typename Number,
7434 typename VectorizedArrayType>
7435inline void
7437 fe_degree,
7438 n_q_points_1d,
7439 n_components_,
7440 Number,
7441 VectorizedArrayType>::reinit(const unsigned int cell_index,
7442 const unsigned int face_number)
7443{
7444 Assert(
7445 this->quad_no <
7446 this->matrix_free->get_mapping_info().face_data_by_cells.size(),
7447 ExcMessage(
7448 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
7451 this->matrix_free->get_mapping_info().cell_type.size());
7452 Assert(this->mapped_geometry == nullptr,
7453 ExcMessage("FEEvaluation was initialized without a matrix-free object."
7454 " Integer indexing is not possible"));
7455 if (this->mapped_geometry != nullptr)
7456 return;
7457 Assert(this->matrix_free != nullptr, ExcNotInitialized());
7458
7459 this->cell_type = this->matrix_free->get_mapping_info()
7460 .faces_by_cells_type[cell_index][face_number];
7461 this->cell = cell_index;
7462 this->subface_index = GeometryInfo<dim>::max_children_per_cell;
7463 this->dof_access_index =
7465
7466 if (this->is_interior_face() == false)
7467 {
7468 // for this case, we need to look into the FaceInfo field that collects
7469 // information from both sides of a face once for the global mesh, and
7470 // pick the face id that is not the local one (cell_this).
7471 for (unsigned int i = 0; i < n_lanes; ++i)
7472 {
7473 // compute actual (non vectorized) cell ID
7474 const unsigned int cell_this = cell_index * n_lanes + i;
7475 // compute face ID
7476 unsigned int face_index =
7477 this->matrix_free->get_cell_and_face_to_plain_faces()(cell_index,
7478 face_number,
7479 i);
7480
7481 this->face_ids[i] = face_index;
7482
7483 if (face_index == numbers::invalid_unsigned_int)
7484 {
7485 this->cell_ids[i] = numbers::invalid_unsigned_int;
7486 this->face_numbers[i] = static_cast<std::uint8_t>(-1);
7487 this->face_orientations[i] = static_cast<std::uint8_t>(-1);
7488 continue; // invalid face ID: no neighbor on boundary
7489 }
7490
7491 const auto &faces =
7492 this->matrix_free->get_face_info(face_index / n_lanes);
7493 // get cell ID on both sides of face
7494 auto cell_m = faces.cells_interior[face_index % n_lanes];
7495 auto cell_p = faces.cells_exterior[face_index % n_lanes];
7496
7497 const bool face_identifies_as_interior = cell_m != cell_this;
7498
7499 Assert(cell_m == cell_this || cell_p == cell_this,
7501
7502 // compare the IDs with the given cell ID
7503 if (face_identifies_as_interior)
7504 {
7505 this->cell_ids[i] = cell_m; // neighbor has the other ID
7506 this->face_numbers[i] = faces.interior_face_no;
7507 }
7508 else
7509 {
7510 this->cell_ids[i] = cell_p;
7511 this->face_numbers[i] = faces.exterior_face_no;
7512 }
7513
7514 const bool orientation_interior_face = faces.face_orientation >= 8;
7515 unsigned int face_orientation = faces.face_orientation % 8;
7516 if (face_identifies_as_interior != orientation_interior_face)
7517 {
7518 constexpr std::array<std::uint8_t, 8> table{
7519 {0, 1, 6, 3, 4, 5, 2, 7}};
7520 face_orientation = table[face_orientation];
7521 }
7522 this->face_orientations[i] = face_orientation;
7523 }
7524 }
7525 else
7526 {
7527 this->face_orientations[0] = 0;
7528 this->face_numbers[0] = face_number;
7529 if (this->matrix_free->n_active_entries_per_cell_batch(this->cell) ==
7530 n_lanes)
7531 {
7533 for (unsigned int i = 0; i < n_lanes; ++i)
7534 this->cell_ids[i] = cell_index * n_lanes + i;
7535 }
7536 else
7537 {
7538 unsigned int i = 0;
7539 for (; i <
7540 this->matrix_free->n_active_entries_per_cell_batch(this->cell);
7541 ++i)
7542 this->cell_ids[i] = cell_index * n_lanes + i;
7543 for (; i < n_lanes; ++i)
7544 this->cell_ids[i] = numbers::invalid_unsigned_int;
7545 }
7546 for (unsigned int i = 0; i < n_lanes; ++i)
7547 this->face_ids[i] =
7548 this->matrix_free->get_cell_and_face_to_plain_faces()(cell_index,
7549 face_number,
7550 i);
7551 }
7552
7553 const unsigned int offsets =
7554 this->matrix_free->get_mapping_info()
7555 .face_data_by_cells[this->quad_no]
7556 .data_index_offsets[cell_index * GeometryInfo<dim>::faces_per_cell +
7557 face_number];
7558 AssertIndexRange(offsets,
7559 this->matrix_free->get_mapping_info()
7560 .face_data_by_cells[this->quad_no]
7561 .JxW_values.size());
7562 this->J_value = &this->matrix_free->get_mapping_info()
7563 .face_data_by_cells[this->quad_no]
7564 .JxW_values[offsets];
7565 this->normal_vectors = &this->matrix_free->get_mapping_info()
7566 .face_data_by_cells[this->quad_no]
7567 .normal_vectors[offsets];
7568 this->jacobian = &this->matrix_free->get_mapping_info()
7569 .face_data_by_cells[this->quad_no]
7570 .jacobians[!this->is_interior_face()][offsets];
7571 this->normal_x_jacobian =
7572 &this->matrix_free->get_mapping_info()
7573 .face_data_by_cells[this->quad_no]
7574 .normals_times_jacobians[!this->is_interior_face()][offsets];
7575 this->jacobian_gradients =
7576 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7577 offsets;
7578 this->jacobian_gradients_non_inverse =
7579 this->mapping_data
7580 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7581 .data() +
7582 offsets;
7583
7584 if (this->matrix_free->get_mapping_info()
7585 .face_data_by_cells[this->quad_no]
7586 .quadrature_point_offsets.empty() == false)
7587 {
7588 const unsigned int index =
7589 this->cell * GeometryInfo<dim>::faces_per_cell + this->face_numbers[0];
7591 this->matrix_free->get_mapping_info()
7592 .face_data_by_cells[this->quad_no]
7593 .quadrature_point_offsets.size());
7594 this->quadrature_points = this->matrix_free->get_mapping_info()
7595 .face_data_by_cells[this->quad_no]
7596 .quadrature_points.data() +
7597 this->matrix_free->get_mapping_info()
7598 .face_data_by_cells[this->quad_no]
7599 .quadrature_point_offsets[index];
7600 }
7601
7602 if constexpr (running_in_debug_mode())
7603 {
7604 this->is_reinitialized = true;
7605 this->dof_values_initialized = false;
7606 this->values_quad_initialized = false;
7607 this->gradients_quad_initialized = false;
7608 this->hessians_quad_initialized = false;
7609 }
7610}
7611
7612
7613
7614template <int dim,
7615 int fe_degree,
7616 int n_q_points_1d,
7617 int n_components_,
7618 typename Number,
7619 typename VectorizedArrayType>
7620inline void
7622 fe_degree,
7623 n_q_points_1d,
7624 n_components_,
7625 Number,
7626 VectorizedArrayType>::
7627 evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
7628{
7629 if constexpr (running_in_debug_mode())
7630 {
7631 Assert(this->dof_values_initialized, ExcNotInitialized());
7632 }
7633
7634 evaluate(this->values_dofs, evaluation_flag);
7635}
7636
7637
7638
7639template <int dim,
7640 int fe_degree,
7641 int n_q_points_1d,
7642 int n_components_,
7643 typename Number,
7644 typename VectorizedArrayType>
7645inline void
7647 fe_degree,
7648 n_q_points_1d,
7649 n_components_,
7650 Number,
7651 VectorizedArrayType>::
7652 evaluate(const VectorizedArrayType *values_array,
7653 const EvaluationFlags::EvaluationFlags evaluation_flag)
7654{
7655 Assert((evaluation_flag &
7658 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
7659 "and EvaluationFlags::hessians are supported."));
7660
7661 const bool hessians_on_general_cells =
7662 evaluation_flag & EvaluationFlags::hessians &&
7663 (this->cell_type > internal::MatrixFreeFunctions::affine);
7664 EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
7665 if (hessians_on_general_cells)
7666 evaluation_flag_actual |= EvaluationFlags::gradients;
7667
7668 if (this->data->element_type ==
7670 evaluation_flag & EvaluationFlags::gradients &&
7671 (this->cell_type > internal::MatrixFreeFunctions::affine))
7672 evaluation_flag_actual |= EvaluationFlags::values;
7673
7674 if constexpr (fe_degree > -1)
7676 template run<fe_degree, n_q_points_1d>(n_components,
7677 evaluation_flag_actual,
7678 values_array,
7679 *this);
7680 else
7682 n_components, evaluation_flag_actual, values_array, *this);
7683
7684 if constexpr (running_in_debug_mode())
7685 {
7686 this->values_quad_initialized =
7687 evaluation_flag_actual & EvaluationFlags::values;
7688 this->gradients_quad_initialized =
7689 evaluation_flag_actual & EvaluationFlags::gradients;
7690 this->hessians_quad_initialized =
7691 evaluation_flag_actual & EvaluationFlags::hessians;
7692 }
7693}
7694
7695
7696
7697template <int dim,
7698 int fe_degree,
7699 int n_q_points_1d,
7700 int n_components_,
7701 typename Number,
7702 typename VectorizedArrayType>
7703inline void
7705 fe_degree,
7706 n_q_points_1d,
7707 n_components_,
7708 Number,
7709 VectorizedArrayType>::
7710 project_to_face(const EvaluationFlags::EvaluationFlags evaluation_flag)
7711{
7712 if constexpr (running_in_debug_mode())
7713 {
7714 Assert(this->dof_values_initialized, ExcNotInitialized());
7715 }
7716
7717 project_to_face(this->values_dofs, evaluation_flag);
7718}
7719
7720
7721
7722template <int dim,
7723 int fe_degree,
7724 int n_q_points_1d,
7725 int n_components_,
7726 typename Number,
7727 typename VectorizedArrayType>
7728inline void
7730 fe_degree,
7731 n_q_points_1d,
7732 n_components_,
7733 Number,
7734 VectorizedArrayType>::
7735 project_to_face(const VectorizedArrayType *values_array,
7736 const EvaluationFlags::EvaluationFlags evaluation_flag)
7737{
7738 Assert((evaluation_flag &
7741 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
7742 "and EvaluationFlags::hessians are supported."));
7743
7744 const bool hessians_on_general_cells =
7745 evaluation_flag & EvaluationFlags::hessians &&
7746 (this->cell_type > internal::MatrixFreeFunctions::affine);
7747 EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
7748 if (hessians_on_general_cells)
7749 evaluation_flag_actual |= EvaluationFlags::gradients;
7750
7751 if (this->data->element_type ==
7753 evaluation_flag & EvaluationFlags::gradients &&
7754 (this->cell_type > internal::MatrixFreeFunctions::affine))
7755 evaluation_flag_actual |= EvaluationFlags::values;
7756
7757 if constexpr (fe_degree > -1)
7759 dim,
7760 VectorizedArrayType>::template run<fe_degree>(n_components,
7761 evaluation_flag_actual,
7762 values_array,
7763 *this);
7764 else
7766 project_to_face(n_components,
7767 evaluation_flag_actual,
7768 values_array,
7769 *this);
7770
7771 // face dofs initialized
7772}
7773
7774
7775
7776template <int dim,
7777 int fe_degree,
7778 int n_q_points_1d,
7779 int n_components_,
7780 typename Number,
7781 typename VectorizedArrayType>
7782inline void
7784 fe_degree,
7785 n_q_points_1d,
7786 n_components_,
7787 Number,
7788 VectorizedArrayType>::
7789 evaluate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag)
7790{
7791 Assert((evaluation_flag &
7794 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
7795 "and EvaluationFlags::hessians are supported."));
7796
7797 const bool hessians_on_general_cells =
7798 evaluation_flag & EvaluationFlags::hessians &&
7799 (this->cell_type > internal::MatrixFreeFunctions::affine);
7800 EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
7801 if (hessians_on_general_cells)
7802 evaluation_flag_actual |= EvaluationFlags::gradients;
7803
7804 if (this->data->element_type ==
7806 evaluation_flag & EvaluationFlags::gradients &&
7807 (this->cell_type > internal::MatrixFreeFunctions::affine))
7808 evaluation_flag_actual |= EvaluationFlags::values;
7809
7810 if constexpr (fe_degree > -1)
7812 dim,
7813 VectorizedArrayType>::template run<fe_degree>(n_components,
7814 evaluation_flag_actual,
7815 *this);
7816 else
7818 evaluate_in_face(n_components, evaluation_flag_actual, *this);
7819
7820 if constexpr (running_in_debug_mode())
7821 {
7822 this->values_quad_initialized =
7823 evaluation_flag_actual & EvaluationFlags::values;
7824 this->gradients_quad_initialized =
7825 evaluation_flag_actual & EvaluationFlags::gradients;
7826 this->hessians_quad_initialized =
7827 evaluation_flag_actual & EvaluationFlags::hessians;
7828 }
7829}
7830
7831
7832
7833template <int dim,
7834 int fe_degree,
7835 int n_q_points_1d,
7836 int n_components_,
7837 typename Number,
7838 typename VectorizedArrayType>
7839inline void
7841 fe_degree,
7842 n_q_points_1d,
7843 n_components_,
7844 Number,
7845 VectorizedArrayType>::
7846 integrate(const EvaluationFlags::EvaluationFlags integration_flag,
7847 const bool sum_into_values)
7848{
7849 integrate(integration_flag, this->values_dofs, sum_into_values);
7850
7851 if constexpr (running_in_debug_mode())
7852 {
7853 this->dof_values_initialized = true;
7854 }
7855}
7856
7857
7858
7859template <int dim,
7860 int fe_degree,
7861 int n_q_points_1d,
7862 int n_components_,
7863 typename Number,
7864 typename VectorizedArrayType>
7865inline void
7867 fe_degree,
7868 n_q_points_1d,
7869 n_components_,
7870 Number,
7871 VectorizedArrayType>::
7872 integrate(const EvaluationFlags::EvaluationFlags integration_flag,
7873 VectorizedArrayType *values_array,
7874 const bool sum_into_values)
7875{
7876 Assert((integration_flag &
7879 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
7880 "and EvaluationFlags::hessians are supported."));
7881
7882 EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
7883 if (integration_flag & EvaluationFlags::hessians &&
7884 (this->cell_type > internal::MatrixFreeFunctions::affine))
7885 {
7886 unsigned int size = n_components * dim * n_q_points;
7887 if ((integration_flag & EvaluationFlags::gradients) != 0u)
7888 {
7889 for (unsigned int i = 0; i < size; ++i)
7890 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7891 }
7892 else
7893 {
7894 for (unsigned int i = 0; i < size; ++i)
7895 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7896 integration_flag_actual |= EvaluationFlags::gradients;
7897 }
7898 }
7899
7900 if (this->data->element_type ==
7902 integration_flag & EvaluationFlags::gradients &&
7903 this->cell_type > internal::MatrixFreeFunctions::affine &&
7904 this->divergence_is_requested == false)
7905 {
7906 unsigned int size = n_components * n_q_points;
7907 if ((integration_flag & EvaluationFlags::values) != 0u)
7908 {
7909 for (unsigned int i = 0; i < size; ++i)
7910 this->values_quad[i] += this->values_from_gradients_quad[i];
7911 }
7912 else
7913 {
7914 for (unsigned int i = 0; i < size; ++i)
7915 this->values_quad[i] = this->values_from_gradients_quad[i];
7916 integration_flag_actual |= EvaluationFlags::values;
7917 }
7918 }
7919
7920 if constexpr (fe_degree > -1)
7922 template run<fe_degree, n_q_points_1d>(n_components,
7923 integration_flag_actual,
7924 values_array,
7925 *this,
7926 sum_into_values);
7927 else
7929 n_components,
7930 integration_flag_actual,
7931 values_array,
7932 *this,
7933 sum_into_values);
7934}
7935
7936
7937
7938template <int dim,
7939 int fe_degree,
7940 int n_q_points_1d,
7941 int n_components_,
7942 typename Number,
7943 typename VectorizedArrayType>
7944inline void
7946 fe_degree,
7947 n_q_points_1d,
7948 n_components_,
7949 Number,
7950 VectorizedArrayType>::
7951 integrate_in_face(const EvaluationFlags::EvaluationFlags integration_flag)
7952{
7953 Assert((integration_flag &
7956 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
7957 "and EvaluationFlags::hessians are supported."));
7958
7959 EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
7960 if (integration_flag & EvaluationFlags::hessians &&
7961 (this->cell_type > internal::MatrixFreeFunctions::affine))
7962 {
7963 unsigned int size = n_components * dim * n_q_points;
7964 if ((integration_flag & EvaluationFlags::gradients) != 0u)
7965 {
7966 for (unsigned int i = 0; i < size; ++i)
7967 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7968 }
7969 else
7970 {
7971 for (unsigned int i = 0; i < size; ++i)
7972 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7973 integration_flag_actual |= EvaluationFlags::gradients;
7974 }
7975 }
7976
7977 if (this->data->element_type ==
7979 integration_flag & EvaluationFlags::gradients &&
7980 this->cell_type > internal::MatrixFreeFunctions::affine &&
7981 this->divergence_is_requested == false)
7982 {
7983 unsigned int size = n_components * n_q_points;
7984 if ((integration_flag & EvaluationFlags::values) != 0u)
7985 {
7986 for (unsigned int i = 0; i < size; ++i)
7987 this->values_quad[i] += this->values_from_gradients_quad[i];
7988 }
7989 else
7990 {
7991 for (unsigned int i = 0; i < size; ++i)
7992 this->values_quad[i] = this->values_from_gradients_quad[i];
7993 integration_flag_actual |= EvaluationFlags::values;
7994 }
7995 }
7996
7997 if constexpr (fe_degree > -1)
7999 dim,
8000 VectorizedArrayType>::template run<fe_degree>(n_components,
8001 integration_flag_actual,
8002 *this);
8003 else
8005 integrate_in_face(n_components, integration_flag_actual, *this);
8006
8007 // face dofs initialized
8008}
8009
8010
8011
8012template <int dim,
8013 int fe_degree,
8014 int n_q_points_1d,
8015 int n_components_,
8016 typename Number,
8017 typename VectorizedArrayType>
8018inline void
8020 fe_degree,
8021 n_q_points_1d,
8022 n_components_,
8023 Number,
8024 VectorizedArrayType>::
8025 collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag,
8026 const bool sum_into_values)
8027{
8028 collect_from_face(integration_flag, this->values_dofs, sum_into_values);
8029
8030 if constexpr (running_in_debug_mode())
8031 {
8032 this->dof_values_initialized = true;
8033 }
8034}
8035
8036
8037
8038template <int dim,
8039 int fe_degree,
8040 int n_q_points_1d,
8041 int n_components_,
8042 typename Number,
8043 typename VectorizedArrayType>
8044inline void
8046 fe_degree,
8047 n_q_points_1d,
8048 n_components_,
8049 Number,
8050 VectorizedArrayType>::
8051 collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag,
8052 VectorizedArrayType *values_array,
8053 const bool sum_into_values)
8054{
8055 Assert((integration_flag &
8058 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
8059 "and EvaluationFlags::hessians are supported."));
8060
8061 EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
8062 if (integration_flag & EvaluationFlags::hessians &&
8063 (this->cell_type > internal::MatrixFreeFunctions::affine))
8064 integration_flag_actual |= EvaluationFlags::gradients;
8065
8066 if (this->data->element_type ==
8068 integration_flag & EvaluationFlags::gradients &&
8069 this->cell_type > internal::MatrixFreeFunctions::affine &&
8070 this->divergence_is_requested == false)
8071 integration_flag_actual |= EvaluationFlags::values;
8072
8073 if constexpr (fe_degree > -1)
8075 dim,
8076 VectorizedArrayType>::template run<fe_degree>(n_components,
8077 integration_flag_actual,
8078 values_array,
8079 *this,
8080 sum_into_values);
8081 else
8083 collect_from_face(n_components,
8084 integration_flag_actual,
8085 values_array,
8086 *this,
8087 sum_into_values);
8088}
8089
8090
8091
8092template <int dim,
8093 int fe_degree,
8094 int n_q_points_1d,
8095 int n_components_,
8096 typename Number,
8097 typename VectorizedArrayType>
8098template <typename VectorType>
8099inline void
8101 fe_degree,
8102 n_q_points_1d,
8103 n_components_,
8104 Number,
8105 VectorizedArrayType>::
8106 gather_evaluate(const VectorType &input_vector,
8107 const EvaluationFlags::EvaluationFlags evaluation_flag)
8108{
8109 Assert((evaluation_flag &
8112 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
8113 "and EvaluationFlags::hessians are supported."));
8114
8115 const auto shared_vector_data = internal::get_shared_vector_data(
8116 &input_vector,
8117 this->dof_access_index ==
8119 this->active_fe_index,
8120 this->dof_info);
8121
8122 if (this->data->data.front().fe_degree > 0 &&
8123 fast_evaluation_supported(this->data->data.front().fe_degree,
8124 this->data->data.front().n_q_points_1d) &&
8126 dim,
8127 typename VectorType::value_type,
8128 VectorizedArrayType>::
8129 supports(evaluation_flag,
8130 *this->data,
8131 internal::get_beginning<typename VectorType::value_type>(
8132 input_vector),
8133 this->dof_info->index_storage_variants[this->dof_access_index]
8134 [this->cell]))
8135 {
8136 if constexpr (fe_degree > -1)
8137 {
8139 dim,
8140 typename VectorType::value_type,
8141 VectorizedArrayType>::template run<fe_degree,
8142 n_q_points_1d>(
8143 n_components,
8144 evaluation_flag,
8145 internal::get_beginning<typename VectorType::value_type>(
8146 input_vector),
8147 shared_vector_data,
8148 *this);
8149 }
8150 else
8151 {
8153 dim,
8154 typename VectorType::value_type,
8155 VectorizedArrayType>::evaluate(n_components,
8156 evaluation_flag,
8157 internal::get_beginning<
8158 typename VectorType::value_type>(
8159 input_vector),
8160 shared_vector_data,
8161 *this);
8162 }
8163 }
8164 else
8165 {
8166 this->read_dof_values(input_vector);
8167 this->evaluate(evaluation_flag);
8168 }
8169
8170 if constexpr (running_in_debug_mode())
8171 {
8172 this->values_quad_initialized = evaluation_flag & EvaluationFlags::values;
8173 this->gradients_quad_initialized =
8174 evaluation_flag & EvaluationFlags::gradients;
8175 this->hessians_quad_initialized =
8176 evaluation_flag & EvaluationFlags::hessians;
8177 }
8178}
8179
8180
8181
8182template <int dim,
8183 int fe_degree,
8184 int n_q_points_1d,
8185 int n_components_,
8186 typename Number,
8187 typename VectorizedArrayType>
8188template <typename VectorType>
8189inline void
8191 dim,
8192 fe_degree,
8193 n_q_points_1d,
8194 n_components_,
8195 Number,
8196 VectorizedArrayType>::integrate_scatter(const bool integrate_values,
8197 const bool integrate_gradients,
8198 VectorType &destination)
8199{
8201 ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8202 ((integrate_gradients) ? EvaluationFlags::gradients :
8204
8205 integrate_scatter(flag, destination);
8206}
8207
8208
8209
8210template <int dim,
8211 int fe_degree,
8212 int n_q_points_1d,
8213 int n_components_,
8214 typename Number,
8215 typename VectorizedArrayType>
8216template <typename VectorType>
8217inline void
8219 fe_degree,
8220 n_q_points_1d,
8221 n_components_,
8222 Number,
8223 VectorizedArrayType>::
8224 integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
8225 VectorType &destination)
8226{
8227 Assert((this->dof_access_index ==
8229 this->is_interior_face() == false) == false,
8231
8232 const auto shared_vector_data = internal::get_shared_vector_data(
8233 &destination,
8234 this->dof_access_index ==
8236 this->active_fe_index,
8237 this->dof_info);
8238
8239 if (this->data->data.front().fe_degree > 0 &&
8240 fast_evaluation_supported(this->data->data.front().fe_degree,
8241 this->data->data.front().n_q_points_1d) &&
8243 dim,
8244 typename VectorType::value_type,
8245 VectorizedArrayType>::
8246 supports(integration_flag,
8247 *this->data,
8248 internal::get_beginning<typename VectorType::value_type>(
8249 destination),
8250 this->dof_info->index_storage_variants[this->dof_access_index]
8251 [this->cell]))
8252 {
8253 if constexpr (fe_degree > -1)
8254 {
8256 dim,
8257 typename VectorType::value_type,
8258 VectorizedArrayType>::template run<fe_degree,
8259 n_q_points_1d>(
8260 n_components,
8261 integration_flag,
8262 internal::get_beginning<typename VectorType::value_type>(
8263 destination),
8264 shared_vector_data,
8265 *this);
8266 }
8267 else
8268 {
8270 dim,
8271 typename VectorType::value_type,
8272 VectorizedArrayType>::integrate(n_components,
8273 integration_flag,
8274 internal::get_beginning<
8275 typename VectorType::value_type>(
8276 destination),
8277 shared_vector_data,
8278 *this);
8279 }
8280 }
8281 else
8282 {
8283 integrate(integration_flag);
8284 this->distribute_local_to_global(destination);
8285 }
8286}
8287
8288
8289
8290template <int dim,
8291 int fe_degree,
8292 int n_q_points_1d,
8293 int n_components_,
8294 typename Number,
8295 typename VectorizedArrayType>
8298 fe_degree,
8299 n_q_points_1d,
8300 n_components_,
8301 Number,
8302 VectorizedArrayType>::dof_indices() const
8303{
8305 0U, dofs_per_cell);
8306}
8307
8308
8309
8310template <int dim,
8311 int fe_degree,
8312 int n_q_points_1d,
8313 int n_components_,
8314 typename Number,
8315 typename VectorizedArrayType>
8316bool
8317FEEvaluation<dim,
8318 fe_degree,
8319 n_q_points_1d,
8320 n_components_,
8321 Number,
8322 VectorizedArrayType>::
8323 fast_evaluation_supported(const unsigned int given_degree,
8324 const unsigned int given_n_q_points_1d)
8325{
8326 return fe_degree == -1 ?
8328 fast_evaluation_supported(given_degree, given_n_q_points_1d) :
8329 true;
8330}
8331
8332
8333
8334template <int dim,
8335 int fe_degree,
8336 int n_q_points_1d,
8337 int n_components_,
8338 typename Number,
8339 typename VectorizedArrayType>
8340bool
8342 fe_degree,
8343 n_q_points_1d,
8344 n_components_,
8345 Number,
8346 VectorizedArrayType>::
8347 fast_evaluation_supported(const unsigned int given_degree,
8348 const unsigned int given_n_q_points_1d)
8349{
8350 return fe_degree == -1 ?
8352 fast_evaluation_supported(given_degree, given_n_q_points_1d) :
8353 true;
8354}
8355
8356
8357
8358template <int dim,
8359 int fe_degree,
8360 int n_q_points_1d,
8361 int n_components_,
8362 typename Number,
8363 typename VectorizedArrayType>
8364bool
8366 fe_degree,
8367 n_q_points_1d,
8368 n_components_,
8369 Number,
8370 VectorizedArrayType>::at_boundary() const
8371{
8372 Assert(this->dof_access_index !=
8375
8376 if (this->is_interior_face() == false)
8377 return false;
8378 else if (this->cell < this->matrix_free->n_inner_face_batches())
8379 return false;
8380 else if (this->cell < (this->matrix_free->n_inner_face_batches() +
8381 this->matrix_free->n_boundary_face_batches()))
8382 return true;
8383 else
8384 return false;
8385}
8386
8387
8388
8389template <int dim,
8390 int fe_degree,
8391 int n_q_points_1d,
8392 int n_components_,
8393 typename Number,
8394 typename VectorizedArrayType>
8397 fe_degree,
8398 n_q_points_1d,
8399 n_components_,
8400 Number,
8401 VectorizedArrayType>::boundary_id() const
8402{
8403 Assert(this->dof_access_index !=
8406
8407 if (at_boundary())
8408 return this->matrix_free->get_boundary_id(this->cell);
8409 else
8411}
8412
8413
8414
8415template <int dim,
8416 int fe_degree,
8417 int n_q_points_1d,
8418 int n_components_,
8419 typename Number,
8420 typename VectorizedArrayType>
8421unsigned int
8423 dim,
8424 fe_degree,
8425 n_q_points_1d,
8426 n_components_,
8427 Number,
8428 VectorizedArrayType>::get_dofs_per_component_projected_to_face()
8429{
8430 return this->data->dofs_per_component_on_face;
8431}
8432
8433
8434
8435template <int dim,
8436 int fe_degree,
8437 int n_q_points_1d,
8438 int n_components_,
8439 typename Number,
8440 typename VectorizedArrayType>
8441unsigned int
8443 fe_degree,
8444 n_q_points_1d,
8445 n_components_,
8446 Number,
8447 VectorizedArrayType>::get_dofs_projected_to_face()
8448{
8449 return this->data->dofs_per_component_on_face * n_components_;
8450}
8451
8452
8453
8454/*------------------------- end FEFaceEvaluation ------------------------- */
8455
8456
8457#endif // ifndef DOXYGEN
8458
8459
8461
8462#endif
value_type get_dof_value(const unsigned int dof) const
void read_write_operation_global(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors) const
value_type get_laplacian(const unsigned int q_point) const
AlignedVector< VectorizedArrayType > * scratch_data_array
gradient_type get_gradient(const unsigned int q_point) const
void submit_gradient(const Tensor< 2, 1, VectorizedArrayType > val_in, const unsigned int q_point)
void submit_normal_hessian(const value_type normal_hessian_in, const unsigned int q_point)
static constexpr unsigned int dimension
void read_write_operation(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< n_lanes > &mask, const bool apply_constraints=true) const
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
void submit_value(const value_type val_in, const unsigned int q_point)
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
void submit_dof_value(const value_type val_in, const unsigned int dof)
std::vector< types::global_dof_index > local_dof_indices
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
std::conditional_t< n_components_==1, Tensor< 2, dim, VectorizedArrayType >, std::conditional_t< n_components_==dim, Tensor< 3, dim, VectorizedArrayType >, Tensor< 1, n_components_, Tensor< 2, dim, VectorizedArrayType > > > > hessian_type
std::conditional_t< n_components_==1, Tensor< 1, dim, VectorizedArrayType >, std::conditional_t< n_components_==dim, Tensor< 2, dim, VectorizedArrayType >, Tensor< 1, n_components_, Tensor< 1, dim, VectorizedArrayType > > > > gradient_type
void read_dof_values_plain(const VectorType &src, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip())
VectorizedArrayType get_divergence(const unsigned int q_point) const
hessian_type get_hessian(const unsigned int q_point) const
FEEvaluationBase(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
FEEvaluationBase & operator=(const FEEvaluationBase &other)
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
static constexpr unsigned int n_components
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_value(const Tensor< 1, 1, VectorizedArrayType > val_in, const unsigned int q_point)
std::conditional_t< n_components_==1, VectorizedArrayType, Tensor< 1, n_components_, VectorizedArrayType > > value_type
const MatrixFree< dim, Number, VectorizedArrayType > & get_matrix_free() const
void apply_hanging_node_constraints(const bool transpose) const
void set_dof_values_plain(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
FEEvaluationBase(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, const unsigned int active_quad_index, const unsigned int face_type)
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
void read_write_operation_contiguous(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< n_lanes > &mask) const
value_type integrate_value() const
FEEvaluationBase(const FEEvaluationBase &other)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
static constexpr unsigned int n_lanes
const MatrixFree< dim, Number, VectorizedArrayType > * matrix_free
value_type get_normal_derivative(const unsigned int q_point) const
value_type get_value(const unsigned int q_point) const
value_type get_normal_hessian(const unsigned int q_point) const
void read_dof_values(const VectorType &src, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip())
const MappingInfoStorageType::QuadratureDescriptor * descriptor
const MappingInfoStorageType * mapping_data
const ShapeInfoType * data
std::shared_ptr< internal::MatrixFreeFunctions::MappingDataOnTheFly< dim, Number > > mapped_geometry
const DoFInfo * dof_info
FEEvaluationData & operator=(const FEEvaluationData &other)
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
const unsigned int dofs_per_component
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
FEEvaluation(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
const unsigned int n_q_points
FEEvaluation(const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
void reinit(const typename Triangulation< dim >::cell_iterator &cell)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int given_n_q_points_1d)
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
void reinit(const unsigned int cell_batch_index)
void reinit(const std::array< unsigned int, n_lanes > &cell_ids)
FEEvaluation(const FiniteElement< dim > &fe, const FEEvaluationData< dim, VectorizedArrayType, false > &other, const unsigned int first_selected_component=0)
void reinit(const TriaIterator< DoFCellAccessor< dim, dim, level_dof_access > > &cell)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
static constexpr unsigned int tensor_dofs_per_cell
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
typename BaseClass::gradient_type gradient_type
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
static constexpr unsigned int dimension
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int static_dofs_per_cell
typename BaseClass::value_type value_type
FEEvaluation & operator=(const FEEvaluation &other)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag)
static constexpr unsigned int n_components
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int)
static constexpr unsigned int static_n_q_points
static constexpr unsigned int n_lanes
FEEvaluation(const FEEvaluation &other)
void check_template_arguments(const unsigned int fe_no, const unsigned int first_selected_component)
static constexpr unsigned int static_dofs_per_component
typename BaseClass::value_type value_type
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
void collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
bool at_boundary() const
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
static constexpr unsigned int static_n_q_points_cell
void project_to_face(const EvaluationFlags::EvaluationFlags evaluation_flag)
static constexpr unsigned int tensor_dofs_per_cell
unsigned int get_dofs_per_component_projected_to_face()
void project_to_face(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_component
void reinit(const unsigned int face_batch_number)
const unsigned int n_q_points
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
void reinit(const unsigned int cell_batch_number, const unsigned int face_number)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int given_n_q_points_1d)
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int n_components
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
void collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag, const bool sum_into_values=false)
static constexpr unsigned int static_dofs_per_component
static constexpr unsigned int n_lanes
static constexpr unsigned int static_n_q_points
unsigned int get_dofs_projected_to_face()
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
static constexpr unsigned int dimension
typename BaseClass::gradient_type gradient_type
types::boundary_id boundary_id() const
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, const bool sum_into_values=false)
void evaluate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag)
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
static constexpr unsigned int static_dofs_per_cell
void integrate_in_face(const EvaluationFlags::EvaluationFlags integration_flag)
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
unsigned int element_multiplicity(const unsigned int index) const
Abstract base class for mapping classes.
Definition mapping.h:320
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
const Number * constraint_pool_begin(const unsigned int pool_index) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
bool indices_initialized() const
const Number * constraint_pool_end(const unsigned int pool_index) const
unsigned int n_components() const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
const internal::MatrixFreeFunctions::ShapeInfo< Number > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
unsigned int n_base_elements(const unsigned int dof_handler_index) const
Definition point.h:113
DEAL_II_HOST constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
#define DEAL_II_ALWAYS_INLINE
Definition config.h:164
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:211
#define DEAL_II_DEPRECATED
Definition config.h:283
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define DEAL_II_NOT_IMPLEMENTED()
Point< 2 > second
Definition grid_out.cc:4633
Point< 2 > first
Definition grid_out.cc:4632
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcAccessToUninitializedField()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMatrixFreeAccessToUninitializedMappingField(std::string arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
UpdateFlags
@ update_jacobian_grads
Gradient of volume element.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
Definition mpi.cc:746
std::size_t size
Definition mpi.cc:745
The namespace for the EvaluationFlags enum.
EvaluationFlags
The EvaluationFlags enum.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:193
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:466
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
void check_vector_compatibility(const VectorType &vec, const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const internal::MatrixFreeFunctions::DoFInfo &dof_info)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
constexpr types::boundary_id internal_face_boundary_id
Definition types.h:323
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:45
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
static void apply(const unsigned int n_components, const unsigned int fe_degree, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const bool transpose, const std::array< MatrixFreeFunctions::compressed_constraint_kind, VectorizedArrayType::size()> &c_mask, VectorizedArrayType *values)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData< dim, Number, true > &fe_eval)
static void integrate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, FEEvaluationData< dim, Number, true > &fe_eval)
static void collect_from_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static void project_to_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
std::vector< std::pair< unsigned short, unsigned short > > constraint_indicator
Definition dof_info.h:538
std::vector< std::pair< unsigned int, unsigned int > > row_starts
Definition dof_info.h:497
std::vector< std::vector< unsigned int > > component_dof_indices_offset
Definition dof_info.h:690
std::vector< std::vector< bool > > hanging_node_constraint_masks_comp
Definition dof_info.h:520
unsigned int fe_index_from_degree(const unsigned int first_selected_component, const unsigned int fe_degree) const
std::vector< unsigned int > dof_indices
Definition dof_info.h:514
std::vector< compressed_constraint_kind > hanging_node_constraint_masks
Definition dof_info.h:526
std::array< std::vector< unsigned int >, 3 > dof_indices_interleave_strides
Definition dof_info.h:574
std::array< std::vector< std::pair< unsigned int, unsigned int > >, 3 > dof_indices_contiguous_sm
Definition dof_info.h:564
std::vector< unsigned int > row_starts_plain_indices
Definition dof_info.h:637
std::vector< unsigned int > component_to_base_index
Definition dof_info.h:677
std::array< std::vector< unsigned int >, 3 > dof_indices_contiguous
Definition dof_info.h:553
std::vector< unsigned int > plain_dof_indices
Definition dof_info.h:647
std::array< std::vector< unsigned char >, 3 > n_vectorization_lanes_filled
Definition dof_info.h:585
std::vector< unsigned int > dof_indices_interleaved
Definition dof_info.h:543
std::array< std::vector< IndexStorageVariants >, 3 > index_storage_variants
Definition dof_info.h:489
unsigned int quad_index_from_n_q_points(const unsigned int n_q_points) const
std::vector< unsigned int > face_partition_data
Definition task_info.h:496
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)