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