15#ifndef dealii_fe_point_evaluation_h
16#define dealii_fe_point_evaluation_h
47 <<
"You are requesting information from an FEPointEvaluationBase "
48 <<
"object for which this kind of information has not been computed. "
49 <<
"What information these objects compute is determined by the update_* "
50 <<
"flags you pass to MappingInfo() in the Constructor. Here, "
51 <<
"the operation you are attempting requires the <" << arg1
52 <<
"> flag to be set, but it was apparently not specified "
53 <<
"upon initialization.");
63 typename Enable =
void>
69 typename ::internal::VectorizedArrayTrait<
78 n_components == spacedim,
90 const unsigned int component,
94 result[component] = vector_entry;
108 for (
unsigned int c = 0; c < n_components; ++c)
109 result_scalar[c] = result[c].sum();
111 return result_scalar;
119 return result[component].
sum();
124 const unsigned int vector_lane,
127 for (
unsigned int i = 0; i < n_components; ++i)
128 for (
unsigned int d = 0; d < dim; ++d)
131 value[d][i], vector_lane);
136 const unsigned int vector_lane,
139 for (
unsigned int i = 0; i < n_components; ++i)
140 for (
unsigned int d = 0; d < dim; ++d)
142 value[d][i], vector_lane) = result[i][d];
147 const unsigned int vector_lane,
150 for (
unsigned int i = 0; i < n_components; ++i)
151 for (
unsigned int d = 0; d < dim; ++d)
153 value[d][i], vector_lane) = result[i][d];
158 const unsigned int vector_lane)
160 for (
unsigned int i = 0; i < n_components; ++i)
161 for (
unsigned int d = 0; d < spacedim; ++d)
168 const unsigned int vector_lane,
171 for (
unsigned int i = 0; i < n_components; ++i)
172 result[i] =
value[i][vector_lane];
185 const unsigned int vector_lane,
188 for (
unsigned int i = 0; i < n_components; ++i)
189 value[i][vector_lane] = result[i];
203 for (
unsigned int i = 0; i < n_components; ++i)
210 const unsigned int vector_lane,
211 const unsigned int component,
215 vector_lane) += shape_value;
220 const unsigned int vector_lane,
221 const unsigned int component)
229 const unsigned int vector_lane,
230 const unsigned int component,
233 for (
unsigned int d = 0; d < spacedim; ++d)
241 const unsigned int vector_lane,
242 const unsigned int component)
245 for (
unsigned int d = 0; d < spacedim; ++d)
253 template <
int dim,
int spacedim,
typename Number>
259 typename ::internal::VectorizedArrayTrait<
276 result = vector_entry;
299 const unsigned int vector_lane,
302 for (
unsigned int d = 0; d < dim; ++d)
303 result[d] =
value[d][vector_lane];
316 const unsigned int vector_lane,
319 for (
unsigned int d = 0; d < dim; ++d)
320 value[d][vector_lane] = result[d];
333 const unsigned int vector_lane)
335 for (
unsigned int d = 0; d < spacedim; ++d)
342 const unsigned int vector_lane,
345 result =
value[vector_lane];
358 const unsigned int vector_lane,
361 value[vector_lane] = result;
380 const unsigned int vector_lane,
390 const unsigned int vector_lane,
398 const unsigned int vector_lane,
402 for (
unsigned int d = 0; d < spacedim; ++d)
409 const unsigned int vector_lane,
413 for (
unsigned int d = 0; d < spacedim; ++d)
420 template <
int dim,
typename Number>
425 std::enable_if_t<dim != 1>>
430 typename ::internal::VectorizedArrayTrait<
444 const unsigned int component,
448 result[component] = vector_entry;
462 for (
unsigned int c = 0; c < dim; ++c)
463 result_scalar[c] = result[c].sum();
465 return result_scalar;
473 return result[component].
sum();
478 const unsigned int vector_lane,
481 for (
unsigned int i = 0; i < dim; ++i)
482 for (
unsigned int d = 0; d < dim; ++d)
485 value[d][i], vector_lane);
490 const unsigned int vector_lane,
493 for (
unsigned int i = 0; i < dim; ++i)
494 for (
unsigned int d = 0; d < dim; ++d)
496 value[d][i], vector_lane) = result[i][d];
501 const unsigned int vector_lane)
503 for (
unsigned int i = 0; i < dim; ++i)
504 for (
unsigned int d = 0; d < dim; ++d)
511 const unsigned int vector_lane,
514 for (
unsigned int i = 0; i < dim; ++i)
515 result[i] =
value[i][vector_lane];
528 const unsigned int vector_lane,
531 for (
unsigned int i = 0; i < dim; ++i)
532 value[i][vector_lane] = result[i];
546 for (
unsigned int i = 0; i < dim; ++i)
553 const unsigned int vector_lane,
554 const unsigned int component,
558 vector_lane) += shape_value;
563 const unsigned int vector_lane,
564 const unsigned int component)
572 const unsigned int vector_lane,
573 const unsigned int component,
576 for (
unsigned int d = 0; d < dim; ++d)
584 const unsigned int vector_lane,
585 const unsigned int component)
588 for (
unsigned int d = 0; d < dim; ++d)
596 template <
int dim,
int spacedim>
599 const unsigned int base_element_number);
601 template <
int dim,
int spacedim>
605 template <
int dim,
int spacedim>
606 std::vector<Polynomials::Polynomial<double>>
619template <
int n_components_,
622 typename Number =
double>
635 using ETT =
typename internal::FEPointEvaluation::
636 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
642 typename ETT::interface_vectorized_unit_gradient_type;
666 const unsigned int first_selected_component = 0);
690 const unsigned int first_selected_component = 0,
776 JxW(
const unsigned int point_index)
const;
836 setup(
const unsigned int first_selected_component);
843 template <
bool is_face,
bool is_linear>
876 std::vector<Polynomials::Polynomial<double>>
poly;
1012 std::unique_ptr<NonMatching::MappingInfo<dim, spacedim, Number>>
1094template <
int n_components_,
1097 typename Number =
double>
1111 using ETT =
typename internal::FEPointEvaluation::
1112 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1119 typename ETT::interface_vectorized_unit_gradient_type;
1142 const unsigned int first_selected_component = 0);
1163 const unsigned int first_selected_component = 0);
1206 template <std::
size_t str
ide_view>
1249 template <std::
size_t str
ide_view>
1253 const bool sum_into_values =
false);
1280 const bool sum_into_values =
false);
1308 template <std::
size_t str
ide_view>
1313 const bool sum_into_values =
false);
1344 const bool sum_into_values =
false);
1366 template <
bool is_linear, std::
size_t str
ide_view>
1375 template <
bool is_linear, std::
size_t str
ide_view>
1380 const unsigned int n_shapes,
1381 const unsigned int qb,
1388 template <
bool is_linear, std::
size_t str
ide_view>
1397 template <std::
size_t str
ide_view>
1408 template <
bool is_linear>
1412 const unsigned int n_shapes,
1413 const unsigned int qb,
1423 template <
bool is_linear, std::
size_t str
ide_view>
1428 const bool sum_into_values);
1433 template <
bool do_JxW,
bool is_linear, std::
size_t str
ide_view>
1438 const bool sum_into_values);
1443 template <
bool do_JxW, std::
size_t str
ide_view>
1448 const bool sum_into_values);
1453 template <
bool do_JxW, std::
size_t str
ide_view>
1458 const bool sum_into_values);
1480template <
int n_components_,
1483 typename Number =
double>
1497 using ETT =
typename internal::FEPointEvaluation::
1498 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1505 typename ETT::interface_vectorized_unit_gradient_type;
1514 const unsigned int first_selected_component = 0);
1528 reinit(
const unsigned int face_index);
1541 template <std::
size_t str
ide_view>
1584 template <std::
size_t str
ide_view>
1588 const bool sum_into_values =
false);
1615 const bool sum_into_values =
false);
1639 template <std::
size_t str
ide_view>
1644 const bool sum_into_values =
false);
1671 const bool sum_into_values =
false);
1679 template <
int str
ide_face_dof = VectorizedArrayType::size()>
1690 template <
int str
ide_face_dof = VectorizedArrayType::size()>
1694 const bool sum_into_values =
false);
1712 template <
bool is_linear, std::
size_t str
ide_view>
1718 template <
bool do_JxW,
bool is_linear, std::
size_t str
ide_view>
1723 const bool sum_into_values);
1729 template <
bool is_linear,
int str
ide_face_dof>
1738 template <
bool do_JxW,
bool is_linear,
int str
ide_face_dof>
1743 const bool sum_into_values);
1751template <
int n_components_,
int dim,
int spacedim,
typename Number>
1756 const unsigned int first_selected_component)
1757 : n_q_batches(
numbers::invalid_unsigned_int)
1758 , n_q_points(
numbers::invalid_unsigned_int)
1759 , n_q_points_scalar(
numbers::invalid_unsigned_int)
1763 , update_flags(update_flags)
1764 , mapping_info_on_the_fly(
1765 std::make_unique<
NonMatching::MappingInfo<dim, spacedim, Number>>(
1768 , mapping_info(mapping_info_on_the_fly.get())
1769 , current_cell_index(
numbers::invalid_unsigned_int)
1770 , current_face_number(
numbers::invalid_unsigned_int)
1771 , is_reinitialized(false)
1774 setup(first_selected_component);
1779template <
int n_components_,
int dim,
int spacedim,
typename Number>
1784 const unsigned int first_selected_component,
1785 const bool is_interior)
1786 : n_q_batches(
numbers::invalid_unsigned_int)
1787 , n_q_points(
numbers::invalid_unsigned_int)
1788 , n_q_points_scalar(
numbers::invalid_unsigned_int)
1789 , mapping(&mapping_info.get_mapping())
1792 , update_flags(mapping_info.get_update_flags())
1793 , mapping_info(&mapping_info)
1794 , current_cell_index(
numbers::invalid_unsigned_int)
1795 , current_face_number(
numbers::invalid_unsigned_int)
1796 , is_reinitialized(false)
1797 , is_interior(is_interior)
1799 setup(first_selected_component);
1806template <
int n_components_,
int dim,
int spacedim,
typename Number>
1810 : n_q_batches(other.n_q_batches)
1811 , n_q_points(other.n_q_points)
1812 , n_q_points_scalar(other.n_q_points_scalar)
1813 , mapping(other.mapping)
1816 , use_linear_path(other.use_linear_path)
1817 , renumber(other.renumber)
1818 , solution_renumbered(other.solution_renumbered)
1819 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1820 , values(other.values)
1821 , gradients(other.gradients)
1822 , dofs_per_component(other.dofs_per_component)
1823 , dofs_per_component_face(other.dofs_per_component_face)
1824 , component_in_base_element(other.component_in_base_element)
1825 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1826 , update_flags(other.update_flags)
1827 , fe_values(other.fe_values)
1828 , mapping_info_on_the_fly(
1829 other.mapping_info_on_the_fly ?
1834 , mapping_info(other.mapping_info)
1835 , current_cell_index(other.current_cell_index)
1836 , current_face_number(other.current_face_number)
1837 , fast_path(other.fast_path)
1838 , is_reinitialized(false)
1839 , shapes(other.shapes)
1840 , shapes_faces(other.shapes_faces)
1841 , is_interior(other.is_interior)
1843 connection_is_reinitialized = mapping_info->connect_is_reinitialized(
1844 [
this]() { this->is_reinitialized =
false; });
1849template <
int n_components_,
int dim,
int spacedim,
typename Number>
1852 : n_q_batches(other.n_q_batches)
1853 , n_q_points(other.n_q_points)
1854 , n_q_points_scalar(other.n_q_points_scalar)
1855 , mapping(other.mapping)
1858 , use_linear_path(other.use_linear_path)
1859 , renumber(other.renumber)
1860 , solution_renumbered(other.solution_renumbered)
1861 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1862 , values(other.values)
1863 , gradients(other.gradients)
1864 , dofs_per_component(other.dofs_per_component)
1865 , dofs_per_component_face(other.dofs_per_component_face)
1866 , component_in_base_element(other.component_in_base_element)
1867 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1868 , update_flags(other.update_flags)
1869 , fe_values(other.fe_values)
1870 , mapping_info_on_the_fly(std::move(other.mapping_info_on_the_fly))
1871 , mapping_info(other.mapping_info)
1872 , current_cell_index(other.current_cell_index)
1873 , current_face_number(other.current_face_number)
1874 , fast_path(other.fast_path)
1875 , is_reinitialized(
false)
1876 , shapes(other.shapes)
1877 , shapes_faces(other.shapes_faces)
1878 , is_interior(other.is_interior)
1880 connection_is_reinitialized = mapping_info->connect_is_reinitialized(
1881 [
this]() { this->is_reinitialized =
false; });
1886template <
int n_components_,
int dim,
int spacedim,
typename Number>
1890 connection_is_reinitialized.disconnect();
1895template <
int n_components_,
int dim,
int spacedim,
typename Number>
1898 const unsigned int first_selected_component)
1901 fe->n_components() + 1);
1903 shapes.reserve(100);
1905 bool same_base_element =
true;
1906 unsigned int base_element_number = 0;
1907 component_in_base_element = 0;
1908 unsigned int component = 0;
1909 for (; base_element_number < fe->n_base_elements(); ++base_element_number)
1910 if (component + fe->element_multiplicity(base_element_number) >
1911 first_selected_component)
1913 if (first_selected_component + n_components >
1914 component + fe->element_multiplicity(base_element_number))
1915 same_base_element =
false;
1916 component_in_base_element = first_selected_component - component;
1920 component += fe->element_multiplicity(base_element_number);
1924 *fe, base_element_number) &&
1927 shape_info.reinit(
QMidpoint<1>(), *fe, base_element_number);
1928 renumber = shape_info.lexicographic_numbering;
1929 dofs_per_component = shape_info.dofs_per_component_on_cell;
1930 dofs_per_component_face = shape_info.dofs_per_component_on_face;
1932 fe->base_element(base_element_number));
1934 bool is_lexicographic =
true;
1935 for (
unsigned int i = 0; i < renumber.size(); ++i)
1936 if (i != renumber[i])
1937 is_lexicographic =
false;
1939 if (is_lexicographic)
1942 use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. &&
1943 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1944 poly[1].value(1.) == 1.) &&
1945 (fe->n_components() == n_components);
1947 const unsigned int size_face = 3 * dofs_per_component_face * n_components;
1948 const unsigned int size_cell = dofs_per_component * n_components;
1949 scratch_data_scalar.resize(size_face + size_cell);
1951 solution_renumbered.resize(dofs_per_component);
1952 solution_renumbered_vectorized.resize(dofs_per_component);
1958 nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
1959 for (
unsigned int d = 0; d < n_components; ++d)
1961 const unsigned int component = first_selected_component + d;
1962 for (
unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1964 const bool is_primitive =
1965 fe->is_primitive() || fe->is_primitive(i);
1967 nonzero_shape_function_component[i][d] =
1968 (component == fe->system_to_component_index(i).first);
1970 nonzero_shape_function_component[i][d] =
1971 (fe->get_nonzero_components(i)[component] ==
true);
1981template <
int n_components_,
int dim,
int spacedim,
typename Number>
1982template <
bool is_face,
bool is_linear>
1986 const unsigned int geometry_index =
1987 mapping_info->template compute_geometry_index_offset<is_face>(
1988 current_cell_index, current_face_number);
1990 cell_type = mapping_info->get_cell_type(geometry_index);
1992 const_cast<unsigned int &
>(n_q_points_scalar) =
1993 mapping_info->get_n_q_points_unvectorized(geometry_index);
1996 const_cast<unsigned int &
>(n_q_batches) =
1997 (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
1999 const unsigned int n_q_points_before = n_q_points;
2001 const_cast<unsigned int &
>(n_q_points) =
2002 (stride == 1) ? n_q_batches : n_q_points_scalar;
2004 if (n_q_points != n_q_points_before)
2007 values.resize(n_q_points);
2009 gradients.resize(n_q_points);
2012 if (n_q_points == 0)
2014 is_reinitialized =
true;
2019 const unsigned int unit_point_offset =
2020 mapping_info->compute_unit_point_index_offset(geometry_index);
2023 unit_point_faces_ptr =
2024 mapping_info->get_unit_point_faces(unit_point_offset);
2026 unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
2029 const unsigned int data_offset =
2030 mapping_info->compute_data_index_offset(geometry_index);
2031 const unsigned int compressed_data_offset =
2032 mapping_info->compute_compressed_data_index_offset(geometry_index);
2035 mapping_info->get_update_flags_mapping();
2037 real_point_ptr = mapping_info->get_real_point(data_offset);
2040 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2042 inverse_jacobian_ptr =
2043 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2045 normal_ptr = mapping_info->get_normal_vector(data_offset);
2047 JxW_ptr = mapping_info->get_JxW(data_offset);
2049 real_point_ptr = mapping_info->get_real_point(data_offset);
2051 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2052 inverse_jacobian_ptr =
2053 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2054 normal_ptr = mapping_info->get_normal_vector(data_offset);
2055 JxW_ptr = mapping_info->get_JxW(data_offset);
2058 if (!is_linear && fast_path)
2060 const std::size_t n_shapes = poly.size();
2062 for (
unsigned int qb = 0; qb < n_q_batches; ++qb)
2067 shapes_faces.resize_fast(n_q_batches * n_shapes);
2069 shapes_faces.data() + qb * n_shapes,
2071 unit_point_faces_ptr[qb],
2077 shapes.resize_fast(n_q_batches * n_shapes);
2079 shapes.data() + qb * n_shapes,
2086 is_reinitialized =
true;
2091template <
int n_components_,
int dim,
int spacedim,
typename Number>
2095 Number>::value_type &
2097 const unsigned int point_index)
const
2100 return values[point_index];
2105template <
int n_components_,
int dim,
int spacedim,
typename Number>
2109 Number>::gradient_type &
2111 const unsigned int point_index)
const
2114 return gradients[point_index];
2119template <
int n_components_,
int dim,
int spacedim,
typename Number>
2123 const unsigned int point_index)
2126 values[point_index] = value;
2131template <
int n_components_,
int dim,
int spacedim,
typename Number>
2135 const unsigned int point_index)
2138 gradients[point_index] = gradient;
2143template <
int n_components_,
int dim,
int spacedim,
typename Number>
2146 const unsigned int point_index)
const
2149 Assert(jacobian_ptr !=
nullptr,
2151 ExcFEPointEvaluationAccessToUninitializedMappingField(
2152 "update_jacobians"));
2153 return jacobian_ptr[cell_type <= ::internal::MatrixFreeFunctions::
2154 GeometryType::affine ?
2161template <
int n_components_,
int dim,
int spacedim,
typename Number>
2164 const unsigned int point_index)
const
2167 Assert(inverse_jacobian_ptr !=
nullptr,
2169 ExcFEPointEvaluationAccessToUninitializedMappingField(
2170 "update_inverse_jacobians"));
2171 return inverse_jacobian_ptr
2180template <
int n_components_,
int dim,
int spacedim,
typename Number>
2183 const unsigned int point_index)
const
2186 Assert(JxW_ptr !=
nullptr,
2188 ExcFEPointEvaluationAccessToUninitializedMappingField(
2189 "update_JxW_values"));
2190 return JxW_ptr[point_index];
2195template <
int n_components_,
int dim,
int spacedim,
typename Number>
2198 const unsigned int point_index)
const
2200 return quadrature_point(point_index);
2206template <
int n_components_,
int dim,
int spacedim,
typename Number>
2209 const unsigned int point_index)
const
2212 Assert(real_point_ptr !=
nullptr,
2214 ExcFEPointEvaluationAccessToUninitializedMappingField(
2215 "update_quadrature_points"));
2216 return real_point_ptr[point_index];
2221template <
int n_components_,
int dim,
int spacedim,
typename Number>
2224 const unsigned int point_index)
const
2227 Assert(unit_point_ptr !=
nullptr,
ExcMessage(
"unit_point_ptr is not set!"));
2229 for (
unsigned int d = 0; d < dim; ++d)
2231 unit_point_ptr[point_index / stride][d], point_index % stride);
2237template <
int n_components_,
int dim,
int spacedim,
typename Number>
2242 return {0U, n_q_points};
2247template <
int n_components_,
int dim,
int spacedim,
typename Number>
2251 const unsigned int first_selected_component)
2255 first_selected_component)
2260template <
int n_components_,
int dim,
int spacedim,
typename Number>
2265 const unsigned int first_selected_component)
2270 first_selected_component)
2275template <
int n_components_,
int dim,
int spacedim,
typename Number>
2282 if (this->use_linear_path)
2283 this->
template do_reinit<false, true>();
2285 this->
template do_reinit<false, false>();
2290template <
int n_components_,
int dim,
int spacedim,
typename Number>
2297 AssertThrow(this->mapping_info_on_the_fly.get() !=
nullptr,
2300 this->mapping_info->reinit(cell, unit_points);
2302 if (!this->fast_path)
2304 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2308 std::vector<
Point<dim>>(unit_points.begin(), unit_points.end())),
2309 this->update_flags);
2310 this->fe_values->reinit(cell);
2313 if (this->use_linear_path)
2314 this->
template do_reinit<false, true>();
2316 this->
template do_reinit<false, false>();
2321template <
int n_components_,
int dim,
int spacedim,
typename Number>
2329 if (this->use_linear_path)
2330 this->
template do_reinit<false, true>();
2332 this->
template do_reinit<false, false>();
2334 if (!this->fast_path)
2336 std::vector<Point<dim>> unit_points(this->n_q_points_scalar);
2338 for (
unsigned int v = 0; v < this->n_q_points_scalar; ++v)
2339 for (
unsigned int d = 0; d < dim; ++d)
2341 this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
2343 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2347 std::vector<
Point<dim>>(unit_points.begin(), unit_points.end())),
2348 this->update_flags);
2350 this->fe_values->reinit(
2351 this->mapping_info->get_cell_iterator(this->current_cell_index));
2357template <
int n_components_,
int dim,
int spacedim,
typename Number>
2358template <std::
size_t str
ide_view>
2364 if (!this->is_reinitialized)
2367 if (this->n_q_points == 0)
2377 if (this->fast_path)
2379 if (this->use_linear_path)
2380 evaluate_fast<true>(solution_values, evaluation_flags);
2382 evaluate_fast<false>(solution_values, evaluation_flags);
2385 evaluate_slow(solution_values, evaluation_flags);
2390template <
int n_components_,
int dim,
int spacedim,
typename Number>
2397 solution_values.
size()),
2403template <
int n_components_,
int dim,
int spacedim,
typename Number>
2404template <std::
size_t str
ide_view>
2409 const bool sum_into_values)
2411 do_integrate<true>(solution_values, integration_flags, sum_into_values);
2416template <
int n_components_,
int dim,
int spacedim,
typename Number>
2421 const bool sum_into_values)
2424 solution_values.
size()),
2431template <
int n_components_,
int dim,
int spacedim,
typename Number>
2439 for (
const auto point_index : this->quadrature_point_indices())
2440 return_value += values[point_index] * this->JxW(point_index);
2442 return ETT::sum_value(return_value);
2447template <
int n_components_,
int dim,
int spacedim,
typename Number>
2448template <std::
size_t str
ide_view>
2453 const bool sum_into_values)
2455 do_integrate<false>(solution_values, integration_flags, sum_into_values);
2460template <
int n_components_,
int dim,
int spacedim,
typename Number>
2465 const bool sum_into_values)
2468 solution_values.
size()),
2475template <
int n_components_,
int dim,
int spacedim,
typename Number>
2476template <
bool is_linear, std::
size_t str
ide_view>
2481 const unsigned int dofs_per_comp =
2484 for (
unsigned int comp = 0; comp < n_components; ++comp)
2486 const std::size_t offset =
2487 (this->component_in_base_element + comp) * dofs_per_comp;
2489 if ((is_linear && n_components == 1) || this->renumber.empty())
2491 for (
unsigned int i = 0; i < dofs_per_comp; ++i)
2492 ETT::read_value(solution_values[i + offset],
2494 this->solution_renumbered[i]);
2498 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2499 for (
unsigned int i = 0; i < dofs_per_comp; ++i)
2500 ETT::read_value(solution_values[renumber_ptr[i]],
2502 this->solution_renumbered[i]);
2509template <
int n_components_,
int dim,
int spacedim,
typename Number>
2510template <
bool is_linear, std::
size_t str
ide_view>
2515 const unsigned int n_shapes,
2516 const unsigned int qb,
2522 std::array<vectorized_value_type, dim + 1> result;
2523 if constexpr (is_linear)
2525 if constexpr (n_components == 1)
2532 stride_view>(solution_values.
data(), this->unit_point_ptr[qb]);
2536 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2544 false>(this->shapes.data() + qb * n_shapes,
2546 this->solution_renumbered.data());
2547 gradient[0] = result[0];
2549 gradient[1] = result[1];
2551 gradient[2] = result[2];
2552 value = result[dim];
2556 if constexpr (is_linear)
2558 if constexpr (n_components == 1)
2563 stride_view>(solution_values.
data(), this->unit_point_ptr[qb]);
2566 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2574 this->shapes.data() + qb * n_shapes,
2576 this->solution_renumbered.data());
2582template <
int n_components_,
int dim,
int spacedim,
typename Number>
2583template <
bool is_linear, std::
size_t str
ide_view>
2589 if (!(is_linear && n_components == 1))
2590 prepare_evaluate_fast<is_linear>(solution_values);
2593 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2595 for (
unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2600 compute_evaluate_fast<is_linear>(
2601 solution_values, evaluation_flags, n_shapes, qb, value, gradient);
2605 for (
unsigned int v = 0, offset = qb * stride;
2606 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2608 ETT::set_value(value, v, this->values[offset]);
2616 for (
unsigned int v = 0, offset = qb * stride;
2617 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2621 ETT::set_gradient(gradient, v, unit_gradient);
2622 this->gradients[offset] =
2626 this->inverse_jacobian_ptr[0].
transpose(), unit_gradient) :
2629 ->inverse_jacobian_ptr[this->cell_type <=
2631 GeometryType::affine ?
2643template <
int n_components_,
int dim,
int spacedim,
typename Number>
2644template <std::
size_t str
ide_view>
2651 Assert(this->fe_values.get() !=
nullptr,
2653 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2655 const std::size_t n_points = this->fe_values->get_quadrature().size();
2659 this->values.resize(this->n_q_points);
2660 std::fill(this->values.begin(), this->values.end(),
value_type());
2661 for (
unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2664 for (
unsigned int d = 0; d < n_components; ++d)
2665 if (this->nonzero_shape_function_component[i][d] &&
2666 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2667 for (
unsigned int qb = 0, q = 0; q < n_points;
2668 ++qb, q += n_lanes_user_interface)
2669 for (
unsigned int v = 0;
2670 v < n_lanes_user_interface && q + v < n_points;
2672 ETT::access(this->values[qb],
2675 this->fe_values->shape_value(i, q + v) * value);
2676 else if (this->nonzero_shape_function_component[i][d])
2677 for (
unsigned int qb = 0, q = 0; q < n_points;
2678 ++qb, q += n_lanes_user_interface)
2679 for (
unsigned int v = 0;
2680 v < n_lanes_user_interface && q + v < n_points;
2682 ETT::access(this->values[qb],
2685 this->fe_values->shape_value_component(i,
2694 this->gradients.resize(this->n_q_points);
2695 std::fill(this->gradients.begin(),
2696 this->gradients.end(),
2698 for (
unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2701 for (
unsigned int d = 0; d < n_components; ++d)
2702 if (this->nonzero_shape_function_component[i][d] &&
2703 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2704 for (
unsigned int qb = 0, q = 0; q < n_points;
2705 ++qb, q += n_lanes_user_interface)
2706 for (
unsigned int v = 0;
2707 v < n_lanes_user_interface && q + v < n_points;
2709 ETT::access(this->gradients[qb],
2712 this->fe_values->shape_grad(i, q + v) * value);
2713 else if (this->nonzero_shape_function_component[i][d])
2714 for (
unsigned int qb = 0, q = 0; q < n_points;
2715 ++qb, q += n_lanes_user_interface)
2716 for (
unsigned int v = 0;
2717 v < n_lanes_user_interface && q + v < n_points;
2720 this->gradients[qb],
2723 this->fe_values->shape_grad_component(i, q + v, d) * value);
2730template <
int n_components_,
int dim,
int spacedim,
typename Number>
2731template <
bool is_linear>
2735 const unsigned int n_shapes,
2736 const unsigned int qb,
2751 solution_values_vectorized_linear :
2752 this->solution_renumbered_vectorized.data(),
2753 this->unit_point_ptr[qb],
2760 this->shapes.data() + qb * n_shapes,
2763 is_linear ? solution_values_vectorized_linear :
2764 this->solution_renumbered_vectorized.data(),
2765 this->unit_point_ptr[qb],
2771template <
int n_components_,
int dim,
int spacedim,
typename Number>
2772template <
bool is_linear, std::
size_t str
ide_view>
2777 const bool sum_into_values)
2779 if (!sum_into_values && this->fe->n_components() > n_components)
2780 for (
unsigned int i = 0; i < solution_values.
size(); ++i)
2781 solution_values[i] = 0;
2783 const unsigned int dofs_per_comp =
2786 for (
unsigned int comp = 0; comp < n_components; ++comp)
2788 const std::size_t offset =
2789 (this->component_in_base_element + comp) * dofs_per_comp;
2791 if (is_linear || this->renumber.empty())
2793 for (
unsigned int i = 0; i < dofs_per_comp; ++i)
2794 if (sum_into_values)
2795 solution_values[i + offset] +=
2796 ETT::sum_value(comp,
2798 *(solution_values_vectorized_linear + i) :
2799 this->solution_renumbered_vectorized[i]);
2801 solution_values[i + offset] =
2802 ETT::sum_value(comp,
2804 *(solution_values_vectorized_linear + i) :
2805 this->solution_renumbered_vectorized[i]);
2809 const unsigned int *renumber_ptr = this->renumber.
data() + offset;
2810 for (
unsigned int i = 0; i < dofs_per_comp; ++i)
2811 if (sum_into_values)
2812 solution_values[renumber_ptr[i]] +=
2813 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2815 solution_values[renumber_ptr[i]] =
2816 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2823template <
int n_components_,
int dim,
int spacedim,
typename Number>
2824template <
bool do_JxW,
bool is_linear, std::
size_t str
ide_view>
2829 const bool sum_into_values)
2832 if constexpr (stride == 1)
2833 if (
const unsigned int n_filled_lanes =
2834 this->n_q_points_scalar & (n_lanes_internal - 1);
2838 for (
unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2839 ETT::set_zero_value(this->values.back(), v);
2841 for (
unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2842 ETT::set_zero_gradient(this->gradients.back(), v);
2846 solution_values_vectorized_linear = {};
2849 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2851 const bool cartesian_cell =
2853 const bool affine_cell =
2855 for (
unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2861 for (
unsigned int v = 0, offset = qb * stride;
2862 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2864 ETT::get_value(value,
2866 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
2867 this->values[offset]);
2870 for (
unsigned int v = 0, offset = qb * stride;
2871 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2875 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
2876 this->gradients[offset];
2884 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
2888 compute_integrate_fast<is_linear>(
2894 solution_values_vectorized_linear.data());
2898 finish_integrate_fast<is_linear>(solution_values,
2899 solution_values_vectorized_linear.
data(),
2905template <
int n_components_,
int dim,
int spacedim,
typename Number>
2906template <
bool do_JxW, std::
size_t str
ide_view>
2911 const bool sum_into_values)
2914 Assert(this->fe_values.get() !=
nullptr,
2916 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2917 if (!sum_into_values)
2918 for (
unsigned int i = 0; i < solution_values.
size(); ++i)
2919 solution_values[i] = 0;
2921 const std::size_t n_points = this->fe_values->get_quadrature().
size();
2926 for (
unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2928 for (
unsigned int d = 0; d < n_components; ++d)
2929 if (this->nonzero_shape_function_component[i][d] &&
2930 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2931 for (
unsigned int qb = 0, q = 0; q < n_points;
2932 ++qb, q += n_lanes_user_interface)
2933 for (
unsigned int v = 0;
2934 v < n_lanes_user_interface && q + v < n_points;
2936 solution_values[i] +=
2937 this->fe_values->shape_value(i, q + v) *
2938 ETT::access(this->values[qb], v, d) *
2939 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2940 else if (this->nonzero_shape_function_component[i][d])
2941 for (
unsigned int qb = 0, q = 0; q < n_points;
2942 ++qb, q += n_lanes_user_interface)
2943 for (
unsigned int v = 0;
2944 v < n_lanes_user_interface && q + v < n_points;
2946 solution_values[i] +=
2947 this->fe_values->shape_value_component(i, q + v, d) *
2948 ETT::access(this->values[qb], v, d) *
2949 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2956 for (
unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2958 for (
unsigned int d = 0; d < n_components; ++d)
2959 if (this->nonzero_shape_function_component[i][d] &&
2960 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2961 for (
unsigned int qb = 0, q = 0; q < n_points;
2962 ++qb, q += n_lanes_user_interface)
2963 for (
unsigned int v = 0;
2964 v < n_lanes_user_interface && q + v < n_points;
2966 solution_values[i] +=
2967 this->fe_values->shape_grad(i, q + v) *
2968 ETT::access(this->gradients[qb], v, d) *
2969 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2970 else if (this->nonzero_shape_function_component[i][d])
2971 for (
unsigned int qb = 0, q = 0; q < n_points;
2972 ++qb, q += n_lanes_user_interface)
2973 for (
unsigned int v = 0;
2974 v < n_lanes_user_interface && q + v < n_points;
2976 solution_values[i] +=
2977 this->fe_values->shape_grad_component(i, q + v, d) *
2978 ETT::access(this->gradients[qb], v, d) *
2979 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2986template <
int n_components_,
int dim,
int spacedim,
typename Number>
2987template <
bool do_JxW, std::
size_t str
ide_view>
2992 const bool sum_into_values)
2994 if (!this->is_reinitialized)
2999 if (this->n_q_points == 0 ||
3001 (integration_flags &
3004 if (!sum_into_values)
3005 for (
unsigned int i = 0; i < solution_values.
size(); ++i)
3006 solution_values[i] = 0;
3011 !do_JxW || this->JxW_ptr !=
nullptr,
3013 "JxW pointer is not set! If you do not want to integrate() use test_and_sum()"));
3016 if (this->fast_path)
3018 if (this->use_linear_path)
3019 integrate_fast<do_JxW, true>(solution_values,
3023 integrate_fast<do_JxW, false>(solution_values,
3028 integrate_slow<do_JxW>(solution_values, integration_flags, sum_into_values);
3033template <
int n_components_,
int dim,
int spacedim,
typename Number>
3036 const unsigned int point_index)
const
3039 Assert(this->normal_ptr !=
nullptr,
3041 ExcFEPointEvaluationAccessToUninitializedMappingField(
3042 "update_normal_vectors"));
3043 if (this->is_interior)
3044 return this->normal_ptr[point_index];
3046 return -this->normal_ptr[point_index];
3051template <
int n_components_,
int dim,
int spacedim,
typename Number>
3056 const bool is_interior,
3057 const unsigned int first_selected_component)
3061 first_selected_component,
3067template <
int n_components_,
int dim,
int spacedim,
typename Number>
3071 const unsigned int face_number)
3074 this->current_face_number = face_number;
3076 if (this->use_linear_path)
3077 this->
template do_reinit<true, true>();
3079 this->
template do_reinit<true, false>();
3084template <
int n_components_,
int dim,
int spacedim,
typename Number>
3087 const unsigned int face_index)
3089 this->current_cell_index = face_index;
3090 this->current_face_number =
3091 this->mapping_info->get_face_number(face_index, this->is_interior);
3093 if (this->use_linear_path)
3094 this->
template do_reinit<true, true>();
3096 this->
template do_reinit<true, false>();
3101template <
int n_components_,
int dim,
int spacedim,
typename Number>
3102template <std::
size_t str
ide_view>
3110 if (this->n_q_points == 0)
3121 if (this->use_linear_path)
3122 do_evaluate<true>(solution_values, evaluation_flags);
3124 do_evaluate<false>(solution_values, evaluation_flags);
3129template <
int n_components_,
int dim,
int spacedim,
typename Number>
3136 solution_values.
size()),
3142template <
int n_components_,
int dim,
int spacedim,
typename Number>
3143template <
bool is_linear, std::
size_t str
ide_view>
3149 const unsigned int dofs_per_comp =
3153 if (stride_view == 1 && this->component_in_base_element == 0 &&
3154 (is_linear || this->renumber.empty()))
3155 input = solution_values.
data();
3158 for (
unsigned int comp = 0; comp < n_components; ++comp)
3160 const std::size_t offset =
3161 (this->component_in_base_element + comp) * dofs_per_comp;
3163 if (is_linear || this->renumber.empty())
3165 for (
unsigned int i = 0; i < dofs_per_comp; ++i)
3166 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3167 solution_values[i + offset];
3171 const unsigned int *renumber_ptr = this->renumber.
data() + offset;
3172 for (
unsigned int i = 0; i < dofs_per_comp; ++i)
3173 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3174 solution_values[renumber_ptr[i]];
3177 input = this->scratch_data_scalar.
data();
3181 this->scratch_data_scalar.begin() + dofs_per_comp * n_components;
3184 template interpolate<true, false>(n_components,
3189 this->current_face_number);
3191 do_evaluate_in_face<is_linear, 1>(output, evaluation_flags);
3196template <
int n_components_,
int dim,
int spacedim,
typename Number>
3197template <std::
size_t str
ide_view>
3202 const bool sum_into_values)
3208 if (this->n_q_points == 0 ||
3210 (integration_flags &
3213 if (!sum_into_values)
3214 for (
unsigned int i = 0; i < solution_values.
size(); ++i)
3215 solution_values[i] = 0;
3221 if (this->use_linear_path)
3222 do_integrate<true, true>(solution_values,
3226 do_integrate<true, false>(solution_values,
3233template <
int n_components_,
int dim,
int spacedim,
typename Number>
3238 const bool sum_into_values)
3241 solution_values.
size()),
3248template <
int n_components_,
int dim,
int spacedim,
typename Number>
3249template <std::
size_t str
ide_view>
3254 const bool sum_into_values)
3260 if (this->n_q_points == 0 ||
3262 (integration_flags &
3265 if (!sum_into_values)
3266 for (
unsigned int i = 0; i < solution_values.
size(); ++i)
3267 solution_values[i] = 0;
3273 if (this->use_linear_path)
3274 do_integrate<false, true>(solution_values,
3278 do_integrate<false, false>(solution_values,
3285template <
int n_components_,
int dim,
int spacedim,
typename Number>
3290 const bool sum_into_values)
3293 solution_values.
size()),
3300template <
int n_components_,
int dim,
int spacedim,
typename Number>
3301template <
bool do_JxW,
bool is_linear, std::
size_t str
ide_view>
3306 const bool sum_into_values)
3308 if (!sum_into_values && this->fe->n_components() > n_components)
3309 for (
unsigned int i = 0; i < solution_values.
size(); ++i)
3310 solution_values[i] = 0;
3312 do_integrate_in_face<do_JxW, is_linear, 1>(this->scratch_data_scalar.begin(),
3316 ScalarNumber *input = this->scratch_data_scalar.begin();
3318 if (stride_view == 1 && this->component_in_base_element == 0 &&
3319 (is_linear || this->renumber.empty()))
3321 if (sum_into_values)
3323 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1,
ScalarNumber>::
3324 template interpolate<false, true>(n_components,
3328 solution_values.
data(),
3329 this->current_face_number);
3332 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1,
ScalarNumber>::
3333 template interpolate<false, false>(n_components,
3337 solution_values.
data(),
3338 this->current_face_number);
3342 const unsigned int dofs_per_comp_face =
3343 is_linear ?
Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3345 const unsigned int size_input = 3 * dofs_per_comp_face * n_components;
3349 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1,
ScalarNumber>::
3350 template interpolate<false, false>(n_components,
3355 this->current_face_number);
3357 const unsigned int dofs_per_comp =
3360 for (
unsigned int comp = 0; comp < n_components; ++comp)
3362 const std::size_t offset =
3363 (this->component_in_base_element + comp) * dofs_per_comp;
3365 if (is_linear || this->renumber.empty())
3367 for (
unsigned int i = 0; i < dofs_per_comp; ++i)
3368 if (sum_into_values)
3369 solution_values[i + offset] +=
3370 output[i + comp * dofs_per_comp];
3372 solution_values[i + offset] =
3373 output[i + comp * dofs_per_comp];
3377 const unsigned int *renumber_ptr = this->renumber.
data() + offset;
3378 for (
unsigned int i = 0; i < dofs_per_comp; ++i)
3379 if (sum_into_values)
3380 solution_values[renumber_ptr[i]] +=
3381 output[i + comp * dofs_per_comp];
3383 solution_values[renumber_ptr[i]] =
3384 output[i + comp * dofs_per_comp];
3392template <
int n_components_,
int dim,
int spacedim,
typename Number>
3393template <
int str
ide_face_dof>
3399 if (this->use_linear_path)
3400 do_evaluate_in_face<true, stride_face_dof>(face_dof_values,
3403 do_evaluate_in_face<false, stride_face_dof>(face_dof_values,
3409template <
int n_components_,
int dim,
int spacedim,
typename Number>
3410template <
bool is_linear,
int str
ide_face_dof>
3417 if constexpr (n_components == 1)
3418 face_dof_values_ptr = face_dof_values;
3421 const unsigned int dofs_per_comp_face =
3422 is_linear ?
Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3423 for (
unsigned int comp = 0; comp < n_components; ++comp)
3424 for (
unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3425 ETT::read_value(face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3428 this->solution_renumbered[i]);
3430 face_dof_values_ptr = this->solution_renumbered.data();
3433 constexpr int stride_face_dof_actual =
3434 n_components == 1 ? stride_face_dof : 1;
3437 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3439 for (
unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3446 const std::array<vectorized_value_type, dim + 1> interpolated_value =
3453 stride_face_dof_actual>(face_dof_values_ptr,
3454 this->unit_point_faces_ptr[qb]) :
3461 stride_face_dof_actual>(this->shapes_faces.data() +
3464 face_dof_values_ptr);
3466 value = interpolated_value[dim - 1];
3469 if (this->current_face_number / 2 == 0)
3471 gradient[0] = interpolated_value[dim];
3473 gradient[1] = interpolated_value[0];
3475 gradient[2] = interpolated_value[1];
3477 else if (this->current_face_number / 2 == 1)
3480 gradient[1] = interpolated_value[dim];
3483 gradient[0] = interpolated_value[1];
3484 gradient[2] = interpolated_value[0];
3487 gradient[0] = interpolated_value[0];
3491 else if (this->current_face_number / 2 == 2)
3495 gradient[0] = interpolated_value[0];
3496 gradient[1] = interpolated_value[1];
3497 gradient[2] = interpolated_value[dim];
3512 stride_face_dof_actual>(face_dof_values_ptr,
3513 this->unit_point_faces_ptr[qb]) :
3519 stride_face_dof_actual>(this->shapes_faces.data() +
3522 face_dof_values_ptr);
3527 for (
unsigned int v = 0, offset = qb * stride;
3528 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3530 ETT::set_value(value, v, this->values[offset]);
3538 for (
unsigned int v = 0, offset = qb * stride;
3539 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3543 ETT::set_gradient(gradient, v, unit_gradient);
3544 this->gradients[offset] =
3548 this->inverse_jacobian_ptr[0].
transpose(), unit_gradient) :
3551 ->inverse_jacobian_ptr[this->cell_type <=
3553 GeometryType::affine ?
3565template <
int n_components_,
int dim,
int spacedim,
typename Number>
3566template <
int str
ide_face_dof>
3571 const bool sum_into_values)
3573 if (this->use_linear_path)
3574 do_integrate_in_face<true, true, stride_face_dof>(face_dof_values,
3578 do_integrate_in_face<true, false, stride_face_dof>(face_dof_values,
3585template <
int n_components_,
int dim,
int spacedim,
typename Number>
3586template <
bool do_JxW,
bool is_linear,
int str
ide_face_dof>
3592 const bool sum_into_values)
3595 if constexpr (stride == 1)
3596 if (
const unsigned int n_filled_lanes =
3597 this->n_q_points_scalar & (n_lanes_internal - 1);
3601 for (
unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3602 ETT::set_zero_value(this->values.back(), v);
3604 for (
unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3605 ETT::set_zero_gradient(this->gradients.back(), v);
3610 solution_values_vectorized_linear = {};
3613 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3615 const bool cartesian_cell =
3617 const bool affine_cell =
3619 for (
unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3625 for (
unsigned int v = 0, offset = qb * stride;
3626 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3628 ETT::get_value(value,
3630 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
3631 this->values[offset]);
3634 for (
unsigned int v = 0, offset = qb * stride;
3635 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3639 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
3640 this->gradients[offset];
3648 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
3654 std::array<vectorized_value_type, 2> value_face = {};
3657 value_face[0] = value;
3660 if (this->current_face_number / 2 == 0)
3662 value_face[1] = gradient[0];
3664 gradient_in_face[0] = gradient[1];
3666 gradient_in_face[1] = gradient[2];
3668 else if (this->current_face_number / 2 == 1)
3671 value_face[1] = gradient[1];
3674 gradient_in_face[0] = gradient[2];
3675 gradient_in_face[1] = gradient[0];
3678 gradient_in_face[0] = gradient[0];
3682 else if (this->current_face_number / 2 == 2)
3686 value_face[1] = gradient[2];
3687 gradient_in_face[0] = gradient[0];
3688 gradient_in_face[1] = gradient[1];
3701 2>(this->shapes_faces.data() + qb * n_shapes,
3705 is_linear ? solution_values_vectorized_linear.data() :
3706 this->solution_renumbered_vectorized.data(),
3707 this->unit_point_faces_ptr[qb],
3715 this->shapes_faces.data() + qb * n_shapes,
3718 is_linear ? solution_values_vectorized_linear.data() :
3719 this->solution_renumbered_vectorized.data(),
3720 this->unit_point_faces_ptr[qb],
3724 const unsigned int dofs_per_comp_face =
3725 is_linear ?
Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3727 for (
unsigned int comp = 0; comp < n_components; ++comp)
3728 for (
unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3729 if (sum_into_values)
3730 face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3732 ETT::sum_value(comp,
3734 *(solution_values_vectorized_linear.data() + i) :
3735 this->solution_renumbered_vectorized[i]);
3737 face_dof_values[(i + comp * 3 * dofs_per_comp_face) * stride_face_dof] =
3738 ETT::sum_value(comp,
3740 *(solution_values_vectorized_linear.data() + i) :
3741 this->solution_renumbered_vectorized[i]);
3746template <
int n_components_,
int dim,
int spacedim,
typename Number>
3749 const unsigned int point_index)
const
3752 Assert(this->normal_ptr !=
nullptr,
3754 ExcFEPointEvaluationAccessToUninitializedMappingField(
3755 "update_normal_vectors"));
3759 for (
unsigned int d = 0; d < dim; ++d)
3763 if (this->is_interior)
3770 if (this->is_interior)
3771 return this->normal_ptr[point_index];
3773 return -(this->normal_ptr[point_index]);
value_type * data() const noexcept
typename ETT::vectorized_value_type vectorized_value_type
static constexpr std::size_t n_lanes_internal
void do_integrate_in_face(ScalarNumber *face_dof_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
void integrate_in_face(ScalarNumber *face_dof_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
typename ETT::scalar_value_type scalar_value_type
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
void evaluate(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
void do_evaluate_in_face(const ScalarNumber *face_dof_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
static constexpr std::size_t stride
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
void do_evaluate(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
static constexpr unsigned int n_components
static constexpr std::size_t n_lanes_user_interface
typename ETT::real_gradient_type gradient_type
void evaluate_in_face(const ScalarNumber *face_dof_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
void reinit(const unsigned int cell_index, const unsigned int face_number)
void test_and_sum(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
FEFacePointEvaluation(NonMatching::MappingInfo< dim, spacedim, Number > &mapping_info, const FiniteElement< dim, spacedim > &fe, const bool is_interior=true, const unsigned int first_selected_component=0)
typename ETT::unit_gradient_type unit_gradient_type
void do_integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
void integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
static constexpr unsigned int dimension
Tensor< 1, spacedim, Number > normal_vector(const unsigned int point_index) const
typename ETT::value_type value_type
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
const DerivativeForm< 1, dim, spacedim, Number > * jacobian_ptr
std::unique_ptr< NonMatching::MappingInfo< dim, spacedim, Number > > mapping_info_on_the_fly
std::vector< gradient_type > gradients
unsigned int dofs_per_component
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
unsigned int dofs_per_component_face
Number JxW(const unsigned int point_index) const
const UpdateFlags update_flags
static constexpr std::size_t n_lanes_user_interface
static constexpr std::size_t stride
internal::MatrixFreeFunctions::GeometryType cell_type
std::vector< Polynomials::Polynomial< double > > poly
Point< spacedim, Number > real_point(const unsigned int point_index) const
const unsigned int n_q_points_scalar
const Point< dim, VectorizedArrayType > * unit_point_ptr
FEPointEvaluationBase(FEPointEvaluationBase &&other) noexcept
AlignedVector< ScalarNumber > scratch_data_scalar
std_cxx20::ranges::iota_view< unsigned int, unsigned int > quadrature_point_indices() const
const value_type & get_value(const unsigned int point_index) const
const Point< dim - 1, VectorizedArrayType > * unit_point_faces_ptr
std::vector< scalar_value_type > solution_renumbered
SmartPointer< NonMatching::MappingInfo< dim, spacedim, Number > > mapping_info
const unsigned int n_q_batches
std::vector< std::array< bool, n_components > > nonzero_shape_function_component
Point< spacedim, Number > quadrature_point(const unsigned int point_index) const
const gradient_type & get_gradient(const unsigned int point_index) const
std::shared_ptr< FEValues< dim, spacedim > > fe_values
FEPointEvaluationBase(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
SmartPointer< const FiniteElement< dim, spacedim > > fe
std::vector< unsigned int > renumber
std::vector< value_type > values
AlignedVector<::ndarray< VectorizedArrayType, 2, dim - 1 > > shapes_faces
Point< dim, Number > unit_point(const unsigned int point_index) const
const unsigned int n_q_points
boost::signals2::connection connection_is_reinitialized
const Tensor< 1, spacedim, Number > * normal_ptr
DerivativeForm< 1, spacedim, dim, Number > inverse_jacobian(const unsigned int point_index) const
AlignedVector< vectorized_value_type > solution_renumbered_vectorized
static constexpr std::size_t n_lanes_internal
const DerivativeForm< 1, spacedim, dim, Number > * inverse_jacobian_ptr
AlignedVector<::ndarray< VectorizedArrayType, 2, dim > > shapes
typename ETT::value_type value_type
unsigned int current_cell_index
scalar_value_type integrate_value() const
void submit_gradient(const gradient_type &, const unsigned int point_index)
void setup(const unsigned int first_selected_component)
FEPointEvaluationBase(NonMatching::MappingInfo< dim, spacedim, Number > &mapping_info, const FiniteElement< dim, spacedim > &fe, const unsigned int first_selected_component=0, const bool is_interior=true)
typename ETT::scalar_value_type scalar_value_type
static constexpr unsigned int dimension
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
typename ETT::vectorized_value_type vectorized_value_type
static constexpr unsigned int n_components
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
unsigned int current_face_number
FEPointEvaluationBase(FEPointEvaluationBase &other) noexcept
DerivativeForm< 1, dim, spacedim, Number > jacobian(const unsigned int point_index) const
typename ETT::real_gradient_type gradient_type
const Point< spacedim, Number > * real_point_ptr
unsigned int component_in_base_element
SmartPointer< const Mapping< dim, spacedim > > mapping
internal::MatrixFreeFunctions::ShapeInfo< ScalarNumber > shape_info
void submit_value(const value_type &value, const unsigned int point_index)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static constexpr std::size_t stride
void integrate_slow(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
static constexpr unsigned int dimension
typename ETT::value_type value_type
void do_integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
void evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
typename ETT::scalar_value_type scalar_value_type
typename ETT::unit_gradient_type unit_gradient_type
void evaluate(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
void prepare_evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values)
void test_and_sum(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
void integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
typename ETT::vectorized_value_type vectorized_value_type
void compute_integrate_fast(const EvaluationFlags::EvaluationFlags &integration_flags, const unsigned int n_shapes, const unsigned int qb, const vectorized_value_type value, const interface_vectorized_unit_gradient_type gradient, vectorized_value_type *solution_values_vectorized_linear)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
typename ETT::real_gradient_type gradient_type
FEPointEvaluation(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
void evaluate_slow(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
static constexpr std::size_t n_lanes_internal
static constexpr unsigned int n_components
void finish_integrate_fast(const StridedArrayView< ScalarNumber, stride_view > &solution_values, vectorized_value_type *solution_values_vectorized_linear, const bool sum_into_values)
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
void compute_evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags, const unsigned int n_shapes, const unsigned int qb, vectorized_value_type &value, interface_vectorized_unit_gradient_type &gradient)
void integrate_fast(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
static constexpr std::size_t n_lanes_user_interface
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
Tensor< 1, spacedim, Number > normal_vector(const unsigned int point_index) const
Abstract base class for mapping classes.
value_type * data() const noexcept
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcFEPointEvaluationAccessToUninitializedMappingField(std::string arg1)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcNotInitialized()
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_ASSERT_UNREACHABLE()
EvaluationFlags
The EvaluationFlags enum.
constexpr T pow(const T base, const int iexp)
std::vector< Polynomials::Polynomial< double > > get_polynomial_space(const FiniteElement< dim, spacedim > &fe)
bool is_fast_path_supported(const FiniteElement< dim, spacedim > &fe, const unsigned int base_element_number)
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
void integrate_tensor_product_value_and_gradient(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values, const Point< dim, Number > &p, const bool do_add)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
void compute_values_of_array(::ndarray< Number, 2, dim > *shapes, const std::vector< Polynomials::Polynomial< double > > &poly, const Point< dim, Number > &p, const unsigned int derivative=1)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_linear(const Number *values, const Point< dim, Number2 > &p)
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_linear(const Number *values, const Point< dim, Number2 > &p)
void integrate_tensor_product_value(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 &value, Number2 *values, const Point< dim, Number > &p, const bool do_add)
static const unsigned int invalid_unsigned_int
boost::integer_range< IncrementableType > iota_view
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
static void set_zero_gradient(unit_gradient_type &value, const unsigned int vector_lane)
static void access(value_type &value, const unsigned int vector_lane, const unsigned int component, const ScalarNumber &shape_value)
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static Tensor< 1, dim, ScalarNumber > access(const real_gradient_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void set_zero_value(value_type &value, const unsigned int vector_lane)
static void set_gradient(const interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, unit_gradient_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int component, scalar_value_type &result)
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int component, const Tensor< 1, dim, ScalarNumber > &shape_gradient)
static scalar_value_type sum_value(const vectorized_value_type &result)
static scalar_value_type sum_value(const scalar_value_type &result)
static ScalarNumber sum_value(const unsigned int component, const vectorized_value_type &result)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const unit_gradient_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int)
static void get_gradient(vectorized_unit_gradient_type &value, const unsigned int, const vectorized_unit_gradient_type &result)
static scalar_value_type sum_value(const scalar_value_type &result)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void set_gradient(const vectorized_unit_gradient_type &value, const unsigned int vector_lane, scalar_unit_gradient_type &result)
static void set_zero_value(value_type &value, const unsigned int vector_lane)
static Tensor< 1, spacedim, ScalarNumber > access(const real_gradient_type &value, const unsigned int vector_lane, const unsigned int)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
VectorizedArrayType vectorized_value_type
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
static scalar_value_type sum_value(const vectorized_value_type &result)
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int, const Tensor< 1, spacedim, ScalarNumber > &shape_gradient)
static void set_gradient(const vectorized_unit_gradient_type &value, const unsigned int, vectorized_unit_gradient_type &result)
ScalarNumber scalar_value_type
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static void set_zero_gradient(real_gradient_type &value, const unsigned int vector_lane)
static void access(value_type &value, const unsigned int vector_lane, const unsigned int, const ScalarNumber &shape_value)
static void get_gradient(vectorized_unit_gradient_type &value, const unsigned int vector_lane, const scalar_unit_gradient_type &result)
static ScalarNumber sum_value(const unsigned int, const vectorized_value_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int, scalar_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static Tensor< 1, spacedim, ScalarNumber > access(const real_gradient_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const DerivativeForm< 1, dim, n_components, Number > &result)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
static scalar_value_type sum_value(const scalar_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static scalar_value_type sum_value(const vectorized_value_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int component, scalar_value_type &result)
Tensor< 1, n_components, ScalarNumber > scalar_value_type
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int component)
static ScalarNumber sum_value(const unsigned int component, const vectorized_value_type &result)
Tensor< 1, n_components, Tensor< 1, dim, VectorizedArrayType > > vectorized_unit_gradient_type
static void set_zero_value(value_type &value, const unsigned int vector_lane)
static void set_gradient(const interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, unit_gradient_type &result)
static void set_zero_gradient(real_gradient_type &value, const unsigned int vector_lane)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const unit_gradient_type &result)
std::conditional_t< n_components==spacedim, Tensor< 2, spacedim, Number >, Tensor< 1, n_components, Tensor< 1, spacedim, Number > > > real_gradient_type
Tensor< 1, n_components, VectorizedArrayType > vectorized_value_type
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int component, const Tensor< 1, spacedim, ScalarNumber > &shape_gradient)
Tensor< 1, n_components, Tensor< 1, dim, Number > > unit_gradient_type
static void access(value_type &value, const unsigned int vector_lane, const unsigned int component, const ScalarNumber &shape_value)
static constexpr std::size_t width()
static constexpr std::size_t stride()
static value_type & get(value_type &value, unsigned int c)