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