Reference documentation for deal.II version 8.1.0
fe_evaluation.h
1 // ---------------------------------------------------------------------
2 // @f$Id: fe_evaluation.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2011 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 
18 #ifndef __deal2__matrix_free_fe_evaluation_h
19 #define __deal2__matrix_free_fe_evaluation_h
20 
21 
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>
28 
29 
31 
32 namespace parallel
33 {
34  namespace distributed
35  {
36  template <typename> class Vector;
37  }
38 }
39 
40 
41 
42 namespace internal
43 {
44  DeclException0 (ExcAccessToUninitializedField);
45 
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);
50 }
51 
52 
53 
85 template <int dim, int dofs_per_cell_, int n_q_points_,
86  int n_components_, typename Number>
88 {
89 public:
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_;
97 
109  void reinit (const unsigned int cell);
110 
119  unsigned int get_cell_data_number() const;
120 
127  internal::MatrixFreeFunctions::CellType get_cell_type() const;
128 
130 
151  template <typename VectorType>
152  void read_dof_values (const VectorType &src);
153 
162  template <typename VectorType>
163  void read_dof_values (const std::vector<VectorType> &src,
164  const unsigned int first_index=0);
165 
170  template <typename VectorType>
171  void read_dof_values (const std::vector<VectorType *> &src,
172  const unsigned int first_index=0);
173 
186  template <typename VectorType>
187  void read_dof_values_plain (const VectorType &src);
188 
196  template <typename VectorType>
197  void read_dof_values_plain (const std::vector<VectorType> &src,
198  const unsigned int first_index=0);
199 
205  template <typename VectorType>
206  void read_dof_values_plain (const std::vector<VectorType *> &src,
207  const unsigned int first_index=0);
208 
216  template<typename VectorType>
217  void distribute_local_to_global (VectorType &dst) const;
218 
228  template<typename VectorType>
229  void distribute_local_to_global (std::vector<VectorType> &dst,
230  const unsigned int first_index=0) const;
231 
236  template<typename VectorType>
237  void distribute_local_to_global (std::vector<VectorType *> &dst,
238  const unsigned int first_index=0) const;
239 
247  template<typename VectorType>
248  void set_dof_values (VectorType &dst) const;
249 
259  template<typename VectorType>
260  void set_dof_values (std::vector<VectorType> &dst,
261  const unsigned int first_index=0) const;
262 
267  template<typename VectorType>
268  void set_dof_values (std::vector<VectorType *> &dst,
269  const unsigned int first_index=0) const;
270 
272 
290  value_type get_dof_value (const unsigned int dof) const;
291 
302  void submit_dof_value (const value_type val_in,
303  const unsigned int dof);
304 
316  value_type get_value (const unsigned int q_point) const;
317 
329  void submit_value (const value_type val_in,
330  const unsigned int q_point);
331 
341  gradient_type get_gradient (const unsigned int q_point) const;
342 
355  void submit_gradient(const gradient_type grad_in,
356  const unsigned int q_point);
357 
369  get_hessian (const unsigned int q_point) const;
370 
379  gradient_type get_hessian_diagonal (const unsigned int q_point) const;
380 
391  value_type get_laplacian (const unsigned int q_point) const;
392 
401  value_type integrate_value () const;
402 
404 
418 
428 
439  const VectorizedArray<Number> *begin_values () const;
440 
452 
465 
478 
491  const VectorizedArray<Number> *begin_hessians () const;
492 
506 
508 
509 protected:
510 
518  FEEvaluationBase (const MatrixFree<dim,Number> &matrix_free,
519  const unsigned int fe_no = 0,
520  const unsigned int quad_no = 0);
521 
528  template<typename VectorType, typename VectorOperation>
529  void read_write_operation (const VectorOperation &operation,
530  VectorType *vectors[]) const;
531 
539  template<typename VectorType>
540  void read_dof_values_plain (const VectorType *src_data[]);
541 
554  VectorizedArray<Number> values_dofs[n_components][dofs_per_cell>0?dofs_per_cell:1];
555 
562  VectorizedArray<Number> values_quad[n_components][n_q_points>0?n_q_points:1];
563 
571  VectorizedArray<Number> gradients_quad[n_components][dim][n_q_points>0?n_q_points:1];
572 
578  VectorizedArray<Number> hessians_quad[n_components][(dim*(dim+1))/2][n_q_points>0?n_q_points:1];
579 
583  const unsigned int quad_no;
584 
589  const unsigned int n_fe_components;
590 
595  const unsigned int active_fe_index;
596 
601  const unsigned int active_quad_index;
602 
607 
613  const internal::MatrixFreeFunctions::DoFInfo &dof_info;
614 
622 
630 
636 
642 
650 
655 
660 
667 
674 
679  unsigned int cell;
680 
687  internal::MatrixFreeFunctions::CellType cell_type;
688 
692  unsigned int cell_data_number;
693 
700 
707 
714 
721 
728 
735 };
736 
737 
738 
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>
751 {
752 public:
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_;
760  typedef FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,
761  Number> BaseClass;
762 
763 protected:
771  FEEvaluationAccess (const MatrixFree<dim,Number> &matrix_free,
772  const unsigned int fe_no = 0,
773  const unsigned int quad_no = 0);
774 };
775 
776 
777 
778 
787 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
788 class FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number> :
789  public FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,1,Number>
790 {
791 public:
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_;
799 
809  value_type get_dof_value (const unsigned int dof) const;
810 
815  void submit_dof_value (const value_type val_in,
816  const unsigned int dof);
817 
825  value_type get_value (const unsigned int q_point) const;
826 
834  void submit_value (const value_type val_in,
835  const unsigned int q_point);
836 
842  gradient_type get_gradient (const unsigned int q_point) const;
843 
852  void submit_gradient(const gradient_type grad_in,
853  const unsigned int q_point);
854 
862  get_hessian (unsigned int q_point) const;
863 
868  gradient_type get_hessian_diagonal (const unsigned int q_point) const;
869 
874  value_type get_laplacian (const unsigned int q_point) const;
875 
884  value_type integrate_value () const;
885 
886 protected:
894  FEEvaluationAccess (const MatrixFree<dim,Number> &matrix_free,
895  const unsigned int fe_no = 0,
896  const unsigned int quad_no = 0);
897 };
898 
899 
900 
910 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
911 class FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number> :
912  public FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,dim,Number>
913 {
914 public:
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_;
923 
928  gradient_type get_gradient (const unsigned int q_point) const;
929 
934  VectorizedArray<Number> get_divergence (const unsigned int q_point) const;
935 
943  get_symmetric_gradient (const unsigned int q_point) const;
944 
950  get_curl (const unsigned int q_point) const;
951 
959  get_hessian (const unsigned int q_point) const;
960 
965  gradient_type get_hessian_diagonal (const unsigned int q_point) const;
966 
975  void submit_gradient(const gradient_type grad_in,
976  const unsigned int q_point);
977 
986  void submit_gradient(const Tensor<1,dim,Tensor<1,dim,VectorizedArray<Number> > > grad_in,
987  const unsigned int q_point);
988 
997  void submit_divergence (const VectorizedArray<Number> div_in,
998  const unsigned int q_point);
999 
1008  void submit_symmetric_gradient(const SymmetricTensor<2,dim,VectorizedArray<Number> > grad_in,
1009  const unsigned int q_point);
1010 
1015  void submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray<Number> > curl_in,
1016  const unsigned int q_point);
1017 
1018 protected:
1026  FEEvaluationAccess (const MatrixFree<dim,Number> &matrix_free,
1027  const unsigned int fe_no = 0,
1028  const unsigned int quad_no = 0);
1029 };
1030 
1031 
1032 
1072 template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
1073  int n_components_ = 1, typename Number = double >
1075  public FEEvaluationAccess<dim,
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>
1079 {
1080 public:
1081  typedef FEEvaluationAccess<dim,
1084  n_components_, Number> BaseClass;
1085  typedef Number number_type;
1086  typedef typename BaseClass::value_type value_type;
1087  typedef typename BaseClass::gradient_type gradient_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;
1092 
1099  FEEvaluationGeneral (const MatrixFree<dim,Number> &matrix_free,
1100  const unsigned int fe_no = 0,
1101  const unsigned int quad_no = 0);
1102 
1110  void evaluate (const bool evaluate_val,
1111  const bool evaluate_grad,
1112  const bool evaluate_hess = false);
1113 
1121  void integrate (const bool integrate_val,
1122  const bool integrate_grad);
1123 
1128  quadrature_point (const unsigned int q_point) const;
1129 
1130 protected:
1131 
1140  template <int direction, bool dof_to_quad, bool add>
1141  void apply_values (const VectorizedArray<Number> in [],
1142  VectorizedArray<Number> out []);
1143 
1152  template <int direction, bool dof_to_quad, bool add>
1153  void apply_gradients (const VectorizedArray<Number> in [],
1154  VectorizedArray<Number> out []);
1155 
1165  template <int direction, bool dof_to_quad, bool add>
1166  void apply_hessians (const VectorizedArray<Number> in [],
1167  VectorizedArray<Number> out []);
1168 
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);
1176 };
1177 
1178 
1179 
1221 template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
1222  int n_components_ = 1, typename Number = double >
1224  public FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
1225 {
1226 public:
1228  typedef Number number_type;
1229  typedef typename BaseClass::value_type value_type;
1230  typedef typename BaseClass::gradient_type gradient_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;
1235 
1242  FEEvaluation (const MatrixFree<dim,Number> &matrix_free,
1243  const unsigned int fe_no = 0,
1244  const unsigned int quad_no = 0);
1245 
1254  void evaluate (const bool evaluate_val,
1255  const bool evaluate_grad,
1256  const bool evaluate_hess = false);
1257 
1265  void integrate (const bool integrate_val,
1266  const bool integrate_grad);
1267 
1268 protected:
1269 
1278  template <int direction, bool dof_to_quad, bool add>
1279  void apply_values (const VectorizedArray<Number> in [],
1280  VectorizedArray<Number> out []);
1281 
1290  template <int direction, bool dof_to_quad, bool add>
1291  void apply_gradients (const VectorizedArray<Number> in [],
1292  VectorizedArray<Number> out []);
1293 
1303  template <int direction, bool dof_to_quad, bool add>
1304  void apply_hessians (const VectorizedArray<Number> in [],
1305  VectorizedArray<Number> out []);
1306 
1307  VectorizedArray<Number> shape_val_evenodd[fe_degree+1][(n_q_points_1d+1)/2];
1308  VectorizedArray<Number> shape_gra_evenodd[fe_degree+1][(n_q_points_1d+1)/2];
1309  VectorizedArray<Number> shape_hes_evenodd[fe_degree+1][(n_q_points_1d+1)/2];
1310 
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);
1318 };
1319 
1320 
1321 
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>
1362 {
1363 public:
1365  typedef Number number_type;
1366  typedef typename BaseClass::value_type value_type;
1367  typedef typename BaseClass::gradient_type gradient_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;
1372 
1379  FEEvaluationGL (const MatrixFree<dim,Number> &matrix_free,
1380  const unsigned int fe_no = 0,
1381  const unsigned int quad_no = 0);
1382 
1391  void evaluate (const bool evaluate_val,
1392  const bool evaluate_grad,
1393  const bool evaluate_lapl = false);
1394 
1402  void integrate (const bool integrate_val,
1403  const bool integrate_grad);
1404 
1405 protected:
1414  template <int direction, bool dof_to_quad, bool add>
1415  void apply_gradients (const VectorizedArray<Number> in [],
1416  VectorizedArray<Number> out []);
1417 };
1418 
1419 
1420 
1421 
1422 /*----------------------- Inline functions ----------------------------------*/
1423 
1424 #ifndef DOXYGEN
1425 
1426 
1427 /*----------------------- FEEvaluationBase ----------------------------------*/
1428 
1429 template <int dim, int dofs_per_cell_, int n_q_points_,
1430  int n_components_, typename Number>
1431 inline
1434  const unsigned int fe_no_in,
1435  const unsigned int quad_no_in)
1436  :
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)),
1450  cartesian_data (0),
1451  jacobian (0),
1452  J_value (0),
1453  quadrature_weights (mapping_info.mapping_data_gen[quad_no].
1454  quadrature_weights[active_quad_index].begin()),
1455  quadrature_points (0),
1456  jacobian_grad (0),
1457  jacobian_grad_upper(0),
1458  cell (numbers::invalid_unsigned_int),
1459  cell_type (internal::MatrixFreeFunctions::undefined),
1460  cell_data_number (0)
1461 {
1462  Assert (matrix_info.mapping_initialized() == true,
1463  ExcNotInitialized());
1464  AssertDimension (matrix_info.get_size_info().vectorization_length,
1466  Assert (n_fe_components == 1 ||
1467  n_components == 1 ||
1468  n_components == n_fe_components,
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."));
1472 
1473 
1474  // do not check for correct dimensions of data fields here, should be done
1475  // in derived classes
1476 }
1477 
1478 
1479 
1480 template <int dim, int dofs_per_cell_, int n_q_points_,
1481  int n_components_, typename Number>
1482 inline
1483 void
1485 ::reinit (const unsigned int cell_in)
1486 {
1487  AssertIndexRange (cell_in, dof_info.row_starts.size()-1);
1488  AssertDimension (((dof_info.cell_active_fe_index.size() > 0) ?
1489  dof_info.cell_active_fe_index[cell_in] : 0),
1490  active_fe_index);
1491  cell = cell_in;
1492  cell_type = mapping_info.get_cell_type(cell);
1493  cell_data_number = mapping_info.get_cell_data_index(cell);
1494 
1495  if (mapping_info.quadrature_points_initialized == true)
1496  {
1497  AssertIndexRange (cell_data_number, mapping_info.
1498  mapping_data_gen[quad_no].rowstart_q_points.size());
1499  const unsigned int index = mapping_info.mapping_data_gen[quad_no].
1500  rowstart_q_points[cell];
1501  AssertIndexRange (index, mapping_info.mapping_data_gen[quad_no].
1502  quadrature_points.size());
1503  quadrature_points =
1504  &mapping_info.mapping_data_gen[quad_no].quadrature_points[index];
1505  }
1506 
1507  if (cell_type == internal::MatrixFreeFunctions::cartesian)
1508  {
1509  cartesian_data = &mapping_info.cartesian_data[cell_data_number].first;
1510  J_value = &mapping_info.cartesian_data[cell_data_number].second;
1511  }
1512  else if (cell_type == internal::MatrixFreeFunctions::affine)
1513  {
1514  jacobian = &mapping_info.affine_data[cell_data_number].first;
1515  J_value = &mapping_info.affine_data[cell_data_number].second;
1516  }
1517  else
1518  {
1519  const unsigned int rowstart = mapping_info.
1520  mapping_data_gen[quad_no].rowstart_jacobians[cell_data_number];
1521  AssertIndexRange (rowstart, mapping_info.
1522  mapping_data_gen[quad_no].jacobians.size());
1523  jacobian =
1524  &mapping_info.mapping_data_gen[quad_no].jacobians[rowstart];
1525  if (mapping_info.JxW_values_initialized == true)
1526  {
1527  AssertIndexRange (rowstart, mapping_info.
1528  mapping_data_gen[quad_no].JxW_values.size());
1529  J_value = &(mapping_info.mapping_data_gen[quad_no].
1530  JxW_values[rowstart]);
1531  }
1532  if (mapping_info.second_derivatives_initialized == true)
1533  {
1534  AssertIndexRange(rowstart, mapping_info.
1535  mapping_data_gen[quad_no].jacobians_grad_diag.size());
1536  jacobian_grad = &mapping_info.mapping_data_gen[quad_no].
1537  jacobians_grad_diag[rowstart];
1538  AssertIndexRange(rowstart, mapping_info.
1539  mapping_data_gen[quad_no].jacobians_grad_upper.size());
1540  jacobian_grad_upper = &mapping_info.mapping_data_gen[quad_no].
1541  jacobians_grad_upper[rowstart];
1542  }
1543  }
1544 
1545 #ifdef DEBUG
1546  dof_values_initialized = false;
1547  values_quad_initialized = false;
1548  gradients_quad_initialized = false;
1549  hessians_quad_initialized = false;
1550 #endif
1551 }
1552 
1553 
1554 
1555 template <int dim, int dofs_per_cell_, int n_q_points_,
1556  int n_components_, typename Number>
1557 inline
1558 unsigned int
1560 ::get_cell_data_number () const
1561 {
1563  return cell_data_number;
1564 }
1565 
1566 
1567 
1568 template <int dim, int dofs_per_cell_, int n_q_points_,
1569  int n_components_, typename Number>
1570 inline
1571 internal::MatrixFreeFunctions::CellType
1573 ::get_cell_type () const
1574 {
1576  return cell_type;
1577 }
1578 
1579 
1580 
1581 namespace internal
1582 {
1583  // write access to generic vectors that have operator ().
1584  template <typename VectorType>
1585  inline
1586  typename VectorType::value_type &
1587  vector_access (VectorType &vec,
1588  const unsigned int entry)
1589  {
1590  return vec(entry);
1591  }
1592 
1593 
1594 
1595  // read access to generic vectors that have operator ().
1596  template <typename VectorType>
1597  inline
1598  typename VectorType::value_type
1599  vector_access (const VectorType &vec,
1600  const unsigned int entry)
1601  {
1602  return vec(entry);
1603  }
1604 
1605 
1606 
1607  // write access to distributed MPI vectors that have a local_element(uint)
1608  // method to access data in local index space, which is what we use in
1609  // DoFInfo and hence in read_dof_values etc.
1610  template <typename Number>
1611  inline
1612  Number &
1613  vector_access (parallel::distributed::Vector<Number> &vec,
1614  const unsigned int entry)
1615  {
1616  return vec.local_element(entry);
1617  }
1618 
1619 
1620 
1621  // read access to distributed MPI vectors that have a local_element(uint)
1622  // method to access data in local index space, which is what we use in
1623  // DoFInfo and hence in read_dof_values etc.
1624  template <typename Number>
1625  inline
1626  Number
1627  vector_access (const parallel::distributed::Vector<Number> &vec,
1628  const unsigned int entry)
1629  {
1630  return vec.local_element(entry);
1631  }
1632 
1633 
1634 
1635  // this is to make sure that the parallel partitioning in the
1636  // parallel::distributed::Vector is really the same as stored in MatrixFree
1637  template <typename VectorType>
1638  inline
1639  void check_vector_compatibility (const VectorType &vec,
1640  const internal::MatrixFreeFunctions::DoFInfo &dof_info)
1641  {
1642  AssertDimension (vec.size(),
1643  dof_info.vector_partitioner->size());
1644  }
1645 
1646  template <typename Number>
1647  inline
1648  void check_vector_compatibility (const parallel::distributed::Vector<Number> &vec,
1649  const internal::MatrixFreeFunctions::DoFInfo &dof_info)
1650  {
1651  Assert (vec.partitioners_are_compatible(*dof_info.vector_partitioner),
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."));
1656  }
1657 
1658  // A class to use the same code to read from and write to vector
1659  template <typename Number>
1660  struct VectorReader
1661  {
1662  template <typename VectorType>
1663  void process_dof (const unsigned int index,
1664  VectorType &vec,
1665  Number &res) const
1666  {
1667  res = vector_access (const_cast<const VectorType &>(vec), index);
1668  }
1669 
1670  void pre_constraints (const Number &,
1671  Number &res) const
1672  {
1673  res = Number();
1674  }
1675 
1676  template <typename VectorType>
1677  void process_constraint (const unsigned int index,
1678  const Number weight,
1679  VectorType &vec,
1680  Number &res) const
1681  {
1682  res += weight * vector_access (const_cast<const VectorType &>(vec), index);
1683  }
1684 
1685  void post_constraints (const Number &sum,
1686  Number &write_pos) const
1687  {
1688  write_pos = sum;
1689  }
1690 
1691  void process_empty (Number &res) const
1692  {
1693  res = Number();
1694  }
1695  };
1696 
1697  // A class to use the same code to read from and write to vector
1698  template <typename Number>
1699  struct VectorDistributorLocalToGlobal
1700  {
1701  template <typename VectorType>
1702  void process_dof (const unsigned int index,
1703  VectorType &vec,
1704  Number &res) const
1705  {
1706  vector_access (vec, index) += res;
1707  }
1708 
1709  void pre_constraints (const Number &input,
1710  Number &res) const
1711  {
1712  res = input;
1713  }
1714 
1715  template <typename VectorType>
1716  void process_constraint (const unsigned int index,
1717  const Number weight,
1718  VectorType &vec,
1719  Number &res) const
1720  {
1721  vector_access (vec, index) += weight * res;
1722  }
1723 
1724  void post_constraints (const Number &,
1725  Number &) const
1726  {
1727  }
1728 
1729  void process_empty (Number &) const
1730  {
1731  }
1732  };
1733 
1734 
1735  // A class to use the same code to read from and write to vector
1736  template <typename Number>
1737  struct VectorSetter
1738  {
1739  template <typename VectorType>
1740  void process_dof (const unsigned int index,
1741  VectorType &vec,
1742  Number &res) const
1743  {
1744  vector_access (vec, index) = res;
1745  }
1746 
1747  void pre_constraints (const Number &,
1748  Number &) const
1749  {
1750  }
1751 
1752  template <typename VectorType>
1753  void process_constraint (const unsigned int,
1754  const Number,
1755  VectorType &,
1756  Number &) const
1757  {
1758  }
1759 
1760  void post_constraints (const Number &,
1761  Number &) const
1762  {
1763  }
1764 
1765  void process_empty (Number &) const
1766  {
1767  }
1768  };
1769 
1770  // allows to select between block vectors and non-block vectors, which
1771  // allows to use a unified interface for extracting blocks on block vectors
1772  // and doing nothing on usual vectors
1773  template <typename VectorType, bool>
1774  struct BlockVectorSelector {};
1775 
1776  template <typename VectorType>
1777  struct BlockVectorSelector<VectorType,true>
1778  {
1779  typedef typename VectorType::BlockType BaseVectorType;
1780 
1781  static BaseVectorType *get_vector_component (VectorType &vec,
1782  const unsigned int component)
1783  {
1784  AssertIndexRange (component, vec.n_blocks());
1785  return &vec.block(component);
1786  }
1787  };
1788 
1789  template <typename VectorType>
1790  struct BlockVectorSelector<VectorType,false>
1791  {
1792  typedef VectorType BaseVectorType;
1793 
1794  static BaseVectorType *get_vector_component (VectorType &vec,
1795  const unsigned int)
1796  {
1797  return &vec;
1798  }
1799  };
1800 }
1801 
1802 
1803 
1804 
1805 template <int dim, int dofs_per_cell_, int n_q_points_,
1806  int n_components_, typename Number>
1807 template<typename VectorType, typename VectorOperation>
1808 inline
1809 void
1811 ::read_write_operation (const VectorOperation &operation,
1812  VectorType *src[]) const
1813 {
1814  // This functions processes all the functions read_dof_values,
1815  // distribute_local_to_global, and set_dof_values with the same code. The
1816  // distinction between these three cases is made by the input
1817  // VectorOperation that either reads values from a vector and puts the data
1818  // into the local data field or write local data into the vector. Certain
1819  // operations are no-ops for the given use case.
1820 
1821  Assert (matrix_info.indices_initialized() == true,
1822  ExcNotInitialized());
1824 
1825  // loop over all local dofs. ind_local holds local number on cell, index
1826  // iterates over the elements of index_local_to_global and dof_indices
1827  // points to the global indices stored in index_local_to_global
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;
1834 
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;
1837 
1838  // scalar case (or case when all components have the same degrees of freedom
1839  // and sit on a different vector each)
1840  if (n_fe_components == 1)
1841  {
1842  const unsigned int n_local_dofs =
1844  for (unsigned int comp=0; comp<n_components; ++comp)
1845  internal::check_vector_compatibility (*src[comp], dof_info);
1846  Number *local_data [n_components];
1847  for (unsigned int comp=0; comp<n_components; ++comp)
1848  local_data[comp] =
1849  const_cast<Number *>(&values_dofs[comp][0][0]);
1850 
1851  // standard case where there are sufficiently many cells to fill all
1852  // vectors
1853  if (at_irregular_cell == false)
1854  {
1855  // check whether there is any constraint on the current cell
1856  if (indicators != indicators_end)
1857  {
1858  for ( ; indicators != indicators_end; ++indicators)
1859  {
1860  // run through values up to next constraint
1861  for (unsigned int j=0; j<indicators->first; ++j)
1862  for (unsigned int comp=0; comp<n_components; ++comp)
1863  operation.process_dof (dof_indices[j], *src[comp],
1864  local_data[comp][ind_local+j]);
1865 
1866  ind_local += indicators->first;
1867  dof_indices += indicators->first;
1868 
1869  // constrained case: build the local value as a linear
1870  // combination of the global value according to constraints
1871  Number value [n_components];
1872  for (unsigned int comp=0; comp<n_components; ++comp)
1873  operation.pre_constraints (local_data[comp][ind_local],
1874  value[comp]);
1875 
1876  const Number *data_val =
1877  matrix_info.constraint_pool_begin(indicators->second);
1878  const Number *end_pool =
1879  matrix_info.constraint_pool_end(indicators->second);
1880  for ( ; data_val != end_pool; ++data_val, ++dof_indices)
1881  for (unsigned int comp=0; comp<n_components; ++comp)
1882  operation.process_constraint (*dof_indices, *data_val,
1883  *src[comp], value[comp]);
1884 
1885  for (unsigned int comp=0; comp<n_components; ++comp)
1886  operation.post_constraints (value[comp],
1887  local_data[comp][ind_local]);
1888 
1889  ind_local++;
1890  }
1891 
1892  // get the dof values past the last constraint
1893  for (; ind_local < n_local_dofs; ++dof_indices, ++ind_local)
1894  {
1895  for (unsigned int comp=0; comp<n_components; ++comp)
1896  operation.process_dof (*dof_indices, *src[comp],
1897  local_data[comp][ind_local]);
1898  }
1899  }
1900  else
1901  {
1902  // no constraint at all: loop bounds are known, compiler can
1903  // unroll without checks
1904  AssertDimension (dof_info.end_indices(cell)-dof_indices,
1905  static_cast<int>(n_local_dofs));
1906  for (unsigned int j=0; j<n_local_dofs; ++j)
1907  for (unsigned int comp=0; comp<n_components; ++comp)
1908  operation.process_dof (dof_indices[j], *src[comp],
1909  local_data[comp][j]);
1910  }
1911  }
1912 
1913  // non-standard case: need to fill in zeros for those components that
1914  // are not present (a bit more expensive), but there is not more than
1915  // one such cell
1916  else
1917  {
1918  Assert (n_irreg_components_filled > 0, ExcInternalError());
1919  for ( ; indicators != indicators_end; ++indicators)
1920  {
1921  for (unsigned int j=0; j<indicators->first; ++j)
1922  {
1923  // non-constrained case: copy the data from the global
1924  // vector, src, to the local one, local_src.
1925  for (unsigned int comp=0; comp<n_components; ++comp)
1926  operation.process_dof (dof_indices[j], *src[comp],
1927  local_data[comp][ind_local]);
1928 
1929  // here we jump over all the components that are artificial
1930  ++ind_local;
1932  >= n_irreg_components_filled)
1933  {
1934  for (unsigned int comp=0; comp<n_components; ++comp)
1935  operation.process_empty (local_data[comp][ind_local]);
1936  ++ind_local;
1937  }
1938  }
1939  dof_indices += indicators->first;
1940 
1941  // constrained case: build the local value as a linear
1942  // combination of the global value according to constraint
1943  Number value [n_components];
1944  for (unsigned int comp=0; comp<n_components; ++comp)
1945  operation.pre_constraints (local_data[comp][ind_local],
1946  value[comp]);
1947 
1948  const Number *data_val =
1949  matrix_info.constraint_pool_begin(indicators->second);
1950  const Number *end_pool =
1951  matrix_info.constraint_pool_end(indicators->second);
1952 
1953  for ( ; data_val != end_pool; ++data_val, ++dof_indices)
1954  for (unsigned int comp=0; comp<n_components; ++comp)
1955  operation.process_constraint (*dof_indices, *data_val,
1956  *src[comp], value[comp]);
1957 
1958  for (unsigned int comp=0; comp<n_components; ++comp)
1959  operation.post_constraints (value[comp],
1960  local_data[comp][ind_local]);
1961  ind_local++;
1963  >= n_irreg_components_filled)
1964  {
1965  for (unsigned int comp=0; comp<n_components; ++comp)
1966  operation.process_empty (local_data[comp][ind_local]);
1967  ++ind_local;
1968  }
1969  }
1970  for (; ind_local<n_local_dofs; ++dof_indices)
1971  {
1972  Assert (dof_indices != dof_info.end_indices(cell),
1973  ExcInternalError());
1974 
1975  // non-constrained case: copy the data from the global vector,
1976  // src, to the local one, local_dst.
1977  for (unsigned int comp=0; comp<n_components; ++comp)
1978  operation.process_dof (*dof_indices, *src[comp],
1979  local_data[comp][ind_local]);
1980  ++ind_local;
1982  >= n_irreg_components_filled)
1983  {
1984  for (unsigned int comp=0; comp<n_components; ++comp)
1985  operation.process_empty(local_data[comp][ind_local]);
1986  ++ind_local;
1987  }
1988  }
1989  }
1990  }
1991  else
1992  // case with vector-valued finite elements where all components are
1993  // included in one single vector. Assumption: first come all entries to
1994  // the first component, then all entries to the second one, and so
1995  // on. This is ensured by the way MatrixFree reads out the indices.
1996  {
1997  internal::check_vector_compatibility (*src[0], dof_info);
1998  Assert (n_fe_components == n_components_, ExcNotImplemented());
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)
2004  {
2005  // check whether there is any constraint on the current cell
2006  if (indicators != indicators_end)
2007  {
2008  for ( ; indicators != indicators_end; ++indicators)
2009  {
2010  // run through values up to next constraint
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;
2016 
2017  // constrained case: build the local value as a linear
2018  // combination of the global value according to constraints
2019  Number value;
2020  operation.pre_constraints (local_data[ind_local], value);
2021 
2022  const Number *data_val =
2023  matrix_info.constraint_pool_begin(indicators->second);
2024  const Number *end_pool =
2025  matrix_info.constraint_pool_end(indicators->second);
2026 
2027  for ( ; data_val != end_pool; ++data_val, ++dof_indices)
2028  operation.process_constraint (*dof_indices, *data_val,
2029  *src[0], value);
2030 
2031  operation.post_constraints (value, local_data[ind_local]);
2032  ind_local++;
2033  }
2034 
2035  // get the dof values past the last constraint
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),
2040  ExcInternalError());
2041  }
2042  else
2043  {
2044  // no constraint at all: loop bounds are known, compiler can
2045  // unroll without checks
2046  AssertDimension (dof_info.end_indices(cell)-dof_indices,
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],
2050  local_data[j]);
2051  }
2052  }
2053 
2054  // non-standard case: need to fill in zeros for those components that
2055  // are not present (a bit more expensive), but there is not more than
2056  // one such cell
2057  else
2058  {
2059  Assert (n_irreg_components_filled > 0, ExcInternalError());
2060  for ( ; indicators != indicators_end; ++indicators)
2061  {
2062  for (unsigned int j=0; j<indicators->first; ++j)
2063  {
2064  // non-constrained case: copy the data from the global
2065  // vector, src, to the local one, local_src.
2066  operation.process_dof (dof_indices[j], *src[0],
2067  local_data[ind_local]);
2068 
2069  // here we jump over all the components that are artificial
2070  ++ind_local;
2072  >= n_irreg_components_filled)
2073  {
2074  operation.process_empty (local_data[ind_local]);
2075  ++ind_local;
2076  }
2077  }
2078  dof_indices += indicators->first;
2079 
2080  // constrained case: build the local value as a linear
2081  // combination of the global value according to constraint
2082  Number value;
2083  operation.pre_constraints (local_data[ind_local], value);
2084 
2085  const Number *data_val =
2086  matrix_info.constraint_pool_begin(indicators->second);
2087  const Number *end_pool =
2088  matrix_info.constraint_pool_end(indicators->second);
2089 
2090  for ( ; data_val != end_pool; ++data_val, ++dof_indices)
2091  operation.process_constraint (*dof_indices, *data_val,
2092  *src[0], value);
2093 
2094  operation.post_constraints (value, local_data[ind_local]);
2095  ind_local++;
2097  >= n_irreg_components_filled)
2098  {
2099  operation.process_empty (local_data[ind_local]);
2100  ++ind_local;
2101  }
2102  }
2103  for (; ind_local<n_local_dofs; ++dof_indices)
2104  {
2105  Assert (dof_indices != dof_info.end_indices(cell),
2106  ExcInternalError());
2107 
2108  // non-constrained case: copy the data from the global vector,
2109  // src, to the local one, local_dst.
2110  operation.process_dof (*dof_indices, *src[0],
2111  local_data[ind_local]);
2112  ++ind_local;
2114  >= n_irreg_components_filled)
2115  {
2116  operation.process_empty (local_data[ind_local]);
2117  ++ind_local;
2118  }
2119  }
2120  }
2121  }
2122 }
2123 
2124 
2125 
2126 template <int dim, int dofs_per_cell_, int n_q_points_,
2127  int n_components_, typename Number>
2128 template<typename VectorType>
2129 inline
2130 void
2132 ::read_dof_values (const VectorType &src)
2133 {
2134  // select between block vectors and non-block vectors. Note that the number
2135  // of components is checked in the internal data
2136  typename internal::BlockVectorSelector<VectorType,
2137  IsBlockVector<VectorType>::value>::BaseVectorType *src_data[n_components];
2138  for (unsigned int d=0; d<n_components; ++d)
2139  src_data[d] = internal::BlockVectorSelector<VectorType, IsBlockVector<VectorType>::value>::get_vector_component(const_cast<VectorType &>(src), d);
2140 
2141  internal::VectorReader<Number> reader;
2142  read_write_operation (reader, src_data);
2143 
2144 #ifdef DEBUG
2145  dof_values_initialized = true;
2146 #endif
2147 }
2148 
2149 
2150 
2151 template <int dim, int dofs_per_cell_, int n_q_points_,
2152  int n_components_, typename Number>
2153 template<typename VectorType>
2154 inline
2155 void
2157 ::read_dof_values (const std::vector<VectorType> &src,
2158  const unsigned int first_index)
2159 {
2160  AssertIndexRange (first_index, src.size());
2161  Assert (n_fe_components == 1, ExcNotImplemented());
2162  Assert ((n_fe_components == 1 ?
2163  (first_index+n_components <= src.size()) : true),
2164  ExcIndexRange (first_index + n_components_, 0, src.size()));
2165 
2166  VectorType *src_data [n_components];
2167  for (unsigned int comp=0; comp<n_components; ++comp)
2168  src_data[comp] = const_cast<VectorType *>(&src[comp+first_index]);
2169 
2170  internal::VectorReader<Number> reader;
2171  read_write_operation (reader, src_data);
2172 
2173 #ifdef DEBUG
2174  dof_values_initialized = true;
2175 #endif
2176 }
2177 
2178 
2179 
2180 template <int dim, int dofs_per_cell_, int n_q_points_,
2181  int n_components_, typename Number>
2182 template<typename VectorType>
2183 inline
2184 void
2186 ::read_dof_values (const std::vector<VectorType *> &src,
2187  const unsigned int first_index)
2188 {
2189  AssertIndexRange (first_index, src.size());
2190  Assert (n_fe_components == 1, ExcNotImplemented());
2191  Assert ((n_fe_components == 1 ?
2192  (first_index+n_components <= src.size()) : true),
2193  ExcIndexRange (first_index + n_components_, 0, src.size()));
2194 
2195  const VectorType *src_data [n_components];
2196  for (unsigned int comp=0; comp<n_components; ++comp)
2197  src_data[comp] = const_cast<VectorType *>(src[comp+first_index]);
2198 
2199  internal::VectorReader<Number> reader;
2200  read_write_operation (reader, src_data);
2201 
2202 #ifdef DEBUG
2203  dof_values_initialized = true;
2204 #endif
2205 }
2206 
2207 
2208 
2209 template <int dim, int dofs_per_cell_, int n_q_points_,
2210  int n_components_, typename Number>
2211 template<typename VectorType>
2212 inline
2213 void
2215 ::read_dof_values_plain (const VectorType &src)
2216 {
2217  // select between block vectors and non-block vectors. Note that the number
2218  // of components is checked in the internal data
2219  const typename internal::BlockVectorSelector<VectorType,
2220  IsBlockVector<VectorType>::value>::BaseVectorType *src_data[n_components];
2221  for (unsigned int d=0; d<n_components; ++d)
2222  src_data[d] = internal::BlockVectorSelector<VectorType, IsBlockVector<VectorType>::value>::get_vector_component(const_cast<VectorType &>(src), d);
2223 
2224  read_dof_values_plain (src_data);
2225 }
2226 
2227 
2228 
2229 template <int dim, int dofs_per_cell_, int n_q_points_,
2230  int n_components_, typename Number>
2231 template<typename VectorType>
2232 inline
2233 void
2235 ::read_dof_values_plain (const std::vector<VectorType> &src,
2236  const unsigned int first_index)
2237 {
2238  AssertIndexRange (first_index, src.size());
2239  Assert (n_fe_components == 1, ExcNotImplemented());
2240  Assert ((n_fe_components == 1 ?
2241  (first_index+n_components <= src.size()) : true),
2242  ExcIndexRange (first_index + n_components_, 0, src.size()));
2243  const VectorType *src_data [n_components];
2244  for (unsigned int comp=0; comp<n_components; ++comp)
2245  src_data[comp] = &src[comp+first_index];
2246  read_dof_values_plain (src_data);
2247 }
2248 
2249 
2250 
2251 template <int dim, int dofs_per_cell_, int n_q_points_,
2252  int n_components_, typename Number>
2253 template<typename VectorType>
2254 inline
2255 void
2257 ::read_dof_values_plain (const std::vector<VectorType *> &src,
2258  const unsigned int first_index)
2259 {
2260  AssertIndexRange (first_index, src.size());
2261  Assert (n_fe_components == 1, ExcNotImplemented());
2262  Assert ((n_fe_components == 1 ?
2263  (first_index+n_components <= src.size()) : true),
2264  ExcIndexRange (first_index + n_components_, 0, src.size()));
2265  const VectorType *src_data [n_components];
2266  for (unsigned int comp=0; comp<n_components; ++comp)
2267  src_data[comp] = src[comp+first_index];
2268  read_dof_values_plain (src_data);
2269 }
2270 
2271 
2272 
2273 template <int dim, int dofs_per_cell_, int n_q_points_,
2274  int n_components_, typename Number>
2275 template<typename VectorType>
2276 inline
2277 void
2279 ::distribute_local_to_global (VectorType &dst) const
2280 {
2281  Assert (dof_values_initialized==true,
2282  internal::ExcAccessToUninitializedField());
2283 
2284  // select between block vectors and non-block vectors. Note that the number
2285  // of components is checked in the internal data
2286  typename internal::BlockVectorSelector<VectorType,
2287  IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
2288  for (unsigned int d=0; d<n_components; ++d)
2289  dst_data[d] = internal::BlockVectorSelector<VectorType, IsBlockVector<VectorType>::value>::get_vector_component(dst, d);
2290 
2291  internal::VectorDistributorLocalToGlobal<Number> distributor;
2292  read_write_operation (distributor, dst_data);
2293 }
2294 
2295 
2296 
2297 template <int dim, int dofs_per_cell_, int n_q_points_,
2298  int n_components_, typename Number>
2299 template<typename VectorType>
2300 inline
2301 void
2303 ::distribute_local_to_global (std::vector<VectorType> &dst,
2304  const unsigned int first_index) const
2305 {
2306  AssertIndexRange (first_index, dst.size());
2307  Assert (n_fe_components == 1, ExcNotImplemented());
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());
2313 
2314  VectorType *dst_data [n_components];
2315  for (unsigned int comp=0; comp<n_components; ++comp)
2316  dst_data[comp] = &dst[comp+first_index];
2317 
2318  internal::VectorDistributorLocalToGlobal<Number> distributor;
2319  read_write_operation (distributor, dst_data);
2320 }
2321 
2322 
2323 
2324 template <int dim, int dofs_per_cell_, int n_q_points_,
2325  int n_components_, typename Number>
2326 template<typename VectorType>
2327 inline
2328 void
2330 ::distribute_local_to_global (std::vector<VectorType *> &dst,
2331  const unsigned int first_index) const
2332 {
2333  AssertIndexRange (first_index, dst.size());
2334  Assert (n_fe_components == 1, ExcNotImplemented());
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());
2340 
2341  VectorType *dst_data [n_components];
2342  for (unsigned int comp=0; comp<n_components; ++comp)
2343  dst_data[comp] = dst[comp+first_index];
2344 
2345  internal::VectorDistributorLocalToGlobal<Number> distributor;
2346  read_write_operation (distributor, dst_data);
2347 }
2348 
2349 
2350 
2351 template <int dim, int dofs_per_cell_, int n_q_points_,
2352  int n_components_, typename Number>
2353 template<typename VectorType>
2354 inline
2355 void
2357 ::set_dof_values (VectorType &dst) const
2358 {
2359  Assert (dof_values_initialized==true,
2360  internal::ExcAccessToUninitializedField());
2361 
2362  // select between block vectors and non-block vectors. Note that the number
2363  // of components is checked in the internal data
2364  typename internal::BlockVectorSelector<VectorType,
2365  IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
2366  for (unsigned int d=0; d<n_components; ++d)
2367  dst_data[d] = internal::BlockVectorSelector<VectorType, IsBlockVector<VectorType>::value>::get_vector_component(dst, d);
2368 
2369  internal::VectorSetter<Number> setter;
2370  read_write_operation (setter, dst_data);
2371 }
2372 
2373 
2374 
2375 template <int dim, int dofs_per_cell_, int n_q_points_,
2376  int n_components_, typename Number>
2377 template<typename VectorType>
2378 inline
2379 void
2381 ::set_dof_values (std::vector<VectorType> &dst,
2382  const unsigned int first_index) const
2383 {
2384  AssertIndexRange (first_index, dst.size());
2385  Assert (n_fe_components == 1, ExcNotImplemented());
2386  Assert ((n_fe_components == 1 ?
2387  (first_index+n_components <= dst.size()) : true),
2388  ExcIndexRange (first_index + n_components_, 0, dst.size()));
2389 
2390  Assert (dof_values_initialized==true,
2391  internal::ExcAccessToUninitializedField());
2392 
2393  VectorType *dst_data [n_components];
2394  for (unsigned int comp=0; comp<n_components; ++comp)
2395  dst_data[comp] = &dst[comp+first_index];
2396 
2397  internal::VectorSetter<Number> setter;
2398  read_write_operation (setter, dst_data);
2399 }
2400 
2401 
2402 
2403 template <int dim, int dofs_per_cell_, int n_q_points_,
2404  int n_components_, typename Number>
2405 template<typename VectorType>
2406 inline
2407 void
2409 ::set_dof_values (std::vector<VectorType *> &dst,
2410  const unsigned int first_index) const
2411 {
2412  AssertIndexRange (first_index, dst.size());
2413  Assert (n_fe_components == 1, ExcNotImplemented());
2414  Assert ((n_fe_components == 1 ?
2415  (first_index+n_components <= dst.size()) : true),
2416  ExcIndexRange (first_index + n_components_, 0, dst.size()));
2417 
2418  Assert (dof_values_initialized==true,
2419  internal::ExcAccessToUninitializedField());
2420 
2421  VectorType *dst_data [n_components];
2422  for (unsigned int comp=0; comp<n_components; ++comp)
2423  dst_data[comp] = dst[comp+first_index];
2424 
2425  internal::VectorSetter<Number> setter;
2426  read_write_operation (setter, dst_data);
2427 }
2428 
2429 
2430 
2431 template <int dim, int dofs_per_cell_, int n_q_points_,
2432  int n_components_, typename Number>
2433 template<typename VectorType>
2434 inline
2435 void
2437 ::read_dof_values_plain (const VectorType *src[])
2438 {
2439  // this is different from the other three operations because we do not use
2440  // constraints here, so this is a separate function.
2441  Assert (matrix_info.indices_initialized() == true,
2442  ExcNotInitialized());
2444  Assert (dof_info.store_plain_indices == true, ExcNotInitialized());
2445 
2446  // loop over all local dofs. ind_local holds local number on cell, index
2447  // iterates over the elements of index_local_to_global and dof_indices
2448  // points to the global indices stored in index_local_to_global
2449  const unsigned int *dof_indices = dof_info.begin_indices_plain(cell);
2450 
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;
2453 
2454  // scalar case (or case when all components have the same degrees of freedom
2455  // and sit on a different vector each)
2456  if (n_fe_components == 1)
2457  {
2458  const unsigned int n_local_dofs =
2460  for (unsigned int comp=0; comp<n_components; ++comp)
2461  internal::check_vector_compatibility (*src[comp], dof_info);
2462  Number *local_src_number [n_components];
2463  for (unsigned int comp=0; comp<n_components; ++comp)
2464  local_src_number[comp] = &values_dofs[comp][0][0];
2465 
2466  // standard case where there are sufficiently many cells to fill all
2467  // vectors
2468  if (at_irregular_cell == false)
2469  {
2470  for (unsigned int j=0; j<n_local_dofs; ++j)
2471  for (unsigned int comp=0; comp<n_components; ++comp)
2472  local_src_number[comp][j] =
2473  internal::vector_access (*src[comp], dof_indices[j]);
2474  }
2475 
2476  // non-standard case: need to fill in zeros for those components that
2477  // are not present (a bit more expensive), but there is not more than
2478  // one such cell
2479  else
2480  {
2481  Assert (n_irreg_components_filled > 0, ExcInternalError());
2482  for (unsigned int ind_local=0; ind_local<n_local_dofs;
2483  ++dof_indices)
2484  {
2485  // non-constrained case: copy the data from the global vector,
2486  // src, to the local one, local_dst.
2487  for (unsigned int comp=0; comp<n_components; ++comp)
2488  local_src_number[comp][ind_local] =
2489  internal::vector_access (*src[comp], *dof_indices);
2490  ++ind_local;
2491  while (ind_local % VectorizedArray<Number>::n_array_elements >= n_irreg_components_filled)
2492  {
2493  for (unsigned int comp=0; comp<n_components; ++comp)
2494  local_src_number[comp][ind_local] = 0.;
2495  ++ind_local;
2496  }
2497  }
2498  }
2499  }
2500  else
2501  // case with vector-valued finite elements where all components are
2502  // included in one single vector. Assumption: first come all entries to
2503  // the first component, then all entries to the second one, and so
2504  // on. This is ensured by the way MatrixFree reads out the indices.
2505  {
2506  internal::check_vector_compatibility (*src[0], dof_info);
2507  Assert (n_fe_components == n_components_, ExcNotImplemented());
2508  const unsigned int n_local_dofs =
2510  Number *local_src_number = &values_dofs[0][0][0];
2511  if (at_irregular_cell == false)
2512  {
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]);
2516  }
2517 
2518  // non-standard case: need to fill in zeros for those components that
2519  // are not present (a bit more expensive), but there is not more than
2520  // one such cell
2521  else
2522  {
2523  Assert (n_irreg_components_filled > 0, ExcInternalError());
2524  for (unsigned int ind_local=0; ind_local<n_local_dofs; ++dof_indices)
2525  {
2526  // non-constrained case: copy the data from the global vector,
2527  // src, to the local one, local_dst.
2528  local_src_number[ind_local] =
2529  internal::vector_access (*src[0], *dof_indices);
2530  ++ind_local;
2531  while (ind_local % VectorizedArray<Number>::n_array_elements >= n_irreg_components_filled)
2532  {
2533  local_src_number[ind_local] = 0.;
2534  ++ind_local;
2535  }
2536  }
2537  }
2538  }
2539 
2540 #ifdef DEBUG
2541  dof_values_initialized = true;
2542 #endif
2543 }
2544 
2545 
2546 
2547 
2548 /*------------------------------ access to data fields ----------------------*/
2549 
2550 
2551 template <int dim, int dofs_per_cell_, int n_q_points_,
2552  int n_components, typename Number>
2553 inline
2556 begin_dof_values () const
2557 {
2558  return &values_dofs[0][0];
2559 }
2560 
2561 
2562 
2563 template <int dim, int dofs_per_cell_, int n_q_points_,
2564  int n_components, typename Number>
2565 inline
2569 {
2570 #ifdef DEBUG
2571  dof_values_initialized = true;
2572 #endif
2573  return &values_dofs[0][0];
2574 }
2575 
2576 
2577 
2578 template <int dim, int dofs_per_cell_, int n_q_points_,
2579  int n_components, typename Number>
2580 inline
2583 begin_values () const
2584 {
2585  Assert (values_quad_initialized || values_quad_submitted,
2586  ExcNotInitialized());
2587  return &values_quad[0][0];
2588 }
2589 
2590 
2591 
2592 template <int dim, int dofs_per_cell_, int n_q_points_,
2593  int n_components, typename Number>
2594 inline
2597 begin_values ()
2598 {
2599 #ifdef DEBUG
2600  values_quad_submitted = true;
2601 #endif
2602  return &values_quad[0][0];
2603 }
2604 
2605 
2606 
2607 template <int dim, int dofs_per_cell_, int n_q_points_,
2608  int n_components, typename Number>
2609 inline
2612 begin_gradients () const
2613 {
2614  Assert (gradients_quad_initialized || gradients_quad_submitted,
2615  ExcNotInitialized());
2616  return &gradients_quad[0][0][0];
2617 }
2618 
2619 
2620 
2621 template <int dim, int dofs_per_cell_, int n_q_points_,
2622  int n_components, typename Number>
2623 inline
2626 begin_gradients ()
2627 {
2628 #ifdef DEBUG
2629  gradients_quad_submitted = true;
2630 #endif
2631  return &gradients_quad[0][0][0];
2632 }
2633 
2634 
2635 
2636 template <int dim, int dofs_per_cell_, int n_q_points_,
2637  int n_components, typename Number>
2638 inline
2641 begin_hessians () const
2642 {
2643  Assert (hessians_quad_initialized, ExcNotInitialized());
2644  return &hessians_quad[0][0][0];
2645 }
2646 
2647 
2648 
2649 template <int dim, int dofs_per_cell_, int n_q_points_,
2650  int n_components, typename Number>
2651 inline
2654 begin_hessians ()
2655 {
2656  return &hessians_quad[0][0][0];
2657 }
2658 
2659 
2660 
2661 template <int dim, int dofs_per_cell_, int n_q_points_,
2662  int n_components_, typename Number>
2663 inline
2666 ::get_dof_value (const unsigned int dof) const
2667 {
2668  AssertIndexRange (dof, dofs_per_cell);
2669  Tensor<1,n_components_,VectorizedArray<Number> > return_value (false);
2670  for (unsigned int comp=0; comp<n_components; comp++)
2671  return_value[comp] = this->values_dofs[comp][dof];
2672  return return_value;
2673 }
2674 
2675 
2676 
2677 template <int dim, int dofs_per_cell_, int n_q_points_,
2678  int n_components_, typename Number>
2679 inline
2682 ::get_value (const unsigned int q_point) const
2683 {
2684  Assert (this->values_quad_initialized==true,
2685  internal::ExcAccessToUninitializedField());
2686  AssertIndexRange (q_point, n_q_points);
2687  Tensor<1,n_components_,VectorizedArray<Number> > return_value (false);
2688  for (unsigned int comp=0; comp<n_components; comp++)
2689  return_value[comp] = this->values_quad[comp][q_point];
2690  return return_value;
2691 }
2692 
2693 
2694 
2695 template <int dim, int dofs_per_cell_, int n_q_points_,
2696  int n_components_, typename Number>
2697 inline
2700 ::get_gradient (const unsigned int q_point) const
2701 {
2702  Assert (this->gradients_quad_initialized==true,
2703  internal::ExcAccessToUninitializedField());
2704  AssertIndexRange (q_point, n_q_points);
2705 
2707 
2708  // Cartesian cell
2709  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
2710  {
2711  for (unsigned int comp=0; comp<n_components; comp++)
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]);
2715  }
2716  // cell with general/affine Jacobian
2717  else
2718  {
2720  this->cell_type == internal::MatrixFreeFunctions::general ?
2721  jacobian[q_point] : jacobian[0];
2722  for (unsigned int comp=0; comp<n_components; comp++)
2723  {
2724  for (unsigned int d=0; d<dim; ++d)
2725  {
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]);
2731  }
2732  }
2733  }
2734  return grad_out;
2735 }
2736 
2737 
2738 
2739 namespace internal
2740 {
2741  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
2742  // store the lower diagonal because of symmetry
2743  template <int dim, int n_q_points, typename Number>
2744  inline
2745  void
2746  hessian_unit_times_jac (const Tensor<2,dim,VectorizedArray<Number> > &jac,
2747  const VectorizedArray<Number> hessians_quad[][n_q_points],
2748  const unsigned int q_point,
2749  VectorizedArray<Number> tmp[dim][dim])
2750  {
2751  for (unsigned int d=0; d<dim; ++d)
2752  {
2753  switch (dim)
2754  {
2755  case 1:
2756  tmp[0][0] = jac[0][0] * hessians_quad[0][q_point];
2757  break;
2758  case 2:
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]);
2763  break;
2764  case 3:
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]);
2774  break;
2775  default:
2776  Assert (false, ExcNotImplemented());
2777  }
2778  }
2779  }
2780 }
2781 
2782 
2783 
2784 template <int dim, int dofs_per_cell_, int n_q_points_,
2785  int n_components_, typename Number>
2786 inline
2789 ::get_hessian (const unsigned int q_point) const
2790 {
2791  Assert (this->hessians_quad_initialized==true,
2792  internal::ExcAccessToUninitializedField());
2793  AssertIndexRange (q_point, n_q_points);
2794 
2796 
2797  // Cartesian cell
2798  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
2799  {
2800  const Tensor<1,dim,VectorizedArray<Number> > &jac = cartesian_data[0];
2801  for (unsigned int comp=0; comp<n_components; comp++)
2802  for (unsigned int d=0; d<dim; ++d)
2803  {
2804  hessian_out[comp][d][d] = (this->hessians_quad[comp][d][q_point] *
2805  jac[d] * jac[d]);
2806  switch (dim)
2807  {
2808  case 1:
2809  break;
2810  case 2:
2811  hessian_out[comp][0][1] = (this->hessians_quad[comp][2][q_point] *
2812  jac[0] * jac[1]);
2813  break;
2814  case 3:
2815  hessian_out[comp][0][1] = (this->hessians_quad[comp][3][q_point] *
2816  jac[0] * jac[1]);
2817  hessian_out[comp][0][2] = (this->hessians_quad[comp][4][q_point] *
2818  jac[0] * jac[2]);
2819  hessian_out[comp][1][2] = (this->hessians_quad[comp][5][q_point] *
2820  jac[1] * jac[2]);
2821  break;
2822  default:
2823  Assert (false, ExcNotImplemented());
2824  }
2825  for (unsigned int e=d+1; e<dim; ++e)
2826  hessian_out[comp][e][d] = hessian_out[comp][d][e];
2827  }
2828  }
2829  // cell with general Jacobian
2830  else if (this->cell_type == internal::MatrixFreeFunctions::general)
2831  {
2832  Assert (this->mapping_info.second_derivatives_initialized == true,
2833  ExcNotInitialized());
2834  const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[q_point];
2835  const Tensor<2,dim,VectorizedArray<Number> > &jac_grad = jacobian_grad[q_point];
2836  const Tensor<1,(dim>1?dim*(dim-1)/2:1),
2838  & jac_grad_UT = jacobian_grad_upper[q_point];
2839  for (unsigned int comp=0; comp<n_components; comp++)
2840  {
2841  // compute laplacian before the gradient because it needs to access
2842  // unscaled gradient data
2843  VectorizedArray<Number> tmp[dim][dim];
2844  internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
2845  q_point, tmp);
2846 
2847  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
2848  for (unsigned int d=0; d<dim; ++d)
2849  for (unsigned int e=d; e<dim; ++e)
2850  {
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];
2854  }
2855 
2856  // add diagonal part of J' * grad(u)
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]);
2861 
2862  // add off-diagonal part of J' * grad(u)
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]);
2868 
2869  // take symmetric part
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];
2873  }
2874  }
2875  // cell with general Jacobian, but constant within the cell
2876  else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
2877  {
2878  const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[0];
2879  for (unsigned int comp=0; comp<n_components; comp++)
2880  {
2881  // compute laplacian before the gradient because it needs to access
2882  // unscaled gradient data
2883  VectorizedArray<Number> tmp[dim][dim];
2884  internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
2885  q_point, tmp);
2886 
2887  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
2888  for (unsigned int d=0; d<dim; ++d)
2889  for (unsigned int e=d; e<dim; ++e)
2890  {
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];
2894  }
2895 
2896  // no J' * grad(u) part here because the Jacobian is constant
2897  // throughout the cell and hence, its derivative is zero
2898 
2899  // take symmetric part
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];
2903  }
2904  }
2906 }
2907 
2908 
2909 
2910 template <int dim, int dofs_per_cell_, int n_q_points_,
2911  int n_components_, typename Number>
2912 inline
2915 ::get_hessian_diagonal (const unsigned int q_point) const
2916 {
2917  Assert (this->hessians_quad_initialized==true,
2918  internal::ExcAccessToUninitializedField());
2919  AssertIndexRange (q_point, n_q_points);
2920 
2922 
2923  // Cartesian cell
2924  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
2925  {
2926  const Tensor<1,dim,VectorizedArray<Number> > &jac = cartesian_data[0];
2927  for (unsigned int comp=0; comp<n_components; comp++)
2928  for (unsigned int d=0; d<dim; ++d)
2929  hessian_out[comp][d] = (this->hessians_quad[comp][d][q_point] *
2930  jac[d] * jac[d]);
2931  }
2932  // cell with general Jacobian
2933  else if (this->cell_type == internal::MatrixFreeFunctions::general)
2934  {
2935  Assert (this->mapping_info.second_derivatives_initialized == true,
2936  ExcNotInitialized());
2937  const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[q_point];
2938  const Tensor<2,dim,VectorizedArray<Number> > &jac_grad = jacobian_grad[q_point];
2939  for (unsigned int comp=0; comp<n_components; comp++)
2940  {
2941  // compute laplacian before the gradient because it needs to access
2942  // unscaled gradient data
2943  VectorizedArray<Number> tmp[dim][dim];
2944  internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
2945  q_point, tmp);
2946 
2947  // compute only the trace part of hessian, J * tmp = J *
2948  // hess_unit(u) * J^T
2949  for (unsigned int d=0; d<dim; ++d)
2950  {
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];
2954  }
2955 
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]);
2960  }
2961  }
2962  // cell with general Jacobian, but constant within the cell
2963  else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
2964  {
2965  const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[0];
2966  for (unsigned int comp=0; comp<n_components; comp++)
2967  {
2968  // compute laplacian before the gradient because it needs to access
2969  // unscaled gradient data
2970  VectorizedArray<Number> tmp[dim][dim];
2971  internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
2972  q_point, tmp);
2973 
2974  // compute only the trace part of hessian, J * tmp = J *
2975  // hess_unit(u) * J^T
2976  for (unsigned int d=0; d<dim; ++d)
2977  {
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];
2981  }
2982  }
2983  }
2984  return hessian_out;
2985 }
2986 
2987 
2988 
2989 template <int dim, int dofs_per_cell_, int n_q_points_,
2990  int n_components_, typename Number>
2991 inline
2994 ::get_laplacian (const unsigned int q_point) const
2995 {
2996  Assert (this->hessians_quad_initialized==true,
2997  internal::ExcAccessToUninitializedField());
2998  AssertIndexRange (q_point, n_q_points);
2999  Tensor<1,n_components_,VectorizedArray<Number> > laplacian_out (false);
3001  = get_hessian_diagonal(q_point);
3002  for (unsigned int comp=0; comp<n_components; ++comp)
3003  {
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];
3007  }
3008  return laplacian_out;
3009 }
3010 
3011 
3012 
3013 template <int dim, int dofs_per_cell_, int n_q_points_,
3014  int n_components_, typename Number>
3015 inline
3016 void
3018 ::submit_dof_value (const Tensor<1,n_components_,VectorizedArray<Number> > val_in,
3019  const unsigned int dof)
3020 {
3021 #ifdef DEBUG
3022  this->dof_values_initialized = true;
3023 #endif
3024  AssertIndexRange (dof, dofs_per_cell);
3025  for (unsigned int comp=0; comp<n_components; comp++)
3026  this->values_dofs[comp][dof] = val_in[comp];
3027 }
3028 
3029 
3030 
3031 template <int dim, int dofs_per_cell_, int n_q_points_,
3032  int n_components_, typename Number>
3033 inline
3034 void
3036 ::submit_value (const Tensor<1,n_components_,VectorizedArray<Number> > val_in,
3037  const unsigned int q_point)
3038 {
3039 #ifdef DEBUG
3041  AssertIndexRange (q_point, n_q_points);
3042  this->values_quad_submitted = true;
3043 #endif
3044  if (this->cell_type == internal::MatrixFreeFunctions::general)
3045  {
3046  const VectorizedArray<Number> JxW = J_value[q_point];
3047  for (unsigned int comp=0; comp<n_components; ++comp)
3048  this->values_quad[comp][q_point] = val_in[comp] * JxW;
3049  }
3050  else //if (this->cell_type < internal::MatrixFreeFunctions::general)
3051  {
3052  const VectorizedArray<Number> JxW = J_value[0] * quadrature_weights[q_point];
3053  for (unsigned int comp=0; comp<n_components; ++comp)
3054  this->values_quad[comp][q_point] = val_in[comp] * JxW;
3055  }
3056 }
3057 
3058 
3059 
3060 template <int dim, int dofs_per_cell_, int n_q_points_,
3061  int n_components_, typename Number>
3062 inline
3063 void
3065 ::submit_gradient (const Tensor<1,n_components_,
3066  Tensor<1,dim,VectorizedArray<Number> > >grad_in,
3067  const unsigned int q_point)
3068 {
3069 #ifdef DEBUG
3071  AssertIndexRange (q_point, n_q_points);
3072  this->gradients_quad_submitted = true;
3073 #endif
3074  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
3075  {
3076  const VectorizedArray<Number> JxW = J_value[0] * quadrature_weights[q_point];
3077  for (unsigned int comp=0; comp<n_components; comp++)
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);
3081  }
3082  else
3083  {
3085  this->cell_type == internal::MatrixFreeFunctions::general ?
3086  jacobian[q_point] : jacobian[0];
3087  const VectorizedArray<Number> JxW =
3088  this->cell_type == internal::MatrixFreeFunctions::general ?
3089  J_value[q_point] : J_value[0] * quadrature_weights[q_point];
3090  for (unsigned int comp=0; comp<n_components; ++comp)
3091  for (unsigned int d=0; d<dim; ++d)
3092  {
3093  VectorizedArray<Number> new_val = jac[0][d] * grad_in[comp][0];
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;
3097  }
3098  }
3099 }
3100 
3101 
3102 
3103 template <int dim, int dofs_per_cell_, int n_q_points_,
3104  int n_components_, typename Number>
3105 inline
3108 ::integrate_value () const
3109 {
3110 #ifdef DEBUG
3112  Assert (this->values_quad_submitted == true,
3113  internal::ExcAccessToUninitializedField());
3114 #endif
3115  Tensor<1,n_components_,VectorizedArray<Number> > return_value (false);
3116  for (unsigned int comp=0; comp<n_components; ++comp)
3117  return_value[comp] = this->values_quad[comp][0];
3118  for (unsigned int q=1; q<n_q_points; ++q)
3119  for (unsigned int comp=0; comp<n_components; ++comp)
3120  return_value[comp] += this->values_quad[comp][q];
3121  return (return_value);
3122 }
3123 
3124 
3125 
3126 /*----------------------- FEEvaluationAccess --------------------------------*/
3127 
3128 
3129 template <int dim, int dofs_per_cell_, int n_q_points_,
3130  int n_components_, typename Number>
3131 inline
3134  const unsigned int fe_no,
3135  const unsigned int quad_no_in)
3136  :
3137  FEEvaluationBase <dim,dofs_per_cell_,n_q_points_,n_components_,Number>
3138  (data_in, fe_no, quad_no_in)
3139 {}
3140 
3141 
3142 
3143 
3144 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
3145 
3146 
3147 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3148 inline
3151  const unsigned int fe_no,
3152  const unsigned int quad_no_in)
3153  :
3154  FEEvaluationBase <dim,dofs_per_cell_,n_q_points_,1,Number>
3155  (data_in, fe_no, quad_no_in)
3156 {}
3157 
3158 
3159 
3160 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3161 inline
3164 ::get_dof_value (const unsigned int dof) const
3165 {
3166  AssertIndexRange (dof, dofs_per_cell);
3167  return this->values_dofs[0][dof];
3168 }
3169 
3170 
3171 
3172 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3173 inline
3176 ::get_value (const unsigned int q_point) const
3177 {
3178  Assert (this->values_quad_initialized==true,
3179  internal::ExcAccessToUninitializedField());
3180  AssertIndexRange (q_point, n_q_points);
3181  return this->values_quad[0][q_point];
3182 }
3183 
3184 
3185 
3186 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3187 inline
3190 ::get_gradient (const unsigned int q_point) const
3191 {
3192  // could use the base class gradient, but that involves too many inefficient
3193  // initialization operations on tensors
3194 
3195  Assert (this->gradients_quad_initialized==true,
3196  internal::ExcAccessToUninitializedField());
3197  AssertIndexRange (q_point, n_q_points);
3198 
3199  Tensor<1,dim,VectorizedArray<Number> > grad_out (false);
3200 
3201  // Cartesian cell
3202  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
3203  {
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]);
3207  }
3208  // cell with general/constant Jacobian
3209  else
3210  {
3212  this->cell_type == internal::MatrixFreeFunctions::general ?
3213  this->jacobian[q_point] : this->jacobian[0];
3214  for (unsigned int d=0; d<dim; ++d)
3215  {
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]);
3219  }
3220  }
3221  return grad_out;
3222 }
3223 
3224 
3225 
3226 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3227 inline
3230 ::get_hessian (const unsigned int q_point) const
3231 {
3232  return BaseClass::get_hessian(q_point)[0];
3233 }
3234 
3235 
3236 
3237 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3238 inline
3241 ::get_hessian_diagonal (const unsigned int q_point) const
3242 {
3243  return BaseClass::get_hessian_diagonal(q_point)[0];
3244 }
3245 
3246 
3247 
3248 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3249 inline
3252 ::get_laplacian (const unsigned int q_point) const
3253 {
3254  return BaseClass::get_laplacian(q_point)[0];
3255 }
3256 
3257 
3258 
3259 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3260 inline
3261 void
3264  const unsigned int dof)
3265 {
3266 #ifdef DEBUG
3267  this->dof_values_initialized = true;
3268  AssertIndexRange (dof, dofs_per_cell);
3269 #endif
3270  this->values_dofs[0][dof] = val_in;
3271 }
3272 
3273 
3274 
3275 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3276 inline
3277 void
3280  const unsigned int q_point)
3281 {
3282 #ifdef DEBUG
3284  AssertIndexRange (q_point, n_q_points);
3285  this->values_quad_submitted = true;
3286 #endif
3287  if (this->cell_type == internal::MatrixFreeFunctions::general)
3288  {
3289  const VectorizedArray<Number> JxW = this->J_value[q_point];
3290  this->values_quad[0][q_point] = val_in * JxW;
3291  }
3292  else //if (this->cell_type < internal::MatrixFreeFunctions::general)
3293  {
3294  const VectorizedArray<Number> JxW = this->J_value[0] * this->quadrature_weights[q_point];
3295  this->values_quad[0][q_point] = val_in * JxW;
3296  }
3297 }
3298 
3299 
3300 
3301 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3302 inline
3303 void
3305 ::submit_gradient (const Tensor<1,dim,VectorizedArray<Number> > grad_in,
3306  const unsigned int q_point)
3307 {
3308 #ifdef DEBUG
3310  AssertIndexRange (q_point, n_q_points);
3311  this->gradients_quad_submitted = true;
3312 #endif
3313  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
3314  {
3315  const VectorizedArray<Number> JxW = this->J_value[0] * this->quadrature_weights[q_point];
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] *
3319  JxW);
3320  }
3321  // general/affine cell type
3322  else
3323  {
3325  this->cell_type == internal::MatrixFreeFunctions::general ?
3326  this->jacobian[q_point] : this->jacobian[0];
3327  const VectorizedArray<Number> JxW =
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)
3331  {
3332  VectorizedArray<Number> new_val = jac[0][d] * grad_in[0];
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;
3336  }
3337  }
3338 }
3339 
3340 
3341 
3342 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3343 inline
3346 ::integrate_value () const
3347 {
3348  return BaseClass::integrate_value()[0];
3349 }
3350 
3351 
3352 
3353 
3354 /*----------------- FEEvaluationAccess vector-valued ------------------------*/
3355 
3356 
3357 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3358 inline
3361  const unsigned int fe_no,
3362  const unsigned int quad_no_in)
3363  :
3364  FEEvaluationBase <dim,dofs_per_cell_,n_q_points_,dim,Number>
3365  (data_in, fe_no, quad_no_in)
3366 {}
3367 
3368 
3369 
3370 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3371 inline
3374 ::get_gradient (const unsigned int q_point) const
3375 {
3376  return BaseClass::get_gradient (q_point);
3377 }
3378 
3379 
3380 
3381 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3382 inline
3385 ::get_divergence (const unsigned int q_point) const
3386 {
3387  Assert (this->gradients_quad_initialized==true,
3388  internal::ExcAccessToUninitializedField());
3389  AssertIndexRange (q_point, n_q_points);
3390 
3391  VectorizedArray<Number> divergence;
3392 
3393  // Cartesian cell
3394  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
3395  {
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]);
3401  }
3402  // cell with general/constant Jacobian
3403  else
3404  {
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]);
3414  }
3415  return divergence;
3416 }
3417 
3418 
3419 
3420 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3421 inline
3424 ::get_symmetric_gradient (const unsigned int q_point) const
3425 {
3426  // copy from generic function into dim-specialization function
3427  const Tensor<2,dim,VectorizedArray<Number> > grad = get_gradient(q_point);
3428  VectorizedArray<Number> symmetrized [(dim*dim+dim)/2];
3430  for (unsigned int d=0; d<dim; ++d)
3431  symmetrized[d] = grad[d][d];
3432  switch (dim)
3433  {
3434  case 1:
3435  break;
3436  case 2:
3437  symmetrized[2] = grad[0][1] + grad[1][0];
3438  symmetrized[2] *= half;
3439  break;
3440  case 3:
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;
3447  break;
3448  default:
3449  Assert (false, ExcNotImplemented());
3450  }
3451  return SymmetricTensor<2,dim,VectorizedArray<Number> > (symmetrized);
3452 }
3453 
3454 
3455 
3456 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3457 inline
3460 ::get_curl (const unsigned int q_point) const
3461 {
3462  // copy from generic function into dim-specialization function
3463  const Tensor<2,dim,VectorizedArray<Number> > grad = get_gradient(q_point);
3465  switch (dim)
3466  {
3467  case 1:
3468  Assert (false,
3469  ExcMessage("Computing the curl in 1d is not a useful operation"));
3470  break;
3471  case 2:
3472  curl[0] = grad[1][0] - grad[0][1];
3473  break;
3474  case 3:
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];
3478  break;
3479  default:
3480  Assert (false, ExcNotImplemented());
3481  }
3482  return curl;
3483 }
3484 
3485 
3486 
3487 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3488 inline
3491 ::get_hessian_diagonal (const unsigned int q_point) const
3492 {
3493  Assert (this->hessians_quad_initialized==true,
3494  internal::ExcAccessToUninitializedField());
3495  AssertIndexRange (q_point, n_q_points);
3496 
3497  return BaseClass::get_hessian_diagonal (q_point);
3498 }
3499 
3500 
3501 
3502 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3503 inline
3506 ::get_hessian (const unsigned int q_point) const
3507 {
3508  Assert (this->hessians_quad_initialized==true,
3509  internal::ExcAccessToUninitializedField());
3510  AssertIndexRange (q_point, n_q_points);
3511  return BaseClass::get_hessian(q_point);
3512 }
3513 
3514 
3515 
3516 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3517 inline
3518 void
3520 ::submit_gradient (const Tensor<2,dim,VectorizedArray<Number> > grad_in,
3521  const unsigned int q_point)
3522 {
3523  BaseClass::submit_gradient (grad_in, q_point);
3524 }
3525 
3526 
3527 
3528 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3529 inline
3530 void
3533  grad_in,
3534  const unsigned int q_point)
3535 {
3536  BaseClass::submit_gradient(grad_in, q_point);
3537 }
3538 
3539 template <int dim, int dofs_per_cell_, int n_q_points_,
3540  typename Number>
3541 inline
3542 void
3545  const unsigned int q_point)
3546 {
3547 #ifdef DEBUG
3549  AssertIndexRange (q_point, n_q_points);
3550  this->gradients_quad_submitted = true;
3551 #endif
3552  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
3553  {
3554  const VectorizedArray<Number> fac = this->J_value[0] *
3555  this->quadrature_weights[q_point] * div_in;
3556  for (unsigned int d=0; d<dim; ++d)
3557  {
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)
3561  {
3562  this->gradients_quad[d][e][q_point] = VectorizedArray<Number>();
3563  this->gradients_quad[e][d][q_point] = VectorizedArray<Number>();
3564  }
3565  }
3566  }
3567  else
3568  {
3570  this->cell_type == internal::MatrixFreeFunctions::general ?
3571  this->jacobian[q_point] : this->jacobian[0];
3572  const VectorizedArray<Number> fac =
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)
3577  {
3578  for (unsigned int e=0; e<dim; ++e)
3579  this->gradients_quad[d][e][q_point] = jac[d][e] * fac;
3580  }
3581  }
3582 }
3583 
3584 
3585 
3586 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3587 inline
3588 void
3591  sym_grad,
3592  const unsigned int q_point)
3593 {
3594  // could have used base class operator, but that involves some overhead
3595  // which is inefficient. it is nice to have the symmetric tensor because
3596  // that saves some operations
3597 #ifdef DEBUG
3599  AssertIndexRange (q_point, n_q_points);
3600  this->gradients_quad_submitted = true;
3601 #endif
3602  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
3603  {
3604  const VectorizedArray<Number> JxW = this->J_value[0] * this->quadrature_weights[q_point];
3605  for (unsigned int d=0; d<dim; ++d)
3606  this->gradients_quad[d][d][q_point] = (sym_grad.access_raw_entry(d) *
3607  JxW *
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)
3611  {
3612  const VectorizedArray<Number> value = sym_grad.access_raw_entry(counter) * JxW;
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]);
3617  }
3618  }
3619  // general/affine cell type
3620  else
3621  {
3622  const VectorizedArray<Number> JxW =
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];
3628  VectorizedArray<Number> weighted [dim][dim];
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)
3633  {
3634  const VectorizedArray<Number> value = sym_grad.access_raw_entry(counter) * JxW;
3635  weighted[i][j] = value;
3636  weighted[j][i] = value;
3637  }
3638  for (unsigned int comp=0; comp<dim; ++comp)
3639  for (unsigned int d=0; d<dim; ++d)
3640  {
3641  VectorizedArray<Number> new_val = jac[0][d] * weighted[comp][0];
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;
3645  }
3646  }
3647 }
3648 
3649 
3650 
3651 template <int dim, int dofs_per_cell_, int n_q_points_, typename Number>
3652 inline
3653 void
3655 ::submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray<Number> > curl,
3656  const unsigned int q_point)
3657 {
3659  switch (dim)
3660  {
3661  case 1:
3662  Assert (false,
3663  ExcMessage("Testing by the curl in 1d is not a useful operation"));
3664  break;
3665  case 2:
3666  grad[1][0] = curl[0];
3667  grad[0][1] = -curl[0];
3668  break;
3669  case 3:
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];
3676  break;
3677  default:
3678  Assert (false, ExcNotImplemented());
3679  }
3680  submit_gradient (grad, q_point);
3681 }
3682 
3683 
3684 
3685 /*----------------------- FEEvaluationGeneral -------------------------------*/
3686 
3687 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
3688  typename Number>
3689 inline
3692  const unsigned int fe_no,
3693  const unsigned int quad_no)
3694  :
3695  BaseClass (data_in, fe_no, quad_no)
3696 {
3697 #ifdef DEBUG
3698  // print error message when the dimensions do not match. Propose a possible
3699  // fix
3700  if (dofs_per_cell != this->data.dofs_per_cell ||
3701  n_q_points != this->data.n_q_points)
3702  {
3703  std::string message =
3704  "-------------------------------------------------------\n";
3705  message += "Illegal arguments in constructor/wrong template arguments!\n";
3706  message += " Called --> FEEvaluation<dim,";
3707  message += Utilities::int_to_string(fe_degree) + ",";
3708  message += Utilities::int_to_string(n_q_points_1d);
3709  message += "," + Utilities::int_to_string(n_components);
3710  message += ",Number>(data, ";
3711  message += Utilities::int_to_string(fe_no) + ", ";
3712  message += Utilities::int_to_string(quad_no) + ")\n";
3713 
3714  // check whether some other vector component has the correct number of
3715  // points
3716  unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
3717  proposed_quad_comp = numbers::invalid_unsigned_int;
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;
3720  else
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]
3723  == dofs_per_cell)
3724  {
3725  proposed_dof_comp = no;
3726  break;
3727  }
3728  if (n_q_points ==
3729  this->mapping_info.mapping_data_gen[quad_no].n_q_points[this->active_quad_index])
3730  proposed_quad_comp = quad_no;
3731  else
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]
3734  == n_q_points)
3735  {
3736  proposed_quad_comp = no;
3737  break;
3738  }
3739  if (proposed_dof_comp != numbers::invalid_unsigned_int &&
3740  proposed_quad_comp != numbers::invalid_unsigned_int)
3741  {
3742  if (proposed_dof_comp != fe_no)
3743  message += "Wrong vector component selection:\n";
3744  else
3745  message += "Wrong quadrature formula selection:\n";
3746  message += " Did you mean FEEvaluation<dim,";
3747  message += Utilities::int_to_string(fe_degree) + ",";
3748  message += Utilities::int_to_string(n_q_points_1d);
3749  message += "," + Utilities::int_to_string(n_components);
3750  message += ",Number>(data, ";
3751  message += Utilities::int_to_string(proposed_dof_comp) + ", ";
3752  message += Utilities::int_to_string(proposed_quad_comp) + ")?\n";
3753  std::string correct_pos;
3754  if (proposed_dof_comp != fe_no)
3755  correct_pos = " ^ ";
3756  else
3757  correct_pos = " ";
3758  if (proposed_quad_comp != quad_no)
3759  correct_pos += " ^\n";
3760  else
3761  correct_pos += " \n";
3762  message += " " + correct_pos;
3763  }
3764  // ok, did not find the numbers specified by the template arguments in
3765  // the given list. Suggest correct template arguments
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,";
3770  message += Utilities::int_to_string(proposed_fe_degree) + ",";
3771  message += Utilities::int_to_string(proposed_n_q_points_1d);
3772  message += "," + Utilities::int_to_string(n_components);
3773  message += ",Number>(data, ";
3774  message += Utilities::int_to_string(fe_no) + ", ";
3775  message += Utilities::int_to_string(quad_no) + ")?\n";
3776  std::string correct_pos;
3777  if (proposed_fe_degree != fe_degree)
3778  correct_pos = " ^";
3779  else
3780  correct_pos = " ";
3781  if (proposed_n_q_points_1d != n_q_points_1d)
3782  correct_pos += " ^\n";
3783  else
3784  correct_pos += " \n";
3785  message += " " + correct_pos;
3786 
3787  Assert (dofs_per_cell == this->data.dofs_per_cell &&
3788  n_q_points == this->data.n_q_points,
3789  ExcMessage(message));
3790  }
3791  AssertDimension (n_q_points,
3792  this->mapping_info.mapping_data_gen[this->quad_no].
3793  n_q_points[this->active_quad_index]);
3794  AssertDimension (dofs_per_cell * this->n_fe_components,
3795  this->dof_info.dofs_per_cell[this->active_fe_index]);
3796 #endif
3797 }
3798 
3799 
3800 
3801 namespace internal
3802 {
3803  // evaluates the given shape data in 1d-3d using the tensor product
3804  // form. does not use a particular layout of entries in the matrices
3805  // like the functions below and corresponds to a usual matrix-matrix
3806  // product
3807  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
3808  int direction, bool dof_to_quad, bool add>
3809  inline
3810  void
3811  apply_tensor_product (const Number *shape_data,
3812  const Number in [],
3813  Number out [])
3814  {
3815  AssertIndexRange (direction, dim);
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);
3818 
3819  const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
3820  const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
3822 
3823  for (int i2=0; i2<n_blocks2; ++i2)
3824  {
3825  for (int i1=0; i1<n_blocks1; ++i1)
3826  {
3827  for (int col=0; col<nn; ++col)
3828  {
3829  Number val0;
3830  if (dof_to_quad == true)
3831  val0 = shape_data[col];
3832  else
3833  val0 = shape_data[col*n_q_points_1d];
3834  Number res0 = val0 * in[0];
3835  for (int ind=1; ind<mm; ++ind)
3836  {
3837  if (dof_to_quad == true)
3838  val0 = shape_data[ind*n_q_points_1d+col];
3839  else
3840  val0 = shape_data[col*n_q_points_1d+ind];
3841  res0 += val0 * in[stride*ind];
3842  }
3843  if (add == false)
3844  out[stride*col] = res0;
3845  else
3846  out[stride*col] += res0;
3847  }
3848 
3849  // increment: in regular case, just go to the next point in
3850  // x-direction. If we are at the end of one chunk in x-dir, need
3851  // to jump over to the next layer in z-direction
3852  switch (direction)
3853  {
3854  case 0:
3855  in += mm;
3856  out += nn;
3857  break;
3858  case 1:
3859  case 2:
3860  ++in;
3861  ++out;
3862  break;
3863  default:
3864  Assert (false, ExcNotImplemented());
3865  }
3866  }
3867  if (direction == 1)
3868  {
3869  in += nn*(mm-1);
3870  out += nn*(nn-1);
3871  }
3872  }
3873  }
3874 
3875 
3876 
3877  // This method applies the tensor product operation to produce face values
3878  // out from cell values. As opposed to the apply_tensor_product method, this
3879  // method assumes that the directions orthogonal to the face have
3880  // fe_degree+1 degrees of freedom per direction and not n_q_points_1d for
3881  // those directions lower than the one currently applied
3882  template <int dim, int fe_degree, typename Number, int face_direction,
3883  bool dof_to_quad, bool add>
3884  inline
3885  void
3886  apply_tensor_product_face (const Number *shape_data,
3887  const Number in [],
3888  Number out [])
3889  {
3890  const int n_blocks1 = dim > 1 ? (fe_degree+1) : 1;
3891  const int n_blocks2 = dim > 2 ? (fe_degree+1) : 1;
3892 
3893  AssertIndexRange (face_direction, dim);
3894  const int mm = dof_to_quad ? (fe_degree+1) : 1,
3895  nn = dof_to_quad ? 1 : (fe_degree+1);
3896 
3898 
3899  for (int i2=0; i2<n_blocks2; ++i2)
3900  {
3901  for (int i1=0; i1<n_blocks1; ++i1)
3902  {
3903  if (dof_to_quad == true)
3904  {
3905  Number res0 = shape_data[0] * in[0];
3906  for (int ind=1; ind<mm; ++ind)
3907  res0 += shape_data[ind] * in[stride*ind];
3908  if (add == false)
3909  out[0] = res0;
3910  else
3911  out[0] += res0;
3912  }
3913  else
3914  {
3915  for (int col=0; col<nn; ++col)
3916  if (add == false)
3917  out[col*stride] = shape_data[col] * in[0];
3918  else
3919  out[col*stride] += shape_data[col] * in[0];
3920  }
3921 
3922  // increment: in regular case, just go to the next point in
3923  // x-direction. If we are at the end of one chunk in x-dir, need
3924  // to jump over to the next layer in z-direction
3925  switch (face_direction)
3926  {
3927  case 0:
3928  in += mm;
3929  out += nn;
3930  break;
3931  case 1:
3932  case 2:
3933  ++in;
3934  ++out;
3935  break;
3936  default:
3937  Assert (false, ExcNotImplemented());
3938  }
3939  }
3940  if (face_direction == 1)
3941  {
3942  in += mm*(mm-1);
3943  out += nn*(nn-1);
3944  }
3945  }
3946  }
3947 
3948 
3949 
3950  // This method specializes the general application of tensor-product based
3951  // elements for "symmetric" finite elements, i.e., when the shape functions
3952  // are symmetric about 0.5 and the quadrature points are, too. In that case,
3953  // the 1D shape values read (sorted lexicographically, rows run over 1D
3954  // dofs, columns over quadrature points):
3955  // Q2 --> [ 0.687 0 -0.087 ]
3956  // [ 0.4 1 0.4 ]
3957  // [-0.087 0 0.687 ]
3958  // Q3 --> [ 0.66 0.003 0.002 0.049 ]
3959  // [ 0.521 1.005 -0.01 -0.230 ]
3960  // [-0.230 -0.01 1.005 0.521 ]
3961  // [ 0.049 0.002 0.003 0.66 ]
3962  // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
3963  // [ 0.608 1.059 0 0.039 0.176 ]
3964  // [-0.409 -0.113 1 -0.113 -0.409 ]
3965  // [ 0.176 0.039 0 1.059 0.608 ]
3966  // [-0.032 -0.007 0 0.022 0.658 ]
3967  //
3968  // In these matrices, we want to use avoid computations involving zeros and
3969  // ones and in addition use the symmetry in entries to reduce the number of
3970  // read operations.
3971  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
3972  int direction, bool dof_to_quad, bool add>
3973  inline
3974  void
3975  apply_tensor_product_values (const Number *shape_values,
3976  const Number in [],
3977  Number out [])
3978  {
3979  AssertIndexRange (direction, dim);
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;
3984 
3985  const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
3986  const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
3988 
3989  for (int i2=0; i2<n_blocks2; ++i2)
3990  {
3991  for (int i1=0; i1<n_blocks1; ++i1)
3992  {
3993  for (int col=0; col<n_cols; ++col)
3994  {
3995  Number val0, val1, in0, in1, res0, res1;
3996  if (dof_to_quad == true)
3997  {
3998  val0 = shape_values[col];
3999  val1 = shape_values[nn-1-col];
4000  }
4001  else
4002  {
4003  val0 = shape_values[col*n_q_points_1d];
4004  val1 = shape_values[(col+1)*n_q_points_1d-1];
4005  }
4006  if (mid > 0)
4007  {
4008  in0 = in[0];
4009  in1 = in[stride*(mm-1)];
4010  res0 = val0 * in0;
4011  res1 = val1 * in0;
4012  res0 += val1 * in1;
4013  res1 += val0 * in1;
4014  for (int ind=1; ind<mid; ++ind)
4015  {
4016  if (dof_to_quad == true)
4017  {
4018  val0 = shape_values[ind*n_q_points_1d+col];
4019  val1 = shape_values[ind*n_q_points_1d+nn-1-col];
4020  }
4021  else
4022  {
4023  val0 = shape_values[col*n_q_points_1d+ind];
4024  val1 = shape_values[(col+1)*n_q_points_1d-1-ind];
4025  }
4026  in0 = in[stride*ind];
4027  in1 = in[stride*(mm-1-ind)];
4028  res0 += val0 * in0;
4029  res1 += val1 * in0;
4030  res0 += val1 * in1;
4031  res1 += val0 * in1;
4032  }
4033  }
4034  else
4035  res0 = res1 = Number();
4036  if (dof_to_quad == true)
4037  {
4038  if (mm % 2 == 1)
4039  {
4040  val0 = shape_values[mid*n_q_points_1d+col];
4041  val1 = val0 * in[stride*mid];
4042  res0 += val1;
4043  res1 += val1;
4044  }
4045  }
4046  else
4047  {
4048  if (mm % 2 == 1 && nn % 2 == 0)
4049  {
4050  val0 = shape_values[col*n_q_points_1d+mid];
4051  val1 = val0 * in[stride*mid];
4052  res0 += val1;
4053  res1 += val1;
4054  }
4055  }
4056  if (add == false)
4057  {
4058  out[stride*col] = res0;
4059  out[stride*(nn-1-col)] = res1;
4060  }
4061  else
4062  {
4063  out[stride*col] += res0;
4064  out[stride*(nn-1-col)] += res1;
4065  }
4066  }
4067  if ( dof_to_quad == true && nn%2==1 && mm%2==1 )
4068  {
4069  if (add==false)
4070  out[stride*n_cols] = in[stride*mid];
4071  else
4072  out[stride*n_cols] += in[stride*mid];
4073  }
4074  else if (dof_to_quad == true && nn%2==1)
4075  {
4076  Number res0;
4077  Number val0 = shape_values[n_cols];
4078  if (mid > 0)
4079  {
4080  res0 = in[0] + in[stride*(mm-1)];
4081  res0 *= val0;
4082  for (int ind=1; ind<mid; ++ind)
4083  {
4084  val0 = shape_values[ind*n_q_points_1d+n_cols];
4085  Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
4086  val1 *= val0;
4087  res0 += val1;
4088  }
4089  }
4090  else
4091  res0 = Number();
4092  if (add == false)
4093  out[stride*n_cols] = res0;
4094  else
4095  out[stride*n_cols] += res0;
4096  }
4097  else if (dof_to_quad == false && nn%2 == 1)
4098  {
4099  Number res0;
4100  if (mid > 0)
4101  {
4102  Number val0 = shape_values[n_cols*n_q_points_1d];
4103  res0 = in[0] + in[stride*(mm-1)];
4104  res0 *= val0;
4105  for (int ind=1; ind<mid; ++ind)
4106  {
4107  val0 = shape_values[n_cols*n_q_points_1d+ind];
4108  Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
4109  val1 *= val0;
4110  res0 += val1;
4111  }
4112  if (mm % 2)
4113  res0 += in[stride*mid];
4114  }
4115  else
4116  res0 = in[0];
4117  if (add == false)
4118  out[stride*n_cols] = res0;
4119  else
4120  out[stride*n_cols] += res0;
4121  }
4122 
4123  // increment: in regular case, just go to the next point in
4124  // x-direction. If we are at the end of one chunk in x-dir, need to
4125  // jump over to the next layer in z-direction
4126  switch (direction)
4127  {
4128  case 0:
4129  in += mm;
4130  out += nn;
4131  break;
4132  case 1:
4133  case 2:
4134  ++in;
4135  ++out;
4136  break;
4137  default:
4138  Assert (false, ExcNotImplemented());
4139  }
4140  }
4141  if (direction == 1)
4142  {
4143  in += nn*(mm-1);
4144  out += nn*(nn-1);
4145  }
4146  }
4147  }
4148 
4149 
4150 
4151  // evaluates the given shape data in 1d-3d using the tensor product
4152  // form assuming the symmetries of unit cell shape gradients for
4153  // finite elements in FEEvaluation
4154 
4155  // For the specialized loop used for the gradient computation in
4156  // here, the 1D shape values read (sorted lexicographically, rows
4157  // run over 1D dofs, columns over quadrature points):
4158  // Q2 --> [-2.549 -1 0.549 ]
4159  // [ 3.098 0 -3.098 ]
4160  // [-0.549 1 2.549 ]
4161  // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
4162  // [ 6.07 -1.44 -2.97 2.196 ]
4163  // [-2.196 2.97 1.44 -6.07 ]
4164  // [ 0.44 -0.5 1.03 4.315 ]
4165  // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
4166  // [10.111 -2.76 -2.667 2.066 -2.306 ]
4167  // [-5.688 5.773 0 -5.773 5.688 ]
4168  // [ 2.306 -2.066 2.667 2.76 -10.111 ]
4169  // [-0.413 0.353 -0.333 -0.353 0.413 ]
4170  //
4171  // In these matrices, we want to use avoid computations involving
4172  // zeros and ones and in addition use the symmetry in entries to
4173  // reduce the number of read operations.
4174  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
4175  int direction, bool dof_to_quad, bool add>
4176  inline
4177  void
4178  apply_tensor_product_gradients (const Number *shape_gradients,
4179  const Number in [],
4180  Number out [])
4181  {
4182  AssertIndexRange (direction, dim);
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;
4187 
4188  const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
4189  const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
4191 
4192  for (int i2=0; i2<n_blocks2; ++i2)
4193  {
4194  for (int i1=0; i1<n_blocks1; ++i1)
4195  {
4196  for (int col=0; col<n_cols; ++col)
4197  {
4198  Number val0, val1, in0, in1, res0, res1;
4199  if (dof_to_quad == true)
4200  {
4201  val0 = shape_gradients[col];
4202  val1 = shape_gradients[nn-1-col];
4203  }
4204  else
4205  {
4206  val0 = shape_gradients[col*n_q_points_1d];
4207  val1 = shape_gradients[(nn-col-1)*n_q_points_1d];
4208  }
4209  if (mid > 0)
4210  {
4211  in0 = in[0];
4212  in1 = in[stride*(mm-1)];
4213  res0 = val0 * in0;
4214  res1 = val1 * in0;
4215  res0 -= val1 * in1;
4216  res1 -= val0 * in1;
4217  for (int ind=1; ind<mid; ++ind)
4218  {
4219  if (dof_to_quad == true)
4220  {
4221  val0 = shape_gradients[ind*n_q_points_1d+col];
4222  val1 = shape_gradients[ind*n_q_points_1d+nn-1-col];
4223  }
4224  else
4225  {
4226  val0 = shape_gradients[col*n_q_points_1d+ind];
4227  val1 = shape_gradients[(nn-col-1)*n_q_points_1d+ind];
4228  }
4229  in0 = in[stride*ind];
4230  in1 = in[stride*(mm-1-ind)];
4231  res0 += val0 * in0;
4232  res1 += val1 * in0;
4233  res0 -= val1 * in1;
4234  res1 -= val0 * in1;
4235  }
4236  }
4237  else
4238  res0 = res1 = Number();
4239  if (mm % 2 == 1)
4240  {
4241  if (dof_to_quad == true)
4242  val0 = shape_gradients[mid*n_q_points_1d+col];
4243  else
4244  val0 = shape_gradients[col*n_q_points_1d+mid];
4245  val1 = val0 * in[stride*mid];
4246  res0 += val1;
4247  res1 -= val1;
4248  }
4249  if (add == false)
4250  {
4251  out[stride*col] = res0;
4252  out[stride*(nn-1-col)] = res1;
4253  }
4254  else
4255  {
4256  out[stride*col] += res0;
4257  out[stride*(nn-1-col)] += res1;
4258  }
4259  }
4260  if ( nn%2 == 1 )
4261  {
4262  Number val0, res0;
4263  if (dof_to_quad == true)
4264  val0 = shape_gradients[n_cols];
4265  else
4266  val0 = shape_gradients[n_cols*n_q_points_1d];
4267  res0 = in[0] - in[stride*(mm-1)];
4268  res0 *= val0;
4269  for (int ind=1; ind<mid; ++ind)
4270  {
4271  if (dof_to_quad == true)
4272  val0 = shape_gradients[ind*n_q_points_1d+n_cols];
4273  else
4274  val0 = shape_gradients[n_cols*n_q_points_1d+ind];
4275  Number val1 = in[stride*ind] - in[stride*(mm-1-ind)];
4276  val1 *= val0;
4277  res0 += val1;
4278  }
4279  if (add == false)
4280  out[stride*n_cols] = res0;
4281  else
4282  out[stride*n_cols] += res0;
4283  }
4284 
4285  // increment: in regular case, just go to the next point in
4286  // x-direction. for y-part in 3D and if we are at the end of one
4287  // chunk in x-dir, need to jump over to the next layer in
4288  // z-direction
4289  switch (direction)
4290  {
4291  case 0:
4292  in += mm;
4293  out += nn;
4294  break;
4295  case 1:
4296  case 2:
4297  ++in;
4298  ++out;
4299  break;
4300  default:
4301  Assert (false, ExcNotImplemented());
4302  }
4303  }
4304 
4305  if (direction == 1)
4306  {
4307  in += nn * (mm-1);
4308  out += nn * (nn-1);
4309  }
4310  }
4311  }
4312 
4313 
4314 
4315  // evaluates the given shape data in 1d-3d using the tensor product
4316  // form assuming the symmetries of unit cell shape hessians for
4317  // finite elements in FEEvaluation
4318  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
4319  int direction, bool dof_to_quad, bool add>
4320  inline
4321  void
4322  apply_tensor_product_hessians (const Number *shape_hessians,
4323  const Number in [],
4324  Number out [])
4325  {
4326  AssertIndexRange (direction, dim);
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;
4331 
4332  const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
4333  const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
4335 
4336  for (int i2=0; i2<n_blocks2; ++i2)
4337  {
4338  for (int i1=0; i1<n_blocks1; ++i1)
4339  {
4340  for (int col=0; col<n_cols; ++col)
4341  {
4342  Number val0, val1, in0, in1, res0, res1;
4343  if (dof_to_quad == true)
4344  {
4345  val0 = shape_hessians[col];
4346  val1 = shape_hessians[nn-1-col];
4347  }
4348  else
4349  {
4350  val0 = shape_hessians[col*n_q_points_1d];
4351  val1 = shape_hessians[(col+1)*n_q_points_1d-1];
4352  }
4353  if (mid > 0)
4354  {
4355  in0 = in[0];
4356  in1 = in[stride*(mm-1)];
4357  res0 = val0 * in0;
4358  res1 = val1 * in0;
4359  res0 += val1 * in1;
4360  res1 += val0 * in1;
4361  for (int ind=1; ind<mid; ++ind)
4362  {
4363  if (dof_to_quad == true)
4364  {
4365  val0 = shape_hessians[ind*n_q_points_1d+col];
4366  val1 = shape_hessians[ind*n_q_points_1d+nn-1-col];
4367  }
4368  else
4369  {
4370  val0 = shape_hessians[col*n_q_points_1d+ind];
4371  val1 = shape_hessians[(col+1)*n_q_points_1d-1-ind];
4372  }
4373  in0 = in[stride*ind];
4374  in1 = in[stride*(mm-1-ind)];
4375  res0 += val0 * in0;
4376  res1 += val1 * in0;
4377  res0 += val1 * in1;
4378  res1 += val0 * in1;
4379  }
4380  }
4381  else
4382  res0 = res1 = Number();
4383  if (mm % 2 == 1)
4384  {
4385  if (dof_to_quad == true)
4386  val0 = shape_hessians[mid*n_q_points_1d+col];
4387  else
4388  val0 = shape_hessians[col*n_q_points_1d+mid];
4389  val1 = val0 * in[stride*mid];
4390  res0 += val1;
4391  res1 += val1;
4392  }
4393  if (add == false)
4394  {
4395  out[stride*col] = res0;
4396  out[stride*(nn-1-col)] = res1;
4397  }
4398  else
4399  {
4400  out[stride*col] += res0;
4401  out[stride*(nn-1-col)] += res1;
4402  }
4403  }
4404  if ( nn%2 == 1 )
4405  {
4406  Number val0, res0;
4407  if (dof_to_quad == true)
4408  val0 = shape_hessians[n_cols];
4409  else
4410  val0 = shape_hessians[n_cols*n_q_points_1d];
4411  if (mid > 0)
4412  {
4413  res0 = in[0] + in[stride*(mm-1)];
4414  res0 *= val0;
4415  for (int ind=1; ind<mid; ++ind)
4416  {
4417  if (dof_to_quad == true)
4418  val0 = shape_hessians[ind*n_q_points_1d+n_cols];
4419  else
4420  val0 = shape_hessians[n_cols*n_q_points_1d+ind];
4421  Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
4422  val1 *= val0;
4423  res0 += val1;
4424  }
4425  }
4426  else
4427  res0 = Number();
4428  if (mm % 2 == 1)
4429  {
4430  if (dof_to_quad == true)
4431  val0 = shape_hessians[mid*n_q_points_1d+n_cols];
4432  else
4433  val0 = shape_hessians[n_cols*n_q_points_1d+mid];
4434  res0 += val0 * in[stride*mid];
4435  }
4436  if (add == false)
4437  out[stride*n_cols] = res0;
4438  else
4439  out[stride*n_cols] += res0;
4440  }
4441 
4442  // increment: in regular case, just go to the next point in
4443  // x-direction. If we are at the end of one chunk in x-dir, need to
4444  // jump over to the next layer in z-direction
4445  switch (direction)
4446  {
4447  case 0:
4448  in += mm;
4449  out += nn;
4450  break;
4451  case 1:
4452  case 2:
4453  ++in;
4454  ++out;
4455  break;
4456  default:
4457  Assert (false, ExcNotImplemented());
4458  }
4459  }
4460  if (direction == 1)
4461  {
4462  in += nn*(mm-1);
4463  out += nn*(nn-1);
4464  }
4465  }
4466  }
4467 
4468 
4469 
4470  // This method implements a different approach to the symmetric case for
4471  // values, gradients, and Hessians also treated with the above functions: It
4472  // is possible to reduce the cost per dimension from N^2 to N^2/2, where N
4473  // is the number of 1D dofs (there are only N^2/2 different entries in the
4474  // shape matrix, so this is plausible). The approach is based on the idea of
4475  // applying the operator on the even and odd part of the input vectors
4476  // separately, given that the shape functions evaluated on quadrature points
4477  // are symmetric. This method is presented e.g. in the book "Implementing
4478  // Spectral Methods for Partial Differential Equations" by David A. Kopriva,
4479  // Springer, 2009, section 3.5.3 (Even-Odd-Decomposition). Even though the
4480  // experiments in the book say that the method is not efficient for N<20, it
4481  // is more efficient in the context where the loop bounds are compile-time
4482  // constants (templates).
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>
4485  inline
4486  void
4487  apply_tensor_product_evenodd (const Number shapes [][(n_q_points_1d+1)/2],
4488  const Number in [],
4489  Number out [])
4490  {
4491  AssertIndexRange (type, 3);
4492  AssertIndexRange (direction, dim);
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;
4497 
4498  const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
4499  const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
4501 
4502  // this code may look very inefficient at first sight due to the many
4503  // different cases with if's at the innermost loop part, but all of the
4504  // conditionals can be evaluated at compile time because they are
4505  // templates, so the compiler should optimize everything away
4506  for (int i2=0; i2<n_blocks2; ++i2)
4507  {
4508  for (int i1=0; i1<n_blocks1; ++i1)
4509  {
4510  Number xp[mid], xm[mid];
4511  for (int i=0; i<mid; ++i)
4512  {
4513  if (dof_to_quad == true && type == 1)
4514  {
4515  xp[i] = in[stride*i] - in[stride*(mm-1-i)];
4516  xm[i] = in[stride*i] + in[stride*(mm-1-i)];
4517  }
4518  else
4519  {
4520  xp[i] = in[stride*i] + in[stride*(mm-1-i)];
4521  xm[i] = in[stride*i] - in[stride*(mm-1-i)];
4522  }
4523  }
4524  for (int col=0; col<n_cols; ++col)
4525  {
4526  Number r0, r1;
4527  if (mid > 0)
4528  {
4529  if (dof_to_quad == true)
4530  {
4531  r0 = shapes[0][col] * xp[0];
4532  r1 = shapes[fe_degree][col] * xm[0];
4533  }
4534  else
4535  {
4536  r0 = shapes[col][0] * xp[0];
4537  r1 = shapes[fe_degree-col][0] * xm[0];
4538  }
4539  for (int ind=1; ind<mid; ++ind)
4540  {
4541  if (dof_to_quad == true)
4542  {
4543  r0 += shapes[ind][col] * xp[ind];
4544  r1 += shapes[fe_degree-ind][col] * xm[ind];
4545  }
4546  else
4547  {
4548  r0 += shapes[col][ind] * xp[ind];
4549  r1 += shapes[fe_degree-col][ind] * xm[ind];
4550  }
4551  }
4552  }
4553  else
4554  r0 = r1 = Number();
4555  if (mm % 2 == 1 && dof_to_quad == true)
4556  {
4557  if (type == 1)
4558  r1 += shapes[mid][col] * in[stride*mid];
4559  else
4560  r0 += shapes[mid][col] * in[stride*mid];
4561  }
4562  else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0))
4563  r0 += shapes[col][mid] * in[stride*mid];
4564 
4565  if (add == false)
4566  {
4567  out[stride*col] = r0 + r1;
4568  if (type == 1 && dof_to_quad == false)
4569  out[stride*(nn-1-col)] = r1 - r0;
4570  else
4571  out[stride*(nn-1-col)] = r0 - r1;
4572  }
4573  else
4574  {
4575  out[stride*col] += r0 + r1;
4576  if (type == 1 && dof_to_quad == false)
4577  out[stride*(nn-1-col)] += r1 - r0;
4578  else
4579  out[stride*(nn-1-col)] += r0 - r1;
4580  }
4581  }
4582  if ( type == 0 && dof_to_quad == true && nn%2==1 && mm%2==1 )
4583  {
4584  if (add==false)
4585  out[stride*n_cols] = in[stride*mid];
4586  else
4587  out[stride*n_cols] += in[stride*mid];
4588  }
4589  else if (dof_to_quad == true && nn%2==1)
4590  {
4591  Number r0;
4592  if (mid > 0)
4593  {
4594  r0 = shapes[0][n_cols] * xp[0];
4595  for (int ind=1; ind<mid; ++ind)
4596  r0 += shapes[ind][n_cols] * xp[ind];
4597  }
4598  else
4599  r0 = Number();
4600  if (type != 1 && mm % 2 == 1)
4601  r0 += shapes[mid][n_cols] * in[stride*mid];
4602 
4603  if (add == false)
4604  out[stride*n_cols] = r0;
4605  else
4606  out[stride*n_cols] += r0;
4607  }
4608  else if (dof_to_quad == false && nn%2 == 1)
4609  {
4610  Number r0;
4611  if (mid > 0)
4612  {
4613  if (type == 1)
4614  {
4615  r0 = shapes[n_cols][0] * xm[0];
4616  for (int ind=1; ind<mid; ++ind)
4617  r0 += shapes[n_cols][ind] * xm[ind];
4618  }
4619  else
4620  {
4621  r0 = shapes[n_cols][0] * xp[0];
4622  for (int ind=1; ind<mid; ++ind)
4623  r0 += shapes[n_cols][ind] * xp[ind];
4624  }
4625  }
4626  else
4627  r0 = Number();
4628 
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];
4633 
4634  if (add == false)
4635  out[stride*n_cols] = r0;
4636  else
4637  out[stride*n_cols] += r0;
4638  }
4639 
4640  // increment: in regular case, just go to the next point in
4641  // x-direction. If we are at the end of one chunk in x-dir, need to
4642  // jump over to the next layer in z-direction
4643  switch (direction)
4644  {
4645  case 0:
4646  in += mm;
4647  out += nn;
4648  break;
4649  case 1:
4650  case 2:
4651  ++in;
4652  ++out;
4653  break;
4654  default:
4655  Assert (false, ExcNotImplemented());
4656  }
4657  }
4658  if (direction == 1)
4659  {
4660  in += nn*(mm-1);
4661  out += nn*(nn-1);
4662  }
4663  }
4664  }
4665 
4666 
4667 
4668  // evaluates the given shape data in 1d-3d using the tensor product
4669  // form assuming the symmetries of unit cell shape gradients for
4670  // finite elements in FEEvaluationGL
4671 
4672  // This function specializes the application of the tensor product loop for
4673  // Gauss-Lobatto elements which are symmetric about 0.5 just as the general
4674  // class of elements treated by FEEvaluation, have diagonal shape matrices
4675  // for the values and have the following gradient matrices (notice the zeros
4676  // on the diagonal in the interior points, which is due to the construction
4677  // of Legendre polynomials):
4678  // Q2 --> [-3 -1 1 ]
4679  // [ 4 0 -4 ]
4680  // [-1 1 3 ]
4681  // Q3 --> [-6 -1.618 0.618 -1 ]
4682  // [ 8.09 0 -2.236 3.09 ]
4683  // [-3.09 2.236 0 -8.09 ]
4684  // [ 1 -0.618 1.618 6 ]
4685  // Q4 --> [-10 -2.482 0.75 -0.518 1 ]
4686  // [ 13.51 0 -2.673 1.528 -2.82 ]
4687  // [-5.333 3.491 0 -3.491 5.333 ]
4688  // [ 2.82 -1.528 2.673 0 -13.51 ]
4689  // [-1 0.518 -0.75 2.482 10 ]
4690  template <int dim, int fe_degree, typename Number,
4691  int direction, bool dof_to_quad, bool add>
4692  inline
4693  void
4694  apply_tensor_product_gradients_gl (const Number *shape_gradients,
4695  const Number in [],
4696  Number out [])
4697  {
4698  AssertIndexRange (direction, dim);
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;
4703 
4704  const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
4705  const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
4707 
4708  for (int i2=0; i2<n_blocks2; ++i2)
4709  {
4710  for (int i1=0; i1<n_blocks1; ++i1)
4711  {
4712  for (int col=0; col<n_cols; ++col)
4713  {
4714  Number val0, val1, in0, in1, res0, res1;
4715  if (mid > 0)
4716  {
4717  if (dof_to_quad == true)
4718  {
4719  val0 = shape_gradients[col];
4720  val1 = shape_gradients[nn-1-col];
4721  }
4722  else
4723  {
4724  val0 = shape_gradients[col*mm];
4725  val1 = shape_gradients[(nn-col-1)*mm];
4726  }
4727  in0 = in[0];
4728  in1 = in[stride*(mm-1)];
4729  if (col == 0)
4730  {
4731  if ((mm+dof_to_quad)%2 == 1)
4732  {
4733  res0 = val0 * in0;
4734  res1 = -in0;
4735  res0 += in1;
4736  res1 -= val0 * in1;
4737  }
4738  else
4739  {
4740  res0 = val0 * in0;
4741  res0 -= in1;
4742  res1 = in0;
4743  res1 -= val0 * in1;
4744  }
4745  }
4746  else
4747  {
4748  res0 = val0 * in0;
4749  res1 = val1 * in0;
4750  res0 -= val1 * in1;
4751  res1 -= val0 * in1;
4752  }
4753  for (int ind=1; ind<mid; ++ind)
4754  {
4755  if (dof_to_quad == true)
4756  {
4757  val0 = shape_gradients[ind*mm+col];
4758  val1 = shape_gradients[ind*mm+nn-1-col];
4759  }
4760  else
4761  {
4762  val0 = shape_gradients[col*mm+ind];
4763  val1 = shape_gradients[(nn-col-1)*mm+ind];
4764  }
4765 
4766  // at inner points, the gradient is zero for ind==col
4767  in0 = in[stride*ind];
4768  in1 = in[stride*(mm-1-ind)];
4769  if (ind == col)
4770  {
4771  res1 += val1 * in0;
4772  res0 -= val1 * in1;
4773  }
4774  else
4775  {
4776  res0 += val0 * in0;
4777  res1 += val1 * in0;
4778  res0 -= val1 * in1;
4779  res1 -= val0 * in1;
4780  }
4781  }
4782  }
4783  else
4784  res0 = res1 = Number();
4785  if (mm % 2 == 1)
4786  {
4787  if (dof_to_quad == true)
4788  val0 = shape_gradients[mid*mm+col];
4789  else
4790  val0 = shape_gradients[col*mm+mid];
4791  val1 = val0 * in[stride*mid];
4792  res0 += val1;
4793  res1 -= val1;
4794  }
4795  if (add == false)
4796  {
4797  out[stride*col] = res0;
4798  out[stride*(nn-1-col)] = res1;
4799  }
4800  else
4801  {
4802  out[stride*col] += res0;
4803  out[stride*(nn-1-col)] += res1;
4804  }
4805  }
4806  if ( nn%2 == 1 )
4807  {
4808  Number val0, res0;
4809  if (dof_to_quad == true)
4810  val0 = shape_gradients[n_cols];
4811  else
4812  val0 = shape_gradients[n_cols*mm];
4813  if (mid > 0)
4814  {
4815  res0 = in[0] - in[stride*(mm-1)];
4816  res0 *= val0;
4817  for (int ind=1; ind<mid; ++ind)
4818  {
4819  if (dof_to_quad == true)
4820  val0 = shape_gradients[ind*mm+n_cols];
4821  else
4822  val0 = shape_gradients[n_cols*mm+ind];
4823  Number val1 = in[stride*ind] - in[stride*(mm-1-ind)];
4824  val1 *= val0;
4825  res0 += val1;
4826  }
4827  }
4828  else
4829  res0 = Number();
4830  if (add == false)
4831  out[stride*n_cols] = res0;
4832  else
4833  out[stride*n_cols] += res0;
4834  }
4835 
4836  // increment: in regular case, just go to the next point in
4837  // x-direction. for y-part in 3D and if we are at the end of one
4838  // chunk in x-dir, need to jump over to the next layer in
4839  // z-direction
4840  switch (direction)
4841  {
4842  case 0:
4843  in += mm;
4844  out += nn;
4845  break;
4846  case 1:
4847  case 2:
4848  ++in;
4849  ++out;
4850  break;
4851  default:
4852  Assert (false, ExcNotImplemented());
4853  }
4854  }
4855 
4856  if (direction == 1)
4857  {
4858  in += nn * (mm-1);
4859  out += nn * (nn-1);
4860  }
4861  }
4862  }
4863 
4864 
4865 
4866  // This performs the evaluation of function values, gradients and Hessians
4867  // for tensor-product finite elements. The operation is used for both
4868  // FEEvaluationGeneral and FEEvaluation, which provide different functions
4869  // apply_values, apply_gradients in the individual coordinate directions
4870  template <typename FEEval>
4871  inline
4872  void
4873  do_evaluate (FEEval &fe_eval,
4874  const bool evaluate_val,
4875  const bool evaluate_grad,
4876  const bool evaluate_lapl)
4877  {
4878  Assert (fe_eval.cell != numbers::invalid_unsigned_int,
4879  ExcNotInitialized());
4880  Assert (fe_eval.dof_values_initialized == true,
4881  internal::ExcAccessToUninitializedField());
4882 
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;
4887 
4888  for (unsigned int c=0; c<n_components; c++)
4889  {
4892 
4893  switch (dim)
4894  {
4895  case 3:
4896 
4897  if (evaluate_grad == true)
4898  {
4899  // grad x
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>
4903  (temp1, temp2);
4904  fe_eval.template apply_values<2,true,false>
4905  (temp2, fe_eval.gradients_quad[c][0]);
4906  }
4907 
4908  if (evaluate_lapl == true)
4909  {
4910  // grad xz
4911  if (evaluate_grad == false)
4912  {
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>
4916  (temp1, temp2);
4917  }
4918  fe_eval.template apply_gradients<2,true,false>
4919  (temp2, fe_eval.hessians_quad[c][4]);
4920 
4921  // grad xy
4922  fe_eval.template apply_gradients<1,true,false>
4923  (temp1, temp2);
4924  fe_eval.template apply_values<2,true,false>
4925  (temp2, fe_eval.hessians_quad[c][3]);
4926 
4927  // grad xx
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>
4931  (temp1, temp2);
4932  fe_eval.template apply_values<2,true,false>
4933  (temp2, fe_eval.hessians_quad[c][0]);
4934  }
4935 
4936  // grad y
4937  fe_eval.template apply_values<0,true,false>
4938  (fe_eval.values_dofs[c], temp1);
4939  if (evaluate_grad == true)
4940  {
4941  fe_eval.template apply_gradients<1,true,false>
4942  (temp1, temp2);
4943  fe_eval.template apply_values<2,true,false>
4944  (temp2, fe_eval.gradients_quad[c][1]);
4945  }
4946 
4947  if (evaluate_lapl == true)
4948  {
4949  // grad yz
4950  if (evaluate_grad == false)
4951  fe_eval.template apply_gradients<1,true,false>
4952  (temp1, temp2);
4953  fe_eval.template apply_gradients<2,true,false>
4954  (temp2, fe_eval.hessians_quad[c][5]);
4955 
4956  // grad yy
4957  fe_eval.template apply_hessians<1,true,false>
4958  (temp1, temp2);
4959  fe_eval.template apply_values<2,true,false>
4960  (temp2, fe_eval.hessians_quad[c][1]);
4961  }
4962 
4963  // grad z: can use the values applied in x direction stored in temp1
4964  fe_eval.template apply_values<1,true,false>
4965  (temp1, temp2);
4966  if (evaluate_grad == true)
4967  fe_eval.template apply_gradients<2,true,false>
4968  (temp2, fe_eval.gradients_quad[c][2]);
4969 
4970  // grad zz: can use the values applied in x and y direction stored
4971  // in temp2
4972  if (evaluate_lapl == true)
4973  fe_eval.template apply_hessians<2,true,false>
4974  (temp2, fe_eval.hessians_quad[c][2]);
4975 
4976  // val: can use the values applied in x & y direction stored in temp2
4977  if (evaluate_val == true)
4978  fe_eval.template apply_values<2,true,false>
4979  (temp2, fe_eval.values_quad[c]);
4980 
4981  break;
4982 
4983  case 2:
4984 
4985  // grad x
4986  if (evaluate_grad == true)
4987  {
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]);
4992  }
4993  if (evaluate_lapl == true)
4994  {
4995  // grad xy
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]);
5001 
5002  // grad xx
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]);
5007  }
5008 
5009  // grad y
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]);
5015 
5016  // grad yy
5017  if (evaluate_lapl == true)
5018  fe_eval.template apply_hessians<1,true,false>
5019  (temp1, fe_eval.hessians_quad[c][1]);
5020 
5021  // val: can use values applied in x
5022  if (evaluate_val == true)
5023  fe_eval.template apply_values<1,true,false>
5024  (temp1, fe_eval.values_quad[c]);
5025 
5026  break;
5027 
5028  case 1:
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]);
5038  break;
5039 
5040  default:
5041  Assert (false, ExcNotImplemented());
5042  }
5043  }
5044 
5045 #ifdef DEBUG
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;
5052 #endif
5053  }
5054 
5055 
5056 
5057  template <typename FEEval>
5058  inline
5059  void
5060  do_integrate (FEEval &fe_eval,
5061  const bool integrate_val,
5062  const bool integrate_grad)
5063  {
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());
5071 
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;
5076 
5077 
5078  for (unsigned int c=0; c<n_components; c++)
5079  {
5082 
5083  switch (dim)
5084  {
5085  case 3:
5086 
5087  if (integrate_val == true)
5088  {
5089  // val
5090  fe_eval.template apply_values<0,false,false>
5091  (fe_eval.values_quad[c], temp1);
5092  }
5093  if (integrate_grad == true)
5094  {
5095  // grad x: can sum to temporary value in temp1
5096  if (integrate_val == true)
5097  fe_eval.template apply_gradients<0,false,true>
5098  (fe_eval.gradients_quad[c][0], temp1);
5099  else
5100  fe_eval.template apply_gradients<0,false,false>
5101  (fe_eval.gradients_quad[c][0], temp1);
5102  }
5103  fe_eval.template apply_values<1,false,false>
5104  (temp1, temp2);
5105  if (integrate_grad == true)
5106  {
5107  // grad y: can sum to temporary x value in temp2
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>
5111  (temp1, temp2);
5112  }
5113  fe_eval.template apply_values<2,false,false>
5114  (temp2, fe_eval.values_dofs[c]);
5115  if (integrate_grad == true)
5116  {
5117  // grad z: can sum to temporary x and y value in output
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>
5121  (temp1, temp2);
5122  fe_eval.template apply_gradients<2,false,true>
5123  (temp2, fe_eval.values_dofs[c]);
5124  }
5125 
5126  break;
5127 
5128  case 2:
5129 
5130  // val
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)
5135  {
5136  //grad x
5137  if (integrate_val == true)
5138  fe_eval.template apply_gradients<0,false,true>
5139  (fe_eval.gradients_quad[c][0], temp1);
5140  else
5141  fe_eval.template apply_gradients<0,false,false>
5142  (fe_eval.gradients_quad[c][0], temp1);
5143  }
5144  fe_eval.template apply_values<1,false,false>
5145  (temp1, fe_eval.values_dofs[c]);
5146  if (integrate_grad == true)
5147  {
5148  // grad y
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]);
5153  }
5154 
5155  break;
5156 
5157  case 1:
5158 
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)
5163  {
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]);
5167  else
5168  fe_eval.template apply_values<0,false,false>
5169  (fe_eval.values_quad[c], fe_eval.values_dofs[c]);
5170  }
5171  break;
5172 
5173  default:
5174  Assert (false, ExcNotImplemented());
5175  }
5176  }
5177 
5178 #ifdef DEBUG
5179  fe_eval.dof_values_initialized = true;
5180 #endif
5181  }
5182 
5183 } // end of namespace internal
5184 
5185 
5186 
5187 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
5188  typename Number>
5189 inline
5190 void
5192 ::evaluate (const bool evaluate_val,
5193  const bool evaluate_grad,
5194  const bool evaluate_lapl)
5195 {
5196  internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl);
5197 }
5198 
5199 
5200 
5201 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
5202  typename Number>
5203 inline
5204 void
5206 ::integrate (const bool integrate_val,
5207  const bool integrate_grad)
5208 {
5209  internal::do_integrate (*this, integrate_val, integrate_grad);
5210 }
5211 
5212 
5213 
5214 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
5215  typename Number>
5216 inline
5219 ::quadrature_point (const unsigned int q) const
5220 {
5221  Assert (this->mapping_info.quadrature_points_initialized == true,
5222  ExcNotInitialized());
5223  AssertIndexRange (q, n_q_points);
5224 
5225  // Cartesian mesh: not all quadrature points are stored, only the
5226  // diagonal. Hence, need to find the tensor product index and retrieve the
5227  // value from that
5228  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5229  {
5230  Point<dim,VectorizedArray<Number> > point (false);
5231  switch (dim)
5232  {
5233  case 1:
5234  return this->quadrature_points[q];
5235  case 2:
5236  point[0] = this->quadrature_points[q%n_q_points_1d][0];
5237  point[1] = this->quadrature_points[q/n_q_points_1d][1];
5238  return point;
5239  case 3:
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];
5243  return point;
5244  default:
5245  Assert (false, ExcNotImplemented());
5246  return point;
5247  }
5248  }
5249  // all other cases: just return the respective data as it is fully stored
5250  else
5251  return this->quadrature_points[q];
5252 }
5253 
5254 
5255 
5256 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
5257  typename Number>
5258 template <int direction, bool dof_to_quad, bool add>
5259 inline
5260 void
5263  VectorizedArray<Number> out [])
5264 {
5265  internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,
5266  VectorizedArray<Number>, direction, dof_to_quad, add>
5267  (this->data.shape_values.begin(), in, out);
5268 }
5269 
5270 
5271 
5272 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
5273  typename Number>
5274 template <int direction, bool dof_to_quad, bool add>
5275 inline
5276 void
5278 ::apply_gradients(const VectorizedArray<Number> in [],
5279  VectorizedArray<Number> out [])
5280 {
5281  internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,
5282  VectorizedArray<Number>, direction, dof_to_quad, add>
5283  (this->data.shape_gradients.begin(), in, out);
5284 }
5285 
5286 
5287 
5288 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
5289  typename Number>
5290 template <int direction, bool dof_to_quad, bool add>
5291 inline
5292 void
5294 ::apply_hessians(const VectorizedArray<Number> in [],
5295  VectorizedArray<Number> out [])
5296 {
5297  internal::apply_tensor_product<dim,fe_degree,n_q_points_1d,
5298  VectorizedArray<Number>, direction, dof_to_quad, add>
5299  (this->data.shape_hessians.begin(), in, out);
5300 }
5301 
5302 
5303 /*-------------------------- FEEvaluation -----------------------------------*/
5304 
5305 
5306 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
5307  typename Number>
5308 inline
5310 ::FEEvaluation (const MatrixFree<dim,Number> &data_in,
5311  const unsigned int fe_no,
5312  const unsigned int quad_no)
5313  :
5314  BaseClass (data_in, fe_no, quad_no)
5315 {
5316  // check whether element is appropriate
5317 #ifdef DEBUG
5318  const double zero_tol =
5319  types_are_equal<Number,double>::value==true?1e-8:1e-7;
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!";
5324 
5325  // symmetry for values
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)
5329  Assert (std::fabs(this->data.shape_values[i*n_q_points_1d+j][0] -
5330  this->data.shape_values[(n_dofs_1d-i)*n_q_points_1d
5331  -j-1][0]) < zero_tol,
5332  ExcMessage(error_message));
5333 
5334  // shape values should be zero at for all basis functions except for one
5335  // where they are one in the middle
5336  if (n_q_points_1d%2 == 1 && n_dofs_1d%2 == 1)
5337  {
5338  for (int i=0; i<static_cast<int>(n_dofs_1d/2); ++i)
5339  Assert (std::fabs(this->data.shape_values[i*n_q_points_1d+
5340  n_q_points_1d/2][0]) < zero_tol,
5341  ExcMessage(error_message));
5342  Assert (std::fabs(this->data.shape_values[(n_dofs_1d/2)*n_q_points_1d+
5343  n_q_points_1d/2][0]-1.)< zero_tol,
5344  ExcMessage(error_message));
5345  }
5346 
5347  // skew-symmetry for gradient, zero of middle basis function in middle
5348  // quadrature point
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)
5351  Assert (std::fabs(this->data.shape_gradients[i*n_q_points_1d+j][0] +
5352  this->data.shape_gradients[(n_dofs_1d-i)*n_q_points_1d-
5353  j-1][0]) < zero_tol,
5354  ExcMessage(error_message));
5355  if (n_dofs_1d%2 == 1 && n_q_points_1d%2 == 1)
5356  Assert (std::fabs(this->data.shape_gradients[(n_dofs_1d/2)*n_q_points_1d+
5357  (n_q_points_1d/2)][0]) < zero_tol,
5358  ExcMessage(error_message));
5359 
5360 
5361  // symmetry for Laplacian
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)
5364  Assert (std::fabs(this->data.shape_hessians[i*n_q_points_1d+j][0] -
5365  this->data.shape_hessians[(n_dofs_1d-i)*n_q_points_1d-
5366  j-1][0]) < zero_tol,
5367  ExcMessage(error_message));
5368 #endif
5369 
5370  // Compute symmetric and skew-symmetric part of shape values for even-odd
5371  // decomposition
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)
5374  {
5375  shape_val_evenodd[i][q] =
5376  0.5 * (this->data.shape_values[i*n_q_points_1d+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] =
5379  0.5 * (this->data.shape_values[i*n_q_points_1d+q] -
5380  this->data.shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
5381 
5382  shape_gra_evenodd[i][q] =
5383  0.5 * (this->data.shape_gradients[i*n_q_points_1d+q] +
5384  this->data.shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
5385  shape_gra_evenodd[fe_degree-i][q] =
5386  0.5 * (this->data.shape_gradients[i*n_q_points_1d+q] -
5387  this->data.shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
5388 
5389  shape_hes_evenodd[i][q] =
5390  0.5 * (this->data.shape_hessians[i*n_q_points_1d+q] +
5391  this->data.shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
5392  shape_hes_evenodd[fe_degree-i][q] =
5393  0.5 * (this->data.shape_hessians[i*n_q_points_1d+q] -
5394  this->data.shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
5395  }
5396  if (fe_degree % 2 == 0)
5397  for (unsigned int q=0; q<(n_q_points_1d+1)/2; ++q)
5398  {
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] =
5402  this->data.shape_gradients[(fe_degree/2)*n_q_points_1d+q];
5403  shape_hes_evenodd[fe_degree/2][q] =
5404  this->data.shape_hessians[(fe_degree/2)*n_q_points_1d+q];
5405  }
5406 }
5407 
5408 
5409 
5410 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
5411  typename Number>
5412 inline
5413 void
5415 ::evaluate (const bool evaluate_val,
5416  const bool evaluate_grad,
5417  const bool evaluate_lapl)
5418 {
5419  internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl);
5420 }
5421 
5422 
5423 
5424 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
5425  typename Number>
5426 inline
5427 void
5429 ::integrate (bool integrate_val,bool integrate_grad)
5430 {
5431  internal::do_integrate (*this, integrate_val, integrate_grad);
5432 }
5433 
5434 
5435 
5436 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
5437  typename Number>
5438 template <int direction, bool dof_to_quad, bool add>
5439 inline
5440 void
5442 ::apply_values (const VectorizedArray<Number> in [],
5443  VectorizedArray<Number> out [])
5444 {
5445  // for linear elements, the even-odd decomposition is slower than the plain
5446  // evaluation (additions to create the symmetric and anti-symmetric part),
5447  // for all other orders, we choose even-odd
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);
5452  else
5453  internal::apply_tensor_product_values<dim,fe_degree,n_q_points_1d,
5454  VectorizedArray<Number>, direction, dof_to_quad, add>
5455  (this->data.shape_values.begin(), in, out);
5456 }
5457 
5458 
5459 
5460 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
5461  typename Number>
5462 template <int direction, bool dof_to_quad, bool add>
5463 inline
5464 void
5466 ::apply_gradients (const VectorizedArray<Number> in [],
5467  VectorizedArray<Number> out [])
5468 {
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);
5473  else
5474  internal::apply_tensor_product_gradients<dim,fe_degree,n_q_points_1d,
5475  VectorizedArray<Number>, direction, dof_to_quad, add>
5476  (this->data.shape_gradients.begin(), in, out);
5477 }
5478 
5479 
5480 
5481 // Laplacian operator application. Very similar to value application because
5482 // the same symmetry relations hold. However, it is not possible to omit some
5483 // values that are zero for the values
5484 template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
5485  typename Number>
5486 template <int direction, bool dof_to_quad, bool add>
5487 inline
5488 void
5490 ::apply_hessians (const VectorizedArray<Number> in [],
5491  VectorizedArray<Number> out [])
5492 {
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);
5497  else
5498  internal::apply_tensor_product_hessians<dim,fe_degree,n_q_points_1d,
5499  VectorizedArray<Number>, direction, dof_to_quad, add>
5500  (this->data.shape_hessians.begin(), in, out);
5501 }
5502 
5503 
5504 
5505 /*------------------------- FEEvaluationGL ----------------------------------*/
5506 
5507 template <int dim, int fe_degree, int n_components_, typename Number>
5508 inline
5511  const unsigned int fe_no,
5512  const unsigned int quad_no)
5513  :
5514  BaseClass (data_in, fe_no, quad_no)
5515 {
5516 #ifdef DEBUG
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!";
5522 
5523  const double zero_tol =
5524  types_are_equal<Number,double>::value==true?1e-9:1e-7;
5525 
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)
5529  if (i!=j)
5530  {
5531  Assert (std::fabs(this->data.shape_values[i*n_points_1d+j][0])<zero_tol,
5532  ExcMessage (error_mess.c_str()));
5533  }
5534  else
5535  {
5536  Assert (std::fabs(this->data.shape_values[i*n_points_1d+
5537  j][0]-1.)<zero_tol,
5538  ExcMessage (error_mess.c_str()));
5539  }
5540  for (unsigned int i=1; i<n_points_1d-1; ++i)
5541  Assert (std::fabs(this->data.shape_gradients[i*n_points_1d+i][0])<zero_tol,
5542  ExcMessage (error_mess.c_str()));
5543  Assert (std::fabs(this->data.shape_gradients[n_points_1d-1][0]-
5544  (n_points_1d%2==0 ? -1. : 1.)) < zero_tol,
5545  ExcMessage (error_mess.c_str()));
5546 #endif
5547 }
5548 
5549 
5550 
5551 template <int dim, int fe_degree, int n_components_, typename Number>
5552 inline
5553 void
5555 ::evaluate (const bool evaluate_val,
5556  const bool evaluate_grad,
5557  const bool evaluate_lapl)
5558 {
5559  Assert (this->cell != numbers::invalid_unsigned_int,
5560  ExcNotInitialized());
5561  Assert (this->dof_values_initialized == true,
5562  internal::ExcAccessToUninitializedField());
5563 
5564  if (evaluate_val == true)
5565  {
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]));
5569 #ifdef DEBUG
5570  this->values_quad_initialized = true;
5571 #endif
5572  }
5573  // separate implementation here compared to the general case because the
5574  // values are an identity operation
5575  if (evaluate_grad == true)
5576  {
5577  for (unsigned int comp=0; comp<n_components; comp++)
5578  {
5579  if (dim == 3)
5580  {
5581  // grad x
5582  apply_gradients<0,true,false> (this->values_dofs[comp],
5583  this->gradients_quad[comp][0]);
5584  // grad y
5585  apply_gradients<1,true,false> (this->values_dofs[comp],
5586  this->gradients_quad[comp][1]);
5587  // grad y
5588  apply_gradients<2,true,false> (this->values_dofs[comp],
5589  this->gradients_quad[comp][2]);
5590  }
5591  else if (dim == 2)
5592  {
5593  // grad x
5594  apply_gradients<0,true,false> (this->values_dofs[comp],
5595  this->gradients_quad[comp][0]);
5596  // grad y
5597  apply_gradients<1,true,false> (this->values_dofs[comp],
5598  this->gradients_quad[comp][1]);
5599  }
5600  else if (dim == 1)
5601  apply_gradients<0,true,false> (this->values_dofs[comp],
5602  this->gradients_quad[comp][0]);
5603  }
5604 #ifdef DEBUG
5605  this->gradients_quad_initialized = true;
5606 #endif
5607  }
5608  if (evaluate_lapl == true)
5609  {
5610  for (unsigned int comp=0; comp<n_components; comp++)
5611  {
5612  if (dim == 3)
5613  {
5614  // grad x
5615  this->template apply_hessians<0,true,false> (this->values_dofs[comp],
5616  this->hessians_quad[comp][0]);
5617  // grad y
5618  this->template apply_hessians<1,true,false> (this->values_dofs[comp],
5619  this->hessians_quad[comp][1]);
5620  // grad y
5621  this->template apply_hessians<2,true,false> (this->values_dofs[comp],
5622  this->hessians_quad[comp][2]);
5623 
5624  VectorizedArray<Number> temp1[n_q_points];
5625  // grad xy
5626  apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
5627  apply_gradients<1,true,false> (temp1, this->hessians_quad[comp][3]);
5628  // grad xz
5629  apply_gradients<2,true,false> (temp1, this->hessians_quad[comp][4]);
5630  // grad yz
5631  apply_gradients<1,true,false> (this->values_dofs[comp], temp1);
5632  apply_gradients<2,true,false> (temp1, this->hessians_quad[comp][5]);
5633  }
5634  else if (dim == 2)
5635  {
5636  // grad x
5637  this->template apply_hessians<0,true,false> (this->values_dofs[comp],
5638  this->hessians_quad[comp][0]);
5639  // grad y
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];
5643  // grad xy
5644  apply_gradients<0,true,false> (this->values_dofs[comp], temp1);
5645  apply_gradients<1,true,false> (temp1, this->hessians_quad[comp][2]);
5646  }
5647  else if (dim == 1)
5648  this->template apply_hessians<0,true,false> (this->values_dofs[comp],
5649  this->hessians_quad[comp][0]);
5650  }
5651 #ifdef DEBUG
5652  this->hessians_quad_initialized = true;
5653 #endif
5654  }
5655 }
5656 
5657 
5658 
5659 template <int dim, int fe_degree, int n_components_, typename Number>
5660 inline
5661 void
5663 ::integrate (const bool integrate_val, const bool integrate_grad)
5664 {
5665  Assert (this->cell != numbers::invalid_unsigned_int,
5666  ExcNotInitialized());
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)
5678  {
5679  for (unsigned int comp=0; comp<n_components; comp++)
5680  {
5681  if (dim == 3)
5682  {
5683  // grad x: If integrate_val == true we have to add to the previous output
5684  if (integrate_val == true)
5685  apply_gradients<0, false, true> (this->gradients_quad[comp][0],
5686  this->values_dofs[comp]);
5687  else
5688  apply_gradients<0, false, false> (this->gradients_quad[comp][0],
5689  this->values_dofs[comp]);
5690 
5691  // grad y: can sum to temporary x value in temp2
5692  apply_gradients<1, false, true> (this->gradients_quad[comp][1],
5693  this->values_dofs[comp]);
5694 
5695  // grad z: can sum to temporary x and y value in output
5696  apply_gradients<2, false, true> (this->gradients_quad[comp][2],
5697  this->values_dofs[comp]);
5698  }
5699  else if (dim == 2)
5700  {
5701  // grad x: If integrate_val == true we have to add to the previous output
5702  if (integrate_val == true)
5703  apply_gradients<0, false, true> (this->gradients_quad[comp][0],
5704  this->values_dofs[comp]);
5705  else
5706  apply_gradients<0, false, false> (this->gradients_quad[comp][0],
5707  this->values_dofs[comp]);
5708 
5709  // grad y: can sum to temporary x value in temp2
5710  apply_gradients<1, false, true> (this->gradients_quad[comp][1],
5711  this->values_dofs[comp]);
5712  }
5713  else if (dim == 1)
5714  {
5715  if (integrate_val == true)
5716  apply_gradients<0, false, true> (this->gradients_quad[comp][0],
5717  this->values_dofs[comp]);
5718  else
5719  apply_gradients<0, false, false> (this->gradients_quad[comp][0],
5720  this->values_dofs[comp]);
5721 
5722  }
5723  }
5724  }
5725 
5726 #ifdef DEBUG
5727  this->dof_values_initialized = true;
5728 #endif
5729 }
5730 
5731 
5732 
5733 template <int dim, int fe_degree, int n_components_, typename Number>
5734 template <int direction, bool dof_to_quad, bool add>
5735 inline
5736 void
5738 ::apply_gradients (const VectorizedArray<Number> in [],
5739  VectorizedArray<Number> out [])
5740 {
5741  internal::apply_tensor_product_gradients_gl<dim,fe_degree,
5742  VectorizedArray<Number>, direction, dof_to_quad, add>
5743  (this->data.shape_gradients.begin(), in, out);
5744 }
5745 
5746 
5747 #endif // ifndef DOXYGEN
5748 
5749 
5750 DEAL_II_NAMESPACE_CLOSE
5751 
5752 #endif
const Tensor< 2, dim, VectorizedArray< Number > > * jacobian_grad
const unsigned int quad_no
unsigned int cell
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
Definition: types.h:191
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
AlignedVector< std::pair< Tensor< 1, dim, VectorizedArray< Number > >, VectorizedArray< Number > > > cartesian_data
Definition: mapping_info.h:139
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
Definition: mapping_info.h:154
bool mapping_initialized() const
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)
Definition: exceptions.h:888
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
Definition: mapping_info.h:284
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)
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
Definition: shape_info.h:107
void apply_hessians(const VectorizedArray< Number > in[], VectorizedArray< Number > out[])
CellType get_cell_type(const unsigned int cell_chunk_no) const
Definition: mapping_info.h:350
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)
Definition: exceptions.h:299
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)
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)
Definition: exceptions.h:505
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)
void evaluate(const bool evaluate_val, const bool evaluate_grad, const bool evaluate_hess=false)
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
Definition: mapping_info.h:363
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)
iterator begin()
void read_dof_values_plain(const VectorType &src)
Definition: types.h:173
Number local_element(const size_type local_index) const
bool partitioners_are_compatible(const Utilities::MPI::Partitioner &part) const
::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
Definition: shape_info.h:81
const VectorizedArray< Number > * begin_hessians() const
const Number * constraint_pool_begin(const unsigned int pool_index) const
AlignedVector< VectorizedArray< Number > > shape_gradients
Definition: shape_info.h:94
bool values_quad_initialized
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
void read_write_operation(const VectorOperation &operation, VectorType *vectors[]) const
void apply_gradients(const VectorizedArray< Number > in[], VectorizedArray< Number > out[])
::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)