Flow123d  master-f44eb46
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(0), DimEvalPoints(1), DimEvalPoints(2), DimEvalPoints(3)}), max_size_(0)
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 
37  dim_eval_points_[dim].add_local_points<dim>( quad.get_points() );
38  uint i_subset = dim_eval_points_[dim].add_subset();
39  std::shared_ptr<BulkIntegral> bulk_integral = std::make_shared<BulkIntegral>(shared_from_this(), dim, i_subset);
40  this->set_max_size();
41  return bulk_integral;
42 }
43 
44 template <>
45 std::shared_ptr<BulkIntegral> EvalPoints::add_bulk<0>(const Quadrature &quad)
46 {
47  ASSERT_EQ(0, quad.dim());
48  uint i_subset = dim_eval_points_[0].add_subset();
49  std::shared_ptr<BulkIntegral> bulk_integral = std::make_shared<BulkIntegral>(shared_from_this(), 0, i_subset);
50  this->set_max_size();
51  return bulk_integral;
52 }
53 
54 template <unsigned int dim>
55 std::shared_ptr<EdgeIntegral> EvalPoints::add_edge(const Quadrature &quad)
56 {
57  ASSERT_EQ(dim, quad.dim()+1);
58 
59  for (unsigned int i=0; i<dim+1; ++i) { // sides
60  Quadrature high_dim_q = quad.make_from_side<dim>(i);
61  dim_eval_points_[dim].add_local_points<dim>( high_dim_q.get_points() );
62  }
63  uint i_subset = dim_eval_points_[dim].add_subset();
64  std::shared_ptr<EdgeIntegral> edge_integral = std::make_shared<EdgeIntegral>(shared_from_this(), dim, i_subset);
65 
66  this->set_max_size();
67  return edge_integral;
68 }
69 
70 template <unsigned int dim>
71 std::shared_ptr<CouplingIntegral> EvalPoints::add_coupling(const Quadrature &quad) {
72  ASSERT_EQ(dim, quad.dim()+1);
73 
74  std::shared_ptr<BulkIntegral> bulk_integral = this->add_bulk<dim-1>(quad);
75  DebugOut() << "coupling bulk subset" << bulk_integral->get_subset_idx();
76  std::shared_ptr<EdgeIntegral> edge_integral = this->add_edge<dim>(quad);
77  DebugOut() << "coupling edge subset" << edge_integral->get_subset_idx();
78  return std::make_shared<CouplingIntegral>(edge_integral, bulk_integral);
79 }
80 
81 template <unsigned int dim>
82 std::shared_ptr<BoundaryIntegral> EvalPoints::add_boundary(const Quadrature &quad) {
83  ASSERT_EQ(dim, quad.dim()+1);
84 
85  std::shared_ptr<BulkIntegral> bulk_integral = this->add_bulk<dim-1>(quad);
86  DebugOut() << "boundary bulk subset: " << bulk_integral->get_subset_idx()
87  << "begin: " << subset_begin(dim-1, bulk_integral->get_subset_idx());
88  std::shared_ptr<EdgeIntegral> edge_integral = this->add_edge<dim>(quad);
89  DebugOut() << "boundary edge subset" << edge_integral->get_subset_idx();
90  return std::make_shared<BoundaryIntegral>(edge_integral, bulk_integral);
91 }
92 
94 : local_points_(dim), n_subsets_(0), dim_(dim)
95 {
96  subset_starts_[0] = 0;
98 }
99 
100 template <unsigned int dim>
102  ASSERT_GT(dim, 0).error("Dimension 0 not supported!\n");
103  unsigned int local_points_old_size = local_points_.size();
104  local_points_.resize(quad_points.size() + local_points_old_size);
105  for (unsigned int i=0; i<quad_points.size(); ++i) {
106  //DebugOut() << "add i: " << i << " : " << quad_points.vec<dim>(i);
107 
108  local_points_.set(i+local_points_old_size) = quad_points.vec<dim>(i);
109  }
110 }
111 
112 
114  ASSERT_LT(n_subsets_, EvalPoints::max_subsets).error("Maximal number of subsets exceeded!\n");
115 
116  n_subsets_++;
117  subset_starts_[n_subsets_] = this->size();
118  return n_subsets_ - 1;
119 }
120 
121 
122 template std::shared_ptr<BulkIntegral> EvalPoints::add_bulk<0>(const Quadrature &);
123 template std::shared_ptr<BulkIntegral> EvalPoints::add_bulk<1>(const Quadrature &);
124 template std::shared_ptr<BulkIntegral> EvalPoints::add_bulk<2>(const Quadrature &);
125 template std::shared_ptr<BulkIntegral> EvalPoints::add_bulk<3>(const Quadrature &);
126 template std::shared_ptr<EdgeIntegral> EvalPoints::add_edge<1>(const Quadrature &);
127 template std::shared_ptr<EdgeIntegral> EvalPoints::add_edge<2>(const Quadrature &);
128 template std::shared_ptr<EdgeIntegral> EvalPoints::add_edge<3>(const Quadrature &);
129 template std::shared_ptr<CouplingIntegral> EvalPoints::add_coupling<2>(const Quadrature &);
130 template std::shared_ptr<CouplingIntegral> EvalPoints::add_coupling<3>(const Quadrature &);
131 template std::shared_ptr<BoundaryIntegral> EvalPoints::add_boundary<1>(const Quadrature &);
132 template std::shared_ptr<BoundaryIntegral> EvalPoints::add_boundary<2>(const Quadrature &);
133 template std::shared_ptr<BoundaryIntegral> EvalPoints::add_boundary<3>(const Quadrature &);
134 template void EvalPoints::DimEvalPoints::add_local_points<1>(const Armor::Array<double> &);
135 template void EvalPoints::DimEvalPoints::add_local_points<2>(const Armor::Array<double> &);
136 template void EvalPoints::DimEvalPoints::add_local_points<3>(const Armor::Array<double> &);
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:301
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
#define ASSERT_GT(a, b)
Definition of comparative assert macro (Greater Than) only for debug mode.
Definition: asserts.hh:317
void reinit(uint size)
Definition: armor.hh:698
ArmaVec< Type, nr > vec(uint mat_index) const
Definition: armor.hh:821
unsigned int size() const
Definition: armor.hh:728
Subobject holds evaluation points data of one dimension (0,1,2,3)
Definition: eval_points.hh:120
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:174
DimEvalPoints(unsigned int dim)
Constructor.
Definition: eval_points.cc:93
uint add_subset()
Adds new subset and its end size to subset_starts_ array.
Definition: eval_points.cc:113
Armor::Array< double > local_points_
Local coords of points vector.
Definition: eval_points.hh:173
void add_local_points(const Armor::Array< double > &quad_points)
Adds set of local point to local_points_ (bulk or side).
Definition: eval_points.cc:101
int subset_begin(unsigned int dim, unsigned int idx) const
Return begin index of appropriate subset data.
Definition: eval_points.hh:70
void set_max_size()
Definition: eval_points.hh:179
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:71
static const unsigned int undefined_dim
Undefined dimension of new (empty) object.
Definition: eval_points.hh:46
unsigned int size(unsigned int dim) const
Return size of evaluation points object (number of points).
Definition: eval_points.hh:58
std::shared_ptr< BoundaryIntegral > add_boundary(const Quadrature &)
The same as add_bulk but for edge points on boundary sides.
Definition: eval_points.cc:82
static const unsigned int max_subset_points
Maximal average number of points hold in subset.
Definition: eval_points.hh:52
std::shared_ptr< BulkIntegral > add_bulk(const Quadrature &)
Definition: eval_points.cc:33
std::shared_ptr< EdgeIntegral > add_edge(const Quadrature &)
The same as add_bulk but for edge points on sides.
Definition: eval_points.cc:55
std::array< DimEvalPoints, 4 > dim_eval_points_
Sub objects of dimensions 0,1,2,3.
Definition: eval_points.hh:184
static constexpr unsigned int max_subsets
Maximal number of hold subsets.
Definition: eval_points.hh:49
EvalPoints()
Constructor.
Definition: eval_points.cc:28
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
const Armor::Array< double > & get_points() const
Return a reference to the whole array of quadrature points.
Definition: quadrature.hh:99
unsigned int dim() const
Definition: quadrature.hh:72
#define DebugOut()
Macro defining 'debug' record of log.
Definition: logger.hh:284
unsigned int uint
Basic definitions of numerical quadrature rules.
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.