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