Flow123d  DF_patch_fe_data_tables-47aef41
patch_fe_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_fe_values.hh
15  * @brief Class FEValues calculates finite element data on the actual
16  * cells such as shape function values, gradients, Jacobian of
17  * the mapping from the reference cell etc.
18  * @author Jan Stebel, David Flanderka
19  */
20 
21 #ifndef PATCH_FE_VALUES_HH_
22 #define PATCH_FE_VALUES_HH_
23 
24 
25 #include <string.h> // for memcpy
26 #include <algorithm> // for swap
27 #include <new> // for operator new[]
28 #include <string> // for operator<<
29 #include <vector> // for vector
30 #include "fem/element_values.hh" // for ElementValues
31 #include "fem/fe_values.hh" // for FEValuesBase
32 #include "fem/fe_values_views.hh" // for FEValuesViews
33 #include "fem/fe_system.hh" // for FESystem
34 #include "fem/eigen_tools.hh"
36 #include "fem/mapping_p1.hh"
37 #include "mesh/ref_element.hh" // for RefElement
38 #include "mesh/accessors.hh"
39 #include "fem/update_flags.hh" // for UpdateFlags
41 #include "fields/eval_subset.hh"
42 #include "system/arena_resource.hh"
43 
44 template<unsigned int spacedim> class PatchFEValues;
45 
46 
47 
48 typedef typename std::vector< std::array<uint, 3> > DimPointTable; ///< Holds triplet (dim; bulk/side; idx of point in subtable)
49 
50 
51 template <class ValueType>
52 class ElQ {
53 public:
54  /// Forbidden default constructor
55  ElQ() = delete;
56 
57  /// Constructor
58  ElQ(PatchPointValues<3> &patch_point_vals, unsigned int begin)
59  : patch_point_vals_(patch_point_vals), begin_(begin) {}
60 
61  ValueType operator()(FMT_UNUSED const BulkPoint &point);
62 
63  ValueType operator()(FMT_UNUSED const SidePoint &point);
64 
65 private:
66  PatchPointValues<3> &patch_point_vals_; ///< Reference to PatchPointValues
67  unsigned int begin_; /// Index of the first component of the bulk Quantity. Size is given by ValueType
68 };
69 
70 
71 template <class ValueType>
72 class FeQ {
73 public:
74  /// Forbidden default constructor
75  FeQ() = delete;
76 
77  // Class similar to current FeView
78  FeQ(PatchPointValues<3> &patch_point_vals, unsigned int begin, unsigned int n_dofs)
79  : patch_point_vals_(patch_point_vals), begin_(begin), n_dofs_(n_dofs) {}
80 
81 
82  ValueType operator()(FMT_UNUSED unsigned int shape_idx, FMT_UNUSED const BulkPoint &point);
83 
84  ValueType operator()(FMT_UNUSED unsigned int shape_idx, FMT_UNUSED const SidePoint &point);
85 
86  // Implementation for EdgePoint, SidePoint, and JoinPoint shoud have a common implementation
87  // resolving to side values
88 
89 private:
90  PatchPointValues<3> &patch_point_vals_; ///< Reference to PatchPointValues
91  unsigned int begin_; ///< Index of the first component of the Quantity. Size is given by ValueType
92  unsigned int n_dofs_; ///< Number of DOFs
93 };
94 
95 
96 template <class ValueType>
98 public:
99  /// Default constructor
101  : patch_point_vals_bulk_(nullptr), patch_point_vals_side_(nullptr), join_idx_(-1) {}
102 
103  /**
104  * Constructor
105  *
106  * @param patch_point_vals_bulk Pointer to PatchPointValues bulk object.
107  * @param patch_point_vals_side Pointer to PatchPointValues side object.
108  * @param begin Index of the first component of the bulk Quantity.
109  * @param begin_side Index of the first component of the side Quantity.
110  * @param n_dofs_bulk Number of DOFs of bulk (lower-dim) element.
111  * @param n_dofs_side Number of DOFs of side (higher-dim) element.
112  * @param join_idx Index function.
113  */
114  JoinShapeAccessor(PatchPointValues<3> *patch_point_vals_bulk, PatchPointValues<3> *patch_point_vals_side,
115  unsigned int begin, unsigned int begin_side, unsigned int n_dofs_bulk, unsigned int n_dofs_side, unsigned int join_idx)
116  : patch_point_vals_bulk_(patch_point_vals_bulk), patch_point_vals_side_(patch_point_vals_side), begin_(begin),
117  begin_side_(begin_side), n_dofs_high_(n_dofs_side), n_dofs_low_(n_dofs_bulk), join_idx_(join_idx) {
118  //ASSERT( (patch_point_vals_bulk->dim()==2) || (patch_point_vals_bulk->dim()==3) )(patch_point_vals_bulk->dim() ).error("Invalid dimension, must be 2 or 3!");
119  }
120 
121  /// Return global index of DOF
122  inline unsigned int join_idx() const {
123  return join_idx_;
124  }
125 
126  /// Return local index of DOF (on low / high-dim) - should be private method
127  inline unsigned int local_idx() const {
128  if (this->is_high_dim()) return (join_idx_ - n_dofs_low_);
129  else return join_idx_;
130  }
131 
132  inline unsigned int n_dofs_low() const {
133  return n_dofs_low_;
134  }
135 
136  inline unsigned int n_dofs_high() const {
137  return n_dofs_high_;
138  }
139 
140  inline unsigned int n_dofs_both() const {
141  return n_dofs_high_ + n_dofs_low_;
142  }
143 
144  inline bool is_high_dim() const {
145  return (join_idx_ >= n_dofs_low_);
146  }
147 
148  /// Iterates to next item.
149  inline void inc() {
150  join_idx_++;
151  }
152 
153  /// Comparison of accessors.
154  bool operator==(const JoinShapeAccessor<ValueType>& other) const {
155  return (join_idx_ == other.join_idx_);
156  }
157 
158 
159  ValueType operator()(const BulkPoint &point);
160 
161  ValueType operator()(const SidePoint &point);
162 
163 private:
164  // attributes:
165  PatchPointValues<3> *patch_point_vals_bulk_; ///< Pointer to bulk PatchPointValues
166  PatchPointValues<3> *patch_point_vals_side_; ///< Pointer to side PatchPointValues
167  unsigned int begin_; ///< Index of the first component of the bulk Quantity. Size is given by ValueType
168  unsigned int begin_side_; ///< Index of the first component of the side Quantity. Size is given by ValueType
169  unsigned int n_dofs_high_; ///< Number of DOFs on high-dim element
170  unsigned int n_dofs_low_; ///< Number of DOFs on low-dim element
171  unsigned int join_idx_; ///< Index of processed DOF
172 };
173 
174 
175 template<unsigned int dim>
177 {
178 protected:
179  // Default constructor
181  {}
182 
183  /// Return FiniteElement of \p component_idx for FESystem or \p fe for other types
184  template<unsigned int FE_dim>
185  std::shared_ptr<FiniteElement<FE_dim>> fe_comp(std::shared_ptr< FiniteElement<FE_dim> > fe, uint component_idx) {
186  FESystem<FE_dim> *fe_sys = dynamic_cast<FESystem<FE_dim>*>( fe.get() );
187  if (fe_sys != nullptr) {
188  return fe_sys->fe()[component_idx];
189  } else {
190  ASSERT_EQ(component_idx, 0).warning("Non-zero component_idx can only be used for FESystem.");
191  return fe;
192  }
193  }
194 
195  /**
196  * @brief Precomputed values of basis functions at the bulk quadrature points.
197  *
198  * Dimensions: (no. of quadrature points)
199  * x (no. of dofs)
200  * x (no. of components in ref. cell)
201  */
202  template<unsigned int FE_dim>
204  std::vector<std::vector<arma::vec> > ref_shape_vals( q->size(), vector<arma::vec>(fe->n_dofs()) );
205 
206  arma::mat shape_values(fe->n_dofs(), fe->n_components());
207  for (unsigned int i=0; i<q->size(); i++)
208  {
209  for (unsigned int j=0; j<fe->n_dofs(); j++)
210  {
211  for (unsigned int c=0; c<fe->n_components(); c++)
212  shape_values(j,c) = fe->shape_value(j, q->point<FE_dim>(i), c);
213 
214  ref_shape_vals[i][j] = trans(shape_values.row(j));
215  }
216  }
217 
218  return ref_shape_vals;
219  }
220 
221  /**
222  * @brief Precomputed values of basis functions at the side quadrature points.
223  *
224  * Dimensions: (sides)
225  * x (no. of quadrature points)
226  * x (no. of dofs)
227  * x (no. of components in ref. cell)
228  */
229  template<unsigned int FE_dim>
232 
233  arma::mat shape_values(fe->n_dofs(), fe->n_components());
234 
235  for (unsigned int sid=0; sid<FE_dim+1; sid++) {
236  auto quad = q->make_from_side<dim>(sid);
237  for (unsigned int i=0; i<quad.size(); i++)
238  {
239  for (unsigned int j=0; j<fe->n_dofs(); j++)
240  {
241  for (unsigned int c=0; c<fe->n_components(); c++) {
242  shape_values(j,c) = fe->shape_value(j, quad.template point<FE_dim>(i), c);
243  }
244 
245  ref_shape_vals[sid][i][j] = trans(shape_values.row(j));
246  }
247  }
248  }
249 
250  return ref_shape_vals;
251  }
252 
253 };
254 
255 template<unsigned int dim>
256 class BulkValues : public BaseValues<dim>
257 {
258 public:
259  /// Constructor
261  : BaseValues<dim>(), patch_point_vals_(patch_point_vals) {
262  ASSERT_EQ(patch_point_vals.dim(), dim);
263  fe_ = fe[Dim<dim>{}];
264  }
265 
266  /**
267  * @brief Register the product of Jacobian determinant and the quadrature
268  * weight at bulk quadrature points.
269  *
270  * @param quad Quadrature.
271  */
272  inline ElQ<Scalar> JxW()
273  {
275  return ElQ<Scalar>(patch_point_vals_, begin);
276  }
277 
278  /// Create bulk accessor of coords entity
280  {
282  return ElQ<Vector>(patch_point_vals_, begin);
283  }
284 
285 // inline ElQ<Tensor> jacobian(std::initializer_list<Quadrature *> quad_list)
286 // {}
287 
288  /// Create bulk accessor of jac determinant entity
290  {
292  return ElQ<Scalar>(patch_point_vals_, begin);
293  }
294 
295  /**
296  * @brief Return the value of the @p function_no-th shape function at
297  * the @p p bulk quadrature point.
298  *
299  * @param component_idx Number of the shape function.
300  */
301  inline FeQ<Scalar> scalar_shape(uint component_idx = 0)
302  {
303  auto fe_component = this->fe_comp(fe_, component_idx);
304  ASSERT_EQ(fe_component->fe_type(), FEType::FEScalar).error("Type of FiniteElement of scalar_shape accessor must be FEScalar!\n");
305 
306  // use lambda reinit function
307  std::vector< std::vector<double> > shape_values( patch_point_vals_.get_quadrature()->size(), vector<double>(fe_component->n_dofs()) );
308  auto ref_shape_vals = this->ref_shape_values_bulk(patch_point_vals_.get_quadrature(), fe_component);
309  for (unsigned int i = 0; i < patch_point_vals_.get_quadrature()->size(); i++)
310  for (unsigned int j = 0; j < fe_component->n_dofs(); j++) {
311  shape_values[i][j] = ref_shape_vals[i][j][0];
312  }
313  uint scalar_shape_op_idx = patch_point_vals_.operations_.size(); // index in operations_ vector
314  auto lambda_scalar_shape = [shape_values, scalar_shape_op_idx](std::vector<ElOp<3>> &operations, TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
315  bulk_reinit::ptop_scalar_shape(operations, op_results, shape_values, scalar_shape_op_idx);
316  };
317  auto &scalar_shape_bulk_op = patch_point_vals_.make_fe_op({1}, lambda_scalar_shape, {}, fe_component->n_dofs());
318  uint begin = scalar_shape_bulk_op.result_row();
319 
320  return FeQ<Scalar>(patch_point_vals_, begin, fe_component->n_dofs());
321  }
322 
323 // inline FeQ<Vector> vector_shape(uint component_idx = 0)
324 // {}
325 
326 // inline FeQ<Tensor> tensor_shape(uint component_idx = 0)
327 // {}
328 
329  /**
330  * @brief Return the value of the @p function_no-th gradient shape function at
331  * the @p p bulk quadrature point.
332  *
333  * @param component_idx Number of the shape function.
334  */
335  inline FeQ<Vector> grad_scalar_shape(uint component_idx=0)
336  {
337  auto fe_component = this->fe_comp(fe_, component_idx);
338  ASSERT_EQ(fe_component->fe_type(), FEType::FEScalar).error("Type of FiniteElement of grad_scalar_shape accessor must be FEScalar!\n");
339 
340  // use lambda reinit function
341  auto ref_shape_grads = this->ref_shape_gradients(fe_component);
342  uint scalar_shape_grads_op_idx = patch_point_vals_.operations_.size(); // index in operations_ vector
343  auto lambda_scalar_shape_grad = [ref_shape_grads, scalar_shape_grads_op_idx](std::vector<ElOp<3>> &operations, TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
344  bulk_reinit::ptop_scalar_shape_grads<dim>(operations, op_results, ref_shape_grads, scalar_shape_grads_op_idx);
345  };
346  auto &grad_scalar_shape_bulk_op = patch_point_vals_.make_fe_op({3}, lambda_scalar_shape_grad, {FeBulk::BulkOps::opInvJac}, fe_component->n_dofs());
347  uint begin = grad_scalar_shape_bulk_op.result_row();
348 
349  return FeQ<Vector>(patch_point_vals_, begin, fe_component->n_dofs());
350  }
351 
352 // inline FeQ<Tensor> grad_vector_shape(std::initializer_list<Quadrature *> quad_list, unsigned int i_comp=0)
353 // {}
354 
355 private:
356  /**
357  * @brief Precomputed gradients of basis functions at the quadrature points.
358  *
359  * Dimensions: (no. of quadrature points)
360  * x (no. of dofs)
361  * x ((dim_ of. ref. cell)x(no. of components in ref. cell))
362  */
365  std::vector<std::vector<arma::mat> > ref_shape_grads( q->size(), vector<arma::mat>(fe->n_dofs()) );
366 
367  arma::mat grad(dim, fe->n_components());
368  for (unsigned int i_pt=0; i_pt<q->size(); i_pt++)
369  {
370  for (unsigned int i_dof=0; i_dof<fe->n_dofs(); i_dof++)
371  {
372  grad.zeros();
373  for (unsigned int c=0; c<fe->n_components(); c++)
374  grad.col(c) += fe->shape_grad(i_dof, q->point<dim>(i_pt), c);
375 
376  ref_shape_grads[i_pt][i_dof] = grad;
377  }
378  }
379 
380  return ref_shape_grads;
381  }
382 
384  std::shared_ptr< FiniteElement<dim> > fe_;
385 };
386 
387 
388 template<unsigned int dim>
389 class SideValues : public BaseValues<dim>
390 {
391 public:
392  /// Constructor
394  : BaseValues<dim>(), patch_point_vals_(patch_point_vals) {
395  ASSERT_EQ(patch_point_vals.dim(), dim);
396  fe_ = fe[Dim<dim>{}];
397  }
398 
399  /// Same as BulkValues::JxW but register at side quadrature points.
400  inline ElQ<Scalar> JxW()
401  {
403  return ElQ<Scalar>(patch_point_vals_, begin);
404  }
405 
406  /**
407  * @brief Register the normal vector to a side at side quadrature points.
408  *
409  * @param quad Quadrature.
410  */
412  {
414  return ElQ<Vector>(patch_point_vals_, begin);
415  }
416 
417  /// Create side accessor of coords entity
419  {
421  return ElQ<Vector>(patch_point_vals_, begin);
422  }
423 
424  /// Create bulk accessor of jac determinant entity
426  {
428  return ElQ<Scalar>(patch_point_vals_, begin);
429  }
430 
431  /// Same as BulkValues::scalar_shape but register at side quadrature points.
432  inline FeQ<Scalar> scalar_shape(uint component_idx = 0)
433  {
434  auto fe_component = this->fe_comp(fe_, component_idx);
435  ASSERT_EQ(fe_component->fe_type(), FEType::FEScalar).error("Type of FiniteElement of scalar_shape accessor must be FEScalar!\n");
436 
437  // use lambda reinit function
439  dim+1,
441  );
442  auto ref_shape_vals = this->ref_shape_values_side(patch_point_vals_.get_quadrature(), fe_component);
443  for (unsigned int s=0; s<dim+1; ++s)
444  for (unsigned int i = 0; i < patch_point_vals_.get_quadrature()->size(); i++)
445  for (unsigned int j = 0; j < fe_component->n_dofs(); j++) {
446  shape_values[s][i][j] = ref_shape_vals[s][i][j][0];
447  }
448  uint scalar_shape_op_idx = patch_point_vals_.operations_.size(); // index in operations_ vector
449  auto lambda_scalar_shape = [shape_values, scalar_shape_op_idx](std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
450  side_reinit::ptop_scalar_shape(operations, op_results, el_table, shape_values, scalar_shape_op_idx);
451  };
452  auto &scalar_shape_bulk_op = patch_point_vals_.make_fe_op({1}, lambda_scalar_shape, {}, fe_component->n_dofs());
453  uint begin = scalar_shape_bulk_op.result_row();
454 
455  return FeQ<Scalar>(patch_point_vals_, begin, fe_component->n_dofs());
456  }
457 
458  /// Same as BulkValues::grad_scalar_shape but register at side quadrature points.
459  inline FeQ<Vector> grad_scalar_shape(uint component_idx=0)
460  {
461  auto fe_component = this->fe_comp(fe_, component_idx);
462  ASSERT_EQ(fe_component->fe_type(), FEType::FEScalar).error("Type of FiniteElement of grad_scalar_shape accessor must be FEScalar!\n");
463 
464  // use lambda reinit function
465  auto ref_shape_grads = this->ref_shape_gradients(fe_component);
466  uint scalar_shape_grads_op_idx = patch_point_vals_.operations_.size(); // index in operations_ vector
467  auto lambda_scalar_shape_grad = [ref_shape_grads, scalar_shape_grads_op_idx](std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
468  side_reinit::ptop_scalar_shape_grads<dim>(operations, op_results, el_table, ref_shape_grads, scalar_shape_grads_op_idx);
469  };
470  auto &grad_scalar_shape_side_op = patch_point_vals_.make_fe_op({3}, lambda_scalar_shape_grad, {FeSide::SideOps::opElInvJac}, fe_component->n_dofs());
471  uint begin = grad_scalar_shape_side_op.result_row();
472 
473  return FeQ<Vector>(patch_point_vals_, begin, fe_component->n_dofs());
474  }
475 
476 private:
477  /**
478  * @brief Precomputed gradients of basis functions at the quadrature points.
479  *
480  * Dimensions: (sides)
481  * x (no. of quadrature points)
482  * x (no. of dofs)
483  * x ((dim_ of. ref. cell)x(no. of components in ref. cell))
484  */
488 
489  arma::mat grad(dim, fe->n_components());
490  for (unsigned int sid=0; sid<dim+1; sid++) {
491  auto quad = q->make_from_side<dim>(sid);
492  for (unsigned int i_pt=0; i_pt<quad.size(); i_pt++)
493  {
494  for (unsigned int i_dof=0; i_dof<fe->n_dofs(); i_dof++)
495  {
496  grad.zeros();
497  for (unsigned int c=0; c<fe->n_components(); c++)
498  grad.col(c) += fe->shape_grad(i_dof, quad.template point<dim>(i_pt), c);
499 
500  ref_shape_grads[sid][i_pt][i_dof] = grad;
501  }
502  }
503  }
504 
505  return ref_shape_grads;
506  }
507 
509  std::shared_ptr< FiniteElement<dim> > fe_;
510 };
511 
512 
513 template<unsigned int dim>
514 class JoinValues : public BaseValues<dim>
515 {
516 public:
517  /// Constructor
518  JoinValues(PatchPointValues<3> *patch_point_vals_bulk, PatchPointValues<3> *patch_point_vals_side, MixedPtr<FiniteElement> fe)
519  : BaseValues<dim>(), patch_point_vals_bulk_(patch_point_vals_bulk), patch_point_vals_side_(patch_point_vals_side) {
520  ASSERT_EQ(patch_point_vals_bulk->dim(), dim-1);
521  ASSERT_EQ(patch_point_vals_side->dim(), dim);
522  fe_high_dim_ = fe[Dim<dim>{}];
523  fe_low_dim_ = fe[Dim<dim-1>{}];
524  }
525 
527  {
528  // element of lower dim (bulk points)
529  auto fe_component_low = this->fe_comp(fe_low_dim_, component_idx);
530  ASSERT_EQ(fe_component_low->fe_type(), FEType::FEScalar).error("Type of FiniteElement of scalar_shape accessor must be FEScalar!\n");
531  // use lambda reinit function
532  std::vector< std::vector<double> > shape_values_bulk( patch_point_vals_bulk_->get_quadrature()->size(), vector<double>(fe_component_low->n_dofs()) );
533  auto ref_shape_vals_bulk = this->ref_shape_values_bulk(patch_point_vals_bulk_->get_quadrature(), fe_component_low);
534  for (unsigned int i = 0; i < patch_point_vals_bulk_->get_quadrature()->size(); i++)
535  for (unsigned int j = 0; j < fe_component_low->n_dofs(); j++) {
536  shape_values_bulk[i][j] = ref_shape_vals_bulk[i][j][0];
537  }
538  uint scalar_shape_op_idx_bulk = patch_point_vals_bulk_->operations_.size(); // index in operations_ vector
539  auto lambda_scalar_shape_bulk = [shape_values_bulk, scalar_shape_op_idx_bulk](std::vector<ElOp<3>> &operations, TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
540  bulk_reinit::ptop_scalar_shape(operations, op_results, shape_values_bulk, scalar_shape_op_idx_bulk);
541  };
542  auto &grad_scalar_shape_bulk_op = patch_point_vals_bulk_->make_fe_op({1}, lambda_scalar_shape_bulk, {}, fe_component_low->n_dofs());
543  uint begin_bulk = grad_scalar_shape_bulk_op.result_row();
544 
545  // element of higher dim (side points)
546  auto fe_component_high = this->fe_comp(fe_high_dim_, component_idx);
547  ASSERT_EQ(fe_component_high->fe_type(), FEType::FEScalar).error("Type of FiniteElement of scalar_shape accessor must be FEScalar!\n");
548  // use lambda reinit function
549  std::vector< std::vector< std::vector<double> > > shape_values_side(
550  dim+1,
552  );
553  auto ref_shape_vals_side = this->ref_shape_values_side(patch_point_vals_side_->get_quadrature(), fe_component_high);
554  for (unsigned int s=0; s<dim+1; ++s)
555  for (unsigned int i = 0; i < patch_point_vals_side_->get_quadrature()->size(); i++)
556  for (unsigned int j = 0; j < fe_component_high->n_dofs(); j++) {
557  shape_values_side[s][i][j] = ref_shape_vals_side[s][i][j][0];
558  }
559  uint scalar_shape_op_idx_side = patch_point_vals_side_->operations_.size(); // index in operations_ vector
560  auto lambda_scalar_shape_side = [shape_values_side, scalar_shape_op_idx_side](std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
561  side_reinit::ptop_scalar_shape(operations, op_results, el_table, shape_values_side, scalar_shape_op_idx_side);
562  };
563  auto &grad_scalar_shape_side_op = patch_point_vals_side_->make_fe_op({1}, lambda_scalar_shape_side, {}, fe_component_high->n_dofs());
564  uint begin_side = grad_scalar_shape_side_op.result_row();
565 
566  auto bgn_it = make_iter<JoinShapeAccessor<Scalar>>( JoinShapeAccessor<Scalar>(patch_point_vals_bulk_, patch_point_vals_side_,
567  begin_bulk, begin_side, fe_component_low->n_dofs(), fe_component_high->n_dofs(), 0) );
568  unsigned int end_idx = fe_component_low->n_dofs() + fe_component_high->n_dofs();
569  auto end_it = make_iter<JoinShapeAccessor<Scalar>>( JoinShapeAccessor<Scalar>(patch_point_vals_bulk_, patch_point_vals_side_,
570  begin_bulk, begin_side, fe_component_low->n_dofs(), fe_component_high->n_dofs(), end_idx) );
571  return Range<JoinShapeAccessor<Scalar>>(bgn_it, end_it);
572  }
573 
574 private:
577  std::shared_ptr< FiniteElement<dim> > fe_high_dim_;
578  std::shared_ptr< FiniteElement<dim-1> > fe_low_dim_;
579 };
580 
581 /// Template specialization of dim = 1
582 template <>
583 class JoinValues<1> : public BaseValues<1>
584 {
585 public:
586  /// Constructor
588  : BaseValues<1>() {}
589 
591  {
593  make_iter<JoinShapeAccessor<Scalar>>(JoinShapeAccessor<Scalar>()),
595  }
596 };
597 
598 
599 template<unsigned int spacedim = 3>
601 public:
602  /// Struct for pre-computing number of elements, sides, bulk points and side points on each dimension.
603  struct TableSizes {
604  public:
605  /// Constructor
609  }
610 
611  /// Set all values to zero
612  void reset() {
613  std::fill(elem_sizes_[0].begin(), elem_sizes_[0].end(), 0);
614  std::fill(elem_sizes_[1].begin(), elem_sizes_[1].end(), 0);
615  std::fill(point_sizes_[0].begin(), point_sizes_[0].end(), 0);
616  std::fill(point_sizes_[1].begin(), point_sizes_[1].end(), 0);
617  }
618 
619  /**
620  * Holds number of elements and sides on each dimension
621  * Format:
622  * { {n_elements_1D, n_elements_2D, n_elements_3D },
623  * {n_sides_1D, n_sides_2D, n_sides_3D } }
624  */
626 
627  /**
628  * Holds number of bulk and side points on each dimension
629  * Format:
630  * { {n_bulk_points_1D, n_bulk_points_2D, n_bulk_points_3D },
631  * {n_side_points_1D, n_side_points_2D, n_side_points_3D } }
632  */
634  };
635 
638  patch_point_vals_side_{ {FeSide::PatchPointValues(1, 0), FeSide::PatchPointValues(2, 0), FeSide::PatchPointValues(3, 0)} },
639  asm_arena_(1024 * 1024, 256),
640  patch_arena_(1024 * 1024, 256)
641  {
642  used_quads_[0] = false; used_quads_[1] = false;
643  }
644 
645  PatchFEValues(unsigned int quad_order, MixedPtr<FiniteElement> fe)
646  : patch_point_vals_bulk_{ {FeBulk::PatchPointValues(1, quad_order),
647  FeBulk::PatchPointValues(2, quad_order),
648  FeBulk::PatchPointValues(3, quad_order)} },
649  patch_point_vals_side_{ {FeSide::PatchPointValues(1, quad_order),
650  FeSide::PatchPointValues(2, quad_order),
651  FeSide::PatchPointValues(3, quad_order)} },
652  fe_(fe),
653  asm_arena_(1024 * 1024, 256),
654  patch_arena_(1024 * 1024, 256) {
655  used_quads_[0] = false; used_quads_[1] = false;
656  }
657 
658 
659  /// Destructor
661  {}
662 
663  /// Return bulk or side quadrature of given dimension
664  Quadrature *get_quadrature(uint dim, bool is_bulk) const {
665  if (is_bulk) return patch_point_vals_bulk_[dim-1].get_quadrature();
666  else return patch_point_vals_side_[dim-1].get_quadrature();
667  }
668 
669  /**
670  * @brief Initialize structures and calculates cell-independent data.
671  *
672  * @param _quadrature The quadrature rule for the cell associated
673  * to given finite element or for the cell side.
674  * @param _flags The update flags.
675  */
676  template<unsigned int DIM>
677  void initialize(Quadrature &_quadrature)
678  {
679  if ( _quadrature.dim() == DIM ) {
680  used_quads_[0] = true;
681  patch_point_vals_bulk_[DIM-1].initialize(3); // bulk
682  } else {
683  used_quads_[1] = true;
684  patch_point_vals_side_[DIM-1].initialize(4); // side
685  }
686  }
687 
688  /// Reset PatchpointValues structures
689  void reset()
690  {
691  for (unsigned int i=0; i<3; ++i) {
692  if (used_quads_[0]) patch_point_vals_bulk_[i].reset();
693  if (used_quads_[1]) patch_point_vals_side_[i].reset();
694  }
695  }
696 
697  /// Reinit data.
699  {
700  for (unsigned int i=0; i<3; ++i) {
701  if (used_quads_[0]) patch_point_vals_bulk_[i].reinit_patch();
702  if (used_quads_[1]) patch_point_vals_side_[i].reinit_patch();
703  }
704  }
705 
706  /**
707  * @brief Returns the number of shape functions.
708  */
709  template<unsigned int dim>
710  inline unsigned int n_dofs() const {
711  ASSERT((dim>=0) && (dim<=3))(dim).error("Dimension must be 0, 1, 2 or 3.");
712  return fe_[Dim<dim>{}]->n_dofs();
713  }
714 
715  /// Return BulkValue object of dimension given by template parameter
716  template<unsigned int dim>
718  ASSERT((dim>0) && (dim<=3))(dim).error("Dimension must be 1, 2 or 3.");
719  return BulkValues<dim>(patch_point_vals_bulk_[dim-1], fe_);
720  }
721 
722  /// Return SideValue object of dimension given by template parameter
723  template<unsigned int dim>
725  ASSERT((dim>0) && (dim<=3))(dim).error("Dimension must be 1, 2 or 3.");
726  return SideValues<dim>(patch_point_vals_side_[dim-1], fe_);
727  }
728 
729  /// Return JoinValue object of dimension given by template parameter
730  template<unsigned int dim>
732  //ASSERT((dim>1) && (dim<=3))(dim).error("Dimension must be 2 or 3.");
733  return JoinValues<dim>(&patch_point_vals_bulk_[dim-2], &patch_point_vals_side_[dim-1], fe_);
734  }
735 
736  /** Following methods are used during update of patch. **/
737 
738  /// Resize tables of patch_point_vals_
739  void resize_tables(TableSizes table_sizes) {
740  for (uint i=0; i<3; ++i) {
741  if (used_quads_[0]) patch_point_vals_bulk_[i].resize_tables(table_sizes.elem_sizes_[0][i], table_sizes.point_sizes_[0][i]);
742  if (used_quads_[1]) patch_point_vals_side_[i].resize_tables(table_sizes.elem_sizes_[1][i], table_sizes.point_sizes_[1][i]);
743  }
744  }
745 
746  /// Register element to patch_point_vals_ table by dimension of element
747  uint register_element(DHCellAccessor cell, uint element_patch_idx) {
748  arma::mat coords;
749  switch (cell.dim()) {
750  case 1:
751  coords = MappingP1<1,spacedim>::element_map(cell.elm());
752  return patch_point_vals_bulk_[0].register_element(coords, element_patch_idx);
753  break;
754  case 2:
755  coords = MappingP1<2,spacedim>::element_map(cell.elm());
756  return patch_point_vals_bulk_[1].register_element(coords, element_patch_idx);
757  break;
758  case 3:
759  coords = MappingP1<3,spacedim>::element_map(cell.elm());
760  return patch_point_vals_bulk_[2].register_element(coords, element_patch_idx);
761  break;
762  default:
763  ASSERT(false);
764  return 0;
765  break;
766  }
767  }
768 
769  /// Register side to patch_point_vals_ table by dimension of side
771  arma::mat side_coords(spacedim, cell_side.dim());
772  for (unsigned int n=0; n<cell_side.dim(); n++)
773  for (unsigned int c=0; c<spacedim; c++)
774  side_coords(c,n) = (*cell_side.side().node(n))[c];
775 
776  arma::mat elm_coords;
777  DHCellAccessor cell = cell_side.cell();
778  switch (cell.dim()) {
779  case 1:
780  elm_coords = MappingP1<1,spacedim>::element_map(cell.elm());
781  return patch_point_vals_side_[0].register_side(elm_coords, side_coords);
782  break;
783  case 2:
784  elm_coords = MappingP1<2,spacedim>::element_map(cell.elm());
785  return patch_point_vals_side_[1].register_side(elm_coords, side_coords);
786  break;
787  case 3:
788  elm_coords = MappingP1<3,spacedim>::element_map(cell.elm());
789  return patch_point_vals_side_[2].register_side(elm_coords, side_coords);
790  break;
791  default:
792  ASSERT(false);
793  return 0;
794  break;
795  }
796  }
797 
798  /// Register bulk point to patch_point_vals_ table by dimension of element
799  uint register_bulk_point(DHCellAccessor cell, uint elem_table_row, uint value_patch_idx) {
800  return patch_point_vals_bulk_[cell.dim()-1].register_bulk_point(elem_table_row, value_patch_idx, cell.elm_idx());
801  }
802 
803  /// Register side point to patch_point_vals_ table by dimension of side
804  uint register_side_point(DHCellSide cell_side, uint elem_table_row, uint value_patch_idx) {
805  return patch_point_vals_side_[cell_side.dim()-1].register_side_point(elem_table_row, value_patch_idx, cell_side.elem_idx(), cell_side.side_idx());
806  }
807 
808  /// Temporary development method
809  void print_data_tables(ostream& stream, bool points, bool ints, bool only_bulk=true) const {
810  stream << endl << "Table of patch FE data:" << endl;
811  for (uint i=0; i<3; ++i) {
812  stream << std::setfill('-') << setw(100) << "" << endl;
813  stream << "Bulk, dimension " << (i+1) << endl;
814  patch_point_vals_bulk_[i].print_data_tables(stream, points, ints);
815  }
816  if (!only_bulk)
817  for (uint i=0; i<3; ++i) {
818  stream << std::setfill('-') << setw(100) << "" << endl;
819  stream << "Side, dimension " << (i+1) << endl;
820  patch_point_vals_side_[i].print_data_tables(stream, points, ints);
821  }
822  stream << std::setfill('=') << setw(100) << "" << endl;
823  }
824 
825  /// Temporary development method
826  void print_operations(ostream& stream) const {
827  stream << endl << "Table of patch FE operations:" << endl;
828  for (uint i=0; i<3; ++i) {
829  stream << std::setfill('-') << setw(100) << "" << endl;
830  stream << "Bulk, dimension " << (i+1) << ", n_rows " << patch_point_vals_bulk_[i].n_rows() << endl;
831  patch_point_vals_bulk_[i].print_operations(stream, 0);
832  }
833  for (uint i=0; i<3; ++i) {
834  stream << std::setfill('-') << setw(100) << "" << endl;
835  stream << "Side, dimension " << (i+1) << ", n_rows " << patch_point_vals_side_[i].n_rows() << endl;
836  patch_point_vals_side_[i].print_operations(stream, 1);
837  }
838  stream << std::setfill('=') << setw(100) << "" << endl;
839  }
840 
841 private:
842  std::array<FeBulk::PatchPointValues<spacedim>, 3> patch_point_vals_bulk_; ///< Sub objects of bulk data of dimensions 1,2,3
843  std::array<FeSide::PatchPointValues<spacedim>, 3> patch_point_vals_side_; ///< Sub objects of side data of dimensions 1,2,3
844 
845  MixedPtr<FiniteElement> fe_; ///< Mixed of shared pointers of FiniteElement object
846  bool used_quads_[2]; ///< Pair of flags signs holds info if bulk and side quadratures are used
849 
850  template <class ValueType>
851  friend class ElQ;
852  template <class ValueType>
853  friend class FeQ;
854  template <class ValueType>
855  friend class JoinShapeAccessor;
856 };
857 
858 
859 template <class ValueType>
860 ValueType ElQ<ValueType>::operator()(const BulkPoint &point) {
861  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
862  return patch_point_vals_.scalar_val(begin_, value_cache_idx);
863 }
864 
865 template <>
867  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
868  return patch_point_vals_.vector_val(begin_, value_cache_idx);
869 }
870 
871 template <>
873  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
874  return patch_point_vals_.tensor_val(begin_, value_cache_idx);
875 }
876 
877 template <class ValueType>
878 ValueType ElQ<ValueType>::operator()(const SidePoint &point) {
879  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
880  return patch_point_vals_.scalar_val(begin_, value_cache_idx);
881 }
882 
883 template <>
885  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
886  return patch_point_vals_.vector_val(begin_, value_cache_idx);
887 }
888 
889 template <>
891  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
892  return patch_point_vals_.tensor_val(begin_, value_cache_idx);
893 }
894 
895 template <class ValueType>
896 ValueType FeQ<ValueType>::operator()(unsigned int shape_idx, const BulkPoint &point) {
897  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
898  return patch_point_vals_.scalar_val(begin_+shape_idx, value_cache_idx);
899 }
900 
901 template <>
902 inline Vector FeQ<Vector>::operator()(unsigned int shape_idx, const BulkPoint &point) {
903  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
904  return patch_point_vals_.vector_val(begin_+3*shape_idx, value_cache_idx);
905 }
906 
907 template <>
908 inline Tensor FeQ<Tensor>::operator()(FMT_UNUSED unsigned int shape_idx, FMT_UNUSED const BulkPoint &point) {
909  Tensor tens; tens.zeros();
910  return tens;
911 }
912 
913 template <class ValueType>
914 ValueType FeQ<ValueType>::operator()(unsigned int shape_idx, const SidePoint &point) {
915  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
916  return patch_point_vals_.scalar_val(begin_+shape_idx, value_cache_idx);
917 }
918 
919 template <>
920 inline Vector FeQ<Vector>::operator()(unsigned int shape_idx, const SidePoint &point) {
921  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
922  return patch_point_vals_.vector_val(begin_+3*shape_idx, value_cache_idx);
923 }
924 
925 template <>
926 inline Tensor FeQ<Tensor>::operator()(FMT_UNUSED unsigned int shape_idx, FMT_UNUSED const SidePoint &point) {
927  Tensor tens; tens.zeros();
928  return tens;
929 }
930 
931 
932 template <class ValueType>
934  if (this->is_high_dim()) {
935  return 0.0;
936  } else {
937  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
938  return patch_point_vals_bulk_->scalar_val(begin_+this->local_idx(), value_cache_idx);
939  }
940 }
941 
942 template <>
944  Vector vect; vect.zeros();
945  return vect;
946 }
947 
948 template <>
950  Tensor tens; tens.zeros();
951  return tens;
952 }
953 
954 template <class ValueType>
956  if (this->is_high_dim()) {
957  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
958  return patch_point_vals_side_->scalar_val(begin_side_+this->local_idx(), value_cache_idx);
959  } else {
960  return 0.0;
961  }
962 }
963 
964 template <>
966  Vector vect; vect.zeros();
967  return vect;
968 }
969 
970 template <>
972  Tensor tens; tens.zeros();
973  return tens;
974 }
975 
976 
977 #endif /* PATCH_FE_VALUES_HH_ */
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
std::vector< std::vector< std::vector< arma::vec > > > ref_shape_values_side(Quadrature *q, std::shared_ptr< FiniteElement< FE_dim >> fe)
Precomputed values of basis functions at the side quadrature points.
std::shared_ptr< FiniteElement< FE_dim > > fe_comp(std::shared_ptr< FiniteElement< FE_dim > > fe, uint component_idx)
Return FiniteElement of component_idx for FESystem or fe for other types.
std::vector< std::vector< arma::vec > > ref_shape_values_bulk(Quadrature *q, std::shared_ptr< FiniteElement< FE_dim >> fe)
Precomputed values of basis functions at the bulk quadrature points.
Base point accessor class.
Definition: eval_subset.hh:55
const ElementCacheMap * elm_cache_map() const
Definition: eval_subset.hh:73
unsigned int elem_patch_idx() const
Definition: eval_subset.hh:79
unsigned int eval_point_idx() const
Return index in EvalPoints object.
Definition: eval_subset.hh:84
ElQ< Scalar > JxW()
Register the product of Jacobian determinant and the quadrature weight at bulk quadrature points.
PatchPointValues< 3 > & patch_point_vals_
ElQ< Vector > coords()
Create bulk accessor of coords entity.
std::vector< std::vector< arma::mat > > ref_shape_gradients(std::shared_ptr< FiniteElement< dim >> fe)
Precomputed gradients of basis functions at the quadrature points.
ElQ< Scalar > determinant()
Create bulk accessor of jac determinant entity.
std::shared_ptr< FiniteElement< dim > > fe_
FeQ< Vector > grad_scalar_shape(uint component_idx=0)
Return the value of the function_no-th gradient shape function at the p bulk quadrature point.
FeQ< Scalar > scalar_shape(uint component_idx=0)
Return the value of the function_no-th shape function at the p bulk quadrature point.
BulkValues(PatchPointValues< 3 > &patch_point_vals, MixedPtr< FiniteElement > fe)
Constructor.
Cell accessor allow iterate over DOF handler cells.
unsigned int dim() const
Return dimension of element appropriate to cell.
unsigned int elm_idx() const
Return serial idx to element of loc_ele_idx_.
ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
Side accessor allows to iterate over sides of DOF handler cell.
unsigned int elem_idx() const
Side side() const
Return Side of given cell and side_idx.
const DHCellAccessor & cell() const
Return DHCellAccessor appropriate to the side.
unsigned int dim() const
Return dimension of element appropriate to the side.
unsigned int side_idx() const
ElQ()=delete
Forbidden default constructor.
PatchPointValues< 3 > & patch_point_vals_
Reference to PatchPointValues.
unsigned int begin_
ValueType operator()(FMT_UNUSED const SidePoint &point)
ValueType operator()(FMT_UNUSED const BulkPoint &point)
ElQ(PatchPointValues< 3 > &patch_point_vals, unsigned int begin)
Constructor.
int element_eval_point(unsigned int i_elem_in_cache, unsigned int i_eval_point) const
Compound finite element on dim dimensional simplex.
Definition: fe_system.hh:102
const std::vector< std::shared_ptr< FiniteElement< dim > > > & fe() const
Definition: fe_system.hh:147
Bulk data specialization, order of item in operations_ vector corresponds to the BulkOps enum.
unsigned int begin_
Index of the first component of the Quantity. Size is given by ValueType.
PatchPointValues< 3 > & patch_point_vals_
Reference to PatchPointValues.
unsigned int n_dofs_
Number of DOFs.
ValueType operator()(FMT_UNUSED unsigned int shape_idx, FMT_UNUSED const BulkPoint &point)
FeQ(PatchPointValues< 3 > &patch_point_vals, unsigned int begin, unsigned int n_dofs)
FeQ()=delete
Forbidden default constructor.
ValueType operator()(FMT_UNUSED unsigned int shape_idx, FMT_UNUSED const SidePoint &point)
Bulk Side specialization, order of item in operations_ vector corresponds to the SideOps enum.
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
unsigned int n_dofs_both() const
unsigned int begin_side_
Index of the first component of the side Quantity. Size is given by ValueType.
unsigned int local_idx() const
Return local index of DOF (on low / high-dim) - should be private method.
bool operator==(const JoinShapeAccessor< ValueType > &other) const
Comparison of accessors.
bool is_high_dim() const
unsigned int begin_
Index of the first component of the bulk Quantity. Size is given by ValueType.
PatchPointValues< 3 > * patch_point_vals_bulk_
Pointer to bulk PatchPointValues.
JoinShapeAccessor(PatchPointValues< 3 > *patch_point_vals_bulk, PatchPointValues< 3 > *patch_point_vals_side, unsigned int begin, unsigned int begin_side, unsigned int n_dofs_bulk, unsigned int n_dofs_side, unsigned int join_idx)
unsigned int n_dofs_high() const
JoinShapeAccessor()
Default constructor.
void inc()
Iterates to next item.
unsigned int join_idx_
Index of processed DOF.
unsigned int n_dofs_high_
Number of DOFs on high-dim element.
unsigned int join_idx() const
Return global index of DOF.
unsigned int n_dofs_low() const
ValueType operator()(const BulkPoint &point)
unsigned int n_dofs_low_
Number of DOFs on low-dim element.
PatchPointValues< 3 > * patch_point_vals_side_
Pointer to side PatchPointValues.
JoinValues(FMT_UNUSED PatchPointValues< 3 > *patch_point_vals_bulk, FMT_UNUSED PatchPointValues< 3 > *patch_point_vals_side, FMT_UNUSED MixedPtr< FiniteElement > fe)
Constructor.
Range< JoinShapeAccessor< Scalar > > scalar_join_shape(FMT_UNUSED uint component_idx=0)
std::shared_ptr< FiniteElement< dim > > fe_high_dim_
JoinValues(PatchPointValues< 3 > *patch_point_vals_bulk, PatchPointValues< 3 > *patch_point_vals_side, MixedPtr< FiniteElement > fe)
Constructor.
PatchPointValues< 3 > * patch_point_vals_bulk_
Range< JoinShapeAccessor< Scalar > > scalar_join_shape(uint component_idx=0)
std::shared_ptr< FiniteElement< dim-1 > > fe_low_dim_
PatchPointValues< 3 > * patch_point_vals_side_
static ElementMap element_map(ElementAccessor< 3 > elm)
Definition: mapping_p1.cc:48
uint register_element(DHCellAccessor cell, uint element_patch_idx)
Register element to patch_point_vals_ table by dimension of element.
~PatchFEValues()
Destructor.
std::array< FeBulk::PatchPointValues< spacedim >, 3 > patch_point_vals_bulk_
Sub objects of bulk data of dimensions 1,2,3.
uint register_bulk_point(DHCellAccessor cell, uint elem_table_row, uint value_patch_idx)
Register bulk point to patch_point_vals_ table by dimension of element.
uint register_side(DHCellSide cell_side)
Register side to patch_point_vals_ table by dimension of side.
void print_data_tables(ostream &stream, bool points, bool ints, bool only_bulk=true) const
Temporary development method.
JoinValues< dim > join_values()
Return JoinValue object of dimension given by template parameter.
Quadrature * get_quadrature(uint dim, bool is_bulk) const
Return bulk or side quadrature of given dimension.
BulkValues< dim > bulk_values()
Return BulkValue object of dimension given by template parameter.
void initialize(Quadrature &_quadrature)
Initialize structures and calculates cell-independent data.
void print_operations(ostream &stream) const
Temporary development method.
std::array< FeSide::PatchPointValues< spacedim >, 3 > patch_point_vals_side_
Sub objects of side data of dimensions 1,2,3.
MixedPtr< FiniteElement > fe_
Mixed of shared pointers of FiniteElement object.
SideValues< dim > side_values()
Return SideValue object of dimension given by template parameter.
AssemblyArena asm_arena_
AssemblyArena patch_arena_
void resize_tables(TableSizes table_sizes)
Resize tables of patch_point_vals_.
PatchFEValues(unsigned int quad_order, MixedPtr< FiniteElement > fe)
void reinit_patch()
Reinit data.
unsigned int n_dofs() const
Returns the number of shape functions.
uint register_side_point(DHCellSide cell_side, uint elem_table_row, uint value_patch_idx)
Register side point to patch_point_vals_ table by dimension of side.
void reset()
Reset PatchpointValues structures.
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 dim() const
Getter for dim_.
Scalar scalar_val(uint result_row, uint point_idx) const
Base class for quadrature rules on simplices in arbitrary dimensions.
Definition: quadrature.hh:48
Quadrature make_from_side(unsigned int sid) const
Definition: quadrature.cc:47
unsigned int size() const
Returns number of quadrature points.
Definition: quadrature.hh:86
unsigned int dim() const
Definition: quadrature.hh:72
Armor::ArmaVec< double, point_dim > point(unsigned int i) const
Returns the ith quadrature point.
Definition: quadrature.hh:151
Range helper class.
General point a+ side_begin_ + ccessor allow iterate over quadrature points of given side defined in ...
Definition: eval_subset.hh:116
unsigned int eval_point_idx() const
Return index in EvalPoints object.
Definition: eval_subset.hh:143
SideValues(PatchPointValues< 3 > &patch_point_vals, MixedPtr< FiniteElement > fe)
Constructor.
ElQ< Vector > normal_vector()
Register the normal vector to a side at side quadrature points.
std::shared_ptr< FiniteElement< dim > > fe_
ElQ< Scalar > JxW()
Same as BulkValues::JxW but register at side quadrature points.
ElQ< Vector > coords()
Create side accessor of coords entity.
FeQ< Scalar > scalar_shape(uint component_idx=0)
Same as BulkValues::scalar_shape but register at side quadrature points.
FeQ< Vector > grad_scalar_shape(uint component_idx=0)
Same as BulkValues::grad_scalar_shape but register at side quadrature points.
std::vector< std::vector< std::vector< arma::mat > > > ref_shape_gradients(std::shared_ptr< FiniteElement< dim >> fe)
Precomputed gradients of basis functions at the quadrature points.
PatchPointValues< 3 > & patch_point_vals_
ElQ< Scalar > determinant()
Create bulk accessor of jac determinant entity.
NodeAccessor< 3 > node(unsigned int i) const
Returns node for given local index i on the side.
Store finite element data on the actual patch such as shape function values, gradients,...
Eigen::Vector< ArrayDbl, Eigen::Dynamic > TableDbl
Definition: eigen_tools.hh:33
Eigen::Vector< ArrayInt, Eigen::Dynamic > TableInt
Definition: eigen_tools.hh:34
Class ElementValues calculates data related to transformation of reference cell to actual cell (Jacob...
Class FESystem for compound finite elements.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
@ FEScalar
Iter< Object > make_iter(Object obj)
Class MappingP1 implements the affine transformation of the unit cell onto the actual cell.
unsigned int uint
ArmaMat< double, N, M > mat
Definition: armor.hh:936
@ opCoords
operations evaluated on quadrature points
@ opInvJac
inverse Jacobian
@ opJxW
JxW value of quadrature point.
@ opNormalVec
normal vector of quadrature point
std::vector< std::array< uint, 3 > > DimPointTable
Holds triplet (dim; bulk/side; idx of point in subtable)
Store finite element data on the actual patch such as shape function values, gradients,...
#define FMT_UNUSED
Definition: posix.h:75
Definitions of particular quadrature rules on simplices.
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.
Definition: mixed.hh:25
Struct for pre-computing number of elements, sides, bulk points and side points on each dimension.
std::vector< std::vector< uint > > point_sizes_
void reset()
Set all values to zero.
std::vector< std::vector< uint > > elem_sizes_
static void ptop_scalar_shape(std::vector< ElOp< 3 >> &operations, TableDbl &op_results, std::vector< std::vector< double > > shape_values, uint scalar_shape_op_idx)
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)
Enum type UpdateFlags indicates which quantities are to be recomputed on each finite element cell.