18 #ifndef __deal2__matrix_free_fe_evaluation_h
19 #define __deal2__matrix_free_fe_evaluation_h
22 #include <deal.II/base/config.h>
24 #include <deal.II/base/template_constraints.h>
25 #include <deal.II/base/symmetric_tensor.h>
26 #include <deal.II/base/vectorization.h>
27 #include <deal.II/matrix_free/matrix_free.h>
46 template <
typename FEEval>
47 void do_evaluate (FEEval &,
const bool,
const bool,
const bool);
48 template <
typename FEEval>
49 void do_integrate (FEEval &,
const bool,
const bool);
85 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
86 int n_components_,
typename Number>
90 typedef Number number_type;
93 static const unsigned int dimension = dim;
94 static const unsigned int n_components = n_components_;
95 static const unsigned int dofs_per_cell = dofs_per_cell_;
96 static const unsigned int n_q_points = n_q_points_;
127 internal::MatrixFreeFunctions::CellType
get_cell_type()
const;
151 template <
typename VectorType>
162 template <
typename VectorType>
164 const unsigned int first_index=0);
170 template <
typename VectorType>
172 const unsigned int first_index=0);
186 template <
typename VectorType>
196 template <
typename VectorType>
198 const unsigned int first_index=0);
205 template <
typename VectorType>
207 const unsigned int first_index=0);
216 template<
typename VectorType>
228 template<
typename VectorType>
230 const unsigned int first_index=0)
const;
236 template<
typename VectorType>
238 const unsigned int first_index=0)
const;
247 template<
typename VectorType>
259 template<
typename VectorType>
261 const unsigned int first_index=0)
const;
267 template<
typename VectorType>
269 const unsigned int first_index=0)
const;
303 const unsigned int dof);
316 value_type
get_value (
const unsigned int q_point)
const;
330 const unsigned int q_point);
341 gradient_type
get_gradient (
const unsigned int q_point)
const;
356 const unsigned int q_point);
519 const unsigned int fe_no = 0,
520 const unsigned int quad_no = 0);
528 template<
typename VectorType,
typename VectorOperation>
530 VectorType *vectors[])
const;
539 template<
typename VectorType>
613 const internal::MatrixFreeFunctions::DoFInfo &
dof_info;
747 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
748 int n_components_,
typename Number>
750 public FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
753 typedef Number number_type;
756 static const unsigned int dimension = dim;
757 static const unsigned int n_components = n_components_;
758 static const unsigned int dofs_per_cell = dofs_per_cell_;
759 static const unsigned int n_q_points = n_q_points_;
772 const unsigned int fe_no = 0,
773 const unsigned int quad_no = 0);
787 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
792 typedef Number number_type;
795 static const unsigned int dimension = dim;
796 static const unsigned int dofs_per_cell = dofs_per_cell_;
797 static const unsigned int n_q_points = n_q_points_;
816 const unsigned int dof);
825 value_type
get_value (
const unsigned int q_point)
const;
835 const unsigned int q_point);
842 gradient_type
get_gradient (
const unsigned int q_point)
const;
853 const unsigned int q_point);
895 const unsigned int fe_no = 0,
896 const unsigned int quad_no = 0);
910 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
915 typedef Number number_type;
918 static const unsigned int dimension = dim;
919 static const unsigned int n_components = dim;
920 static const unsigned int dofs_per_cell = dofs_per_cell_;
921 static const unsigned int n_q_points = n_q_points_;
928 gradient_type
get_gradient (
const unsigned int q_point)
const;
943 get_symmetric_gradient (
const unsigned int q_point)
const;
950 get_curl (
const unsigned int q_point)
const;
976 const unsigned int q_point);
987 const unsigned int q_point);
998 const unsigned int q_point);
1009 const unsigned int q_point);
1016 const unsigned int q_point);
1027 const unsigned int fe_no = 0,
1028 const unsigned int quad_no = 0);
1072 template <
int dim,
int fe_degree,
int n_q_points_1d = fe_degree+1,
1073 int n_components_ = 1,
typename Number =
double >
1076 Utilities::fixed_int_power<fe_degree+1,dim>::value,
1077 Utilities::fixed_int_power<n_q_points_1d,dim>::value,
1078 n_components_,Number>
1085 typedef Number number_type;
1088 static const unsigned int dimension = dim;
1089 static const unsigned int n_components = n_components_;
1090 static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell;
1091 static const unsigned int n_q_points = BaseClass::n_q_points;
1100 const unsigned int fe_no = 0,
1101 const unsigned int quad_no = 0);
1110 void evaluate (
const bool evaluate_val,
1111 const bool evaluate_grad,
1112 const bool evaluate_hess =
false);
1121 void integrate (
const bool integrate_val,
1122 const bool integrate_grad);
1140 template <
int direction,
bool dof_to_quad,
bool add>
1152 template <
int direction,
bool dof_to_quad,
bool add>
1165 template <
int direction,
bool dof_to_quad,
bool add>
1172 template <
typename FEEval>
friend void
1173 internal::do_evaluate (FEEval &,
const bool,
const bool,
const bool);
1174 template <
typename FEEval>
friend void
1175 internal::do_integrate (FEEval &,
const bool,
const bool);
1221 template <
int dim,
int fe_degree,
int n_q_points_1d = fe_degree+1,
1222 int n_components_ = 1,
typename Number =
double >
1228 typedef Number number_type;
1231 static const unsigned int dimension = dim;
1232 static const unsigned int n_components = n_components_;
1233 static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell;
1234 static const unsigned int n_q_points = BaseClass::n_q_points;
1243 const unsigned int fe_no = 0,
1244 const unsigned int quad_no = 0);
1254 void evaluate (
const bool evaluate_val,
1255 const bool evaluate_grad,
1256 const bool evaluate_hess =
false);
1265 void integrate (
const bool integrate_val,
1266 const bool integrate_grad);
1278 template <
int direction,
bool dof_to_quad,
bool add>
1290 template <
int direction,
bool dof_to_quad,
bool add>
1303 template <
int direction,
bool dof_to_quad,
bool add>
1314 template <
typename FEEval>
friend void
1315 internal::do_evaluate (FEEval &,
const bool,
const bool,
const bool);
1316 template <
typename FEEval>
friend void
1317 internal::do_integrate (FEEval &,
const bool,
const bool);
1359 template <
int dim,
int fe_degree,
int n_components_ = 1,
typename Number =
double >
1361 public FEEvaluation<dim,fe_degree,fe_degree+1,n_components_,Number>
1365 typedef Number number_type;
1368 static const unsigned int dimension = dim;
1369 static const unsigned int n_components = n_components_;
1370 static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell;
1371 static const unsigned int n_q_points = BaseClass::n_q_points;
1380 const unsigned int fe_no = 0,
1381 const unsigned int quad_no = 0);
1391 void evaluate (
const bool evaluate_val,
1392 const bool evaluate_grad,
1393 const bool evaluate_lapl =
false);
1402 void integrate (
const bool integrate_val,
1403 const bool integrate_grad);
1414 template <
int direction,
bool dof_to_quad,
bool add>
1429 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
1430 int n_components_,
typename Number>
1434 const unsigned int fe_no_in,
1435 const unsigned int quad_no_in)
1437 quad_no (quad_no_in),
1438 n_fe_components (data_in.get_dof_info(fe_no_in).
n_components),
1439 active_fe_index (data_in.get_dof_info(fe_no_in).fe_index_from_dofs_per_cell
1440 (dofs_per_cell_ * n_fe_components)),
1441 active_quad_index (data_in.get_mapping_info().
1442 mapping_data_gen[quad_no_in].
1443 quad_index_from_n_q_points(n_q_points_)),
1444 matrix_info (data_in),
1445 dof_info (data_in.get_dof_info(fe_no_in)),
1446 mapping_info (data_in.get_mapping_info()),
1447 data (data_in.get_shape_info
1448 (fe_no_in, quad_no_in, active_fe_index,
1449 active_quad_index)),
1453 quadrature_weights (mapping_info.mapping_data_gen[quad_no].
1454 quadrature_weights[active_quad_index].begin()),
1455 quadrature_points (0),
1457 jacobian_grad_upper(0),
1458 cell (
numbers::invalid_unsigned_int),
1459 cell_type (
internal::MatrixFreeFunctions::undefined),
1460 cell_data_number (0)
1466 Assert (n_fe_components == 1 ||
1469 ExcMessage (
"The underlying FE is vector-valued. In this case, the "
1470 "template argument n_components must be a the same "
1471 "as the number of underlying vector components."));
1480 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
1481 int n_components_,
typename Number>
1485 ::reinit (
const unsigned int cell_in)
1489 dof_info.cell_active_fe_index[cell_in] : 0),
1498 mapping_data_gen[quad_no].rowstart_q_points.size());
1500 rowstart_q_points[cell];
1502 quadrature_points.size());
1507 if (cell_type == internal::MatrixFreeFunctions::cartesian)
1509 cartesian_data = &mapping_info.
cartesian_data[cell_data_number].first;
1512 else if (cell_type == internal::MatrixFreeFunctions::affine)
1514 jacobian = &mapping_info.
affine_data[cell_data_number].first;
1515 J_value = &mapping_info.
affine_data[cell_data_number].second;
1519 const unsigned int rowstart = mapping_info.
1520 mapping_data_gen[quad_no].rowstart_jacobians[cell_data_number];
1522 mapping_data_gen[quad_no].jacobians.size());
1528 mapping_data_gen[quad_no].JxW_values.size());
1530 JxW_values[rowstart]);
1535 mapping_data_gen[quad_no].jacobians_grad_diag.size());
1537 jacobians_grad_diag[rowstart];
1539 mapping_data_gen[quad_no].jacobians_grad_upper.size());
1541 jacobians_grad_upper[rowstart];
1546 dof_values_initialized =
false;
1547 values_quad_initialized =
false;
1548 gradients_quad_initialized =
false;
1549 hessians_quad_initialized =
false;
1555 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
1556 int n_components_,
typename Number>
1563 return cell_data_number;
1568 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
1569 int n_components_,
typename Number>
1571 internal::MatrixFreeFunctions::CellType
1584 template <
typename VectorType>
1586 typename VectorType::value_type &
1587 vector_access (VectorType &vec,
1588 const unsigned int entry)
1596 template <
typename VectorType>
1598 typename VectorType::value_type
1599 vector_access (
const VectorType &vec,
1600 const unsigned int entry)
1610 template <
typename Number>
1614 const unsigned int entry)
1624 template <
typename Number>
1628 const unsigned int entry)
1637 template <
typename VectorType>
1639 void check_vector_compatibility (
const VectorType &vec,
1640 const internal::MatrixFreeFunctions::DoFInfo &dof_info)
1643 dof_info.vector_partitioner->size());
1646 template <
typename Number>
1649 const internal::MatrixFreeFunctions::DoFInfo &dof_info)
1652 ExcMessage(
"The parallel layout of the given vector is not "
1653 "compatible with the parallel partitioning in MatrixFree. "
1654 "Use MatrixFree::initialize_dof_vector to get a "
1655 "compatible vector."));
1659 template <
typename Number>
1662 template <
typename VectorType>
1663 void process_dof (
const unsigned int index,
1667 res = vector_access (const_cast<const VectorType &>(vec), index);
1670 void pre_constraints (
const Number &,
1676 template <
typename VectorType>
1677 void process_constraint (
const unsigned int index,
1678 const Number weight,
1682 res += weight * vector_access (const_cast<const VectorType &>(vec), index);
1685 void post_constraints (
const Number &sum,
1686 Number &write_pos)
const
1691 void process_empty (Number &res)
const
1698 template <
typename Number>
1699 struct VectorDistributorLocalToGlobal
1701 template <
typename VectorType>
1702 void process_dof (
const unsigned int index,
1706 vector_access (vec, index) += res;
1709 void pre_constraints (
const Number &input,
1715 template <
typename VectorType>
1716 void process_constraint (
const unsigned int index,
1717 const Number weight,
1721 vector_access (vec, index) += weight * res;
1724 void post_constraints (
const Number &,
1729 void process_empty (Number &)
const
1736 template <
typename Number>
1739 template <
typename VectorType>
1740 void process_dof (
const unsigned int index,
1744 vector_access (vec, index) = res;
1747 void pre_constraints (
const Number &,
1752 template <
typename VectorType>
1753 void process_constraint (
const unsigned int,
1760 void post_constraints (
const Number &,
1765 void process_empty (Number &)
const
1773 template <
typename VectorType,
bool>
1774 struct BlockVectorSelector {};
1776 template <
typename VectorType>
1777 struct BlockVectorSelector<VectorType,true>
1779 typedef typename VectorType::BlockType BaseVectorType;
1781 static BaseVectorType *get_vector_component (VectorType &vec,
1782 const unsigned int component)
1785 return &vec.block(component);
1789 template <
typename VectorType>
1790 struct BlockVectorSelector<VectorType,false>
1792 typedef VectorType BaseVectorType;
1794 static BaseVectorType *get_vector_component (VectorType &vec,
1805 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
1806 int n_components_,
typename Number>
1807 template<
typename VectorType,
typename VectorOperation>
1812 VectorType *src[])
const
1828 const unsigned int *dof_indices = dof_info.begin_indices(cell);
1829 const std::pair<unsigned short,unsigned short> *indicators =
1830 dof_info.begin_indicators(cell);
1831 const std::pair<unsigned short,unsigned short> *indicators_end =
1832 dof_info.end_indicators(cell);
1833 unsigned int ind_local = 0;
1835 const unsigned int n_irreg_components_filled = dof_info.row_starts[cell][2];
1836 const bool at_irregular_cell = n_irreg_components_filled > 0;
1840 if (n_fe_components == 1)
1842 const unsigned int n_local_dofs =
1845 internal::check_vector_compatibility (*src[comp], dof_info);
1849 const_cast<Number *>(&values_dofs[comp][0][0]);
1853 if (at_irregular_cell ==
false)
1856 if (indicators != indicators_end)
1858 for ( ; indicators != indicators_end; ++indicators)
1861 for (
unsigned int j=0; j<indicators->first; ++j)
1863 operation.process_dof (dof_indices[j], *src[comp],
1864 local_data[comp][ind_local+j]);
1866 ind_local += indicators->first;
1867 dof_indices += indicators->first;
1873 operation.pre_constraints (local_data[comp][ind_local],
1876 const Number *data_val =
1878 const Number *end_pool =
1880 for ( ; data_val != end_pool; ++data_val, ++dof_indices)
1882 operation.process_constraint (*dof_indices, *data_val,
1883 *src[comp], value[comp]);
1886 operation.post_constraints (value[comp],
1887 local_data[comp][ind_local]);
1893 for (; ind_local < n_local_dofs; ++dof_indices, ++ind_local)
1896 operation.process_dof (*dof_indices, *src[comp],
1897 local_data[comp][ind_local]);
1905 static_cast<int>(n_local_dofs));
1906 for (
unsigned int j=0; j<n_local_dofs; ++j)
1908 operation.process_dof (dof_indices[j], *src[comp],
1909 local_data[comp][j]);
1919 for ( ; indicators != indicators_end; ++indicators)
1921 for (
unsigned int j=0; j<indicators->first; ++j)
1926 operation.process_dof (dof_indices[j], *src[comp],
1927 local_data[comp][ind_local]);
1932 >= n_irreg_components_filled)
1935 operation.process_empty (local_data[comp][ind_local]);
1939 dof_indices += indicators->first;
1945 operation.pre_constraints (local_data[comp][ind_local],
1948 const Number *data_val =
1950 const Number *end_pool =
1953 for ( ; data_val != end_pool; ++data_val, ++dof_indices)
1955 operation.process_constraint (*dof_indices, *data_val,
1956 *src[comp], value[comp]);
1959 operation.post_constraints (value[comp],
1960 local_data[comp][ind_local]);
1963 >= n_irreg_components_filled)
1966 operation.process_empty (local_data[comp][ind_local]);
1970 for (; ind_local<n_local_dofs; ++dof_indices)
1972 Assert (dof_indices != dof_info.end_indices(cell),
1978 operation.process_dof (*dof_indices, *src[comp],
1979 local_data[comp][ind_local]);
1982 >= n_irreg_components_filled)
1985 operation.process_empty(local_data[comp][ind_local]);
1997 internal::check_vector_compatibility (*src[0], dof_info);
1999 const unsigned int n_local_dofs =
2001 Number *local_data =
2002 const_cast<Number *
>(&values_dofs[0][0][0]);
2003 if (at_irregular_cell ==
false)
2006 if (indicators != indicators_end)
2008 for ( ; indicators != indicators_end; ++indicators)
2011 for (
unsigned int j=0; j<indicators->first; ++j)
2012 operation.process_dof (dof_indices[j], *src[0],
2013 local_data[ind_local+j]);
2014 ind_local += indicators->first;
2015 dof_indices += indicators->first;
2020 operation.pre_constraints (local_data[ind_local], value);
2022 const Number *data_val =
2024 const Number *end_pool =
2027 for ( ; data_val != end_pool; ++data_val, ++dof_indices)
2028 operation.process_constraint (*dof_indices, *data_val,
2031 operation.post_constraints (value, local_data[ind_local]);
2036 for (; ind_local<n_local_dofs; ++dof_indices, ++ind_local)
2037 operation.process_dof (*dof_indices, *src[0],
2038 local_data[ind_local]);
2039 Assert (dof_indices == dof_info.end_indices(cell),
2047 static_cast<int>(n_local_dofs));
2048 for (
unsigned int j=0; j<n_local_dofs; ++j)
2049 operation.process_dof (dof_indices[j], *src[0],
2060 for ( ; indicators != indicators_end; ++indicators)
2062 for (
unsigned int j=0; j<indicators->first; ++j)
2066 operation.process_dof (dof_indices[j], *src[0],
2067 local_data[ind_local]);
2072 >= n_irreg_components_filled)
2074 operation.process_empty (local_data[ind_local]);
2078 dof_indices += indicators->first;
2083 operation.pre_constraints (local_data[ind_local], value);
2085 const Number *data_val =
2087 const Number *end_pool =
2090 for ( ; data_val != end_pool; ++data_val, ++dof_indices)
2091 operation.process_constraint (*dof_indices, *data_val,
2094 operation.post_constraints (value, local_data[ind_local]);
2097 >= n_irreg_components_filled)
2099 operation.process_empty (local_data[ind_local]);
2103 for (; ind_local<n_local_dofs; ++dof_indices)
2105 Assert (dof_indices != dof_info.end_indices(cell),
2110 operation.process_dof (*dof_indices, *src[0],
2111 local_data[ind_local]);
2114 >= n_irreg_components_filled)
2116 operation.process_empty (local_data[ind_local]);
2126 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2127 int n_components_,
typename Number>
2128 template<
typename VectorType>
2136 typename internal::BlockVectorSelector<VectorType,
2141 internal::VectorReader<Number> reader;
2142 read_write_operation (reader, src_data);
2145 dof_values_initialized =
true;
2151 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2152 int n_components_,
typename Number>
2153 template<
typename VectorType>
2158 const unsigned int first_index)
2162 Assert ((n_fe_components == 1 ?
2163 (first_index+n_components <= src.size()) :
true),
2164 ExcIndexRange (first_index + n_components_, 0, src.size()));
2168 src_data[comp] = const_cast<VectorType *>(&src[comp+first_index]);
2170 internal::VectorReader<Number> reader;
2171 read_write_operation (reader, src_data);
2174 dof_values_initialized =
true;
2180 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2181 int n_components_,
typename Number>
2182 template<
typename VectorType>
2187 const unsigned int first_index)
2191 Assert ((n_fe_components == 1 ?
2192 (first_index+n_components <= src.size()) :
true),
2193 ExcIndexRange (first_index + n_components_, 0, src.size()));
2197 src_data[comp] = const_cast<VectorType *>(src[comp+first_index]);
2199 internal::VectorReader<Number> reader;
2200 read_write_operation (reader, src_data);
2203 dof_values_initialized =
true;
2209 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2210 int n_components_,
typename Number>
2211 template<
typename VectorType>
2219 const typename internal::BlockVectorSelector<VectorType,
2224 read_dof_values_plain (src_data);
2229 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2230 int n_components_,
typename Number>
2231 template<
typename VectorType>
2236 const unsigned int first_index)
2240 Assert ((n_fe_components == 1 ?
2241 (first_index+n_components <= src.size()) :
true),
2242 ExcIndexRange (first_index + n_components_, 0, src.size()));
2245 src_data[comp] = &src[comp+first_index];
2246 read_dof_values_plain (src_data);
2251 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2252 int n_components_,
typename Number>
2253 template<
typename VectorType>
2258 const unsigned int first_index)
2262 Assert ((n_fe_components == 1 ?
2263 (first_index+n_components <= src.size()) :
true),
2264 ExcIndexRange (first_index + n_components_, 0, src.size()));
2267 src_data[comp] = src[comp+first_index];
2268 read_dof_values_plain (src_data);
2273 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2274 int n_components_,
typename Number>
2275 template<
typename VectorType>
2281 Assert (dof_values_initialized==
true,
2282 internal::ExcAccessToUninitializedField());
2286 typename internal::BlockVectorSelector<VectorType,
2291 internal::VectorDistributorLocalToGlobal<Number> distributor;
2292 read_write_operation (distributor, dst_data);
2297 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2298 int n_components_,
typename Number>
2299 template<
typename VectorType>
2304 const unsigned int first_index)
const
2308 Assert ((n_fe_components == 1 ?
2309 (first_index+n_components <= dst.size()) :
true),
2310 ExcIndexRange (first_index + n_components_, 0, dst.size()));
2311 Assert (dof_values_initialized==
true,
2312 internal::ExcAccessToUninitializedField());
2316 dst_data[comp] = &dst[comp+first_index];
2318 internal::VectorDistributorLocalToGlobal<Number> distributor;
2319 read_write_operation (distributor, dst_data);
2324 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2325 int n_components_,
typename Number>
2326 template<
typename VectorType>
2331 const unsigned int first_index)
const
2335 Assert ((n_fe_components == 1 ?
2336 (first_index+n_components <= dst.size()) :
true),
2337 ExcIndexRange (first_index + n_components_, 0, dst.size()));
2338 Assert (dof_values_initialized==
true,
2339 internal::ExcAccessToUninitializedField());
2343 dst_data[comp] = dst[comp+first_index];
2345 internal::VectorDistributorLocalToGlobal<Number> distributor;
2346 read_write_operation (distributor, dst_data);
2351 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2352 int n_components_,
typename Number>
2353 template<
typename VectorType>
2359 Assert (dof_values_initialized==
true,
2360 internal::ExcAccessToUninitializedField());
2364 typename internal::BlockVectorSelector<VectorType,
2369 internal::VectorSetter<Number> setter;
2370 read_write_operation (setter, dst_data);
2375 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2376 int n_components_,
typename Number>
2377 template<
typename VectorType>
2382 const unsigned int first_index)
const
2386 Assert ((n_fe_components == 1 ?
2387 (first_index+n_components <= dst.size()) :
true),
2388 ExcIndexRange (first_index + n_components_, 0, dst.size()));
2390 Assert (dof_values_initialized==
true,
2391 internal::ExcAccessToUninitializedField());
2395 dst_data[comp] = &dst[comp+first_index];
2397 internal::VectorSetter<Number> setter;
2398 read_write_operation (setter, dst_data);
2403 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2404 int n_components_,
typename Number>
2405 template<
typename VectorType>
2410 const unsigned int first_index)
const
2414 Assert ((n_fe_components == 1 ?
2415 (first_index+n_components <= dst.size()) :
true),
2416 ExcIndexRange (first_index + n_components_, 0, dst.size()));
2418 Assert (dof_values_initialized==
true,
2419 internal::ExcAccessToUninitializedField());
2423 dst_data[comp] = dst[comp+first_index];
2425 internal::VectorSetter<Number> setter;
2426 read_write_operation (setter, dst_data);
2431 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2432 int n_components_,
typename Number>
2433 template<
typename VectorType>
2449 const unsigned int *dof_indices = dof_info.begin_indices_plain(cell);
2451 const unsigned int n_irreg_components_filled = dof_info.row_starts[cell][2];
2452 const bool at_irregular_cell = n_irreg_components_filled > 0;
2456 if (n_fe_components == 1)
2458 const unsigned int n_local_dofs =
2461 internal::check_vector_compatibility (*src[comp], dof_info);
2464 local_src_number[comp] = &values_dofs[comp][0][0];
2468 if (at_irregular_cell ==
false)
2470 for (
unsigned int j=0; j<n_local_dofs; ++j)
2472 local_src_number[comp][j] =
2473 internal::vector_access (*src[comp], dof_indices[j]);
2482 for (
unsigned int ind_local=0; ind_local<n_local_dofs;
2488 local_src_number[comp][ind_local] =
2489 internal::vector_access (*src[comp], *dof_indices);
2494 local_src_number[comp][ind_local] = 0.;
2506 internal::check_vector_compatibility (*src[0], dof_info);
2508 const unsigned int n_local_dofs =
2510 Number *local_src_number = &values_dofs[0][0][0];
2511 if (at_irregular_cell ==
false)
2513 for (
unsigned int j=0; j<n_local_dofs; ++j)
2514 local_src_number[j] =
2515 internal::vector_access (*src[0], dof_indices[j]);
2524 for (
unsigned int ind_local=0; ind_local<n_local_dofs; ++dof_indices)
2528 local_src_number[ind_local] =
2529 internal::vector_access (*src[0], *dof_indices);
2533 local_src_number[ind_local] = 0.;
2541 dof_values_initialized =
true;
2551 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2558 return &values_dofs[0][0];
2563 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2571 dof_values_initialized =
true;
2573 return &values_dofs[0][0];
2578 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2585 Assert (values_quad_initialized || values_quad_submitted,
2587 return &values_quad[0][0];
2592 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2600 values_quad_submitted =
true;
2602 return &values_quad[0][0];
2607 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2614 Assert (gradients_quad_initialized || gradients_quad_submitted,
2616 return &gradients_quad[0][0][0];
2621 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2629 gradients_quad_submitted =
true;
2631 return &gradients_quad[0][0][0];
2636 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2644 return &hessians_quad[0][0][0];
2649 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2656 return &hessians_quad[0][0][0];
2661 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2662 int n_components_,
typename Number>
2671 return_value[comp] = this->values_dofs[comp][dof];
2672 return return_value;
2677 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2678 int n_components_,
typename Number>
2684 Assert (this->values_quad_initialized==
true,
2685 internal::ExcAccessToUninitializedField());
2689 return_value[comp] = this->values_quad[comp][q_point];
2690 return return_value;
2695 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2696 int n_components_,
typename Number>
2702 Assert (this->gradients_quad_initialized==
true,
2703 internal::ExcAccessToUninitializedField());
2709 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
2712 for (
unsigned int d=0; d<dim; ++d)
2713 grad_out[comp][d] = (this->gradients_quad[comp][d][q_point] *
2714 cartesian_data[0][d]);
2720 this->cell_type == internal::MatrixFreeFunctions::general ?
2721 jacobian[q_point] : jacobian[0];
2724 for (
unsigned int d=0; d<dim; ++d)
2726 grad_out[comp][d] = (jac[d][0] *
2727 this->gradients_quad[comp][0][q_point]);
2728 for (
unsigned int e=1; e<dim; ++e)
2729 grad_out[comp][d] += (jac[d][e] *
2730 this->gradients_quad[comp][e][q_point]);
2743 template <
int dim,
int n_q_po
ints,
typename Number>
2748 const unsigned int q_point,
2751 for (
unsigned int d=0; d<dim; ++d)
2756 tmp[0][0] = jac[0][0] * hessians_quad[0][q_point];
2759 tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
2760 jac[d][1] * hessians_quad[2][q_point]);
2761 tmp[1][d] = (jac[d][0] * hessians_quad[2][q_point] +
2762 jac[d][1] * hessians_quad[1][q_point]);
2765 tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
2766 jac[d][1] * hessians_quad[3][q_point] +
2767 jac[d][2] * hessians_quad[4][q_point]);
2768 tmp[1][d] = (jac[d][0] * hessians_quad[3][q_point] +
2769 jac[d][1] * hessians_quad[1][q_point] +
2770 jac[d][2] * hessians_quad[5][q_point]);
2771 tmp[2][d] = (jac[d][0] * hessians_quad[4][q_point] +
2772 jac[d][1] * hessians_quad[5][q_point] +
2773 jac[d][2] * hessians_quad[2][q_point]);
2784 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2785 int n_components_,
typename Number>
2791 Assert (this->hessians_quad_initialized==
true,
2792 internal::ExcAccessToUninitializedField());
2798 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
2802 for (
unsigned int d=0; d<dim; ++d)
2804 hessian_out[comp][d][d] = (this->hessians_quad[comp][d][q_point] *
2811 hessian_out[comp][0][1] = (this->hessians_quad[comp][2][q_point] *
2815 hessian_out[comp][0][1] = (this->hessians_quad[comp][3][q_point] *
2817 hessian_out[comp][0][2] = (this->hessians_quad[comp][4][q_point] *
2819 hessian_out[comp][1][2] = (this->hessians_quad[comp][5][q_point] *
2825 for (
unsigned int e=d+1; e<dim; ++e)
2826 hessian_out[comp][e][d] = hessian_out[comp][d][e];
2830 else if (this->cell_type == internal::MatrixFreeFunctions::general)
2838 & jac_grad_UT = jacobian_grad_upper[q_point];
2844 internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
2848 for (
unsigned int d=0; d<dim; ++d)
2849 for (
unsigned int e=d; e<dim; ++e)
2851 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
2852 for (
unsigned int f=1; f<dim; ++f)
2853 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
2857 for (
unsigned int d=0; d<dim; ++d)
2858 for (
unsigned int e=0; e<dim; ++e)
2859 hessian_out[comp][d][d] += (jac_grad[d][e] *
2860 this->gradients_quad[comp][e][q_point]);
2863 for (
unsigned int d=0, count=0; d<dim; ++d)
2864 for (
unsigned int e=d+1; e<dim; ++e, ++count)
2865 for (
unsigned int f=0; f<dim; ++f)
2866 hessian_out[comp][d][e] += (jac_grad_UT[count][f] *
2867 this->gradients_quad[comp][f][q_point]);
2870 for (
unsigned int d=0; d<dim; ++d)
2871 for (
unsigned int e=d+1; e<dim; ++e)
2872 hessian_out[comp][e][d] = hessian_out[comp][d][e];
2884 internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
2888 for (
unsigned int d=0; d<dim; ++d)
2889 for (
unsigned int e=d; e<dim; ++e)
2891 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
2892 for (
unsigned int f=1; f<dim; ++f)
2893 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
2900 for (
unsigned int d=0; d<dim; ++d)
2901 for (
unsigned int e=d+1; e<dim; ++e)
2902 hessian_out[comp][e][d] = hessian_out[comp][d][e];
2910 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2911 int n_components_,
typename Number>
2917 Assert (this->hessians_quad_initialized==
true,
2918 internal::ExcAccessToUninitializedField());
2924 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
2928 for (
unsigned int d=0; d<dim; ++d)
2929 hessian_out[comp][d] = (this->hessians_quad[comp][d][q_point] *
2933 else if (this->cell_type == internal::MatrixFreeFunctions::general)
2944 internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
2949 for (
unsigned int d=0; d<dim; ++d)
2951 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
2952 for (
unsigned int f=1; f<dim; ++f)
2953 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
2956 for (
unsigned int d=0; d<dim; ++d)
2957 for (
unsigned int e=0; e<dim; ++e)
2958 hessian_out[comp][d] += (jac_grad[d][e] *
2959 this->gradients_quad[comp][e][q_point]);
2971 internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
2976 for (
unsigned int d=0; d<dim; ++d)
2978 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
2979 for (
unsigned int f=1; f<dim; ++f)
2980 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
2989 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
2990 int n_components_,
typename Number>
2996 Assert (this->hessians_quad_initialized==
true,
2997 internal::ExcAccessToUninitializedField());
3001 = get_hessian_diagonal(q_point);
3004 laplacian_out[comp] = hess_diag[comp][0];
3005 for (
unsigned int d=1; d<dim; ++d)
3006 laplacian_out[comp] += hess_diag[comp][d];
3008 return laplacian_out;
3013 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
3014 int n_components_,
typename Number>
3019 const unsigned int dof)
3022 this->dof_values_initialized =
true;
3026 this->values_dofs[comp][dof] = val_in[comp];
3031 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
3032 int n_components_,
typename Number>
3037 const unsigned int q_point)
3042 this->values_quad_submitted =
true;
3044 if (this->cell_type == internal::MatrixFreeFunctions::general)
3048 this->values_quad[comp][q_point] = val_in[comp] * JxW;
3054 this->values_quad[comp][q_point] = val_in[comp] * JxW;
3060 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
3061 int n_components_,
typename Number>
3067 const unsigned int q_point)
3072 this->gradients_quad_submitted =
true;
3074 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
3078 for (
unsigned int d=0; d<dim; ++d)
3079 this->gradients_quad[comp][d][q_point] = (grad_in[comp][d] *
3080 cartesian_data[0][d] * JxW);
3085 this->cell_type == internal::MatrixFreeFunctions::general ?
3086 jacobian[q_point] : jacobian[0];
3088 this->cell_type == internal::MatrixFreeFunctions::general ?
3089 J_value[q_point] : J_value[0] * quadrature_weights[q_point];
3091 for (
unsigned int d=0; d<dim; ++d)
3094 for (
unsigned int e=1; e<dim; ++e)
3095 new_val += (jac[e][d] * grad_in[comp][e]);
3096 this->gradients_quad[comp][d][q_point] = new_val * JxW;
3103 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
3104 int n_components_,
typename Number>
3112 Assert (this->values_quad_submitted ==
true,
3113 internal::ExcAccessToUninitializedField());
3117 return_value[comp] = this->values_quad[comp][0];
3118 for (
unsigned int q=1; q<n_q_points; ++q)
3120 return_value[comp] += this->values_quad[comp][q];
3121 return (return_value);
3129 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
3130 int n_components_,
typename Number>
3134 const unsigned int fe_no,
3135 const unsigned int quad_no_in)
3138 (data_in, fe_no, quad_no_in)
3147 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3151 const unsigned int fe_no,
3152 const unsigned int quad_no_in)
3155 (data_in, fe_no, quad_no_in)
3160 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3167 return this->values_dofs[0][dof];
3172 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3178 Assert (this->values_quad_initialized==
true,
3179 internal::ExcAccessToUninitializedField());
3181 return this->values_quad[0][q_point];
3186 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3195 Assert (this->gradients_quad_initialized==
true,
3196 internal::ExcAccessToUninitializedField());
3202 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
3204 for (
unsigned int d=0; d<dim; ++d)
3205 grad_out[d] = (this->gradients_quad[0][d][q_point] *
3206 this->cartesian_data[0][d]);
3212 this->cell_type == internal::MatrixFreeFunctions::general ?
3213 this->jacobian[q_point] : this->jacobian[0];
3214 for (
unsigned int d=0; d<dim; ++d)
3216 grad_out[d] = (jac[d][0] * this->gradients_quad[0][0][q_point]);
3217 for (
unsigned int e=1; e<dim; ++e)
3218 grad_out[d] += (jac[d][e] * this->gradients_quad[0][e][q_point]);
3226 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3232 return BaseClass::get_hessian(q_point)[0];
3237 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3243 return BaseClass::get_hessian_diagonal(q_point)[0];
3248 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3254 return BaseClass::get_laplacian(q_point)[0];
3259 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3264 const unsigned int dof)
3267 this->dof_values_initialized =
true;
3270 this->values_dofs[0][dof] = val_in;
3275 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3280 const unsigned int q_point)
3285 this->values_quad_submitted =
true;
3287 if (this->cell_type == internal::MatrixFreeFunctions::general)
3290 this->values_quad[0][q_point] = val_in * JxW;
3295 this->values_quad[0][q_point] = val_in * JxW;
3301 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3306 const unsigned int q_point)
3311 this->gradients_quad_submitted =
true;
3313 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
3316 for (
unsigned int d=0; d<dim; ++d)
3317 this->gradients_quad[0][d][q_point] = (grad_in[d] *
3318 this->cartesian_data[0][d] *
3325 this->cell_type == internal::MatrixFreeFunctions::general ?
3326 this->jacobian[q_point] : this->jacobian[0];
3328 this->cell_type == internal::MatrixFreeFunctions::general ?
3329 this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
3330 for (
unsigned int d=0; d<dim; ++d)
3333 for (
unsigned int e=1; e<dim; ++e)
3334 new_val += jac[e][d] * grad_in[e];
3335 this->gradients_quad[0][d][q_point] = new_val * JxW;
3342 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3348 return BaseClass::integrate_value()[0];
3357 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3361 const unsigned int fe_no,
3362 const unsigned int quad_no_in)
3365 (data_in, fe_no, quad_no_in)
3370 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3376 return BaseClass::get_gradient (q_point);
3381 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3387 Assert (this->gradients_quad_initialized==
true,
3388 internal::ExcAccessToUninitializedField());
3394 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
3396 divergence = (this->gradients_quad[0][0][q_point] *
3397 this->cartesian_data[0][0]);
3398 for (
unsigned int d=1; d<dim; ++d)
3399 divergence += (this->gradients_quad[d][d][q_point] *
3400 this->cartesian_data[0][d]);
3406 this->cell_type == internal::MatrixFreeFunctions::general ?
3407 this->jacobian[q_point] : this->jacobian[0];
3408 divergence = (jac[0][0] * this->gradients_quad[0][0][q_point]);
3409 for (
unsigned int e=1; e<dim; ++e)
3410 divergence += (jac[0][e] * this->gradients_quad[0][e][q_point]);
3411 for (
unsigned int d=1; d<dim; ++d)
3412 for (
unsigned int e=0; e<dim; ++e)
3413 divergence += (jac[d][e] * this->gradients_quad[d][e][q_point]);
3420 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3430 for (
unsigned int d=0; d<dim; ++d)
3431 symmetrized[d] = grad[d][d];
3437 symmetrized[2] = grad[0][1] + grad[1][0];
3438 symmetrized[2] *= half;
3441 symmetrized[3] = grad[0][1] + grad[1][0];
3442 symmetrized[3] *= half;
3443 symmetrized[4] = grad[0][2] + grad[2][0];
3444 symmetrized[4] *= half;
3445 symmetrized[5] = grad[1][2] + grad[2][1];
3446 symmetrized[5] *= half;
3456 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3460 ::get_curl (
const unsigned int q_point)
const
3469 ExcMessage(
"Computing the curl in 1d is not a useful operation"));
3472 curl[0] = grad[1][0] - grad[0][1];
3475 curl[0] = grad[2][1] - grad[1][2];
3476 curl[1] = grad[0][2] - grad[2][0];
3477 curl[2] = grad[1][0] - grad[0][1];
3487 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3493 Assert (this->hessians_quad_initialized==
true,
3494 internal::ExcAccessToUninitializedField());
3497 return BaseClass::get_hessian_diagonal (q_point);
3502 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3508 Assert (this->hessians_quad_initialized==
true,
3509 internal::ExcAccessToUninitializedField());
3511 return BaseClass::get_hessian(q_point);
3516 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3521 const unsigned int q_point)
3523 BaseClass::submit_gradient (grad_in, q_point);
3528 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3534 const unsigned int q_point)
3536 BaseClass::submit_gradient(grad_in, q_point);
3539 template <
int dim,
int dofs_per_cell_,
int n_q_points_,
3545 const unsigned int q_point)
3550 this->gradients_quad_submitted =
true;
3552 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
3555 this->quadrature_weights[q_point] * div_in;
3556 for (
unsigned int d=0; d<dim; ++d)
3558 this->gradients_quad[d][d][q_point] = (fac *
3559 this->cartesian_data[0][d]);
3560 for (
unsigned int e=d+1; e<dim; ++e)
3570 this->cell_type == internal::MatrixFreeFunctions::general ?
3571 this->jacobian[q_point] : this->jacobian[0];
3573 (this->cell_type == internal::MatrixFreeFunctions::general ?
3574 this->J_value[q_point] : this->J_value[0] *
3575 this->quadrature_weights[q_point]) * div_in;
3576 for (
unsigned int d=0; d<dim; ++d)
3578 for (
unsigned int e=0; e<dim; ++e)
3579 this->gradients_quad[d][e][q_point] = jac[d][e] * fac;
3586 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3592 const unsigned int q_point)
3600 this->gradients_quad_submitted =
true;
3602 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
3605 for (
unsigned int d=0; d<dim; ++d)
3606 this->gradients_quad[d][d][q_point] = (sym_grad.access_raw_entry(d) *
3608 this->cartesian_data[0][d]);
3609 for (
unsigned int e=0, counter=dim; e<dim; ++e)
3610 for (
unsigned int d=e+1; d<dim; ++d, ++counter)
3613 this->gradients_quad[e][d][q_point] = (value *
3614 this->cartesian_data[0][d]);
3615 this->gradients_quad[d][e][q_point] = (value *
3616 this->cartesian_data[0][e]);
3623 this->cell_type == internal::MatrixFreeFunctions::general ?
3624 this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
3626 this->cell_type == internal::MatrixFreeFunctions::general ?
3627 this->jacobian[q_point] : this->jacobian[0];
3629 for (
unsigned int i=0; i<dim; ++i)
3630 weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
3631 for (
unsigned int i=0, counter=dim; i<dim; ++i)
3632 for (
unsigned int j=i+1; j<dim; ++j, ++counter)
3635 weighted[i][j] = value;
3636 weighted[j][i] = value;
3638 for (
unsigned int comp=0; comp<dim; ++comp)
3639 for (
unsigned int d=0; d<dim; ++d)
3642 for (
unsigned int e=1; e<dim; ++e)
3643 new_val += jac[e][d] * weighted[comp][e];
3644 this->gradients_quad[comp][d][q_point] = new_val;
3651 template <
int dim,
int dofs_per_cell_,
int n_q_po
ints_,
typename Number>
3656 const unsigned int q_point)
3663 ExcMessage(
"Testing by the curl in 1d is not a useful operation"));
3666 grad[1][0] = curl[0];
3667 grad[0][1] = -curl[0];
3670 grad[2][1] = curl[0];
3671 grad[1][2] = -curl[0];
3672 grad[0][2] = curl[1];
3673 grad[2][0] = -curl[1];
3674 grad[1][0] = curl[2];
3675 grad[0][1] = -curl[2];
3680 submit_gradient (grad, q_point);
3687 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
3692 const unsigned int fe_no,
3693 const unsigned int quad_no)
3695 BaseClass (data_in, fe_no, quad_no)
3701 n_q_points != this->data.n_q_points)
3703 std::string message =
3704 "-------------------------------------------------------\n";
3705 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
3706 message +=
" Called --> FEEvaluation<dim,";
3710 message +=
",Number>(data, ";
3718 if (dofs_per_cell == this->matrix_info.
get_dof_info(fe_no).dofs_per_cell[this->active_fe_index])
3719 proposed_dof_comp = fe_no;
3721 for (
unsigned int no=0; no<this->matrix_info.
n_components(); ++no)
3722 if (this->matrix_info.
get_dof_info(no).dofs_per_cell[this->active_fe_index]
3725 proposed_dof_comp = no;
3729 this->mapping_info.
mapping_data_gen[quad_no].n_q_points[this->active_quad_index])
3730 proposed_quad_comp = quad_no;
3732 for (
unsigned int no=0; no<this->mapping_info.
mapping_data_gen.size(); ++no)
3733 if (this->mapping_info.
mapping_data_gen[no].n_q_points[this->active_quad_index]
3736 proposed_quad_comp = no;
3742 if (proposed_dof_comp != fe_no)
3743 message +=
"Wrong vector component selection:\n";
3745 message +=
"Wrong quadrature formula selection:\n";
3746 message +=
" Did you mean FEEvaluation<dim,";
3750 message +=
",Number>(data, ";
3753 std::string correct_pos;
3754 if (proposed_dof_comp != fe_no)
3755 correct_pos =
" ^ ";
3758 if (proposed_quad_comp != quad_no)
3759 correct_pos +=
" ^\n";
3761 correct_pos +=
" \n";
3762 message +=
" " + correct_pos;
3766 const unsigned int proposed_fe_degree =
static_cast<unsigned int>(std::pow(1.001*this->data.
dofs_per_cell,1./dim))-1;
3767 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(std::pow(1.001*this->data.
n_q_points,1./dim));
3768 message +=
"Wrong template arguments:\n";
3769 message +=
" Did you mean FEEvaluation<dim,";
3773 message +=
",Number>(data, ";
3776 std::string correct_pos;
3777 if (proposed_fe_degree != fe_degree)
3781 if (proposed_n_q_points_1d != n_q_points_1d)
3782 correct_pos +=
" ^\n";
3784 correct_pos +=
" \n";
3785 message +=
" " + correct_pos;
3788 n_q_points == this->data.n_q_points,
3793 n_q_points[this->active_quad_index]);
3795 this->dof_info.dofs_per_cell[this->active_fe_index]);
3807 template <
int dim,
int fe_degree,
int n_q_points_1d,
typename Number,
3808 int direction,
bool dof_to_quad,
bool add>
3811 apply_tensor_product (
const Number *shape_data,
3816 const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
3817 nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
3819 const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
3820 const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
3823 for (
int i2=0; i2<n_blocks2; ++i2)
3825 for (
int i1=0; i1<n_blocks1; ++i1)
3827 for (
int col=0; col<nn; ++col)
3830 if (dof_to_quad ==
true)
3831 val0 = shape_data[col];
3833 val0 = shape_data[col*n_q_points_1d];
3834 Number res0 = val0 * in[0];
3835 for (
int ind=1; ind<mm; ++ind)
3837 if (dof_to_quad ==
true)
3838 val0 = shape_data[ind*n_q_points_1d+col];
3840 val0 = shape_data[col*n_q_points_1d+ind];
3841 res0 += val0 * in[stride*ind];
3844 out[stride*col] = res0;
3846 out[stride*col] += res0;
3882 template <
int dim,
int fe_degree,
typename Number,
int face_direction,
3883 bool dof_to_quad,
bool add>
3886 apply_tensor_product_face (
const Number *shape_data,
3890 const int n_blocks1 = dim > 1 ? (fe_degree+1) : 1;
3891 const int n_blocks2 = dim > 2 ? (fe_degree+1) : 1;
3894 const int mm = dof_to_quad ? (fe_degree+1) : 1,
3895 nn = dof_to_quad ? 1 : (fe_degree+1);
3899 for (
int i2=0; i2<n_blocks2; ++i2)
3901 for (
int i1=0; i1<n_blocks1; ++i1)
3903 if (dof_to_quad ==
true)
3905 Number res0 = shape_data[0] * in[0];
3906 for (
int ind=1; ind<mm; ++ind)
3907 res0 += shape_data[ind] * in[stride*ind];
3915 for (
int col=0; col<nn; ++col)
3917 out[col*stride] = shape_data[col] * in[0];
3919 out[col*stride] += shape_data[col] * in[0];
3925 switch (face_direction)
3940 if (face_direction == 1)
3971 template <
int dim,
int fe_degree,
int n_q_points_1d,
typename Number,
3972 int direction,
bool dof_to_quad,
bool add>
3975 apply_tensor_product_values (
const Number *shape_values,
3980 const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
3981 nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
3982 const int n_cols = nn / 2;
3983 const int mid = mm / 2;
3985 const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
3986 const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
3989 for (
int i2=0; i2<n_blocks2; ++i2)
3991 for (
int i1=0; i1<n_blocks1; ++i1)
3993 for (
int col=0; col<n_cols; ++col)
3995 Number val0, val1, in0, in1, res0, res1;
3996 if (dof_to_quad ==
true)
3998 val0 = shape_values[col];
3999 val1 = shape_values[nn-1-col];
4003 val0 = shape_values[col*n_q_points_1d];
4004 val1 = shape_values[(col+1)*n_q_points_1d-1];
4009 in1 = in[stride*(mm-1)];
4014 for (
int ind=1; ind<mid; ++ind)
4016 if (dof_to_quad ==
true)
4018 val0 = shape_values[ind*n_q_points_1d+col];
4019 val1 = shape_values[ind*n_q_points_1d+nn-1-col];
4023 val0 = shape_values[col*n_q_points_1d+ind];
4024 val1 = shape_values[(col+1)*n_q_points_1d-1-ind];
4026 in0 = in[stride*ind];
4027 in1 = in[stride*(mm-1-ind)];
4035 res0 = res1 = Number();
4036 if (dof_to_quad ==
true)
4040 val0 = shape_values[mid*n_q_points_1d+col];
4041 val1 = val0 * in[stride*mid];
4048 if (mm % 2 == 1 && nn % 2 == 0)
4050 val0 = shape_values[col*n_q_points_1d+mid];
4051 val1 = val0 * in[stride*mid];
4058 out[stride*col] = res0;
4059 out[stride*(nn-1-col)] = res1;
4063 out[stride*col] += res0;
4064 out[stride*(nn-1-col)] += res1;
4067 if ( dof_to_quad ==
true && nn%2==1 && mm%2==1 )
4070 out[stride*n_cols] = in[stride*mid];
4072 out[stride*n_cols] += in[stride*mid];
4074 else if (dof_to_quad ==
true && nn%2==1)
4077 Number val0 = shape_values[n_cols];
4080 res0 = in[0] + in[stride*(mm-1)];
4082 for (
int ind=1; ind<mid; ++ind)
4084 val0 = shape_values[ind*n_q_points_1d+n_cols];
4085 Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
4093 out[stride*n_cols] = res0;
4095 out[stride*n_cols] += res0;
4097 else if (dof_to_quad ==
false && nn%2 == 1)
4102 Number val0 = shape_values[n_cols*n_q_points_1d];
4103 res0 = in[0] + in[stride*(mm-1)];
4105 for (
int ind=1; ind<mid; ++ind)
4107 val0 = shape_values[n_cols*n_q_points_1d+ind];
4108 Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
4113 res0 += in[stride*mid];
4118 out[stride*n_cols] = res0;
4120 out[stride*n_cols] += res0;
4174 template <
int dim,
int fe_degree,
int n_q_points_1d,
typename Number,
4175 int direction,
bool dof_to_quad,
bool add>
4178 apply_tensor_product_gradients (
const Number *shape_gradients,
4183 const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
4184 nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
4185 const int n_cols = nn / 2;
4186 const int mid = mm / 2;
4188 const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
4189 const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
4192 for (
int i2=0; i2<n_blocks2; ++i2)
4194 for (
int i1=0; i1<n_blocks1; ++i1)
4196 for (
int col=0; col<n_cols; ++col)
4198 Number val0, val1, in0, in1, res0, res1;
4199 if (dof_to_quad ==
true)
4201 val0 = shape_gradients[col];
4202 val1 = shape_gradients[nn-1-col];
4206 val0 = shape_gradients[col*n_q_points_1d];
4207 val1 = shape_gradients[(nn-col-1)*n_q_points_1d];
4212 in1 = in[stride*(mm-1)];
4217 for (
int ind=1; ind<mid; ++ind)
4219 if (dof_to_quad ==
true)
4221 val0 = shape_gradients[ind*n_q_points_1d+col];
4222 val1 = shape_gradients[ind*n_q_points_1d+nn-1-col];
4226 val0 = shape_gradients[col*n_q_points_1d+ind];
4227 val1 = shape_gradients[(nn-col-1)*n_q_points_1d+ind];
4229 in0 = in[stride*ind];
4230 in1 = in[stride*(mm-1-ind)];
4238 res0 = res1 = Number();
4241 if (dof_to_quad ==
true)
4242 val0 = shape_gradients[mid*n_q_points_1d+col];
4244 val0 = shape_gradients[col*n_q_points_1d+mid];
4245 val1 = val0 * in[stride*mid];
4251 out[stride*col] = res0;
4252 out[stride*(nn-1-col)] = res1;
4256 out[stride*col] += res0;
4257 out[stride*(nn-1-col)] += res1;
4263 if (dof_to_quad ==
true)
4264 val0 = shape_gradients[n_cols];
4266 val0 = shape_gradients[n_cols*n_q_points_1d];
4267 res0 = in[0] - in[stride*(mm-1)];
4269 for (
int ind=1; ind<mid; ++ind)
4271 if (dof_to_quad ==
true)
4272 val0 = shape_gradients[ind*n_q_points_1d+n_cols];
4274 val0 = shape_gradients[n_cols*n_q_points_1d+ind];
4275 Number val1 = in[stride*ind] - in[stride*(mm-1-ind)];
4280 out[stride*n_cols] = res0;
4282 out[stride*n_cols] += res0;
4318 template <
int dim,
int fe_degree,
int n_q_points_1d,
typename Number,
4319 int direction,
bool dof_to_quad,
bool add>
4322 apply_tensor_product_hessians (
const Number *shape_hessians,
4327 const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
4328 nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
4329 const int n_cols = nn / 2;
4330 const int mid = mm / 2;
4332 const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
4333 const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
4336 for (
int i2=0; i2<n_blocks2; ++i2)
4338 for (
int i1=0; i1<n_blocks1; ++i1)
4340 for (
int col=0; col<n_cols; ++col)
4342 Number val0, val1, in0, in1, res0, res1;
4343 if (dof_to_quad ==
true)
4345 val0 = shape_hessians[col];
4346 val1 = shape_hessians[nn-1-col];
4350 val0 = shape_hessians[col*n_q_points_1d];
4351 val1 = shape_hessians[(col+1)*n_q_points_1d-1];
4356 in1 = in[stride*(mm-1)];
4361 for (
int ind=1; ind<mid; ++ind)
4363 if (dof_to_quad ==
true)
4365 val0 = shape_hessians[ind*n_q_points_1d+col];
4366 val1 = shape_hessians[ind*n_q_points_1d+nn-1-col];
4370 val0 = shape_hessians[col*n_q_points_1d+ind];
4371 val1 = shape_hessians[(col+1)*n_q_points_1d-1-ind];
4373 in0 = in[stride*ind];
4374 in1 = in[stride*(mm-1-ind)];
4382 res0 = res1 = Number();
4385 if (dof_to_quad ==
true)
4386 val0 = shape_hessians[mid*n_q_points_1d+col];
4388 val0 = shape_hessians[col*n_q_points_1d+mid];
4389 val1 = val0 * in[stride*mid];
4395 out[stride*col] = res0;
4396 out[stride*(nn-1-col)] = res1;
4400 out[stride*col] += res0;
4401 out[stride*(nn-1-col)] += res1;
4407 if (dof_to_quad ==
true)
4408 val0 = shape_hessians[n_cols];
4410 val0 = shape_hessians[n_cols*n_q_points_1d];
4413 res0 = in[0] + in[stride*(mm-1)];
4415 for (
int ind=1; ind<mid; ++ind)
4417 if (dof_to_quad ==
true)
4418 val0 = shape_hessians[ind*n_q_points_1d+n_cols];
4420 val0 = shape_hessians[n_cols*n_q_points_1d+ind];
4421 Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
4430 if (dof_to_quad ==
true)
4431 val0 = shape_hessians[mid*n_q_points_1d+n_cols];
4433 val0 = shape_hessians[n_cols*n_q_points_1d+mid];
4434 res0 += val0 * in[stride*mid];
4437 out[stride*n_cols] = res0;
4439 out[stride*n_cols] += res0;
4483 template <
int dim,
int fe_degree,
int n_q_points_1d,
typename Number,
4484 int direction,
bool dof_to_quad,
bool add,
int type>
4487 apply_tensor_product_evenodd (
const Number shapes [][(n_q_points_1d+1)/2],
4493 const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
4494 nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
4495 const int n_cols = nn / 2;
4496 const int mid = mm / 2;
4498 const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
4499 const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
4506 for (
int i2=0; i2<n_blocks2; ++i2)
4508 for (
int i1=0; i1<n_blocks1; ++i1)
4510 Number xp[mid], xm[mid];
4511 for (
int i=0; i<mid; ++i)
4513 if (dof_to_quad ==
true && type == 1)
4515 xp[i] = in[stride*i] - in[stride*(mm-1-i)];
4516 xm[i] = in[stride*i] + in[stride*(mm-1-i)];
4520 xp[i] = in[stride*i] + in[stride*(mm-1-i)];
4521 xm[i] = in[stride*i] - in[stride*(mm-1-i)];
4524 for (
int col=0; col<n_cols; ++col)
4529 if (dof_to_quad ==
true)
4531 r0 = shapes[0][col] * xp[0];
4532 r1 = shapes[fe_degree][col] * xm[0];
4536 r0 = shapes[col][0] * xp[0];
4537 r1 = shapes[fe_degree-col][0] * xm[0];
4539 for (
int ind=1; ind<mid; ++ind)
4541 if (dof_to_quad ==
true)
4543 r0 += shapes[ind][col] * xp[ind];
4544 r1 += shapes[fe_degree-ind][col] * xm[ind];
4548 r0 += shapes[col][ind] * xp[ind];
4549 r1 += shapes[fe_degree-col][ind] * xm[ind];
4555 if (mm % 2 == 1 && dof_to_quad ==
true)
4558 r1 += shapes[mid][col] * in[stride*mid];
4560 r0 += shapes[mid][col] * in[stride*mid];
4562 else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0))
4563 r0 += shapes[col][mid] * in[stride*mid];
4567 out[stride*col] = r0 + r1;
4568 if (type == 1 && dof_to_quad ==
false)
4569 out[stride*(nn-1-col)] = r1 - r0;
4571 out[stride*(nn-1-col)] = r0 - r1;
4575 out[stride*col] += r0 + r1;
4576 if (type == 1 && dof_to_quad ==
false)
4577 out[stride*(nn-1-col)] += r1 - r0;
4579 out[stride*(nn-1-col)] += r0 - r1;
4582 if ( type == 0 && dof_to_quad ==
true && nn%2==1 && mm%2==1 )
4585 out[stride*n_cols] = in[stride*mid];
4587 out[stride*n_cols] += in[stride*mid];
4589 else if (dof_to_quad ==
true && nn%2==1)
4594 r0 = shapes[0][n_cols] * xp[0];
4595 for (
int ind=1; ind<mid; ++ind)
4596 r0 += shapes[ind][n_cols] * xp[ind];
4600 if (type != 1 && mm % 2 == 1)
4601 r0 += shapes[mid][n_cols] * in[stride*mid];
4604 out[stride*n_cols] = r0;
4606 out[stride*n_cols] += r0;
4608 else if (dof_to_quad ==
false && nn%2 == 1)
4615 r0 = shapes[n_cols][0] * xm[0];
4616 for (
int ind=1; ind<mid; ++ind)
4617 r0 += shapes[n_cols][ind] * xm[ind];
4621 r0 = shapes[n_cols][0] * xp[0];
4622 for (
int ind=1; ind<mid; ++ind)
4623 r0 += shapes[n_cols][ind] * xp[ind];
4629 if (type == 0 && mm % 2 == 1)
4630 r0 += in[stride*mid];
4631 else if (type == 2 && mm % 2 == 1)
4632 r0 += shapes[n_cols][mid] * in[stride*mid];
4635 out[stride*n_cols] = r0;
4637 out[stride*n_cols] += r0;
4690 template <
int dim,
int fe_degree,
typename Number,
4691 int direction,
bool dof_to_quad,
bool add>
4694 apply_tensor_product_gradients_gl (
const Number *shape_gradients,
4699 const int mm = fe_degree+1;
4700 const int nn = fe_degree+1;
4701 const int n_cols = nn / 2;
4702 const int mid = mm / 2;
4704 const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
4705 const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
4708 for (
int i2=0; i2<n_blocks2; ++i2)
4710 for (
int i1=0; i1<n_blocks1; ++i1)
4712 for (
int col=0; col<n_cols; ++col)
4714 Number val0, val1, in0, in1, res0, res1;
4717 if (dof_to_quad ==
true)
4719 val0 = shape_gradients[col];
4720 val1 = shape_gradients[nn-1-col];
4724 val0 = shape_gradients[col*mm];
4725 val1 = shape_gradients[(nn-col-1)*mm];
4728 in1 = in[stride*(mm-1)];
4731 if ((mm+dof_to_quad)%2 == 1)
4753 for (
int ind=1; ind<mid; ++ind)
4755 if (dof_to_quad ==
true)
4757 val0 = shape_gradients[ind*mm+col];
4758 val1 = shape_gradients[ind*mm+nn-1-col];
4762 val0 = shape_gradients[col*mm+ind];
4763 val1 = shape_gradients[(nn-col-1)*mm+ind];
4767 in0 = in[stride*ind];
4768 in1 = in[stride*(mm-1-ind)];
4784 res0 = res1 = Number();
4787 if (dof_to_quad ==
true)
4788 val0 = shape_gradients[mid*mm+col];
4790 val0 = shape_gradients[col*mm+mid];
4791 val1 = val0 * in[stride*mid];
4797 out[stride*col] = res0;
4798 out[stride*(nn-1-col)] = res1;
4802 out[stride*col] += res0;
4803 out[stride*(nn-1-col)] += res1;
4809 if (dof_to_quad ==
true)
4810 val0 = shape_gradients[n_cols];
4812 val0 = shape_gradients[n_cols*mm];
4815 res0 = in[0] - in[stride*(mm-1)];
4817 for (
int ind=1; ind<mid; ++ind)
4819 if (dof_to_quad ==
true)
4820 val0 = shape_gradients[ind*mm+n_cols];
4822 val0 = shape_gradients[n_cols*mm+ind];
4823 Number val1 = in[stride*ind] - in[stride*(mm-1-ind)];
4831 out[stride*n_cols] = res0;
4833 out[stride*n_cols] += res0;
4870 template <
typename FEEval>
4873 do_evaluate (FEEval &fe_eval,
4874 const bool evaluate_val,
4875 const bool evaluate_grad,
4876 const bool evaluate_lapl)
4880 Assert (fe_eval.dof_values_initialized ==
true,
4881 internal::ExcAccessToUninitializedField());
4883 const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ?
4884 FEEval::dofs_per_cell : FEEval::n_q_points;
4885 const unsigned int n_components = FEEval::n_components;
4886 const unsigned int dim = FEEval::dimension;
4888 for (
unsigned int c=0; c<n_components; c++)
4897 if (evaluate_grad ==
true)
4900 fe_eval.template apply_gradients<0,true,false>
4901 (fe_eval.values_dofs[c], temp1);
4902 fe_eval.template apply_values<1,true,false>
4904 fe_eval.template apply_values<2,true,false>
4905 (temp2, fe_eval.gradients_quad[c][0]);
4908 if (evaluate_lapl ==
true)
4911 if (evaluate_grad ==
false)
4913 fe_eval.template apply_gradients<0,true,false>
4914 (fe_eval.values_dofs[c], temp1);
4915 fe_eval.template apply_values<1,true,false>
4918 fe_eval.template apply_gradients<2,true,false>
4919 (temp2, fe_eval.hessians_quad[c][4]);
4922 fe_eval.template apply_gradients<1,true,false>
4924 fe_eval.template apply_values<2,true,false>
4925 (temp2, fe_eval.hessians_quad[c][3]);
4928 fe_eval.template apply_hessians<0,true,false>
4929 (fe_eval.values_dofs[c], temp1);
4930 fe_eval.template apply_values<1,true,false>
4932 fe_eval.template apply_values<2,true,false>
4933 (temp2, fe_eval.hessians_quad[c][0]);
4937 fe_eval.template apply_values<0,true,false>
4938 (fe_eval.values_dofs[c], temp1);
4939 if (evaluate_grad ==
true)
4941 fe_eval.template apply_gradients<1,true,false>
4943 fe_eval.template apply_values<2,true,false>
4944 (temp2, fe_eval.gradients_quad[c][1]);
4947 if (evaluate_lapl ==
true)
4950 if (evaluate_grad ==
false)
4951 fe_eval.template apply_gradients<1,true,false>
4953 fe_eval.template apply_gradients<2,true,false>
4954 (temp2, fe_eval.hessians_quad[c][5]);
4957 fe_eval.template apply_hessians<1,true,false>
4959 fe_eval.template apply_values<2,true,false>
4960 (temp2, fe_eval.hessians_quad[c][1]);
4964 fe_eval.template apply_values<1,true,false>
4966 if (evaluate_grad ==
true)
4967 fe_eval.template apply_gradients<2,true,false>
4968 (temp2, fe_eval.gradients_quad[c][2]);
4972 if (evaluate_lapl ==
true)
4973 fe_eval.template apply_hessians<2,true,false>
4974 (temp2, fe_eval.hessians_quad[c][2]);
4977 if (evaluate_val ==
true)
4978 fe_eval.template apply_values<2,true,false>
4979 (temp2, fe_eval.values_quad[c]);
4986 if (evaluate_grad ==
true)
4988 fe_eval.template apply_gradients<0,true,false>
4989 (fe_eval.values_dofs[c], temp1);
4990 fe_eval.template apply_values<1,true,false>
4991 (temp1, fe_eval.gradients_quad[c][0]);
4993 if (evaluate_lapl ==
true)
4996 if (evaluate_grad ==
false)
4997 fe_eval.template apply_gradients<0,true,false>
4998 (fe_eval.values_dofs[c], temp1);
4999 fe_eval.template apply_gradients<1,true,false>
5000 (temp1, fe_eval.hessians_quad[c][2]);
5003 fe_eval.template apply_hessians<0,true,false>
5004 (fe_eval.values_dofs[c], temp1);
5005 fe_eval.template apply_values<1,true,false>
5006 (temp1, fe_eval.hessians_quad[c][0]);
5010 fe_eval.template apply_values<0,true,false>
5011 (fe_eval.values_dofs[c], temp1);
5012 if (evaluate_grad ==
true)
5013 fe_eval.template apply_gradients<1,true,false>
5014 (temp1, fe_eval.gradients_quad[c][1]);
5017 if (evaluate_lapl ==
true)
5018 fe_eval.template apply_hessians<1,true,false>
5019 (temp1, fe_eval.hessians_quad[c][1]);
5022 if (evaluate_val ==
true)
5023 fe_eval.template apply_values<1,true,false>
5024 (temp1, fe_eval.values_quad[c]);
5029 if (evaluate_val ==
true)
5030 fe_eval.template apply_values<0,true,false>
5031 (fe_eval.values_dofs[c], fe_eval.values_quad[c]);
5032 if (evaluate_grad ==
true)
5033 fe_eval.template apply_gradients<0,true,false>
5034 (fe_eval.values_dofs[c], fe_eval.gradients_quad[c][0]);
5035 if (evaluate_lapl ==
true)
5036 fe_eval.template apply_hessians<0,true,false>
5037 (fe_eval.values_dofs[c], fe_eval.hessians_quad[c][0]);
5046 if (evaluate_val ==
true)
5047 fe_eval.values_quad_initialized =
true;
5048 if (evaluate_grad ==
true)
5049 fe_eval.gradients_quad_initialized =
true;
5050 if (evaluate_lapl ==
true)
5051 fe_eval.hessians_quad_initialized =
true;
5057 template <
typename FEEval>
5060 do_integrate (FEEval &fe_eval,
5061 const bool integrate_val,
5062 const bool integrate_grad)
5065 if (integrate_val ==
true)
5066 Assert (fe_eval.values_quad_submitted ==
true,
5067 ExcAccessToUninitializedField());
5068 if (integrate_grad ==
true)
5069 Assert (fe_eval.gradients_quad_submitted ==
true,
5070 ExcAccessToUninitializedField());
5072 const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ?
5073 FEEval::dofs_per_cell : FEEval::n_q_points;
5074 const unsigned int n_components = FEEval::n_components;
5075 const unsigned int dim = FEEval::dimension;
5078 for (
unsigned int c=0; c<n_components; c++)
5087 if (integrate_val ==
true)
5090 fe_eval.template apply_values<0,false,false>
5091 (fe_eval.values_quad[c], temp1);
5093 if (integrate_grad ==
true)
5096 if (integrate_val ==
true)
5097 fe_eval.template apply_gradients<0,false,true>
5098 (fe_eval.gradients_quad[c][0], temp1);
5100 fe_eval.template apply_gradients<0,false,false>
5101 (fe_eval.gradients_quad[c][0], temp1);
5103 fe_eval.template apply_values<1,false,false>
5105 if (integrate_grad ==
true)
5108 fe_eval.template apply_values<0,false,false>
5109 (fe_eval.gradients_quad[c][1], temp1);
5110 fe_eval.template apply_gradients<1,false,true>
5113 fe_eval.template apply_values<2,false,false>
5114 (temp2, fe_eval.values_dofs[c]);
5115 if (integrate_grad ==
true)
5118 fe_eval.template apply_values<0,false,false>
5119 (fe_eval.gradients_quad[c][2], temp1);
5120 fe_eval.template apply_values<1,false,false>
5122 fe_eval.template apply_gradients<2,false,true>
5123 (temp2, fe_eval.values_dofs[c]);
5131 if (integrate_val ==
true)
5132 fe_eval.template apply_values<0,false,false>
5133 (fe_eval.values_quad[c], temp1);
5134 if (integrate_grad ==
true)
5137 if (integrate_val ==
true)
5138 fe_eval.template apply_gradients<0,false,true>
5139 (fe_eval.gradients_quad[c][0], temp1);
5141 fe_eval.template apply_gradients<0,false,false>
5142 (fe_eval.gradients_quad[c][0], temp1);
5144 fe_eval.template apply_values<1,false,false>
5145 (temp1, fe_eval.values_dofs[c]);
5146 if (integrate_grad ==
true)
5149 fe_eval.template apply_values<0,false,false>
5150 (fe_eval.gradients_quad[c][1], temp1);
5151 fe_eval.template apply_gradients<1,false,true>
5152 (temp1, fe_eval.values_dofs[c]);
5159 if (integrate_grad ==
true)
5160 fe_eval.template apply_gradients<0,false,false>
5161 (fe_eval.gradients_quad[c][0], fe_eval.values_dofs[c]);
5162 if (integrate_val ==
true)
5164 if (integrate_grad ==
true)
5165 fe_eval.template apply_values<0,false,true>
5166 (fe_eval.values_quad[c], fe_eval.values_dofs[c]);
5168 fe_eval.template apply_values<0,false,false>
5169 (fe_eval.values_quad[c], fe_eval.values_dofs[c]);
5179 fe_eval.dof_values_initialized =
true;
5187 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
5193 const bool evaluate_grad,
5194 const bool evaluate_lapl)
5196 internal::do_evaluate (*
this, evaluate_val, evaluate_grad, evaluate_lapl);
5201 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
5207 const bool integrate_grad)
5209 internal::do_integrate (*
this, integrate_val, integrate_grad);
5214 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
5228 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5234 return this->quadrature_points[q];
5236 point[0] = this->quadrature_points[q%n_q_points_1d][0];
5237 point[1] = this->quadrature_points[q/n_q_points_1d][1];
5240 point[0] = this->quadrature_points[q%n_q_points_1d][0];
5241 point[1] = this->quadrature_points[(q/n_q_points_1d)%n_q_points_1d][1];
5242 point[2] = this->quadrature_points[q/(n_q_points_1d*n_q_points_1d)][2];
5251 return this->quadrature_points[q];
5256 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
5258 template <
int direction,
bool dof_to_quad,
bool add>
5265 internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,
5272 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
5274 template <
int direction,
bool dof_to_quad,
bool add>
5279 VectorizedArray<Number> out [])
5281 internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,
5282 VectorizedArray<Number>, direction, dof_to_quad, add>
5288 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
5290 template <
int direction,
bool dof_to_quad,
bool add>
5295 VectorizedArray<Number> out [])
5297 internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,
5298 VectorizedArray<Number>, direction, dof_to_quad, add>
5306 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
5311 const unsigned int fe_no,
5312 const unsigned int quad_no)
5314 BaseClass (data_in, fe_no, quad_no)
5318 const double zero_tol =
5320 std::string error_message =
"FEEvaluation not appropriate.\n";
5321 error_message +=
" It assumes symmetry of quadrature points w.r.t. 0.5 \n";
5322 error_message +=
" and the basis functions starting from left and right.\n";
5323 error_message +=
"Try FEEvaluationGeneral<...> instead!";
5326 const unsigned int n_dofs_1d = fe_degree + 1;
5327 for (
unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
5328 for (
unsigned int j=0; j<n_q_points_1d; ++j)
5330 this->data.shape_values[(n_dofs_1d-i)*n_q_points_1d
5331 -j-1][0]) < zero_tol,
5336 if (n_q_points_1d%2 == 1 && n_dofs_1d%2 == 1)
5338 for (
int i=0; i<static_cast<int>(n_dofs_1d/2); ++i)
5340 n_q_points_1d/2][0]) < zero_tol,
5343 n_q_points_1d/2][0]-1.)< zero_tol,
5349 for (
unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
5350 for (
unsigned int j=0; j<n_q_points_1d; ++j)
5352 this->data.shape_gradients[(n_dofs_1d-i)*n_q_points_1d-
5353 j-1][0]) < zero_tol,
5355 if (n_dofs_1d%2 == 1 && n_q_points_1d%2 == 1)
5357 (n_q_points_1d/2)][0]) < zero_tol,
5362 for (
unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
5363 for (
unsigned int j=0; j<n_q_points_1d; ++j)
5365 this->data.shape_hessians[(n_dofs_1d-i)*n_q_points_1d-
5366 j-1][0]) < zero_tol,
5372 for (
unsigned int i=0; i<(fe_degree+1)/2; ++i)
5373 for (
unsigned int q=0; q<(n_q_points_1d+1)/2; ++q)
5375 shape_val_evenodd[i][q] =
5377 this->data.
shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
5378 shape_val_evenodd[fe_degree-i][q] =
5380 this->data.
shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
5382 shape_gra_evenodd[i][q] =
5385 shape_gra_evenodd[fe_degree-i][q] =
5389 shape_hes_evenodd[i][q] =
5392 shape_hes_evenodd[fe_degree-i][q] =
5396 if (fe_degree % 2 == 0)
5397 for (
unsigned int q=0; q<(n_q_points_1d+1)/2; ++q)
5399 shape_val_evenodd[fe_degree/2][q] =
5400 this->data.
shape_values[(fe_degree/2)*n_q_points_1d+q];
5401 shape_gra_evenodd[fe_degree/2][q] =
5403 shape_hes_evenodd[fe_degree/2][q] =
5410 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
5416 const bool evaluate_grad,
5417 const bool evaluate_lapl)
5419 internal::do_evaluate (*
this, evaluate_val, evaluate_grad, evaluate_lapl);
5424 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
5429 ::integrate (
bool integrate_val,
bool integrate_grad)
5431 internal::do_integrate (*
this, integrate_val, integrate_grad);
5436 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
5438 template <
int direction,
bool dof_to_quad,
bool add>
5443 VectorizedArray<Number> out [])
5448 if (fe_degree > 1 || n_q_points_1d > 3)
5449 internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
5450 VectorizedArray<Number>, direction, dof_to_quad, add, 0>
5451 (shape_val_evenodd, in, out);
5453 internal::apply_tensor_product_values<dim,fe_degree,n_q_points_1d,
5454 VectorizedArray<Number>, direction, dof_to_quad, add>
5460 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
5462 template <
int direction,
bool dof_to_quad,
bool add>
5467 VectorizedArray<Number> out [])
5469 if (fe_degree > 1 || n_q_points_1d > 3)
5470 internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
5471 VectorizedArray<Number>, direction, dof_to_quad, add, 1>
5472 (shape_gra_evenodd, in, out);
5474 internal::apply_tensor_product_gradients<dim,fe_degree,n_q_points_1d,
5475 VectorizedArray<Number>, direction, dof_to_quad, add>
5484 template <
int dim,
int fe_degree,
int n_q_points_1d,
int n_components_,
5486 template <
int direction,
bool dof_to_quad,
bool add>
5491 VectorizedArray<Number> out [])
5493 if (fe_degree > 1 || n_q_points_1d > 3)
5494 internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
5495 VectorizedArray<Number>, direction, dof_to_quad, add, 2>
5496 (shape_hes_evenodd, in, out);
5498 internal::apply_tensor_product_hessians<dim,fe_degree,n_q_points_1d,
5499 VectorizedArray<Number>, direction, dof_to_quad, add>
5507 template <
int dim,
int fe_degree,
int n_components_,
typename Number>
5511 const unsigned int fe_no,
5512 const unsigned int quad_no)
5514 BaseClass (data_in, fe_no, quad_no)
5517 std::string error_mess =
"FEEvaluationGL not appropriate. It assumes:\n";
5518 error_mess +=
" - identity operation for shape values\n";
5519 error_mess +=
" - zero diagonal at interior points for gradients\n";
5520 error_mess +=
" - gradient equal to unity at element boundary\n";
5521 error_mess +=
"Try FEEvaluation<...> instead!";
5523 const double zero_tol =
5526 const unsigned int n_points_1d = fe_degree+1;
5527 for (
unsigned int i=0; i<n_points_1d; ++i)
5528 for (
unsigned int j=0; j<n_points_1d; ++j)
5540 for (
unsigned int i=1; i<n_points_1d-1; ++i)
5544 (n_points_1d%2==0 ? -1. : 1.)) < zero_tol,
5551 template <
int dim,
int fe_degree,
int n_components_,
typename Number>
5556 const bool evaluate_grad,
5557 const bool evaluate_lapl)
5561 Assert (this->dof_values_initialized ==
true,
5562 internal::ExcAccessToUninitializedField());
5564 if (evaluate_val ==
true)
5566 std::memcpy (&this->values_quad[0][0], &this->values_dofs[0][0],
5567 dofs_per_cell * n_components *
5568 sizeof (this->values_dofs[0][0]));
5570 this->values_quad_initialized =
true;
5575 if (evaluate_grad ==
true)
5582 apply_gradients<0,true,false> (this->values_dofs[comp],
5583 this->gradients_quad[comp][0]);
5585 apply_gradients<1,true,false> (this->values_dofs[comp],
5586 this->gradients_quad[comp][1]);
5588 apply_gradients<2,true,false> (this->values_dofs[comp],
5589 this->gradients_quad[comp][2]);
5594 apply_gradients<0,true,false> (this->values_dofs[comp],
5595 this->gradients_quad[comp][0]);
5597 apply_gradients<1,true,false> (this->values_dofs[comp],
5598 this->gradients_quad[comp][1]);
5601 apply_gradients<0,true,false> (this->values_dofs[comp],
5602 this->gradients_quad[comp][0]);
5605 this->gradients_quad_initialized =
true;
5608 if (evaluate_lapl ==
true)
5615 this->
template apply_hessians<0,true,false> (this->values_dofs[comp],
5616 this->hessians_quad[comp][0]);
5618 this->
template apply_hessians<1,true,false> (this->values_dofs[comp],
5619 this->hessians_quad[comp][1]);
5621 this->
template apply_hessians<2,true,false> (this->values_dofs[comp],
5622 this->hessians_quad[comp][2]);
5624 VectorizedArray<Number> temp1[n_q_points];
5626 apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
5627 apply_gradients<1,true,false> (temp1, this->hessians_quad[comp][3]);
5629 apply_gradients<2,true,false> (temp1, this->hessians_quad[comp][4]);
5631 apply_gradients<1,true,false> (this->values_dofs[comp], temp1);
5632 apply_gradients<2,true,false> (temp1, this->hessians_quad[comp][5]);
5637 this->
template apply_hessians<0,true,false> (this->values_dofs[comp],
5638 this->hessians_quad[comp][0]);
5640 this->
template apply_hessians<1,true,false> (this->values_dofs[comp],
5641 this->hessians_quad[comp][1]);
5642 VectorizedArray<Number> temp1[n_q_points];
5644 apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
5645 apply_gradients<1,true,false> (temp1, this->hessians_quad[comp][2]);
5648 this->
template apply_hessians<0,true,false> (this->values_dofs[comp],
5649 this->hessians_quad[comp][0]);
5652 this->hessians_quad_initialized =
true;
5659 template <
int dim,
int fe_degree,
int n_components_,
typename Number>
5663 ::integrate (
const bool integrate_val,
const bool integrate_grad)
5667 if (integrate_val ==
true)
5668 Assert (this->values_quad_submitted ==
true,
5669 internal::ExcAccessToUninitializedField());
5670 if (integrate_grad ==
true)
5671 Assert (this->gradients_quad_submitted ==
true,
5672 internal::ExcAccessToUninitializedField());
5673 if (integrate_val ==
true)
5674 std::memcpy (&this->values_dofs[0][0], &this->values_quad[0][0],
5675 dofs_per_cell * n_components *
5676 sizeof (this->values_dofs[0][0]));
5677 if (integrate_grad ==
true)
5684 if (integrate_val ==
true)
5685 apply_gradients<0, false, true> (this->gradients_quad[comp][0],
5686 this->values_dofs[comp]);
5688 apply_gradients<0, false, false> (this->gradients_quad[comp][0],
5689 this->values_dofs[comp]);
5692 apply_gradients<1, false, true> (this->gradients_quad[comp][1],
5693 this->values_dofs[comp]);
5696 apply_gradients<2, false, true> (this->gradients_quad[comp][2],
5697 this->values_dofs[comp]);
5702 if (integrate_val ==
true)
5703 apply_gradients<0, false, true> (this->gradients_quad[comp][0],
5704 this->values_dofs[comp]);
5706 apply_gradients<0, false, false> (this->gradients_quad[comp][0],
5707 this->values_dofs[comp]);
5710 apply_gradients<1, false, true> (this->gradients_quad[comp][1],
5711 this->values_dofs[comp]);
5715 if (integrate_val ==
true)
5716 apply_gradients<0, false, true> (this->gradients_quad[comp][0],
5717 this->values_dofs[comp]);
5719 apply_gradients<0, false, false> (this->gradients_quad[comp][0],
5720 this->values_dofs[comp]);
5727 this->dof_values_initialized =
true;
5733 template <
int dim,
int fe_degree,
int n_components_,
typename Number>
5734 template <
int direction,
bool dof_to_quad,
bool add>
5739 VectorizedArray<Number> out [])
5741 internal::apply_tensor_product_gradients_gl<dim,fe_degree,
5742 VectorizedArray<Number>, direction, dof_to_quad, add>
5747 #endif // ifndef DOXYGEN
5750 DEAL_II_NAMESPACE_CLOSE
const Tensor< 2, dim, VectorizedArray< Number > > * jacobian_grad
const unsigned int quad_no
void apply_gradients(const VectorizedArray< Number > in[], VectorizedArray< Number > out[])
bool gradients_quad_submitted
const unsigned int n_fe_components
static const unsigned int invalid_unsigned_int
#define AssertDimension(dim1, dim2)
AlignedVector< std::pair< Tensor< 1, dim, VectorizedArray< Number > >, VectorizedArray< Number > > > cartesian_data
void integrate(const bool integrate_val, const bool integrate_grad)
void reinit(const unsigned int cell)
AlignedVector< std::pair< Tensor< 2, dim, VectorizedArray< Number > >, VectorizedArray< Number > > > affine_data
bool mapping_initialized() const
bool quadrature_points_initialized
VectorizedArray< Number > make_vectorized_array(const Number &u)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
::ExceptionBase & ExcMessage(std::string arg1)
VectorizedArray< Number > hessians_quad[n_components][(dim *(dim+1))/2][n_q_points >0?n_q_points:1]
#define AssertIndexRange(index, range)
internal::MatrixFreeFunctions::CellType cell_type
const VectorizedArray< Number > * begin_gradients() const
const Point< dim, VectorizedArray< Number > > * quadrature_points
value_type get_dof_value(const unsigned int dof) const
std::vector< MappingInfoDependent > mapping_data_gen
const MatrixFree< dim, Number > & matrix_info
void apply_hessians(const VectorizedArray< Number > in[], VectorizedArray< Number > out[])
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
Tensor< 1, n_components_, Tensor< 2, dim, VectorizedArray< Number > > > get_hessian(const unsigned int q_point) const
void evaluate(const bool evaluate_val, const bool evaluate_grad, const bool evaluate_hess=false)
void integrate(const bool integrate_val, const bool integrate_grad)
bool partitioners_are_compatible(const Utilities::MPI::Partitioner &part) const
const Tensor< 1, dim, VectorizedArray< Number > > * cartesian_data
void distribute_local_to_global(VectorType &dst) const
unsigned int cell_data_number
const internal::MatrixFreeFunctions::ShapeInfo< Number > & data
const Tensor< 2, dim, VectorizedArray< Number > > * jacobian
AlignedVector< VectorizedArray< Number > > shape_hessians
void apply_hessians(const VectorizedArray< Number > in[], VectorizedArray< Number > out[])
CellType get_cell_type(const unsigned int cell_chunk_no) const
FEEvaluationGL(const MatrixFree< dim, Number > &matrix_free, const unsigned int fe_no=0, const unsigned int quad_no=0)
FEEvaluationBase(const MatrixFree< dim, Number > &matrix_free, const unsigned int fe_no=0, const unsigned int quad_no=0)
value_type integrate_value() const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number > & mapping_info
#define Assert(cond, exc)
VectorizedArray< Number > gradients_quad[n_components][dim][n_q_points >0?n_q_points:1]
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int fe_component=0) const
unsigned int get_cell_data_number() const
unsigned int n_components() const
void evaluate(const bool evaluate_val, const bool evaluate_grad, const bool evaluate_lapl=false)
Number local_element(const size_type local_index) const
FEEvaluation(const MatrixFree< dim, Number > &matrix_free, const unsigned int fe_no=0, const unsigned int quad_no=0)
const internal::MatrixFreeFunctions::DoFInfo & dof_info
const VectorizedArray< Number > * begin_dof_values() const
bool indices_initialized() const
void submit_dof_value(const value_type val_in, const unsigned int dof)
#define DeclException0(Exception0)
std::string int_to_string(const unsigned int i, const unsigned int digits=numbers::invalid_unsigned_int)
void apply_gradients(const VectorizedArray< Number > in[], VectorizedArray< Number > out[])
VectorizedArray< Number > values_dofs[n_components][dofs_per_cell >0?dofs_per_cell:1]
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const internal::MatrixFreeFunctions::SizeInfo & get_size_info() const
const Tensor< 1,(dim >1?dim *(dim-1)/2:1), Tensor< 1, dim, VectorizedArray< Number > > > * jacobian_grad_upper
FEEvaluationAccess(const MatrixFree< dim, Number > &matrix_free, const unsigned int fe_no=0, const unsigned int quad_no=0)
void read_dof_values(const VectorType &src)
bool values_quad_submitted
void evaluate(const bool evaluate_val, const bool evaluate_grad, const bool evaluate_hess=false)
bool JxW_values_initialized
VectorizedArray< Number > values_quad[n_components][n_q_points >0?n_q_points:1]
const VectorizedArray< Number > * quadrature_weights
const VectorizedArray< Number > * J_value
const Number * constraint_pool_end(const unsigned int pool_index) const
gradient_type get_gradient(const unsigned int q_point) const
void apply_values(const VectorizedArray< Number > in[], VectorizedArray< Number > out[])
unsigned int get_cell_data_index(const unsigned int cell_chunk_no) const
void integrate(const bool integrate_val, const bool integrate_grad)
bool hessians_quad_initialized
bool gradients_quad_initialized
value_type get_value(const unsigned int q_point) const
internal::MatrixFreeFunctions::CellType get_cell_type() const
value_type get_laplacian(const unsigned int q_point) const
void submit_value(const value_type val_in, const unsigned int q_point)
void read_dof_values_plain(const VectorType &src)
unsigned int dofs_per_cell
bool second_derivatives_initialized
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcNotInitialized()
const unsigned int active_fe_index
const VectorizedArray< Number > * begin_values() const
const unsigned int active_quad_index
AlignedVector< VectorizedArray< Number > > shape_values
const VectorizedArray< Number > * begin_hessians() const
const Number * constraint_pool_begin(const unsigned int pool_index) const
AlignedVector< VectorizedArray< Number > > shape_gradients
bool values_quad_initialized
void read_write_operation(const VectorOperation &operation, VectorType *vectors[]) const
void apply_gradients(const VectorizedArray< Number > in[], VectorizedArray< Number > out[])
bool dof_values_initialized
::ExceptionBase & ExcInternalError()
Point< dim, VectorizedArray< Number > > quadrature_point(const unsigned int q_point) const
void set_dof_values(VectorType &dst) const
void apply_values(const VectorizedArray< Number > in[], VectorizedArray< Number > out[])
FEEvaluationGeneral(const MatrixFree< dim, Number > &matrix_free, const unsigned int fe_no=0, const unsigned int quad_no=0)
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)