Flow123d  DF_patch_fe_data_tables-3df4116
patch_point_values.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file patch_point_values.hh
15  * @brief Store finite element data on the actual patch
16  * such as shape function values, gradients, Jacobian
17  * of the mapping from the reference cell etc.
18  * @author David Flanderka
19  */
20 
21 #ifndef PATCH_POINT_VALUES_HH_
22 #define PATCH_POINT_VALUES_HH_
23 
24 #include <Eigen/Dense>
25 
26 #include "fem/eigen_tools.hh"
27 #include "fem/dh_cell_accessor.hh"
28 #include "fem/element_values.hh"
30 #include "fem/arena_resource.hh"
31 #include "fem/arena_vec.hh"
32 
33 
34 template<unsigned int spacedim> class PatchFEValues;
35 template<unsigned int spacedim> class ElOp;
36 template<unsigned int dim> class BulkValues;
37 template<unsigned int dim> class SideValues;
38 using Scalar = double;
41 
42 
43 
44 
45 /// Type for conciseness
46 using ReinitFunction = std::function<void(std::vector<ElOp<3>> &, TableDbl &, TableInt &)>;
47 
48 
49 namespace FeBulk {
50  /**
51  * Enumeration of element bulk operations
52  *
53  * Operations are stored in fix order. Order in enum is equal to order
54  * in PatchPointVale::operations_ vector. FE operations are added dynamically
55  * by request of user.
56  */
57  enum BulkOps
58  {
59  /// operations evaluated on elements
60  opElCoords, ///< coordinations of all nodes of element
61  opJac, ///< Jacobian of element
62  opInvJac, ///< inverse Jacobian
63  opJacDet, ///< determinant of Jacobian
64  /// operations evaluated on quadrature points
65  opCoords, ///< coordinations of quadrature point
66  opWeights, ///< weight of quadrature point
67  opJxW ///< JxW value of quadrature point
68  };
69 }
70 
71 
72 namespace FeSide {
73  /**
74  * Enumeration of element side operations
75  *
76  * Operations are stored in fix order. Order in enum is equal to order
77  * in PatchPointVale::operations_ vector. FE operations are added dynamically
78  * by request of user.
79  */
80  enum SideOps
81  {
82  /// operations evaluated on elements
83  opElCoords, ///< coordinations of all nodes of element
84  opElJac, ///< Jacobian of element
85  opElInvJac, ///< inverse Jacobian of element
86  /// operations evaluated on sides
87  opSideCoords, ///< coordinations of all nodes of side
88  opSideJac, ///< Jacobian of element
89  opSideJacDet, ///< determinant of Jacobian of side
90  /// operation executed expansion to quadrature point (value of element / side > values on quadrature points)
91 // opExpansionElCoords, ///< expands coordinates on element
92 // opExpansionElJac, ///< expands Jacobian on element
93 // opExpansionElInvJac, ///< expands inverse Jacobian on element
94 // opExpansionSideCoords, ///< expands coordinates on side
95 // opExpansionSideJac, ///< expands Jacobian on side
96 // opExpansionSideJacDet, ///< expands Jacobian determinant on side
97  /// operations evaluated on quadrature points
98  opCoords, ///< coordinations of quadrature point
99  opWeights, ///< weight of quadrature point
100  opJxW, ///< JxW value of quadrature point
101  opNormalVec ///< normal vector of quadrature point
102  };
103 }
104 
105 
106 /// Distinguishes operations by type and size of output rows
108 {
109  elemOp, ///< operation is evaluated on elements or sides
110  pointOp, ///< operation is evaluated on quadrature points
111  fixedSizeOp ///< operation has fixed size and it is filled during initialization
112 };
113 
114 
115 
116 /**
117  * v Class for storing FE data of quadrature points on one patch.
118  *
119  * Store data of bulk or side quadrature points of one dimension.
120  */
121 template<unsigned int spacedim = 3>
123 {
124 public:
125  /**
126  * Constructor
127  *
128  * @param dim Set dimension
129  */
131  : dim_(dim), n_rows_(0), elements_map_(300, 0), points_map_(300, 0), asm_arena_(asm_arena) {
132  reset();
133  }
134 
135  /**
136  * Initialize object, set number of columns (quantities) in tables.
137  *
138  * Number of columns of int_vals_ table is passed by argument \p int_cols, number of columns
139  * of other tables is given by n_rows_ value.
140  */
141  void initialize(uint int_cols) {
142  this->reset();
143 
144  point_vals_.resize(n_rows_);
145  int_vals_.resize(int_cols);
146  }
147 
148  inline void init_finalize(PatchArena *patch_arena) {
149  patch_arena_ = patch_arena;
150  }
151 
152  /// Reset number of columns (points and elements)
153  inline void reset() {
154  n_points_ = 0;
155  n_elems_ = 0;
156  i_elem_ = 0;
157  }
158 
159  /// Getter for dim_
160  inline uint dim() const {
161  return dim_;
162  }
163 
164  /// Getter for n_rows_
165  inline uint n_rows() const {
166  return n_rows_;
167  }
168 
169  /// Getter for n_elems_
170  inline uint n_elems() const {
171  return n_elems_;
172  }
173 
174  /// Getter for n_points_
175  inline uint n_points() const {
176  return n_points_;
177  }
178 
179  /// Getter for quadrature
181  return quad_;
182  }
183 
184  /// Resize data tables. Method is called before reinit of patch.
186  n_elems_ = n_elems;
187  // n_points_ = n_points; // will be set here later
189  for (auto &elOp : operations_)
190  if (elOp.size_type() != fixedSizeOp) {
191  elOp.allocate_result(sizes[elOp.size_type()], *patch_arena_);
192  }
193  for (uint i=0; i<point_vals_.rows(); ++i) {
194  if (row_sizes_[i] == elemOp) {
195  point_vals_(i).resize(n_elems);
196  point_vals_(i).setZero(n_elems,1);
197  } else if (row_sizes_[i] == pointOp) {
198  point_vals_(i).resize(n_points);
199  point_vals_(i).setZero(n_points,1);
200  }
201  }
203  }
204 
205  /**
206  * Register element, add to point_vals_ table
207  *
208  * @param coords Coordinates of element nodes.
209  * @param element_patch_idx Index of element on patch.
210  */
211  uint register_element(arma::mat coords, uint element_patch_idx) {
213  uint res_column = op.result_row();
214  auto &coords_mat = op.result_matrix();
215  std::size_t i_elem = i_elem_;
216  for (uint i_col=0; i_col<coords.n_cols; ++i_col)
217  for (uint i_row=0; i_row<coords.n_rows; ++i_row) {
218  point_vals_(res_column)(i_elem_) = coords(i_row, i_col);
219  coords_mat(i_row, i_col)(i_elem) = coords(i_row, i_col);
220  ++res_column;
221  }
222 
223  elements_map_[element_patch_idx] = i_elem_;
224  return i_elem_++;
225  }
226 
227  /**
228  * Register side, add to point_vals_ table
229  *
230  * @param coords Coordinates of element nodes.
231  * @param side_coords Coordinates of side nodes.
232  */
233  uint register_side(arma::mat elm_coords, arma::mat side_coords) {
234  {
236  uint res_column = op.result_row();
237  auto &coords_mat = op.result_matrix();
238  std::size_t i_elem = i_elem_;
239  for (uint i_col=0; i_col<elm_coords.n_cols; ++i_col)
240  for (uint i_row=0; i_row<elm_coords.n_rows; ++i_row) {
241  point_vals_(res_column)(i_elem_) = elm_coords(i_row, i_col);
242  coords_mat(i_row, i_col)(i_elem) = elm_coords(i_row, i_col);
243  ++res_column;
244  }
245  }
246 
247  {
249  uint res_column = op.result_row();
250  auto &coords_mat = op.result_matrix();
251  std::size_t i_elem = i_elem_;
252  for (uint i_col=0; i_col<side_coords.n_cols; ++i_col)
253  for (uint i_row=0; i_row<side_coords.n_rows; ++i_row) {
254  point_vals_(res_column)(i_elem_) = side_coords(i_row, i_col);
255  coords_mat(i_row, i_col)(i_elem) = side_coords(i_row, i_col);
256  ++res_column;
257  }
258  }
259 
260  return i_elem_++;
261  }
262 
263  /**
264  * Register bulk point, add to int_vals_ table
265  *
266  * @param elem_table_row Index of element in temporary element table.
267  * @param value_patch_idx Index of point in ElementCacheMap.
268  * @param elem_idx Index of element in Mesh.
269  * @param i_point_on_elem Index of point on element
270  */
271  uint register_bulk_point(uint elem_table_row, uint value_patch_idx, uint elem_idx, uint i_point_on_elem) {
272  uint point_pos = i_point_on_elem * n_elems_ + elem_table_row; // index of bulk point on patch
273  int_vals_(0)(point_pos) = value_patch_idx;
274  int_vals_(1)(point_pos) = elem_table_row;
275  int_vals_(2)(point_pos) = elem_idx;
276 
277  points_map_[value_patch_idx] = point_pos;
278  n_points_++;
279  return point_pos;
280  }
281 
282  /**
283  * Register side point, add to int_vals_ table
284  *
285  * @param elem_table_row Index of side in temporary element table.
286  * @param value_patch_idx Index of point in ElementCacheMap.
287  * @param elem_idx Index of element in Mesh.
288  * @param side_idx Index of side on element.
289  */
290  uint register_side_point(uint elem_table_row, uint value_patch_idx, uint elem_idx, uint side_idx) {
291  int_vals_(0)(n_points_) = value_patch_idx;
292  int_vals_(1)(n_points_) = elem_table_row;
293  int_vals_(2)(n_points_) = elem_idx;
294  int_vals_(3)(n_points_) = side_idx;
295 
296  points_map_[value_patch_idx] = n_points_;
297  return n_points_++;
298  }
299 
300  /**
301  * Adds accessor of new operation to operations_ vector
302  *
303  * @param shape Shape of function output
304  * @param reinit_f Reinitialize function
305  * @param input_ops_vec Indices of input operations in operations_ vector.
306  * @param size_type Type of operation by size of rows
307  */
308  ElOp<spacedim> &make_new_op(std::initializer_list<uint> shape, ReinitFunction reinit_f, std::vector<uint> input_ops_vec, OpSizeType size_type = pointOp) {
309  ElOp<spacedim> op_accessor(this->dim_, shape, this->n_rows_, reinit_f, size_type, input_ops_vec);
310  this->n_rows_ += op_accessor.n_comp();
311  row_sizes_.insert(row_sizes_.end(), op_accessor.n_comp(), size_type);
312  operations_.push_back(op_accessor);
313  return operations_[operations_.size()-1];
314  }
315 
316  /**
317  * Adds accessor of new operation with fixed data size (ref data) to operations_ vector
318  *
319  * @param shape Shape of function output
320  * @param reinit_f Reinitialize function
321  */
322  ElOp<spacedim> &make_fixed_op(std::initializer_list<uint> shape, ReinitFunction reinit_f) {
323  return make_new_op(shape, reinit_f, {}, fixedSizeOp);
324  }
325 
326  /**
327  * Adds accessor of expansion operation to operations_ vector
328  *
329  * @param el_op Source operation of expansion.
330  * @param shape Shape of function output
331  * @param reinit_f Reinitialize function
332  */
333  ElOp<spacedim> &make_expansion(ElOp<spacedim> &el_op, std::initializer_list<uint> shape, ReinitFunction reinit_f) {
334  ElOp<spacedim> op_accessor(this->dim_, shape, el_op.result_row(), reinit_f, OpSizeType::pointOp);
335  // shape passed from el_op throws:
336  // C++ exception with description "std::bad_alloc" thrown in the test body.
337  operations_.push_back(op_accessor);
338  return operations_[operations_.size()-1];
339  }
340 
341  /**
342  * Adds accessor of FE operation and adds operation dynamically to operations_ vector
343  *
344  * @param shape Shape of function output
345  * @param reinit_f Reinitialize function
346  * @param input_ops_vec Indices of input operations in operations_ vector.
347  * @param n_dofs Number of DOFs
348  * @param size_type Type of operation by size of rows
349  */
350  ElOp<spacedim> &make_fe_op(std::initializer_list<uint> shape, ReinitFunction reinit_f, std::vector<uint> input_ops_vec, uint n_dofs,
351  OpSizeType size_type = pointOp) {
352  ElOp<spacedim> op_accessor(this->dim_, shape, this->n_rows_, reinit_f, size_type, input_ops_vec, n_dofs);
353  this->n_rows_ += op_accessor.n_comp() * n_dofs;
354  row_sizes_.insert(row_sizes_.end(), op_accessor.n_comp() * n_dofs, size_type);
355  operations_.push_back(op_accessor);
356  return operations_[operations_.size()-1];
357  }
358 
359 
360  /**
361  * Reinitializes patch data.
362  *
363  * Calls reinit functions defined on each operations.
364  */
365  void reinit_patch() {
366  if (n_elems_ == 0) return; // skip if tables are empty
367  for (uint i=0; i<operations_.size(); ++i)
368  operations_[i].reinit_function(operations_, point_vals_, int_vals_);
369  }
370 
371  /**
372  * Returns scalar output value given by index of first row and index of quadrature point.
373  *
374  * @param result_row Row of operation in point_vals_ data table
375  * @param point_idx Index of quadrature point in ElementCacheMap
376  */
377  inline Scalar scalar_val(uint result_row, uint point_idx) const {
378  return point_vals_(result_row)(points_map_[point_idx]);
379  }
380  inline Scalar scalar_value(uint op_idx, uint point_idx, uint i_dof=0) const {
381  return operations_[op_idx].result_matrix()(0)(points_map_[point_idx] + i_dof*n_points_);
382  }
383 
384  /**
385  * Returns vector output value given by index of first row and index of quadrature point.
386  *
387  * @param result_row First row of operation in point_vals_ data table
388  * @param point_idx Index of quadrature point in ElementCacheMap
389  */
390  inline Vector vector_val(uint result_row, uint point_idx) const {
391  Vector val;
392  for (uint i=0; i<3; ++i)
393  val(i) = point_vals_(result_row+i)(points_map_[point_idx]);
394  return val;
395  }
396  inline Vector vector_value(uint op_idx, uint point_idx, uint i_dof=0) const {
397  Vector val;
398  const auto &op_matrix = operations_[op_idx].result_matrix();
399  uint op_matrix_idx = points_map_[point_idx] + i_dof*n_points_;
400  for (uint i=0; i<3; ++i)
401  val(i) = op_matrix(i)(op_matrix_idx);
402  return val;
403  }
404 
405  /**
406  * Returns tensor output value given by index of first row and index of quadrature point.
407  *
408  * @param result_row First row of operation in point_vals_ data table
409  * @param point_idx Index of quadrature point in ElementCacheMap
410  */
411  inline Tensor tensor_val(uint result_row, uint point_idx) const {
412  Tensor val;
413  for (uint i=0; i<3; ++i)
414  for (uint j=0; j<3; ++j)
415  val(i,j) = point_vals_(result_row+3*i+j)(points_map_[point_idx]);
416  return val;
417  }
418  inline Tensor tensor_value(uint op_idx, uint point_idx, uint i_dof=0) const {
419  Tensor val;
420  const auto &op_matrix = operations_[op_idx].result_matrix();
421  uint op_matrix_idx = points_map_[point_idx] + i_dof*n_points_;
422  for (uint i=0; i<3; ++i)
423  for (uint j=0; j<3; ++j)
424  val(i,j) = op_matrix(i,j)(op_matrix_idx);
425  return val;
426  }
427 
428  /**
429  * Performs output of data tables to stream.
430  *
431  * Development method.
432  * @param points Allows switched off output of point table,
433  * @param ints Allows switched off output of int (connectivity to elements) table,
434  */
435  void print_data_tables(ostream& stream, bool points, bool ints) const {
436  if (points) {
437  stream << "Point vals: " << point_vals_.rows() << " - " << point_vals_.cols() << std::endl;
438  for (auto &op : operations_) {
439  const auto &mat = op.result_matrix();
440  for (uint i_mat=0; i_mat<mat.rows()*mat.cols(); ++i_mat) {
441  if (mat(i_mat).data_size()==0) stream << "<empty>";
442  else {
443  const double *vals = mat(i_mat).data_ptr();
444  for (size_t i_val=0; i_val<mat(i_mat).data_size(); ++i_val)
445  stream << vals[i_val] << " ";
446  }
447  stream << std::endl;
448  }
449  }
450  stream << std::endl;
451  }
452  if (ints) {
453  stream << "Int vals: " << int_vals_.rows() << " - " << int_vals_.cols() << std::endl;
454  for (uint i_row=0; i_row<n_points_; ++i_row) {
455  for (uint i_col=0; i_col<3; ++i_col)
456  stream << int_vals_(i_col)(i_row) << " ";
457  stream << std::endl;
458  }
459  stream << std::endl;
460  }
461  }
462 
463  /**
464  * Performs table of fixed operations to stream.
465  *
466  * Development method.
467  * @param bulk_side Needs set 0 (bulk) or 1 (side) for correct output of operation names.
468  */
469  void print_operations(ostream& stream, uint bulk_side) const {
471  {
472  { "el_coords", "jacobian", "inv_jac", "jac_det", "pt_coords", "weights", "JxW", "", "", "", "", "" },
473  { "el_coords", "el_jac", "el_inv_jac", "side_coords", "side_jac", "side_jac_det", "exp_el_coords", "exp_el_jac", "exp_el_inv_jac",
474  "exp_side_coords", "exp_side_jac", "exp_side_jac_det", "pt_coords", "weights", "JxW", "normal_vec", "", "", "", "", "" }
475  };
476  stream << std::setfill(' ') << " Operation" << setw(12) << "" << "Shape" << setw(2) << ""
477  << "Result row" << setw(2) << "" << "n DOFs" << setw(2) << "" << "Input operations" << endl;
478  for (uint i=0; i<operations_.size(); ++i) {
479  stream << " " << std::left << setw(20) << op_names[bulk_side][i] << "" << " " << setw(6) << operations_[i].format_shape() << "" << " "
480  << setw(11) << operations_[i].result_row() << "" << " " << setw(7) << operations_[i].n_dofs() << "" << " ";
481  auto &input_ops = operations_[i].input_ops();
482  for (auto i_o : input_ops) stream << op_names[bulk_side][i_o] << " ";
483  stream << std::endl;
484  }
485  }
486 
487 protected:
488  /**
489  * Store data of bulk or side quadrature points of one dimension
490  *
491  * Number of columns is given by n_rows_, number of used rows by n_points_.
492  */
494  /**
495  * Hold integer values of quadrature points of previous table.
496  *
497  * Table contains following columns:
498  * 0: Index of quadrature point on patch
499  * 1: Row of element/side in point_vals_ table in registration step (before expansion)
500  * 2: Element idx in Mesh
501  * 3: Index of side in element (column is allocated only for side point table)
502  * Number of used rows is given by n_points_.
503  */
505 
506  /// Vector of all defined operations
508 
509  uint dim_; ///< Dimension
510  uint n_rows_; ///< Number of columns of \p point_vals table
511  uint n_points_; ///< Number of points in patch
512  uint n_elems_; ///< Number of elements in patch
513  uint i_elem_; ///< Index of registered element in table, helper value used during patch creating.
514  Quadrature *quad_; ///< Quadrature of given dimension and order passed in constructor.
515 
516  std::vector<uint> elements_map_; ///< Map of element patch indices to el_vals_ table
517  std::vector<uint> points_map_; ///< Map of point patch indices to point_vals_ and int_vals_ tables
518  std::vector<OpSizeType> row_sizes_; ///< hold sizes of rows by type of operation
519  AssemblyArena &asm_arena_; ///< Reference to global assembly arena of PatchFeValues
520  PatchArena *patch_arena_; ///< Pointer to global patch arena of PatchFeValues
521 
522  friend class PatchFEValues<spacedim>;
523  friend class ElOp<spacedim>;
524  template<unsigned int dim>
525  friend class BulkValues;
526  template<unsigned int dim>
527  friend class SideValues;
528  template<unsigned int dim>
529  friend class JoinValues;
530 };
531 
532 /**
533  * @brief Class represents element or FE operations.
534  */
535 template<unsigned int spacedim = 3>
536 class ElOp {
537 public:
538  /**
539  * Constructor
540  *
541  * Set all data members.
542  */
545  {}
546 
547  /// Aligns shape_vec to 2 items (equal to matrix number of dimensions)
548  std::vector<uint> set_shape_vec(std::initializer_list<uint> shape) const {
549  std::vector<uint> shape_vec(shape);
550  if (shape_vec.size() == 1) shape_vec.push_back(1);
551  ASSERT_EQ(shape_vec.size(), 2);
552  return shape_vec;
553  }
554 
555  /**
556  * Return number of operation components
557  *
558  * Value is computed from shape_ vector
559  */
560  inline uint n_comp() const {
561  return shape_[0] * shape_[1];
562  }
563 
564  /// Getter for dimension
565  inline uint dim() const {
566  return dim_;
567  }
568 
569  /// Getter for result_row_
570  inline uint result_row() const {
571  return result_row_;
572  }
573 
574  /// Getter for size_type_
576  return size_type_;
577  }
578 
579  /// Getter for n_dofs_
580  inline uint n_dofs() const {
581  return n_dofs_;
582  }
583 
584  /// Getter for input_ops_
585  inline const std::vector<uint> &input_ops() const {
586  return input_ops_;
587  }
588 
589  /// Getter for shape_
590  inline const std::vector<uint> &shape() const {
591  return shape_;
592  }
593 
594  /**
595  * Format shape to string
596  *
597  * Method is used in output development method.
598  */
599  inline std::string format_shape() const {
600  stringstream ss;
601  ss << shape_[0] << "x" << shape_[1];
602  return ss.str();
603  }
604 
605  /// Call reinit function on element table if function is defined
606  inline void reinit_function(std::vector<ElOp<spacedim>> &operations, TableDbl &data_table, TableInt &int_table) {
607  reinit_func(operations, data_table, int_table);
608  }
609 
610  inline void allocate_result(size_t data_size, PatchArena &arena) {
611  result_ = Eigen::Matrix<ArenaVec<double>, Eigen::Dynamic, Eigen::Dynamic>(shape_[0], shape_[1]);
612  for (uint i=0; i<n_comp(); ++i)
613  result_(i) = ArenaVec<double>(data_size, arena);
614  }
615 
616  /// Return map referenced result as Eigen::Matrix
617  Eigen::Matrix<ArenaVec<double>, Eigen::Dynamic, Eigen::Dynamic> &result_matrix() {
618  return result_;
619  }
620 // Eigen::Map<Eigen::Matrix<ArenaVec<double>, Eigen::Dynamic, Eigen::Dynamic>> result_matrix() {
621 // return Eigen::Map<Eigen::Matrix<ArenaVec<double>, Eigen::Dynamic, Eigen::Dynamic>>(result_.data(), shape_[0], shape_[1]);
622 // }
623 
624  /// Same as previous but return const reference
625  const Eigen::Matrix<ArenaVec<double>, Eigen::Dynamic, Eigen::Dynamic> &result_matrix() const {
626  return result_;
627  }
628 // const Eigen::Map<Eigen::Matrix<ArenaVec<double>, Eigen::Dynamic, Eigen::Dynamic>> result_matrix() const {
629 // return Eigen::Map<Eigen::Matrix<ArenaVec<double>, Eigen::Dynamic, Eigen::Dynamic>>(result_.data(), shape_[0], shape_[1]);
630 // }
631 
632  /// Return map referenced Eigen::Matrix of given dimension
633  template<unsigned int dim1, unsigned int dim2>
634  Eigen::Map<Eigen::Matrix<ArrayDbl, dim1, dim2>> value(TableDbl &op_results, uint i_dof = 0) const {
635  return Eigen::Map<Eigen::Matrix<ArrayDbl, dim1, dim2>>(op_results.data() + result_row_ + i_dof * n_comp(), dim1, dim2);
636  }
637 
638  /// Return map referenced Eigen::Matrix of given dimensions
639  Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> matrix_value(TableDbl &op_results, uint dim1, uint dim2) const {
640  return Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>>(op_results(result_row_).data(), dim1, dim2);
641  }
642 
643  /// Return map referenced Eigen::Matrix of given dimensions
644  Eigen::Map<Eigen::Vector<double, Eigen::Dynamic>> vector_value(TableDbl &op_results) const {
645  return Eigen::Map<Eigen::Vector<double, Eigen::Dynamic>>(op_results(result_row_).data(), op_results(result_row_).rows());
646  }
647 
648 
649 protected:
650  uint dim_; ///< Dimension
651  std::vector<uint> shape_; ///< Shape of stored data (size of vector or number of rows and cols of matrix)
652  uint result_row_; ///< First row to scalar, vector or matrix result TODO replace
653  Eigen::Matrix<ArenaVec<double>, Eigen::Dynamic, Eigen::Dynamic> result_;
654  ///< Result matrix of operation
655  OpSizeType size_type_; ///< Type of operation by size of vector (element, point or fixed size)
656  std::vector<uint> input_ops_; ///< Indices of operations in PatchPointValues::operations_ vector on which ElOp is depended
657  uint n_dofs_; ///< Number of DOFs of FE operations (or 1 in case of element operations)
658 
659  ReinitFunction reinit_func; ///< Pointer to patch reinit function of element data table specialized by operation
660 };
661 
662 
663 /// Defines common functionality of reinit operations.
665  // empty base operation
666  static inline void op_base(FMT_UNUSED std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
667  // empty
668  }
669 
670  // expansion operation
671  static inline void expand_data(ElOp<3> &op, TableDbl &op_results, TableInt &el_table) {
672  uint row_begin = op.result_row();
673  uint row_end = row_begin + op.n_comp();
674  uint size = op_results(row_begin).rows();
675  for (int i_pt=size-1; i_pt>=0; --i_pt) {
676  uint el_table_idx = el_table(1)(i_pt);
677  for (uint i_q=row_begin; i_q<row_end; ++i_q)
678  op_results(i_q)(i_pt) = op_results(i_q)(el_table_idx);
679  }
680  }
681 };
682 
683 /// Defines reinit operations on bulk points.
684 struct bulk_reinit {
685  // element operations
686  template<unsigned int dim>
687  static inline void elop_jac(std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
688  // result matrix(spacedim, dim), input matrix(spacedim, dim+1)
689  auto &op = operations[FeBulk::BulkOps::opJac];
690  auto &jac_value = op.result_matrix();
691  const auto &coords_value = operations[ op.input_ops()[0] ].result_matrix();
692  for (unsigned int i=0; i<3; i++)
693  for (unsigned int j=0; j<dim; j++)
694  jac_value(i,j) = coords_value(i,j+1) - coords_value(i,0);
695  }
696  template<unsigned int dim>
697  static inline void elop_inv_jac(std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
698  // result matrix(spacedim, dim), input matrix(spacedim, dim+1)
699  auto &op = operations[FeBulk::BulkOps::opInvJac];
700  auto &inv_jac_value = op.result_matrix();
701  const auto &jac_value = operations[ op.input_ops()[0] ].result_matrix();
702  inv_jac_value = eigen_arena_tools::inverse<3, dim>(jac_value);
703  }
704  template<unsigned int dim>
705  static inline void elop_jac_det(std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
706  // result double, input matrix(spacedim, dim)
707  auto &op = operations[FeBulk::BulkOps::opJacDet];
708  auto &jac_det_value = op.result_matrix();
709  const auto &jac_value = operations[ op.input_ops()[0] ].result_matrix();
710  jac_det_value(0,0) = eigen_arena_tools::determinant<3, dim>(jac_value).abs();
711  }
712 
713  // point operations
714  static inline void ptop_coords(FMT_UNUSED std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
715  // Implement
716  }
717  static inline void ptop_weights(std::vector<ElOp<3>> &operations, PatchArena *arena, const std::vector<double> &point_weights) {
718  auto &op = operations[FeBulk::BulkOps::opWeights];
719  op.allocate_result(point_weights.size(), *arena);
720  auto &weights_value = op.result_matrix();
721  for (uint i=0; i<point_weights.size(); ++i)
722  weights_value(0,0)(i) = point_weights[i];
723  }
724  static inline void ptop_JxW(std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
725  auto &op = operations[FeBulk::BulkOps::opJxW];
726  auto &weights_value = operations[ op.input_ops()[0] ].result_matrix();
727  auto &jac_det_value = operations[ op.input_ops()[1] ].result_matrix();
728  ArenaOVec<double> weights_ovec( weights_value(0,0), weights_value(0,0).data_size() );
729  ArenaOVec<double> jac_det_ovec( jac_det_value(0,0), jac_det_value(0,0).data_size() );
730  ArenaOVec<double> jxw_ovec = jac_det_ovec * weights_ovec;
731  auto &jxw_value = op.result_matrix();
732  jxw_value(0,0) = jxw_ovec.get_vec();
733  }
734  static inline void ptop_scalar_shape(std::vector<ElOp<3>> &operations,
735  std::vector< std::vector<double> > shape_values, uint scalar_shape_op_idx) {
736  auto &op = operations[scalar_shape_op_idx];
737  uint n_dofs = shape_values.size();
738  uint n_points = shape_values[0].size();
739  uint n_elem = op.result_matrix()(0).data_size() / n_points;
740 
741  auto &shape_matrix = op.result_matrix();
742  ArenaVec<double> ref_vec(n_points * n_dofs, op.result_matrix()(0).arena());
743  ArenaVec<double> elem_vec(n_elem, op.result_matrix()(0).arena());
744  for (uint i=0; i<n_elem; ++i) {
745  elem_vec(i) = 1.0;
746  }
747  ArenaOVec<double> ref_ovec(ref_vec, ref_vec.data_size());
748  ArenaOVec<double> elem_ovec(elem_vec, elem_vec.data_size());
749  for (uint i_dof=0; i_dof<n_dofs; ++i_dof)
750  for (uint i_p=0; i_p<n_points; ++i_p)
751  ref_vec(i_dof * n_points + i_p) = shape_values[i_dof][i_p];
752  ArenaOVec<double> shape_ovec = elem_ovec * ref_ovec;
753  shape_matrix(0) = shape_ovec.get_vec();
754  }
755  template<unsigned int dim>
756  static inline void ptop_scalar_shape_grads(std::vector<ElOp<3>> &operations,
757  std::vector< std::vector<arma::mat> > ref_shape_grads, uint scalar_shape_grads_op_idx) {
758  auto &op = operations[scalar_shape_grads_op_idx];
759  auto &inv_jac_vec = operations[ op.input_ops()[0] ].result_matrix();
760  uint n_points = ref_shape_grads.size();
761  uint n_dofs = ref_shape_grads[0].size();
762 
763  Eigen::Matrix<ArenaVec<double>, dim, 1> ref_grads_vec;
764  Eigen::Matrix<ArenaOVec<double>, dim, 1> ref_grads_ovec;
765  for (uint i=0; i<ref_grads_vec.rows(); ++i) {
766  ref_grads_vec(i) = ArenaVec<double>(n_points * n_dofs, op.result_matrix()(0).arena());
767  for (uint i_dof=0; i_dof<n_dofs; ++i_dof)
768  for (uint i_p=0; i_p<n_points; ++i_p)
769  ref_grads_vec(i)(i_dof * n_points + i_p) = ref_shape_grads[i_dof][i_p](i);
770  ref_grads_ovec(i) = ArenaOVec(ref_grads_vec(i), ref_grads_vec(i).data_size());
771  }
772 
773  Eigen::Matrix<ArenaOVec<double>, dim, 3> inv_jac_ovec;
774  for (uint i=0; i<dim*3; ++i) {
775  inv_jac_ovec(i) = ArenaOVec(inv_jac_vec(i), inv_jac_vec(i).data_size());
776  }
777 
778  auto &result_vec = op.result_matrix();
779  Eigen::Matrix<ArenaOVec<double>, 3, 1> result_ovec = inv_jac_ovec.transpose() * ref_grads_ovec;
780  for (uint i=0; i<3; ++i) {
781  result_vec(i) = result_ovec(i).get_vec();
782  }
783  }
784 };
785 
786 
787 
788 /// Defines reinit operations on side points.
789 struct side_reinit {
790  // element operations
791  template<unsigned int dim>
792  static inline void elop_el_jac(std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
793  // result matrix(spacedim, dim), input matrix(spacedim, dim+1)
794  auto &op = operations[FeSide::SideOps::opElJac];
795  auto &jac_value = op.result_matrix();
796  const auto &coords_value = operations[ op.input_ops()[0] ].result_matrix();
797  for (unsigned int i=0; i<3; i++)
798  for (unsigned int j=0; j<dim; j++)
799  jac_value(i,j) = coords_value(i,j+1) - coords_value(i,0);
800  }
801  template<unsigned int dim>
802  static inline void elop_el_inv_jac(std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
803  auto &op = operations[FeSide::SideOps::opElInvJac];
804  auto &inv_jac_value = op.result_matrix();
805  const auto &jac_value = operations[ op.input_ops()[0] ].result_matrix();
806  inv_jac_value = eigen_arena_tools::inverse<3, dim>(jac_value);
807  }
808  template<unsigned int dim>
809  static inline void elop_sd_jac(std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
810  // result matrix(spacedim, dim), input matrix(spacedim, dim+1)
811  auto &op = operations[FeSide::SideOps::opSideJac];
812  auto &jac_value = op.result_matrix();
813  const auto &coords_value = operations[ op.input_ops()[0] ].result_matrix();
814  for (unsigned int i=0; i<3; i++)
815  for (unsigned int j=0; j<dim-1; j++)
816  jac_value(i,j) = coords_value(i,j+1) - coords_value(i,0);
817  }
818  template<unsigned int dim>
819  static inline void elop_sd_jac_det(std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
820  // result double, input matrix(spacedim, dim)
821  auto &op = operations[FeSide::SideOps::opSideJacDet];
822  auto &jac_det_value = op.result_matrix();
823  const auto &jac_value = operations[ op.input_ops()[0] ].result_matrix();
824  jac_det_value(0,0) = eigen_arena_tools::determinant<3, dim-1>(jac_value).abs();
825  }
826 
827  // expansion operations
828 // static inline void expd_el_coords(std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
829 // auto &op = operations[FeSide::SideOps::opExpansionElCoords];
830 // common_reinit::expand_data(op, op_results, el_table);
831 // }
832 // static inline void expd_el_jac(std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
833 // auto &op = operations[FeSide::SideOps::opExpansionElJac];
834 // common_reinit::expand_data(op, op_results, el_table);
835 // }
836 // static inline void expd_el_inv_jac(std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
837 // auto &op = operations[FeSide::SideOps::opExpansionElInvJac];
838 // common_reinit::expand_data(op, op_results, el_table);
839 // }
840 //
841 // static inline void expd_sd_coords(std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
842 // auto &op = operations[FeSide::SideOps::opExpansionSideCoords];
843 // common_reinit::expand_data(op, op_results, el_table);
844 // }
845 // template<unsigned int dim>
846 // static inline void expd_sd_jac(std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
847 // auto &op = operations[FeSide::SideOps::opExpansionSideJac];
848 // common_reinit::expand_data(op, op_results, el_table);
849 // }
850 // static inline void expd_sd_jac_det(std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
851 // auto &op = operations[FeSide::SideOps::opExpansionSideJacDet];
852 // common_reinit::expand_data(op, op_results, el_table);
853 // }
854 
855  // Point operations
856  static inline void ptop_coords(FMT_UNUSED std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
857  // Implement
858  }
859  static inline void ptop_weights(std::vector<ElOp<3>> &operations, PatchArena *arena, const std::vector<double> &point_weights) {
860  auto &op = operations[FeSide::SideOps::opWeights];
861  op.allocate_result(point_weights.size(), *arena);
862  auto &weights_value = op.result_matrix();
863  for (uint i=0; i<point_weights.size(); ++i)
864  weights_value(0,0)(i) = point_weights[i];
865  }
866  static inline void ptop_JxW(std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
867  auto &op = operations[FeSide::SideOps::opJxW];
868  auto &weights_value = operations[ op.input_ops()[0] ].result_matrix();
869  auto &jac_det_value = operations[ op.input_ops()[1] ].result_matrix();
870  ArenaOVec<double> weights_ovec( weights_value(0,0), weights_value(0,0).data_size() );
871  ArenaOVec<double> jac_det_ovec( jac_det_value(0,0), jac_det_value(0,0).data_size() );
872  ArenaOVec<double> jxw_ovec = jac_det_ovec * weights_ovec;
873  auto &jxw_value = op.result_matrix();
874  jxw_value(0,0) = jxw_ovec.get_vec();
875  }
876  template<unsigned int dim>
877  static inline void ptop_normal_vec(std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
878  auto &op = operations[FeSide::SideOps::opNormalVec];
879  auto normal_value = op.value<3, 1>(op_results);
880  auto inv_jac_mat_value = operations[ op.input_ops()[0] ].value<dim, 3>(op_results);
881  normal_value = inv_jac_mat_value.transpose() * RefElement<dim>::normal_vector_array( el_table(3) );
882 
883  ArrayDbl norm_vec;
884  norm_vec.resize(normal_value(0).rows());
885  Eigen::VectorXd A(3);
886 
887  for (uint i=0; i<normal_value(0).rows(); ++i) {
888  A(0) = normal_value(0)(i);
889  A(1) = normal_value(1)(i);
890  A(2) = normal_value(2)(i);
891  norm_vec(i) = A.norm();
892  }
893  for (uint i=0; i<3; ++i) {
894  normal_value(i) /= norm_vec;
895  }
896  }
897  static inline void ptop_scalar_shape(std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table,
898  std::vector< std::vector< std::vector<double> > > shape_values, uint scalar_shape_op_idx) {
899  auto &op = operations[scalar_shape_op_idx];
900  uint n_points = shape_values[0].size();
901 
902  for (uint i_row=0; i_row<shape_values[0][0].size(); ++i_row) {
903  ArrayDbl &result_row = op_results( op.result_row()+i_row );
904  for (uint i_pt=0; i_pt<result_row.rows(); ++i_pt)
905  result_row(i_pt) = shape_values[el_table(3)(i_pt)][i_pt % n_points][i_row];
906  }
907  }
908  template<unsigned int dim>
909  static inline void ptop_scalar_shape_grads(std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table,
910  std::vector< std::vector< std::vector<arma::mat> > > ref_shape_grads, uint scalar_shape_grads_op_idx) {
911  auto &op = operations[scalar_shape_grads_op_idx];
912  uint n_points = ref_shape_grads[0].size();
913  uint n_dofs = ref_shape_grads[0][0].size();
914 
915  Eigen::Vector<ArrayDbl, Eigen::Dynamic> ref_shape_grads_expd;
916  ref_shape_grads_expd.resize(dim*n_dofs);
917  for (uint i=0; i<ref_shape_grads_expd.rows(); ++i)
918  ref_shape_grads_expd(i).resize(op_results(0).rows());
919 
920  for (uint i_dof=0; i_dof<n_dofs; ++i_dof)
921  for (uint i_c=0; i_c<dim; ++i_c) {
922  ArrayDbl &shape_grad_row = ref_shape_grads_expd(i_dof*dim+i_c);
923  for (uint i_pt=0; i_pt<shape_grad_row.rows(); ++i_pt)
924  shape_grad_row(i_pt) = ref_shape_grads[el_table(3)(i_pt)][i_pt % n_points][i_dof][i_c];
925  }
926 
927  auto inv_jac_value = operations[ op.input_ops()[0] ].value<dim, 3>(op_results);
928  for (uint i_dof=0; i_dof<n_dofs; ++i_dof) {
929  auto shape_grad_value = op.value<3, 1>(op_results, i_dof);
930  Eigen::Map<Eigen::Matrix<ArrayDbl, dim, 1>> ref_shape_grads_dof_value(ref_shape_grads_expd.data() + dim*i_dof, dim, 1);
931  shape_grad_value = inv_jac_value.transpose() * ref_shape_grads_dof_value;
932  }
933  }
934 };
935 
936 
937 // template specialization
938 template<>
939 inline void side_reinit::elop_sd_jac<1>(FMT_UNUSED std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
940 }
941 
942 template<>
943 inline void side_reinit::elop_sd_jac_det<1>(std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
944  auto &op = operations[FeSide::SideOps::opSideJacDet];
945  auto &result_vec = op.result_matrix();
946  for (uint i=0;i<result_vec(0,0).data_size(); ++i) {
947  result_vec(0,0)(i) = 1.0;
948  }
949 }
950 
951 //template<>
952 //inline void side_reinit::expd_sd_jac<1>(FMT_UNUSED std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
953 //}
954 
955 
956 
957 namespace FeBulk {
958 
959  /// Bulk data specialization, order of item in operations_ vector corresponds to the BulkOps enum
960  template<unsigned int spacedim = 3>
961  class PatchPointValues : public ::PatchPointValues<spacedim> {
962  public:
963  /// Constructor
964  PatchPointValues(uint dim, uint quad_order, AssemblyArena &asm_arena)
965  : ::PatchPointValues<spacedim>(dim, asm_arena) {
966  this->quad_ = new QGauss(dim, 2*quad_order);
967  switch (dim) {
968  case 1:
969  init<1>();
970  break;
971  case 2:
972  init<2>();
973  break;
974  case 3:
975  init<3>();
976  break;
977  }
978  }
979 
980  private:
981  /// Initialize operations vector
982  template<unsigned int dim>
983  void init() {
984  // First step: adds element values operations
985  /*auto &el_coords =*/ this->make_new_op( {spacedim, this->dim_+1}, &common_reinit::op_base, {}, OpSizeType::elemOp );
986 
987  /*auto &el_jac =*/ this->make_new_op( {spacedim, this->dim_}, &bulk_reinit::elop_jac<dim>, {BulkOps::opElCoords}, OpSizeType::elemOp );
988 
989  /*auto &el_inv_jac =*/ this->make_new_op( {this->dim_, spacedim}, &bulk_reinit::elop_inv_jac<dim>, {BulkOps::opJac}, OpSizeType::elemOp );
990 
991  /*auto &el_jac_det =*/ this->make_new_op( {1}, &bulk_reinit::elop_jac_det<dim>, {BulkOps::opJac}, OpSizeType::elemOp );
992 
993  // Second step: adds point values operations
994  /*auto &pt_coords =*/ this->make_new_op( {spacedim}, &bulk_reinit::ptop_coords, {} );
995 
996  // use lambda reinit function
997  const std::vector<double> &point_weights_vec = this->quad_->get_weights();
998  auto lambda_weights = [this, point_weights_vec](std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
999  bulk_reinit::ptop_weights(operations, this->patch_arena_, point_weights_vec);
1000  };
1001  /*auto &weights =*/ this->make_fixed_op( {1}, lambda_weights );
1002 
1003  /*auto &JxW =*/ this->make_new_op( {1}, &bulk_reinit::ptop_JxW, {BulkOps::opWeights, BulkOps::opJacDet} );
1004  }
1005  };
1006 
1007 } // closing namespace FeBulk
1008 
1009 
1010 
1011 namespace FeSide {
1012 
1013 /// Bulk Side specialization, order of item in operations_ vector corresponds to the SideOps enum
1014  template<unsigned int spacedim = 3>
1015  class PatchPointValues : public ::PatchPointValues<spacedim> {
1016  public:
1017  /// Constructor
1018  PatchPointValues(uint dim, uint quad_order, AssemblyArena &asm_arena)
1019  : ::PatchPointValues<spacedim>(dim, asm_arena) {
1020  this->quad_ = new QGauss(dim-1, 2*quad_order);
1021  switch (dim) {
1022  case 1:
1023  init<1>();
1024  break;
1025  case 2:
1026  init<2>();
1027  break;
1028  case 3:
1029  init<3>();
1030  break;
1031  }
1032  }
1033 
1034  private:
1035  /// Initialize operations vector
1036  template<unsigned int dim>
1037  void init() {
1038  // First step: adds element values operations
1039  /*auto &el_coords =*/ this->make_new_op( {spacedim, this->dim_+1}, &common_reinit::op_base, {}, OpSizeType::elemOp );
1040 
1041  /*auto &el_jac =*/ this->make_new_op( {spacedim, this->dim_}, &side_reinit::elop_el_jac<dim>, {SideOps::opElCoords}, OpSizeType::elemOp );
1042 
1043  /*auto &el_inv_jac =*/ this->make_new_op( {this->dim_, spacedim}, &side_reinit::elop_el_inv_jac<dim>, {SideOps::opElJac}, OpSizeType::elemOp );
1044 
1045  /*auto &sd_coords =*/ this->make_new_op( {spacedim, this->dim_}, &common_reinit::op_base, {}, OpSizeType::elemOp );
1046 
1047  /*auto &sd_jac =*/ this->make_new_op( {spacedim, this->dim_-1}, &side_reinit::elop_sd_jac<dim>, {SideOps::opSideCoords}, OpSizeType::elemOp );
1048 
1049  /*auto &sd_jac_det =*/ this->make_new_op( {1}, &side_reinit::elop_sd_jac_det<dim>, {SideOps::opSideJac}, OpSizeType::elemOp );
1050 
1051  // Second step: adds expand operations (element values to point values)
1052 // this->make_expansion( el_coords, {spacedim, this->dim_+1}, &side_reinit::expd_el_coords );
1053 //
1054 // this->make_expansion( el_jac, {spacedim, this->dim_}, &side_reinit::expd_el_jac );
1055 //
1056 // this->make_expansion( el_inv_jac, {this->dim_, spacedim}, &side_reinit::expd_el_inv_jac );
1057 //
1058 // this->make_expansion( sd_coords, {spacedim, this->dim_}, &side_reinit::expd_sd_coords );
1059 //
1060 // this->make_expansion( sd_jac, {spacedim, this->dim_-1}, &side_reinit::expd_sd_jac<dim> );
1061 //
1062 // this->make_expansion( sd_jac_det, {1}, &side_reinit::expd_sd_jac_det );
1063 
1064  // Third step: adds point values operations
1065  /*auto &coords =*/ this->make_new_op( {spacedim}, &side_reinit::ptop_coords, {} );
1066 
1067  // use lambda reinit function
1068  const std::vector<double> &point_weights_vec = this->quad_->get_weights();
1069  auto lambda_weights = [this, point_weights_vec](std::vector<ElOp<3>> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
1070  side_reinit::ptop_weights(operations, this->patch_arena_, point_weights_vec);
1071  };
1072  /*auto &weights =*/ this->make_fixed_op( {1}, lambda_weights );
1073 
1075 
1076  /*auto &normal_vec =*/ this->make_new_op( {spacedim}, &side_reinit::ptop_normal_vec<dim>, {SideOps::opElInvJac} );
1077  }
1078  };
1079 
1080 } // closing namespace FeSide
1081 
1082 
1083 
1084 #endif /* PATCH_POINT_VALUES_HH_ */
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
Outer product - only proposal of multi operator.
Definition: arena_vec.hh:237
ArenaVec< T > get_vec() const
Convert ArenaOVec to ArenaVec and its.
Definition: arena_vec.hh:274
size_t data_size() const
Getter for data_size_.
Definition: arena_vec.hh:104
Class represents element or FE operations.
Eigen::Matrix< ArenaVec< double >, Eigen::Dynamic, Eigen::Dynamic > & result_matrix()
Return map referenced result as Eigen::Matrix.
void allocate_result(size_t data_size, PatchArena &arena)
std::string format_shape() const
std::vector< uint > input_ops_
Indices of operations in PatchPointValues::operations_ vector on which ElOp is depended.
std::vector< uint > set_shape_vec(std::initializer_list< uint > shape) const
Aligns shape_vec to 2 items (equal to matrix number of dimensions)
uint dim_
Dimension.
Eigen::Matrix< ArenaVec< double >, Eigen::Dynamic, Eigen::Dynamic > result_
Result matrix of operation.
std::vector< uint > shape_
Shape of stored data (size of vector or number of rows and cols of matrix)
uint result_row() const
Getter for result_row_.
const std::vector< uint > & shape() const
Getter for shape_.
const std::vector< uint > & input_ops() const
Getter for input_ops_.
OpSizeType size_type_
Type of operation by size of vector (element, point or fixed size)
uint result_row_
First row to scalar, vector or matrix result TODO replace.
const Eigen::Matrix< ArenaVec< double >, Eigen::Dynamic, Eigen::Dynamic > & result_matrix() const
Same as previous but return const reference.
uint n_dofs() const
Getter for n_dofs_.
ReinitFunction reinit_func
Pointer to patch reinit function of element data table specialized by operation.
uint n_dofs_
Number of DOFs of FE operations (or 1 in case of element operations)
uint dim() const
Getter for dimension.
Eigen::Map< Eigen::Matrix< ArrayDbl, dim1, dim2 > > value(TableDbl &op_results, uint i_dof=0) const
Return map referenced Eigen::Matrix of given dimension.
Eigen::Map< Eigen::Vector< double, Eigen::Dynamic > > vector_value(TableDbl &op_results) const
Return map referenced Eigen::Matrix of given dimensions.
uint n_comp() const
ElOp(uint dim, std::initializer_list< uint > shape, uint result_row, ReinitFunction reinit_f, OpSizeType size_type, std::vector< uint > input_ops={}, uint n_dofs=1)
Eigen::Map< Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > > matrix_value(TableDbl &op_results, uint dim1, uint dim2) const
Return map referenced Eigen::Matrix of given dimensions.
OpSizeType size_type() const
Getter for size_type_.
void reinit_function(std::vector< ElOp< spacedim >> &operations, TableDbl &data_table, TableInt &int_table)
Call reinit function on element table if function is defined.
Bulk data specialization, order of item in operations_ vector corresponds to the BulkOps enum.
PatchPointValues(uint dim, uint quad_order, AssemblyArena &asm_arena)
Constructor.
void init()
Initialize operations vector.
Bulk Side specialization, order of item in operations_ vector corresponds to the SideOps enum.
PatchPointValues(uint dim, uint quad_order, AssemblyArena &asm_arena)
Constructor.
void init()
Initialize operations vector.
void print_data_tables(ostream &stream, bool points, bool ints) const
PatchPointValues(uint dim, AssemblyArena &asm_arena)
ElOp< spacedim > & make_expansion(ElOp< spacedim > &el_op, std::initializer_list< uint > shape, ReinitFunction reinit_f)
AssemblyArena & asm_arena_
Reference to global assembly arena of PatchFeValues.
std::vector< OpSizeType > row_sizes_
hold sizes of rows by type of operation
void initialize(uint int_cols)
uint n_rows() const
Getter for n_rows_.
std::vector< uint > points_map_
Map of point patch indices to point_vals_ and int_vals_ tables.
uint register_side_point(uint elem_table_row, uint value_patch_idx, uint elem_idx, uint side_idx)
uint n_rows_
Number of columns of point_vals table.
uint register_bulk_point(uint elem_table_row, uint value_patch_idx, uint elem_idx, uint i_point_on_elem)
uint register_element(arma::mat coords, uint element_patch_idx)
void resize_tables(uint n_elems, uint n_points)
Resize data tables. Method is called before reinit of patch.
Tensor tensor_val(uint result_row, uint point_idx) const
uint register_side(arma::mat elm_coords, arma::mat side_coords)
ElOp< spacedim > & make_fe_op(std::initializer_list< uint > shape, ReinitFunction reinit_f, std::vector< uint > input_ops_vec, uint n_dofs, OpSizeType size_type=pointOp)
std::vector< ElOp< spacedim > > operations_
Vector of all defined operations.
Quadrature * get_quadrature() const
Getter for quadrature.
uint n_elems() const
Getter for n_elems_.
uint i_elem_
Index of registered element in table, helper value used during patch creating.
void reset()
Reset number of columns (points and elements)
uint n_points() const
Getter for n_points_.
ElOp< spacedim > & make_fixed_op(std::initializer_list< uint > shape, ReinitFunction reinit_f)
std::vector< uint > elements_map_
Map of element patch indices to el_vals_ table.
uint dim() const
Getter for dim_.
Tensor tensor_value(uint op_idx, uint point_idx, uint i_dof=0) const
Scalar scalar_value(uint op_idx, uint point_idx, uint i_dof=0) const
void init_finalize(PatchArena *patch_arena)
Vector vector_val(uint result_row, uint point_idx) const
uint n_points_
Number of points in patch.
void print_operations(ostream &stream, uint bulk_side) const
Quadrature * quad_
Quadrature of given dimension and order passed in constructor.
Scalar scalar_val(uint result_row, uint point_idx) const
ElOp< spacedim > & make_new_op(std::initializer_list< uint > shape, ReinitFunction reinit_f, std::vector< uint > input_ops_vec, OpSizeType size_type=pointOp)
Vector vector_value(uint op_idx, uint point_idx, uint i_dof=0) const
PatchArena * patch_arena_
Pointer to global patch arena of PatchFeValues.
uint n_elems_
Number of elements in patch.
Symmetric Gauss-Legendre quadrature formulae on simplices.
Base class for quadrature rules on simplices in arbitrary dimensions.
Definition: quadrature.hh:48
const std::vector< double > & get_weights() const
Return a reference to the whole array of weights.
Definition: quadrature.hh:111
static Eigen::Vector< ArrayDbl, dim > normal_vector_array(Eigen::Array< uint, Eigen::Dynamic, 1 > loc_side_idx_array)
Definition: ref_element.cc:279
Store finite element data on the actual patch such as shape function values, gradients,...
Eigen::Vector< ArrayDbl, Eigen::Dynamic > TableDbl
Definition: eigen_tools.hh:36
Eigen::Array< double, Eigen::Dynamic, 1 > ArrayDbl
Definitions of Eigen structures.
Definition: eigen_tools.hh:34
Eigen::Vector< ArrayInt, Eigen::Dynamic > TableInt
Definition: eigen_tools.hh:37
Class ElementValues calculates data related to transformation of reference cell to actual cell (Jacob...
unsigned int uint
ArmaMat< double, N, M > mat
Definition: armor.hh:936
@ opCoords
operations evaluated on quadrature points
@ opInvJac
inverse Jacobian
@ opWeights
weight of quadrature point
@ opJac
Jacobian of element.
@ opElCoords
operations evaluated on elements
@ opJxW
JxW value of quadrature point.
@ opNormalVec
normal vector of quadrature point
@ opSideCoords
operations evaluated on sides
@ opJxW
JxW value of quadrature point.
@ opCoords
operation executed expansion to quadrature point (value of element / side > values on quadrature poin...
@ opElCoords
operations evaluated on elements
@ opElJac
Jacobian of element.
@ opWeights
weight of quadrature point
@ opSideJac
Jacobian of element.
ArenaVec< double > determinant(const Eigen::Matrix< ArenaVec< double >, m, n > &A)
Calculates determinant of a rectangular matrix.
void resize_table(typename Eigen::Vector< ET, Eigen::Dynamic > &table, uint size)
Resize vector of Eigen::Array to given size.
Definition: eigen_tools.hh:226
std::function< void(std::vector< ElOp< 3 > > &, TableDbl &, TableInt &)> ReinitFunction
Type for conciseness.
double Scalar
OpSizeType
Distinguishes operations by type and size of output rows.
@ pointOp
operation is evaluated on quadrature points
@ elemOp
operation is evaluated on elements or sides
@ fixedSizeOp
operation has fixed size and it is filled during initialization
#define FMT_UNUSED
Definition: posix.h:75
Definitions of particular quadrature rules on simplices.
Defines reinit operations on bulk points.
static void elop_jac(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_scalar_shape_grads(std::vector< ElOp< 3 >> &operations, std::vector< std::vector< arma::mat > > ref_shape_grads, uint scalar_shape_grads_op_idx)
static void ptop_JxW(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_coords(FMT_UNUSED std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void elop_jac_det(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_weights(std::vector< ElOp< 3 >> &operations, PatchArena *arena, const std::vector< double > &point_weights)
static void elop_inv_jac(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_scalar_shape(std::vector< ElOp< 3 >> &operations, std::vector< std::vector< double > > shape_values, uint scalar_shape_op_idx)
Defines common functionality of reinit operations.
static void expand_data(ElOp< 3 > &op, TableDbl &op_results, TableInt &el_table)
static void op_base(FMT_UNUSED std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
Defines reinit operations on side points.
static void elop_sd_jac_det(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_normal_vec(std::vector< ElOp< 3 >> &operations, TableDbl &op_results, TableInt &el_table)
static void ptop_weights(std::vector< ElOp< 3 >> &operations, PatchArena *arena, const std::vector< double > &point_weights)
static void ptop_scalar_shape(std::vector< ElOp< 3 >> &operations, TableDbl &op_results, TableInt &el_table, std::vector< std::vector< std::vector< double > > > shape_values, uint scalar_shape_op_idx)
static void ptop_JxW(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_coords(FMT_UNUSED std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void elop_sd_jac(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void ptop_scalar_shape_grads(std::vector< ElOp< 3 >> &operations, TableDbl &op_results, TableInt &el_table, std::vector< std::vector< std::vector< arma::mat > > > ref_shape_grads, uint scalar_shape_grads_op_idx)
static void elop_el_jac(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)
static void elop_el_inv_jac(std::vector< ElOp< 3 >> &operations, FMT_UNUSED TableDbl &op_results, FMT_UNUSED TableInt &el_table)