Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_values_views.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2023 - 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#ifndef dealii_fe_values_views_h
16#define dealii_fe_values_views_h
17
18#include <deal.II/base/config.h>
19
21#include <deal.II/base/lazy.h>
24#include <deal.II/base/tensor.h>
25
28
30
31#include <type_traits>
32#include <vector>
33
35
36// Forward declaration
37#ifndef DOXYGEN
38template <int dim, int spacedim = dim>
39class FEValuesBase;
40#endif
41
42namespace internal
43{
48 template <int dim, typename NumberType = double>
49 struct CurlType;
50
57 template <typename NumberType>
59 {
61 };
62
69 template <typename NumberType>
71 {
73 };
74
81 template <typename NumberType>
83 {
85 };
86} // namespace internal
87
88
89
111namespace FEValuesViews
112{
124 template <int dim, int spacedim = dim>
125 class Scalar
126 {
127 public:
133 using value_type = double;
134
141
148
155
162 template <typename Number>
164
171 template <typename Number>
174
181 template <typename Number>
184
191 template <typename Number>
194
201 template <typename Number>
204
231
235 Scalar();
236
242 Scalar(const FEValuesBase<dim, spacedim> &fe_values_base,
243 const unsigned int component);
244
249 Scalar(const Scalar<dim, spacedim> &) = delete;
250
254 // NOLINTNEXTLINE OSX does not compile with noexcept
256
260 ~Scalar() = default;
261
266 Scalar &
268
272 Scalar &
273 operator=(Scalar<dim, spacedim> &&) noexcept = default;
274
289 value(const unsigned int shape_function, const unsigned int q_point) const;
290
302 gradient(const unsigned int shape_function,
303 const unsigned int q_point) const;
304
316 hessian(const unsigned int shape_function,
317 const unsigned int q_point) const;
318
330 third_derivative(const unsigned int shape_function,
331 const unsigned int q_point) const;
332
350 template <typename Number>
351 void
352 get_function_values(const ReadVector<Number> &fe_function,
353 std::vector<solution_value_type<Number>> &values) const;
354
389 template <class InputVector>
390 void
392 const InputVector &dof_values,
393 std::vector<solution_value_type<typename InputVector::value_type>>
394 &values) const;
395
413 template <typename Number>
414 void
416 const ReadVector<Number> &fe_function,
417 std::vector<solution_gradient_type<Number>> &gradients) const;
418
425 template <class InputVector>
426 void
428 const InputVector &dof_values,
429 std::vector<solution_gradient_type<typename InputVector::value_type>>
430 &gradients) const;
431
449 template <typename Number>
450 void
452 const ReadVector<Number> &fe_function,
453 std::vector<solution_hessian_type<Number>> &hessians) const;
454
461 template <class InputVector>
462 void
464 const InputVector &dof_values,
465 std::vector<solution_hessian_type<typename InputVector::value_type>>
466 &hessians) const;
467
468
487 template <typename Number>
488 void
490 const ReadVector<Number> &fe_function,
491 std::vector<solution_laplacian_type<Number>> &laplacians) const;
492
499 template <class InputVector>
500 void
502 const InputVector &dof_values,
503 std::vector<solution_laplacian_type<typename InputVector::value_type>>
504 &laplacians) const;
505
506
525 template <typename Number>
526 void
528 const ReadVector<Number> &fe_function,
529 std::vector<solution_third_derivative_type<Number>> &third_derivatives)
530 const;
531
538 template <class InputVector>
539 void
541 const InputVector &dof_values,
542 std::vector<
543 solution_third_derivative_type<typename InputVector::value_type>>
544 &third_derivatives) const;
545
546
547 private:
551 SmartPointer<const FEValuesBase<dim, spacedim>> fe_values;
552
557 unsigned int component;
558
563 };
564
565
566
596 template <int dim, int spacedim = dim>
597 class Vector
598 {
599 public:
606
616
628
634 using divergence_type = double;
635
642 using curl_type = typename ::internal::CurlType<spacedim>::type;
643
650
657
664 template <typename Number>
666
673 template <typename Number>
676
683 template <typename Number>
686
693 template <typename Number>
696
703 template <typename Number>
706
713 template <typename Number>
715
722 template <typename Number>
725
732 template <typename Number>
735
741 {
751
761 unsigned int row_index[spacedim];
762
773 };
774
778 Vector();
779
788 Vector(const FEValuesBase<dim, spacedim> &fe_values_base,
789 const unsigned int first_vector_component);
790
795 Vector(const Vector<dim, spacedim> &) = delete;
796
800 // NOLINTNEXTLINE OSX does not compile with noexcept
802
806 ~Vector() = default;
807
812 Vector &
814
818 // NOLINTNEXTLINE OSX does not compile with noexcept
819 Vector &
820 operator=(Vector<dim, spacedim> &&) = default; // NOLINT
821
839 value(const unsigned int shape_function, const unsigned int q_point) const;
840
855 gradient(const unsigned int shape_function,
856 const unsigned int q_point) const;
857
874 symmetric_gradient(const unsigned int shape_function,
875 const unsigned int q_point) const;
876
888 divergence(const unsigned int shape_function,
889 const unsigned int q_point) const;
890
912 curl(const unsigned int shape_function, const unsigned int q_point) const;
913
925 hessian(const unsigned int shape_function,
926 const unsigned int q_point) const;
927
939 third_derivative(const unsigned int shape_function,
940 const unsigned int q_point) const;
941
959 template <typename Number>
960 void
961 get_function_values(const ReadVector<Number> &fe_function,
962 std::vector<solution_value_type<Number>> &values) const;
963
998 template <class InputVector>
999 void
1001 const InputVector &dof_values,
1003 &values) const;
1004
1022 template <typename Number>
1023 void
1025 const ReadVector<Number> &fe_function,
1026 std::vector<solution_gradient_type<Number>> &gradients) const;
1027
1034 template <class InputVector>
1035 void
1037 const InputVector &dof_values,
1039 &gradients) const;
1040
1064 template <typename Number>
1065 void
1066 get_function_symmetric_gradients(
1067 const ReadVector<Number> &fe_function,
1069 &symmetric_gradients) const;
1070
1077 template <class InputVector>
1078 void
1079 get_function_symmetric_gradients_from_local_dof_values(
1080 const InputVector &dof_values,
1081 std::vector<
1083 &symmetric_gradients) const;
1084
1103 template <typename Number>
1104 void
1105 get_function_divergences(
1106 const ReadVector<Number> &fe_function,
1107 std::vector<solution_divergence_type<Number>> &divergences) const;
1108
1115 template <class InputVector>
1116 void
1117 get_function_divergences_from_local_dof_values(
1118 const InputVector &dof_values,
1120 &divergences) const;
1121
1140 template <typename Number>
1141 void
1142 get_function_curls(const ReadVector<Number> &fe_function,
1143 std::vector<solution_curl_type<Number>> &curls) const;
1144
1151 template <class InputVector>
1152 void
1153 get_function_curls_from_local_dof_values(
1154 const InputVector &dof_values,
1156 const;
1157
1175 template <typename Number>
1176 void
1178 const ReadVector<Number> &fe_function,
1179 std::vector<solution_hessian_type<Number>> &hessians) const;
1180
1187 template <class InputVector>
1188 void
1190 const InputVector &dof_values,
1192 &hessians) const;
1193
1212 template <typename Number>
1213 void
1215 const ReadVector<Number> &fe_function,
1216 std::vector<solution_laplacian_type<Number>> &laplacians) const;
1217
1224 template <class InputVector>
1225 void
1227 const InputVector &dof_values,
1229 &laplacians) const;
1230
1249 template <typename Number>
1250 void
1252 const ReadVector<Number> &fe_function,
1253 std::vector<solution_third_derivative_type<Number>> &third_derivatives)
1254 const;
1255
1262 template <class InputVector>
1263 void
1265 const InputVector &dof_values,
1266 std::vector<
1268 &third_derivatives) const;
1269
1270 private:
1275
1281
1285 std::vector<ShapeFunctionData> shape_function_data;
1286 };
1287
1288
1289 template <int rank, int dim, int spacedim = dim>
1291
1314 template <int dim, int spacedim>
1315 class SymmetricTensor<2, dim, spacedim>
1316 {
1317 public:
1325
1336
1343 template <typename Number>
1345
1352 template <typename Number>
1355
1356
1362 {
1371 bool is_nonzero_shape_function_component
1372 [value_type::n_independent_components];
1373
1383 unsigned int row_index[value_type::n_independent_components];
1384
1394
1399 };
1400
1405
1415 SymmetricTensor(const FEValuesBase<dim, spacedim> &fe_values_base,
1416 const unsigned int first_tensor_component);
1417
1423
1427 // NOLINTNEXTLINE OSX does not compile with noexcept
1429
1436
1442
1461 value(const unsigned int shape_function, const unsigned int q_point) const;
1462
1477 divergence(const unsigned int shape_function,
1478 const unsigned int q_point) const;
1479
1497 template <typename Number>
1498 void
1499 get_function_values(const ReadVector<Number> &fe_function,
1500 std::vector<solution_value_type<Number>> &values) const;
1501
1536 template <class InputVector>
1537 void
1538 get_function_values_from_local_dof_values(
1539 const InputVector &dof_values,
1540 std::vector<solution_value_type<typename InputVector::value_type>>
1541 &values) const;
1542
1564 template <typename Number>
1565 void
1566 get_function_divergences(
1567 const ReadVector<Number> &fe_function,
1568 std::vector<solution_divergence_type<Number>> &divergences) const;
1569
1576 template <class InputVector>
1577 void
1578 get_function_divergences_from_local_dof_values(
1579 const InputVector &dof_values,
1580 std::vector<solution_divergence_type<typename InputVector::value_type>>
1581 &divergences) const;
1582
1583 private:
1587 SmartPointer<const FEValuesBase<dim, spacedim>> fe_values;
1588
1593 unsigned int first_tensor_component;
1594
1598 std::vector<ShapeFunctionData> shape_function_data;
1599 };
1600
1601
1602 template <int rank, int dim, int spacedim = dim>
1603 class Tensor;
1604
1623 template <int dim, int spacedim>
1624 class Tensor<2, dim, spacedim>
1625 {
1626 public:
1632
1637
1643
1650 template <typename Number>
1652
1659 template <typename Number>
1662
1669 template <typename Number>
1672
1673
1679 {
1688 bool is_nonzero_shape_function_component
1689 [value_type::n_independent_components];
1690
1700 unsigned int row_index[value_type::n_independent_components];
1701
1711
1716 };
1717
1721 Tensor();
1722
1728
1732 // NOLINTNEXTLINE OSX does not compile with noexcept
1734
1738 ~Tensor() = default;
1739
1749 Tensor(const FEValuesBase<dim, spacedim> &fe_values_base,
1750 const unsigned int first_tensor_component);
1751
1752
1757 Tensor &
1759
1763 Tensor &
1764 operator=(Tensor<2, dim, spacedim> &&) = default; // NOLINT
1765
1783 value(const unsigned int shape_function, const unsigned int q_point) const;
1784
1799 divergence(const unsigned int shape_function,
1800 const unsigned int q_point) const;
1801
1816 gradient(const unsigned int shape_function,
1817 const unsigned int q_point) const;
1818
1836 template <typename Number>
1837 void
1838 get_function_values(const ReadVector<Number> &fe_function,
1839 std::vector<solution_value_type<Number>> &values) const;
1840
1875 template <class InputVector>
1876 void
1877 get_function_values_from_local_dof_values(
1878 const InputVector &dof_values,
1880 &values) const;
1881
1903 template <typename Number>
1904 void
1905 get_function_divergences(
1906 const ReadVector<Number> &fe_function,
1907 std::vector<solution_divergence_type<Number>> &divergences) const;
1908
1915 template <class InputVector>
1916 void
1917 get_function_divergences_from_local_dof_values(
1918 const InputVector &dof_values,
1920 &divergences) const;
1921
1938 template <typename Number>
1939 void
1940 get_function_gradients(
1941 const ReadVector<Number> &fe_function,
1942 std::vector<solution_gradient_type<Number>> &gradients) const;
1943
1950 template <class InputVector>
1951 void
1952 get_function_gradients_from_local_dof_values(
1953 const InputVector &dof_values,
1955 &gradients) const;
1956
1957 private:
1962
1968
1972 std::vector<ShapeFunctionData> shape_function_data;
1973 };
1974
1975} // namespace FEValuesViews
1976
1977
1978namespace internal
1979{
1981 {
1986 template <int dim, int spacedim, typename Extractor>
1988 {};
1989
1997 template <int dim, int spacedim>
1998 struct ViewType<dim, spacedim, FEValuesExtractors::Scalar>
1999 {
2000 using type = typename ::FEValuesViews::Scalar<dim, spacedim>;
2001 };
2002
2010 template <int dim, int spacedim>
2011 struct ViewType<dim, spacedim, FEValuesExtractors::Vector>
2012 {
2013 using type = typename ::FEValuesViews::Vector<dim, spacedim>;
2014 };
2015
2023 template <int dim, int spacedim, int rank>
2024 struct ViewType<dim, spacedim, FEValuesExtractors::Tensor<rank>>
2025 {
2026 using type = typename ::FEValuesViews::Tensor<rank, dim, spacedim>;
2027 };
2028
2036 template <int dim, int spacedim, int rank>
2037 struct ViewType<dim, spacedim, FEValuesExtractors::SymmetricTensor<rank>>
2038 {
2039 using type =
2040 typename ::FEValuesViews::SymmetricTensor<rank, dim, spacedim>;
2041 };
2042
2050 template <int dim, int spacedim>
2051 struct Cache
2052 {
2057 std::vector<Lazy<::FEValuesViews::Scalar<dim, spacedim>>> scalars;
2058 std::vector<Lazy<::FEValuesViews::Vector<dim, spacedim>>> vectors;
2059 std::vector<
2062 std::vector<Lazy<::FEValuesViews::Tensor<2, dim, spacedim>>>
2064
2068 Cache(const FEValuesBase<dim, spacedim> &fe_values);
2069 };
2070 } // namespace FEValuesViews
2071} // namespace internal
2072
2073namespace FEValuesViews
2074{
2079 template <int dim, int spacedim, typename Extractor>
2080 using View = typename ::internal::FEValuesViews::
2081 ViewType<dim, spacedim, Extractor>::type;
2082} // namespace FEValuesViews
2083
2084#ifndef DOXYGEN
2085
2086/*---------------- Inline functions: namespace FEValuesViews -----------------*/
2087
2088namespace FEValuesViews
2089{
2090 template <int dim, int spacedim>
2091 inline typename Scalar<dim, spacedim>::value_type
2092 Scalar<dim, spacedim>::value(const unsigned int shape_function,
2093 const unsigned int q_point) const
2094 {
2095 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2096 Assert(
2097 fe_values->update_flags & update_values,
2099 "update_values"))));
2100
2101 // an adaptation of the FEValuesBase::shape_value_component function
2102 // except that here we know the component as fixed and we have
2103 // pre-computed and cached a bunch of information. See the comments there.
2104 if (shape_function_data[shape_function].is_nonzero_shape_function_component)
2105 return fe_values->finite_element_output.shape_values(
2106 shape_function_data[shape_function].row_index, q_point);
2107 else
2108 return 0;
2109 }
2110
2111
2112
2113 template <int dim, int spacedim>
2114 inline typename Scalar<dim, spacedim>::gradient_type
2115 Scalar<dim, spacedim>::gradient(const unsigned int shape_function,
2116 const unsigned int q_point) const
2117 {
2118 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2119 Assert(fe_values->update_flags & update_gradients,
2121 "update_gradients")));
2122
2123 // an adaptation of the FEValuesBase::shape_grad_component
2124 // function except that here we know the component as fixed and we have
2125 // pre-computed and cached a bunch of information. See the comments there.
2126 if (shape_function_data[shape_function].is_nonzero_shape_function_component)
2127 return fe_values->finite_element_output
2128 .shape_gradients[shape_function_data[shape_function].row_index]
2129 [q_point];
2130 else
2131 return gradient_type();
2132 }
2133
2134
2135
2136 template <int dim, int spacedim>
2137 inline typename Scalar<dim, spacedim>::hessian_type
2138 Scalar<dim, spacedim>::hessian(const unsigned int shape_function,
2139 const unsigned int q_point) const
2140 {
2141 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2142 Assert(fe_values->update_flags & update_hessians,
2144 "update_hessians")));
2145
2146 // an adaptation of the FEValuesBase::shape_hessian_component
2147 // function except that here we know the component as fixed and we have
2148 // pre-computed and cached a bunch of information. See the comments there.
2149 if (shape_function_data[shape_function].is_nonzero_shape_function_component)
2150 return fe_values->finite_element_output
2151 .shape_hessians[shape_function_data[shape_function].row_index][q_point];
2152 else
2153 return hessian_type();
2154 }
2155
2156
2157
2158 template <int dim, int spacedim>
2159 inline typename Scalar<dim, spacedim>::third_derivative_type
2160 Scalar<dim, spacedim>::third_derivative(const unsigned int shape_function,
2161 const unsigned int q_point) const
2162 {
2163 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2164 Assert(fe_values->update_flags & update_3rd_derivatives,
2166 "update_3rd_derivatives")));
2167
2168 // an adaptation of the FEValuesBase::shape_3rdderivative_component
2169 // function except that here we know the component as fixed and we have
2170 // pre-computed and cached a bunch of information. See the comments there.
2171 if (shape_function_data[shape_function].is_nonzero_shape_function_component)
2172 return fe_values->finite_element_output
2173 .shape_3rd_derivatives[shape_function_data[shape_function].row_index]
2174 [q_point];
2175 else
2176 return third_derivative_type();
2177 }
2178
2179
2180
2181 template <int dim, int spacedim>
2182 inline typename Vector<dim, spacedim>::value_type
2183 Vector<dim, spacedim>::value(const unsigned int shape_function,
2184 const unsigned int q_point) const
2185 {
2186 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2187 Assert(fe_values->update_flags & update_values,
2189 "update_values")));
2190
2191 // same as for the scalar case except that we have one more index
2192 const int snc =
2193 shape_function_data[shape_function].single_nonzero_component;
2194 if (snc == -2)
2195 return value_type();
2196 else if (snc != -1)
2197 {
2198 value_type return_value;
2199 return_value[shape_function_data[shape_function]
2200 .single_nonzero_component_index] =
2201 fe_values->finite_element_output.shape_values(snc, q_point);
2202 return return_value;
2203 }
2204 else
2205 {
2206 value_type return_value;
2207 for (unsigned int d = 0; d < dim; ++d)
2208 if (shape_function_data[shape_function]
2209 .is_nonzero_shape_function_component[d])
2210 return_value[d] = fe_values->finite_element_output.shape_values(
2211 shape_function_data[shape_function].row_index[d], q_point);
2212
2213 return return_value;
2214 }
2215 }
2216
2217
2218
2219 template <int dim, int spacedim>
2221 Vector<dim, spacedim>::gradient(const unsigned int shape_function,
2222 const unsigned int q_point) const
2223 {
2224 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2225 Assert(fe_values->update_flags & update_gradients,
2227 "update_gradients")));
2228
2229 // same as for the scalar case except that we have one more index
2230 const int snc =
2231 shape_function_data[shape_function].single_nonzero_component;
2232 if (snc == -2)
2233 return gradient_type();
2234 else if (snc != -1)
2235 {
2236 gradient_type return_value;
2237 return_value[shape_function_data[shape_function]
2238 .single_nonzero_component_index] =
2239 fe_values->finite_element_output.shape_gradients[snc][q_point];
2240 return return_value;
2241 }
2242 else
2243 {
2244 gradient_type return_value;
2245 for (unsigned int d = 0; d < dim; ++d)
2246 if (shape_function_data[shape_function]
2247 .is_nonzero_shape_function_component[d])
2248 return_value[d] =
2249 fe_values->finite_element_output.shape_gradients
2250 [shape_function_data[shape_function].row_index[d]][q_point];
2251
2252 return return_value;
2253 }
2254 }
2255
2256
2257
2258 template <int dim, int spacedim>
2260 Vector<dim, spacedim>::divergence(const unsigned int shape_function,
2261 const unsigned int q_point) const
2262 {
2263 // this function works like in the case above
2264 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2265 Assert(fe_values->update_flags & update_gradients,
2267 "update_gradients")));
2268
2269 // same as for the scalar case except that we have one more index
2270 const int snc =
2271 shape_function_data[shape_function].single_nonzero_component;
2272 if (snc == -2)
2273 return divergence_type();
2274 else if (snc != -1)
2275 return fe_values->finite_element_output
2276 .shape_gradients[snc][q_point][shape_function_data[shape_function]
2277 .single_nonzero_component_index];
2278 else
2279 {
2280 divergence_type return_value = 0;
2281 for (unsigned int d = 0; d < dim; ++d)
2282 if (shape_function_data[shape_function]
2283 .is_nonzero_shape_function_component[d])
2284 return_value +=
2285 fe_values->finite_element_output.shape_gradients
2286 [shape_function_data[shape_function].row_index[d]][q_point][d];
2287
2288 return return_value;
2289 }
2290 }
2291
2292
2293
2294 template <int dim, int spacedim>
2295 inline typename Vector<dim, spacedim>::curl_type
2296 Vector<dim, spacedim>::curl(const unsigned int shape_function,
2297 const unsigned int q_point) const
2298 {
2299 // this function works like in the case above
2300
2301 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2302 Assert(fe_values->update_flags & update_gradients,
2304 "update_gradients")));
2305 // same as for the scalar case except that we have one more index
2306 const int snc =
2307 shape_function_data[shape_function].single_nonzero_component;
2308
2309 if (snc == -2)
2310 return curl_type();
2311
2312 else
2313 switch (dim)
2314 {
2315 case 1:
2316 {
2317 Assert(false,
2318 ExcMessage(
2319 "Computing the curl in 1d is not a useful operation"));
2320 return curl_type();
2321 }
2322
2323 case 2:
2324 {
2325 if (snc != -1)
2326 {
2327 curl_type return_value;
2328
2329 // the single nonzero component can only be zero or one in 2d
2330 if (shape_function_data[shape_function]
2331 .single_nonzero_component_index == 0)
2332 return_value[0] =
2333 -1.0 * fe_values->finite_element_output
2334 .shape_gradients[snc][q_point][1];
2335 else
2336 return_value[0] = fe_values->finite_element_output
2337 .shape_gradients[snc][q_point][0];
2338
2339 return return_value;
2340 }
2341
2342 else
2343 {
2344 curl_type return_value;
2345
2346 return_value[0] = 0.0;
2347
2348 if (shape_function_data[shape_function]
2349 .is_nonzero_shape_function_component[0])
2350 return_value[0] -=
2351 fe_values->finite_element_output
2352 .shape_gradients[shape_function_data[shape_function]
2353 .row_index[0]][q_point][1];
2354
2355 if (shape_function_data[shape_function]
2356 .is_nonzero_shape_function_component[1])
2357 return_value[0] +=
2358 fe_values->finite_element_output
2359 .shape_gradients[shape_function_data[shape_function]
2360 .row_index[1]][q_point][0];
2361
2362 return return_value;
2363 }
2364 }
2365
2366 case 3:
2367 {
2368 if (snc != -1)
2369 {
2370 curl_type return_value;
2371
2372 switch (shape_function_data[shape_function]
2373 .single_nonzero_component_index)
2374 {
2375 case 0:
2376 {
2377 return_value[0] = 0;
2378 return_value[1] = fe_values->finite_element_output
2379 .shape_gradients[snc][q_point][2];
2380 return_value[2] =
2381 -1.0 * fe_values->finite_element_output
2382 .shape_gradients[snc][q_point][1];
2383 return return_value;
2384 }
2385
2386 case 1:
2387 {
2388 return_value[0] =
2389 -1.0 * fe_values->finite_element_output
2390 .shape_gradients[snc][q_point][2];
2391 return_value[1] = 0;
2392 return_value[2] = fe_values->finite_element_output
2393 .shape_gradients[snc][q_point][0];
2394 return return_value;
2395 }
2396
2397 default:
2398 {
2399 return_value[0] = fe_values->finite_element_output
2400 .shape_gradients[snc][q_point][1];
2401 return_value[1] =
2402 -1.0 * fe_values->finite_element_output
2403 .shape_gradients[snc][q_point][0];
2404 return_value[2] = 0;
2405 return return_value;
2406 }
2407 }
2408 }
2409
2410 else
2411 {
2412 curl_type return_value;
2413
2414 for (unsigned int i = 0; i < dim; ++i)
2415 return_value[i] = 0.0;
2416
2417 if (shape_function_data[shape_function]
2418 .is_nonzero_shape_function_component[0])
2419 {
2420 return_value[1] +=
2421 fe_values->finite_element_output
2422 .shape_gradients[shape_function_data[shape_function]
2423 .row_index[0]][q_point][2];
2424 return_value[2] -=
2425 fe_values->finite_element_output
2426 .shape_gradients[shape_function_data[shape_function]
2427 .row_index[0]][q_point][1];
2428 }
2429
2430 if (shape_function_data[shape_function]
2431 .is_nonzero_shape_function_component[1])
2432 {
2433 return_value[0] -=
2434 fe_values->finite_element_output
2435 .shape_gradients[shape_function_data[shape_function]
2436 .row_index[1]][q_point][2];
2437 return_value[2] +=
2438 fe_values->finite_element_output
2439 .shape_gradients[shape_function_data[shape_function]
2440 .row_index[1]][q_point][0];
2441 }
2442
2443 if (shape_function_data[shape_function]
2444 .is_nonzero_shape_function_component[2])
2445 {
2446 return_value[0] +=
2447 fe_values->finite_element_output
2448 .shape_gradients[shape_function_data[shape_function]
2449 .row_index[2]][q_point][1];
2450 return_value[1] -=
2451 fe_values->finite_element_output
2452 .shape_gradients[shape_function_data[shape_function]
2453 .row_index[2]][q_point][0];
2454 }
2455
2456 return return_value;
2457 }
2458 }
2459 }
2460 // should not end up here
2462 return curl_type();
2463 }
2464
2465
2466
2467 template <int dim, int spacedim>
2469 Vector<dim, spacedim>::hessian(const unsigned int shape_function,
2470 const unsigned int q_point) const
2471 {
2472 // this function works like in the case above
2473 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2474 Assert(fe_values->update_flags & update_hessians,
2476 "update_hessians")));
2477
2478 // same as for the scalar case except that we have one more index
2479 const int snc =
2480 shape_function_data[shape_function].single_nonzero_component;
2481 if (snc == -2)
2482 return hessian_type();
2483 else if (snc != -1)
2484 {
2485 hessian_type return_value;
2486 return_value[shape_function_data[shape_function]
2487 .single_nonzero_component_index] =
2488 fe_values->finite_element_output.shape_hessians[snc][q_point];
2489 return return_value;
2490 }
2491 else
2492 {
2493 hessian_type return_value;
2494 for (unsigned int d = 0; d < dim; ++d)
2495 if (shape_function_data[shape_function]
2496 .is_nonzero_shape_function_component[d])
2497 return_value[d] =
2498 fe_values->finite_element_output.shape_hessians
2499 [shape_function_data[shape_function].row_index[d]][q_point];
2500
2501 return return_value;
2502 }
2503 }
2504
2505
2506
2507 template <int dim, int spacedim>
2509 Vector<dim, spacedim>::third_derivative(const unsigned int shape_function,
2510 const unsigned int q_point) const
2511 {
2512 // this function works like in the case above
2513 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2514 Assert(fe_values->update_flags & update_3rd_derivatives,
2516 "update_3rd_derivatives")));
2517
2518 // same as for the scalar case except that we have one more index
2519 const int snc =
2520 shape_function_data[shape_function].single_nonzero_component;
2521 if (snc == -2)
2522 return third_derivative_type();
2523 else if (snc != -1)
2524 {
2525 third_derivative_type return_value;
2526 return_value[shape_function_data[shape_function]
2527 .single_nonzero_component_index] =
2528 fe_values->finite_element_output.shape_3rd_derivatives[snc][q_point];
2529 return return_value;
2530 }
2531 else
2532 {
2533 third_derivative_type return_value;
2534 for (unsigned int d = 0; d < dim; ++d)
2535 if (shape_function_data[shape_function]
2536 .is_nonzero_shape_function_component[d])
2537 return_value[d] =
2538 fe_values->finite_element_output.shape_3rd_derivatives
2539 [shape_function_data[shape_function].row_index[d]][q_point];
2540
2541 return return_value;
2542 }
2543 }
2544
2545
2546
2547 namespace internal
2548 {
2553 inline ::SymmetricTensor<2, 1>
2554 symmetrize_single_row(const unsigned int n, const ::Tensor<1, 1> &t)
2555 {
2556 AssertIndexRange(n, 1);
2557 (void)n;
2558
2559 return {{t[0]}};
2560 }
2561
2562
2563
2564 inline ::SymmetricTensor<2, 2>
2565 symmetrize_single_row(const unsigned int n, const ::Tensor<1, 2> &t)
2566 {
2567 switch (n)
2568 {
2569 case 0:
2570 {
2571 return {{t[0], 0, t[1] / 2}};
2572 }
2573 case 1:
2574 {
2575 return {{0, t[1], t[0] / 2}};
2576 }
2577 default:
2578 {
2579 AssertIndexRange(n, 2);
2580 return {};
2581 }
2582 }
2583 }
2584
2585
2586
2587 inline ::SymmetricTensor<2, 3>
2588 symmetrize_single_row(const unsigned int n, const ::Tensor<1, 3> &t)
2589 {
2590 switch (n)
2591 {
2592 case 0:
2593 {
2594 return {{t[0], 0, 0, t[1] / 2, t[2] / 2, 0}};
2595 }
2596 case 1:
2597 {
2598 return {{0, t[1], 0, t[0] / 2, 0, t[2] / 2}};
2599 }
2600 case 2:
2601 {
2602 return {{0, 0, t[2], 0, t[0] / 2, t[1] / 2}};
2603 }
2604 default:
2605 {
2606 AssertIndexRange(n, 3);
2607 return {};
2608 }
2609 }
2610 }
2611 } // namespace internal
2612
2613
2614
2615 template <int dim, int spacedim>
2617 Vector<dim, spacedim>::symmetric_gradient(const unsigned int shape_function,
2618 const unsigned int q_point) const
2619 {
2620 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2621 Assert(fe_values->update_flags & update_gradients,
2623 "update_gradients")));
2624
2625 // same as for the scalar case except that we have one more index
2626 const int snc =
2627 shape_function_data[shape_function].single_nonzero_component;
2628 if (snc == -2)
2629 return symmetric_gradient_type();
2630 else if (snc != -1)
2631 return internal::symmetrize_single_row(
2632 shape_function_data[shape_function].single_nonzero_component_index,
2633 fe_values->finite_element_output.shape_gradients[snc][q_point]);
2634 else
2635 {
2636 gradient_type return_value;
2637 for (unsigned int d = 0; d < dim; ++d)
2638 if (shape_function_data[shape_function]
2639 .is_nonzero_shape_function_component[d])
2640 return_value[d] =
2641 fe_values->finite_element_output.shape_gradients
2642 [shape_function_data[shape_function].row_index[d]][q_point];
2643
2644 return symmetrize(return_value);
2645 }
2646 }
2647
2648
2649
2650 template <int dim, int spacedim>
2652 SymmetricTensor<2, dim, spacedim>::value(const unsigned int shape_function,
2653 const unsigned int q_point) const
2654 {
2655 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2656 Assert(fe_values->update_flags & update_values,
2658 "update_values")));
2659
2660 // similar to the vector case where we have more then one index and we need
2661 // to convert between unrolled and component indexing for tensors
2662 const int snc =
2663 shape_function_data[shape_function].single_nonzero_component;
2664
2665 if (snc == -2)
2666 {
2667 // shape function is zero for the selected components
2668 return value_type();
2669 }
2670 else if (snc != -1)
2671 {
2672 value_type return_value;
2673 const unsigned int comp =
2674 shape_function_data[shape_function].single_nonzero_component_index;
2675 return_value[value_type::unrolled_to_component_indices(comp)] =
2676 fe_values->finite_element_output.shape_values(snc, q_point);
2677 return return_value;
2678 }
2679 else
2680 {
2681 value_type return_value;
2682 for (unsigned int d = 0; d < value_type::n_independent_components; ++d)
2683 if (shape_function_data[shape_function]
2684 .is_nonzero_shape_function_component[d])
2685 return_value[value_type::unrolled_to_component_indices(d)] =
2686 fe_values->finite_element_output.shape_values(
2687 shape_function_data[shape_function].row_index[d], q_point);
2688 return return_value;
2689 }
2690 }
2691
2692
2693
2694 template <int dim, int spacedim>
2697 const unsigned int shape_function,
2698 const unsigned int q_point) const
2699 {
2700 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2701 Assert(fe_values->update_flags & update_gradients,
2703 "update_gradients")));
2704
2705 const int snc =
2706 shape_function_data[shape_function].single_nonzero_component;
2707
2708 if (snc == -2)
2709 {
2710 // shape function is zero for the selected components
2711 return divergence_type();
2712 }
2713 else if (snc != -1)
2714 {
2715 // we have a single non-zero component when the symmetric tensor is
2716 // represented in unrolled form. this implies we potentially have
2717 // two non-zero components when represented in component form! we
2718 // will only have one non-zero entry if the non-zero component lies on
2719 // the diagonal of the tensor.
2720 //
2721 // the divergence of a second-order tensor is a first order tensor.
2722 //
2723 // assume the second-order tensor is A with components A_{ij}. then
2724 // A_{ij} = A_{ji} and there is only one (if diagonal) or two non-zero
2725 // entries in the tensorial representation. define the
2726 // divergence as:
2727 // b_i \dealcoloneq \dfrac{\partial phi_{ij}}{\partial x_j}.
2728 // (which is incidentally also
2729 // b_j \dealcoloneq \dfrac{\partial phi_{ij}}{\partial x_i}).
2730 // In both cases, a sum is implied.
2731 //
2732 // Now, we know the nonzero component in unrolled form: it is indicated
2733 // by 'snc'. we can figure out which tensor components belong to this:
2734 const unsigned int comp =
2735 shape_function_data[shape_function].single_nonzero_component_index;
2736 const unsigned int ii =
2737 value_type::unrolled_to_component_indices(comp)[0];
2738 const unsigned int jj =
2739 value_type::unrolled_to_component_indices(comp)[1];
2740
2741 // given the form of the divergence above, if ii=jj there is only a
2742 // single nonzero component of the full tensor and the gradient
2743 // equals
2744 // b_ii \dealcoloneq \dfrac{\partial phi_{ii,ii}}{\partial x_ii}.
2745 // all other entries of 'b' are zero
2746 //
2747 // on the other hand, if ii!=jj, then there are two nonzero entries in
2748 // the full tensor and
2749 // b_ii \dealcoloneq \dfrac{\partial phi_{ii,jj}}{\partial x_ii}.
2750 // b_jj \dealcoloneq \dfrac{\partial phi_{ii,jj}}{\partial x_jj}.
2751 // again, all other entries of 'b' are zero
2752 const ::Tensor<1, spacedim> &phi_grad =
2753 fe_values->finite_element_output.shape_gradients[snc][q_point];
2754
2755 divergence_type return_value;
2756 return_value[ii] = phi_grad[jj];
2757
2758 if (ii != jj)
2759 return_value[jj] = phi_grad[ii];
2760
2761 return return_value;
2762 }
2763 else
2764 {
2766 divergence_type return_value;
2767 return return_value;
2768 }
2769 }
2770
2771
2772
2773 template <int dim, int spacedim>
2775 Tensor<2, dim, spacedim>::value(const unsigned int shape_function,
2776 const unsigned int q_point) const
2777 {
2778 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2779 Assert(fe_values->update_flags & update_values,
2781 "update_values")));
2782
2783 // similar to the vector case where we have more then one index and we need
2784 // to convert between unrolled and component indexing for tensors
2785 const int snc =
2786 shape_function_data[shape_function].single_nonzero_component;
2787
2788 if (snc == -2)
2789 {
2790 // shape function is zero for the selected components
2791 return value_type();
2792 }
2793 else if (snc != -1)
2794 {
2795 value_type return_value;
2796 const unsigned int comp =
2797 shape_function_data[shape_function].single_nonzero_component_index;
2798 const TableIndices<2> indices =
2800 return_value[indices] =
2801 fe_values->finite_element_output.shape_values(snc, q_point);
2802 return return_value;
2803 }
2804 else
2805 {
2806 value_type return_value;
2807 for (unsigned int d = 0; d < dim * dim; ++d)
2808 if (shape_function_data[shape_function]
2809 .is_nonzero_shape_function_component[d])
2810 {
2811 const TableIndices<2> indices =
2813 return_value[indices] =
2814 fe_values->finite_element_output.shape_values(
2815 shape_function_data[shape_function].row_index[d], q_point);
2816 }
2817 return return_value;
2818 }
2819 }
2820
2821
2822
2823 template <int dim, int spacedim>
2825 Tensor<2, dim, spacedim>::divergence(const unsigned int shape_function,
2826 const unsigned int q_point) const
2827 {
2828 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2829 Assert(fe_values->update_flags & update_gradients,
2831 "update_gradients")));
2832
2833 const int snc =
2834 shape_function_data[shape_function].single_nonzero_component;
2835
2836 if (snc == -2)
2837 {
2838 // shape function is zero for the selected components
2839 return divergence_type();
2840 }
2841 else if (snc != -1)
2842 {
2843 // we have a single non-zero component when the tensor is
2844 // represented in unrolled form.
2845 //
2846 // the divergence of a second-order tensor is a first order tensor.
2847 //
2848 // assume the second-order tensor is A with components A_{ij},
2849 // then divergence is d_i := \frac{\partial A_{ij}}{\partial x_j}
2850 //
2851 // Now, we know the nonzero component in unrolled form: it is indicated
2852 // by 'snc'. we can figure out which tensor components belong to this:
2853 const unsigned int comp =
2854 shape_function_data[shape_function].single_nonzero_component_index;
2855 const TableIndices<2> indices =
2857 const unsigned int ii = indices[0];
2858 const unsigned int jj = indices[1];
2859
2860 const ::Tensor<1, spacedim> &phi_grad =
2861 fe_values->finite_element_output.shape_gradients[snc][q_point];
2862
2863 divergence_type return_value;
2864 // note that we contract \nabla from the right
2865 return_value[ii] = phi_grad[jj];
2866
2867 return return_value;
2868 }
2869 else
2870 {
2872 divergence_type return_value;
2873 return return_value;
2874 }
2875 }
2876
2877
2878
2879 template <int dim, int spacedim>
2881 Tensor<2, dim, spacedim>::gradient(const unsigned int shape_function,
2882 const unsigned int q_point) const
2883 {
2884 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell());
2885 Assert(fe_values->update_flags & update_gradients,
2887 "update_gradients")));
2888
2889 const int snc =
2890 shape_function_data[shape_function].single_nonzero_component;
2891
2892 if (snc == -2)
2893 {
2894 // shape function is zero for the selected components
2895 return gradient_type();
2896 }
2897 else if (snc != -1)
2898 {
2899 // we have a single non-zero component when the tensor is
2900 // represented in unrolled form.
2901 //
2902 // the gradient of a second-order tensor is a third order tensor.
2903 //
2904 // assume the second-order tensor is A with components A_{ij},
2905 // then gradient is B_{ijk} := \frac{\partial A_{ij}}{\partial x_k}
2906 //
2907 // Now, we know the nonzero component in unrolled form: it is indicated
2908 // by 'snc'. we can figure out which tensor components belong to this:
2909 const unsigned int comp =
2910 shape_function_data[shape_function].single_nonzero_component_index;
2911 const TableIndices<2> indices =
2913 const unsigned int ii = indices[0];
2914 const unsigned int jj = indices[1];
2915
2916 const ::Tensor<1, spacedim> &phi_grad =
2917 fe_values->finite_element_output.shape_gradients[snc][q_point];
2918
2919 gradient_type return_value;
2920 return_value[ii][jj] = phi_grad;
2921
2922 return return_value;
2923 }
2924 else
2925 {
2927 gradient_type return_value;
2928 return return_value;
2929 }
2930 }
2931} // namespace FEValuesViews
2932
2933#endif
2934
2936
2937#endif
void get_function_third_derivatives_from_local_dof_values(const InputVector &dof_values, std::vector< solution_third_derivative_type< typename InputVector::value_type > > &third_derivatives) const
Scalar & operator=(const Scalar< dim, spacedim > &)=delete
void get_function_laplacians(const ReadVector< Number > &fe_function, std::vector< solution_laplacian_type< Number > > &laplacians) const
value_type value(const unsigned int shape_function, const unsigned int q_point) const
typename ProductType< Number, value_type >::type solution_laplacian_type
void get_function_hessians_from_local_dof_values(const InputVector &dof_values, std::vector< solution_hessian_type< typename InputVector::value_type > > &hessians) const
void get_function_values_from_local_dof_values(const InputVector &dof_values, std::vector< solution_value_type< typename InputVector::value_type > > &values) const
gradient_type gradient(const unsigned int shape_function, const unsigned int q_point) const
std::vector< ShapeFunctionData > shape_function_data
void get_function_gradients(const ReadVector< Number > &fe_function, std::vector< solution_gradient_type< Number > > &gradients) const
Scalar & operator=(Scalar< dim, spacedim > &&) noexcept=default
third_derivative_type third_derivative(const unsigned int shape_function, const unsigned int q_point) const
typename ProductType< Number, gradient_type >::type solution_gradient_type
typename ProductType< Number, value_type >::type solution_value_type
void get_function_values(const ReadVector< Number > &fe_function, std::vector< solution_value_type< Number > > &values) const
void get_function_laplacians_from_local_dof_values(const InputVector &dof_values, std::vector< solution_laplacian_type< typename InputVector::value_type > > &laplacians) const
SmartPointer< const FEValuesBase< dim, spacedim > > fe_values
void get_function_gradients_from_local_dof_values(const InputVector &dof_values, std::vector< solution_gradient_type< typename InputVector::value_type > > &gradients) const
Scalar(Scalar< dim, spacedim > &&)=default
typename ProductType< Number, hessian_type >::type solution_hessian_type
Scalar(const Scalar< dim, spacedim > &)=delete
hessian_type hessian(const unsigned int shape_function, const unsigned int q_point) const
void get_function_third_derivatives(const ReadVector< Number > &fe_function, std::vector< solution_third_derivative_type< Number > > &third_derivatives) const
typename ProductType< Number, third_derivative_type >::type solution_third_derivative_type
void get_function_hessians(const ReadVector< Number > &fe_function, std::vector< solution_hessian_type< Number > > &hessians) const
SymmetricTensor(const SymmetricTensor< 2, dim, spacedim > &)=delete
SymmetricTensor(SymmetricTensor< 2, dim, spacedim > &&)=default
SymmetricTensor & operator=(SymmetricTensor< 2, dim, spacedim > &&) noexcept=default
typename ProductType< Number, value_type >::type solution_value_type
SymmetricTensor & operator=(const SymmetricTensor< 2, dim, spacedim > &)=delete
typename ProductType< Number, divergence_type >::type solution_divergence_type
typename ProductType< Number, value_type >::type solution_value_type
divergence_type divergence(const unsigned int shape_function, const unsigned int q_point) const
typename ProductType< Number, divergence_type >::type solution_divergence_type
typename ProductType< Number, gradient_type >::type solution_gradient_type
Tensor(const Tensor< 2, dim, spacedim > &)=delete
SmartPointer< const FEValuesBase< dim, spacedim > > fe_values
Tensor & operator=(const Tensor< 2, dim, spacedim > &)=delete
Tensor & operator=(Tensor< 2, dim, spacedim > &&)=default
value_type value(const unsigned int shape_function, const unsigned int q_point) const
gradient_type gradient(const unsigned int shape_function, const unsigned int q_point) const
Tensor(Tensor< 2, dim, spacedim > &&)=default
std::vector< ShapeFunctionData > shape_function_data
Vector(const Vector< dim, spacedim > &)=delete
unsigned int first_vector_component
hessian_type hessian(const unsigned int shape_function, const unsigned int q_point) const
typename ProductType< Number, value_type >::type solution_laplacian_type
typename ProductType< Number, third_derivative_type >::type solution_third_derivative_type
Vector & operator=(Vector< dim, spacedim > &&)=default
typename ProductType< Number, hessian_type >::type solution_hessian_type
typename ProductType< Number, divergence_type >::type solution_divergence_type
gradient_type gradient(const unsigned int shape_function, const unsigned int q_point) const
divergence_type divergence(const unsigned int shape_function, const unsigned int q_point) const
Vector & operator=(const Vector< dim, spacedim > &)=delete
third_derivative_type third_derivative(const unsigned int shape_function, const unsigned int q_point) const
typename ProductType< Number, value_type >::type solution_value_type
symmetric_gradient_type symmetric_gradient(const unsigned int shape_function, const unsigned int q_point) const
Vector(Vector< dim, spacedim > &&)=default
typename ProductType< Number, gradient_type >::type solution_gradient_type
typename ::internal::CurlType< spacedim >::type curl_type
typename ProductType< Number, curl_type >::type solution_curl_type
typename ProductType< Number, symmetric_gradient_type >::type solution_symmetric_gradient_type
std::vector< ShapeFunctionData > shape_function_data
value_type value(const unsigned int shape_function, const unsigned int q_point) const
curl_type curl(const unsigned int shape_function, const unsigned int q_point) const
SmartPointer< const FEValuesBase< dim, spacedim > > fe_values
Definition lazy.h:134
std::conditional_t< rank_==1, Number, Tensor< rank_ - 1, dim, Number > > value_type
Definition tensor.h:514
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Number value_type
Definition vector.h:138
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_3rd_derivatives
Third derivatives of shape functions.
@ update_gradients
Shape function gradients.
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
typename ::internal::FEValuesViews:: ViewType< dim, spacedim, Extractor >::type View
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
STL namespace.
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
std::vector< Lazy<::FEValuesViews::Vector< dim, spacedim > > > vectors
std::vector< Lazy<::FEValuesViews::SymmetricTensor< 2, dim, spacedim > > > symmetric_second_order_tensors
std::vector< Lazy<::FEValuesViews::Tensor< 2, dim, spacedim > > > second_order_tensors
std::vector< Lazy<::FEValuesViews::Scalar< dim, spacedim > > > scalars
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)