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