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