Flow123d  JS_before_hm-883-gc471082
eval_points.cc
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 eval_points.cc
15  * @brief
16  * @author David Flanderka
17  */
18 
19 #include "fields/eval_points.hh"
20 #include "fields/eval_subset.hh"
21 #include "quadrature/quadrature.hh"
22 #include "mesh/ref_element.hh"
23 #include <memory>
24 
25 
26 const unsigned int EvalPoints::undefined_dim = 10;
27 
29 : dim_eval_points_({DimEvalPoints(1), DimEvalPoints(2), DimEvalPoints(3)})
30 {}
31 
32 template <unsigned int dim>
33 std::shared_ptr<BulkIntegral> EvalPoints::add_bulk(const Quadrature &quad)
34 {
35  ASSERT_EQ(dim, quad.dim());
36  std::shared_ptr<BulkIntegral> bulk_integral = std::make_shared<BulkIntegral>(shared_from_this(), dim);
37  dim_eval_points_[dim-1].add_local_points<dim>( quad.get_points() );
38  dim_eval_points_[dim-1].add_subset();
39  return bulk_integral;
40 }
41 
42 template <unsigned int dim>
43 std::shared_ptr<EdgeIntegral> EvalPoints::add_edge(const Quadrature &quad)
44 {
45  ASSERT_EQ(dim, quad.dim()+1);
46  unsigned int old_data_size=this->size(dim), new_data_size; // interval of side subset data
47  unsigned int points_per_side = quad.make_from_side<dim>(0, 0).get_points().size();
48  unsigned int n_side_permutations = RefElement<dim>::n_side_permutations;
49 
50  std::shared_ptr<EdgeIntegral> edge_integral = std::make_shared<EdgeIntegral>(shared_from_this(), dim, n_side_permutations, points_per_side);
51  unsigned int*** perm_indices = edge_integral->perm_indices_;
52 
53  // permutation 0
54  for (unsigned int i=0; i<dim+1; ++i) { // sides
55  Quadrature high_dim_q = quad.make_from_side<dim>(i, 0);
56  dim_eval_points_[dim-1].add_local_points<dim>( high_dim_q.get_points() );
57  }
58  dim_eval_points_[dim-1].add_subset();
59  new_data_size = this->size(dim);
60  unsigned int i_data=old_data_size;
61  for (unsigned int i_side=0; i_side<dim+1; ++i_side) {
62  for (unsigned int i_point=0; i_point<points_per_side; ++i_point) {
63  perm_indices[i_side][0][i_point] = i_data;
64  ++i_data;
65  }
66  }
67 
68  // permutation 1...N
69  for (unsigned int i_perm=1; i_perm<n_side_permutations; ++i_perm) {
70  for (unsigned int i_side=0; i_side<dim+1; ++i_side) {
71  Quadrature high_dim_q = quad.make_from_side<dim>(i_side, i_perm);
72  const Armor::Array<double> & quad_points = high_dim_q.get_points();
73  for (unsigned int i_point=0; i_point<quad_points.size(); ++i_point) {
74  perm_indices[i_side][i_perm][i_point] = dim_eval_points_[dim-1].find_permute_point<dim>( quad_points.vec<dim>(i_point), old_data_size, new_data_size );
75  }
76  }
77  }
78 
79  return edge_integral;
80 }
81 
82 template <unsigned int dim>
83 std::shared_ptr<CouplingIntegral> EvalPoints::add_coupling(const Quadrature &quad) {
84  ASSERT_EQ(dim, quad.dim()+1);
85 
86  std::shared_ptr<BulkIntegral> bulk_integral = this->add_bulk<dim-1>(quad);
87  std::shared_ptr<EdgeIntegral> edge_integral = this->add_edge<dim>(quad);
88  return std::make_shared<CouplingIntegral>(edge_integral, bulk_integral);
89 }
90 
91 template <unsigned int dim>
92 std::shared_ptr<BoundaryIntegral> EvalPoints::add_boundary(const Quadrature &quad) {
93  ASSERT_EQ(dim, quad.dim()+1);
94  return std::make_shared<BoundaryIntegral>(this->add_edge<dim>(quad));
95 }
96 
98 : local_points_(dim), n_subsets_(0), dim_(dim)
99 {
100  subset_starts_[0] = 0;
102 }
103 
104 template <unsigned int dim>
106  unsigned int local_points_old_size = local_points_.size();
107  local_points_.resize(quad_points.size() + local_points_old_size);
108  for (unsigned int i=0; i<quad_points.size(); ++i) {
109  local_points_.set(i+local_points_old_size) = quad_points.vec<dim>(i);
110  }
111 }
112 
113 template <unsigned int dim>
114 unsigned int EvalPoints::DimEvalPoints::find_permute_point(arma::vec coords, unsigned int data_begin, unsigned int data_end) {
115  for (unsigned int loc_idx=data_begin; loc_idx<data_end; ++loc_idx) {
116  // Check if point exists in local points vector.
117  if ( arma::norm(coords-local_points_.vec<dim>(loc_idx), 2) < 4*std::numeric_limits<double>::epsilon() ) return loc_idx;
118  }
119 
120  ASSERT(false);
121  return 0;
122 }
123 
125  ASSERT_LT_DBG(n_subsets_, EvalPoints::max_subsets).error("Maximal number of subsets exceeded!\n");
126 
127  n_subsets_++;
128  subset_starts_[n_subsets_] = this->size();
129 }
130 
131 
132 template std::shared_ptr<BulkIntegral> EvalPoints::add_bulk<1>(const Quadrature &);
133 template std::shared_ptr<BulkIntegral> EvalPoints::add_bulk<2>(const Quadrature &);
134 template std::shared_ptr<BulkIntegral> EvalPoints::add_bulk<3>(const Quadrature &);
135 template std::shared_ptr<EdgeIntegral> EvalPoints::add_edge<1>(const Quadrature &);
136 template std::shared_ptr<EdgeIntegral> EvalPoints::add_edge<2>(const Quadrature &);
137 template std::shared_ptr<EdgeIntegral> EvalPoints::add_edge<3>(const Quadrature &);
138 template std::shared_ptr<CouplingIntegral> EvalPoints::add_coupling<2>(const Quadrature &);
139 template std::shared_ptr<CouplingIntegral> EvalPoints::add_coupling<3>(const Quadrature &);
140 template std::shared_ptr<BoundaryIntegral> EvalPoints::add_boundary<1>(const Quadrature &);
141 template std::shared_ptr<BoundaryIntegral> EvalPoints::add_boundary<2>(const Quadrature &);
142 template std::shared_ptr<BoundaryIntegral> EvalPoints::add_boundary<3>(const Quadrature &);
143 template void EvalPoints::DimEvalPoints::add_local_points<1>(const Armor::Array<double> &);
144 template void EvalPoints::DimEvalPoints::add_local_points<2>(const Armor::Array<double> &);
145 template void EvalPoints::DimEvalPoints::add_local_points<3>(const Armor::Array<double> &);
146 template unsigned int EvalPoints::DimEvalPoints::find_permute_point<1>(arma::vec, unsigned int, unsigned int);
147 template unsigned int EvalPoints::DimEvalPoints::find_permute_point<2>(arma::vec, unsigned int, unsigned int);
148 template unsigned int EvalPoints::DimEvalPoints::find_permute_point<3>(arma::vec, unsigned int, unsigned int);
static constexpr unsigned int max_subsets
Maximal number of hold subsets.
Definition: eval_points.hh:49
std::shared_ptr< BulkIntegral > add_bulk(const Quadrature &)
Definition: eval_points.cc:33
void resize(uint size)
Definition: armor.hh:700
Armor::Array< double > local_points_
Local coords of points vector.
Definition: eval_points.hh:164
unsigned int size() const
Definition: armor.hh:718
ArmaVec< double, N > vec
Definition: armor.hh:861
unsigned int dim() const
Definition: quadrature.hh:73
unsigned int n_subsets_
Number of subset.
Definition: eval_points.hh:166
const Armor::Array< double > & get_points() const
Return a reference to the whole array of quadrature points.
Definition: quadrature.hh:104
Quadrature make_from_side(unsigned int sid, unsigned int pid) const
Definition: quadrature.cc:47
static const unsigned int max_subset_points
Maximal average number of points hold in subset.
Definition: eval_points.hh:52
DimEvalPoints(unsigned int dim)
Constructor.
Definition: eval_points.cc:97
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library...
Definition: asserts.hh:347
ArmaVec< Type, nr > vec(uint mat_index) const
Definition: armor.hh:807
static const unsigned int undefined_dim
Undefined dimension of new (empty) object.
Definition: eval_points.hh:46
Base class for quadrature rules on simplices in arbitrary dimensions.
Definition: quadrature.hh:48
std::shared_ptr< EdgeIntegral > add_edge(const Quadrature &)
The same as add_bulk but for edge points on sides.
Definition: eval_points.cc:43
void reinit(uint size)
Definition: armor.hh:688
std::array< DimEvalPoints, 3 > dim_eval_points_
Sub objects of dimensions 1,2,3.
Definition: eval_points.hh:171
std::shared_ptr< CouplingIntegral > add_coupling(const Quadrature &)
The same as add_bulk but for points between side points of element of dim and bulk points of element ...
Definition: eval_points.cc:83
unsigned int find_permute_point(arma::vec coords, unsigned int data_begin, unsigned int data_end)
Find position of local point (coords) in subvector of local points given by limits <data_begin...
Definition: eval_points.cc:114
Basic definitions of numerical quadrature rules.
ArrayMatSet set(uint index)
Definition: armor.hh:821
EvalPoints()
Constructor.
Definition: eval_points.cc:28
unsigned int size() const
Return size of evaluation points object (number of points).
Definition: eval_points.hh:119
unsigned int size(unsigned int dim) const
Return size of evaluation points object (number of points).
Definition: eval_points.hh:58
void add_local_points(const Armor::Array< double > &quad_points)
Adds set of local point to local_points_ (bulk or side of given permutation).
Definition: eval_points.cc:105
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc...
std::shared_ptr< BoundaryIntegral > add_boundary(const Quadrature &)
The same as add_bulk but for edge points on boundary sides.
Definition: eval_points.cc:92
std::array< int, EvalPoints::max_subsets+1 > subset_starts_
Indices of subsets data in local_points_ vector, used size is n_subsets_ + 1.
Definition: eval_points.hh:165
unsigned int size() const
Returns number of quadrature points.
Definition: quadrature.hh:87
void add_subset()
Adds new subset and its end size to subset_starts_ array.
Definition: eval_points.cc:124
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual)
Definition: asserts.hh:328
#define ASSERT_LT_DBG(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:300