Flow123d  DF_patch_fe_data_tables-5938e3d
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, unsigned int op_idx)
59  : patch_point_vals_(patch_point_vals), begin_(begin), op_idx_(op_idx) {}
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  unsigned int op_idx_; /// Index of operation in patch_point_vals_.operations vector
69 };
70 
71 
72 template <class ValueType>
73 class FeQ {
74 public:
75  /// Forbidden default constructor
76  FeQ() = delete;
77 
78  // Class similar to current FeView
79  FeQ(PatchPointValues<3> &patch_point_vals, unsigned int begin, unsigned int n_dofs)
80  : patch_point_vals_(patch_point_vals), begin_(begin), n_dofs_(n_dofs) {}
81 
82 
83  ValueType operator()(FMT_UNUSED unsigned int shape_idx, FMT_UNUSED const BulkPoint &point);
84 
85  ValueType operator()(FMT_UNUSED unsigned int shape_idx, FMT_UNUSED const SidePoint &point);
86 
87  // Implementation for EdgePoint, SidePoint, and JoinPoint shoud have a common implementation
88  // resolving to side values
89 
90 private:
91  PatchPointValues<3> &patch_point_vals_; ///< Reference to PatchPointValues
92  unsigned int begin_; ///< Index of the first component of the Quantity. Size is given by ValueType
93  unsigned int n_dofs_; ///< Number of DOFs
94 };
95 
96 
97 template <class ValueType>
99 public:
100  /// Default constructor
102  : patch_point_vals_bulk_(nullptr), patch_point_vals_side_(nullptr), join_idx_(-1) {}
103 
104  /**
105  * Constructor
106  *
107  * @param patch_point_vals_bulk Pointer to PatchPointValues bulk object.
108  * @param patch_point_vals_side Pointer to PatchPointValues side object.
109  * @param begin Index of the first component of the bulk Quantity.
110  * @param begin_side Index of the first component of the side Quantity.
111  * @param n_dofs_bulk Number of DOFs of bulk (lower-dim) element.
112  * @param n_dofs_side Number of DOFs of side (higher-dim) element.
113  * @param join_idx Index function.
114  */
115  JoinShapeAccessor(PatchPointValues<3> *patch_point_vals_bulk, PatchPointValues<3> *patch_point_vals_side,
116  unsigned int begin, unsigned int begin_side, unsigned int n_dofs_bulk, unsigned int n_dofs_side, unsigned int join_idx)
117  : patch_point_vals_bulk_(patch_point_vals_bulk), patch_point_vals_side_(patch_point_vals_side), begin_(begin),
118  begin_side_(begin_side), n_dofs_high_(n_dofs_side), n_dofs_low_(n_dofs_bulk), join_idx_(join_idx) {
119  //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!");
120  }
121 
122  /// Return global index of DOF
123  inline unsigned int join_idx() const {
124  return join_idx_;
125  }
126 
127  /// Return local index of DOF (on low / high-dim) - should be private method
128  inline unsigned int local_idx() const {
129  if (this->is_high_dim()) return (join_idx_ - n_dofs_low_);
130  else return join_idx_;
131  }
132 
133  inline unsigned int n_dofs_low() const {
134  return n_dofs_low_;
135  }
136 
137  inline unsigned int n_dofs_high() const {
138  return n_dofs_high_;
139  }
140 
141  inline unsigned int n_dofs_both() const {
142  return n_dofs_high_ + n_dofs_low_;
143  }
144 
145  inline bool is_high_dim() const {
146  return (join_idx_ >= n_dofs_low_);
147  }
148 
149  /// Iterates to next item.
150  inline void inc() {
151  join_idx_++;
152  }
153 
154  /// Comparison of accessors.
155  bool operator==(const JoinShapeAccessor<ValueType>& other) const {
156  return (join_idx_ == other.join_idx_);
157  }
158 
159 
160  ValueType operator()(const BulkPoint &point);
161 
162  ValueType operator()(const SidePoint &point);
163 
164 private:
165  // attributes:
166  PatchPointValues<3> *patch_point_vals_bulk_; ///< Pointer to bulk PatchPointValues
167  PatchPointValues<3> *patch_point_vals_side_; ///< Pointer to side PatchPointValues
168  unsigned int begin_; ///< Index of the first component of the bulk Quantity. Size is given by ValueType
169  unsigned int begin_side_; ///< Index of the first component of the side Quantity. Size is given by ValueType
170  unsigned int n_dofs_high_; ///< Number of DOFs on high-dim element
171  unsigned int n_dofs_low_; ///< Number of DOFs on low-dim element
172  unsigned int join_idx_; ///< Index of processed DOF
173 };
174 
175 
176 template<unsigned int dim>
178 {
179 protected:
180  // Default constructor
182  {}
183 
184  /// Return FiniteElement of \p component_idx for FESystem or \p fe for other types
185  template<unsigned int FE_dim>
186  std::shared_ptr<FiniteElement<FE_dim>> fe_comp(std::shared_ptr< FiniteElement<FE_dim> > fe, uint component_idx) {
187  FESystem<FE_dim> *fe_sys = dynamic_cast<FESystem<FE_dim>*>( fe.get() );
188  if (fe_sys != nullptr) {
189  return fe_sys->fe()[component_idx];
190  } else {
191  ASSERT_EQ(component_idx, 0).warning("Non-zero component_idx can only be used for FESystem.");
192  return fe;
193  }
194  }
195 
196  /**
197  * @brief Precomputed values of basis functions at the bulk quadrature points.
198  *
199  * Dimensions: (no. of quadrature points)
200  * x (no. of dofs)
201  * x (no. of components in ref. cell)
202  */
203  template<unsigned int FE_dim>
205  std::vector<std::vector<arma::vec> > ref_shape_vals( q->size(), vector<arma::vec>(fe->n_dofs()) );
206 
207  arma::mat shape_values(fe->n_dofs(), fe->n_components());
208  for (unsigned int i=0; i<q->size(); i++)
209  {
210  for (unsigned int j=0; j<fe->n_dofs(); j++)
211  {
212  for (unsigned int c=0; c<fe->n_components(); c++)
213  shape_values(j,c) = fe->shape_value(j, q->point<FE_dim>(i), c);
214 
215  ref_shape_vals[i][j] = trans(shape_values.row(j));
216  }
217  }
218 
219  return ref_shape_vals;
220  }
221 
222  /**
223  * @brief Precomputed values of basis functions at the side quadrature points.
224  *
225  * Dimensions: (sides)
226  * x (no. of quadrature points)
227  * x (no. of dofs)
228  * x (no. of components in ref. cell)
229  */
230  template<unsigned int FE_dim>
233 
234  arma::mat shape_values(fe->n_dofs(), fe->n_components());
235 
236  for (unsigned int sid=0; sid<FE_dim+1; sid++) {
237  auto quad = q->make_from_side<dim>(sid);
238  for (unsigned int i=0; i<quad.size(); i++)
239  {
240  for (unsigned int j=0; j<fe->n_dofs(); j++)
241  {
242  for (unsigned int c=0; c<fe->n_components(); c++) {
243  shape_values(j,c) = fe->shape_value(j, quad.template point<FE_dim>(i), c);
244  }
245 
246  ref_shape_vals[sid][i][j] = trans(shape_values.row(j));
247  }
248  }
249  }
250 
251  return ref_shape_vals;
252  }
253 
254 };
255 
256 template<unsigned int dim>
257 class BulkValues : public BaseValues<dim>
258 {
259 public:
260  /// Constructor
262  : BaseValues<dim>(), patch_point_vals_(patch_point_vals) {
263  ASSERT_EQ(patch_point_vals.dim(), dim);
264  fe_ = fe[Dim<dim>{}];
265  }
266 
267  /**
268  * @brief Register the product of Jacobian determinant and the quadrature
269  * weight at bulk quadrature points.
270  *
271  * @param quad Quadrature.
272  */
273  inline ElQ<Scalar> JxW()
274  {
277  }
278 
279  /// Create bulk accessor of coords entity
281  {
284  }
285 
286 // inline ElQ<Tensor> jacobian(std::initializer_list<Quadrature *> quad_list)
287 // {}
288 
289  /// Create bulk accessor of jac determinant entity
291  {
294  }
295 
296  /**
297  * @brief Return the value of the @p function_no-th shape function at
298  * the @p p bulk quadrature point.
299  *
300  * @param component_idx Number of the shape function.
301  */
302  inline FeQ<Scalar> scalar_shape(uint component_idx = 0)
303  {
304  auto fe_component = this->fe_comp(fe_, component_idx);
305  ASSERT_EQ(fe_component->fe_type(), FEType::FEScalar).error("Type of FiniteElement of scalar_shape accessor must be FEScalar!\n");
306 
307  // use lambda reinit function
308  std::vector< std::vector<double> > shape_values( patch_point_vals_.get_quadrature()->size(), vector<double>(fe_component->n_dofs()) );
309  auto ref_shape_vals = this->ref_shape_values_bulk(patch_point_vals_.get_quadrature(), fe_component);
310  for (unsigned int i = 0; i < patch_point_vals_.get_quadrature()->size(); i++)
311  for (unsigned int j = 0; j < fe_component->n_dofs(); j++) {
312  shape_values[i][j] = ref_shape_vals[i][j][0];
313  }
314  uint scalar_shape_op_idx = patch_point_vals_.operations_.size(); // index in operations_ vector
315  auto lambda_scalar_shape = [shape_values, scalar_shape_op_idx](std::vector<ElOp<3>> &operations, TableDbl &op_results, FMT_UNUSED TableInt &el_table) {
316  bulk_reinit::ptop_scalar_shape(operations, op_results, shape_values, scalar_shape_op_idx);
317  };
318  auto &scalar_shape_bulk_op = patch_point_vals_.make_fe_op({1}, lambda_scalar_shape, {}, fe_component->n_dofs());
319  uint begin = scalar_shape_bulk_op.result_row();
320 
321  return FeQ<Scalar>(patch_point_vals_, begin, fe_component->n_dofs());
322  }
323 
324 // inline FeQ<Vector> vector_shape(uint component_idx = 0)
325 // {}
326 
327 // inline FeQ<Tensor> tensor_shape(uint component_idx = 0)
328 // {}
329 
330  /**
331  * @brief Return the value of the @p function_no-th gradient shape function at
332  * the @p p bulk quadrature point.
333  *
334  * @param component_idx Number of the shape function.
335  */
336  inline FeQ<Vector> grad_scalar_shape(uint component_idx=0)
337  {
338  auto fe_component = this->fe_comp(fe_, component_idx);
339  ASSERT_EQ(fe_component->fe_type(), FEType::FEScalar).error("Type of FiniteElement of grad_scalar_shape accessor must be FEScalar!\n");
340 
341  // use lambda reinit function
342  auto ref_shape_grads = this->ref_shape_gradients(fe_component);
343  uint scalar_shape_grads_op_idx = patch_point_vals_.operations_.size(); // index in operations_ vector
344  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) {
345  bulk_reinit::ptop_scalar_shape_grads<dim>(operations, op_results, ref_shape_grads, scalar_shape_grads_op_idx);
346  };
347  auto &grad_scalar_shape_bulk_op = patch_point_vals_.make_fe_op({3}, lambda_scalar_shape_grad, {FeBulk::BulkOps::opInvJac}, fe_component->n_dofs());
348  uint begin = grad_scalar_shape_bulk_op.result_row();
349 
350  return FeQ<Vector>(patch_point_vals_, begin, fe_component->n_dofs());
351  }
352 
353 // inline FeQ<Tensor> grad_vector_shape(std::initializer_list<Quadrature *> quad_list, unsigned int i_comp=0)
354 // {}
355 
356 private:
357  /**
358  * @brief Precomputed gradients of basis functions at the quadrature points.
359  *
360  * Dimensions: (no. of quadrature points)
361  * x (no. of dofs)
362  * x ((dim_ of. ref. cell)x(no. of components in ref. cell))
363  */
366  std::vector<std::vector<arma::mat> > ref_shape_grads( q->size(), vector<arma::mat>(fe->n_dofs()) );
367 
368  arma::mat grad(dim, fe->n_components());
369  for (unsigned int i_pt=0; i_pt<q->size(); i_pt++)
370  {
371  for (unsigned int i_dof=0; i_dof<fe->n_dofs(); i_dof++)
372  {
373  grad.zeros();
374  for (unsigned int c=0; c<fe->n_components(); c++)
375  grad.col(c) += fe->shape_grad(i_dof, q->point<dim>(i_pt), c);
376 
377  ref_shape_grads[i_pt][i_dof] = grad;
378  }
379  }
380 
381  return ref_shape_grads;
382  }
383 
385  std::shared_ptr< FiniteElement<dim> > fe_;
386 };
387 
388 
389 template<unsigned int dim>
390 class SideValues : public BaseValues<dim>
391 {
392 public:
393  /// Constructor
395  : BaseValues<dim>(), patch_point_vals_(patch_point_vals) {
396  ASSERT_EQ(patch_point_vals.dim(), dim);
397  fe_ = fe[Dim<dim>{}];
398  }
399 
400  /// Same as BulkValues::JxW but register at side quadrature points.
401  inline ElQ<Scalar> JxW()
402  {
405  }
406 
407  /**
408  * @brief Register the normal vector to a side at side quadrature points.
409  *
410  * @param quad Quadrature.
411  */
413  {
416  }
417 
418  /// Create side accessor of coords entity
420  {
423  }
424 
425  /// Create bulk accessor of jac determinant entity
427  {
430  }
431 
432  /// Same as BulkValues::scalar_shape but register at side quadrature points.
433  inline FeQ<Scalar> scalar_shape(uint component_idx = 0)
434  {
435  auto fe_component = this->fe_comp(fe_, component_idx);
436  ASSERT_EQ(fe_component->fe_type(), FEType::FEScalar).error("Type of FiniteElement of scalar_shape accessor must be FEScalar!\n");
437 
438  // use lambda reinit function
440  dim+1,
442  );
443  auto ref_shape_vals = this->ref_shape_values_side(patch_point_vals_.get_quadrature(), fe_component);
444  for (unsigned int s=0; s<dim+1; ++s)
445  for (unsigned int i = 0; i < patch_point_vals_.get_quadrature()->size(); i++)
446  for (unsigned int j = 0; j < fe_component->n_dofs(); j++) {
447  shape_values[s][i][j] = ref_shape_vals[s][i][j][0];
448  }
449  uint scalar_shape_op_idx = patch_point_vals_.operations_.size(); // index in operations_ vector
450  auto lambda_scalar_shape = [shape_values, scalar_shape_op_idx](std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
451  side_reinit::ptop_scalar_shape(operations, op_results, el_table, shape_values, scalar_shape_op_idx);
452  };
453  auto &scalar_shape_bulk_op = patch_point_vals_.make_fe_op({1}, lambda_scalar_shape, {}, fe_component->n_dofs());
454  uint begin = scalar_shape_bulk_op.result_row();
455 
456  return FeQ<Scalar>(patch_point_vals_, begin, fe_component->n_dofs());
457  }
458 
459  /// Same as BulkValues::grad_scalar_shape but register at side quadrature points.
460  inline FeQ<Vector> grad_scalar_shape(uint component_idx=0)
461  {
462  auto fe_component = this->fe_comp(fe_, component_idx);
463  ASSERT_EQ(fe_component->fe_type(), FEType::FEScalar).error("Type of FiniteElement of grad_scalar_shape accessor must be FEScalar!\n");
464 
465  // use lambda reinit function
466  auto ref_shape_grads = this->ref_shape_gradients(fe_component);
467  uint scalar_shape_grads_op_idx = patch_point_vals_.operations_.size(); // index in operations_ vector
468  auto lambda_scalar_shape_grad = [ref_shape_grads, scalar_shape_grads_op_idx](std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
469  side_reinit::ptop_scalar_shape_grads<dim>(operations, op_results, el_table, ref_shape_grads, scalar_shape_grads_op_idx);
470  };
471  auto &grad_scalar_shape_side_op = patch_point_vals_.make_fe_op({3}, lambda_scalar_shape_grad, {FeSide::SideOps::opElInvJac}, fe_component->n_dofs());
472  uint begin = grad_scalar_shape_side_op.result_row();
473 
474  return FeQ<Vector>(patch_point_vals_, begin, fe_component->n_dofs());
475  }
476 
477 private:
478  /**
479  * @brief Precomputed gradients of basis functions at the quadrature points.
480  *
481  * Dimensions: (sides)
482  * x (no. of quadrature points)
483  * x (no. of dofs)
484  * x ((dim_ of. ref. cell)x(no. of components in ref. cell))
485  */
489 
490  arma::mat grad(dim, fe->n_components());
491  for (unsigned int sid=0; sid<dim+1; sid++) {
492  auto quad = q->make_from_side<dim>(sid);
493  for (unsigned int i_pt=0; i_pt<quad.size(); i_pt++)
494  {
495  for (unsigned int i_dof=0; i_dof<fe->n_dofs(); i_dof++)
496  {
497  grad.zeros();
498  for (unsigned int c=0; c<fe->n_components(); c++)
499  grad.col(c) += fe->shape_grad(i_dof, quad.template point<dim>(i_pt), c);
500 
501  ref_shape_grads[sid][i_pt][i_dof] = grad;
502  }
503  }
504  }
505 
506  return ref_shape_grads;
507  }
508 
510  std::shared_ptr< FiniteElement<dim> > fe_;
511 };
512 
513 
514 template<unsigned int dim>
515 class JoinValues : public BaseValues<dim>
516 {
517 public:
518  /// Constructor
519  JoinValues(PatchPointValues<3> *patch_point_vals_bulk, PatchPointValues<3> *patch_point_vals_side, MixedPtr<FiniteElement> fe)
520  : BaseValues<dim>(), patch_point_vals_bulk_(patch_point_vals_bulk), patch_point_vals_side_(patch_point_vals_side) {
521  ASSERT_EQ(patch_point_vals_bulk->dim(), dim-1);
522  ASSERT_EQ(patch_point_vals_side->dim(), dim);
523  fe_high_dim_ = fe[Dim<dim>{}];
524  fe_low_dim_ = fe[Dim<dim-1>{}];
525  }
526 
528  {
529  // element of lower dim (bulk points)
530  auto fe_component_low = this->fe_comp(fe_low_dim_, component_idx);
531  ASSERT_EQ(fe_component_low->fe_type(), FEType::FEScalar).error("Type of FiniteElement of scalar_shape accessor must be FEScalar!\n");
532  // use lambda reinit function
533  std::vector< std::vector<double> > shape_values_bulk( patch_point_vals_bulk_->get_quadrature()->size(), vector<double>(fe_component_low->n_dofs()) );
534  auto ref_shape_vals_bulk = this->ref_shape_values_bulk(patch_point_vals_bulk_->get_quadrature(), fe_component_low);
535  for (unsigned int i = 0; i < patch_point_vals_bulk_->get_quadrature()->size(); i++)
536  for (unsigned int j = 0; j < fe_component_low->n_dofs(); j++) {
537  shape_values_bulk[i][j] = ref_shape_vals_bulk[i][j][0];
538  }
539  uint scalar_shape_op_idx_bulk = patch_point_vals_bulk_->operations_.size(); // index in operations_ vector
540  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) {
541  bulk_reinit::ptop_scalar_shape(operations, op_results, shape_values_bulk, scalar_shape_op_idx_bulk);
542  };
543  auto &grad_scalar_shape_bulk_op = patch_point_vals_bulk_->make_fe_op({1}, lambda_scalar_shape_bulk, {}, fe_component_low->n_dofs());
544  uint begin_bulk = grad_scalar_shape_bulk_op.result_row();
545 
546  // element of higher dim (side points)
547  auto fe_component_high = this->fe_comp(fe_high_dim_, component_idx);
548  ASSERT_EQ(fe_component_high->fe_type(), FEType::FEScalar).error("Type of FiniteElement of scalar_shape accessor must be FEScalar!\n");
549  // use lambda reinit function
550  std::vector< std::vector< std::vector<double> > > shape_values_side(
551  dim+1,
553  );
554  auto ref_shape_vals_side = this->ref_shape_values_side(patch_point_vals_side_->get_quadrature(), fe_component_high);
555  for (unsigned int s=0; s<dim+1; ++s)
556  for (unsigned int i = 0; i < patch_point_vals_side_->get_quadrature()->size(); i++)
557  for (unsigned int j = 0; j < fe_component_high->n_dofs(); j++) {
558  shape_values_side[s][i][j] = ref_shape_vals_side[s][i][j][0];
559  }
560  uint scalar_shape_op_idx_side = patch_point_vals_side_->operations_.size(); // index in operations_ vector
561  auto lambda_scalar_shape_side = [shape_values_side, scalar_shape_op_idx_side](std::vector<ElOp<3>> &operations, TableDbl &op_results, TableInt &el_table) {
562  side_reinit::ptop_scalar_shape(operations, op_results, el_table, shape_values_side, scalar_shape_op_idx_side);
563  };
564  auto &grad_scalar_shape_side_op = patch_point_vals_side_->make_fe_op({1}, lambda_scalar_shape_side, {}, fe_component_high->n_dofs());
565  uint begin_side = grad_scalar_shape_side_op.result_row();
566 
567  auto bgn_it = make_iter<JoinShapeAccessor<Scalar>>( JoinShapeAccessor<Scalar>(patch_point_vals_bulk_, patch_point_vals_side_,
568  begin_bulk, begin_side, fe_component_low->n_dofs(), fe_component_high->n_dofs(), 0) );
569  unsigned int end_idx = fe_component_low->n_dofs() + fe_component_high->n_dofs();
570  auto end_it = make_iter<JoinShapeAccessor<Scalar>>( JoinShapeAccessor<Scalar>(patch_point_vals_bulk_, patch_point_vals_side_,
571  begin_bulk, begin_side, fe_component_low->n_dofs(), fe_component_high->n_dofs(), end_idx) );
572  return Range<JoinShapeAccessor<Scalar>>(bgn_it, end_it);
573  }
574 
575 private:
578  std::shared_ptr< FiniteElement<dim> > fe_high_dim_;
579  std::shared_ptr< FiniteElement<dim-1> > fe_low_dim_;
580 };
581 
582 /// Template specialization of dim = 1
583 template <>
584 class JoinValues<1> : public BaseValues<1>
585 {
586 public:
587  /// Constructor
589  : BaseValues<1>() {}
590 
592  {
594  make_iter<JoinShapeAccessor<Scalar>>(JoinShapeAccessor<Scalar>()),
596  }
597 };
598 
599 
600 template<unsigned int spacedim = 3>
602 public:
603  /// Struct for pre-computing number of elements, sides, bulk points and side points on each dimension.
604  struct TableSizes {
605  public:
606  /// Constructor
610  }
611 
612  /// Set all values to zero
613  void reset() {
614  std::fill(elem_sizes_[0].begin(), elem_sizes_[0].end(), 0);
615  std::fill(elem_sizes_[1].begin(), elem_sizes_[1].end(), 0);
616  std::fill(point_sizes_[0].begin(), point_sizes_[0].end(), 0);
617  std::fill(point_sizes_[1].begin(), point_sizes_[1].end(), 0);
618  }
619 
620  /**
621  * Holds number of elements and sides on each dimension
622  * Format:
623  * { {n_elements_1D, n_elements_2D, n_elements_3D },
624  * {n_sides_1D, n_sides_2D, n_sides_3D } }
625  */
627 
628  /**
629  * Holds number of bulk and side points on each dimension
630  * Format:
631  * { {n_bulk_points_1D, n_bulk_points_2D, n_bulk_points_3D },
632  * {n_side_points_1D, n_side_points_2D, n_side_points_3D } }
633  */
635  };
636 
638  : asm_arena_(1024 * 1024, 256),
639  patch_arena_(1024 * 1024, 256),
643  patch_point_vals_side_{ {FeSide::PatchPointValues(1, 0, patch_arena_),
644  FeSide::PatchPointValues(2, 0, patch_arena_),
645  FeSide::PatchPointValues(3, 0, patch_arena_)} }
646  {
647  used_quads_[0] = false; used_quads_[1] = false;
648  }
649 
650  PatchFEValues(unsigned int quad_order, MixedPtr<FiniteElement> fe)
651  : asm_arena_(1024 * 1024, 256),
652  patch_arena_(1024 * 1024, 256),
653  patch_point_vals_bulk_{ {FeBulk::PatchPointValues(1, quad_order, patch_arena_),
654  FeBulk::PatchPointValues(2, quad_order, patch_arena_),
655  FeBulk::PatchPointValues(3, quad_order, patch_arena_)} },
656  patch_point_vals_side_{ {FeSide::PatchPointValues(1, quad_order, patch_arena_),
657  FeSide::PatchPointValues(2, quad_order, patch_arena_),
658  FeSide::PatchPointValues(3, quad_order, patch_arena_)} },
659  fe_(fe)
660  {
661  used_quads_[0] = false; used_quads_[1] = false;
662  }
663 
664 
665  /// Destructor
667  {}
668 
669  /// Return bulk or side quadrature of given dimension
670  Quadrature *get_quadrature(uint dim, bool is_bulk) const {
671  if (is_bulk) return patch_point_vals_bulk_[dim-1].get_quadrature();
672  else return patch_point_vals_side_[dim-1].get_quadrature();
673  }
674 
675  /**
676  * @brief Initialize structures and calculates cell-independent data.
677  *
678  * @param _quadrature The quadrature rule for the cell associated
679  * to given finite element or for the cell side.
680  * @param _flags The update flags.
681  */
682  template<unsigned int DIM>
683  void initialize(Quadrature &_quadrature)
684  {
685  if ( _quadrature.dim() == DIM ) {
686  used_quads_[0] = true;
687  patch_point_vals_bulk_[DIM-1].initialize(3); // bulk
688  } else {
689  used_quads_[1] = true;
690  patch_point_vals_side_[DIM-1].initialize(4); // side
691  }
692  }
693 
694  /// Reset PatchpointValues structures
695  void reset()
696  {
697  for (unsigned int i=0; i<3; ++i) {
698  if (used_quads_[0]) patch_point_vals_bulk_[i].reset();
699  if (used_quads_[1]) patch_point_vals_side_[i].reset();
700  }
701  }
702 
703  /// Reinit data.
705  {
706  for (unsigned int i=0; i<3; ++i) {
707  if (used_quads_[0]) patch_point_vals_bulk_[i].reinit_patch();
708  if (used_quads_[1]) patch_point_vals_side_[i].reinit_patch();
709  }
710  }
711 
712  /**
713  * @brief Returns the number of shape functions.
714  */
715  template<unsigned int dim>
716  inline unsigned int n_dofs() const {
717  ASSERT((dim>=0) && (dim<=3))(dim).error("Dimension must be 0, 1, 2 or 3.");
718  return fe_[Dim<dim>{}]->n_dofs();
719  }
720 
721  /// Return BulkValue object of dimension given by template parameter
722  template<unsigned int dim>
724  ASSERT((dim>0) && (dim<=3))(dim).error("Dimension must be 1, 2 or 3.");
725  return BulkValues<dim>(patch_point_vals_bulk_[dim-1], fe_);
726  }
727 
728  /// Return SideValue object of dimension given by template parameter
729  template<unsigned int dim>
731  ASSERT((dim>0) && (dim<=3))(dim).error("Dimension must be 1, 2 or 3.");
732  return SideValues<dim>(patch_point_vals_side_[dim-1], fe_);
733  }
734 
735  /// Return JoinValue object of dimension given by template parameter
736  template<unsigned int dim>
738  //ASSERT((dim>1) && (dim<=3))(dim).error("Dimension must be 2 or 3.");
739  return JoinValues<dim>(&patch_point_vals_bulk_[dim-2], &patch_point_vals_side_[dim-1], fe_);
740  }
741 
742  /** Following methods are used during update of patch. **/
743 
744  /// Resize tables of patch_point_vals_
745  void resize_tables(TableSizes table_sizes) {
746  for (uint i=0; i<3; ++i) {
747  if (used_quads_[0]) patch_point_vals_bulk_[i].resize_tables(table_sizes.elem_sizes_[0][i], table_sizes.point_sizes_[0][i]);
748  if (used_quads_[1]) patch_point_vals_side_[i].resize_tables(table_sizes.elem_sizes_[1][i], table_sizes.point_sizes_[1][i]);
749  }
750  }
751 
752  /// Register element to patch_point_vals_ table by dimension of element
753  uint register_element(DHCellAccessor cell, uint element_patch_idx) {
754  arma::mat coords;
755  switch (cell.dim()) {
756  case 1:
757  coords = MappingP1<1,spacedim>::element_map(cell.elm());
758  return patch_point_vals_bulk_[0].register_element(coords, element_patch_idx);
759  break;
760  case 2:
761  coords = MappingP1<2,spacedim>::element_map(cell.elm());
762  return patch_point_vals_bulk_[1].register_element(coords, element_patch_idx);
763  break;
764  case 3:
765  coords = MappingP1<3,spacedim>::element_map(cell.elm());
766  return patch_point_vals_bulk_[2].register_element(coords, element_patch_idx);
767  break;
768  default:
769  ASSERT(false);
770  return 0;
771  break;
772  }
773  }
774 
775  /// Register side to patch_point_vals_ table by dimension of side
777  arma::mat side_coords(spacedim, cell_side.dim());
778  for (unsigned int n=0; n<cell_side.dim(); n++)
779  for (unsigned int c=0; c<spacedim; c++)
780  side_coords(c,n) = (*cell_side.side().node(n))[c];
781 
782  arma::mat elm_coords;
783  DHCellAccessor cell = cell_side.cell();
784  switch (cell.dim()) {
785  case 1:
786  elm_coords = MappingP1<1,spacedim>::element_map(cell.elm());
787  return patch_point_vals_side_[0].register_side(elm_coords, side_coords);
788  break;
789  case 2:
790  elm_coords = MappingP1<2,spacedim>::element_map(cell.elm());
791  return patch_point_vals_side_[1].register_side(elm_coords, side_coords);
792  break;
793  case 3:
794  elm_coords = MappingP1<3,spacedim>::element_map(cell.elm());
795  return patch_point_vals_side_[2].register_side(elm_coords, side_coords);
796  break;
797  default:
798  ASSERT(false);
799  return 0;
800  break;
801  }
802  }
803 
804  /// Register bulk point to patch_point_vals_ table by dimension of element
805  uint register_bulk_point(DHCellAccessor cell, uint elem_table_row, uint value_patch_idx) {
806  return patch_point_vals_bulk_[cell.dim()-1].register_bulk_point(elem_table_row, value_patch_idx, cell.elm_idx());
807  }
808 
809  /// Register side point to patch_point_vals_ table by dimension of side
810  uint register_side_point(DHCellSide cell_side, uint elem_table_row, uint value_patch_idx) {
811  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());
812  }
813 
814  /// Temporary development method
815  void print_data_tables(ostream& stream, bool points, bool ints, bool only_bulk=true) const {
816  stream << endl << "Table of patch FE data:" << endl;
817  for (uint i=0; i<3; ++i) {
818  stream << std::setfill('-') << setw(100) << "" << endl;
819  stream << "Bulk, dimension " << (i+1) << endl;
820  patch_point_vals_bulk_[i].print_data_tables(stream, points, ints);
821  }
822  if (!only_bulk)
823  for (uint i=0; i<3; ++i) {
824  stream << std::setfill('-') << setw(100) << "" << endl;
825  stream << "Side, dimension " << (i+1) << endl;
826  patch_point_vals_side_[i].print_data_tables(stream, points, ints);
827  }
828  stream << std::setfill('=') << setw(100) << "" << endl;
829  }
830 
831  /// Temporary development method
832  void print_operations(ostream& stream) const {
833  stream << endl << "Table of patch FE operations:" << endl;
834  for (uint i=0; i<3; ++i) {
835  stream << std::setfill('-') << setw(100) << "" << endl;
836  stream << "Bulk, dimension " << (i+1) << ", n_rows " << patch_point_vals_bulk_[i].n_rows() << endl;
837  patch_point_vals_bulk_[i].print_operations(stream, 0);
838  }
839  for (uint i=0; i<3; ++i) {
840  stream << std::setfill('-') << setw(100) << "" << endl;
841  stream << "Side, dimension " << (i+1) << ", n_rows " << patch_point_vals_side_[i].n_rows() << endl;
842  patch_point_vals_side_[i].print_operations(stream, 1);
843  }
844  stream << std::setfill('=') << setw(100) << "" << endl;
845  }
846 
847 private:
850  std::array<FeBulk::PatchPointValues<spacedim>, 3> patch_point_vals_bulk_; ///< Sub objects of bulk data of dimensions 1,2,3
851  std::array<FeSide::PatchPointValues<spacedim>, 3> patch_point_vals_side_; ///< Sub objects of side data of dimensions 1,2,3
852 
853  MixedPtr<FiniteElement> fe_; ///< Mixed of shared pointers of FiniteElement object
854  bool used_quads_[2]; ///< Pair of flags signs holds info if bulk and side quadratures are used
855 
856  template <class ValueType>
857  friend class ElQ;
858  template <class ValueType>
859  friend class FeQ;
860  template <class ValueType>
861  friend class JoinShapeAccessor;
862 };
863 
864 
865 template <class ValueType>
866 ValueType ElQ<ValueType>::operator()(const BulkPoint &point) {
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_.scalar_value(op_idx_, 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_.vector_value(op_idx_, value_cache_idx);
875 }
876 
877 template <>
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_.tensor_value(op_idx_, value_cache_idx);
881 }
882 
883 template <class ValueType>
884 ValueType ElQ<ValueType>::operator()(const SidePoint &point) {
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_.scalar_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_.vector_val(begin_, value_cache_idx);
893 }
894 
895 template <>
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_.tensor_val(begin_, value_cache_idx);
899 }
900 
901 template <class ValueType>
902 ValueType FeQ<ValueType>::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_.scalar_val(begin_+shape_idx, value_cache_idx);
905 }
906 
907 template <>
908 inline Vector FeQ<Vector>::operator()(unsigned int shape_idx, const BulkPoint &point) {
909  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
910  return patch_point_vals_.vector_val(begin_+3*shape_idx, value_cache_idx);
911 }
912 
913 template <>
914 inline Tensor FeQ<Tensor>::operator()(FMT_UNUSED unsigned int shape_idx, FMT_UNUSED const BulkPoint &point) {
915  Tensor tens; tens.zeros();
916  return tens;
917 }
918 
919 template <class ValueType>
920 ValueType FeQ<ValueType>::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_.scalar_val(begin_+shape_idx, value_cache_idx);
923 }
924 
925 template <>
926 inline Vector FeQ<Vector>::operator()(unsigned int shape_idx, const SidePoint &point) {
927  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
928  return patch_point_vals_.vector_val(begin_+3*shape_idx, value_cache_idx);
929 }
930 
931 template <>
932 inline Tensor FeQ<Tensor>::operator()(FMT_UNUSED unsigned int shape_idx, FMT_UNUSED const SidePoint &point) {
933  Tensor tens; tens.zeros();
934  return tens;
935 }
936 
937 
938 template <class ValueType>
940  if (this->is_high_dim()) {
941  return 0.0;
942  } else {
943  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
944  return patch_point_vals_bulk_->scalar_val(begin_+this->local_idx(), value_cache_idx);
945  }
946 }
947 
948 template <>
950  Vector vect; vect.zeros();
951  return vect;
952 }
953 
954 template <>
956  Tensor tens; tens.zeros();
957  return tens;
958 }
959 
960 template <class ValueType>
962  if (this->is_high_dim()) {
963  unsigned int value_cache_idx = point.elm_cache_map()->element_eval_point(point.elem_patch_idx(), point.eval_point_idx());
964  return patch_point_vals_side_->scalar_val(begin_side_+this->local_idx(), value_cache_idx);
965  } else {
966  return 0.0;
967  }
968 }
969 
970 template <>
972  Vector vect; vect.zeros();
973  return vect;
974 }
975 
976 template <>
978  Tensor tens; tens.zeros();
979  return tens;
980 }
981 
982 
983 #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(PatchPointValues< 3 > &patch_point_vals, unsigned int begin, unsigned int op_idx)
Constructor.
unsigned int op_idx_
Index of the first component of the bulk Quantity. Size is given by ValueType.
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)
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:35
Eigen::Vector< ArrayInt, Eigen::Dynamic > TableInt
Definition: eigen_tools.hh:36
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.