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