deal.II version GIT relicensing-2220-ge8e5360a96 2024-12-05 09: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
156 template <typename VectorType>
157 void
159 const VectorType &src,
160 const unsigned int first_index = 0,
161 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
162
191 template <typename VectorType>
192 void
194 const VectorType &src,
195 const unsigned int first_index = 0,
196 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
197
229 template <typename VectorType>
230 void
232 VectorType &dst,
233 const unsigned int first_index = 0,
234 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
235
274 template <typename VectorType>
275 void
277 VectorType &dst,
278 const unsigned int first_index = 0,
279 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
280
284 template <typename VectorType>
285 void
287 VectorType &dst,
288 const unsigned int first_index = 0,
289 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
290
313 get_dof_value(const unsigned int dof) const;
314
323 void
324 submit_dof_value(const value_type val_in, const unsigned int dof);
325
337 get_value(const unsigned int q_point) const;
338
354 void
355 submit_value(const value_type val_in, const unsigned int q_point);
356
364 template <int n_components_local = n_components,
365 typename = std::enable_if_t<n_components == n_components_local>>
366 void
368 const unsigned int q_point);
369
380 get_gradient(const unsigned int q_point) const;
381
393 get_normal_derivative(const unsigned int q_point) const;
394
409 void
410 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
411
419 template <int dim_ = dim,
420 typename = std::enable_if_t<dim_ == 1 && n_components == dim_>>
421 void
423 const unsigned int q_point);
424
440 void
442 const unsigned int q_point);
443
452 get_hessian(const unsigned int q_point) const;
453
460 get_hessian_diagonal(const unsigned int q_point) const;
461
470 get_laplacian(const unsigned int q_point) const;
471
482 get_normal_hessian(const unsigned int q_point) const;
483
498 void
499 submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
500
516 void
517 submit_normal_hessian(const value_type normal_hessian_in,
518 const unsigned int q_point);
519
527 template <int dim_ = dim, typename = std::enable_if_t<n_components_ == dim_>>
528 VectorizedArrayType
529 get_divergence(const unsigned int q_point) const;
530
546 template <int dim_ = dim, typename = std::enable_if_t<n_components_ == dim_>>
547 void
548 submit_divergence(const VectorizedArrayType div_in,
549 const unsigned int q_point);
550
559 template <int dim_ = dim, typename = std::enable_if_t<n_components_ == dim_>>
561 get_symmetric_gradient(const unsigned int q_point) const;
562
578 template <int dim_ = dim, typename = std::enable_if_t<n_components_ == dim_>>
579 void
582 const unsigned int q_point);
583
592 template <int dim_ = dim,
593 typename = std::enable_if_t<n_components_ == dim_ && dim_ != 1>>
594 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
595 get_curl(const unsigned int q_point) const;
596
610 template <int dim_ = dim,
611 typename = std::enable_if_t<n_components_ == dim_ && dim != 1>>
612 void
614 const unsigned int q_point);
615
635
643
644protected:
655 const unsigned int dof_no,
656 const unsigned int first_selected_component,
657 const unsigned int quad_no,
658 const unsigned int fe_degree,
659 const unsigned int n_q_points,
660 const bool is_interior_face,
661 const unsigned int active_fe_index,
662 const unsigned int active_quad_index,
663 const unsigned int face_type);
664
702 const Mapping<dim> &mapping,
703 const FiniteElement<dim> &fe,
704 const Quadrature<1> &quadrature,
705 const UpdateFlags update_flags,
706 const unsigned int first_selected_component,
708
716
725
730
737 template <typename VectorType, typename VectorOperation>
738 void
740 const VectorOperation &operation,
741 const std::array<VectorType *, n_components_> &vectors,
742 const std::array<
744 n_components_> &vectors_sm,
745 const std::bitset<n_lanes> &mask,
746 const bool apply_constraints = true) const;
747
755 template <typename VectorType, typename VectorOperation>
756 void
758 const VectorOperation &operation,
759 const std::array<VectorType *, n_components_> &vectors,
760 const std::array<
762 n_components_> &vectors_sm,
763 const std::bitset<n_lanes> &mask) const;
764
772 template <typename VectorType, typename VectorOperation>
773 void
775 const VectorOperation &operation,
776 const std::array<VectorType *, n_components_> &vectors) const;
777
781 void
783
788
793
798 mutable std::vector<types::global_dof_index> local_dof_indices;
799};
800
801
802
803// backward compatibility
804template <int dim,
805 int n_components_,
806 typename Number,
807 bool is_face,
808 typename VectorizedArrayType = VectorizedArray<Number>>
811
1375template <int dim,
1376 int fe_degree,
1377 int n_q_points_1d,
1378 int n_components_,
1379 typename Number,
1380 typename VectorizedArrayType>
1382 n_components_,
1383 Number,
1384 false,
1385 VectorizedArrayType>
1386{
1387 static_assert(
1388 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1389 "Type of Number and of VectorizedArrayType do not match.");
1390
1391public:
1397
1401 using number_type = Number;
1402
1409
1416
1420 static constexpr unsigned int dimension = dim;
1421
1426 static constexpr unsigned int n_components = n_components_;
1427
1431 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1432
1441 static constexpr unsigned int static_n_q_points =
1442 Utilities::pow(n_q_points_1d, dim);
1443
1453 static constexpr unsigned int static_dofs_per_component =
1454 Utilities::pow(fe_degree + 1, dim);
1455
1465 static constexpr unsigned int tensor_dofs_per_cell =
1467
1477 static constexpr unsigned int static_dofs_per_cell =
1479
1516 const unsigned int dof_no = 0,
1517 const unsigned int quad_no = 0,
1518 const unsigned int first_selected_component = 0,
1521
1530 const std::pair<unsigned int, unsigned int> &range,
1531 const unsigned int dof_no = 0,
1532 const unsigned int quad_no = 0,
1533 const unsigned int first_selected_component = 0);
1534
1564 const FiniteElement<dim> &fe,
1565 const Quadrature<1> &quadrature,
1566 const UpdateFlags update_flags,
1567 const unsigned int first_selected_component = 0);
1568
1575 const Quadrature<1> &quadrature,
1576 const UpdateFlags update_flags,
1577 const unsigned int first_selected_component = 0);
1578
1591 const unsigned int first_selected_component = 0);
1592
1600
1607 FEEvaluation &
1609
1618 void
1619 reinit(const unsigned int cell_batch_index);
1620
1627 void
1628 reinit(const std::array<unsigned int, n_lanes> &cell_ids);
1629
1642 template <bool level_dof_access>
1643 void
1645
1656 void
1658
1662 static bool
1663 fast_evaluation_supported(const unsigned int given_degree,
1664 const unsigned int given_n_q_points_1d);
1665
1675 void
1677
1690 void
1691 evaluate(const VectorizedArrayType *values_array,
1692 const EvaluationFlags::EvaluationFlags evaluation_flag);
1693
1707 template <typename VectorType>
1708 void
1709 gather_evaluate(const VectorType &input_vector,
1710 const EvaluationFlags::EvaluationFlags evaluation_flag);
1711
1721 void
1723
1735 void
1737 VectorizedArrayType *values_array,
1738 const bool sum_into_values = false);
1739
1753 template <typename VectorType>
1754 void
1756 VectorType &output_vector);
1757
1765
1772 const unsigned int dofs_per_component;
1773
1780 const unsigned int dofs_per_cell;
1781
1789 const unsigned int n_q_points;
1790
1791private:
1796 void
1797 check_template_arguments(const unsigned int fe_no,
1798 const unsigned int first_selected_component);
1799};
1800
1801
1802
1838template <int dim,
1839 int fe_degree,
1840 int n_q_points_1d = fe_degree + 1,
1841 int n_components_ = 1,
1842 typename Number = double,
1843 typename VectorizedArrayType = VectorizedArray<Number>>
1845 n_components_,
1846 Number,
1847 true,
1848 VectorizedArrayType>
1849{
1850 static_assert(
1851 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1852 "Type of Number and of VectorizedArrayType do not match.");
1853
1854public:
1860
1864 using number_type = Number;
1865
1872
1879
1883 static constexpr unsigned int dimension = dim;
1884
1889 static constexpr unsigned int n_components = n_components_;
1890
1894 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1895
1905 static constexpr unsigned int static_n_q_points =
1906 Utilities::pow(n_q_points_1d, dim - 1);
1907
1916 static constexpr unsigned int static_n_q_points_cell =
1917 Utilities::pow(n_q_points_1d, dim);
1918
1927 static constexpr unsigned int static_dofs_per_component =
1928 Utilities::pow(fe_degree + 1, dim);
1929
1938 static constexpr unsigned int tensor_dofs_per_cell =
1940
1949 static constexpr unsigned int static_dofs_per_cell =
1951
1995 const bool is_interior_face = true,
1996 const unsigned int dof_no = 0,
1997 const unsigned int quad_no = 0,
1998 const unsigned int first_selected_component = 0,
2001 const unsigned int face_type = numbers::invalid_unsigned_int);
2002
2012 const std::pair<unsigned int, unsigned int> &range,
2013 const bool is_interior_face = true,
2014 const unsigned int dof_no = 0,
2015 const unsigned int quad_no = 0,
2016 const unsigned int first_selected_component = 0);
2017
2028 void
2029 reinit(const unsigned int face_batch_number);
2030
2038 void
2039 reinit(const unsigned int cell_batch_number, const unsigned int face_number);
2040
2044 static bool
2045 fast_evaluation_supported(const unsigned int given_degree,
2046 const unsigned int given_n_q_points_1d);
2047
2058 void
2060
2073 void
2074 evaluate(const VectorizedArrayType *values_array,
2075 const EvaluationFlags::EvaluationFlags evaluation_flag);
2076
2081 void
2083
2088 void
2089 project_to_face(const VectorizedArrayType *values_array,
2090 const EvaluationFlags::EvaluationFlags evaluation_flag);
2091
2096 void
2098
2110 template <typename VectorType>
2111 void
2112 gather_evaluate(const VectorType &input_vector,
2113 const EvaluationFlags::EvaluationFlags evaluation_flag);
2114
2124 void
2126 const bool sum_into_values = false);
2127
2137 void
2139 VectorizedArrayType *values_array,
2140 const bool sum_into_values = false);
2141
2148 void
2150
2155 void
2157 const bool sum_into_values = false);
2158
2163 void
2165 VectorizedArrayType *values_array,
2166 const bool sum_into_values = false);
2167
2179 template <typename VectorType>
2180 void
2182 VectorType &output_vector);
2183
2187 template <typename VectorType>
2188 void
2189 integrate_scatter(const bool integrate_values,
2190 const bool integrate_gradients,
2191 VectorType &output_vector);
2192
2200
2205 bool
2207
2222
2227 unsigned int
2229
2234 unsigned int
2236
2243 const unsigned int dofs_per_component;
2244
2251 const unsigned int dofs_per_cell;
2252
2260 const unsigned int n_q_points;
2261};
2262
2263
2264
2265namespace internal
2266{
2267 namespace MatrixFreeFunctions
2268 {
2269 // a helper function to compute the number of DoFs of a DGP element at
2270 // compile time, depending on the degree
2271 template <int dim, int degree>
2273 {
2274 // this division is always without remainder
2275 static constexpr unsigned int value =
2276 (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
2277 };
2278
2279 // base specialization: 1d elements have 'degree+1' degrees of freedom
2280 template <int degree>
2281 struct DGP_dofs_per_component<1, degree>
2282 {
2283 static constexpr unsigned int value = degree + 1;
2284 };
2285 } // namespace MatrixFreeFunctions
2286} // namespace internal
2287
2288
2289/*----------------------- Inline functions ----------------------------------*/
2290
2291#ifndef DOXYGEN
2292
2293
2294namespace internal
2295{
2296 // Extract all internal data pointers and indices in a single function that
2297 // get passed on to the constructor of FEEvaluationData, avoiding to look
2298 // things up multiple times
2299 template <bool is_face,
2300 int dim,
2301 typename Number,
2302 typename VectorizedArrayType>
2304 InitializationData
2305 extract_initialization_data(
2307 const unsigned int dof_no,
2308 const unsigned int first_selected_component,
2309 const unsigned int quad_no,
2310 const unsigned int fe_degree,
2311 const unsigned int n_q_points,
2312 const unsigned int active_fe_index_given,
2313 const unsigned int active_quad_index_given,
2314 const unsigned int face_type)
2315 {
2317 InitializationData init_data;
2318
2319 init_data.dof_info = &matrix_free.get_dof_info(dof_no);
2320 init_data.mapping_data =
2321 &internal::MatrixFreeFunctions::
2322 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2323 matrix_free.get_mapping_info(), quad_no);
2324
2325 init_data.active_fe_index =
2326 fe_degree != numbers::invalid_unsigned_int ?
2327 init_data.dof_info->fe_index_from_degree(first_selected_component,
2328 fe_degree) :
2329 (active_fe_index_given != numbers::invalid_unsigned_int ?
2330 active_fe_index_given :
2331 0);
2332
2333 init_data.active_quad_index =
2334 fe_degree == numbers::invalid_unsigned_int ?
2335 (active_quad_index_given != numbers::invalid_unsigned_int ?
2336 active_quad_index_given :
2337 std::min<unsigned int>(
2338 init_data.active_fe_index,
2339 init_data.mapping_data->descriptor.size() /
2340 (is_face ? std::max<unsigned int>(1, dim - 1) : 1) -
2341 1)) :
2342 init_data.mapping_data->quad_index_from_n_q_points(n_q_points);
2343
2344 init_data.shape_info = &matrix_free.get_shape_info(
2345 dof_no,
2346 quad_no,
2347 init_data.dof_info->component_to_base_index[first_selected_component],
2348 init_data.active_fe_index,
2349 init_data.active_quad_index);
2350 init_data.descriptor =
2351 &init_data.mapping_data->descriptor
2352 [is_face ?
2353 (init_data.active_quad_index * std::max<unsigned int>(1, dim - 1) +
2354 (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
2355 init_data.active_quad_index];
2356
2357 return init_data;
2358 }
2359} // namespace internal
2360
2361
2362
2363/*----------------------- FEEvaluationBase ----------------------------------*/
2364
2365template <int dim,
2366 int n_components_,
2367 typename Number,
2368 bool is_face,
2369 typename VectorizedArrayType>
2370inline FEEvaluationBase<dim,
2371 n_components_,
2372 Number,
2373 is_face,
2374 VectorizedArrayType>::
2375 FEEvaluationBase(
2377 const unsigned int dof_no,
2378 const unsigned int first_selected_component,
2379 const unsigned int quad_no,
2380 const unsigned int fe_degree,
2381 const unsigned int n_q_points,
2382 const bool is_interior_face,
2383 const unsigned int active_fe_index,
2384 const unsigned int active_quad_index,
2385 const unsigned int face_type)
2386 : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2387 internal::extract_initialization_data<is_face>(matrix_free,
2388 dof_no,
2389 first_selected_component,
2390 quad_no,
2391 fe_degree,
2392 n_q_points,
2393 active_fe_index,
2394 active_quad_index,
2395 face_type),
2396 is_interior_face,
2397 quad_no,
2398 first_selected_component)
2399 , scratch_data_array(matrix_free.acquire_scratch_data())
2400 , matrix_free(&matrix_free)
2401{
2402 this->set_data_pointers(scratch_data_array, n_components_);
2403 Assert(
2404 this->dof_info->start_components.back() == 1 ||
2405 static_cast<int>(n_components_) <=
2406 static_cast<int>(
2407 this->dof_info->start_components
2408 [this->dof_info->component_to_base_index[first_selected_component] +
2409 1]) -
2410 first_selected_component,
2411 ExcMessage(
2412 "You tried to construct a vector-valued evaluator with " +
2413 std::to_string(n_components) +
2414 " components. However, "
2415 "the current base element has only " +
2416 std::to_string(
2417 this->dof_info->start_components
2418 [this->dof_info->component_to_base_index[first_selected_component] +
2419 1] -
2420 first_selected_component) +
2421 " components left when starting from local element index " +
2422 std::to_string(
2423 first_selected_component -
2424 this->dof_info->start_components
2425 [this->dof_info->component_to_base_index[first_selected_component]]) +
2426 " (global index " + std::to_string(first_selected_component) + ")"));
2427
2428 // do not check for correct dimensions of data fields here, should be done
2429 // in derived classes
2430}
2431
2432
2433
2434template <int dim,
2435 int n_components_,
2436 typename Number,
2437 bool is_face,
2438 typename VectorizedArrayType>
2439inline FEEvaluationBase<dim,
2440 n_components_,
2441 Number,
2442 is_face,
2443 VectorizedArrayType>::
2444 FEEvaluationBase(
2445 const Mapping<dim> &mapping,
2446 const FiniteElement<dim> &fe,
2447 const Quadrature<1> &quadrature,
2448 const UpdateFlags update_flags,
2449 const unsigned int first_selected_component,
2451 : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2452 other != nullptr &&
2453 other->mapped_geometry->get_quadrature() == quadrature ?
2454 other->mapped_geometry :
2455 std::make_shared<internal::MatrixFreeFunctions::
2456 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2457 mapping,
2458 quadrature,
2459 update_flags),
2460 n_components_,
2461 first_selected_component)
2462 , scratch_data_array(new AlignedVector<VectorizedArrayType>())
2463 , matrix_free(nullptr)
2464{
2465 const unsigned int base_element_number =
2466 fe.component_to_base_index(first_selected_component).first;
2467 Assert(fe.element_multiplicity(base_element_number) == 1 ||
2468 fe.element_multiplicity(base_element_number) -
2469 first_selected_component >=
2470 n_components_,
2471 ExcMessage("The underlying element must at least contain as many "
2472 "components as requested by this class"));
2473 (void)base_element_number;
2474
2475 Assert(this->data == nullptr, ExcInternalError());
2477 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2478 fe,
2479 fe.component_to_base_index(first_selected_component).first);
2480
2481 this->set_data_pointers(scratch_data_array, n_components_);
2482}
2483
2484
2485
2486template <int dim,
2487 int n_components_,
2488 typename Number,
2489 bool is_face,
2490 typename VectorizedArrayType>
2491inline FEEvaluationBase<dim,
2492 n_components_,
2493 Number,
2494 is_face,
2495 VectorizedArrayType>::
2496 FEEvaluationBase(const FEEvaluationBase<dim,
2497 n_components_,
2498 Number,
2499 is_face,
2500 VectorizedArrayType> &other)
2501 : FEEvaluationData<dim, VectorizedArrayType, is_face>(other)
2502 , scratch_data_array(other.matrix_free == nullptr ?
2503 new AlignedVector<VectorizedArrayType>() :
2504 other.matrix_free->acquire_scratch_data())
2505 , matrix_free(other.matrix_free)
2506{
2507 if (other.matrix_free == nullptr)
2508 {
2509 Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
2510 this->data =
2512
2513 // Create deep copy of mapped geometry for use in parallel
2514 this->mapped_geometry =
2515 std::make_shared<internal::MatrixFreeFunctions::
2516 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2517 other.mapped_geometry->get_fe_values().get_mapping(),
2518 other.mapped_geometry->get_quadrature(),
2519 other.mapped_geometry->get_fe_values().get_update_flags());
2520
2521 if constexpr (is_face == false)
2522 this->mapping_data = &this->mapped_geometry->get_data_storage();
2523 else
2524 Assert(false,
2525 ExcNotImplemented("On-the-fly geometry evaluation with "
2526 "face evaluators is not currently "
2527 "implemented!"));
2528
2529 this->cell = 0;
2530
2531 this->jacobian =
2532 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2533 this->J_value =
2534 this->mapped_geometry->get_data_storage().JxW_values.begin();
2535 this->jacobian_gradients =
2536 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2537 this->jacobian_gradients_non_inverse =
2538 this->mapped_geometry->get_data_storage()
2539 .jacobian_gradients_non_inverse[0]
2540 .begin();
2541 this->quadrature_points =
2542 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2543 }
2544
2545 this->set_data_pointers(scratch_data_array, n_components_);
2546}
2547
2548
2549
2550template <int dim,
2551 int n_components_,
2552 typename Number,
2553 bool is_face,
2554 typename VectorizedArrayType>
2555inline FEEvaluationBase<dim,
2556 n_components_,
2557 Number,
2558 is_face,
2559 VectorizedArrayType> &
2561operator=(const FEEvaluationBase<dim,
2562 n_components_,
2563 Number,
2564 is_face,
2565 VectorizedArrayType> &other)
2566{
2567 // release old memory
2568 if (matrix_free == nullptr)
2569 {
2570 delete this->data;
2571 delete scratch_data_array;
2572 }
2573 else
2574 {
2575 matrix_free->release_scratch_data(scratch_data_array);
2576 }
2577
2579
2580 matrix_free = other.matrix_free;
2581
2582 if (other.matrix_free == nullptr)
2583 {
2584 Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
2585 this->data =
2587 scratch_data_array = new AlignedVector<VectorizedArrayType>();
2588
2589 // Create deep copy of mapped geometry for use in parallel
2590 this->mapped_geometry =
2591 std::make_shared<internal::MatrixFreeFunctions::
2592 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2593 other.mapped_geometry->get_fe_values().get_mapping(),
2594 other.mapped_geometry->get_quadrature(),
2595 other.mapped_geometry->get_fe_values().get_update_flags());
2596
2597 if constexpr (is_face == false)
2598 this->mapping_data = &this->mapped_geometry->get_data_storage();
2599 else
2600 Assert(false,
2601 ExcNotImplemented("On-the-fly geometry evaluation with "
2602 "face evaluators is not currently "
2603 "implemented!"));
2604 this->cell = 0;
2605
2606 this->jacobian =
2607 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2608 this->J_value =
2609 this->mapped_geometry->get_data_storage().JxW_values.begin();
2610 this->jacobian_gradients =
2611 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2612 this->jacobian_gradients_non_inverse =
2613 this->mapped_geometry->get_data_storage()
2614 .jacobian_gradients_non_inverse[0]
2615 .begin();
2616 this->quadrature_points =
2617 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2618 }
2619 else
2620 {
2621 scratch_data_array = matrix_free->acquire_scratch_data();
2622 }
2623
2624 this->set_data_pointers(scratch_data_array, n_components_);
2625
2626 return *this;
2627}
2628
2629
2630
2631template <int dim,
2632 int n_components_,
2633 typename Number,
2634 bool is_face,
2635 typename VectorizedArrayType>
2636inline FEEvaluationBase<dim,
2637 n_components_,
2638 Number,
2639 is_face,
2640 VectorizedArrayType>::~FEEvaluationBase()
2641{
2642 if (matrix_free != nullptr)
2643 {
2644 try
2645 {
2646 matrix_free->release_scratch_data(scratch_data_array);
2647 }
2648 catch (...)
2649 {}
2650 }
2651 else
2652 {
2653 delete scratch_data_array;
2654 delete this->data;
2655 }
2656}
2657
2658
2659
2660template <int dim,
2661 int n_components_,
2662 typename Number,
2663 bool is_face,
2664 typename VectorizedArrayType>
2667 get_matrix_free() const
2668{
2669 Assert(matrix_free != nullptr,
2670 ExcMessage(
2671 "FEEvaluation was not initialized with a MatrixFree object!"));
2672 return *matrix_free;
2673}
2674
2675
2676
2677namespace internal
2678{
2679 // given a block vector return the underlying vector type
2680 // including constness (specified by bool)
2681 template <typename VectorType, bool>
2682 struct ConstBlockVectorSelector;
2683
2684 template <typename VectorType>
2685 struct ConstBlockVectorSelector<VectorType, true>
2686 {
2687 using BaseVectorType = const typename VectorType::BlockType;
2688 };
2689
2690 template <typename VectorType>
2691 struct ConstBlockVectorSelector<VectorType, false>
2692 {
2693 using BaseVectorType = typename VectorType::BlockType;
2694 };
2695
2696 // allows to select between block vectors and non-block vectors, which
2697 // allows to use a unified interface for extracting blocks on block vectors
2698 // and doing nothing on usual vectors
2699 template <typename VectorType, bool>
2700 struct BlockVectorSelector;
2701
2702 template <typename VectorType>
2703 struct BlockVectorSelector<VectorType, true>
2704 {
2705 using BaseVectorType = typename ConstBlockVectorSelector<
2706 VectorType,
2707 std::is_const_v<VectorType>>::BaseVectorType;
2708
2709 static BaseVectorType *
2710 get_vector_component(VectorType &vec, const unsigned int component)
2711 {
2712 AssertIndexRange(component, vec.n_blocks());
2713 return &vec.block(component);
2714 }
2715 };
2716
2717 template <typename VectorType>
2718 struct BlockVectorSelector<VectorType, false>
2719 {
2720 using BaseVectorType = VectorType;
2721
2722 static BaseVectorType *
2723 get_vector_component(VectorType &vec, const unsigned int component)
2724 {
2725 // FEEvaluation allows to combine several vectors from a scalar
2726 // FiniteElement into a "vector-valued" FEEvaluation object with
2727 // multiple components. These components can be extracted with the other
2728 // get_vector_component functions. If we do not get a vector of vectors
2729 // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
2730 // must make sure that we do not duplicate the components in input
2731 // and/or duplicate the resulting integrals. In such a case, we should
2732 // only get the zeroth component in the vector contained set nullptr for
2733 // the others which allows us to catch unintended use in
2734 // read_write_operation.
2735 if (component == 0)
2736 return &vec;
2737 else
2738 return nullptr;
2739 }
2740 };
2741
2742 template <typename VectorType>
2743 struct BlockVectorSelector<std::vector<VectorType>, false>
2744 {
2745 using BaseVectorType = VectorType;
2746
2747 static BaseVectorType *
2748 get_vector_component(std::vector<VectorType> &vec,
2749 const unsigned int component)
2750 {
2751 AssertIndexRange(component, vec.size());
2752 return &vec[component];
2753 }
2754 };
2755
2756 template <typename VectorType>
2757 struct BlockVectorSelector<const std::vector<VectorType>, false>
2758 {
2759 using BaseVectorType = const VectorType;
2760
2761 static const BaseVectorType *
2762 get_vector_component(const std::vector<VectorType> &vec,
2763 const unsigned int component)
2764 {
2765 AssertIndexRange(component, vec.size());
2766 return &vec[component];
2767 }
2768 };
2769
2770 template <typename VectorType>
2771 struct BlockVectorSelector<std::vector<VectorType *>, false>
2772 {
2773 using BaseVectorType = VectorType;
2774
2775 static BaseVectorType *
2776 get_vector_component(std::vector<VectorType *> &vec,
2777 const unsigned int component)
2778 {
2779 AssertIndexRange(component, vec.size());
2780 return vec[component];
2781 }
2782 };
2783
2784 template <typename VectorType>
2785 struct BlockVectorSelector<const std::vector<VectorType *>, false>
2786 {
2787 using BaseVectorType = const VectorType;
2788
2789 static const BaseVectorType *
2790 get_vector_component(const std::vector<VectorType *> &vec,
2791 const unsigned int component)
2792 {
2793 AssertIndexRange(component, vec.size());
2794 return vec[component];
2795 }
2796 };
2797
2798 template <typename VectorType, std::size_t N>
2799 struct BlockVectorSelector<std::array<VectorType *, N>, false>
2800 {
2801 using BaseVectorType = VectorType;
2802
2803 static BaseVectorType *
2804 get_vector_component(std::array<VectorType *, N> &vec,
2805 const unsigned int component)
2806 {
2807 AssertIndexRange(component, vec.size());
2808 return vec[component];
2809 }
2810 };
2811} // namespace internal
2812
2813
2814
2815template <int dim,
2816 int n_components_,
2817 typename Number,
2818 bool is_face,
2819 typename VectorizedArrayType>
2820template <typename VectorType, typename VectorOperation>
2821inline void
2824 const VectorOperation &operation,
2825 const std::array<VectorType *, n_components_> &src,
2826 const std::array<
2828 n_components_> &src_sm,
2829 const std::bitset<n_lanes> &mask,
2830 const bool apply_constraints) const
2831{
2832 // Case 1: No MatrixFree object given, simple case because we do not need to
2833 // process constraints and need not care about vectorization -> go to
2834 // separate function
2835 if (this->matrix_free == nullptr)
2836 {
2837 read_write_operation_global(operation, src);
2838 return;
2839 }
2840
2841 Assert(this->dof_info != nullptr, ExcNotInitialized());
2842 const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
2843 Assert(this->matrix_free->indices_initialized() == true, ExcNotInitialized());
2844 if (this->n_fe_components == 1)
2845 for (unsigned int comp = 0; comp < n_components; ++comp)
2846 {
2847 Assert(src[comp] != nullptr,
2848 ExcMessage("The finite element underlying this FEEvaluation "
2849 "object is scalar, but you requested " +
2850 std::to_string(n_components) +
2851 " components via the template argument in "
2852 "FEEvaluation. In that case, you must pass an "
2853 "std::vector<VectorType> or a BlockVector to " +
2854 "read_dof_values and distribute_local_to_global."));
2856 *this->matrix_free,
2857 *this->dof_info);
2858 }
2859 else
2860 {
2862 *this->matrix_free,
2863 *this->dof_info);
2864 }
2865
2866 const bool accesses_exterior_dofs =
2867 this->dof_access_index ==
2869 this->is_interior_face() == false;
2870
2871 // Case 2: contiguous indices which use reduced storage of indices and can
2872 // use vectorized load/store operations -> go to separate function
2873 if (this->cell != numbers::invalid_unsigned_int)
2874 {
2876 this->cell,
2877 dof_info.index_storage_variants[this->dof_access_index].size());
2878
2879 bool is_contiguous = true;
2880 // check if exterior cells are not contiguous (ECL case)
2881 if (accesses_exterior_dofs)
2882 {
2883 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
2884 const unsigned int n_filled_lanes =
2887 [this->cell];
2888 // we have to check all filled lanes which are active in the mask
2889 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2890 if (mask[v] == true &&
2891 dof_info.index_storage_variants
2893 [cells[v] / n_lanes] <
2895 contiguous)
2896 is_contiguous = false;
2897 } // or if cell/face batch is not contiguous
2898 else if (dof_info.index_storage_variants
2899 [is_face ?
2900 this->dof_access_index :
2901 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2902 [this->cell] < internal::MatrixFreeFunctions::DoFInfo::
2903 IndexStorageVariants::contiguous)
2904 {
2905 is_contiguous = false;
2906 }
2907
2908 if (is_contiguous)
2909 {
2910 read_write_operation_contiguous(operation, src, src_sm, mask);
2911 return;
2912 }
2913 }
2914
2915 // Case 3: standard operation with one index per degree of freedom -> go on
2916 // here
2917 std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
2918
2919 const bool masking_is_active = mask.count() < n_lanes;
2920 if (masking_is_active)
2921 for (unsigned int v = 0; v < n_lanes; ++v)
2922 if (mask[v] == false)
2924
2925 bool has_hn_constraints = false;
2926
2927 if (is_face == false)
2928 {
2929 if (!dof_info.hanging_node_constraint_masks.empty() &&
2930 !dof_info.hanging_node_constraint_masks_comp.empty() &&
2931 dof_info
2932 .hanging_node_constraint_masks_comp[this->active_fe_index]
2933 [this->first_selected_component])
2934 for (unsigned int v = 0; v < n_lanes; ++v)
2935 if (cells[v] != numbers::invalid_unsigned_int &&
2936 dof_info.hanging_node_constraint_masks[cells[v]] !=
2939 has_hn_constraints = true;
2940 }
2941
2942 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
2943 vector_selector;
2944
2945 const bool use_vectorized_path =
2946 !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
2947
2948 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
2949 std::array<VectorizedArrayType *, n_components> values_dofs;
2950 for (unsigned int c = 0; c < n_components; ++c)
2951 values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
2952 c * dofs_per_component;
2953
2954 if (this->cell != numbers::invalid_unsigned_int &&
2955 dof_info.index_storage_variants
2956 [is_face ? this->dof_access_index :
2957 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2958 [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
2959 IndexStorageVariants::interleaved &&
2960 use_vectorized_path)
2961 {
2962 const unsigned int *dof_indices =
2963 dof_info.dof_indices_interleaved.data() +
2964 dof_info.row_starts[this->cell * this->n_fe_components * n_lanes]
2965 .first +
2966 this->dof_info
2967 ->component_dof_indices_offset[this->active_fe_index]
2968 [this->first_selected_component] *
2969 n_lanes;
2970
2971 std::array<typename VectorType::value_type *, n_components> src_ptrs;
2972 if (n_components == 1 || this->n_fe_components == 1)
2973 for (unsigned int comp = 0; comp < n_components; ++comp)
2974 src_ptrs[comp] =
2975 const_cast<typename VectorType::value_type *>(src[comp]->begin());
2976 else
2977 src_ptrs[0] =
2978 const_cast<typename VectorType::value_type *>(src[0]->begin());
2979
2980 if (n_components == 1 || this->n_fe_components == 1)
2981 for (unsigned int i = 0; i < dofs_per_component;
2982 ++i, dof_indices += n_lanes)
2983 for (unsigned int comp = 0; comp < n_components; ++comp)
2984 operation.process_dof_gather(dof_indices,
2985 *src[comp],
2986 0,
2987 src_ptrs[comp],
2988 values_dofs[comp][i],
2989 vector_selector);
2990 else
2991 for (unsigned int comp = 0; comp < n_components; ++comp)
2992 for (unsigned int i = 0; i < dofs_per_component;
2993 ++i, dof_indices += n_lanes)
2994 operation.process_dof_gather(dof_indices,
2995 *src[0],
2996 0,
2997 src_ptrs[0],
2998 values_dofs[comp][i],
2999 vector_selector);
3000 return;
3001 }
3002
3003 // Allocate pointers, then initialize all of them to nullptrs and
3004 // below overwrite the ones we actually use:
3005 std::array<const unsigned int *, n_lanes> dof_indices;
3006 dof_indices.fill(nullptr);
3007
3008 // Assign the appropriate cell ids for face/cell case and get the pointers
3009 // to the dof indices of the cells on all lanes
3010
3011 bool has_constraints = false;
3012 const unsigned int n_components_read =
3013 this->n_fe_components > 1 ? n_components : 1;
3014
3015 if (is_face)
3016 {
3017 for (unsigned int v = 0; v < n_lanes; ++v)
3018 {
3019 if (cells[v] == numbers::invalid_unsigned_int)
3020 continue;
3021
3022 Assert(cells[v] < dof_info.row_starts.size() - 1, ExcInternalError());
3023 const std::pair<unsigned int, unsigned int> *my_index_start =
3024 &dof_info.row_starts[cells[v] * this->n_fe_components +
3025 this->first_selected_component];
3026
3027 // check whether any of the SIMD lanes has constraints, i.e., the
3028 // constraint indicator which is the second entry of row_starts
3029 // increments on this cell
3030 if (my_index_start[n_components_read].second !=
3031 my_index_start[0].second)
3032 has_constraints = true;
3033
3034 dof_indices[v] =
3035 dof_info.dof_indices.data() + my_index_start[0].first;
3036 }
3037 }
3038 else
3039 {
3040 for (unsigned int v = 0; v < n_lanes; ++v)
3041 {
3042 if (cells[v] == numbers::invalid_unsigned_int)
3043 continue;
3044
3045 const std::pair<unsigned int, unsigned int> *my_index_start =
3046 &dof_info.row_starts[cells[v] * this->n_fe_components +
3047 this->first_selected_component];
3048 if (my_index_start[n_components_read].second !=
3049 my_index_start[0].second)
3050 has_constraints = true;
3051
3052 if (dof_info.hanging_node_constraint_masks.size() > 0 &&
3053 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
3054 dof_info.hanging_node_constraint_masks[cells[v]] !=
3057 dof_info.hanging_node_constraint_masks_comp
3058 [this->active_fe_index][this->first_selected_component])
3059 has_hn_constraints = true;
3060
3061 Assert(my_index_start[n_components_read].first ==
3062 my_index_start[0].first ||
3063 my_index_start[0].first < dof_info.dof_indices.size(),
3064 ExcIndexRange(0,
3065 my_index_start[0].first,
3066 dof_info.dof_indices.size()));
3067 dof_indices[v] =
3068 dof_info.dof_indices.data() + my_index_start[0].first;
3069 }
3070 }
3071
3072 if (std::count_if(cells.begin(), cells.end(), [](const auto i) {
3073 return i != numbers::invalid_unsigned_int;
3074 }) < n_lanes)
3075 for (unsigned int comp = 0; comp < n_components; ++comp)
3076 for (unsigned int i = 0; i < dofs_per_component; ++i)
3077 operation.process_empty(values_dofs[comp][i]);
3078
3079 // Case where we have no constraints throughout the whole cell: Can go
3080 // through the list of DoFs directly
3081 if (!has_constraints && apply_constraints)
3082 {
3083 if (n_components == 1 || this->n_fe_components == 1)
3084 {
3085 for (unsigned int v = 0; v < n_lanes; ++v)
3086 {
3087 if (cells[v] == numbers::invalid_unsigned_int)
3088 continue;
3089
3090 for (unsigned int i = 0; i < dofs_per_component; ++i)
3091 for (unsigned int comp = 0; comp < n_components; ++comp)
3092 operation.process_dof(dof_indices[v][i],
3093 *src[comp],
3094 values_dofs[comp][i][v]);
3095 }
3096 }
3097 else
3098 {
3099 for (unsigned int comp = 0; comp < n_components; ++comp)
3100 for (unsigned int v = 0; v < n_lanes; ++v)
3101 {
3102 if (cells[v] == numbers::invalid_unsigned_int)
3103 continue;
3104
3105 for (unsigned int i = 0; i < dofs_per_component; ++i)
3106 operation.process_dof(
3107 dof_indices[v][comp * dofs_per_component + i],
3108 *src[0],
3109 values_dofs[comp][i][v]);
3110 }
3111 }
3112 return;
3113 }
3114
3115 // In the case where there are some constraints to be resolved, loop over
3116 // all vector components that are filled and then over local dofs. ind_local
3117 // holds local number on cell, index iterates over the elements of
3118 // index_local_to_global and dof_indices points to the global indices stored
3119 // in index_local_to_global
3120
3121 for (unsigned int v = 0; v < n_lanes; ++v)
3122 {
3123 if (cells[v] == numbers::invalid_unsigned_int)
3124 continue;
3125
3126 const unsigned int cell_index = cells[v];
3127 const unsigned int cell_dof_index =
3128 cell_index * this->n_fe_components + this->first_selected_component;
3129 const unsigned int n_components_read =
3130 this->n_fe_components > 1 ? n_components : 1;
3131 unsigned int index_indicators =
3132 dof_info.row_starts[cell_dof_index].second;
3133 unsigned int next_index_indicators =
3134 dof_info.row_starts[cell_dof_index + 1].second;
3135
3136 // For read_dof_values_plain, redirect the dof_indices field to the
3137 // unconstrained indices
3138 if (apply_constraints == false &&
3139 (dof_info.row_starts[cell_dof_index].second !=
3140 dof_info.row_starts[cell_dof_index + n_components_read].second ||
3141 ((dof_info.hanging_node_constraint_masks.size() > 0 &&
3142 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
3146 dof_info.hanging_node_constraint_masks_comp
3147 [this->active_fe_index][this->first_selected_component])))
3148 {
3152 dof_indices[v] =
3153 dof_info.plain_dof_indices.data() +
3154 this->dof_info
3155 ->component_dof_indices_offset[this->active_fe_index]
3156 [this->first_selected_component] +
3158 next_index_indicators = index_indicators;
3159 }
3160
3161 if (n_components == 1 || this->n_fe_components == 1)
3162 {
3163 unsigned int ind_local = 0;
3164 for (; index_indicators != next_index_indicators; ++index_indicators)
3165 {
3166 const std::pair<unsigned short, unsigned short> indicator =
3167 dof_info.constraint_indicator[index_indicators];
3168 // run through values up to next constraint
3169 for (unsigned int j = 0; j < indicator.first; ++j)
3170 for (unsigned int comp = 0; comp < n_components; ++comp)
3171 operation.process_dof(dof_indices[v][j],
3172 *src[comp],
3173 values_dofs[comp][ind_local + j][v]);
3174
3175 ind_local += indicator.first;
3176 dof_indices[v] += indicator.first;
3177
3178 // constrained case: build the local value as a linear
3179 // combination of the global value according to constraints
3180 Number value[n_components];
3181 for (unsigned int comp = 0; comp < n_components; ++comp)
3182 operation.pre_constraints(values_dofs[comp][ind_local][v],
3183 value[comp]);
3184
3185 const Number *data_val =
3186 this->matrix_free->constraint_pool_begin(indicator.second);
3187 const Number *end_pool =
3188 this->matrix_free->constraint_pool_end(indicator.second);
3189 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3190 for (unsigned int comp = 0; comp < n_components; ++comp)
3191 operation.process_constraint(*dof_indices[v],
3192 *data_val,
3193 *src[comp],
3194 value[comp]);
3195
3196 for (unsigned int comp = 0; comp < n_components; ++comp)
3197 operation.post_constraints(value[comp],
3198 values_dofs[comp][ind_local][v]);
3199 ++ind_local;
3200 }
3201
3202 AssertIndexRange(ind_local, dofs_per_component + 1);
3203
3204 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3205 for (unsigned int comp = 0; comp < n_components; ++comp)
3206 operation.process_dof(*dof_indices[v],
3207 *src[comp],
3208 values_dofs[comp][ind_local][v]);
3209 }
3210 else
3211 {
3212 // case with vector-valued finite elements where all components are
3213 // included in one single vector. Assumption: first come all entries
3214 // to the first component, then all entries to the second one, and
3215 // so on. This is ensured by the way MatrixFree reads out the
3216 // indices.
3217 for (unsigned int comp = 0; comp < n_components; ++comp)
3218 {
3219 unsigned int ind_local = 0;
3220
3221 // check whether there is any constraint on the current cell
3222 for (; index_indicators != next_index_indicators;
3223 ++index_indicators)
3224 {
3225 const std::pair<unsigned short, unsigned short> indicator =
3226 dof_info.constraint_indicator[index_indicators];
3227
3228 // run through values up to next constraint
3229 for (unsigned int j = 0; j < indicator.first; ++j)
3230 operation.process_dof(dof_indices[v][j],
3231 *src[0],
3232 values_dofs[comp][ind_local + j][v]);
3233 ind_local += indicator.first;
3234 dof_indices[v] += indicator.first;
3235
3236 // constrained case: build the local value as a linear
3237 // combination of the global value according to constraints
3238 Number value;
3239 operation.pre_constraints(values_dofs[comp][ind_local][v],
3240 value);
3241
3242 const Number *data_val =
3243 this->matrix_free->constraint_pool_begin(indicator.second);
3244 const Number *end_pool =
3245 this->matrix_free->constraint_pool_end(indicator.second);
3246
3247 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3248 operation.process_constraint(*dof_indices[v],
3249 *data_val,
3250 *src[0],
3251 value);
3252
3253 operation.post_constraints(value,
3254 values_dofs[comp][ind_local][v]);
3255 ++ind_local;
3256 }
3257
3258 AssertIndexRange(ind_local, dofs_per_component + 1);
3259
3260 // get the dof values past the last constraint
3261 for (; ind_local < dofs_per_component;
3262 ++dof_indices[v], ++ind_local)
3263 {
3264 AssertIndexRange(*dof_indices[v], src[0]->size());
3265 operation.process_dof(*dof_indices[v],
3266 *src[0],
3267 values_dofs[comp][ind_local][v]);
3268 }
3269
3270 if (apply_constraints == true && comp + 1 < n_components)
3271 next_index_indicators =
3272 dof_info.row_starts[cell_dof_index + comp + 2].second;
3273 }
3274 }
3275 }
3276}
3277
3278
3279
3280template <int dim,
3281 int n_components_,
3282 typename Number,
3283 bool is_face,
3284 typename VectorizedArrayType>
3285template <typename VectorType, typename VectorOperation>
3286inline void
3289 const VectorOperation &operation,
3290 const std::array<VectorType *, n_components_> &src) const
3291{
3292 Assert(!local_dof_indices.empty(), ExcNotInitialized());
3293
3294 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3295 unsigned int index = this->first_selected_component * dofs_per_component;
3296 for (unsigned int comp = 0; comp < n_components; ++comp)
3297 {
3298 for (unsigned int i = 0; i < dofs_per_component; ++i, ++index)
3299 {
3300 operation.process_empty(
3301 this->values_dofs[comp * dofs_per_component + i]);
3302 operation.process_dof_global(
3303 local_dof_indices[this->data->lexicographic_numbering[index]],
3304 *src[0],
3305 this->values_dofs[comp * dofs_per_component + i][0]);
3306 }
3307 }
3308}
3309
3310
3311
3312template <int dim,
3313 int n_components_,
3314 typename Number,
3315 bool is_face,
3316 typename VectorizedArrayType>
3317template <typename VectorType, typename VectorOperation>
3318inline void
3321 const VectorOperation &operation,
3322 const std::array<VectorType *, n_components_> &src,
3323 const std::array<
3325 n_components_> &vectors_sm,
3326 const std::bitset<n_lanes> &mask) const
3327{
3328 // This functions processes the functions read_dof_values,
3329 // distribute_local_to_global, and set_dof_values with the same code for
3330 // contiguous cell indices (DG case). The distinction between these three
3331 // cases is made by the input VectorOperation that either reads values from
3332 // a vector and puts the data into the local data field or write local data
3333 // into the vector. Certain operations are no-ops for the given use case.
3334
3335 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
3336 vector_selector;
3338 is_face ? this->dof_access_index :
3340 const unsigned int n_active_lanes = mask.count();
3341
3342 const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
3343 const std::vector<unsigned int> &dof_indices_cont =
3344 dof_info.dof_indices_contiguous[ind];
3345
3346 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3347 std::array<VectorizedArrayType *, n_components> values_dofs{{nullptr}};
3348 for (unsigned int c = 0; c < n_components; ++c)
3349 values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3350 c * dofs_per_component;
3351
3353
3354 const bool accesses_exterior_dofs =
3355 this->dof_access_index ==
3357 this->is_interior_face() == false;
3358
3359 // Simple case: We have contiguous storage, so we can simply copy out the
3360 // data
3361 if (dof_info.index_storage_variants[ind][this->cell] ==
3363 interleaved_contiguous &&
3364 n_active_lanes == n_lanes && !accesses_exterior_dofs)
3365 {
3366 const unsigned int dof_index =
3367 dof_indices_cont[this->cell * n_lanes] +
3368 this->dof_info
3369 ->component_dof_indices_offset[this->active_fe_index]
3370 [this->first_selected_component] *
3371 n_lanes;
3372 if (n_components == 1 || this->n_fe_components == 1)
3373 for (unsigned int comp = 0; comp < n_components; ++comp)
3374 operation.process_dofs_vectorized(dofs_per_component,
3375 dof_index,
3376 *src[comp],
3377 values_dofs[comp],
3378 vector_selector);
3379 else
3380 operation.process_dofs_vectorized(dofs_per_component * n_components,
3381 dof_index,
3382 *src[0],
3383 values_dofs[0],
3384 vector_selector);
3385 return;
3386 }
3387
3388 const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3389
3390 // More general case: Must go through the components one by one and apply
3391 // some transformations
3392 const unsigned int n_filled_lanes =
3393 dof_info.n_vectorization_lanes_filled[ind][this->cell];
3394
3395 const bool use_vectorized_path = n_filled_lanes == n_lanes &&
3396 n_active_lanes == n_lanes &&
3397 !accesses_exterior_dofs;
3398
3399 if (vectors_sm[0] != nullptr)
3400 {
3401 const auto compute_vector_ptrs = [&](const unsigned int comp) {
3402 std::array<typename VectorType::value_type *, n_lanes> vector_ptrs{
3403 {nullptr}};
3404
3405 const auto upper_bound =
3406 std::min<unsigned int>(n_filled_lanes, n_lanes);
3407 for (unsigned int v = 0; v < upper_bound; ++v)
3408 {
3409 if (mask[v] == false)
3410 {
3411 vector_ptrs[v] = nullptr;
3412 continue;
3413 }
3414
3417 Assert(ind < dof_info.dof_indices_contiguous_sm.size(),
3418 ExcIndexRange(ind,
3419 0,
3420 dof_info.dof_indices_contiguous_sm.size()));
3421 Assert(
3422 cells[v] < dof_info.dof_indices_contiguous_sm[ind].size(),
3423 ExcIndexRange(cells[v],
3424 0,
3425 dof_info.dof_indices_contiguous_sm[ind].size()));
3426
3427 const auto &temp =
3428 dof_info.dof_indices_contiguous_sm[ind][cells[v]];
3429
3430 if (temp.first != numbers::invalid_unsigned_int)
3431 vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
3432 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3434 [this->active_fe_index][this->first_selected_component]);
3435 else
3436 vector_ptrs[v] = nullptr;
3437 }
3438 for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
3439 vector_ptrs[v] = nullptr;
3440
3441 return vector_ptrs;
3442 };
3443
3444 if (use_vectorized_path)
3445 {
3446 if (n_components == 1 || this->n_fe_components == 1)
3447 {
3448 for (unsigned int comp = 0; comp < n_components; ++comp)
3449 {
3450 auto vector_ptrs = compute_vector_ptrs(comp);
3451 operation.process_dofs_vectorized_transpose(
3452 dofs_per_component,
3453 vector_ptrs,
3454 values_dofs[comp],
3455 vector_selector);
3456 }
3457 }
3458 else
3459 {
3460 auto vector_ptrs = compute_vector_ptrs(0);
3461 operation.process_dofs_vectorized_transpose(dofs_per_component *
3462 n_components,
3463 vector_ptrs,
3464 &values_dofs[0][0],
3465 vector_selector);
3466 }
3467 }
3468 else
3469 for (unsigned int comp = 0; comp < n_components; ++comp)
3470 {
3471 auto vector_ptrs = compute_vector_ptrs(
3472 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3473
3474 for (unsigned int i = 0; i < dofs_per_component; ++i)
3475 operation.process_empty(values_dofs[comp][i]);
3476
3477 if (n_components == 1 || this->n_fe_components == 1)
3478 {
3479 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3480 if (mask[v] == true)
3481 for (unsigned int i = 0; i < dofs_per_component; ++i)
3482 operation.process_dof(vector_ptrs[v][i],
3483 values_dofs[comp][i][v]);
3484 }
3485 else
3486 {
3487 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3488 if (mask[v] == true)
3489 for (unsigned int i = 0; i < dofs_per_component; ++i)
3490 operation.process_dof(
3491 vector_ptrs[v][i + comp * dofs_per_component],
3492 values_dofs[comp][i][v]);
3493 }
3494 }
3495 return;
3496 }
3497
3498 std::array<unsigned int, n_lanes> dof_indices{
3500 Assert(n_filled_lanes <= n_lanes, ExcInternalError());
3501 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3502 {
3503 Assert(mask[v] == false || cells[v] != numbers::invalid_unsigned_int,
3505 if (mask[v] == true)
3506 dof_indices[v] =
3507 dof_indices_cont[cells[v]] +
3508 this->dof_info
3509 ->component_dof_indices_offset[this->active_fe_index]
3510 [this->first_selected_component] *
3511 dof_info.dof_indices_interleave_strides[ind][cells[v]];
3512 }
3513
3514 // In the case with contiguous cell indices, we know that there are no
3515 // constraints and that the indices within each element are contiguous
3516 if (use_vectorized_path)
3517 {
3518 if (dof_info.index_storage_variants[ind][this->cell] ==
3520 contiguous)
3521 {
3522 if (n_components == 1 || this->n_fe_components == 1)
3523 for (unsigned int comp = 0; comp < n_components; ++comp)
3524 operation.process_dofs_vectorized_transpose(dofs_per_component,
3525 dof_indices.data(),
3526 *src[comp],
3527 values_dofs[comp],
3528 vector_selector);
3529 else
3530 operation.process_dofs_vectorized_transpose(dofs_per_component *
3531 n_components,
3532 dof_indices.data(),
3533 *src[0],
3534 &values_dofs[0][0],
3535 vector_selector);
3536 }
3537 else if (dof_info.index_storage_variants[ind][this->cell] ==
3539 interleaved_contiguous_strided)
3540 {
3541 std::array<typename VectorType::value_type *, n_components> src_ptrs{
3542 {nullptr}};
3543 if (n_components == 1 || this->n_fe_components == 1)
3544 for (unsigned int comp = 0; comp < n_components; ++comp)
3545 src_ptrs[comp] = const_cast<typename VectorType::value_type *>(
3546 src[comp]->begin());
3547 else
3548 src_ptrs[0] =
3549 const_cast<typename VectorType::value_type *>(src[0]->begin());
3550
3551 if (n_components == 1 || this->n_fe_components == 1)
3552 for (unsigned int i = 0; i < dofs_per_component; ++i)
3553 {
3554 for (unsigned int comp = 0; comp < n_components; ++comp)
3555 operation.process_dof_gather(dof_indices.data(),
3556 *src[comp],
3557 i * n_lanes,
3558 src_ptrs[comp] + i * n_lanes,
3559 values_dofs[comp][i],
3560 vector_selector);
3561 }
3562 else
3563 for (unsigned int comp = 0; comp < n_components; ++comp)
3564 for (unsigned int i = 0; i < dofs_per_component; ++i)
3565 {
3566 operation.process_dof_gather(
3567 dof_indices.data(),
3568 *src[0],
3569 (comp * dofs_per_component + i) * n_lanes,
3570 src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
3571 values_dofs[comp][i],
3572 vector_selector);
3573 }
3574 }
3575 else
3576 {
3577 Assert(dof_info.index_storage_variants[ind][this->cell] ==
3579 IndexStorageVariants::interleaved_contiguous_mixed_strides,
3581 std::array<typename VectorType::value_type *, n_components> src_ptrs{
3582 {nullptr}};
3583 if (n_components == 1 || this->n_fe_components == 1)
3584 for (unsigned int comp = 0; comp < n_components; ++comp)
3585 src_ptrs[comp] = const_cast<typename VectorType::value_type *>(
3586 src[comp]->begin());
3587 else
3588 src_ptrs[0] =
3589 const_cast<typename VectorType::value_type *>(src[0]->begin());
3590
3591 const unsigned int *offsets =
3592 &dof_info.dof_indices_interleave_strides[ind][n_lanes * this->cell];
3593 if (n_components == 1 || this->n_fe_components == 1)
3594 for (unsigned int i = 0; i < dofs_per_component; ++i)
3595 {
3596 for (unsigned int comp = 0; comp < n_components; ++comp)
3597 operation.process_dof_gather(dof_indices.data(),
3598 *src[comp],
3599 0,
3600 src_ptrs[comp],
3601 values_dofs[comp][i],
3602 vector_selector);
3604 for (unsigned int v = 0; v < n_lanes; ++v)
3605 dof_indices[v] += offsets[v];
3606 }
3607 else
3608 for (unsigned int comp = 0; comp < n_components; ++comp)
3609 for (unsigned int i = 0; i < dofs_per_component; ++i)
3610 {
3611 operation.process_dof_gather(dof_indices.data(),
3612 *src[0],
3613 0,
3614 src_ptrs[0],
3615 values_dofs[comp][i],
3616 vector_selector);
3618 for (unsigned int v = 0; v < n_lanes; ++v)
3619 dof_indices[v] += offsets[v];
3620 }
3621 }
3622 }
3623 else
3624 for (unsigned int comp = 0; comp < n_components; ++comp)
3625 {
3626 for (unsigned int i = 0; i < dofs_per_component; ++i)
3627 operation.process_empty(values_dofs[comp][i]);
3628 if (accesses_exterior_dofs)
3629 {
3630 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3631 if (mask[v] == true)
3632 {
3633 if (dof_info.index_storage_variants
3634 [ind][cells[v] / VectorizedArrayType::size()] ==
3637 {
3638 if (n_components == 1 || this->n_fe_components == 1)
3639 {
3640 for (unsigned int i = 0; i < dofs_per_component; ++i)
3641 operation.process_dof(dof_indices[v] + i,
3642 *src[comp],
3643 values_dofs[comp][i][v]);
3644 }
3645 else
3646 {
3647 for (unsigned int i = 0; i < dofs_per_component; ++i)
3648 operation.process_dof(dof_indices[v] + i +
3649 comp * dofs_per_component,
3650 *src[0],
3651 values_dofs[comp][i][v]);
3652 }
3653 }
3654 else
3655 {
3656 const unsigned int offset =
3657 dof_info.dof_indices_interleave_strides[ind][cells[v]];
3658 AssertIndexRange(offset, VectorizedArrayType::size() + 1);
3659 if (n_components == 1 || this->n_fe_components == 1)
3660 {
3661 for (unsigned int i = 0; i < dofs_per_component; ++i)
3662 operation.process_dof(dof_indices[v] + i * offset,
3663 *src[comp],
3664 values_dofs[comp][i][v]);
3665 }
3666 else
3667 {
3668 for (unsigned int i = 0; i < dofs_per_component; ++i)
3669 operation.process_dof(
3670 dof_indices[v] +
3671 (i + comp * dofs_per_component) * offset,
3672 *src[0],
3673 values_dofs[comp][i][v]);
3674 }
3675 }
3676 }
3677 }
3678 else
3679 {
3680 if (dof_info.index_storage_variants[ind][this->cell] ==
3682 contiguous)
3683 {
3684 if (n_components == 1 || this->n_fe_components == 1)
3685 {
3686 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3687 if (mask[v] == true)
3688 for (unsigned int i = 0; i < dofs_per_component; ++i)
3689 operation.process_dof(dof_indices[v] + i,
3690 *src[comp],
3691 values_dofs[comp][i][v]);
3692 }
3693 else
3694 {
3695 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3696 if (mask[v] == true)
3697 for (unsigned int i = 0; i < dofs_per_component; ++i)
3698 operation.process_dof(dof_indices[v] + i +
3699 comp * dofs_per_component,
3700 *src[0],
3701 values_dofs[comp][i][v]);
3702 }
3703 }
3704 else
3705 {
3706 const unsigned int *offsets =
3708 [ind][VectorizedArrayType::size() * this->cell];
3709 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3710 AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
3711 if (n_components == 1 || this->n_fe_components == 1)
3712 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3713 {
3714 if (mask[v] == true)
3715 for (unsigned int i = 0; i < dofs_per_component; ++i)
3716 operation.process_dof(dof_indices[v] + i * offsets[v],
3717 *src[comp],
3718 values_dofs[comp][i][v]);
3719 }
3720 else
3721 {
3722 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3723 if (mask[v] == true)
3724 for (unsigned int i = 0; i < dofs_per_component; ++i)
3725 operation.process_dof(
3726 dof_indices[v] +
3727 (i + comp * dofs_per_component) * offsets[v],
3728 *src[0],
3729 values_dofs[comp][i][v]);
3730 }
3731 }
3732 }
3733 }
3734}
3735
3736namespace internal
3737{
3738 template <
3739 typename Number,
3740 typename VectorType,
3741 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType> * = nullptr>
3742 decltype(std::declval<VectorType>().begin())
3743 get_beginning(VectorType &vec)
3744 {
3745 return vec.begin();
3746 }
3747
3748 template <
3749 typename Number,
3750 typename VectorType,
3751 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * = nullptr>
3752 typename VectorType::value_type *
3753 get_beginning(VectorType &)
3754 {
3755 return nullptr;
3756 }
3757
3758 template <typename VectorType,
3759 std::enable_if_t<has_shared_vector_data<VectorType>, VectorType> * =
3760 nullptr>
3761 const std::vector<ArrayView<const typename VectorType::value_type>> *
3762 get_shared_vector_data(VectorType *vec,
3763 const bool is_valid_mode_for_sm,
3764 const unsigned int active_fe_index,
3766 {
3767 // note: no hp is supported
3768 if (is_valid_mode_for_sm &&
3769 dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
3770 .size() > 0 &&
3771 active_fe_index == 0)
3772 return &vec->shared_vector_data();
3773 else
3774 return nullptr;
3775 }
3776
3777 template <typename VectorType,
3778 std::enable_if_t<!has_shared_vector_data<VectorType>, VectorType>
3779 * = nullptr>
3780 const std::vector<ArrayView<const typename VectorType::value_type>> *
3781 get_shared_vector_data(VectorType *,
3782 const bool,
3783 const unsigned int,
3785 {
3786 return nullptr;
3787 }
3788
3789 template <int n_components, typename VectorType>
3790 std::pair<
3791 std::array<typename internal::BlockVectorSelector<
3792 VectorType,
3793 IsBlockVector<VectorType>::value>::BaseVectorType *,
3794 n_components>,
3795 std::array<
3796 const std::vector<ArrayView<const typename internal::BlockVectorSelector<
3797 VectorType,
3798 IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
3799 n_components>>
3800 get_vector_data(VectorType &src,
3801 const unsigned int first_index,
3802 const bool is_valid_mode_for_sm,
3803 const unsigned int active_fe_index,
3805 {
3806 // select between block vectors and non-block vectors. Note that the number
3807 // of components is checked in the internal data
3808 std::pair<
3809 std::array<typename internal::BlockVectorSelector<
3810 VectorType,
3811 IsBlockVector<VectorType>::value>::BaseVectorType *,
3812 n_components>,
3813 std::array<
3814 const std::vector<
3815 ArrayView<const typename internal::BlockVectorSelector<
3816 VectorType,
3817 IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
3818 n_components>>
3819 src_data;
3820
3821 for (unsigned int d = 0; d < n_components; ++d)
3822 src_data.first[d] = internal::BlockVectorSelector<
3823 VectorType,
3824 IsBlockVector<VectorType>::value>::get_vector_component(src,
3825 d +
3826 first_index);
3827
3828 for (unsigned int d = 0; d < n_components; ++d)
3829 src_data.second[d] = get_shared_vector_data(
3830 const_cast<typename internal::BlockVectorSelector<
3831 std::remove_const_t<VectorType>,
3833 *>(src_data.first[d]),
3834 is_valid_mode_for_sm,
3835 active_fe_index,
3836 dof_info);
3837
3838 return src_data;
3839 }
3840} // namespace internal
3841
3842
3843
3844template <int dim,
3845 int n_components_,
3846 typename Number,
3847 bool is_face,
3848 typename VectorizedArrayType>
3849inline void
3852{
3853 if (this->dof_info == nullptr ||
3854 this->dof_info->hanging_node_constraint_masks.empty() ||
3855 this->dof_info->hanging_node_constraint_masks_comp.empty() ||
3856 this->dof_info->hanging_node_constraint_masks_comp
3857 [this->active_fe_index][this->first_selected_component] == false)
3858 return; // nothing to do with faces
3859
3860 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
3861 constraint_mask{{internal::MatrixFreeFunctions::
3862 unconstrained_compressed_constraint_kind}};
3863
3864 bool hn_available = false;
3865
3866 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3867
3868 for (unsigned int v = 0; v < n_lanes; ++v)
3869 {
3870 if (cells[v] == numbers::invalid_unsigned_int)
3871 {
3872 constraint_mask[v] = internal::MatrixFreeFunctions::
3874 continue;
3875 }
3876
3877 const unsigned int cell_index = cells[v];
3878 const auto mask =
3880 constraint_mask[v] = mask;
3881
3882 hn_available |= (mask != internal::MatrixFreeFunctions::
3884 }
3885
3886 if (hn_available == false)
3887 return; // no hanging node on cell batch -> nothing to do
3888
3890 apply(n_components,
3891 this->data->data.front().fe_degree,
3892 this->get_shape_info(),
3893 transpose,
3894 constraint_mask,
3895 this->values_dofs);
3896}
3897
3898
3899
3900template <int dim,
3901 int n_components_,
3902 typename Number,
3903 bool is_face,
3904 typename VectorizedArrayType>
3905template <typename VectorType>
3906inline void
3908 read_dof_values(const VectorType &src,
3909 const unsigned int first_index,
3910 const std::bitset<n_lanes> &mask)
3911{
3912 const auto src_data = internal::get_vector_data<n_components_>(
3913 src,
3914 first_index,
3915 this->dof_access_index ==
3917 this->active_fe_index,
3918 this->dof_info);
3919
3921 read_write_operation(reader, src_data.first, src_data.second, mask, true);
3922
3923 apply_hanging_node_constraints(false);
3924
3925# ifdef DEBUG
3926 this->dof_values_initialized = true;
3927# endif
3928}
3929
3930
3931
3932template <int dim,
3933 int n_components_,
3934 typename Number,
3935 bool is_face,
3936 typename VectorizedArrayType>
3937template <typename VectorType>
3938inline void
3940 read_dof_values_plain(const VectorType &src,
3941 const unsigned int first_index,
3942 const std::bitset<n_lanes> &mask)
3943{
3944 const auto src_data = internal::get_vector_data<n_components_>(
3945 src,
3946 first_index,
3947 this->dof_access_index ==
3949 this->active_fe_index,
3950 this->dof_info);
3951
3953 read_write_operation(reader, src_data.first, src_data.second, mask, false);
3954
3955# ifdef DEBUG
3956 this->dof_values_initialized = true;
3957# endif
3958}
3959
3960
3961
3962template <int dim,
3963 int n_components_,
3964 typename Number,
3965 bool is_face,
3966 typename VectorizedArrayType>
3967template <typename VectorType>
3968inline void
3970 distribute_local_to_global(VectorType &dst,
3971 const unsigned int first_index,
3972 const std::bitset<n_lanes> &mask) const
3973{
3974# ifdef DEBUG
3975 Assert(this->dof_values_initialized == true,
3977# endif
3978
3979 apply_hanging_node_constraints(true);
3980
3981 const auto dst_data = internal::get_vector_data<n_components_>(
3982 dst,
3983 first_index,
3984 this->dof_access_index ==
3986 this->active_fe_index,
3987 this->dof_info);
3988
3990 distributor;
3991 read_write_operation(distributor, dst_data.first, dst_data.second, mask);
3992}
3993
3994
3995
3996template <int dim,
3997 int n_components_,
3998 typename Number,
3999 bool is_face,
4000 typename VectorizedArrayType>
4001template <typename VectorType>
4002inline void
4004 set_dof_values(VectorType &dst,
4005 const unsigned int first_index,
4006 const std::bitset<n_lanes> &mask) const
4007{
4008# ifdef DEBUG
4009 Assert(this->dof_values_initialized == true,
4011# endif
4012
4013 const auto dst_data = internal::get_vector_data<n_components_>(
4014 dst,
4015 first_index,
4016 this->dof_access_index ==
4018 this->active_fe_index,
4019 this->dof_info);
4020
4022 read_write_operation(setter, dst_data.first, dst_data.second, mask);
4023}
4024
4025
4026
4027template <int dim,
4028 int n_components_,
4029 typename Number,
4030 bool is_face,
4031 typename VectorizedArrayType>
4032template <typename VectorType>
4033inline void
4035 set_dof_values_plain(VectorType &dst,
4036 const unsigned int first_index,
4037 const std::bitset<n_lanes> &mask) const
4038{
4039# ifdef DEBUG
4040 Assert(this->dof_values_initialized == true,
4042# endif
4043
4044 const auto dst_data = internal::get_vector_data<n_components_>(
4045 dst,
4046 first_index,
4047 this->dof_access_index ==
4049 this->active_fe_index,
4050 this->dof_info);
4051
4053 read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
4054}
4055
4056
4057
4058/*------------------------------ access to data fields ----------------------*/
4059
4060
4061
4062template <int dim,
4063 int n_components_,
4064 typename Number,
4065 bool is_face,
4066 typename VectorizedArrayType>
4068 typename FEEvaluationBase<dim,
4069 n_components_,
4070 Number,
4071 is_face,
4072 VectorizedArrayType>::value_type
4074 get_dof_value(const unsigned int dof) const
4075{
4076 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4077 if constexpr (n_components == 1)
4078 return this->values_dofs[dof];
4079 else
4080 {
4081 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4082 Tensor<1, n_components_, VectorizedArrayType> return_value;
4083 for (unsigned int comp = 0; comp < n_components; ++comp)
4084 return_value[comp] = this->values_dofs[comp * dofs + dof];
4085 return return_value;
4086 }
4087}
4088
4089
4090
4091template <int dim,
4092 int n_components_,
4093 typename Number,
4094 bool is_face,
4095 typename VectorizedArrayType>
4097 typename FEEvaluationBase<dim,
4098 n_components_,
4099 Number,
4100 is_face,
4101 VectorizedArrayType>::value_type
4103 get_value(const unsigned int q_point) const
4104{
4105# ifdef DEBUG
4106 Assert(this->values_quad_initialized == true,
4108# endif
4109
4110 AssertIndexRange(q_point, this->n_quadrature_points);
4111 if constexpr (n_components == 1)
4112 return this->values_quad[q_point];
4113 else
4114 {
4115 if (n_components == dim &&
4116 this->data->element_type ==
4118 {
4119 // Piola transform is required
4120# ifdef DEBUG
4121 Assert(this->values_quad_initialized == true,
4123# endif
4124
4125 AssertIndexRange(q_point, this->n_quadrature_points);
4126 Assert(this->J_value != nullptr,
4128 "update_values"));
4129 const std::size_t nqp = this->n_quadrature_points;
4131
4132 if (!is_face &&
4134 {
4135 // Cartesian cell
4136 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
4137 const VectorizedArrayType inv_det =
4138 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4139 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4140 this->jacobian[0][2][2];
4141
4142 // J * u * det(J^-1)
4143 for (unsigned int comp = 0; comp < n_components; ++comp)
4144 value_out[comp] = this->values_quad[comp * nqp + q_point] *
4145 jac[comp][comp] * inv_det;
4146 }
4147 else
4148 {
4149 // Affine or general cell
4150 const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
4151 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
4152 this->jacobian[q_point] :
4153 this->jacobian[0];
4155 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
4156 transpose(invert(inv_t_jac)) :
4157 this->jacobian[1];
4158
4159 // Derivatives are reordered for faces. Need to take this into
4160 // account
4161 const VectorizedArrayType inv_det =
4162 (is_face && dim == 2 && this->get_face_no() < 2) ?
4163 -determinant(inv_t_jac) :
4164 determinant(inv_t_jac);
4165 // J * u * det(J^-1)
4166 for (unsigned int comp = 0; comp < n_components; ++comp)
4167 {
4168 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
4169 for (unsigned int e = 1; e < dim; ++e)
4170 value_out[comp] +=
4171 this->values_quad[e * nqp + q_point] * jac[comp][e];
4172 value_out[comp] *= inv_det;
4173 }
4174 }
4175 return value_out;
4176 }
4177 else
4178 {
4179 const std::size_t nqp = this->n_quadrature_points;
4181 for (unsigned int comp = 0; comp < n_components; ++comp)
4182 return_value[comp] = this->values_quad[comp * nqp + q_point];
4183 return return_value;
4184 }
4185 }
4186}
4187
4188
4189
4190template <int dim,
4191 int n_components_,
4192 typename Number,
4193 bool is_face,
4194 typename VectorizedArrayType>
4196 typename FEEvaluationBase<dim,
4197 n_components_,
4198 Number,
4199 is_face,
4200 VectorizedArrayType>::gradient_type
4202 get_gradient(const unsigned int q_point) const
4203{
4204# ifdef DEBUG
4205 Assert(this->gradients_quad_initialized == true,
4207# endif
4208
4209 AssertIndexRange(q_point, this->n_quadrature_points);
4210 Assert(this->jacobian != nullptr,
4212 "update_gradients"));
4213 const std::size_t nqp = this->n_quadrature_points;
4214
4215 if constexpr (n_components == dim && dim > 1)
4216 {
4217 if (this->data->element_type ==
4219 {
4220 // Piola transform is required
4221# ifdef DEBUG
4222 Assert(this->gradients_quad_initialized == true,
4224# endif
4225
4226 AssertIndexRange(q_point, this->n_quadrature_points);
4227 Assert(this->jacobian != nullptr,
4229 "update_gradients"));
4230 const std::size_t nqp = this->n_quadrature_points;
4231 const std::size_t nqp_d = nqp * dim;
4233 const VectorizedArrayType *gradients =
4234 this->gradients_quad + q_point * dim;
4235
4236
4237 if (!is_face &&
4239 {
4240 // Cartesian cell
4241 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
4242 this->jacobian[0];
4244 this->jacobian[1];
4245 const VectorizedArrayType inv_det =
4246 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4247 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4248 this->jacobian[0][2][2];
4249
4250 // J * grad_quad * J^-1 * det(J^-1)
4251 for (unsigned int d = 0; d < dim; ++d)
4252 for (unsigned int comp = 0; comp < n_components; ++comp)
4253 grad_out[comp][d] = gradients[comp * nqp_d + d] *
4254 inv_t_jac[d][d] *
4255 (jac[comp][comp] * inv_det);
4256 }
4257 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4258 {
4259 // Affine cell
4260 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
4261 this->jacobian[0];
4263 this->jacobian[1];
4264
4265 // Derivatives are reordered for faces. Need to take this into
4266 // account
4267 const VectorizedArrayType inv_det =
4268 (is_face && dim == 2 && this->get_face_no() < 2) ?
4269 -determinant(inv_t_jac) :
4270 determinant(inv_t_jac);
4271
4272 VectorizedArrayType tmp[dim][dim];
4273 // J * grad_quad * J^-1 * det(J^-1)
4274 for (unsigned int d = 0; d < dim; ++d)
4275 for (unsigned int e = 0; e < dim; ++e)
4276 {
4277 tmp[d][e] = inv_t_jac[d][0] * gradients[e * nqp_d + 0];
4278 for (unsigned int f = 1; f < dim; ++f)
4279 tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
4280 }
4281 for (unsigned int comp = 0; comp < n_components; ++comp)
4282 for (unsigned int d = 0; d < dim; ++d)
4283 {
4284 VectorizedArrayType res = jac[comp][0] * tmp[d][0];
4285 for (unsigned int f = 1; f < dim; ++f)
4286 res += jac[comp][f] * tmp[d][f];
4287
4288 grad_out[comp][d] = res * inv_det;
4289 }
4290 }
4291 else
4292 {
4293 // General cell
4294
4295 // This assert could be removed if we make sure that this is
4296 // updated even though update_hessians or update_jacobian_grads is
4297 // not passed, i.e make the necessary changes in
4298 // MatrixFreeFunctions::MappingInfoStorage::compute_update_flags
4299 Assert(this->jacobian_gradients_non_inverse != nullptr,
4301 "update_hessians"));
4302
4303 const auto jac_grad =
4304 this->jacobian_gradients_non_inverse[q_point];
4305 const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
4306 this->jacobian[q_point];
4307
4308 // Derivatives are reordered for faces. Need to take this into
4309 // account
4310 const VectorizedArrayType inv_det =
4311 (is_face && dim == 2 && this->get_face_no() < 2) ?
4312 -determinant(inv_t_jac) :
4313 determinant(inv_t_jac);
4315 invert(inv_t_jac);
4316
4317 // (J * grad_quad) * J^-1 * det(J^-1), part in braces
4318 VectorizedArrayType tmp[dim][dim];
4319 for (unsigned int d = 0; d < dim; ++d)
4320 for (unsigned int e = 0; e < dim; ++e)
4321 {
4322 tmp[e][d] = t_jac[0][d] * gradients[0 * nqp_d + e];
4323 for (unsigned int f = 1; f < dim; ++f)
4324 tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
4325 }
4326
4327 // Add (jac_grad * values) * J^{-1} * det(J^{-1}), combine terms
4328 // outside braces with gradient part from above
4329 for (unsigned int d = 0; d < dim; ++d)
4330 {
4331 for (unsigned int e = 0; e < dim; ++e)
4332 tmp[e][d] +=
4333 jac_grad[e][d] * this->values_quad[e * nqp + q_point];
4334 for (unsigned int f = 0, r = dim; f < dim; ++f)
4335 for (unsigned int k = f + 1; k < dim; ++k, ++r)
4336 {
4337 tmp[k][d] +=
4338 jac_grad[r][d] * this->values_quad[f * nqp + q_point];
4339 tmp[f][d] +=
4340 jac_grad[r][d] * this->values_quad[k * nqp + q_point];
4341 }
4342 }
4343
4344 // Apply J^{-1} appearing in both terms outside braces above
4345 for (unsigned int d = 0; d < dim; ++d)
4346 for (unsigned int e = 0; e < dim; ++e)
4347 {
4348 VectorizedArrayType res = tmp[0][d] * inv_t_jac[e][0];
4349 for (unsigned int f = 1; f < dim; ++f)
4350 res += tmp[f][d] * inv_t_jac[e][f];
4351 grad_out[d][e] = res;
4352 }
4353
4354 // Add -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1})),
4355 // which can be expressed as a rank-1 update tmp[d] * tmp4[e],
4356 // where tmp = J * values and tmp4 = (J^{-T} * jac_grad * J^{-1})
4357 VectorizedArrayType tmp3[dim], tmp4[dim];
4358 for (unsigned int d = 0; d < dim; ++d)
4359 {
4360 tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0];
4361 for (unsigned int e = 1; e < dim; ++e)
4362 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
4363 }
4364 for (unsigned int e = 0, k = dim; e < dim; ++e)
4365 for (unsigned int f = e + 1; f < dim; ++k, ++f)
4366 for (unsigned int d = 0; d < dim; ++d)
4367 {
4368 tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d];
4369 tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d];
4370 }
4371 for (unsigned int d = 0; d < dim; ++d)
4372 {
4373 tmp4[d] = tmp3[0] * inv_t_jac[d][0];
4374 for (unsigned int e = 1; e < dim; ++e)
4375 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
4376 }
4377
4378 VectorizedArrayType tmp2[dim];
4379 for (unsigned int d = 0; d < dim; ++d)
4380 {
4381 tmp2[d] = t_jac[0][d] * this->values_quad[q_point];
4382 for (unsigned e = 1; e < dim; ++e)
4383 tmp2[d] +=
4384 t_jac[e][d] * this->values_quad[e * nqp + q_point];
4385 }
4386
4387 for (unsigned int d = 0; d < dim; ++d)
4388 for (unsigned int e = 0; e < dim; ++e)
4389 {
4390 grad_out[d][e] -= tmp4[e] * tmp2[d];
4391
4392 // finally multiply by det(J^{-1}) necessary in all
4393 // contributions above
4394 grad_out[d][e] *= inv_det;
4395 }
4396 }
4397 return grad_out;
4398 }
4399 }
4401
4402 // Cartesian cell
4403 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4404 {
4405 for (unsigned int comp = 0; comp < n_components; ++comp)
4406 for (unsigned int d = 0; d < dim; ++d)
4407 grad_out[comp][d] =
4408 this->gradients_quad[(comp * nqp + q_point) * dim + d] *
4409 this->jacobian[0][d][d];
4410 }
4411 // cell with general/affine Jacobian
4412 else
4413 {
4415 this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
4416 q_point :
4417 0];
4418 for (unsigned int comp = 0; comp < n_components; ++comp)
4419 for (unsigned int d = 0; d < dim; ++d)
4420 {
4421 grad_out[comp][d] =
4422 jac[d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4423 for (unsigned int e = 1; e < dim; ++e)
4424 grad_out[comp][d] +=
4425 jac[d][e] *
4426 this->gradients_quad[(comp * nqp + q_point) * dim + e];
4427 }
4428 }
4429 if constexpr (n_components == 1)
4430 return grad_out[0];
4431 else
4432 return grad_out;
4433}
4434
4435
4436
4437template <int dim,
4438 int n_components_,
4439 typename Number,
4440 bool is_face,
4441 typename VectorizedArrayType>
4443 typename FEEvaluationBase<dim,
4444 n_components_,
4445 Number,
4446 is_face,
4447 VectorizedArrayType>::value_type
4449 get_normal_derivative(const unsigned int q_point) const
4450{
4451 AssertIndexRange(q_point, this->n_quadrature_points);
4452# ifdef DEBUG
4453 Assert(this->gradients_quad_initialized == true,
4455# endif
4456
4457 Assert(this->normal_x_jacobian != nullptr,
4459 "update_gradients"));
4460
4461 const std::size_t nqp = this->n_quadrature_points;
4463
4464 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4465 for (unsigned int comp = 0; comp < n_components; ++comp)
4466 grad_out[comp] =
4467 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4468 (this->normal_x_jacobian[0][dim - 1]);
4469 else
4470 {
4471 const std::size_t index =
4472 this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
4473 for (unsigned int comp = 0; comp < n_components; ++comp)
4474 {
4475 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4476 this->normal_x_jacobian[index][0];
4477 for (unsigned int d = 1; d < dim; ++d)
4478 grad_out[comp] +=
4479 this->gradients_quad[(comp * nqp + q_point) * dim + d] *
4480 this->normal_x_jacobian[index][d];
4481 }
4482 }
4483 if constexpr (n_components == 1)
4484 return grad_out[0];
4485 else
4486 return grad_out;
4487}
4488
4489
4490
4491namespace internal
4492{
4493 // compute tmp = hess_unit(u) * J^T. do this manually because we do not
4494 // store the lower diagonal because of symmetry
4495 template <typename VectorizedArrayType>
4496 inline void
4497 hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
4498 const VectorizedArrayType *const hessians,
4499 const unsigned int,
4500 VectorizedArrayType (&tmp)[1][1])
4501 {
4502 tmp[0][0] = jac[0][0] * hessians[0];
4503 }
4504
4505 template <typename VectorizedArrayType>
4506 inline void
4507 hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
4508 const VectorizedArrayType *const hessians,
4509 const unsigned int nqp,
4510 VectorizedArrayType (&tmp)[2][2])
4511 {
4512 for (unsigned int d = 0; d < 2; ++d)
4513 {
4514 tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
4515 tmp[1][d] =
4516 (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
4517 }
4518 }
4519
4520 template <typename VectorizedArrayType>
4521 inline void
4522 hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
4523 const VectorizedArrayType *const hessians,
4524 const unsigned int nqp,
4525 VectorizedArrayType (&tmp)[3][3])
4526 {
4527 for (unsigned int d = 0; d < 3; ++d)
4528 {
4529 tmp[0][d] =
4530 (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
4531 jac[d][2] * hessians[4 * nqp]);
4532 tmp[1][d] =
4533 (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
4534 jac[d][2] * hessians[5 * nqp]);
4535 tmp[2][d] =
4536 (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
4537 jac[d][2] * hessians[2 * nqp]);
4538 }
4539 }
4540} // namespace internal
4541
4542
4543
4544template <int dim,
4545 int n_components_,
4546 typename Number,
4547 bool is_face,
4548 typename VectorizedArrayType>
4549inline typename FEEvaluationBase<dim,
4550 n_components_,
4551 Number,
4552 is_face,
4553 VectorizedArrayType>::hessian_type
4555 get_hessian(const unsigned int q_point) const
4556{
4557# ifdef DEBUG
4558 Assert(this->hessians_quad_initialized == true,
4560# endif
4561 AssertIndexRange(q_point, this->n_quadrature_points);
4562
4563 Assert(this->jacobian != nullptr,
4565 "update_hessian"));
4567 this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4568 0 :
4569 q_point];
4570
4572
4573 const std::size_t nqp = this->n_quadrature_points;
4574 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4575
4576 // Cartesian cell
4577 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4578 {
4579 for (unsigned int comp = 0; comp < n_components; ++comp)
4580 {
4581 for (unsigned int d = 0; d < dim; ++d)
4582 hessian_out[comp][d][d] =
4583 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4584 (jac[d][d] * jac[d][d]);
4585 switch (dim)
4586 {
4587 case 1:
4588 break;
4589 case 2:
4590 hessian_out[comp][0][1] =
4591 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4592 (jac[0][0] * jac[1][1]);
4593 break;
4594 case 3:
4595 hessian_out[comp][0][1] =
4596 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4597 (jac[0][0] * jac[1][1]);
4598 hessian_out[comp][0][2] =
4599 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4600 (jac[0][0] * jac[2][2]);
4601 hessian_out[comp][1][2] =
4602 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4603 (jac[1][1] * jac[2][2]);
4604 break;
4605 default:
4607 }
4608 for (unsigned int d = 0; d < dim; ++d)
4609 for (unsigned int e = d + 1; e < dim; ++e)
4610 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4611 }
4612 }
4613 // cell with general Jacobian, but constant within the cell
4614 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4615 {
4616 for (unsigned int comp = 0; comp < n_components; ++comp)
4617 {
4618 VectorizedArrayType tmp[dim][dim];
4619 internal::hessian_unit_times_jac(
4620 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4621
4622 // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4623 for (unsigned int d = 0; d < dim; ++d)
4624 for (unsigned int e = d; e < dim; ++e)
4625 {
4626 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4627 for (unsigned int f = 1; f < dim; ++f)
4628 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4629 }
4630
4631 // no J' * grad(u) part here because the Jacobian is constant
4632 // throughout the cell and hence, its derivative is zero
4633
4634 // take symmetric part
4635 for (unsigned int d = 0; d < dim; ++d)
4636 for (unsigned int e = d + 1; e < dim; ++e)
4637 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4638 }
4639 }
4640 // cell with general Jacobian
4641 else
4642 {
4643 const auto &jac_grad = this->jacobian_gradients[q_point];
4644 for (unsigned int comp = 0; comp < n_components; ++comp)
4645 {
4646 VectorizedArrayType tmp[dim][dim];
4647 internal::hessian_unit_times_jac(
4648 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4649
4650 // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4651 for (unsigned int d = 0; d < dim; ++d)
4652 for (unsigned int e = d; e < dim; ++e)
4653 {
4654 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4655 for (unsigned int f = 1; f < dim; ++f)
4656 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4657 }
4658
4659 // add diagonal part of J' * grad(u)
4660 for (unsigned int d = 0; d < dim; ++d)
4661 for (unsigned int e = 0; e < dim; ++e)
4662 hessian_out[comp][d][d] +=
4663 jac_grad[d][e] *
4664 this->gradients_quad[(comp * nqp + q_point) * dim + e];
4665
4666 // add off-diagonal part of J' * grad(u)
4667 for (unsigned int d = 0, count = dim; d < dim; ++d)
4668 for (unsigned int e = d + 1; e < dim; ++e, ++count)
4669 for (unsigned int f = 0; f < dim; ++f)
4670 hessian_out[comp][d][e] +=
4671 jac_grad[count][f] *
4672 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4673
4674 // take symmetric part
4675 for (unsigned int d = 0; d < dim; ++d)
4676 for (unsigned int e = d + 1; e < dim; ++e)
4677 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4678 }
4679 }
4680 if constexpr (n_components == 1)
4681 return hessian_out[0];
4682 else
4683 return hessian_out;
4684}
4685
4686
4687
4688template <int dim,
4689 int n_components_,
4690 typename Number,
4691 bool is_face,
4692 typename VectorizedArrayType>
4693inline typename FEEvaluationBase<dim,
4694 n_components_,
4695 Number,
4696 is_face,
4697 VectorizedArrayType>::gradient_type
4699 get_hessian_diagonal(const unsigned int q_point) const
4700{
4701 Assert(!is_face, ExcNotImplemented());
4702# ifdef DEBUG
4703 Assert(this->hessians_quad_initialized == true,
4705# endif
4706 AssertIndexRange(q_point, this->n_quadrature_points);
4707
4708 Assert(this->jacobian != nullptr, ExcNotImplemented());
4710 this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4711 0 :
4712 q_point];
4713
4714 const std::size_t nqp = this->n_quadrature_points;
4715 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4717
4718 // Cartesian cell
4719 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4720 {
4721 for (unsigned int comp = 0; comp < n_components; ++comp)
4722 for (unsigned int d = 0; d < dim; ++d)
4723 hessian_out[comp][d] =
4724 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4725 (jac[d][d] * jac[d][d]);
4726 }
4727 // cell with general Jacobian, but constant within the cell
4728 else if (this->cell_type == internal::MatrixFreeFunctions::affine)
4729 {
4730 for (unsigned int comp = 0; comp < n_components; ++comp)
4731 {
4732 // compute laplacian before the gradient because it needs to access
4733 // unscaled gradient data
4734 VectorizedArrayType tmp[dim][dim];
4735 internal::hessian_unit_times_jac(
4736 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4737
4738 // compute only the trace part of hessian, J * tmp = J *
4739 // hess_unit(u) * J^T
4740 for (unsigned int d = 0; d < dim; ++d)
4741 {
4742 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4743 for (unsigned int f = 1; f < dim; ++f)
4744 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4745 }
4746 }
4747 }
4748 // cell with general Jacobian
4749 else
4750 {
4751 const auto &jac_grad = this->jacobian_gradients[q_point];
4752 for (unsigned int comp = 0; comp < n_components; ++comp)
4753 {
4754 // compute laplacian before the gradient because it needs to access
4755 // unscaled gradient data
4756 VectorizedArrayType tmp[dim][dim];
4757 internal::hessian_unit_times_jac(
4758 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4759
4760 // compute only the trace part of hessian, J * tmp = J *
4761 // hess_unit(u) * J^T
4762 for (unsigned int d = 0; d < dim; ++d)
4763 {
4764 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4765 for (unsigned int f = 1; f < dim; ++f)
4766 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4767 }
4768
4769 for (unsigned int d = 0; d < dim; ++d)
4770 for (unsigned int e = 0; e < dim; ++e)
4771 hessian_out[comp][d] +=
4772 jac_grad[d][e] *
4773 this->gradients_quad[(comp * nqp + q_point) * dim + e];
4774 }
4775 }
4776
4777 if constexpr (n_components == 1)
4778 return hessian_out[0];
4779 else
4780 return hessian_out;
4781}
4782
4783
4784
4785template <int dim,
4786 int n_components_,
4787 typename Number,
4788 bool is_face,
4789 typename VectorizedArrayType>
4790inline typename FEEvaluationBase<dim,
4791 n_components_,
4792 Number,
4793 is_face,
4794 VectorizedArrayType>::value_type
4796 get_laplacian(const unsigned int q_point) const
4797{
4798 Assert(is_face == false, ExcNotImplemented());
4799# ifdef DEBUG
4800 Assert(this->hessians_quad_initialized == true,
4802# endif
4803 AssertIndexRange(q_point, this->n_quadrature_points);
4804
4805 const gradient_type hess_diag = get_hessian_diagonal(q_point);
4806 if constexpr (n_components == 1)
4807 {
4808 VectorizedArrayType sum = hess_diag[0];
4809 for (unsigned int d = 1; d < dim; ++d)
4810 sum += hess_diag[d];
4811 return sum;
4812 }
4813 else
4814 {
4816 for (unsigned int comp = 0; comp < n_components; ++comp)
4817 {
4818 laplacian_out[comp] = hess_diag[comp][0];
4819 for (unsigned int d = 1; d < dim; ++d)
4820 laplacian_out[comp] += hess_diag[comp][d];
4821 }
4822 return laplacian_out;
4823 }
4824}
4825
4826
4827
4828template <int dim,
4829 int n_components_,
4830 typename Number,
4831 bool is_face,
4832 typename VectorizedArrayType>
4833inline typename FEEvaluationBase<dim,
4834 n_components_,
4835 Number,
4836 is_face,
4837 VectorizedArrayType>::value_type
4839 get_normal_hessian(const unsigned int q_point) const
4840{
4841# ifdef DEBUG
4842 Assert(this->hessians_quad_initialized == true,
4844# endif
4845 AssertIndexRange(q_point, this->n_quadrature_points);
4846
4847 Assert(this->normal_x_jacobian != nullptr,
4849 "update_hessians"));
4850
4852
4853 const std::size_t nqp = this->n_quadrature_points;
4854 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4855
4856 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4857 {
4858 const auto nxj = this->normal_x_jacobian[0];
4859
4860 for (unsigned int comp = 0; comp < n_components; ++comp)
4861 {
4862 for (unsigned int d = 0; d < dim; ++d)
4863 hessian_out[comp] +=
4864 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4865 (nxj[d]) * (nxj[d]);
4866
4867 switch (dim)
4868 {
4869 case 1:
4870 break;
4871 case 2:
4872 hessian_out[comp] +=
4873 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4874 (nxj[0] * nxj[1]);
4875 break;
4876 case 3:
4877 hessian_out[comp] +=
4878 2. * this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4879 (nxj[0] * nxj[1]);
4880 hessian_out[comp] +=
4881 2. * this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4882 (nxj[0] * nxj[2]);
4883 hessian_out[comp] +=
4884 2. * this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4885 (nxj[1] * nxj[2]);
4886 break;
4887 default:
4889 }
4890 }
4891 }
4892 // cell with general Jacobian
4893 else
4894 {
4895 const auto normal = this->normal_vector(q_point);
4896 const auto hessian = get_hessian(q_point);
4897
4898 if constexpr (n_components == 1)
4899 hessian_out[0] = hessian * normal * normal;
4900 else
4901 for (unsigned int comp = 0; comp < n_components; ++comp)
4902 hessian_out[comp] = hessian[comp] * normal * normal;
4903 }
4904 if constexpr (n_components == 1)
4905 return hessian_out[0];
4906 else
4907 return hessian_out;
4908}
4909
4910
4911
4912template <int dim,
4913 int n_components_,
4914 typename Number,
4915 bool is_face,
4916 typename VectorizedArrayType>
4917inline DEAL_II_ALWAYS_INLINE void
4919 submit_dof_value(const value_type val_in, const unsigned int dof)
4920{
4921# ifdef DEBUG
4922 this->dof_values_initialized = true;
4923# endif
4924 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4925 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4926 for (unsigned int comp = 0; comp < n_components; ++comp)
4927 if constexpr (n_components == 1)
4928 this->values_dofs[comp * dofs + dof] = val_in;
4929 else
4930 this->values_dofs[comp * dofs + dof] = val_in[comp];
4931}
4932
4933
4934
4935template <int dim,
4936 int n_components_,
4937 typename Number,
4938 bool is_face,
4939 typename VectorizedArrayType>
4940inline DEAL_II_ALWAYS_INLINE void
4942 submit_value(const value_type val_in, const unsigned int q_point)
4943{
4944# ifdef DEBUG
4945 Assert(this->is_reinitialized, ExcNotInitialized());
4946# endif
4947 AssertIndexRange(q_point, this->n_quadrature_points);
4948 Assert(this->J_value != nullptr,
4950 "update_values"));
4951# ifdef DEBUG
4952 this->values_quad_submitted = true;
4953# endif
4954
4955 const std::size_t nqp = this->n_quadrature_points;
4956 VectorizedArrayType *values = this->values_quad + q_point;
4957
4958 const VectorizedArrayType JxW =
4959 this->cell_type <= internal::MatrixFreeFunctions::affine ?
4960 this->J_value[0] * this->quadrature_weights[q_point] :
4961 this->J_value[q_point];
4962 if constexpr (n_components == 1)
4963 values[0] = val_in * JxW;
4964 else
4965 {
4966 if (n_components == dim &&
4967 this->data->element_type ==
4969 {
4970 // Piola transform is required
4971 AssertIndexRange(q_point, this->n_quadrature_points);
4972 Assert(this->J_value != nullptr,
4974 "update_value"));
4975# ifdef DEBUG
4976 Assert(this->is_reinitialized, ExcNotInitialized());
4977 this->values_quad_submitted = true;
4978# endif
4979
4980 VectorizedArrayType *values = this->values_quad + q_point;
4981 const std::size_t nqp = this->n_quadrature_points;
4982
4983 if (!is_face &&
4985 {
4986 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
4987 const VectorizedArrayType weight =
4988 this->quadrature_weights[q_point];
4989
4990 for (unsigned int comp = 0; comp < n_components; ++comp)
4991 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
4992 }
4993 else
4994 {
4995 // Affine or general cell
4996 const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
4997 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
4998 this->jacobian[q_point] :
4999 this->jacobian[0];
5000
5001 // Derivatives are reordered for faces. Need to take this into
5002 // account and 1/inv_det != J_value for faces
5003 const VectorizedArrayType fac =
5004 (!is_face) ?
5005 this->quadrature_weights[q_point] :
5006 (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
5007 this->J_value[q_point] :
5008 this->J_value[0] * this->quadrature_weights[q_point]) *
5009 ((dim == 2 && this->get_face_no() < 2) ?
5010 -determinant(inv_t_jac) :
5011 determinant(inv_t_jac)));
5013 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5014 transpose(invert(inv_t_jac)) :
5015 this->jacobian[1];
5016
5017 // J^T * u * factor
5018 for (unsigned int comp = 0; comp < n_components; ++comp)
5019 {
5020 values[comp * nqp] = val_in[0] * jac[0][comp];
5021 for (unsigned int e = 1; e < dim; ++e)
5022 values[comp * nqp] += val_in[e] * jac[e][comp];
5023 values[comp * nqp] *= fac;
5024 }
5025 }
5026 }
5027 else
5028 for (unsigned int comp = 0; comp < n_components; ++comp)
5029 values[comp * nqp] = val_in[comp] * JxW;
5030 }
5031}
5032
5033
5034
5035template <int dim,
5036 int n_components_,
5037 typename Number,
5038 bool is_face,
5039 typename VectorizedArrayType>
5040template <int, typename>
5041inline DEAL_II_ALWAYS_INLINE void
5044 const unsigned int q_point)
5045{
5046 static_assert(n_components == 1,
5047 "Do not try to modify the default template parameters used for"
5048 " selectively enabling this function via std::enable_if!");
5049 submit_value(val_in[0], q_point);
5050}
5051
5052
5053
5054template <int dim,
5055 int n_components_,
5056 typename Number,
5057 bool is_face,
5058 typename VectorizedArrayType>
5059inline DEAL_II_ALWAYS_INLINE void
5061 submit_gradient(const gradient_type grad_in, const unsigned int q_point)
5062{
5063# ifdef DEBUG
5064 Assert(this->is_reinitialized, ExcNotInitialized());
5065# endif
5066 AssertIndexRange(q_point, this->n_quadrature_points);
5067 Assert(this->J_value != nullptr,
5069 "update_gradients"));
5070 Assert(this->jacobian != nullptr,
5072 "update_gradients"));
5073# ifdef DEBUG
5074 this->gradients_quad_submitted = true;
5075# endif
5076
5077 if constexpr (dim > 1 && n_components == dim)
5078 {
5079 if (this->data->element_type ==
5081 {
5082 // Piola transform is required
5083
5084# ifdef DEBUG
5085 Assert(this->is_reinitialized, ExcNotInitialized());
5086# endif
5087 AssertIndexRange(q_point, this->n_quadrature_points);
5088 Assert(this->J_value != nullptr,
5090 "update_gradients"));
5091 Assert(this->jacobian != nullptr,
5093 "update_gradients"));
5094# ifdef DEBUG
5095 this->gradients_quad_submitted = true;
5096# endif
5097
5098 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5099 VectorizedArrayType *values =
5100 this->values_from_gradients_quad + q_point;
5101 const std::size_t nqp = this->n_quadrature_points;
5102 const std::size_t nqp_d = nqp * dim;
5103
5104 if (!is_face &&
5106 {
5107 // Cartesian cell
5108 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5109 this->jacobian[0];
5111 this->jacobian[1];
5112 const VectorizedArrayType weight =
5113 this->quadrature_weights[q_point];
5114 for (unsigned int d = 0; d < dim; ++d)
5115 for (unsigned int comp = 0; comp < n_components; ++comp)
5116 gradients[comp * nqp_d + d] = grad_in[comp][d] *
5117 inv_t_jac[d][d] *
5118 (jac[comp][comp] * weight);
5119 }
5120 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5121 {
5122 // Affine cell
5123 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5124 this->jacobian[0];
5126 this->jacobian[1];
5127
5128 // Derivatives are reordered for faces. Need to take this into
5129 // account and 1/inv_det != J_value for faces
5130 const VectorizedArrayType fac =
5131 (!is_face) ?
5132 this->quadrature_weights[q_point] :
5133 this->J_value[0] * this->quadrature_weights[q_point] *
5134 ((dim == 2 && this->get_face_no() < 2) ?
5135 -determinant(inv_t_jac) :
5136 determinant(inv_t_jac));
5137
5138 // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
5139 VectorizedArrayType tmp[dim][dim];
5140 for (unsigned int d = 0; d < dim; ++d)
5141 for (unsigned int e = 0; e < dim; ++e)
5142 {
5143 tmp[d][e] = inv_t_jac[0][d] * grad_in[e][0];
5144 for (unsigned int f = 1; f < dim; ++f)
5145 tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
5146 }
5147 for (unsigned int comp = 0; comp < n_components; ++comp)
5148 for (unsigned int d = 0; d < dim; ++d)
5149 {
5150 VectorizedArrayType res = jac[0][comp] * tmp[d][0];
5151 for (unsigned int f = 1; f < dim; ++f)
5152 res += jac[f][comp] * tmp[d][f];
5153
5154 gradients[comp * nqp_d + d] = res * fac;
5155 }
5156 }
5157 else
5158 {
5159 // General cell
5160
5161 const auto jac_grad =
5162 this->jacobian_gradients_non_inverse[q_point];
5163 const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
5164 this->jacobian[q_point];
5165
5166 // Derivatives are reordered for faces. Need to take this into
5167 // account and 1/inv_det != J_value for faces
5168 const VectorizedArrayType fac =
5169 (!is_face) ? this->quadrature_weights[q_point] :
5170 this->J_value[q_point] *
5171 ((dim == 2 && this->get_face_no() < 2) ?
5172 -determinant(inv_t_jac) :
5173 determinant(inv_t_jac));
5175 invert(inv_t_jac);
5176
5177 // Start evaluation for values part below to enable the compiler
5178 // to possibly re-use the same computation in get_gradient()
5179 // without interfering with stores to 'gradients'
5180 VectorizedArrayType tmp3[dim], tmp4[dim];
5181 for (unsigned int d = 0; d < dim; ++d)
5182 {
5183 tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0];
5184 for (unsigned int e = 1; e < dim; ++e)
5185 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
5186 }
5187 for (unsigned int e = 0, k = dim; e < dim; ++e)
5188 for (unsigned int f = e + 1; f < dim; ++k, ++f)
5189 for (unsigned int d = 0; d < dim; ++d)
5190 {
5191 tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d];
5192 tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d];
5193 }
5194 for (unsigned int d = 0; d < dim; ++d)
5195 {
5196 tmp4[d] = tmp3[0] * inv_t_jac[d][0];
5197 for (unsigned int e = 1; e < dim; ++e)
5198 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
5199 }
5200
5201 const Tensor<2, dim, VectorizedArrayType> grad_in_scaled =
5202 fac * grad_in;
5203
5204 VectorizedArrayType tmp[dim][dim];
5205
5206 // J * (J^{-1} * (grad_in * factor))
5207 for (unsigned int d = 0; d < dim; ++d)
5208 for (unsigned int e = 0; e < dim; ++e)
5209 {
5210 tmp[d][e] = inv_t_jac[0][d] * grad_in_scaled[e][0];
5211 for (unsigned int f = 1; f < dim; ++f)
5212 tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
5213 }
5214
5215 for (unsigned int d = 0; d < dim; ++d)
5216 for (unsigned int e = 0; e < dim; ++e)
5217 {
5218 VectorizedArrayType res = t_jac[d][0] * tmp[e][0];
5219 for (unsigned int f = 1; f < dim; ++f)
5220 res += t_jac[d][f] * tmp[e][f];
5221
5222 gradients[d * nqp_d + e] = res;
5223 }
5224
5225 // jac_grad * (J^{-1} * (grad_in * factor)), re-use part in braces
5226 // as 'tmp' from above
5227 VectorizedArrayType value[dim];
5228 for (unsigned int d = 0; d < dim; ++d)
5229 {
5230 value[d] = tmp[d][0] * jac_grad[d][0];
5231 for (unsigned int e = 1; e < dim; ++e)
5232 value[d] += tmp[d][e] * jac_grad[d][e];
5233 }
5234 for (unsigned int e = 0, k = dim; e < dim; ++e)
5235 for (unsigned int f = e + 1; f < dim; ++k, ++f)
5236 for (unsigned int d = 0; d < dim; ++d)
5237 {
5238 value[e] += tmp[f][d] * jac_grad[k][d];
5239 value[f] += tmp[e][d] * jac_grad[k][d];
5240 }
5241
5242 // -(grad_in * factor) * J * (J^{-T} * jac_grad * J^{-1})
5243 // = -(grad_in * factor) * J * ( \------- tmp4 ---------/ )
5244 for (unsigned int d = 0; d < dim; ++d)
5245 {
5246 VectorizedArrayType tmp2 = grad_in_scaled[d][0] * tmp4[0];
5247 for (unsigned int e = 1; e < dim; ++e)
5248 tmp2 += grad_in_scaled[d][e] * tmp4[e];
5249 for (unsigned int e = 0; e < dim; ++e)
5250 value[e] -= t_jac[e][d] * tmp2;
5251 }
5252
5253 for (unsigned int d = 0; d < dim; ++d)
5254 values[d * nqp] = value[d];
5255 }
5256 return;
5257 }
5258 }
5259
5260 const std::size_t nqp_d = this->n_quadrature_points * dim;
5261 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5262
5263 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5264 {
5265 const VectorizedArrayType JxW =
5266 this->J_value[0] * this->quadrature_weights[q_point];
5267
5268 // Load all entries before starting to write back to make sure the
5269 // compiler sees opportunity of loads in a possibly nearby
5270 // get_gradient() function (i.e., the compiler should not think that
5271 // 'jacobian' could alias with 'gradients').
5272 std::array<VectorizedArrayType, dim> jac;
5273 for (unsigned int d = 0; d < dim; ++d)
5274 jac[d] = this->jacobian[0][d][d];
5275
5276 for (unsigned int d = 0; d < dim; ++d)
5277 {
5278 const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW;
5279 if constexpr (n_components == 1)
5280 gradients[d] = grad_in[d] * factor;
5281 else
5282 for (unsigned int comp = 0; comp < n_components; ++comp)
5283 gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5284 }
5285 }
5286 else
5287 {
5289 this->cell_type > internal::MatrixFreeFunctions::affine ?
5290 this->jacobian[q_point] :
5291 this->jacobian[0];
5292 const VectorizedArrayType JxW =
5293 this->cell_type > internal::MatrixFreeFunctions::affine ?
5294 this->J_value[q_point] :
5295 this->J_value[0] * this->quadrature_weights[q_point];
5296 if constexpr (n_components == 1)
5297 for (unsigned int d = 0; d < dim; ++d)
5298 {
5299 VectorizedArrayType new_val = jac[0][d] * grad_in[0];
5300 for (unsigned int e = 1; e < dim; ++e)
5301 new_val += (jac[e][d] * grad_in[e]);
5302 gradients[d] = new_val * JxW;
5303 }
5304 else
5305 for (unsigned int comp = 0; comp < n_components; ++comp)
5306 for (unsigned int d = 0; d < dim; ++d)
5307 {
5308 VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
5309 for (unsigned int e = 1; e < dim; ++e)
5310 new_val += (jac[e][d] * grad_in[comp][e]);
5311 gradients[comp * nqp_d + d] = new_val * JxW;
5312 }
5313 }
5314}
5315
5316
5317
5318template <int dim,
5319 int n_components_,
5320 typename Number,
5321 bool is_face,
5322 typename VectorizedArrayType>
5323template <int, typename>
5324inline DEAL_II_ALWAYS_INLINE void
5327 const unsigned int q_point)
5328{
5329 static_assert(n_components == 1 && dim == 1,
5330 "Do not try to modify the default template parameters used for"
5331 " selectively enabling this function via std::enable_if!");
5332 submit_gradient(grad_in[0], q_point);
5333}
5334
5335
5336
5337template <int dim,
5338 int n_components_,
5339 typename Number,
5340 bool is_face,
5341 typename VectorizedArrayType>
5342inline DEAL_II_ALWAYS_INLINE void
5344 submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
5345{
5346 AssertIndexRange(q_point, this->n_quadrature_points);
5347 Assert(this->normal_x_jacobian != nullptr,
5349 "update_gradients"));
5350# ifdef DEBUG
5351 this->gradients_quad_submitted = true;
5352# endif
5353
5354 const std::size_t nqp_d = this->n_quadrature_points * dim;
5355 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5356
5357 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5358 {
5359 const VectorizedArrayType JxW_jac = this->J_value[0] *
5360 this->quadrature_weights[q_point] *
5361 this->normal_x_jacobian[0][dim - 1];
5362 for (unsigned int comp = 0; comp < n_components; ++comp)
5363 {
5364 for (unsigned int d = 0; d < dim - 1; ++d)
5365 gradients[comp * nqp_d + d] = VectorizedArrayType();
5366 if constexpr (n_components == 1)
5367 gradients[dim - 1] = grad_in * JxW_jac;
5368 else
5369 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5370 }
5371 }
5372 else
5373 {
5374 const unsigned int index =
5375 this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5377 this->normal_x_jacobian[index];
5378 const VectorizedArrayType JxW =
5379 (this->cell_type <= internal::MatrixFreeFunctions::affine) ?
5380 this->J_value[index] * this->quadrature_weights[q_point] :
5381 this->J_value[index];
5382 for (unsigned int comp = 0; comp < n_components; ++comp)
5383 for (unsigned int d = 0; d < dim; ++d)
5384 if constexpr (n_components == 1)
5385 gradients[d] = (grad_in * JxW) * jac[d];
5386 else
5387 gradients[comp * nqp_d + d] = (grad_in[comp] * JxW) * jac[d];
5388 }
5389}
5390
5391
5392
5393template <int dim,
5394 int n_components_,
5395 typename Number,
5396 bool is_face,
5397 typename VectorizedArrayType>
5398inline DEAL_II_ALWAYS_INLINE void
5400 submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
5401{
5402# ifdef DEBUG
5403 Assert(this->is_reinitialized, ExcNotInitialized());
5404# endif
5405 AssertIndexRange(q_point, this->n_quadrature_points);
5406 Assert(this->J_value != nullptr,
5408 "update_hessians"));
5409 Assert(this->jacobian != nullptr,
5411 "update_hessians"));
5412# ifdef DEBUG
5413 this->hessians_quad_submitted = true;
5414# endif
5415
5416 // compute hessian_unit = J^T * hessian_in(u) * J
5417 const std::size_t nqp = this->n_quadrature_points;
5418 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5419 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5420 {
5421 const VectorizedArrayType JxW =
5422 this->J_value[0] * this->quadrature_weights[q_point];
5423
5424 // diagonal part
5425 for (unsigned int d = 0; d < dim; ++d)
5426 {
5427 const auto jac_d = this->jacobian[0][d][d];
5428 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5429 for (unsigned int comp = 0; comp < n_components; ++comp)
5430 if constexpr (n_components == 1)
5431 this->hessians_quad[d * nqp + q_point] =
5432 hessian_in[d][d] * factor;
5433 else
5434 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5435 hessian_in[comp][d][d] * factor;
5436 }
5437
5438 // off diagonal part
5439 for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5440 for (unsigned int e = 0; e < d; ++e, ++off_dia)
5441 {
5442 const auto jac_d = this->jacobian[0][d][d];
5443 const auto jac_e = this->jacobian[0][e][e];
5444 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5445 for (unsigned int comp = 0; comp < n_components; ++comp)
5446 if constexpr (n_components == 1)
5447 this->hessians_quad[off_dia * nqp + q_point] =
5448 (hessian_in[d][e] + hessian_in[e][d]) * factor;
5449 else
5450 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5451 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5452 }
5453 }
5454 // cell with general Jacobian, but constant within the cell
5455 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5456 {
5457 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
5458 const VectorizedArrayType JxW =
5459 this->J_value[0] * this->quadrature_weights[q_point];
5460 for (unsigned int comp = 0; comp < n_components; ++comp)
5461 {
5463 if constexpr (n_components == 1)
5464 hessian_c = hessian_in;
5465 else
5466 hessian_c = hessian_in[comp];
5467
5468 // 1. tmp = hessian(u) * J
5469 VectorizedArrayType tmp[dim][dim];
5470 for (unsigned int i = 0; i < dim; ++i)
5471 for (unsigned int j = 0; j < dim; ++j)
5472 {
5473 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5474 for (unsigned int k = 1; k < dim; ++k)
5475 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5476 }
5477
5478 // 2. hessian_unit = J^T * tmp
5479 VectorizedArrayType tmp2[dim][dim];
5480 for (unsigned int i = 0; i < dim; ++i)
5481 for (unsigned int j = 0; j < dim; ++j)
5482 {
5483 tmp2[i][j] = jac[0][i] * tmp[0][j];
5484 for (unsigned int k = 1; k < dim; ++k)
5485 tmp2[i][j] += jac[k][i] * tmp[k][j];
5486 }
5487
5488 // diagonal part
5489 for (unsigned int d = 0; d < dim; ++d)
5490 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5491 tmp2[d][d] * JxW;
5492
5493 // off diagonal part
5494 for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5495 for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5496 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5497 (tmp2[d][e] + tmp2[e][d]) * JxW;
5498 }
5499 }
5500 else
5501 {
5502 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
5503 const VectorizedArrayType JxW = this->J_value[q_point];
5504 const auto &jac_grad = this->jacobian_gradients[q_point];
5505 for (unsigned int comp = 0; comp < n_components; ++comp)
5506 {
5508 if constexpr (n_components == 1)
5509 hessian_c = hessian_in;
5510 else
5511 hessian_c = hessian_in[comp];
5512
5513 // 1. tmp = hessian(u) * J
5514 VectorizedArrayType tmp[dim][dim];
5515 for (unsigned int i = 0; i < dim; ++i)
5516 for (unsigned int j = 0; j < dim; ++j)
5517 {
5518 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5519 for (unsigned int k = 1; k < dim; ++k)
5520 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5521 }
5522
5523 // 2. hessian_unit = J^T * tmp
5524 VectorizedArrayType tmp2[dim][dim];
5525 for (unsigned int i = 0; i < dim; ++i)
5526 for (unsigned int j = 0; j < dim; ++j)
5527 {
5528 tmp2[i][j] = jac[0][i] * tmp[0][j];
5529 for (unsigned int k = 1; k < dim; ++k)
5530 tmp2[i][j] += jac[k][i] * tmp[k][j];
5531 }
5532
5533 // diagonal part
5534 for (unsigned int d = 0; d < dim; ++d)
5535 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5536 tmp2[d][d] * JxW;
5537
5538 // off diagonal part
5539 for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5540 for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5541 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5542 (tmp2[d][e] + tmp2[e][d]) * JxW;
5543
5544 // 3. gradient_unit = J' * hessian
5545 for (unsigned int d = 0; d < dim; ++d)
5546 {
5547 VectorizedArrayType sum = 0;
5548 for (unsigned int e = 0; e < dim; ++e)
5549 sum += hessian_c[e][e] * jac_grad[e][d];
5550 for (unsigned int e = 0, count = dim; e < dim; ++e)
5551 for (unsigned int f = e + 1; f < dim; ++f, ++count)
5552 sum +=
5553 (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][d];
5554 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5555 d] = sum * JxW;
5556 }
5557 }
5558 }
5559}
5560
5561
5562
5563template <int dim,
5564 int n_components_,
5565 typename Number,
5566 bool is_face,
5567 typename VectorizedArrayType>
5568inline DEAL_II_ALWAYS_INLINE void
5570 submit_normal_hessian(const value_type normal_hessian_in,
5571 const unsigned int q_point)
5572{
5573# ifdef DEBUG
5574 Assert(this->is_reinitialized, ExcNotInitialized());
5575# endif
5576 AssertIndexRange(q_point, this->n_quadrature_points);
5577 Assert(this->J_value != nullptr,
5579 "update_hessians"));
5580 Assert(this->jacobian != nullptr,
5582 "update_hessians"));
5583# ifdef DEBUG
5584 this->hessians_quad_submitted = true;
5585# endif
5586
5587 // compute hessian_unit = J^T * hessian_in(u) * J
5588 const std::size_t nqp = this->n_quadrature_points;
5589 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5590 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5591 {
5592 const VectorizedArrayType JxW =
5593 this->J_value[0] * this->quadrature_weights[q_point];
5594
5595 const auto nxj = this->normal_x_jacobian[0];
5596
5597 // diagonal part
5598 for (unsigned int d = 0; d < dim; ++d)
5599 {
5600 const auto nxj_d = nxj[d];
5601 const VectorizedArrayType factor = nxj_d * nxj_d * JxW;
5602 for (unsigned int comp = 0; comp < n_components; ++comp)
5603 if constexpr (n_components == 1)
5604 this->hessians_quad[d * nqp + q_point] =
5605 normal_hessian_in * factor;
5606 else
5607 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5608 normal_hessian_in[comp] * factor;
5609 }
5610
5611 // off diagonal part
5612 for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5613 for (unsigned int e = 0; e < d; ++e, ++off_dia)
5614 {
5615 const auto jac_d = nxj[d];
5616 const auto jac_e = nxj[e];
5617 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5618 for (unsigned int comp = 0; comp < n_components; ++comp)
5619 if constexpr (n_components == 1)
5620 this->hessians_quad[off_dia * nqp + q_point] =
5621 2. * normal_hessian_in * factor;
5622 else
5623 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5624 2. * normal_hessian_in[comp] * factor;
5625 }
5626 }
5627 else
5628 {
5629 const auto normal = this->normal_vector(q_point);
5630 const auto normal_projector = outer_product(normal, normal);
5631 if constexpr (n_components == 1)
5632 submit_hessian(normal_hessian_in * normal_projector, q_point);
5633 else
5634 {
5635 hessian_type tmp;
5636 for (unsigned int comp = 0; comp < n_components; ++comp)
5637 tmp[comp] = normal_hessian_in[comp] * normal_projector;
5638 submit_hessian(tmp, q_point);
5639 }
5640 }
5641}
5642
5643
5644
5645template <int dim,
5646 int n_components_,
5647 typename Number,
5648 bool is_face,
5649 typename VectorizedArrayType>
5650inline typename FEEvaluationBase<dim,
5651 n_components_,
5652 Number,
5653 is_face,
5654 VectorizedArrayType>::value_type
5656 integrate_value() const
5657{
5658# ifdef DEBUG
5659 Assert(this->is_reinitialized, ExcNotInitialized());
5660 Assert(this->values_quad_submitted == true,
5662# endif
5663
5665 const std::size_t nqp = this->n_quadrature_points;
5666 for (unsigned int q = 0; q < nqp; ++q)
5667 for (unsigned int comp = 0; comp < n_components; ++comp)
5668 return_value[comp] += this->values_quad[comp * nqp + q];
5669 if constexpr (n_components == 1)
5670 return return_value[0];
5671 else
5672 return return_value;
5673}
5674
5675
5676
5677template <int dim,
5678 int n_components_,
5679 typename Number,
5680 bool is_face,
5681 typename VectorizedArrayType>
5682template <int, typename>
5683inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5685 get_divergence(const unsigned int q_point) const
5686{
5687 static_assert(n_components == dim,
5688 "Do not try to modify the default template parameters used for"
5689 " selectively enabling this function via std::enable_if!");
5690
5691# ifdef DEBUG
5692 Assert(this->gradients_quad_initialized == true,
5694# endif
5695 AssertIndexRange(q_point, this->n_quadrature_points);
5696 Assert(this->jacobian != nullptr,
5698 "update_gradients"));
5699
5700 VectorizedArrayType divergence;
5701 const std::size_t nqp = this->n_quadrature_points;
5702
5703 if (dim > 1 &&
5704 this->data->element_type ==
5706 {
5707 VectorizedArrayType inv_det =
5708 (!is_face &&
5709 this->cell_type == internal::MatrixFreeFunctions::cartesian) ?
5710 this->jacobian[0][0][0] *
5711 ((dim == 2) ? this->jacobian[0][1][1] :
5712 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
5713 determinant(this->jacobian[this->cell_type >
5714 internal::MatrixFreeFunctions::affine ?
5715 q_point :
5716 0]);
5717
5718 // on faces in 2d, the determinant has the wrong sign due to ordering of
5719 // derivatives
5720 if (is_face && dim == 2 && this->get_face_no() < 2)
5721 inv_det = -inv_det;
5722
5723 // div * det(J^-1)
5724 divergence = this->gradients_quad[q_point * dim];
5725 for (unsigned int d = 1; d < dim; ++d)
5726 divergence += this->gradients_quad[(d * nqp + q_point) * dim + d];
5727 divergence *= inv_det;
5728 }
5729 else
5730 {
5731 if (!is_face &&
5733 {
5734 // Cartesian cell
5735 divergence =
5736 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
5737 for (unsigned int d = 1; d < dim; ++d)
5738 divergence += this->gradients_quad[(d * nqp + q_point) * dim + d] *
5739 this->jacobian[0][d][d];
5740 }
5741 else
5742 {
5743 // cell with general/constant Jacobian
5745 this->cell_type == internal::MatrixFreeFunctions::general ?
5746 this->jacobian[q_point] :
5747 this->jacobian[0];
5748 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
5749 for (unsigned int e = 1; e < dim; ++e)
5750 divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
5751 for (unsigned int d = 1; d < dim; ++d)
5752 for (unsigned int e = 0; e < dim; ++e)
5753 divergence +=
5754 jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim + e];
5755 }
5756 }
5757 return divergence;
5758}
5759
5760
5761
5762template <int dim,
5763 int n_components_,
5764 typename Number,
5765 bool is_face,
5766 typename VectorizedArrayType>
5767template <int, typename>
5770 get_symmetric_gradient(const unsigned int q_point) const
5771{
5772 static_assert(n_components == dim,
5773 "Do not try to modify the default template parameters used for"
5774 " selectively enabling this function via std::enable_if!");
5775
5776 // copy from generic function into dim-specialization function
5777 const auto grad = get_gradient(q_point);
5778 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5779 VectorizedArrayType half = Number(0.5);
5780 for (unsigned int d = 0; d < dim; ++d)
5781 symmetrized[d] = grad[d][d];
5782 switch (dim)
5783 {
5784 case 1:
5785 break;
5786 case 2:
5787 symmetrized[2] = grad[0][1] + grad[1][0];
5788 symmetrized[2] *= half;
5789 break;
5790 case 3:
5791 symmetrized[3] = grad[0][1] + grad[1][0];
5792 symmetrized[3] *= half;
5793 symmetrized[4] = grad[0][2] + grad[2][0];
5794 symmetrized[4] *= half;
5795 symmetrized[5] = grad[1][2] + grad[2][1];
5796 symmetrized[5] *= half;
5797 break;
5798 default:
5800 }
5802}
5803
5804
5805
5806template <int dim,
5807 int n_components_,
5808 typename Number,
5809 bool is_face,
5810 typename VectorizedArrayType>
5811template <int, typename>
5813 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
5815 get_curl(const unsigned int q_point) const
5816{
5817 static_assert(dim > 1 && n_components == dim,
5818 "Do not try to modify the default template parameters used for"
5819 " selectively enabling this function via std::enable_if!");
5820
5821 // copy from generic function into dim-specialization function
5822 const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
5823 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
5824 switch (dim)
5825 {
5826 case 2:
5827 curl[0] = grad[1][0] - grad[0][1];
5828 break;
5829 case 3:
5830 curl[0] = grad[2][1] - grad[1][2];
5831 curl[1] = grad[0][2] - grad[2][0];
5832 curl[2] = grad[1][0] - grad[0][1];
5833 break;
5834 default:
5836 }
5837 return curl;
5838}
5839
5840
5841
5842template <int dim,
5843 int n_components_,
5844 typename Number,
5845 bool is_face,
5846 typename VectorizedArrayType>
5847template <int, typename>
5848inline DEAL_II_ALWAYS_INLINE void
5850 submit_divergence(const VectorizedArrayType div_in,
5851 const unsigned int q_point)
5852{
5853 static_assert(n_components == dim,
5854 "Do not try to modify the default template parameters used for"
5855 " selectively enabling this function via std::enable_if!");
5856
5857# ifdef DEBUG
5858 Assert(this->is_reinitialized, ExcNotInitialized());
5859# endif
5860 AssertIndexRange(q_point, this->n_quadrature_points);
5861 Assert(this->J_value != nullptr,
5863 "update_gradients"));
5864 Assert(this->jacobian != nullptr,
5866 "update_gradients"));
5867# ifdef DEBUG
5868 this->gradients_quad_submitted = true;
5869# endif
5870
5871 const std::size_t nqp_d = this->n_quadrature_points * dim;
5872 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5873
5874 if (this->data->element_type ==
5876 {
5877 // General cell
5878
5879 // Derivatives are reordered for faces. Need to take this into account
5880 // and 1/inv_det != J_value for faces
5881 const VectorizedArrayType fac =
5882 (!is_face) ?
5883 this->quadrature_weights[q_point] * div_in :
5884 (this->cell_type > internal::MatrixFreeFunctions::affine ?
5885 this->J_value[q_point] :
5886 this->J_value[0] * this->quadrature_weights[q_point]) *
5887 div_in *
5889 this->jacobian[this->cell_type >
5890 internal::MatrixFreeFunctions::affine ?
5891 q_point :
5892 0]) *
5893 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
5894
5895 for (unsigned int d = 0; d < dim; ++d)
5896 {
5897 for (unsigned int e = 0; e < dim; ++e)
5898 gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
5899 }
5900 this->divergence_is_requested = true;
5901 }
5902 else
5903 {
5904 if (!is_face &&
5906 {
5907 const VectorizedArrayType fac =
5908 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
5909 for (unsigned int d = 0; d < dim; ++d)
5910 {
5911 const VectorizedArrayType jac_dd = this->jacobian[0][d][d];
5912 for (unsigned int e = 0; e < dim; ++e)
5913 gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
5914 }
5915 }
5916 else
5917 {
5919 this->cell_type == internal::MatrixFreeFunctions::general ?
5920 this->jacobian[q_point] :
5921 this->jacobian[0];
5922 const VectorizedArrayType fac =
5923 (this->cell_type == internal::MatrixFreeFunctions::general ?
5924 this->J_value[q_point] :
5925 this->J_value[0] * this->quadrature_weights[q_point]) *
5926 div_in;
5927 for (unsigned int d = 0; d < dim; ++d)
5928 {
5929 for (unsigned int e = 0; e < dim; ++e)
5930 gradients[d * nqp_d + e] = jac[d][e] * fac;
5931 }
5932 }
5933 }
5934}
5935
5936
5937
5938template <int dim,
5939 int n_components_,
5940 typename Number,
5941 bool is_face,
5942 typename VectorizedArrayType>
5943template <int, typename>
5944inline DEAL_II_ALWAYS_INLINE void
5948 const unsigned int q_point)
5949{
5950 static_assert(n_components == dim,
5951 "Do not try to modify the default template parameters used for"
5952 " selectively enabling this function via std::enable_if!");
5953
5955 this->data->element_type !=
5958
5959 // could have used base class operator, but that involves some overhead
5960 // which is inefficient. it is nice to have the symmetric tensor because
5961 // that saves some operations
5962# ifdef DEBUG
5963 Assert(this->is_reinitialized, ExcNotInitialized());
5964# endif
5965 AssertIndexRange(q_point, this->n_quadrature_points);
5966 Assert(this->J_value != nullptr,
5968 "update_gradients"));
5969 Assert(this->jacobian != nullptr,
5971 "update_gradients"));
5972# ifdef DEBUG
5973 this->gradients_quad_submitted = true;
5974# endif
5975
5976 const std::size_t nqp_d = this->n_quadrature_points * dim;
5977 VectorizedArrayType *gradients = this->gradients_quad + dim * q_point;
5978 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5979 {
5980 const VectorizedArrayType JxW =
5981 this->J_value[0] * this->quadrature_weights[q_point];
5982 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
5983 for (unsigned int d = 0; d < dim; ++d)
5984 gradients[d * nqp_d + d] =
5985 (sym_grad.access_raw_entry(d) * JxW * jac[d][d]);
5986 for (unsigned int e = 0, counter = dim; e < dim; ++e)
5987 for (unsigned int d = e + 1; d < dim; ++d, ++counter)
5988 {
5989 const VectorizedArrayType value =
5990 sym_grad.access_raw_entry(counter) * JxW;
5991 gradients[e * nqp_d + d] = value * jac[d][d];
5992 gradients[d * nqp_d + e] = value * jac[e][e];
5993 }
5994 }
5995 // general/affine cell type
5996 else
5997 {
5998 const VectorizedArrayType JxW =
5999 this->cell_type == internal::MatrixFreeFunctions::general ?
6000 this->J_value[q_point] :
6001 this->J_value[0] * this->quadrature_weights[q_point];
6003 this->cell_type == internal::MatrixFreeFunctions::general ?
6004 this->jacobian[q_point] :
6005 this->jacobian[0];
6006 VectorizedArrayType weighted[dim][dim];
6007 for (unsigned int i = 0; i < dim; ++i)
6008 weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6009 for (unsigned int i = 0, counter = dim; i < dim; ++i)
6010 for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6011 {
6012 const VectorizedArrayType value =
6013 sym_grad.access_raw_entry(counter) * JxW;
6014 weighted[i][j] = value;
6015 weighted[j][i] = value;
6016 }
6017 for (unsigned int comp = 0; comp < dim; ++comp)
6018 for (unsigned int d = 0; d < dim; ++d)
6019 {
6020 VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6021 for (unsigned int e = 1; e < dim; ++e)
6022 new_val += jac[e][d] * weighted[comp][e];
6023 gradients[comp * nqp_d + d] = new_val;
6024 }
6025 }
6026}
6027
6028
6029
6030template <int dim,
6031 int n_components_,
6032 typename Number,
6033 bool is_face,
6034 typename VectorizedArrayType>
6035template <int, typename>
6036inline DEAL_II_ALWAYS_INLINE void
6039 const unsigned int q_point)
6040{
6041 static_assert(n_components == dim,
6042 "Do not try to modify the default template parameters used for"
6043 " selectively enabling this function via std::enable_if!");
6044
6046 switch (dim)
6047 {
6048 case 2:
6049 grad[1][0] = curl[0];
6050 grad[0][1] = -curl[0];
6051 break;
6052 case 3:
6053 grad[2][1] = curl[0];
6054 grad[1][2] = -curl[0];
6055 grad[0][2] = curl[1];
6056 grad[2][0] = -curl[1];
6057 grad[1][0] = curl[2];
6058 grad[0][1] = -curl[2];
6059 break;
6060 default:
6062 }
6063 submit_gradient(grad, q_point);
6064}
6065
6066
6067
6068/*-------------------------- FEEvaluation -----------------------------------*/
6069
6070
6071template <int dim,
6072 int fe_degree,
6073 int n_q_points_1d,
6074 int n_components_,
6075 typename Number,
6076 typename VectorizedArrayType>
6077inline FEEvaluation<dim,
6078 fe_degree,
6079 n_q_points_1d,
6080 n_components_,
6081 Number,
6082 VectorizedArrayType>::
6083 FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
6084 const unsigned int fe_no,
6085 const unsigned int quad_no,
6086 const unsigned int first_selected_component,
6087 const unsigned int active_fe_index,
6088 const unsigned int active_quad_index)
6089 : BaseClass(matrix_free,
6090 fe_no,
6091 first_selected_component,
6092 quad_no,
6093 fe_degree,
6094 static_n_q_points,
6095 true /*note: this is not a face*/,
6096 active_fe_index,
6097 active_quad_index,
6098 numbers::invalid_unsigned_int /*face_type*/)
6099 , dofs_per_component(this->data->dofs_per_component_on_cell)
6100 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6101 , n_q_points(this->data->n_q_points)
6102{
6103 check_template_arguments(fe_no, 0);
6104}
6105
6106
6107
6108template <int dim,
6109 int fe_degree,
6110 int n_q_points_1d,
6111 int n_components_,
6112 typename Number,
6113 typename VectorizedArrayType>
6114inline FEEvaluation<dim,
6115 fe_degree,
6116 n_q_points_1d,
6117 n_components_,
6118 Number,
6119 VectorizedArrayType>::
6120 FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
6121 const std::pair<unsigned int, unsigned int> &range,
6122 const unsigned int dof_no,
6123 const unsigned int quad_no,
6124 const unsigned int first_selected_component)
6125 : FEEvaluation(matrix_free,
6126 dof_no,
6127 quad_no,
6128 first_selected_component,
6129 matrix_free.get_cell_active_fe_index(range))
6130{}
6131
6132
6133
6134template <int dim,
6135 int fe_degree,
6136 int n_q_points_1d,
6137 int n_components_,
6138 typename Number,
6139 typename VectorizedArrayType>
6140inline FEEvaluation<dim,
6141 fe_degree,
6142 n_q_points_1d,
6143 n_components_,
6144 Number,
6145 VectorizedArrayType>::
6146 FEEvaluation(const Mapping<dim> &mapping,
6147 const FiniteElement<dim> &fe,
6148 const Quadrature<1> &quadrature,
6149 const UpdateFlags update_flags,
6150 const unsigned int first_selected_component)
6151 : BaseClass(mapping,
6152 fe,
6153 quadrature,
6154 update_flags,
6155 first_selected_component,
6156 nullptr)
6157 , dofs_per_component(this->data->dofs_per_component_on_cell)
6158 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6159 , n_q_points(this->data->n_q_points)
6160{
6161 check_template_arguments(numbers::invalid_unsigned_int, 0);
6162}
6163
6164
6165
6166template <int dim,
6167 int fe_degree,
6168 int n_q_points_1d,
6169 int n_components_,
6170 typename Number,
6171 typename VectorizedArrayType>
6172inline FEEvaluation<dim,
6173 fe_degree,
6174 n_q_points_1d,
6175 n_components_,
6176 Number,
6177 VectorizedArrayType>::
6178 FEEvaluation(const FiniteElement<dim> &fe,
6179 const Quadrature<1> &quadrature,
6180 const UpdateFlags update_flags,
6181 const unsigned int first_selected_component)
6182 : BaseClass(StaticMappingQ1<dim>::mapping,
6183 fe,
6184 quadrature,
6185 update_flags,
6186 first_selected_component,
6187 nullptr)
6188 , dofs_per_component(this->data->dofs_per_component_on_cell)
6189 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6190 , n_q_points(this->data->n_q_points)
6191{
6192 check_template_arguments(numbers::invalid_unsigned_int, 0);
6193}
6194
6195
6196
6197template <int dim,
6198 int fe_degree,
6199 int n_q_points_1d,
6200 int n_components_,
6201 typename Number,
6202 typename VectorizedArrayType>
6203inline FEEvaluation<dim,
6204 fe_degree,
6205 n_q_points_1d,
6206 n_components_,
6207 Number,
6208 VectorizedArrayType>::
6209 FEEvaluation(const FiniteElement<dim> &fe,
6211 const unsigned int first_selected_component)
6212 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6213 fe,
6214 other.mapped_geometry->get_quadrature(),
6215 other.mapped_geometry->get_fe_values().get_update_flags(),
6216 first_selected_component,
6217 &other)
6218 , dofs_per_component(this->data->dofs_per_component_on_cell)
6219 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6220 , n_q_points(this->data->n_q_points)
6221{
6222 check_template_arguments(numbers::invalid_unsigned_int, 0);
6223}
6224
6225
6226
6227template <int dim,
6228 int fe_degree,
6229 int n_q_points_1d,
6230 int n_components_,
6231 typename Number,
6232 typename VectorizedArrayType>
6233inline FEEvaluation<dim,
6234 fe_degree,
6235 n_q_points_1d,
6236 n_components_,
6237 Number,
6238 VectorizedArrayType>::FEEvaluation(const FEEvaluation
6239 &other)
6240 : BaseClass(other)
6241 , dofs_per_component(this->data->dofs_per_component_on_cell)
6242 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6243 , n_q_points(this->data->n_q_points)
6244{
6245 check_template_arguments(numbers::invalid_unsigned_int, 0);
6246}
6247
6248
6249
6250template <int dim,
6251 int fe_degree,
6252 int n_q_points_1d,
6253 int n_components_,
6254 typename Number,
6255 typename VectorizedArrayType>
6256inline FEEvaluation<dim,
6257 fe_degree,
6258 n_q_points_1d,
6259 n_components_,
6260 Number,
6261 VectorizedArrayType> &
6262FEEvaluation<dim,
6263 fe_degree,
6264 n_q_points_1d,
6265 n_components_,
6266 Number,
6267 VectorizedArrayType>::operator=(const FEEvaluation &other)
6268{
6269 BaseClass::operator=(other);
6270 check_template_arguments(numbers::invalid_unsigned_int, 0);
6271 return *this;
6272}
6273
6274
6275
6276template <int dim,
6277 int fe_degree,
6278 int n_q_points_1d,
6279 int n_components_,
6280 typename Number,
6281 typename VectorizedArrayType>
6282inline void
6283FEEvaluation<dim,
6284 fe_degree,
6285 n_q_points_1d,
6286 n_components_,
6287 Number,
6288 VectorizedArrayType>::
6289 check_template_arguments(const unsigned int dof_no,
6290 const unsigned int first_selected_component)
6291{
6292 (void)dof_no;
6293 (void)first_selected_component;
6294
6295 Assert(
6296 this->data->dofs_per_component_on_cell > 0,
6297 ExcMessage(
6298 "There is nothing useful you can do with an FEEvaluation object with "
6299 "FE_Nothing, i.e., without DoFs! If you have passed to "
6300 "MatrixFree::reinit() a collection of finite elements also containing "
6301 "FE_Nothing, please check - before creating FEEvaluation - the category "
6302 "of the current range by calling either "
6303 "MatrixFree::get_cell_range_category(range) or "
6304 "MatrixFree::get_face_range_category(range). The returned category "
6305 "is the index of the active FE, which you can use to exclude "
6306 "FE_Nothing."));
6307
6308# ifdef DEBUG
6309 // print error message when the dimensions do not match. Propose a possible
6310 // fix
6311 if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
6312 static_cast<unsigned int>(fe_degree) !=
6313 this->data->data.front().fe_degree) ||
6314 n_q_points != this->n_quadrature_points)
6315 {
6316 std::string message =
6317 "-------------------------------------------------------\n";
6318 message += "Illegal arguments in constructor/wrong template arguments!\n";
6319 message += " Called --> FEEvaluation<dim,";
6320 message += Utilities::int_to_string(fe_degree) + ",";
6321 message += Utilities::int_to_string(n_q_points_1d);
6322 message += "," + Utilities::int_to_string(n_components);
6323 message += ",Number>(data";
6324 if (first_selected_component != numbers::invalid_unsigned_int)
6325 {
6326 message += ", " + Utilities::int_to_string(dof_no) + ", ";
6327 message += Utilities::int_to_string(this->quad_no) + ", ";
6328 message += Utilities::int_to_string(first_selected_component);
6329 }
6330 message += ")\n";
6331
6332 // check whether some other vector component has the correct number of
6333 // points
6334 unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
6335 proposed_fe_comp = numbers::invalid_unsigned_int,
6336 proposed_quad_comp = numbers::invalid_unsigned_int;
6337 if (dof_no != numbers::invalid_unsigned_int)
6338 {
6339 if (static_cast<unsigned int>(fe_degree) ==
6340 this->data->data.front().fe_degree)
6341 {
6342 proposed_dof_comp = dof_no;
6343 proposed_fe_comp = first_selected_component;
6344 }
6345 else
6346 for (unsigned int no = 0; no < this->matrix_free->n_components();
6347 ++no)
6348 for (unsigned int nf = 0;
6349 nf < this->matrix_free->n_base_elements(no);
6350 ++nf)
6351 if (this->matrix_free
6352 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6353 .data.front()
6354 .fe_degree == static_cast<unsigned int>(fe_degree))
6355 {
6356 proposed_dof_comp = no;
6357 proposed_fe_comp = nf;
6358 break;
6359 }
6360 if (n_q_points ==
6361 this->mapping_data->descriptor[this->active_quad_index]
6362 .n_q_points)
6363 proposed_quad_comp = this->quad_no;
6364 else
6365 for (unsigned int no = 0;
6366 no < this->matrix_free->get_mapping_info().cell_data.size();
6367 ++no)
6368 if (this->matrix_free->get_mapping_info()
6369 .cell_data[no]
6370 .descriptor[this->active_quad_index]
6371 .n_q_points == n_q_points)
6372 {
6373 proposed_quad_comp = no;
6374 break;
6375 }
6376 }
6377 if (proposed_dof_comp != numbers::invalid_unsigned_int &&
6378 proposed_quad_comp != numbers::invalid_unsigned_int)
6379 {
6380 if (proposed_dof_comp != first_selected_component)
6381 message += "Wrong vector component selection:\n";
6382 else
6383 message += "Wrong quadrature formula selection:\n";
6384 message += " Did you mean FEEvaluation<dim,";
6385 message += Utilities::int_to_string(fe_degree) + ",";
6386 message += Utilities::int_to_string(n_q_points_1d);
6387 message += "," + Utilities::int_to_string(n_components);
6388 message += ",Number>(data";
6389 if (dof_no != numbers::invalid_unsigned_int)
6390 {
6391 message +=
6392 ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
6393 message += Utilities::int_to_string(proposed_quad_comp) + ", ";
6394 message += Utilities::int_to_string(proposed_fe_comp);
6395 }
6396 message += ")?\n";
6397 std::string correct_pos;
6398 if (proposed_dof_comp != dof_no)
6399 correct_pos = " ^ ";
6400 else
6401 correct_pos = " ";
6402 if (proposed_quad_comp != this->quad_no)
6403 correct_pos += " ^ ";
6404 else
6405 correct_pos += " ";
6406 if (proposed_fe_comp != first_selected_component)
6407 correct_pos += " ^\n";
6408 else
6409 correct_pos += " \n";
6410 message += " " +
6411 correct_pos;
6412 }
6413 // ok, did not find the numbers specified by the template arguments in
6414 // the given list. Suggest correct template arguments
6415 const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
6416 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6417 message += "Wrong template arguments:\n";
6418 message += " Did you mean FEEvaluation<dim,";
6419 message +=
6420 Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
6421 message += Utilities::int_to_string(proposed_n_q_points_1d);
6422 message += "," + Utilities::int_to_string(n_components);
6423 message += ",Number>(data";
6424 if (dof_no != numbers::invalid_unsigned_int)
6425 {
6426 message += ", " + Utilities::int_to_string(dof_no) + ", ";
6427 message += Utilities::int_to_string(this->quad_no);
6428 message += ", " + Utilities::int_to_string(first_selected_component);
6429 }
6430 message += ")?\n";
6431 std::string correct_pos;
6432 if (this->data->data.front().fe_degree !=
6433 static_cast<unsigned int>(fe_degree))
6434 correct_pos = " ^";
6435 else
6436 correct_pos = " ";
6437 if (proposed_n_q_points_1d != n_q_points_1d)
6438 correct_pos += " ^\n";
6439 else
6440 correct_pos += " \n";
6441 message += " " + correct_pos;
6442
6443 Assert(static_cast<unsigned int>(fe_degree) ==
6444 this->data->data.front().fe_degree &&
6445 n_q_points == this->n_quadrature_points,
6446 ExcMessage(message));
6447 }
6448 if (dof_no != numbers::invalid_unsigned_int)
6450 n_q_points,
6451 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6452# endif
6453}
6454
6455
6456
6457template <int dim,
6458 int fe_degree,
6459 int n_q_points_1d,
6460 int n_components_,
6461 typename Number,
6462 typename VectorizedArrayType>
6463inline void
6464FEEvaluation<dim,
6465 fe_degree,
6466 n_q_points_1d,
6467 n_components_,
6468 Number,
6469 VectorizedArrayType>::reinit(const unsigned int cell_index)
6470{
6471 Assert(this->matrix_free != nullptr,
6472 ExcMessage("FEEvaluation was initialized without a matrix-free object."
6473 " Integer indexing is not possible."));
6474
6475 Assert(this->dof_info != nullptr, ExcNotInitialized());
6476 Assert(this->mapping_data != nullptr, ExcNotInitialized());
6477 this->cell = cell_index;
6478 this->cell_type =
6479 this->matrix_free->get_mapping_info().get_cell_type(cell_index);
6480
6481 const unsigned int offsets =
6482 this->mapping_data->data_index_offsets[cell_index];
6483 this->jacobian = &this->mapping_data->jacobians[0][offsets];
6484 this->J_value = &this->mapping_data->JxW_values[offsets];
6485 if (!this->mapping_data->jacobian_gradients[0].empty())
6486 {
6487 this->jacobian_gradients =
6488 this->mapping_data->jacobian_gradients[0].data() + offsets;
6489 this->jacobian_gradients_non_inverse =
6490 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
6491 }
6492
6493 if (this->matrix_free->n_active_entries_per_cell_batch(this->cell) == n_lanes)
6494 {
6496 for (unsigned int i = 0; i < n_lanes; ++i)
6497 this->cell_ids[i] = cell_index * n_lanes + i;
6498 }
6499 else
6500 {
6501 unsigned int i = 0;
6502 for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell);
6503 ++i)
6504 this->cell_ids[i] = cell_index * n_lanes + i;
6505 for (; i < n_lanes; ++i)
6506 this->cell_ids[i] = numbers::invalid_unsigned_int;
6507 }
6508
6509 if (this->mapping_data->quadrature_points.empty() == false)
6510 this->quadrature_points =
6511 &this->mapping_data->quadrature_points
6512 [this->mapping_data->quadrature_point_offsets[this->cell]];
6513
6514# ifdef DEBUG
6515 this->is_reinitialized = true;
6516 this->dof_values_initialized = false;
6517 this->values_quad_initialized = false;
6518 this->gradients_quad_initialized = false;
6519 this->hessians_quad_initialized = false;
6520# endif
6521}
6522
6523
6524
6525template <int dim,
6526 int fe_degree,
6527 int n_q_points_1d,
6528 int n_components_,
6529 typename Number,
6530 typename VectorizedArrayType>
6531inline void
6532FEEvaluation<dim,
6533 fe_degree,
6534 n_q_points_1d,
6535 n_components_,
6536 Number,
6537 VectorizedArrayType>::reinit(const std::array<unsigned int,
6538 n_lanes> &cell_ids)
6539{
6540 Assert(this->dof_info != nullptr, ExcNotInitialized());
6541 Assert(this->mapping_data != nullptr, ExcNotInitialized());
6542
6543 this->cell = numbers::invalid_unsigned_int;
6544 this->cell_ids = cell_ids;
6545
6546 // determine type of cell batch
6548
6549 for (unsigned int v = 0; v < n_lanes; ++v)
6550 {
6551 const unsigned int cell_index = cell_ids[v];
6552
6554 continue;
6555
6556 this->cell_type =
6557 std::max(this->cell_type,
6558 this->matrix_free->get_mapping_info().get_cell_type(
6559 cell_index / n_lanes));
6560 }
6561
6562 // allocate memory for internal data storage
6563 if (this->mapped_geometry == nullptr)
6564 this->mapped_geometry =
6565 std::make_shared<internal::MatrixFreeFunctions::
6566 MappingDataOnTheFly<dim, VectorizedArrayType>>();
6567
6568 auto &mapping_storage = this->mapped_geometry->get_data_storage();
6569
6570 auto &this_jacobian_data = mapping_storage.jacobians[0];
6571 auto &this_J_value_data = mapping_storage.JxW_values;
6572 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
6573 auto &this_jacobian_gradients_non_inverse_data =
6574 mapping_storage.jacobian_gradients_non_inverse[0];
6575 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
6576
6578 {
6579 if (this_jacobian_data.size() != 2)
6580 this_jacobian_data.resize_fast(2);
6581
6582 if (this_J_value_data.size() != 1)
6583 this_J_value_data.resize_fast(1);
6584
6585 const auto &update_flags_cells =
6586 this->matrix_free->get_mapping_info().update_flags_cells;
6587
6588 if (update_flags_cells & update_jacobian_grads &&
6589 this_jacobian_gradients_data.size() != 1)
6590 {
6591 this_jacobian_gradients_data.resize_fast(1);
6592 this_jacobian_gradients_non_inverse_data.resize_fast(1);
6593 }
6594
6595 if (update_flags_cells & update_quadrature_points &&
6596 this_quadrature_points_data.size() != 1)
6597 this_quadrature_points_data.resize_fast(1);
6598 }
6599 else
6600 {
6601 if (this_jacobian_data.size() != this->n_quadrature_points)
6602 this_jacobian_data.resize_fast(this->n_quadrature_points);
6603
6604 if (this_J_value_data.size() != this->n_quadrature_points)
6605 this_J_value_data.resize_fast(this->n_quadrature_points);
6606
6607 const auto &update_flags_cells =
6608 this->matrix_free->get_mapping_info().update_flags_cells;
6609
6610 if (update_flags_cells & update_jacobian_grads &&
6611 this_jacobian_gradients_data.size() != this->n_quadrature_points)
6612 {
6613 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
6614 this_jacobian_gradients_non_inverse_data.resize_fast(
6615 this->n_quadrature_points);
6616 }
6617
6618 if (update_flags_cells & update_quadrature_points &&
6619 this_quadrature_points_data.size() != this->n_quadrature_points)
6620 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
6621 }
6622
6623 // set pointers to internal data storage
6624 this->jacobian = this_jacobian_data.data();
6625 this->J_value = this_J_value_data.data();
6626 this->jacobian_gradients = this_jacobian_gradients_data.data();
6627 this->jacobian_gradients_non_inverse =
6628 this_jacobian_gradients_non_inverse_data.data();
6629 this->quadrature_points = this_quadrature_points_data.data();
6630
6631 // fill internal data storage lane by lane
6632 for (unsigned int v = 0; v < n_lanes; ++v)
6633 {
6634 const unsigned int cell_index = cell_ids[v];
6635
6637 continue;
6638
6639 const unsigned int cell_batch_index = cell_index / n_lanes;
6640 const unsigned int offsets =
6641 this->mapping_data->data_index_offsets[cell_batch_index];
6642 const unsigned int lane = cell_index % n_lanes;
6643
6644 if (this->cell_type <=
6646 {
6647 // case that all cells are Cartesian or affine
6648 for (unsigned int q = 0; q < 2; ++q)
6649 for (unsigned int i = 0; i < dim; ++i)
6650 for (unsigned int j = 0; j < dim; ++j)
6651 this_jacobian_data[q][i][j][v] =
6652 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
6653
6654 const unsigned int q = 0;
6655
6656 this_J_value_data[q][v] =
6657 this->mapping_data->JxW_values[offsets + q][lane];
6658
6659 const auto &update_flags_cells =
6660 this->matrix_free->get_mapping_info().update_flags_cells;
6661
6662 if (update_flags_cells & update_jacobian_grads)
6663 {
6664 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6665 for (unsigned int j = 0; j < dim; ++j)
6666 this_jacobian_gradients_data[q][i][j][v] =
6667 this->mapping_data
6668 ->jacobian_gradients[0][offsets + q][i][j][lane];
6669
6670 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6671 for (unsigned int j = 0; j < dim; ++j)
6672 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6673 this->mapping_data
6674 ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
6675 [lane];
6676 }
6677
6678 if (update_flags_cells & update_quadrature_points)
6679 for (unsigned int i = 0; i < dim; ++i)
6680 this_quadrature_points_data[q][i][v] =
6681 this->mapping_data->quadrature_points
6682 [this->mapping_data
6683 ->quadrature_point_offsets[cell_batch_index] +
6684 q][i][lane];
6685 }
6686 else
6687 {
6688 // general case that at least one cell is not Cartesian or affine
6689 const auto cell_type =
6690 this->matrix_free->get_mapping_info().get_cell_type(
6691 cell_batch_index);
6692
6693 for (unsigned int q = 0; q < this->n_quadrature_points; ++q)
6694 {
6695 const unsigned int q_src =
6696 (cell_type <=
6698 0 :
6699 q;
6700
6701 this_J_value_data[q][v] =
6702 this->mapping_data->JxW_values[offsets + q_src][lane];
6703
6704 for (unsigned int i = 0; i < dim; ++i)
6705 for (unsigned int j = 0; j < dim; ++j)
6706 this_jacobian_data[q][i][j][v] =
6707 this->mapping_data
6708 ->jacobians[0][offsets + q_src][i][j][lane];
6709
6710 const auto &update_flags_cells =
6711 this->matrix_free->get_mapping_info().update_flags_cells;
6712
6713 if (update_flags_cells & update_jacobian_grads)
6714 {
6715 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6716 for (unsigned int j = 0; j < dim; ++j)
6717 this_jacobian_gradients_data[q][i][j][v] =
6718 this->mapping_data
6719 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
6720
6721 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6722 for (unsigned int j = 0; j < dim; ++j)
6723 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6724 this->mapping_data
6725 ->jacobian_gradients_non_inverse[0][offsets + q_src]
6726 [i][j][lane];
6727 }
6728
6729 if (update_flags_cells & update_quadrature_points)
6730 {
6731 if (cell_type <=
6733 {
6734 // affine case: quadrature points are not available but
6735 // have to be computed from the corner point and the
6736 // Jacobian
6738 this->mapping_data->quadrature_points
6739 [this->mapping_data
6740 ->quadrature_point_offsets[cell_batch_index] +
6741 0];
6742
6744 this->mapping_data->jacobians[0][offsets + 1];
6746 for (unsigned int d = 0; d < dim; ++d)
6747 point[d] +=
6748 jac[d][d] *
6749 static_cast<Number>(
6750 this->descriptor->quadrature.point(q)[d]);
6751 else
6752 for (unsigned int d = 0; d < dim; ++d)
6753 for (unsigned int e = 0; e < dim; ++e)
6754 point[d] +=
6755 jac[d][e] *
6756 static_cast<Number>(
6757 this->descriptor->quadrature.point(q)[e]);
6758
6759 for (unsigned int i = 0; i < dim; ++i)
6760 this_quadrature_points_data[q][i][v] = point[i][lane];
6761 }
6762 else
6763 {
6764 // general case: quadrature points are available
6765 for (unsigned int i = 0; i < dim; ++i)
6766 this_quadrature_points_data[q][i][v] =
6767 this->mapping_data->quadrature_points
6768 [this->mapping_data
6769 ->quadrature_point_offsets[cell_batch_index] +
6770 q][i][lane];
6771 }
6772 }
6773 }
6774 }
6775 }
6776
6777# ifdef DEBUG
6778 this->is_reinitialized = true;
6779 this->dof_values_initialized = false;
6780 this->values_quad_initialized = false;
6781 this->gradients_quad_initialized = false;
6782 this->hessians_quad_initialized = false;
6783# endif
6784}
6785
6786
6787
6788template <int dim,
6789 int fe_degree,
6790 int n_q_points_1d,
6791 int n_components_,
6792 typename Number,
6793 typename VectorizedArrayType>
6794template <bool level_dof_access>
6795inline void
6796FEEvaluation<dim,
6797 fe_degree,
6798 n_q_points_1d,
6799 n_components_,
6800 Number,
6801 VectorizedArrayType>::
6803{
6804 Assert(this->matrix_free == nullptr,
6805 ExcMessage("Cannot use initialization from cell iterator if "
6806 "initialized from MatrixFree object. Use variant for "
6807 "on the fly computation with arguments as for FEValues "
6808 "instead"));
6809 Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
6810 this->mapped_geometry->reinit(
6811 static_cast<typename Triangulation<dim>::cell_iterator>(cell));
6812 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6813 if (level_dof_access)
6814 cell->get_mg_dof_indices(this->local_dof_indices);
6815 else
6816 cell->get_dof_indices(this->local_dof_indices);
6817
6818# ifdef DEBUG
6819 this->is_reinitialized = true;
6820# endif
6821}
6822
6823
6824
6825template <int dim,
6826 int fe_degree,
6827 int n_q_points_1d,
6828 int n_components_,
6829 typename Number,
6830 typename VectorizedArrayType>
6831inline void
6832FEEvaluation<dim,
6833 fe_degree,
6834 n_q_points_1d,
6835 n_components_,
6836 Number,
6837 VectorizedArrayType>::
6838 reinit(const typename Triangulation<dim>::cell_iterator &cell)
6839{
6840 Assert(this->matrix_free == 0,
6841 ExcMessage("Cannot use initialization from cell iterator if "
6842 "initialized from MatrixFree object. Use variant for "
6843 "on the fly computation with arguments as for FEValues "
6844 "instead"));
6845 Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
6846 this->mapped_geometry->reinit(cell);
6847
6848# ifdef DEBUG
6849 this->is_reinitialized = true;
6850# endif
6851}
6852
6853
6854
6855template <int dim,
6856 int fe_degree,
6857 int n_q_points_1d,
6858 int n_components_,
6859 typename Number,
6860 typename VectorizedArrayType>
6861inline void
6862FEEvaluation<dim,
6863 fe_degree,
6864 n_q_points_1d,
6865 n_components_,
6866 Number,
6867 VectorizedArrayType>::
6868 evaluate(const EvaluationFlags::EvaluationFlags evaluation_flags)
6869{
6870# ifdef DEBUG
6871 Assert(this->dof_values_initialized == true,
6873# endif
6874 evaluate(this->values_dofs, evaluation_flags);
6875}
6876
6877
6878
6879template <int dim,
6880 int fe_degree,
6881 int n_q_points_1d,
6882 int n_components_,
6883 typename Number,
6884 typename VectorizedArrayType>
6885inline void
6886FEEvaluation<dim,
6887 fe_degree,
6888 n_q_points_1d,
6889 n_components_,
6890 Number,
6891 VectorizedArrayType>::
6892 evaluate(const VectorizedArrayType *values_array,
6893 const EvaluationFlags::EvaluationFlags evaluation_flag)
6894{
6895 const bool hessians_on_general_cells =
6896 evaluation_flag & EvaluationFlags::hessians &&
6897 (this->cell_type > internal::MatrixFreeFunctions::affine);
6898 EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
6899 if (hessians_on_general_cells)
6900 evaluation_flag_actual |= EvaluationFlags::gradients;
6901
6902 if (this->data->element_type ==
6904 evaluation_flag & EvaluationFlags::gradients &&
6905 (this->cell_type > internal::MatrixFreeFunctions::affine))
6906 evaluation_flag_actual |= EvaluationFlags::values;
6907
6908 if constexpr (fe_degree > -1)
6909 {
6911 template run<fe_degree, n_q_points_1d>(n_components,
6912 evaluation_flag_actual,
6913 values_array,
6914 *this);
6915 }
6916 else
6917 {
6919 n_components,
6920 evaluation_flag_actual,
6921 const_cast<VectorizedArrayType *>(values_array),
6922 *this);
6923 }
6924
6925# ifdef DEBUG
6926 this->values_quad_initialized =
6927 evaluation_flag_actual & EvaluationFlags::values;
6928 this->gradients_quad_initialized =
6929 evaluation_flag_actual & EvaluationFlags::gradients;
6930 this->hessians_quad_initialized =
6931 evaluation_flag_actual & EvaluationFlags::hessians;
6932# endif
6933}
6934
6935
6936namespace internal
6937{
6941 template <typename Number,
6942 typename VectorizedArrayType,
6943 typename VectorType,
6944 typename EvaluatorType,
6945 std::enable_if_t<internal::has_begin<VectorType> &&
6947 VectorType> * = nullptr>
6948 VectorizedArrayType *
6949 check_vector_access_inplace(const EvaluatorType &fe_eval, VectorType &vector)
6950 {
6951 // for user-defined cell batches this functionality is not supported
6952 if (fe_eval.get_current_cell_index() == numbers::invalid_unsigned_int)
6953 return nullptr;
6954
6955 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
6956 const auto &dof_info = fe_eval.get_dof_info();
6957
6958 // If the index storage is interleaved and contiguous and the vector
6959 // storage has the correct alignment, we can directly pass the pointer
6960 // into the vector to the evaluate() and integrate() calls, without
6961 // reading the vector entries into a separate data field. This saves some
6962 // operations.
6963 if (std::is_same_v<typename VectorType::value_type, Number> &&
6964 dof_info.index_storage_variants
6967 interleaved_contiguous &&
6968 reinterpret_cast<std::size_t>(
6969 vector.begin() +
6970 dof_info.dof_indices_contiguous
6971 [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
6972 [cell * VectorizedArrayType::size()]) %
6973 sizeof(VectorizedArrayType) ==
6974 0)
6975 {
6976 return reinterpret_cast<VectorizedArrayType *>(
6977 vector.begin() +
6978 dof_info.dof_indices_contiguous
6980 [cell * VectorizedArrayType::size()] +
6982 [fe_eval.get_active_fe_index()]
6983 [fe_eval.get_first_selected_component()] *
6984 VectorizedArrayType::size());
6985 }
6986 else
6987 return nullptr;
6988 }
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 &, VectorType &)
7002 {
7003 return nullptr;
7004 }
7005} // namespace internal
7006
7007
7008
7009template <int dim,
7010 int fe_degree,
7011 int n_q_points_1d,
7012 int n_components_,
7013 typename Number,
7014 typename VectorizedArrayType>
7015template <typename VectorType>
7016inline void
7017FEEvaluation<dim,
7018 fe_degree,
7019 n_q_points_1d,
7020 n_components_,
7021 Number,
7022 VectorizedArrayType>::
7023 gather_evaluate(const VectorType &input_vector,
7024 const EvaluationFlags::EvaluationFlags evaluation_flag)
7025{
7026 const VectorizedArrayType *src_ptr =
7027 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7028 *this, input_vector);
7029 if (src_ptr != nullptr)
7030 evaluate(src_ptr, evaluation_flag);
7031 else
7032 {
7033 this->read_dof_values(input_vector);
7034 evaluate(this->begin_dof_values(), evaluation_flag);
7035 }
7036}
7037
7038
7039
7040template <int dim,
7041 int fe_degree,
7042 int n_q_points_1d,
7043 int n_components_,
7044 typename Number,
7045 typename VectorizedArrayType>
7046inline void
7047FEEvaluation<dim,
7048 fe_degree,
7049 n_q_points_1d,
7050 n_components_,
7051 Number,
7052 VectorizedArrayType>::
7053 integrate(const EvaluationFlags::EvaluationFlags integration_flag)
7054{
7055 integrate(integration_flag, this->values_dofs);
7056
7057# ifdef DEBUG
7058 this->dof_values_initialized = true;
7059# endif
7060}
7061
7062
7063
7064template <int dim,
7065 int fe_degree,
7066 int n_q_points_1d,
7067 int n_components_,
7068 typename Number,
7069 typename VectorizedArrayType>
7070inline void
7071FEEvaluation<dim,
7072 fe_degree,
7073 n_q_points_1d,
7074 n_components_,
7075 Number,
7076 VectorizedArrayType>::
7077 integrate(const EvaluationFlags::EvaluationFlags integration_flag,
7078 VectorizedArrayType *values_array,
7079 const bool sum_into_values_array)
7080{
7081# ifdef DEBUG
7082 if (integration_flag & EvaluationFlags::values)
7083 Assert(this->values_quad_submitted == true,
7085 if (integration_flag & EvaluationFlags::gradients)
7086 Assert(this->gradients_quad_submitted == true,
7088 if ((integration_flag & EvaluationFlags::hessians) != 0u)
7089 Assert(this->hessians_quad_submitted == true,
7091# endif
7092 Assert(this->matrix_free != nullptr ||
7093 this->mapped_geometry->is_initialized(),
7095
7096 Assert(
7097 (integration_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients |
7099 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7100 "EvaluationFlags::hessians are supported."));
7101
7102 EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
7103 if (integration_flag & EvaluationFlags::hessians &&
7104 (this->cell_type > internal::MatrixFreeFunctions::affine))
7105 {
7106 unsigned int size = n_components * dim * n_q_points;
7107 if ((integration_flag & EvaluationFlags::gradients) != 0u)
7108 {
7109 for (unsigned int i = 0; i < size; ++i)
7110 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7111 }
7112 else
7113 {
7114 for (unsigned int i = 0; i < size; ++i)
7115 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7116 integration_flag_actual |= EvaluationFlags::gradients;
7117 }
7118 }
7119
7120 if (n_components == dim &&
7121 this->data->element_type ==
7123 integration_flag & EvaluationFlags::gradients &&
7124 this->cell_type > internal::MatrixFreeFunctions::affine &&
7125 this->divergence_is_requested == false)
7126 {
7127 unsigned int size = n_components * n_q_points;
7128 if ((integration_flag & EvaluationFlags::values) != 0u)
7129 {
7130 for (unsigned int i = 0; i < size; ++i)
7131 this->values_quad[i] += this->values_from_gradients_quad[i];
7132 }
7133 else
7134 {
7135 for (unsigned int i = 0; i < size; ++i)
7136 this->values_quad[i] = this->values_from_gradients_quad[i];
7137 integration_flag_actual |= EvaluationFlags::values;
7138 }
7139 }
7140
7141 if constexpr (fe_degree > -1)
7142 {
7144 template run<fe_degree, n_q_points_1d>(n_components,
7145 integration_flag_actual,
7146 values_array,
7147 *this,
7148 sum_into_values_array);
7149 }
7150 else
7151 {
7153 n_components,
7154 integration_flag_actual,
7155 values_array,
7156 *this,
7157 sum_into_values_array);
7158 }
7159
7160# ifdef DEBUG
7161 this->dof_values_initialized = true;
7162# endif
7163}
7164
7165
7166
7167template <int dim,
7168 int fe_degree,
7169 int n_q_points_1d,
7170 int n_components_,
7171 typename Number,
7172 typename VectorizedArrayType>
7173template <typename VectorType>
7174inline void
7175FEEvaluation<dim,
7176 fe_degree,
7177 n_q_points_1d,
7178 n_components_,
7179 Number,
7180 VectorizedArrayType>::
7181 integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
7182 VectorType &destination)
7183{
7184 VectorizedArrayType *dst_ptr =
7185 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
7186 *this, destination);
7187 if (dst_ptr != nullptr)
7188 integrate(integration_flag, dst_ptr, true);
7189 else
7190 {
7191 integrate(integration_flag, this->begin_dof_values());
7192 this->distribute_local_to_global(destination);
7193 }
7194}
7195
7196
7197
7198template <int dim,
7199 int fe_degree,
7200 int n_q_points_1d,
7201 int n_components_,
7202 typename Number,
7203 typename VectorizedArrayType>
7205FEEvaluation<dim,
7206 fe_degree,
7207 n_q_points_1d,
7208 n_components_,
7209 Number,
7210 VectorizedArrayType>::dof_indices() const
7211{
7212 return {0U, dofs_per_cell};
7213}
7214
7215
7216
7217/*-------------------------- FEFaceEvaluation ---------------------------*/
7218
7219
7220
7221template <int dim,
7222 int fe_degree,
7223 int n_q_points_1d,
7224 int n_components_,
7225 typename Number,
7226 typename VectorizedArrayType>
7227inline FEFaceEvaluation<dim,
7228 fe_degree,
7229 n_q_points_1d,
7230 n_components_,
7231 Number,
7232 VectorizedArrayType>::
7233 FEFaceEvaluation(
7235 const bool is_interior_face,
7236 const unsigned int dof_no,
7237 const unsigned int quad_no,
7238 const unsigned int first_selected_component,
7239 const unsigned int active_fe_index,
7240 const unsigned int active_quad_index,
7241 const unsigned int face_type)
7242 : BaseClass(matrix_free,
7243 dof_no,
7244 first_selected_component,
7245 quad_no,
7246 fe_degree,
7247 static_n_q_points,
7248 is_interior_face,
7249 active_fe_index,
7250 active_quad_index,
7251 face_type)
7252 , dofs_per_component(this->data->dofs_per_component_on_cell)
7253 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7254 , n_q_points(this->n_quadrature_points)
7255{}
7256
7257
7258
7259template <int dim,
7260 int fe_degree,
7261 int n_q_points_1d,
7262 int n_components_,
7263 typename Number,
7264 typename VectorizedArrayType>
7265inline FEFaceEvaluation<dim,
7266 fe_degree,
7267 n_q_points_1d,
7268 n_components_,
7269 Number,
7270 VectorizedArrayType>::
7271 FEFaceEvaluation(
7273 const std::pair<unsigned int, unsigned int> &rang