deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+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\}}\)
Loading...
Searching...
No Matches
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
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
741 const std::array<VectorType *, n_components_> &vectors,
742 const std::array<
745 const std::bitset<n_lanes> &mask,
746 const bool apply_constraints = true) const;
747
755 template <typename VectorType, typename VectorOperation>
756 void
759 const std::array<VectorType *, n_components_> &vectors,
760 const std::array<
763 const std::bitset<n_lanes> &mask) const;
764
772 template <typename VectorType, typename VectorOperation>
773 void
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,
1693
1707 template <typename VectorType>
1708 void
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
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
1878 using gradient_type = typename BaseClass::gradient_type;
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,
2076
2081 void
2083
2088 void
2089 project_to_face(const VectorizedArrayType *values_array,
2091
2096 void
2098
2110 template <typename VectorType>
2111 void
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>
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 {
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) :
2331 0);
2332
2333 init_data.active_quad_index =
2334 fe_degree == numbers::invalid_unsigned_int ?
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,
2372 Number,
2373 is_face,
2374 VectorizedArrayType>::
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,
2441 Number,
2442 is_face,
2443 VectorizedArrayType>::
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),
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 >=
2471 ExcMessage("The underlying element must at least contain as many "
2472 "components as requested by this class"));
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,
2493 Number,
2494 is_face,
2495 VectorizedArrayType>::
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,
2557 Number,
2558 is_face,
2559 VectorizedArrayType> &
2561operator=(const FEEvaluationBase<dim,
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,
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>
2683
2684 template <typename VectorType>
2686 {
2687 using BaseVectorType = const typename VectorType::BlockType;
2688 };
2689
2690 template <typename VectorType>
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>
2704 {
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>
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
2825 const std::array<VectorType *, n_components_> &src,
2826 const std::array<
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)
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;
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]] !=
2940 }
2941
2942 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
2944
2945 const bool use_vectorized_path =
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 &&
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],
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],
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
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];
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])
3060
3062 my_index_start[0].first ||
3063 my_index_start[0].first < dof_info.dof_indices.size(),
3064 ExcIndexRange(0,
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
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] +
3159 }
3160
3161 if (n_components == 1 || this->n_fe_components == 1)
3162 {
3163 unsigned int ind_local = 0;
3165 {
3166 const std::pair<unsigned short, unsigned short> indicator =
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
3224 {
3225 const std::pair<unsigned short, unsigned short> indicator =
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)
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
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
3322 const std::array<VectorType *, n_components_> &src,
3323 const std::array<
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>
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 &&
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],
3379 else
3380 operation.process_dofs_vectorized(dofs_per_component * n_components,
3381 dof_index,
3382 *src[0],
3383 values_dofs[0],
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 &&
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(),
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
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
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],
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],
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
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],
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],
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],
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],
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],
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],
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]);
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
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,
3794 n_components>,
3795 std::array<
3796 const std::vector<ArrayView<const typename internal::BlockVectorSelector<
3797 VectorType,
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,
3812 n_components>,
3813 std::array<
3814 const std::vector<
3815 ArrayView<const typename internal::BlockVectorSelector<
3816 VectorType,
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,
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]),
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
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,
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# ifdef DEBUG
3927 this->dof_values_initialized = true;
3928# endif
3929}
3930
3931
3932
3933template <int dim,
3934 int n_components_,
3935 typename Number,
3936 bool is_face,
3937 typename VectorizedArrayType>
3938template <typename VectorType>
3939inline void
3941 read_dof_values_plain(const VectorType &src,
3942 const unsigned int first_index,
3943 const std::bitset<n_lanes> &mask)
3944{
3945 const auto src_data = internal::get_vector_data<n_components_>(
3946 src,
3948 this->dof_access_index ==
3950 this->active_fe_index,
3951 this->dof_info);
3952
3954 read_write_operation(reader, src_data.first, src_data.second, mask, false);
3955
3956# ifdef DEBUG
3957 this->dof_values_initialized = true;
3958# endif
3959}
3960
3961
3962
3963template <int dim,
3964 int n_components_,
3965 typename Number,
3966 bool is_face,
3967 typename VectorizedArrayType>
3968template <typename VectorType>
3969inline void
3971 distribute_local_to_global(VectorType &dst,
3972 const unsigned int first_index,
3973 const std::bitset<n_lanes> &mask) const
3974{
3975# ifdef DEBUG
3976 Assert(this->dof_values_initialized == true,
3978# endif
3979
3980 apply_hanging_node_constraints(true);
3981
3982 const auto dst_data = internal::get_vector_data<n_components_>(
3983 dst,
3985 this->dof_access_index ==
3987 this->active_fe_index,
3988 this->dof_info);
3989
3992 read_write_operation(distributor, dst_data.first, dst_data.second, mask);
3993}
3994
3995
3996
3997template <int dim,
3998 int n_components_,
3999 typename Number,
4000 bool is_face,
4001 typename VectorizedArrayType>
4002template <typename VectorType>
4003inline void
4005 set_dof_values(VectorType &dst,
4006 const unsigned int first_index,
4007 const std::bitset<n_lanes> &mask) const
4008{
4009# ifdef DEBUG
4010 Assert(this->dof_values_initialized == true,
4012# endif
4013
4014 const auto dst_data = internal::get_vector_data<n_components_>(
4015 dst,
4017 this->dof_access_index ==
4019 this->active_fe_index,
4020 this->dof_info);
4021
4023 read_write_operation(setter, dst_data.first, dst_data.second, mask);
4024}
4025
4026
4027
4028template <int dim,
4029 int n_components_,
4030 typename Number,
4031 bool is_face,
4032 typename VectorizedArrayType>
4033template <typename VectorType>
4034inline void
4036 set_dof_values_plain(VectorType &dst,
4037 const unsigned int first_index,
4038 const std::bitset<n_lanes> &mask) const
4039{
4040# ifdef DEBUG
4041 Assert(this->dof_values_initialized == true,
4043# endif
4044
4045 const auto dst_data = internal::get_vector_data<n_components_>(
4046 dst,
4048 this->dof_access_index ==
4050 this->active_fe_index,
4051 this->dof_info);
4052
4054 read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
4055}
4056
4057
4058
4059/*------------------------------ access to data fields ----------------------*/
4060
4061
4062
4063template <int dim,
4064 int n_components_,
4065 typename Number,
4066 bool is_face,
4067 typename VectorizedArrayType>
4069 typename FEEvaluationBase<dim,
4071 Number,
4072 is_face,
4073 VectorizedArrayType>::value_type
4075 get_dof_value(const unsigned int dof) const
4076{
4077 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4078 if constexpr (n_components == 1)
4079 return this->values_dofs[dof];
4080 else
4081 {
4082 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4083 Tensor<1, n_components_, VectorizedArrayType> return_value;
4084 for (unsigned int comp = 0; comp < n_components; ++comp)
4085 return_value[comp] = this->values_dofs[comp * dofs + dof];
4086 return return_value;
4087 }
4088}
4089
4090
4091
4092template <int dim,
4093 int n_components_,
4094 typename Number,
4095 bool is_face,
4096 typename VectorizedArrayType>
4098 typename FEEvaluationBase<dim,
4100 Number,
4101 is_face,
4102 VectorizedArrayType>::value_type
4104 get_value(const unsigned int q_point) const
4105{
4106# ifdef DEBUG
4107 Assert(this->values_quad_initialized == true,
4109# endif
4110
4111 AssertIndexRange(q_point, this->n_quadrature_points);
4112 if constexpr (n_components == 1)
4113 return this->values_quad[q_point];
4114 else
4115 {
4116 if (n_components == dim &&
4117 this->data->element_type ==
4119 {
4120 // Piola transform is required
4121# ifdef DEBUG
4122 Assert(this->values_quad_initialized == true,
4124# endif
4125
4126 AssertIndexRange(q_point, this->n_quadrature_points);
4127 Assert(this->J_value != nullptr,
4129 "update_values"));
4130 const std::size_t nqp = this->n_quadrature_points;
4132
4133 if (!is_face &&
4135 {
4136 // Cartesian cell
4137 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
4138 const VectorizedArrayType inv_det =
4139 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4140 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4141 this->jacobian[0][2][2];
4142
4143 // J * u * det(J^-1)
4144 for (unsigned int comp = 0; comp < n_components; ++comp)
4145 value_out[comp] = this->values_quad[comp * nqp + q_point] *
4146 jac[comp][comp] * inv_det;
4147 }
4148 else
4149 {
4150 // Affine or general cell
4152 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
4153 this->jacobian[q_point] :
4154 this->jacobian[0];
4156 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
4158 this->jacobian[1];
4159
4160 // Derivatives are reordered for faces. Need to take this into
4161 // account
4162 const VectorizedArrayType inv_det =
4163 (is_face && dim == 2 && this->get_face_no() < 2) ?
4166 // J * u * det(J^-1)
4167 for (unsigned int comp = 0; comp < n_components; ++comp)
4168 {
4169 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
4170 for (unsigned int e = 1; e < dim; ++e)
4171 value_out[comp] +=
4172 this->values_quad[e * nqp + q_point] * jac[comp][e];
4174 }
4175 }
4176 return value_out;
4177 }
4178 else
4179 {
4180 const std::size_t nqp = this->n_quadrature_points;
4182 for (unsigned int comp = 0; comp < n_components; ++comp)
4183 return_value[comp] = this->values_quad[comp * nqp + q_point];
4184 return return_value;
4185 }
4186 }
4187}
4188
4189
4190
4191template <int dim,
4192 int n_components_,
4193 typename Number,
4194 bool is_face,
4195 typename VectorizedArrayType>
4197 typename FEEvaluationBase<dim,
4199 Number,
4200 is_face,
4201 VectorizedArrayType>::gradient_type
4203 get_gradient(const unsigned int q_point) const
4204{
4205# ifdef DEBUG
4206 Assert(this->gradients_quad_initialized == true,
4208# endif
4209
4210 AssertIndexRange(q_point, this->n_quadrature_points);
4211 Assert(this->jacobian != nullptr,
4213 "update_gradients"));
4214 const std::size_t nqp = this->n_quadrature_points;
4215
4216 if constexpr (n_components == dim && dim > 1)
4217 {
4218 if (this->data->element_type ==
4220 {
4221 // Piola transform is required
4222# ifdef DEBUG
4223 Assert(this->gradients_quad_initialized == true,
4225# endif
4226
4227 AssertIndexRange(q_point, this->n_quadrature_points);
4228 Assert(this->jacobian != nullptr,
4230 "update_gradients"));
4231 const std::size_t nqp = this->n_quadrature_points;
4232 const std::size_t nqp_d = nqp * dim;
4234 const VectorizedArrayType *gradients =
4235 this->gradients_quad + q_point * dim;
4236
4237
4238 if (!is_face &&
4240 {
4241 // Cartesian cell
4243 this->jacobian[0];
4245 this->jacobian[1];
4246 const VectorizedArrayType inv_det =
4247 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4248 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4249 this->jacobian[0][2][2];
4250
4251 // J * grad_quad * J^-1 * det(J^-1)
4252 for (unsigned int d = 0; d < dim; ++d)
4253 for (unsigned int comp = 0; comp < n_components; ++comp)
4254 grad_out[comp][d] = gradients[comp * nqp_d + d] *
4255 inv_t_jac[d][d] *
4256 (jac[comp][comp] * inv_det);
4257 }
4258 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4259 {
4260 // Affine cell
4262 this->jacobian[0];
4264 this->jacobian[1];
4265
4266 // Derivatives are reordered for faces. Need to take this into
4267 // account
4268 const VectorizedArrayType inv_det =
4269 (is_face && dim == 2 && this->get_face_no() < 2) ?
4272
4273 VectorizedArrayType tmp[dim][dim];
4274 // J * grad_quad * J^-1 * det(J^-1)
4275 for (unsigned int d = 0; d < dim; ++d)
4276 for (unsigned int e = 0; e < dim; ++e)
4277 {
4278 tmp[d][e] = inv_t_jac[d][0] * gradients[e * nqp_d + 0];
4279 for (unsigned int f = 1; f < dim; ++f)
4280 tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
4281 }
4282 for (unsigned int comp = 0; comp < n_components; ++comp)
4283 for (unsigned int d = 0; d < dim; ++d)
4284 {
4285 VectorizedArrayType res = jac[comp][0] * tmp[d][0];
4286 for (unsigned int f = 1; f < dim; ++f)
4287 res += jac[comp][f] * tmp[d][f];
4288
4289 grad_out[comp][d] = res * inv_det;
4290 }
4291 }
4292 else
4293 {
4294 // General cell
4295
4296 // This assert could be removed if we make sure that this is
4297 // updated even though update_hessians or update_jacobian_grads is
4298 // not passed, i.e make the necessary changes in
4299 // MatrixFreeFunctions::MappingInfoStorage::compute_update_flags
4300 Assert(this->jacobian_gradients_non_inverse != nullptr,
4302 "update_hessians"));
4303
4304 const auto jac_grad =
4305 this->jacobian_gradients_non_inverse[q_point];
4307 this->jacobian[q_point];
4308
4309 // Derivatives are reordered for faces. Need to take this into
4310 // account
4311 const VectorizedArrayType inv_det =
4312 (is_face && dim == 2 && this->get_face_no() < 2) ?
4317
4318 // (J * grad_quad) * J^-1 * det(J^-1), part in braces
4319 VectorizedArrayType tmp[dim][dim];
4320 for (unsigned int d = 0; d < dim; ++d)
4321 for (unsigned int e = 0; e < dim; ++e)
4322 {
4323 tmp[e][d] = t_jac[0][d] * gradients[0 * nqp_d + e];
4324 for (unsigned int f = 1; f < dim; ++f)
4325 tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
4326 }
4327
4328 // Add (jac_grad * values) * J^{-1} * det(J^{-1}), combine terms
4329 // outside braces with gradient part from above
4330 for (unsigned int d = 0; d < dim; ++d)
4331 {
4332 for (unsigned int e = 0; e < dim; ++e)
4333 tmp[e][d] +=
4334 jac_grad[e][d] * this->values_quad[e * nqp + q_point];
4335 for (unsigned int f = 0, r = dim; f < dim; ++f)
4336 for (unsigned int k = f + 1; k < dim; ++k, ++r)
4337 {
4338 tmp[k][d] +=
4339 jac_grad[r][d] * this->values_quad[f * nqp + q_point];
4340 tmp[f][d] +=
4341 jac_grad[r][d] * this->values_quad[k * nqp + q_point];
4342 }
4343 }
4344
4345 // Apply J^{-1} appearing in both terms outside braces above
4346 for (unsigned int d = 0; d < dim; ++d)
4347 for (unsigned int e = 0; e < dim; ++e)
4348 {
4349 VectorizedArrayType res = tmp[0][d] * inv_t_jac[e][0];
4350 for (unsigned int f = 1; f < dim; ++f)
4351 res += tmp[f][d] * inv_t_jac[e][f];
4352 grad_out[d][e] = res;
4353 }
4354
4355 // Add -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1})),
4356 // which can be expressed as a rank-1 update tmp[d] * tmp4[e],
4357 // where tmp = J * values and tmp4 = (J^{-T} * jac_grad * J^{-1})
4358 VectorizedArrayType tmp3[dim], tmp4[dim];
4359 for (unsigned int d = 0; d < dim; ++d)
4360 {
4361 tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0];
4362 for (unsigned int e = 1; e < dim; ++e)
4363 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
4364 }
4365 for (unsigned int e = 0, k = dim; e < dim; ++e)
4366 for (unsigned int f = e + 1; f < dim; ++k, ++f)
4367 for (unsigned int d = 0; d < dim; ++d)
4368 {
4369 tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d];
4370 tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d];
4371 }
4372 for (unsigned int d = 0; d < dim; ++d)
4373 {
4374 tmp4[d] = tmp3[0] * inv_t_jac[d][0];
4375 for (unsigned int e = 1; e < dim; ++e)
4376 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
4377 }
4378
4379 VectorizedArrayType tmp2[dim];
4380 for (unsigned int d = 0; d < dim; ++d)
4381 {
4382 tmp2[d] = t_jac[0][d] * this->values_quad[q_point];
4383 for (unsigned e = 1; e < dim; ++e)
4384 tmp2[d] +=
4385 t_jac[e][d] * this->values_quad[e * nqp + q_point];
4386 }
4387
4388 for (unsigned int d = 0; d < dim; ++d)
4389 for (unsigned int e = 0; e < dim; ++e)
4390 {
4391 grad_out[d][e] -= tmp4[e] * tmp2[d];
4392
4393 // finally multiply by det(J^{-1}) necessary in all
4394 // contributions above
4395 grad_out[d][e] *= inv_det;
4396 }
4397 }
4398 return grad_out;
4399 }
4400 }
4402
4403 // Cartesian cell
4404 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4405 {
4406 for (unsigned int comp = 0; comp < n_components; ++comp)
4407 for (unsigned int d = 0; d < dim; ++d)
4408 grad_out[comp][d] =
4409 this->gradients_quad[(comp * nqp + q_point) * dim + d] *
4410 this->jacobian[0][d][d];
4411 }
4412 // cell with general/affine Jacobian
4413 else
4414 {
4416 this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
4417 q_point :
4418 0];
4419 for (unsigned int comp = 0; comp < n_components; ++comp)
4420 for (unsigned int d = 0; d < dim; ++d)
4421 {
4422 grad_out[comp][d] =
4423 jac[d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4424 for (unsigned int e = 1; e < dim; ++e)
4425 grad_out[comp][d] +=
4426 jac[d][e] *
4427 this->gradients_quad[(comp * nqp + q_point) * dim + e];
4428 }
4429 }
4430 if constexpr (n_components == 1)
4431 return grad_out[0];
4432 else
4433 return grad_out;
4434}
4435
4436
4437
4438template <int dim,
4439 int n_components_,
4440 typename Number,
4441 bool is_face,
4442 typename VectorizedArrayType>
4444 typename FEEvaluationBase<dim,
4446 Number,
4447 is_face,
4448 VectorizedArrayType>::value_type
4450 get_normal_derivative(const unsigned int q_point) const
4451{
4452 AssertIndexRange(q_point, this->n_quadrature_points);
4453# ifdef DEBUG
4454 Assert(this->gradients_quad_initialized == true,
4456# endif
4457
4458 Assert(this->normal_x_jacobian != nullptr,
4460 "update_gradients"));
4461
4462 const std::size_t nqp = this->n_quadrature_points;
4464
4465 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4466 for (unsigned int comp = 0; comp < n_components; ++comp)
4467 grad_out[comp] =
4468 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4469 (this->normal_x_jacobian[0][dim - 1]);
4470 else
4471 {
4472 const std::size_t index =
4473 this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
4474 for (unsigned int comp = 0; comp < n_components; ++comp)
4475 {
4476 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4477 this->normal_x_jacobian[index][0];
4478 for (unsigned int d = 1; d < dim; ++d)
4479 grad_out[comp] +=
4480 this->gradients_quad[(comp * nqp + q_point) * dim + d] *
4481 this->normal_x_jacobian[index][d];
4482 }
4483 }
4484 if constexpr (n_components == 1)
4485 return grad_out[0];
4486 else
4487 return grad_out;
4488}
4489
4490
4491
4492namespace internal
4493{
4494 // compute tmp = hess_unit(u) * J^T. do this manually because we do not
4495 // store the lower diagonal because of symmetry
4496 template <typename VectorizedArrayType>
4497 inline void
4499 const VectorizedArrayType *const hessians,
4500 const unsigned int,
4501 VectorizedArrayType (&tmp)[1][1])
4502 {
4503 tmp[0][0] = jac[0][0] * hessians[0];
4504 }
4505
4506 template <typename VectorizedArrayType>
4507 inline void
4509 const VectorizedArrayType *const hessians,
4510 const unsigned int nqp,
4511 VectorizedArrayType (&tmp)[2][2])
4512 {
4513 for (unsigned int d = 0; d < 2; ++d)
4514 {
4515 tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
4516 tmp[1][d] =
4517 (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
4518 }
4519 }
4520
4521 template <typename VectorizedArrayType>
4522 inline void
4524 const VectorizedArrayType *const hessians,
4525 const unsigned int nqp,
4526 VectorizedArrayType (&tmp)[3][3])
4527 {
4528 for (unsigned int d = 0; d < 3; ++d)
4529 {
4530 tmp[0][d] =
4531 (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
4532 jac[d][2] * hessians[4 * nqp]);
4533 tmp[1][d] =
4534 (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
4535 jac[d][2] * hessians[5 * nqp]);
4536 tmp[2][d] =
4537 (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
4538 jac[d][2] * hessians[2 * nqp]);
4539 }
4540 }
4541} // namespace internal
4542
4543
4544
4545template <int dim,
4546 int n_components_,
4547 typename Number,
4548 bool is_face,
4549 typename VectorizedArrayType>
4550inline typename FEEvaluationBase<dim,
4552 Number,
4553 is_face,
4554 VectorizedArrayType>::hessian_type
4556 get_hessian(const unsigned int q_point) const
4557{
4558# ifdef DEBUG
4559 Assert(this->hessians_quad_initialized == true,
4561# endif
4562 AssertIndexRange(q_point, this->n_quadrature_points);
4563
4564 Assert(this->jacobian != nullptr,
4566 "update_hessian"));
4568 this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4569 0 :
4570 q_point];
4571
4573
4574 const std::size_t nqp = this->n_quadrature_points;
4575 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4576
4577 // Cartesian cell
4578 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4579 {
4580 for (unsigned int comp = 0; comp < n_components; ++comp)
4581 {
4582 for (unsigned int d = 0; d < dim; ++d)
4583 hessian_out[comp][d][d] =
4584 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4585 (jac[d][d] * jac[d][d]);
4586 switch (dim)
4587 {
4588 case 1:
4589 break;
4590 case 2:
4591 hessian_out[comp][0][1] =
4592 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4593 (jac[0][0] * jac[1][1]);
4594 break;
4595 case 3:
4596 hessian_out[comp][0][1] =
4597 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4598 (jac[0][0] * jac[1][1]);
4599 hessian_out[comp][0][2] =
4600 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4601 (jac[0][0] * jac[2][2]);
4602 hessian_out[comp][1][2] =
4603 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4604 (jac[1][1] * jac[2][2]);
4605 break;
4606 default:
4608 }
4609 for (unsigned int d = 0; d < dim; ++d)
4610 for (unsigned int e = d + 1; e < dim; ++e)
4611 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4612 }
4613 }
4614 // cell with general Jacobian, but constant within the cell
4615 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4616 {
4617 for (unsigned int comp = 0; comp < n_components; ++comp)
4618 {
4619 VectorizedArrayType tmp[dim][dim];
4620 internal::hessian_unit_times_jac(
4621 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4622
4623 // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4624 for (unsigned int d = 0; d < dim; ++d)
4625 for (unsigned int e = d; e < dim; ++e)
4626 {
4627 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4628 for (unsigned int f = 1; f < dim; ++f)
4629 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4630 }
4631
4632 // no J' * grad(u) part here because the Jacobian is constant
4633 // throughout the cell and hence, its derivative is zero
4634
4635 // take symmetric part
4636 for (unsigned int d = 0; d < dim; ++d)
4637 for (unsigned int e = d + 1; e < dim; ++e)
4638 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4639 }
4640 }
4641 // cell with general Jacobian
4642 else
4643 {
4644 const auto &jac_grad = this->jacobian_gradients[q_point];
4645 for (unsigned int comp = 0; comp < n_components; ++comp)
4646 {
4647 VectorizedArrayType tmp[dim][dim];
4648 internal::hessian_unit_times_jac(
4649 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4650
4651 // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4652 for (unsigned int d = 0; d < dim; ++d)
4653 for (unsigned int e = d; e < dim; ++e)
4654 {
4655 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4656 for (unsigned int f = 1; f < dim; ++f)
4657 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4658 }
4659
4660 // add diagonal part of J' * grad(u)
4661 for (unsigned int d = 0; d < dim; ++d)
4662 for (unsigned int e = 0; e < dim; ++e)
4663 hessian_out[comp][d][d] +=
4664 jac_grad[d][e] *
4665 this->gradients_quad[(comp * nqp + q_point) * dim + e];
4666
4667 // add off-diagonal part of J' * grad(u)
4668 for (unsigned int d = 0, count = dim; d < dim; ++d)
4669 for (unsigned int e = d + 1; e < dim; ++e, ++count)
4670 for (unsigned int f = 0; f < dim; ++f)
4671 hessian_out[comp][d][e] +=
4672 jac_grad[count][f] *
4673 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4674
4675 // take symmetric part
4676 for (unsigned int d = 0; d < dim; ++d)
4677 for (unsigned int e = d + 1; e < dim; ++e)
4678 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4679 }
4680 }
4681 if constexpr (n_components == 1)
4682 return hessian_out[0];
4683 else
4684 return hessian_out;
4685}
4686
4687
4688
4689template <int dim,
4690 int n_components_,
4691 typename Number,
4692 bool is_face,
4693 typename VectorizedArrayType>
4694inline typename FEEvaluationBase<dim,
4696 Number,
4697 is_face,
4698 VectorizedArrayType>::gradient_type
4700 get_hessian_diagonal(const unsigned int q_point) const
4701{
4702 Assert(!is_face, ExcNotImplemented());
4703# ifdef DEBUG
4704 Assert(this->hessians_quad_initialized == true,
4706# endif
4707 AssertIndexRange(q_point, this->n_quadrature_points);
4708
4709 Assert(this->jacobian != nullptr, ExcNotImplemented());
4711 this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4712 0 :
4713 q_point];
4714
4715 const std::size_t nqp = this->n_quadrature_points;
4716 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4718
4719 // Cartesian cell
4720 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4721 {
4722 for (unsigned int comp = 0; comp < n_components; ++comp)
4723 for (unsigned int d = 0; d < dim; ++d)
4724 hessian_out[comp][d] =
4725 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4726 (jac[d][d] * jac[d][d]);
4727 }
4728 // cell with general Jacobian, but constant within the cell
4729 else if (this->cell_type == internal::MatrixFreeFunctions::affine)
4730 {
4731 for (unsigned int comp = 0; comp < n_components; ++comp)
4732 {
4733 // compute laplacian before the gradient because it needs to access
4734 // unscaled gradient data
4735 VectorizedArrayType tmp[dim][dim];
4736 internal::hessian_unit_times_jac(
4737 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4738
4739 // compute only the trace part of hessian, J * tmp = J *
4740 // hess_unit(u) * J^T
4741 for (unsigned int d = 0; d < dim; ++d)
4742 {
4743 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4744 for (unsigned int f = 1; f < dim; ++f)
4745 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4746 }
4747 }
4748 }
4749 // cell with general Jacobian
4750 else
4751 {
4752 const auto &jac_grad = this->jacobian_gradients[q_point];
4753 for (unsigned int comp = 0; comp < n_components; ++comp)
4754 {
4755 // compute laplacian before the gradient because it needs to access
4756 // unscaled gradient data
4757 VectorizedArrayType tmp[dim][dim];
4758 internal::hessian_unit_times_jac(
4759 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4760
4761 // compute only the trace part of hessian, J * tmp = J *
4762 // hess_unit(u) * J^T
4763 for (unsigned int d = 0; d < dim; ++d)
4764 {
4765 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4766 for (unsigned int f = 1; f < dim; ++f)
4767 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4768 }
4769
4770 for (unsigned int d = 0; d < dim; ++d)
4771 for (unsigned int e = 0; e < dim; ++e)
4772 hessian_out[comp][d] +=
4773 jac_grad[d][e] *
4774 this->gradients_quad[(comp * nqp + q_point) * dim + e];
4775 }
4776 }
4777
4778 if constexpr (n_components == 1)
4779 return hessian_out[0];
4780 else
4781 return hessian_out;
4782}
4783
4784
4785
4786template <int dim,
4787 int n_components_,
4788 typename Number,
4789 bool is_face,
4790 typename VectorizedArrayType>
4791inline typename FEEvaluationBase<dim,
4793 Number,
4794 is_face,
4795 VectorizedArrayType>::value_type
4797 get_laplacian(const unsigned int q_point) const
4798{
4799 Assert(is_face == false, ExcNotImplemented());
4800# ifdef DEBUG
4801 Assert(this->hessians_quad_initialized == true,
4803# endif
4804 AssertIndexRange(q_point, this->n_quadrature_points);
4805
4806 const gradient_type hess_diag = get_hessian_diagonal(q_point);
4807 if constexpr (n_components == 1)
4808 {
4809 VectorizedArrayType sum = hess_diag[0];
4810 for (unsigned int d = 1; d < dim; ++d)
4811 sum += hess_diag[d];
4812 return sum;
4813 }
4814 else
4815 {
4817 for (unsigned int comp = 0; comp < n_components; ++comp)
4818 {
4820 for (unsigned int d = 1; d < dim; ++d)
4822 }
4823 return laplacian_out;
4824 }
4825}
4826
4827
4828
4829template <int dim,
4830 int n_components_,
4831 typename Number,
4832 bool is_face,
4833 typename VectorizedArrayType>
4834inline typename FEEvaluationBase<dim,
4836 Number,
4837 is_face,
4838 VectorizedArrayType>::value_type
4840 get_normal_hessian(const unsigned int q_point) const
4841{
4842# ifdef DEBUG
4843 Assert(this->hessians_quad_initialized == true,
4845# endif
4846 AssertIndexRange(q_point, this->n_quadrature_points);
4847
4848 Assert(this->normal_x_jacobian != nullptr,
4850 "update_hessians"));
4851
4853
4854 const std::size_t nqp = this->n_quadrature_points;
4855 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4856
4857 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4858 {
4859 const auto nxj = this->normal_x_jacobian[0];
4860
4861 for (unsigned int comp = 0; comp < n_components; ++comp)
4862 {
4863 for (unsigned int d = 0; d < dim; ++d)
4864 hessian_out[comp] +=
4865 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4866 (nxj[d]) * (nxj[d]);
4867
4868 switch (dim)
4869 {
4870 case 1:
4871 break;
4872 case 2:
4873 hessian_out[comp] +=
4874 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4875 (nxj[0] * nxj[1]);
4876 break;
4877 case 3:
4878 hessian_out[comp] +=
4879 2. * this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4880 (nxj[0] * nxj[1]);
4881 hessian_out[comp] +=
4882 2. * this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4883 (nxj[0] * nxj[2]);
4884 hessian_out[comp] +=
4885 2. * this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4886 (nxj[1] * nxj[2]);
4887 break;
4888 default:
4890 }
4891 }
4892 }
4893 // cell with general Jacobian
4894 else
4895 {
4896 const auto normal = this->normal_vector(q_point);
4897 const auto hessian = get_hessian(q_point);
4898
4899 if constexpr (n_components == 1)
4900 hessian_out[0] = hessian * normal * normal;
4901 else
4902 for (unsigned int comp = 0; comp < n_components; ++comp)
4903 hessian_out[comp] = hessian[comp] * normal * normal;
4904 }
4905 if constexpr (n_components == 1)
4906 return hessian_out[0];
4907 else
4908 return hessian_out;
4909}
4910
4911
4912
4913template <int dim,
4914 int n_components_,
4915 typename Number,
4916 bool is_face,
4917 typename VectorizedArrayType>
4918inline DEAL_II_ALWAYS_INLINE void
4920 submit_dof_value(const value_type val_in, const unsigned int dof)
4921{
4922# ifdef DEBUG
4923 this->dof_values_initialized = true;
4924# endif
4925 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4926 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4927 for (unsigned int comp = 0; comp < n_components; ++comp)
4928 if constexpr (n_components == 1)
4929 this->values_dofs[comp * dofs + dof] = val_in;
4930 else
4931 this->values_dofs[comp * dofs + dof] = val_in[comp];
4932}
4933
4934
4935
4936template <int dim,
4937 int n_components_,
4938 typename Number,
4939 bool is_face,
4940 typename VectorizedArrayType>
4941inline DEAL_II_ALWAYS_INLINE void
4943 submit_value(const value_type val_in, const unsigned int q_point)
4944{
4945# ifdef DEBUG
4946 Assert(this->is_reinitialized, ExcNotInitialized());
4947# endif
4948 AssertIndexRange(q_point, this->n_quadrature_points);
4949 Assert(this->J_value != nullptr,
4951 "update_values"));
4952# ifdef DEBUG
4953 this->values_quad_submitted = true;
4954# endif
4955
4956 const std::size_t nqp = this->n_quadrature_points;
4957 VectorizedArrayType *values = this->values_quad + q_point;
4958
4959 const VectorizedArrayType JxW =
4960 this->cell_type <= internal::MatrixFreeFunctions::affine ?
4961 this->J_value[0] * this->quadrature_weights[q_point] :
4962 this->J_value[q_point];
4963 if constexpr (n_components == 1)
4964 values[0] = val_in * JxW;
4965 else
4966 {
4967 if (n_components == dim &&
4968 this->data->element_type ==
4970 {
4971 // Piola transform is required
4972 AssertIndexRange(q_point, this->n_quadrature_points);
4973 Assert(this->J_value != nullptr,
4975 "update_value"));
4976# ifdef DEBUG
4977 Assert(this->is_reinitialized, ExcNotInitialized());
4978 this->values_quad_submitted = true;
4979# endif
4980
4981 VectorizedArrayType *values = this->values_quad + q_point;
4982 const std::size_t nqp = this->n_quadrature_points;
4983
4984 if (!is_face &&
4986 {
4987 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
4988 const VectorizedArrayType weight =
4989 this->quadrature_weights[q_point];
4990
4991 for (unsigned int comp = 0; comp < n_components; ++comp)
4992 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
4993 }
4994 else
4995 {
4996 // Affine or general cell
4998 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
4999 this->jacobian[q_point] :
5000 this->jacobian[0];
5001
5002 // Derivatives are reordered for faces. Need to take this into
5003 // account and 1/inv_det != J_value for faces
5004 const VectorizedArrayType fac =
5005 (!is_face) ?
5006 this->quadrature_weights[q_point] :
5007 (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
5008 this->J_value[q_point] :
5009 this->J_value[0] * this->quadrature_weights[q_point]) *
5010 ((dim == 2 && this->get_face_no() < 2) ?
5014 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5016 this->jacobian[1];
5017
5018 // J^T * u * factor
5019 for (unsigned int comp = 0; comp < n_components; ++comp)
5020 {
5021 values[comp * nqp] = val_in[0] * jac[0][comp];
5022 for (unsigned int e = 1; e < dim; ++e)
5023 values[comp * nqp] += val_in[e] * jac[e][comp];
5024 values[comp * nqp] *= fac;
5025 }
5026 }
5027 }
5028 else
5029 for (unsigned int comp = 0; comp < n_components; ++comp)
5030 values[comp * nqp] = val_in[comp] * JxW;
5031 }
5032}
5033
5034
5035
5036template <int dim,
5037 int n_components_,
5038 typename Number,
5039 bool is_face,
5040 typename VectorizedArrayType>
5041template <int, typename>
5042inline DEAL_II_ALWAYS_INLINE void
5045 const unsigned int q_point)
5046{
5047 static_assert(n_components == 1,
5048 "Do not try to modify the default template parameters used for"
5049 " selectively enabling this function via std::enable_if!");
5050 submit_value(val_in[0], q_point);
5051}
5052
5053
5054
5055template <int dim,
5056 int n_components_,
5057 typename Number,
5058 bool is_face,
5059 typename VectorizedArrayType>
5060inline DEAL_II_ALWAYS_INLINE void
5062 submit_gradient(const gradient_type grad_in, const unsigned int q_point)
5063{
5064# ifdef DEBUG
5065 Assert(this->is_reinitialized, ExcNotInitialized());
5066# endif
5067 AssertIndexRange(q_point, this->n_quadrature_points);
5068 Assert(this->J_value != nullptr,
5070 "update_gradients"));
5071 Assert(this->jacobian != nullptr,
5073 "update_gradients"));
5074# ifdef DEBUG
5075 this->gradients_quad_submitted = true;
5076# endif
5077
5078 if constexpr (dim > 1 && n_components == dim)
5079 {
5080 if (this->data->element_type ==
5082 {
5083 // Piola transform is required
5084
5085# ifdef DEBUG
5086 Assert(this->is_reinitialized, ExcNotInitialized());
5087# endif
5088 AssertIndexRange(q_point, this->n_quadrature_points);
5089 Assert(this->J_value != nullptr,
5091 "update_gradients"));
5092 Assert(this->jacobian != nullptr,
5094 "update_gradients"));
5095# ifdef DEBUG
5096 this->gradients_quad_submitted = true;
5097# endif
5098
5099 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5100 VectorizedArrayType *values =
5101 this->values_from_gradients_quad + q_point;
5102 const std::size_t nqp = this->n_quadrature_points;
5103 const std::size_t nqp_d = nqp * dim;
5104
5105 if (!is_face &&
5107 {
5108 // Cartesian cell
5110 this->jacobian[0];
5112 this->jacobian[1];
5113 const VectorizedArrayType weight =
5114 this->quadrature_weights[q_point];
5115 for (unsigned int d = 0; d < dim; ++d)
5116 for (unsigned int comp = 0; comp < n_components; ++comp)
5117 gradients[comp * nqp_d + d] = grad_in[comp][d] *
5118 inv_t_jac[d][d] *
5119 (jac[comp][comp] * weight);
5120 }
5121 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5122 {
5123 // Affine cell
5125 this->jacobian[0];
5127 this->jacobian[1];
5128
5129 // Derivatives are reordered for faces. Need to take this into
5130 // account and 1/inv_det != J_value for faces
5131 const VectorizedArrayType fac =
5132 (!is_face) ?
5133 this->quadrature_weights[q_point] :
5134 this->J_value[0] * this->quadrature_weights[q_point] *
5135 ((dim == 2 && this->get_face_no() < 2) ?
5138
5139 // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
5140 VectorizedArrayType tmp[dim][dim];
5141 for (unsigned int d = 0; d < dim; ++d)
5142 for (unsigned int e = 0; e < dim; ++e)
5143 {
5144 tmp[d][e] = inv_t_jac[0][d] * grad_in[e][0];
5145 for (unsigned int f = 1; f < dim; ++f)
5146 tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
5147 }
5148 for (unsigned int comp = 0; comp < n_components; ++comp)
5149 for (unsigned int d = 0; d < dim; ++d)
5150 {
5151 VectorizedArrayType res = jac[0][comp] * tmp[d][0];
5152 for (unsigned int f = 1; f < dim; ++f)
5153 res += jac[f][comp] * tmp[d][f];
5154
5155 gradients[comp * nqp_d + d] = res * fac;
5156 }
5157 }
5158 else
5159 {
5160 // General cell
5161
5162 const auto jac_grad =
5163 this->jacobian_gradients_non_inverse[q_point];
5165 this->jacobian[q_point];
5166
5167 // Derivatives are reordered for faces. Need to take this into
5168 // account and 1/inv_det != J_value for faces
5169 const VectorizedArrayType fac =
5170 (!is_face) ? this->quadrature_weights[q_point] :
5171 this->J_value[q_point] *
5172 ((dim == 2 && this->get_face_no() < 2) ?
5177
5178 // Start evaluation for values part below to enable the compiler
5179 // to possibly re-use the same computation in get_gradient()
5180 // without interfering with stores to 'gradients'
5181 VectorizedArrayType tmp3[dim], tmp4[dim];
5182 for (unsigned int d = 0; d < dim; ++d)
5183 {
5184 tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0];
5185 for (unsigned int e = 1; e < dim; ++e)
5186 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
5187 }
5188 for (unsigned int e = 0, k = dim; e < dim; ++e)
5189 for (unsigned int f = e + 1; f < dim; ++k, ++f)
5190 for (unsigned int d = 0; d < dim; ++d)
5191 {
5192 tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d];
5193 tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d];
5194 }
5195 for (unsigned int d = 0; d < dim; ++d)
5196 {
5197 tmp4[d] = tmp3[0] * inv_t_jac[d][0];
5198 for (unsigned int e = 1; e < dim; ++e)
5199 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
5200 }
5201
5203 fac * grad_in;
5204
5205 VectorizedArrayType tmp[dim][dim];
5206
5207 // J * (J^{-1} * (grad_in * factor))
5208 for (unsigned int d = 0; d < dim; ++d)
5209 for (unsigned int e = 0; e < dim; ++e)
5210 {
5211 tmp[d][e] = inv_t_jac[0][d] * grad_in_scaled[e][0];
5212 for (unsigned int f = 1; f < dim; ++f)
5213 tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
5214 }
5215
5216 for (unsigned int d = 0; d < dim; ++d)
5217 for (unsigned int e = 0; e < dim; ++e)
5218 {
5219 VectorizedArrayType res = t_jac[d][0] * tmp[e][0];
5220 for (unsigned int f = 1; f < dim; ++f)
5221 res += t_jac[d][f] * tmp[e][f];
5222
5223 gradients[d * nqp_d + e] = res;
5224 }
5225
5226 // jac_grad * (J^{-1} * (grad_in * factor)), re-use part in braces
5227 // as 'tmp' from above
5228 VectorizedArrayType value[dim];
5229 for (unsigned int d = 0; d < dim; ++d)
5230 {
5231 value[d] = tmp[d][0] * jac_grad[d][0];
5232 for (unsigned int e = 1; e < dim; ++e)
5233 value[d] += tmp[d][e] * jac_grad[d][e];
5234 }
5235 for (unsigned int e = 0, k = dim; e < dim; ++e)
5236 for (unsigned int f = e + 1; f < dim; ++k, ++f)
5237 for (unsigned int d = 0; d < dim; ++d)
5238 {
5239 value[e] += tmp[f][d] * jac_grad[k][d];
5240 value[f] += tmp[e][d] * jac_grad[k][d];
5241 }
5242
5243 // -(grad_in * factor) * J * (J^{-T} * jac_grad * J^{-1})
5244 // = -(grad_in * factor) * J * ( \------- tmp4 ---------/ )
5245 for (unsigned int d = 0; d < dim; ++d)
5246 {
5247 VectorizedArrayType tmp2 = grad_in_scaled[d][0] * tmp4[0];
5248 for (unsigned int e = 1; e < dim; ++e)
5249 tmp2 += grad_in_scaled[d][e] * tmp4[e];
5250 for (unsigned int e = 0; e < dim; ++e)
5251 value[e] -= t_jac[e][d] * tmp2;
5252 }
5253
5254 for (unsigned int d = 0; d < dim; ++d)
5255 values[d * nqp] = value[d];
5256 }
5257 return;
5258 }
5259 }
5260
5261 const std::size_t nqp_d = this->n_quadrature_points * dim;
5262 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5263
5264 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5265 {
5266 const VectorizedArrayType JxW =
5267 this->J_value[0] * this->quadrature_weights[q_point];
5268
5269 // Load all entries before starting to write back to make sure the
5270 // compiler sees opportunity of loads in a possibly nearby
5271 // get_gradient() function (i.e., the compiler should not think that
5272 // 'jacobian' could alias with 'gradients').
5273 std::array<VectorizedArrayType, dim> jac;
5274 for (unsigned int d = 0; d < dim; ++d)
5275 jac[d] = this->jacobian[0][d][d];
5276
5277 for (unsigned int d = 0; d < dim; ++d)
5278 {
5279 const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW;
5280 if constexpr (n_components == 1)
5281 gradients[d] = grad_in[d] * factor;
5282 else
5283 for (unsigned int comp = 0; comp < n_components; ++comp)
5284 gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5285 }
5286 }
5287 else
5288 {
5290 this->cell_type > internal::MatrixFreeFunctions::affine ?
5291 this->jacobian[q_point] :
5292 this->jacobian[0];
5293 const VectorizedArrayType JxW =
5294 this->cell_type > internal::MatrixFreeFunctions::affine ?
5295 this->J_value[q_point] :
5296 this->J_value[0] * this->quadrature_weights[q_point];
5297 if constexpr (n_components == 1)
5298 for (unsigned int d = 0; d < dim; ++d)
5299 {
5300 VectorizedArrayType new_val = jac[0][d] * grad_in[0];
5301 for (unsigned int e = 1; e < dim; ++e)
5302 new_val += (jac[e][d] * grad_in[e]);
5303 gradients[d] = new_val * JxW;
5304 }
5305 else
5306 for (unsigned int comp = 0; comp < n_components; ++comp)
5307 for (unsigned int d = 0; d < dim; ++d)
5308 {
5309 VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
5310 for (unsigned int e = 1; e < dim; ++e)
5311 new_val += (jac[e][d] * grad_in[comp][e]);
5312 gradients[comp * nqp_d + d] = new_val * JxW;
5313 }
5314 }
5315}
5316
5317
5318
5319template <int dim,
5320 int n_components_,
5321 typename Number,
5322 bool is_face,
5323 typename VectorizedArrayType>
5324template <int, typename>
5325inline DEAL_II_ALWAYS_INLINE void
5328 const unsigned int q_point)
5329{
5330 static_assert(n_components == 1 && dim == 1,
5331 "Do not try to modify the default template parameters used for"
5332 " selectively enabling this function via std::enable_if!");
5333 submit_gradient(grad_in[0], q_point);
5334}
5335
5336
5337
5338template <int dim,
5339 int n_components_,
5340 typename Number,
5341 bool is_face,
5342 typename VectorizedArrayType>
5343inline DEAL_II_ALWAYS_INLINE void
5345 submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
5346{
5347 AssertIndexRange(q_point, this->n_quadrature_points);
5348 Assert(this->normal_x_jacobian != nullptr,
5350 "update_gradients"));
5351# ifdef DEBUG
5352 this->gradients_quad_submitted = true;
5353# endif
5354
5355 const std::size_t nqp_d = this->n_quadrature_points * dim;
5356 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5357
5358 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5359 {
5360 const VectorizedArrayType JxW_jac = this->J_value[0] *
5361 this->quadrature_weights[q_point] *
5362 this->normal_x_jacobian[0][dim - 1];
5363 for (unsigned int comp = 0; comp < n_components; ++comp)
5364 {
5365 for (unsigned int d = 0; d < dim - 1; ++d)
5366 gradients[comp * nqp_d + d] = VectorizedArrayType();
5367 if constexpr (n_components == 1)
5368 gradients[dim - 1] = grad_in * JxW_jac;
5369 else
5370 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5371 }
5372 }
5373 else
5374 {
5375 const unsigned int index =
5376 this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5378 this->normal_x_jacobian[index];
5379 const VectorizedArrayType JxW =
5380 (this->cell_type <= internal::MatrixFreeFunctions::affine) ?
5381 this->J_value[index] * this->quadrature_weights[q_point] :
5382 this->J_value[index];
5383 for (unsigned int comp = 0; comp < n_components; ++comp)
5384 for (unsigned int d = 0; d < dim; ++d)
5385 if constexpr (n_components == 1)
5386 gradients[d] = (grad_in * JxW) * jac[d];
5387 else
5388 gradients[comp * nqp_d + d] = (grad_in[comp] * JxW) * jac[d];
5389 }
5390}
5391
5392
5393
5394template <int dim,
5395 int n_components_,
5396 typename Number,
5397 bool is_face,
5398 typename VectorizedArrayType>
5399inline DEAL_II_ALWAYS_INLINE void
5401 submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
5402{
5403# ifdef DEBUG
5404 Assert(this->is_reinitialized, ExcNotInitialized());
5405# endif
5406 AssertIndexRange(q_point, this->n_quadrature_points);
5407 Assert(this->J_value != nullptr,
5409 "update_hessians"));
5410 Assert(this->jacobian != nullptr,
5412 "update_hessians"));
5413# ifdef DEBUG
5414 this->hessians_quad_submitted = true;
5415# endif
5416
5417 // compute hessian_unit = J^T * hessian_in(u) * J
5418 const std::size_t nqp = this->n_quadrature_points;
5419 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5420 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5421 {
5422 const VectorizedArrayType JxW =
5423 this->J_value[0] * this->quadrature_weights[q_point];
5424
5425 // diagonal part
5426 for (unsigned int d = 0; d < dim; ++d)
5427 {
5428 const auto jac_d = this->jacobian[0][d][d];
5429 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5430 for (unsigned int comp = 0; comp < n_components; ++comp)
5431 if constexpr (n_components == 1)
5432 this->hessians_quad[d * nqp + q_point] =
5433 hessian_in[d][d] * factor;
5434 else
5435 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5436 hessian_in[comp][d][d] * factor;
5437 }
5438
5439 // off diagonal part
5440 for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5441 for (unsigned int e = 0; e < d; ++e, ++off_dia)
5442 {
5443 const auto jac_d = this->jacobian[0][d][d];
5444 const auto jac_e = this->jacobian[0][e][e];
5445 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5446 for (unsigned int comp = 0; comp < n_components; ++comp)
5447 if constexpr (n_components == 1)
5448 this->hessians_quad[off_dia * nqp + q_point] =
5449 (hessian_in[d][e] + hessian_in[e][d]) * factor;
5450 else
5451 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5452 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5453 }
5454 }
5455 // cell with general Jacobian, but constant within the cell
5456 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5457 {
5458 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
5459 const VectorizedArrayType JxW =
5460 this->J_value[0] * this->quadrature_weights[q_point];
5461 for (unsigned int comp = 0; comp < n_components; ++comp)
5462 {
5464 if constexpr (n_components == 1)
5466 else
5468
5469 // 1. tmp = hessian(u) * J
5470 VectorizedArrayType tmp[dim][dim];
5471 for (unsigned int i = 0; i < dim; ++i)
5472 for (unsigned int j = 0; j < dim; ++j)
5473 {
5474 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5475 for (unsigned int k = 1; k < dim; ++k)
5476 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5477 }
5478
5479 // 2. hessian_unit = J^T * tmp
5480 VectorizedArrayType tmp2[dim][dim];
5481 for (unsigned int i = 0; i < dim; ++i)
5482 for (unsigned int j = 0; j < dim; ++j)
5483 {
5484 tmp2[i][j] = jac[0][i] * tmp[0][j];
5485 for (unsigned int k = 1; k < dim; ++k)
5486 tmp2[i][j] += jac[k][i] * tmp[k][j];
5487 }
5488
5489 // diagonal part
5490 for (unsigned int d = 0; d < dim; ++d)
5491 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5492 tmp2[d][d] * JxW;
5493
5494 // off diagonal part
5495 for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5496 for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5497 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5498 (tmp2[d][e] + tmp2[e][d]) * JxW;
5499 }
5500 }
5501 else
5502 {
5503 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
5504 const VectorizedArrayType JxW = this->J_value[q_point];
5505 const auto &jac_grad = this->jacobian_gradients[q_point];
5506 for (unsigned int comp = 0; comp < n_components; ++comp)
5507 {
5509 if constexpr (n_components == 1)
5511 else
5513
5514 // 1. tmp = hessian(u) * J
5515 VectorizedArrayType tmp[dim][dim];
5516 for (unsigned int i = 0; i < dim; ++i)
5517 for (unsigned int j = 0; j < dim; ++j)
5518 {
5519 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5520 for (unsigned int k = 1; k < dim; ++k)
5521 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5522 }
5523
5524 // 2. hessian_unit = J^T * tmp
5525 VectorizedArrayType tmp2[dim][dim];
5526 for (unsigned int i = 0; i < dim; ++i)
5527 for (unsigned int j = 0; j < dim; ++j)
5528 {
5529 tmp2[i][j] = jac[0][i] * tmp[0][j];
5530 for (unsigned int k = 1; k < dim; ++k)
5531 tmp2[i][j] += jac[k][i] * tmp[k][j];
5532 }
5533
5534 // diagonal part
5535 for (unsigned int d = 0; d < dim; ++d)
5536 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5537 tmp2[d][d] * JxW;
5538
5539 // off diagonal part
5540 for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5541 for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5542 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5543 (tmp2[d][e] + tmp2[e][d]) * JxW;
5544
5545 // 3. gradient_unit = J' * hessian
5546 for (unsigned int d = 0; d < dim; ++d)
5547 {
5548 VectorizedArrayType sum = 0;
5549 for (unsigned int e = 0; e < dim; ++e)
5550 sum += hessian_c[e][e] * jac_grad[e][d];
5551 for (unsigned int e = 0, count = dim; e < dim; ++e)
5552 for (unsigned int f = e + 1; f < dim; ++f, ++count)
5553 sum +=
5554 (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][d];
5555 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5556 d] = sum * JxW;
5557 }
5558 }
5559 }
5560}
5561
5562
5563
5564template <int dim,
5565 int n_components_,
5566 typename Number,
5567 bool is_face,
5568 typename VectorizedArrayType>
5569inline DEAL_II_ALWAYS_INLINE void
5572 const unsigned int q_point)
5573{
5574# ifdef DEBUG
5575 Assert(this->is_reinitialized, ExcNotInitialized());
5576# endif
5577 AssertIndexRange(q_point, this->n_quadrature_points);
5578 Assert(this->J_value != nullptr,
5580 "update_hessians"));
5581 Assert(this->jacobian != nullptr,
5583 "update_hessians"));
5584# ifdef DEBUG
5585 this->hessians_quad_submitted = true;
5586# endif
5587
5588 // compute hessian_unit = J^T * hessian_in(u) * J
5589 const std::size_t nqp = this->n_quadrature_points;
5590 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5591 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5592 {
5593 const VectorizedArrayType JxW =
5594 this->J_value[0] * this->quadrature_weights[q_point];
5595
5596 const auto nxj = this->normal_x_jacobian[0];
5597
5598 // diagonal part
5599 for (unsigned int d = 0; d < dim; ++d)
5600 {
5601 const auto nxj_d = nxj[d];
5602 const VectorizedArrayType factor = nxj_d * nxj_d * JxW;
5603 for (unsigned int comp = 0; comp < n_components; ++comp)
5604 if constexpr (n_components == 1)
5605 this->hessians_quad[d * nqp + q_point] =
5606 normal_hessian_in * factor;
5607 else
5608 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5609 normal_hessian_in[comp] * factor;
5610 }
5611
5612 // off diagonal part
5613 for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5614 for (unsigned int e = 0; e < d; ++e, ++off_dia)
5615 {
5616 const auto jac_d = nxj[d];
5617 const auto jac_e = nxj[e];
5618 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5619 for (unsigned int comp = 0; comp < n_components; ++comp)
5620 if constexpr (n_components == 1)
5621 this->hessians_quad[off_dia * nqp + q_point] =
5622 2. * normal_hessian_in * factor;
5623 else
5624 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5625 2. * normal_hessian_in[comp] * factor;
5626 }
5627 }
5628 else
5629 {
5630 const auto normal = this->normal_vector(q_point);
5631 const auto normal_projector = outer_product(normal, normal);
5632 if constexpr (n_components == 1)
5633 submit_hessian(normal_hessian_in * normal_projector, q_point);
5634 else
5635 {
5636 hessian_type tmp;
5637 for (unsigned int comp = 0; comp < n_components; ++comp)
5639 submit_hessian(tmp, q_point);
5640 }
5641 }
5642}
5643
5644
5645
5646template <int dim,
5647 int n_components_,
5648 typename Number,
5649 bool is_face,
5650 typename VectorizedArrayType>
5651inline typename FEEvaluationBase<dim,
5653 Number,
5654 is_face,
5655 VectorizedArrayType>::value_type
5657 integrate_value() const
5658{
5659# ifdef DEBUG
5660 Assert(this->is_reinitialized, ExcNotInitialized());
5661 Assert(this->values_quad_submitted == true,
5663# endif
5664
5666 const std::size_t nqp = this->n_quadrature_points;
5667 for (unsigned int q = 0; q < nqp; ++q)
5668 for (unsigned int comp = 0; comp < n_components; ++comp)
5669 return_value[comp] += this->values_quad[comp * nqp + q];
5670 if constexpr (n_components == 1)
5671 return return_value[0];
5672 else
5673 return return_value;
5674}
5675
5676
5677
5678template <int dim,
5679 int n_components_,
5680 typename Number,
5681 bool is_face,
5682 typename VectorizedArrayType>
5683template <int, typename>
5684inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5686 get_divergence(const unsigned int q_point) const
5687{
5688 static_assert(n_components == dim,
5689 "Do not try to modify the default template parameters used for"
5690 " selectively enabling this function via std::enable_if!");
5691
5692# ifdef DEBUG
5693 Assert(this->gradients_quad_initialized == true,
5695# endif
5696 AssertIndexRange(q_point, this->n_quadrature_points);
5697 Assert(this->jacobian != nullptr,
5699 "update_gradients"));
5700
5701 VectorizedArrayType divergence;
5702 const std::size_t nqp = this->n_quadrature_points;
5703
5704 if (dim > 1 &&
5705 this->data->element_type ==
5707 {
5708 VectorizedArrayType inv_det =
5709 (!is_face &&
5710 this->cell_type == internal::MatrixFreeFunctions::cartesian) ?
5711 this->jacobian[0][0][0] *
5712 ((dim == 2) ? this->jacobian[0][1][1] :
5713 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
5714 determinant(this->jacobian[this->cell_type >
5715 internal::MatrixFreeFunctions::affine ?
5716 q_point :
5717 0]);
5718
5719 // on faces in 2d, the determinant has the wrong sign due to ordering of
5720 // derivatives
5721 if (is_face && dim == 2 && this->get_face_no() < 2)
5722 inv_det = -inv_det;
5723
5724 // div * det(J^-1)
5725 divergence = this->gradients_quad[q_point * dim];
5726 for (unsigned int d = 1; d < dim; ++d)
5727 divergence += this->gradients_quad[(d * nqp + q_point) * dim + d];
5728 divergence *= inv_det;
5729 }
5730 else
5731 {
5732 if (!is_face &&
5734 {
5735 // Cartesian cell
5736 divergence =
5737 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
5738 for (unsigned int d = 1; d < dim; ++d)
5739 divergence += this->gradients_quad[(d * nqp + q_point) * dim + d] *
5740 this->jacobian[0][d][d];
5741 }
5742 else
5743 {
5744 // cell with general/constant Jacobian
5746 this->cell_type == internal::MatrixFreeFunctions::general ?
5747 this->jacobian[q_point] :
5748 this->jacobian[0];
5749 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
5750 for (unsigned int e = 1; e < dim; ++e)
5751 divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
5752 for (unsigned int d = 1; d < dim; ++d)
5753 for (unsigned int e = 0; e < dim; ++e)
5754 divergence +=
5755 jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim + e];
5756 }
5757 }
5758 return divergence;
5759}
5760
5761
5762
5763template <int dim,
5764 int n_components_,
5765 typename Number,
5766 bool is_face,
5767 typename VectorizedArrayType>
5768template <int, typename>
5771 get_symmetric_gradient(const unsigned int q_point) const
5772{
5773 static_assert(n_components == dim,
5774 "Do not try to modify the default template parameters used for"
5775 " selectively enabling this function via std::enable_if!");
5776
5777 // copy from generic function into dim-specialization function
5778 const auto grad = get_gradient(q_point);
5779 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5780 VectorizedArrayType half = Number(0.5);
5781 for (unsigned int d = 0; d < dim; ++d)
5782 symmetrized[d] = grad[d][d];
5783 switch (dim)
5784 {
5785 case 1:
5786 break;
5787 case 2:
5788 symmetrized[2] = grad[0][1] + grad[1][0];
5789 symmetrized[2] *= half;
5790 break;
5791 case 3:
5792 symmetrized[3] = grad[0][1] + grad[1][0];
5793 symmetrized[3] *= half;
5794 symmetrized[4] = grad[0][2] + grad[2][0];
5795 symmetrized[4] *= half;
5796 symmetrized[5] = grad[1][2] + grad[2][1];
5797 symmetrized[5] *= half;
5798 break;
5799 default:
5801 }
5803}
5804
5805
5806
5807template <int dim,
5808 int n_components_,
5809 typename Number,
5810 bool is_face,
5811 typename VectorizedArrayType>
5812template <int, typename>
5814 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
5816 get_curl(const unsigned int q_point) const
5817{
5818 static_assert(dim > 1 && n_components == dim,
5819 "Do not try to modify the default template parameters used for"
5820 " selectively enabling this function via std::enable_if!");
5821
5822 // copy from generic function into dim-specialization function
5823 const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
5824 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
5825 switch (dim)
5826 {
5827 case 2:
5828 curl[0] = grad[1][0] - grad[0][1];
5829 break;
5830 case 3:
5831 curl[0] = grad[2][1] - grad[1][2];
5832 curl[1] = grad[0][2] - grad[2][0];
5833 curl[2] = grad[1][0] - grad[0][1];
5834 break;
5835 default:
5837 }
5838 return curl;
5839}
5840
5841
5842
5843template <int dim,
5844 int n_components_,
5845 typename Number,
5846 bool is_face,
5847 typename VectorizedArrayType>
5848template <int, typename>
5849inline DEAL_II_ALWAYS_INLINE void
5851 submit_divergence(const VectorizedArrayType div_in,
5852 const unsigned int q_point)
5853{
5854 static_assert(n_components == dim,
5855 "Do not try to modify the default template parameters used for"
5856 " selectively enabling this function via std::enable_if!");
5857
5858# ifdef DEBUG
5859 Assert(this->is_reinitialized, ExcNotInitialized());
5860# endif
5861 AssertIndexRange(q_point, this->n_quadrature_points);
5862 Assert(this->J_value != nullptr,
5864 "update_gradients"));
5865 Assert(this->jacobian != nullptr,
5867 "update_gradients"));
5868# ifdef DEBUG
5869 this->gradients_quad_submitted = true;
5870# endif
5871
5872 const std::size_t nqp_d = this->n_quadrature_points * dim;
5873 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5874
5875 if (this->data->element_type ==
5877 {
5878 // General cell
5879
5880 // Derivatives are reordered for faces. Need to take this into account
5881 // and 1/inv_det != J_value for faces
5882 const VectorizedArrayType fac =
5883 (!is_face) ?
5884 this->quadrature_weights[q_point] * div_in :
5885 (this->cell_type > internal::MatrixFreeFunctions::affine ?
5886 this->J_value[q_point] :
5887 this->J_value[0] * this->quadrature_weights[q_point]) *
5888 div_in *
5890 this->jacobian[this->cell_type >
5891 internal::MatrixFreeFunctions::affine ?
5892 q_point :
5893 0]) *
5894 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
5895
5896 for (unsigned int d = 0; d < dim; ++d)
5897 {
5898 for (unsigned int e = 0; e < dim; ++e)
5899 gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
5900 }
5901 this->divergence_is_requested = true;
5902 }
5903 else
5904 {
5905 if (!is_face &&
5907 {
5908 const VectorizedArrayType fac =
5909 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
5910 for (unsigned int d = 0; d < dim; ++d)
5911 {
5912 const VectorizedArrayType jac_dd = this->jacobian[0][d][d];
5913 for (unsigned int e = 0; e < dim; ++e)
5914 gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
5915 }
5916 }
5917 else
5918 {
5920 this->cell_type == internal::MatrixFreeFunctions::general ?
5921 this->jacobian[q_point] :
5922 this->jacobian[0];
5923 const VectorizedArrayType fac =
5924 (this->cell_type == internal::MatrixFreeFunctions::general ?
5925 this->J_value[q_point] :
5926 this->J_value[0] * this->quadrature_weights[q_point]) *
5927 div_in;
5928 for (unsigned int d = 0; d < dim; ++d)
5929 {
5930 for (unsigned int e = 0; e < dim; ++e)
5931 gradients[d * nqp_d + e] = jac[d][e] * fac;
5932 }
5933 }
5934 }
5935}
5936
5937
5938
5939template <int dim,
5940 int n_components_,
5941 typename Number,
5942 bool is_face,
5943 typename VectorizedArrayType>
5944template <int, typename>
5945inline DEAL_II_ALWAYS_INLINE void
5949 const unsigned int q_point)
5950{
5951 static_assert(n_components == dim,
5952 "Do not try to modify the default template parameters used for"
5953 " selectively enabling this function via std::enable_if!");
5954
5956 this->data->element_type !=
5959
5960 // could have used base class operator, but that involves some overhead
5961 // which is inefficient. it is nice to have the symmetric tensor because
5962 // that saves some operations
5963# ifdef DEBUG
5964 Assert(this->is_reinitialized, ExcNotInitialized());
5965# endif
5966 AssertIndexRange(q_point, this->n_quadrature_points);
5967 Assert(this->J_value != nullptr,
5969 "update_gradients"));
5970 Assert(this->jacobian != nullptr,
5972 "update_gradients"));
5973# ifdef DEBUG
5974 this->gradients_quad_submitted = true;
5975# endif
5976
5977 const std::size_t nqp_d = this->n_quadrature_points * dim;
5978 VectorizedArrayType *gradients = this->gradients_quad + dim * q_point;
5979 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5980 {
5981 const VectorizedArrayType JxW =
5982 this->J_value[0] * this->quadrature_weights[q_point];
5983 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
5984 for (unsigned int d = 0; d < dim; ++d)
5985 gradients[d * nqp_d + d] =
5986 (sym_grad.access_raw_entry(d) * JxW * jac[d][d]);
5987 for (unsigned int e = 0, counter = dim; e < dim; ++e)
5988 for (unsigned int d = e + 1; d < dim; ++d, ++counter)
5989 {
5990 const VectorizedArrayType value =
5991 sym_grad.access_raw_entry(counter) * JxW;
5992 gradients[e * nqp_d + d] = value * jac[d][d];
5993 gradients[d * nqp_d + e] = value * jac[e][e];
5994 }
5995 }
5996 // general/affine cell type
5997 else
5998 {
5999 const VectorizedArrayType JxW =
6000 this->cell_type == internal::MatrixFreeFunctions::general ?
6001 this->J_value[q_point] :
6002 this->J_value[0] * this->quadrature_weights[q_point];
6004 this->cell_type == internal::MatrixFreeFunctions::general ?
6005 this->jacobian[q_point] :
6006 this->jacobian[0];
6007 VectorizedArrayType weighted[dim][dim];
6008 for (unsigned int i = 0; i < dim; ++i)
6009 weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6010 for (unsigned int i = 0, counter = dim; i < dim; ++i)
6011 for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6012 {
6013 const VectorizedArrayType value =
6014 sym_grad.access_raw_entry(counter) * JxW;
6015 weighted[i][j] = value;
6016 weighted[j][i] = value;
6017 }
6018 for (unsigned int comp = 0; comp < dim; ++comp)
6019 for (unsigned int d = 0; d < dim; ++d)
6020 {
6021 VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6022 for (unsigned int e = 1; e < dim; ++e)
6023 new_val += jac[e][d] * weighted[comp][e];
6024 gradients[comp * nqp_d + d] = new_val;
6025 }
6026 }
6027}
6028
6029
6030
6031template <int dim,
6032 int n_components_,
6033 typename Number,
6034 bool is_face,
6035 typename VectorizedArrayType>
6036template <int, typename>
6037inline DEAL_II_ALWAYS_INLINE void
6040 const unsigned int q_point)
6041{
6042 static_assert(n_components == dim,
6043 "Do not try to modify the default template parameters used for"
6044 " selectively enabling this function via std::enable_if!");
6045
6047 switch (dim)
6048 {
6049 case 2:
6050 grad[1][0] = curl[0];
6051 grad[0][1] = -curl[0];
6052 break;
6053 case 3:
6054 grad[2][1] = curl[0];
6055 grad[1][2] = -curl[0];
6056 grad[0][2] = curl[1];
6057 grad[2][0] = -curl[1];
6058 grad[1][0] = curl[2];
6059 grad[0][1] = -curl[2];
6060 break;
6061 default:
6063 }
6064 submit_gradient(grad, q_point);
6065}
6066
6067
6068
6069/*-------------------------- FEEvaluation -----------------------------------*/
6070
6071
6072template <int dim,
6073 int fe_degree,
6074 int n_q_points_1d,
6075 int n_components_,
6076 typename Number,
6077 typename VectorizedArrayType>
6078inline FEEvaluation<dim,
6079 fe_degree,
6080 n_q_points_1d,
6082 Number,
6083 VectorizedArrayType>::
6085 const unsigned int fe_no,
6086 const unsigned int quad_no,
6087 const unsigned int first_selected_component,
6088 const unsigned int active_fe_index,
6089 const unsigned int active_quad_index)
6090 : BaseClass(matrix_free,
6091 fe_no,
6092 first_selected_component,
6093 quad_no,
6094 fe_degree,
6095 static_n_q_points,
6096 true /*note: this is not a face*/,
6097 active_fe_index,
6098 active_quad_index,
6099 numbers::invalid_unsigned_int /*face_type*/)
6100 , dofs_per_component(this->data->dofs_per_component_on_cell)
6101 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6102 , n_q_points(this->data->n_q_points)
6103{
6104 check_template_arguments(fe_no, 0);
6105}
6106
6107
6108
6109template <int dim,
6110 int fe_degree,
6111 int n_q_points_1d,
6112 int n_components_,
6113 typename Number,
6114 typename VectorizedArrayType>
6115inline FEEvaluation<dim,
6116 fe_degree,
6117 n_q_points_1d,
6119 Number,
6120 VectorizedArrayType>::
6122 const std::pair<unsigned int, unsigned int> &range,
6123 const unsigned int dof_no,
6124 const unsigned int quad_no,
6125 const unsigned int first_selected_component)
6126 : FEEvaluation(matrix_free,
6127 dof_no,
6128 quad_no,
6129 first_selected_component,
6130 matrix_free.get_cell_active_fe_index(range))
6131{}
6132
6133
6134
6135template <int dim,
6136 int fe_degree,
6137 int n_q_points_1d,
6138 int n_components_,
6139 typename Number,
6140 typename VectorizedArrayType>
6141inline FEEvaluation<dim,
6142 fe_degree,
6143 n_q_points_1d,
6145 Number,
6146 VectorizedArrayType>::
6147 FEEvaluation(const Mapping<dim> &mapping,
6148 const FiniteElement<dim> &fe,
6149 const Quadrature<1> &quadrature,
6150 const UpdateFlags update_flags,
6151 const unsigned int first_selected_component)
6152 : BaseClass(mapping,
6153 fe,
6154 quadrature,
6155 update_flags,
6156 first_selected_component,
6157 nullptr)
6158 , dofs_per_component(this->data->dofs_per_component_on_cell)
6159 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6160 , n_q_points(this->data->n_q_points)
6161{
6162 check_template_arguments(numbers::invalid_unsigned_int, 0);
6163}
6164
6165
6166
6167template <int dim,
6168 int fe_degree,
6169 int n_q_points_1d,
6170 int n_components_,
6171 typename Number,
6172 typename VectorizedArrayType>
6173inline FEEvaluation<dim,
6174 fe_degree,
6175 n_q_points_1d,
6177 Number,
6178 VectorizedArrayType>::
6180 const Quadrature<1> &quadrature,
6181 const UpdateFlags update_flags,
6182 const unsigned int first_selected_component)
6183 : BaseClass(StaticMappingQ1<dim>::mapping,
6184 fe,
6185 quadrature,
6186 update_flags,
6187 first_selected_component,
6188 nullptr)
6189 , dofs_per_component(this->data->dofs_per_component_on_cell)
6190 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6191 , n_q_points(this->data->n_q_points)
6192{
6193 check_template_arguments(numbers::invalid_unsigned_int, 0);
6194}
6195
6196
6197
6198template <int dim,
6199 int fe_degree,
6200 int n_q_points_1d,
6201 int n_components_,
6202 typename Number,
6203 typename VectorizedArrayType>
6204inline FEEvaluation<dim,
6205 fe_degree,
6206 n_q_points_1d,
6208 Number,
6209 VectorizedArrayType>::
6212 const unsigned int first_selected_component)
6213 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6214 fe,
6215 other.mapped_geometry->get_quadrature(),
6216 other.mapped_geometry->get_fe_values().get_update_flags(),
6217 first_selected_component,
6218 &other)
6219 , dofs_per_component(this->data->dofs_per_component_on_cell)
6220 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6221 , n_q_points(this->data->n_q_points)
6222{
6223 check_template_arguments(numbers::invalid_unsigned_int, 0);
6224}
6225
6226
6227
6228template <int dim,
6229 int fe_degree,
6230 int n_q_points_1d,
6231 int n_components_,
6232 typename Number,
6233 typename VectorizedArrayType>
6234inline FEEvaluation<dim,
6235 fe_degree,
6236 n_q_points_1d,
6238 Number,
6239 VectorizedArrayType>::FEEvaluation(const FEEvaluation
6240 &other)
6241 : BaseClass(other)
6242 , dofs_per_component(this->data->dofs_per_component_on_cell)
6243 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6244 , n_q_points(this->data->n_q_points)
6245{
6246 check_template_arguments(numbers::invalid_unsigned_int, 0);
6247}
6248
6249
6250
6251template <int dim,
6252 int fe_degree,
6253 int n_q_points_1d,
6254 int n_components_,
6255 typename Number,
6256 typename VectorizedArrayType>
6257inline FEEvaluation<dim,
6258 fe_degree,
6259 n_q_points_1d,
6261 Number,
6262 VectorizedArrayType> &
6263FEEvaluation<dim,
6264 fe_degree,
6265 n_q_points_1d,
6267 Number,
6268 VectorizedArrayType>::operator=(const FEEvaluation &other)
6269{
6270 BaseClass::operator=(other);
6271 check_template_arguments(numbers::invalid_unsigned_int, 0);
6272 return *this;
6273}
6274
6275
6276
6277template <int dim,
6278 int fe_degree,
6279 int n_q_points_1d,
6280 int n_components_,
6281 typename Number,
6282 typename VectorizedArrayType>
6283inline void
6284FEEvaluation<dim,
6285 fe_degree,
6286 n_q_points_1d,
6288 Number,
6289 VectorizedArrayType>::
6290 check_template_arguments(const unsigned int dof_no,
6291 const unsigned int first_selected_component)
6292{
6293 (void)dof_no;
6294 (void)first_selected_component;
6295
6296 Assert(
6297 this->data->dofs_per_component_on_cell > 0,
6298 ExcMessage(
6299 "There is nothing useful you can do with an FEEvaluation object with "
6300 "FE_Nothing, i.e., without DoFs! If you have passed to "
6301 "MatrixFree::reinit() a collection of finite elements also containing "
6302 "FE_Nothing, please check - before creating FEEvaluation - the category "
6303 "of the current range by calling either "
6304 "MatrixFree::get_cell_range_category(range) or "
6305 "MatrixFree::get_face_range_category(range). The returned category "
6306 "is the index of the active FE, which you can use to exclude "
6307 "FE_Nothing."));
6308
6309# ifdef DEBUG
6310 // print error message when the dimensions do not match. Propose a possible
6311 // fix
6312 if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
6313 static_cast<unsigned int>(fe_degree) !=
6314 this->data->data.front().fe_degree) ||
6315 n_q_points != this->n_quadrature_points)
6316 {
6317 std::string message =
6318 "-------------------------------------------------------\n";
6319 message += "Illegal arguments in constructor/wrong template arguments!\n";
6320 message += " Called --> FEEvaluation<dim,";
6321 message += Utilities::int_to_string(fe_degree) + ",";
6322 message += Utilities::int_to_string(n_q_points_1d);
6323 message += "," + Utilities::int_to_string(n_components);
6324 message += ",Number>(data";
6325 if (first_selected_component != numbers::invalid_unsigned_int)
6326 {
6327 message += ", " + Utilities::int_to_string(dof_no) + ", ";
6328 message += Utilities::int_to_string(this->quad_no) + ", ";
6329 message += Utilities::int_to_string(first_selected_component);
6330 }
6331 message += ")\n";
6332
6333 // check whether some other vector component has the correct number of
6334 // points
6339 {
6340 if (static_cast<unsigned int>(fe_degree) ==
6341 this->data->data.front().fe_degree)
6342 {
6343 proposed_dof_comp = dof_no;
6344 proposed_fe_comp = first_selected_component;
6345 }
6346 else
6347 for (unsigned int no = 0; no < this->matrix_free->n_components();
6348 ++no)
6349 for (unsigned int nf = 0;
6350 nf < this->matrix_free->n_base_elements(no);
6351 ++nf)
6352 if (this->matrix_free
6353 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6354 .data.front()
6355 .fe_degree == static_cast<unsigned int>(fe_degree))
6356 {
6359 break;
6360 }
6361 if (n_q_points ==
6362 this->mapping_data->descriptor[this->active_quad_index]
6363 .n_q_points)
6364 proposed_quad_comp = this->quad_no;
6365 else
6366 for (unsigned int no = 0;
6367 no < this->matrix_free->get_mapping_info().cell_data.size();
6368 ++no)
6369 if (this->matrix_free->get_mapping_info()
6370 .cell_data[no]
6371 .descriptor[this->active_quad_index]
6372 .n_q_points == n_q_points)
6373 {
6375 break;
6376 }
6377 }
6380 {
6381 if (proposed_dof_comp != first_selected_component)
6382 message += "Wrong vector component selection:\n";
6383 else
6384 message += "Wrong quadrature formula selection:\n";
6385 message += " Did you mean FEEvaluation<dim,";
6386 message += Utilities::int_to_string(fe_degree) + ",";
6387 message += Utilities::int_to_string(n_q_points_1d);
6388 message += "," + Utilities::int_to_string(n_components);
6389 message += ",Number>(data";
6391 {
6392 message +=
6396 }
6397 message += ")?\n";
6398 std::string correct_pos;
6400 correct_pos = " ^ ";
6401 else
6402 correct_pos = " ";
6403 if (proposed_quad_comp != this->quad_no)
6404 correct_pos += " ^ ";
6405 else
6406 correct_pos += " ";
6407 if (proposed_fe_comp != first_selected_component)
6408 correct_pos += " ^\n";
6409 else
6410 correct_pos += " \n";
6411 message += " " +
6413 }
6414 // ok, did not find the numbers specified by the template arguments in
6415 // the given list. Suggest correct template arguments
6416 const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
6417 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6418 message += "Wrong template arguments:\n";
6419 message += " Did you mean FEEvaluation<dim,";
6420 message +=
6421 Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
6423 message += "," + Utilities::int_to_string(n_components);
6424 message += ",Number>(data";
6426 {
6427 message += ", " + Utilities::int_to_string(dof_no) + ", ";
6428 message += Utilities::int_to_string(this->quad_no);
6429 message += ", " + Utilities::int_to_string(first_selected_component);
6430 }
6431 message += ")?\n";
6432 std::string correct_pos;
6433 if (this->data->data.front().fe_degree !=
6434 static_cast<unsigned int>(fe_degree))
6435 correct_pos = " ^";
6436 else
6437 correct_pos = " ";
6438 if (proposed_n_q_points_1d != n_q_points_1d)
6439 correct_pos += " ^\n";
6440 else
6441 correct_pos += " \n";
6442 message += " " + correct_pos;
6443
6444 Assert(static_cast<unsigned int>(fe_degree) ==
6445 this->data->data.front().fe_degree &&
6446 n_q_points == this->n_quadrature_points,
6448 }
6451 n_q_points,
6452 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6453# endif
6454}
6455
6456
6457
6458template <int dim,
6459 int fe_degree,
6460 int n_q_points_1d,
6461 int n_components_,
6462 typename Number,
6463 typename VectorizedArrayType>
6464inline void
6465FEEvaluation<dim,
6466 fe_degree,
6467 n_q_points_1d,
6469 Number,
6470 VectorizedArrayType>::reinit(const unsigned int cell_index)
6471{
6472 Assert(this->matrix_free != nullptr,
6473 ExcMessage("FEEvaluation was initialized without a matrix-free object."
6474 " Integer indexing is not possible."));
6475
6476 Assert(this->dof_info != nullptr, ExcNotInitialized());
6477 Assert(this->mapping_data != nullptr, ExcNotInitialized());
6478 this->cell = cell_index;
6479 this->cell_type =
6480 this->matrix_free->get_mapping_info().get_cell_type(cell_index);
6481
6482 const unsigned int offsets =
6483 this->mapping_data->data_index_offsets[cell_index];
6484 this->jacobian = &this->mapping_data->jacobians[0][offsets];
6485 this->J_value = &this->mapping_data->JxW_values[offsets];
6486 if (!this->mapping_data->jacobian_gradients[0].empty())
6487 {
6488 this->jacobian_gradients =
6489 this->mapping_data->jacobian_gradients[0].data() + offsets;
6490 this->jacobian_gradients_non_inverse =
6491 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
6492 }
6493
6494 if (this->matrix_free->n_active_entries_per_cell_batch(this->cell) == n_lanes)
6495 {
6497 for (unsigned int i = 0; i < n_lanes; ++i)
6498 this->cell_ids[i] = cell_index * n_lanes + i;
6499 }
6500 else
6501 {
6502 unsigned int i = 0;
6503 for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell);
6504 ++i)
6505 this->cell_ids[i] = cell_index * n_lanes + i;
6506 for (; i < n_lanes; ++i)
6507 this->cell_ids[i] = numbers::invalid_unsigned_int;
6508 }
6509
6510 if (this->mapping_data->quadrature_points.empty() == false)
6511 this->quadrature_points =
6512 &this->mapping_data->quadrature_points
6513 [this->mapping_data->quadrature_point_offsets[this->cell]];
6514
6515# ifdef DEBUG
6516 this->is_reinitialized = true;
6517 this->dof_values_initialized = false;
6518 this->values_quad_initialized = false;
6519 this->gradients_quad_initialized = false;
6520 this->hessians_quad_initialized = false;
6521# endif
6522}
6523
6524
6525
6526template <int dim,
6527 int fe_degree,
6528 int n_q_points_1d,
6529 int n_components_,
6530 typename Number,
6531 typename VectorizedArrayType>
6532inline void
6533FEEvaluation<dim,
6534 fe_degree,
6535 n_q_points_1d,
6537 Number,
6538 VectorizedArrayType>::reinit(const std::array<unsigned int,
6539 n_lanes> &cell_ids)
6540{
6541 Assert(this->dof_info != nullptr, ExcNotInitialized());
6542 Assert(this->mapping_data != nullptr, ExcNotInitialized());
6543
6544 this->cell = numbers::invalid_unsigned_int;
6545 this->cell_ids = cell_ids;
6546
6547 // determine type of cell batch
6549
6550 for (unsigned int v = 0; v < n_lanes; ++v)
6551 {
6552 const unsigned int cell_index = cell_ids[v];
6553
6555 continue;
6556
6557 this->cell_type =
6558 std::max(this->cell_type,
6559 this->matrix_free->get_mapping_info().get_cell_type(
6560 cell_index / n_lanes));
6561 }
6562
6563 // allocate memory for internal data storage
6564 if (this->mapped_geometry == nullptr)
6565 this->mapped_geometry =
6566 std::make_shared<internal::MatrixFreeFunctions::
6567 MappingDataOnTheFly<dim, VectorizedArrayType>>();
6568
6569 auto &mapping_storage = this->mapped_geometry->get_data_storage();
6570
6571 auto &this_jacobian_data = mapping_storage.jacobians[0];
6572 auto &this_J_value_data = mapping_storage.JxW_values;
6573 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
6575 mapping_storage.jacobian_gradients_non_inverse[0];
6576 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
6577
6579 {
6580 if (this_jacobian_data.size() != 2)
6581 this_jacobian_data.resize_fast(2);
6582
6583 if (this_J_value_data.size() != 1)
6584 this_J_value_data.resize_fast(1);
6585
6586 const auto &update_flags_cells =
6587 this->matrix_free->get_mapping_info().update_flags_cells;
6588
6589 if (update_flags_cells & update_jacobian_grads &&
6590 this_jacobian_gradients_data.size() != 1)
6591 {
6592 this_jacobian_gradients_data.resize_fast(1);
6594 }
6595
6596 if (update_flags_cells & update_quadrature_points &&
6597 this_quadrature_points_data.size() != 1)
6598 this_quadrature_points_data.resize_fast(1);
6599 }
6600 else
6601 {
6602 if (this_jacobian_data.size() != this->n_quadrature_points)
6603 this_jacobian_data.resize_fast(this->n_quadrature_points);
6604
6605 if (this_J_value_data.size() != this->n_quadrature_points)
6606 this_J_value_data.resize_fast(this->n_quadrature_points);
6607
6608 const auto &update_flags_cells =
6609 this->matrix_free->get_mapping_info().update_flags_cells;
6610
6611 if (update_flags_cells & update_jacobian_grads &&
6612 this_jacobian_gradients_data.size() != this->n_quadrature_points)
6613 {
6614 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
6616 this->n_quadrature_points);
6617 }
6618
6619 if (update_flags_cells & update_quadrature_points &&
6620 this_quadrature_points_data.size() != this->n_quadrature_points)
6621 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
6622 }
6623
6624 // set pointers to internal data storage
6625 this->jacobian = this_jacobian_data.data();
6626 this->J_value = this_J_value_data.data();
6627 this->jacobian_gradients = this_jacobian_gradients_data.data();
6628 this->jacobian_gradients_non_inverse =
6631
6632 // fill internal data storage lane by lane
6633 for (unsigned int v = 0; v < n_lanes; ++v)
6634 {
6635 const unsigned int cell_index = cell_ids[v];
6636
6638 continue;
6639
6640 const unsigned int cell_batch_index = cell_index / n_lanes;
6641 const unsigned int offsets =
6642 this->mapping_data->data_index_offsets[cell_batch_index];
6643 const unsigned int lane = cell_index % n_lanes;
6644
6645 if (this->cell_type <=
6647 {
6648 // case that all cells are Cartesian or affine
6649 for (unsigned int q = 0; q < 2; ++q)
6650 for (unsigned int i = 0; i < dim; ++i)
6651 for (unsigned int j = 0; j < dim; ++j)
6652 this_jacobian_data[q][i][j][v] =
6653 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
6654
6655 const unsigned int q = 0;
6656
6657 this_J_value_data[q][v] =
6658 this->mapping_data->JxW_values[offsets + q][lane];
6659
6660 const auto &update_flags_cells =
6661 this->matrix_free->get_mapping_info().update_flags_cells;
6662
6663 if (update_flags_cells & update_jacobian_grads)
6664 {
6665 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6666 for (unsigned int j = 0; j < dim; ++j)
6668 this->mapping_data
6669 ->jacobian_gradients[0][offsets + q][i][j][lane];
6670
6671 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6672 for (unsigned int j = 0; j < dim; ++j)
6674 this->mapping_data
6675 ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
6676 [lane];
6677 }
6678
6679 if (update_flags_cells & update_quadrature_points)
6680 for (unsigned int i = 0; i < dim; ++i)
6682 this->mapping_data->quadrature_points
6683 [this->mapping_data
6684 ->quadrature_point_offsets[cell_batch_index] +
6685 q][i][lane];
6686 }
6687 else
6688 {
6689 // general case that at least one cell is not Cartesian or affine
6690 const auto cell_type =
6691 this->matrix_free->get_mapping_info().get_cell_type(
6693
6694 for (unsigned int q = 0; q < this->n_quadrature_points; ++q)
6695 {
6696 const unsigned int q_src =
6697 (cell_type <=
6699 0 :
6700 q;
6701
6702 this_J_value_data[q][v] =
6703 this->mapping_data->JxW_values[offsets + q_src][lane];
6704
6705 for (unsigned int i = 0; i < dim; ++i)
6706 for (unsigned int j = 0; j < dim; ++j)
6707 this_jacobian_data[q][i][j][v] =
6708 this->mapping_data
6709 ->jacobians[0][offsets + q_src][i][j][lane];
6710
6711 const auto &update_flags_cells =
6712 this->matrix_free->get_mapping_info().update_flags_cells;
6713
6714 if (update_flags_cells & update_jacobian_grads)
6715 {
6716 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6717 for (unsigned int j = 0; j < dim; ++j)
6719 this->mapping_data
6720 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
6721
6722 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6723 for (unsigned int j = 0; j < dim; ++j)
6725 this->mapping_data
6726 ->jacobian_gradients_non_inverse[0][offsets + q_src]
6727 [i][j][lane];
6728 }
6729
6730 if (update_flags_cells & update_quadrature_points)
6731 {
6732 if (cell_type <=
6734 {
6735 // affine case: quadrature points are not available but
6736 // have to be computed from the corner point and the
6737 // Jacobian
6739 this->mapping_data->quadrature_points
6740 [this->mapping_data
6741 ->quadrature_point_offsets[cell_batch_index] +
6742 0];
6743
6745 this->mapping_data->jacobians[0][offsets + 1];
6747 for (unsigned int d = 0; d < dim; ++d)
6748 point[d] +=
6749 jac[d][d] *
6750 static_cast<Number>(
6751 this->descriptor->quadrature.point(q)[d]);
6752 else
6753 for (unsigned int d = 0; d < dim; ++d)
6754 for (unsigned int e = 0; e < dim; ++e)
6755 point[d] +=
6756 jac[d][e] *
6757 static_cast<Number>(
6758 this->descriptor->quadrature.point(q)[e]);
6759
6760 for (unsigned int i = 0; i < dim; ++i)
6761 this_quadrature_points_data[q][i][v] = point[i][lane];
6762 }
6763 else
6764 {
6765 // general case: quadrature points are available
6766 for (unsigned int i = 0; i < dim; ++i)
6768 this->mapping_data->quadrature_points
6769 [this->mapping_data
6770 ->quadrature_point_offsets[cell_batch_index] +
6771 q][i][lane];
6772 }
6773 }
6774 }
6775 }
6776 }
6777
6778# ifdef DEBUG
6779 this->is_reinitialized = true;
6780 this->dof_values_initialized = false;
6781 this->values_quad_initialized = false;
6782 this->gradients_quad_initialized = false;
6783 this->hessians_quad_initialized = false;
6784# endif
6785}
6786
6787
6788
6789template <int dim,
6790 int fe_degree,
6791 int n_q_points_1d,
6792 int n_components_,
6793 typename Number,
6794 typename VectorizedArrayType>
6795template <bool level_dof_access>
6796inline void
6797FEEvaluation<dim,
6798 fe_degree,
6799 n_q_points_1d,
6801 Number,
6802 VectorizedArrayType>::
6804{
6805 Assert(this->matrix_free == nullptr,
6806 ExcMessage("Cannot use initialization from cell iterator if "
6807 "initialized from MatrixFree object. Use variant for "
6808 "on the fly computation with arguments as for FEValues "
6809 "instead"));
6810 Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
6811 this->mapped_geometry->reinit(
6812 static_cast<typename Triangulation<dim>::cell_iterator>(cell));
6813 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6814 if (level_dof_access)
6815 cell->get_mg_dof_indices(this->local_dof_indices);
6816 else
6817 cell->get_dof_indices(this->local_dof_indices);
6818
6819# ifdef DEBUG
6820 this->is_reinitialized = true;
6821# endif
6822}
6823
6824
6825
6826template <int dim,
6827 int fe_degree,
6828 int n_q_points_1d,
6829 int n_components_,
6830 typename Number,
6831 typename VectorizedArrayType>
6832inline void
6833FEEvaluation<dim,
6834 fe_degree,
6835 n_q_points_1d,
6837 Number,
6838 VectorizedArrayType>::
6839 reinit(const typename Triangulation<dim>::cell_iterator &cell)
6840{
6841 Assert(this->matrix_free == 0,
6842 ExcMessage("Cannot use initialization from cell iterator if "
6843 "initialized from MatrixFree object. Use variant for "
6844 "on the fly computation with arguments as for FEValues "
6845 "instead"));
6846 Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
6847 this->mapped_geometry->reinit(cell);
6848
6849# ifdef DEBUG
6850 this->is_reinitialized = true;
6851# endif
6852}
6853
6854
6855
6856template <int dim,
6857 int fe_degree,
6858 int n_q_points_1d,
6859 int n_components_,
6860 typename Number,
6861 typename VectorizedArrayType>
6862inline void
6863FEEvaluation<dim,
6864 fe_degree,
6865 n_q_points_1d,
6867 Number,
6868 VectorizedArrayType>::
6869 evaluate(const EvaluationFlags::EvaluationFlags evaluation_flags)
6870{
6871# ifdef DEBUG
6872 Assert(this->dof_values_initialized == true,
6874# endif
6875 evaluate(this->values_dofs, evaluation_flags);
6876}
6877
6878
6879
6880template <int dim,
6881 int fe_degree,
6882 int n_q_points_1d,
6883 int n_components_,
6884 typename Number,
6885 typename VectorizedArrayType>
6886inline void
6887FEEvaluation<dim,
6888 fe_degree,
6889 n_q_points_1d,
6891 Number,
6892 VectorizedArrayType>::
6893 evaluate(const VectorizedArrayType *values_array,
6895{
6896 const bool hessians_on_general_cells =
6898 (this->cell_type > internal::MatrixFreeFunctions::affine);
6902
6903 if (this->data->element_type ==
6905 evaluation_flag & EvaluationFlags::gradients &&
6908
6909 if constexpr (fe_degree > -1)
6910 {
6915 *this);
6916 }
6917 else
6918 {
6920 n_components,
6922 const_cast<VectorizedArrayType *>(values_array),
6923 *this);
6924 }
6925
6926# ifdef DEBUG
6927 this->values_quad_initialized =
6929 this->gradients_quad_initialized =
6931 this->hessians_quad_initialized =
6933# endif
6934}
6935
6936
6937namespace internal
6938{
6942 template <typename Number,
6943 typename VectorizedArrayType,
6944 typename VectorType,
6945 typename EvaluatorType,
6946 std::enable_if_t<internal::has_begin<VectorType> &&
6948 VectorType> * = nullptr>
6949 VectorizedArrayType *
6950 check_vector_access_inplace(const EvaluatorType &fe_eval, VectorType &vector)
6951 {
6952 // for user-defined cell batches this functionality is not supported
6953 if (fe_eval.get_current_cell_index() == numbers::invalid_unsigned_int)
6954 return nullptr;
6955
6956 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
6957 const auto &dof_info = fe_eval.get_dof_info();
6958
6959 // If the index storage is interleaved and contiguous and the vector
6960 // storage has the correct alignment, we can directly pass the pointer
6961 // into the vector to the evaluate() and integrate() calls, without
6962 // reading the vector entries into a separate data field. This saves some
6963 // operations.
6964 if (std::is_same_v<typename VectorType::value_type, Number> &&
6965 dof_info.index_storage_variants
6968 interleaved_contiguous &&
6969 reinterpret_cast<std::size_t>(
6970 vector.begin() +
6971 dof_info.dof_indices_contiguous
6972 [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
6973 [cell * VectorizedArrayType::size()]) %
6974 sizeof(VectorizedArrayType) ==
6975 0)
6976 {
6977 return reinterpret_cast<VectorizedArrayType *>(
6978 vector.begin() +
6979 dof_info.dof_indices_contiguous
6981 [cell * VectorizedArrayType::size()] +
6983 [fe_eval.get_active_fe_index()]
6984 [fe_eval.get_first_selected_component()] *
6985 VectorizedArrayType::size());
6986 }
6987 else
6988 return nullptr;
6989 }
6990
6994 template <typename Number,
6995 typename VectorizedArrayType,
6996 typename VectorType,
6997 typename EvaluatorType,
6998 std::enable_if_t<!internal::has_begin<VectorType> ||
7000 VectorType> * = nullptr>
7001 VectorizedArrayType *
7002 check_vector_access_inplace(const EvaluatorType &, VectorType &)
7003 {
7004 return nullptr;
7005 }
7006} // namespace internal
7007
7008
7009
7010template <int dim,
7011 int fe_degree,
7012 int n_q_points_1d,
7013 int n_components_,
7014 typename Number,
7015 typename VectorizedArrayType>
7016template <typename VectorType>
7017inline void
7018FEEvaluation<dim,
7019 fe_degree,
7020 n_q_points_1d,
7022 Number,
7023 VectorizedArrayType>::
7024 gather_evaluate(const VectorType &input_vector,
7026{
7027 const VectorizedArrayType *src_ptr =
7028 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7029 *this, input_vector);
7030 if (src_ptr != nullptr)
7031 evaluate(src_ptr, evaluation_flag);
7032 else
7033 {
7034 this->read_dof_values(input_vector);
7035 evaluate(this->begin_dof_values(), evaluation_flag);
7036 }
7037}
7038
7039
7040
7041template <int dim,
7042 int fe_degree,
7043 int n_q_points_1d,
7044 int n_components_,
7045 typename Number,
7046 typename VectorizedArrayType>
7047inline void
7048FEEvaluation<dim,
7049 fe_degree,
7050 n_q_points_1d,
7052 Number,
7053 VectorizedArrayType>::
7055{
7056 integrate(integration_flag, this->values_dofs);
7057
7058# ifdef DEBUG
7059 this->dof_values_initialized = true;
7060# endif
7061}
7062
7063
7064
7065template <int dim,
7066 int fe_degree,
7067 int n_q_points_1d,
7068 int n_components_,
7069 typename Number,
7070 typename VectorizedArrayType>
7071inline void
7072FEEvaluation<dim,
7073 fe_degree,
7074 n_q_points_1d,
7076 Number,
7077 VectorizedArrayType>::
7079 VectorizedArrayType *values_array,
7080 const bool sum_into_values_array)
7081{
7082# ifdef DEBUG
7084 Assert(this->values_quad_submitted == true,
7087 Assert(this->gradients_quad_submitted == true,
7090 Assert(this->hessians_quad_submitted == true,
7092# endif
7093 Assert(this->matrix_free != nullptr ||
7094 this->mapped_geometry->is_initialized(),
7096
7097 Assert(
7100 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7101 "EvaluationFlags::hessians are supported."));
7102
7105 (this->cell_type > internal::MatrixFreeFunctions::affine))
7106 {
7107 unsigned int size = n_components * dim * n_q_points;
7109 {
7110 for (unsigned int i = 0; i < size; ++i)
7111 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7112 }
7113 else
7114 {
7115 for (unsigned int i = 0; i < size; ++i)
7116 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7118 }
7119 }
7120
7121 if (n_components == dim &&
7122 this->data->element_type ==
7124 integration_flag & EvaluationFlags::gradients &&
7126 this->divergence_is_requested == false)
7127 {
7128 unsigned int size = n_components * n_q_points;
7130 {
7131 for (unsigned int i = 0; i < size; ++i)
7132 this->values_quad[i] += this->values_from_gradients_quad[i];
7133 }
7134 else
7135 {
7136 for (unsigned int i = 0; i < size; ++i)
7137 this->values_quad[i] = this->values_from_gradients_quad[i];
7139 }
7140 }
7141
7142 if constexpr (fe_degree > -1)
7143 {
7148 *this,
7150 }
7151 else
7152 {
7154 n_components,
7157 *this,
7159 }
7160
7161# ifdef DEBUG
7162 this->dof_values_initialized = true;
7163# endif
7164}
7165
7166
7167
7168template <int dim,
7169 int fe_degree,
7170 int n_q_points_1d,
7171 int n_components_,
7172 typename Number,
7173 typename VectorizedArrayType>
7174template <typename VectorType>
7175inline void
7176FEEvaluation<dim,
7177 fe_degree,
7178 n_q_points_1d,
7180 Number,
7181 VectorizedArrayType>::
7183 VectorType &destination)
7184{
7185 VectorizedArrayType *dst_ptr =
7186 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
7187 *this, destination);
7188 if (dst_ptr != nullptr)
7189 integrate(integration_flag, dst_ptr, true);
7190 else
7191 {
7192 integrate(integration_flag, this->begin_dof_values());
7193 this->distribute_local_to_global(destination);
7194 }
7195}
7196
7197
7198
7199template <int dim,
7200 int fe_degree,
7201 int n_q_points_1d,
7202 int n_components_,
7203 typename Number,
7204 typename VectorizedArrayType>
7206FEEvaluation<dim,
7207 fe_degree,
7208 n_q_points_1d,
7210 Number,
7211 VectorizedArrayType>::dof_indices() const
7212{
7213 return {0U, dofs_per_cell};
7214}
7215
7216
7217
7218/*-------------------------- FEFaceEvaluation ---------------------------*/
7219
7220
7221
7222template <int dim,
7223 int fe_degree,
7224 int n_q_points_1d,
7225 int n_components_,
7226 typename Number,
7227 typename VectorizedArrayType>
7228inline FEFaceEvaluation<dim,
7229 fe_degree,
7230 n_q_points_1d,
7232 Number,
7233 VectorizedArrayType>::
7236 const bool is_interior_face,
7237 const unsigned int dof_no,
7238 const unsigned int quad_no,
7239 const unsigned int first_selected_component,
7240 const unsigned int active_fe_index,
7241 const unsigned int active_quad_index,
7242 const unsigned int face_type)
7243 : BaseClass(matrix_free,
7244 dof_no,
7245 first_selected_component,
7246 quad_no,
7247 fe_degree,
7248 static_n_q_points,
7249 is_interior_face,
7250 active_fe_index,
7251 active_quad_index,
7252 face_type)
7253 , dofs_per_component(this->data->dofs_per_component_on_cell)
7254 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7255 , n_q_points(this->n_quadrature_points)
7256{}
7257
7258
7259
7260template <int dim,
7261 int fe_degree,
7262 int n_q_points_1d,
7263 int n_components_,
7264 typename Number,
7265 typename VectorizedArrayType>
7266inline FEFaceEvaluation<dim,
7267 fe_degree,
7268 n_q_points_1d,
7270 Number,
7271 VectorizedArrayType>::
7274 const std::pair<unsigned int, unsigned int> &range,
7275 const bool is_interior_face,
7276 const unsigned int dof_no,
7277 const unsigned int quad_no,
7278 const unsigned int first_selected_component)
7279 : FEFaceEvaluation(matrix_free,
7280 is_interior_face,
7281 dof_no,
7282 quad_no,
7283 first_selected_component,
7284 matrix_free.get_face_active_fe_index(range,
7285 is_interior_face),
7287 matrix_free.get_face_info(range.first).face_type)
7288{}
7289
7290
7291
7292template <int dim,
7293 int fe_degree,
7294 int n_q_points_1d,
7295 int n_components_,
7296 typename Number,
7297 typename VectorizedArrayType>
7298inline void
7300 fe_degree,
7301 n_q_points_1d,
7303 Number,
7304 VectorizedArrayType>::reinit(const unsigned int face_index)
7305{
7306 Assert(this->mapped_geometry == nullptr,
7307 ExcMessage("FEEvaluation was initialized without a matrix-free object."
7308 " Integer indexing is not possible"));
7309 if (this->mapped_geometry != nullptr)
7310 return;
7311
7312 this->cell = face_index;
7313 this->dof_access_index =
7314 this->is_interior_face() ?
7317 Assert(this->mapping_data != nullptr, ExcNotInitialized());
7318
7319 if (face_index >=
7320 this->matrix_free->get_task_info().face_partition_data.back() &&
7321 face_index <
7322 this->matrix_free->get_task_info().boundary_partition_data.back())
7323 Assert(this->is_interior_face(),
7324 ExcMessage(
7325 "Boundary faces do not have a neighbor. When looping over "
7326 "boundary faces use FEFaceEvaluation with the parameter "
7327 "is_interior_face set to true. "));
7328
7329 this->reinit_face(this->matrix_free->get_face_info(face_index));
7330
7331 unsigned int i = 0;
7332 for (; i < this->matrix_free->n_active_entries_per_face_batch(this->cell);
7333 ++i)
7334 this->face_ids[i] = face_index * n_lanes + i;
7335 for (; i < n_lanes; ++i)
7336 this->face_ids[i] = numbers::invalid_unsigned_int;
7337
7338 this->cell_type = this->matrix_free->get_mapping_info().face_type[face_index];
7339 const unsigned int offsets =
7340 this->mapping_data->data_index_offsets[face_index];
7341 this->J_value = &this->mapping_data->JxW_values[offsets];
7342 this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
7343 this->jacobian =
7344 &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
7345 this->normal_x_jacobian =
7346 &this->mapping_data
7347 ->normals_times_jacobians[!this->is_interior_face()][offsets];
7348 this->jacobian_gradients =
7349 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7350 offsets;
7351 this->jacobian_gradients_non_inverse =
7352 this->mapping_data
7353 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7354 .data() +
7355 offsets;
7356
7357 if (this->mapping_data->quadrature_point_offsets.empty() == false)
7358 {
7359 AssertIndexRange(this->cell,
7360 this->mapping_data->quadrature_point_offsets.size());
7361 this->quadrature_points =
7362 this->mapping_data->quadrature_points.data() +
7363 this->mapping_data->quadrature_point_offsets[this->cell];
7364 }
7365
7366# ifdef DEBUG
7367 this->is_reinitialized = true;
7368 this->dof_values_initialized = false;
7369 this->values_quad_initialized = false;
7370 this->gradients_quad_initialized = false;
7371 this->hessians_quad_initialized = false;
7372# endif
7373}
7374
7375
7376
7377template <int dim,
7378 int fe_degree,
7379 int n_q_points_1d,
7380 int n_components_,
7381 typename Number,
7382 typename VectorizedArrayType>
7383inline void
7385 fe_degree,
7386 n_q_points_1d,
7388 Number,
7389 VectorizedArrayType>::reinit(const unsigned int cell_index,
7390 const unsigned int face_number)
7391{
7392 Assert(
7393 this->quad_no <
7394 this->matrix_free->get_mapping_info().face_data_by_cells.size(),
7395 ExcMessage(
7396 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
7399 this->matrix_free->get_mapping_info().cell_type.size());
7400 Assert(this->mapped_geometry == nullptr,
7401 ExcMessage("FEEvaluation was initialized without a matrix-free object."
7402 " Integer indexing is not possible"));
7403 if (this->mapped_geometry != nullptr)
7404 return;
7405 Assert(this->matrix_free != nullptr, ExcNotInitialized());
7406
7407 this->cell_type = this->matrix_free->get_mapping_info()
7408 .faces_by_cells_type[cell_index][face_number];
7409 this->cell = cell_index;
7410 this->subface_index = GeometryInfo<dim>::max_children_per_cell;
7411 this->dof_access_index =
7413
7414 if (this->is_interior_face() == false)
7415 {
7416 // for this case, we need to look into the FaceInfo field that collects
7417 // information from both sides of a face once for the global mesh, and
7418 // pick the face id that is not the local one (cell_this).
7419 for (unsigned int i = 0; i < n_lanes; ++i)
7420 {
7421 // compute actual (non vectorized) cell ID
7422 const unsigned int cell_this = cell_index * n_lanes + i;
7423 // compute face ID
7424 unsigned int face_index =
7425 this->matrix_free->get_cell_and_face_to_plain_faces()(cell_index,
7426 face_number,
7427 i);
7428
7429 this->face_ids[i] = face_index;
7430
7431 if (face_index == numbers::invalid_unsigned_int)
7432 {
7433 this->cell_ids[i] = numbers::invalid_unsigned_int;
7434 this->face_numbers[i] = static_cast<std::uint8_t>(-1);
7435 this->face_orientations[i] = static_cast<std::uint8_t>(-1);
7436 continue; // invalid face ID: no neighbor on boundary
7437 }
7438
7439 const auto &faces =
7440 this->matrix_free->get_face_info(face_index / n_lanes);
7441 // get cell ID on both sides of face
7442 auto cell_m = faces.cells_interior[face_index % n_lanes];
7443 auto cell_p = faces.cells_exterior[face_index % n_lanes];
7444
7446
7449
7450 // compare the IDs with the given cell ID
7452 {
7453 this->cell_ids[i] = cell_m; // neighbor has the other ID
7454 this->face_numbers[i] = faces.interior_face_no;
7455 }
7456 else
7457 {
7458 this->cell_ids[i] = cell_p;
7459 this->face_numbers[i] = faces.exterior_face_no;
7460 }
7461
7462 const bool orientation_interior_face = faces.face_orientation >= 8;
7463 unsigned int face_orientation = faces.face_orientation % 8;
7465 {
7466 constexpr std::array<std::uint8_t, 8> table{
7467 {0, 1, 6, 3, 4, 5, 2, 7}};
7468 face_orientation = table[face_orientation];
7469 }
7470 this->face_orientations[i] = face_orientation;
7471 }
7472 }
7473 else
7474 {
7475 this->face_orientations[0] = 0;
7476 this->face_numbers[0] = face_number;
7477 if (this->matrix_free->n_active_entries_per_cell_batch(this->cell) ==
7478 n_lanes)
7479 {
7481 for (unsigned int i = 0; i < n_lanes; ++i)
7482 this->cell_ids[i] = cell_index * n_lanes + i;
7483 }
7484 else
7485 {
7486 unsigned int i = 0;
7487 for (; i <
7488 this->matrix_free->n_active_entries_per_cell_batch(this->cell);
7489 ++i)
7490 this->cell_ids[i] = cell_index * n_lanes + i;
7491 for (; i < n_lanes; ++i)
7492 this->cell_ids[i] = numbers::invalid_unsigned_int;
7493 }
7494 for (unsigned int i = 0; i < n_lanes; ++i)
7495 this->face_ids[i] =
7496 this->matrix_free->get_cell_and_face_to_plain_faces()(cell_index,
7497 face_number,
7498 i);
7499 }
7500
7501 const unsigned int offsets =
7502 this->matrix_free->get_mapping_info()
7503 .face_data_by_cells[this->quad_no]
7504 .data_index_offsets[cell_index * GeometryInfo<dim>::faces_per_cell +
7505 face_number];
7507 this->matrix_free->get_mapping_info()
7508 .face_data_by_cells[this->quad_no]
7509 .JxW_values.size());
7510 this->J_value = &this->matrix_free->get_mapping_info()
7511 .face_data_by_cells[this->quad_no]
7512 .JxW_values[offsets];
7513 this->normal_vectors = &this->matrix_free->get_mapping_info()
7514 .face_data_by_cells[this->quad_no]
7515 .normal_vectors[offsets];
7516 this->jacobian = &this->matrix_free->get_mapping_info()
7517 .face_data_by_cells[this->quad_no]
7518 .jacobians[!this->is_interior_face()][offsets];
7519 this->normal_x_jacobian =
7520 &this->matrix_free->get_mapping_info()
7521 .face_data_by_cells[this->quad_no]
7522 .normals_times_jacobians[!this->is_interior_face()][offsets];
7523 this->jacobian_gradients =
7524 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7525 offsets;
7526 this->jacobian_gradients_non_inverse =
7527 this->mapping_data
7528 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7529 .data() +
7530 offsets;
7531
7532 if (this->matrix_free->get_mapping_info()
7533 .face_data_by_cells[this->quad_no]
7534 .quadrature_point_offsets.empty() == false)
7535 {
7536 const unsigned int index =
7537 this->cell * GeometryInfo<dim>::faces_per_cell + this->face_numbers[0];
7539 this->matrix_free->get_mapping_info()
7540 .face_data_by_cells[this->quad_no]
7541 .quadrature_point_offsets.size());
7542 this->quadrature_points = this->matrix_free->get_mapping_info()
7543 .face_data_by_cells[this->quad_no]
7544 .quadrature_points.data() +
7545 this->matrix_free->get_mapping_info()
7546 .face_data_by_cells[this->quad_no]
7547 .quadrature_point_offsets[index];
7548 }
7549
7550# ifdef DEBUG
7551 this->is_reinitialized = true;
7552 this->dof_values_initialized = false;
7553 this->values_quad_initialized = false;
7554 this->gradients_quad_initialized = false;
7555 this->hessians_quad_initialized = false;
7556# endif
7557}
7558
7559
7560
7561template <int dim,
7562 int fe_degree,
7563 int n_q_points_1d,
7564 int n_components_,
7565 typename Number,
7566 typename VectorizedArrayType>
7567inline void
7569 fe_degree,
7570 n_q_points_1d,
7572 Number,
7573 VectorizedArrayType>::
7575{
7576# ifdef DEBUG
7577 Assert(this->dof_values_initialized, ExcNotInitialized());
7578# endif
7579
7580 evaluate(this->values_dofs, evaluation_flag);
7581}
7582
7583
7584
7585template <int dim,
7586 int fe_degree,
7587 int n_q_points_1d,
7588 int n_components_,
7589 typename Number,
7590 typename VectorizedArrayType>
7591inline void
7593 fe_degree,
7594 n_q_points_1d,
7596 Number,
7597 VectorizedArrayType>::
7598 evaluate(const VectorizedArrayType *values_array,
7600{
7604 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
7605 "and EvaluationFlags::hessians are supported."));
7606
7607 const bool hessians_on_general_cells =
7609 (this->cell_type > internal::MatrixFreeFunctions::affine);
7613
7614 if (this->data->element_type ==
7616 evaluation_flag & EvaluationFlags::gradients &&
7619
7620 if constexpr (fe_degree > -1)
7622 template run<fe_degree, n_q_points_1d>(n_components,
7625 *this);
7626 else
7628 n_components, evaluation_flag_actual, values_array, *this);
7629
7630# ifdef DEBUG
7631 this->values_quad_initialized =
7633 this->gradients_quad_initialized =
7635 this->hessians_quad_initialized =
7637# endif
7638}
7639
7640
7641
7642template <int dim,
7643 int fe_degree,
7644 int n_q_points_1d,
7645 int n_components_,
7646 typename Number,
7647 typename VectorizedArrayType>
7648inline void
7650 fe_degree,
7651 n_q_points_1d,
7653 Number,
7654 VectorizedArrayType>::
7656{
7657# ifdef DEBUG
7658 Assert(this->dof_values_initialized, ExcNotInitialized());
7659# endif
7660
7661 project_to_face(this->values_dofs, evaluation_flag);
7662}
7663
7664
7665
7666template <int dim,
7667 int fe_degree,
7668 int n_q_points_1d,
7669 int n_components_,
7670 typename Number,
7671 typename VectorizedArrayType>
7672inline void
7674 fe_degree,
7675 n_q_points_1d,
7677 Number,
7678 VectorizedArrayType>::
7679 project_to_face(const VectorizedArrayType *values_array,
7681{
7685 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
7686 "and EvaluationFlags::hessians are supported."));
7687
7688 const bool hessians_on_general_cells =
7690 (this->cell_type > internal::MatrixFreeFunctions::affine);
7694
7695 if (this->data->element_type ==
7697 evaluation_flag & EvaluationFlags::gradients &&
7700
7701 if constexpr (fe_degree > -1)
7703 dim,
7704 VectorizedArrayType>::template run<fe_degree>(n_components,
7707 *this);
7708 else
7710 project_to_face(n_components,
7713 *this);
7714
7715 // face dofs initialized
7716}
7717
7718
7719
7720template <int dim,
7721 int fe_degree,
7722 int n_q_points_1d,
7723 int n_components_,
7724 typename Number,
7725 typename VectorizedArrayType>
7726inline void
7728 fe_degree,
7729 n_q_points_1d,
7731 Number,
7732 VectorizedArrayType>::
7734{
7738 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
7739 "and EvaluationFlags::hessians are supported."));
7740
7741 const bool hessians_on_general_cells =
7743 (this->cell_type > internal::MatrixFreeFunctions::affine);
7747
7748 if (this->data->element_type ==
7750 evaluation_flag & EvaluationFlags::gradients &&
7753
7754 if constexpr (fe_degree > -1)
7756 dim,
7757 VectorizedArrayType>::template run<fe_degree>(n_components,
7759 *this);
7760 else
7762 evaluate_in_face(n_components, evaluation_flag_actual, *this);
7763
7764# ifdef DEBUG
7765 this->values_quad_initialized =
7767 this->gradients_quad_initialized =
7769 this->hessians_quad_initialized =
7771# endif
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,
7787 Number,
7788 VectorizedArrayType>::
7790 const bool sum_into_values)
7791{
7792 integrate(integration_flag, this->values_dofs, sum_into_values);
7793
7794# ifdef DEBUG
7795 this->dof_values_initialized = true;
7796# endif
7797}
7798
7799
7800
7801template <int dim,
7802 int fe_degree,
7803 int n_q_points_1d,
7804 int n_components_,
7805 typename Number,
7806 typename VectorizedArrayType>
7807inline void
7809 fe_degree,
7810 n_q_points_1d,
7812 Number,
7813 VectorizedArrayType>::
7815 VectorizedArrayType *values_array,
7816 const bool sum_into_values)
7817{
7821 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
7822 "and EvaluationFlags::hessians are supported."));
7823
7826 (this->cell_type > internal::MatrixFreeFunctions::affine))
7827 {
7828 unsigned int size = n_components * dim * n_q_points;
7830 {
7831 for (unsigned int i = 0; i < size; ++i)
7832 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7833 }
7834 else
7835 {
7836 for (unsigned int i = 0; i < size; ++i)
7837 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7839 }
7840 }
7841
7842 if (this->data->element_type ==
7844 integration_flag & EvaluationFlags::gradients &&
7846 this->divergence_is_requested == false)
7847 {
7848 unsigned int size = n_components * n_q_points;
7850 {
7851 for (unsigned int i = 0; i < size; ++i)
7852 this->values_quad[i] += this->values_from_gradients_quad[i];
7853 }
7854 else
7855 {
7856 for (unsigned int i = 0; i < size; ++i)
7857 this->values_quad[i] = this->values_from_gradients_quad[i];
7859 }
7860 }
7861
7862 if constexpr (fe_degree > -1)
7864 template run<fe_degree, n_q_points_1d>(n_components,
7867 *this,
7869 else
7871 n_components,
7874 *this,
7876}
7877
7878
7879
7880template <int dim,
7881 int fe_degree,
7882 int n_q_points_1d,
7883 int n_components_,
7884 typename Number,
7885 typename VectorizedArrayType>
7886inline void
7888 fe_degree,
7889 n_q_points_1d,
7891 Number,
7892 VectorizedArrayType>::
7894{
7898 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
7899 "and EvaluationFlags::hessians are supported."));
7900
7903 (this->cell_type > internal::MatrixFreeFunctions::affine))
7904 {
7905 unsigned int size = n_components * dim * n_q_points;
7907 {
7908 for (unsigned int i = 0; i < size; ++i)
7909 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7910 }
7911 else
7912 {
7913 for (unsigned int i = 0; i < size; ++i)
7914 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7916 }
7917 }
7918
7919 if (this->data->element_type ==
7921 integration_flag & EvaluationFlags::gradients &&
7923 this->divergence_is_requested == false)
7924 {
7925 unsigned int size = n_components * n_q_points;
7927 {
7928 for (unsigned int i = 0; i < size; ++i)
7929 this->values_quad[i] += this->values_from_gradients_quad[i];
7930 }
7931 else
7932 {
7933 for (unsigned int i = 0; i < size; ++i)
7934 this->values_quad[i] = this->values_from_gradients_quad[i];
7936 }
7937 }
7938
7939 if constexpr (fe_degree > -1)
7941 dim,
7942 VectorizedArrayType>::template run<fe_degree>(n_components,
7944 *this);
7945 else
7947 integrate_in_face(n_components, integration_flag_actual, *this);
7948
7949 // face dofs initialized
7950}
7951
7952
7953
7954template <int dim,
7955 int fe_degree,
7956 int n_q_points_1d,
7957 int n_components_,
7958 typename Number,
7959 typename VectorizedArrayType>
7960inline void
7962 fe_degree,
7963 n_q_points_1d,
7965 Number,
7966 VectorizedArrayType>::
7968 const bool sum_into_values)
7969{
7970 collect_from_face(integration_flag, this->values_dofs, sum_into_values);
7971
7972# ifdef DEBUG
7973 this->dof_values_initialized = true;
7974# endif
7975}
7976
7977
7978
7979template <int dim,
7980 int fe_degree,
7981 int n_q_points_1d,
7982 int n_components_,
7983 typename Number,
7984 typename VectorizedArrayType>
7985inline void
7987 fe_degree,
7988 n_q_points_1d,
7990 Number,
7991 VectorizedArrayType>::
7993 VectorizedArrayType *values_array,
7994 const bool sum_into_values)
7995{
7999 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
8000 "and EvaluationFlags::hessians are supported."));
8001
8004 (this->cell_type > internal::MatrixFreeFunctions::affine))
8006
8007 if (this->data->element_type ==
8009 integration_flag & EvaluationFlags::gradients &&
8011 this->divergence_is_requested == false)
8013
8014 if constexpr (fe_degree > -1)
8016 dim,
8017 VectorizedArrayType>::template run<fe_degree>(n_components,
8020 *this,
8022 else
8024 collect_from_face(n_components,
8027 *this,
8029}
8030
8031
8032
8033template <int dim,
8034 int fe_degree,
8035 int n_q_points_1d,
8036 int n_components_,
8037 typename Number,
8038 typename VectorizedArrayType>
8039template <typename VectorType>
8040inline void
8042 fe_degree,
8043 n_q_points_1d,
8045 Number,
8046 VectorizedArrayType>::
8047 gather_evaluate(const VectorType &input_vector,
8049{
8053 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
8054 "and EvaluationFlags::hessians are supported."));
8055
8056 const auto shared_vector_data = internal::get_shared_vector_data(
8057 &input_vector,
8058 this->dof_access_index ==
8060 this->active_fe_index,
8061 this->dof_info);
8062
8063 if (this->data->data.front().fe_degree > 0 &&
8064 fast_evaluation_supported(this->data->data.front().fe_degree,
8065 this->data->data.front().n_q_points_1d) &&
8067 dim,
8068 typename VectorType::value_type,
8069 VectorizedArrayType>::
8070 supports(evaluation_flag,
8071 *this->data,
8073 input_vector),
8074 this->dof_info->index_storage_variants[this->dof_access_index]
8075 [this->cell]))
8076 {
8077 if constexpr (fe_degree > -1)
8078 {
8080 dim,
8081 typename VectorType::value_type,
8082 VectorizedArrayType>::template run<fe_degree,
8083 n_q_points_1d>(
8084 n_components,
8086 internal::get_beginning<typename VectorType::value_type>(
8087 input_vector),
8088 shared_vector_data,
8089 *this);
8090 }
8091 else
8092 {
8094 dim,
8095 typename VectorType::value_type,
8096 VectorizedArrayType>::evaluate(n_components,
8098 internal::get_beginning<
8099 typename VectorType::value_type>(
8100 input_vector),
8101 shared_vector_data,
8102 *this);
8103 }
8104 }
8105 else
8106 {
8107 this->read_dof_values(input_vector);
8108 this->evaluate(evaluation_flag);
8109 }
8110
8111# ifdef DEBUG
8112 this->values_quad_initialized = evaluation_flag & EvaluationFlags::values;
8113 this->gradients_quad_initialized =
8115 this->hessians_quad_initialized = evaluation_flag & EvaluationFlags::hessians;
8116# endif
8117}
8118
8119
8120
8121template <int dim,
8122 int fe_degree,
8123 int n_q_points_1d,
8124 int n_components_,
8125 typename Number,
8126 typename VectorizedArrayType>
8127template <typename VectorType>
8128inline void
8130 dim,
8131 fe_degree,
8132 n_q_points_1d,
8134 Number,
8135 VectorizedArrayType>::integrate_scatter(const bool integrate_values,
8136 const bool integrate_gradients,
8137 VectorType &destination)
8138{
8140 ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8141 ((integrate_gradients) ? EvaluationFlags::gradients :
8143
8144 integrate_scatter(flag, destination);
8145}
8146
8147
8148
8149template <int dim,
8150 int fe_degree,
8151 int n_q_points_1d,
8152 int n_components_,
8153 typename Number,
8154 typename VectorizedArrayType>
8155template <typename VectorType>
8156inline void
8158 fe_degree,
8159 n_q_points_1d,
8161 Number,
8162 VectorizedArrayType>::
8164 VectorType &destination)
8165{
8166 Assert((this->dof_access_index ==
8168 this->is_interior_face() == false) == false,
8170
8171 const auto shared_vector_data = internal::get_shared_vector_data(
8172 &destination,
8173 this->dof_access_index ==
8175 this->active_fe_index,
8176 this->dof_info);
8177
8178 if (this->data->data.front().fe_degree > 0 &&
8179 fast_evaluation_supported(this->data->data.front().fe_degree,
8180 this->data->data.front().n_q_points_1d) &&
8182 dim,
8183 typename VectorType::value_type,
8184 VectorizedArrayType>::
8185 supports(integration_flag,
8186 *this->data,
8188 destination),
8189 this->dof_info->index_storage_variants[this->dof_access_index]
8190 [this->cell]))
8191 {
8192 if constexpr (fe_degree > -1)
8193 {
8195 dim,
8196 typename VectorType::value_type,
8197 VectorizedArrayType>::template run<fe_degree,
8198 n_q_points_1d>(
8199 n_components,
8201 internal::get_beginning<typename VectorType::value_type>(
8202 destination),
8203 shared_vector_data,
8204 *this);
8205 }
8206 else
8207 {
8209 dim,
8210 typename VectorType::value_type,
8211 VectorizedArrayType>::integrate(n_components,
8213 internal::get_beginning<
8214 typename VectorType::value_type>(
8215 destination),
8216 shared_vector_data,
8217 *this);
8218 }
8219 }
8220 else
8221 {
8222 integrate(integration_flag);
8223 this->distribute_local_to_global(destination);
8224 }
8225}
8226
8227
8228
8229template <int dim,
8230 int fe_degree,
8231 int n_q_points_1d,
8232 int n_components_,
8233 typename Number,
8234 typename VectorizedArrayType>
8237 fe_degree,
8238 n_q_points_1d,
8240 Number,
8241 VectorizedArrayType>::dof_indices() const
8242{
8243 return {0U, dofs_per_cell};
8244}
8245
8246
8247
8248template <int dim,
8249 int fe_degree,
8250 int n_q_points_1d,
8251 int n_components_,
8252 typename Number,
8253 typename VectorizedArrayType>
8254bool
8255FEEvaluation<dim,
8256 fe_degree,
8257 n_q_points_1d,
8259 Number,
8260 VectorizedArrayType>::
8261 fast_evaluation_supported(const unsigned int given_degree,
8262 const unsigned int given_n_q_points_1d)
8263{
8264 return fe_degree == -1 ?
8267 true;
8268}
8269
8270
8271
8272template <int dim,
8273 int fe_degree,
8274 int n_q_points_1d,
8275 int n_components_,
8276 typename Number,
8277 typename VectorizedArrayType>
8278bool
8280 fe_degree,
8281 n_q_points_1d,
8283 Number,
8284 VectorizedArrayType>::
8285 fast_evaluation_supported(const unsigned int given_degree,
8286 const unsigned int given_n_q_points_1d)
8287{
8288 return fe_degree == -1 ?
8291 true;
8292}
8293
8294
8295
8296template <int dim,
8297 int fe_degree,
8298 int n_q_points_1d,
8299 int n_components_,
8300 typename Number,
8301 typename VectorizedArrayType>
8302bool
8304 fe_degree,
8305 n_q_points_1d,
8307 Number,
8308 VectorizedArrayType>::at_boundary() const
8309{
8310 Assert(this->dof_access_index !=
8313
8314 if (this->is_interior_face() == false)
8315 return false;
8316 else if (this->cell < this->matrix_free->n_inner_face_batches())
8317 return false;
8318 else if (this->cell < (this->matrix_free->n_inner_face_batches() +
8319 this->matrix_free->n_boundary_face_batches()))
8320 return true;
8321 else
8322 return false;
8323}
8324
8325
8326
8327template <int dim,
8328 int fe_degree,
8329 int n_q_points_1d,
8330 int n_components_,
8331 typename Number,
8332 typename VectorizedArrayType>
8335 fe_degree,
8336 n_q_points_1d,
8338 Number,
8339 VectorizedArrayType>::boundary_id() const
8340{
8341 Assert(this->dof_access_index !=
8344
8345 if (at_boundary())
8346 return this->matrix_free->get_boundary_id(this->cell);
8347 else
8349}
8350
8351
8352
8353template <int dim,
8354 int fe_degree,
8355 int n_q_points_1d,
8356 int n_components_,
8357 typename Number,
8358 typename VectorizedArrayType>
8359unsigned int
8361 dim,
8362 fe_degree,
8363 n_q_points_1d,
8365 Number,
8366 VectorizedArrayType>::get_dofs_per_component_projected_to_face()
8367{
8368 return this->data->dofs_per_component_on_face;
8369}
8370
8371
8372
8373template <int dim,
8374 int fe_degree,
8375 int n_q_points_1d,
8376 int n_components_,
8377 typename Number,
8378 typename VectorizedArrayType>
8379unsigned int
8381 fe_degree,
8382 n_q_points_1d,
8384 Number,
8385 VectorizedArrayType>::get_dofs_projected_to_face()
8386{
8387 return this->data->dofs_per_component_on_face * n_components_;
8388}
8389
8390
8391
8392/*------------------------- end FEFaceEvaluation ------------------------- */
8393
8394
8395#endif // ifndef DOXYGEN
8396
8397
8399
8400#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())
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::conditional_t< rank_==1, Number, Tensor< rank_ - 1, dim, Number > > value_type
Definition tensor.h:506
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:156
#define DEAL_II_DEPRECATED
Definition config.h:228
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Point< 2 > second
Definition grid_out.cc:4630
Point< 2 > first
Definition grid_out.cc:4629
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.
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< index_type > data
Definition mpi.cc:740
std::size_t size
Definition mpi.cc:739
The namespace for the EvaluationFlags enum.
EvaluationFlags
The EvaluationFlags enum.
constexpr char U
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:474
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
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::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
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)