Flow123d  master-f44eb46
assembly_observe.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 assembly_observe.hh
15  * @brief
16  */
17 
18 #ifndef ASSEMBLY_OBSERVE_HH_
19 #define ASSEMBLY_OBSERVE_HH_
20 
21 #include <unordered_map>
22 
25 #include "fem/dofhandler.hh"
27 #include "io/observe.hh"
28 
29 
30 /**
31  * @brief Generic class of observe output assemblation.
32  *
33  * Class
34  * - holds assemblation structures (EvalPoints, Integral objects, Integral data table).
35  * - associates assemblation objects specified by dimension
36  * - provides general assemble method
37  * - provides methods that allow construction of element patches
38  */
39 template < template<IntDim...> class DimAssembly>
41 {
42 public:
43  /// Constructor
44  GenericAssemblyObserve( typename DimAssembly<1>::EqFields *eq_fields, const std::unordered_set<string> &observe_fields_list,
45  std::shared_ptr<Observe> observe)
46  : multidim_assembly_(eq_fields, observe_fields_list, observe.get()), observe_(observe), bulk_integral_data_(20, 10)
47  {
48  eval_points_ = std::make_shared<EvalPoints>();
49  multidim_assembly_[1_d]->create_observe_integrals(eval_points_, integrals_);
50  multidim_assembly_[2_d]->create_observe_integrals(eval_points_, integrals_);
51  multidim_assembly_[3_d]->create_observe_integrals(eval_points_, integrals_);
53  multidim_assembly_[1_d]->initialize(&element_cache_map_);
54  multidim_assembly_[2_d]->initialize(&element_cache_map_);
55  multidim_assembly_[3_d]->initialize(&element_cache_map_);
56  }
57 
58  /// Getter to set of assembly objects
60  return multidim_assembly_;
61  }
62 
63  /**
64  * @brief General assemble methods.
65  *
66  * Loops through local cells and calls assemble methods of assembly
67  * object of each cells over space dimension.
68  */
69  void assemble(std::shared_ptr<DOFHandlerMultiDim> dh) override {
70  START_TIMER( DimAssembly<1>::name() );
71 
72  unsigned int i_ep, subset_begin, subset_idx;
73  auto &patch_point_data = observe_->patch_point_data();
74  for(auto & p_data : patch_point_data) {
75  subset_idx = integrals_.bulk_[p_data.i_quad]->get_subset_idx();
76  subset_begin = eval_points_->subset_begin(p_data.i_quad+1, subset_idx);
77  i_ep = subset_begin + p_data.i_quad_point;
78  DHCellAccessor dh_cell = dh->cell_accessor_from_element(p_data.elem_idx);
79  bulk_integral_data_.emplace_back(dh_cell, p_data.i_quad_point);
80  element_cache_map_.eval_point_data_.emplace_back(p_data.i_reg, p_data.elem_idx, i_ep, 0);
81  }
84 
85  this->reallocate_cache();
87  multidim_assembly_[1_d]->eq_fields_->cache_update(element_cache_map_);
88 
89  multidim_assembly_[1_d]->assemble_cell_integrals(bulk_integral_data_);
90  multidim_assembly_[2_d]->assemble_cell_integrals(bulk_integral_data_);
91  multidim_assembly_[3_d]->assemble_cell_integrals(bulk_integral_data_);
94  END_TIMER( DimAssembly<1>::name() );
95  }
96 
97 
98 private:
99  /// Calls cache_reallocate method on set of used fields
100  inline void reallocate_cache() {
101  multidim_assembly_[1_d]->eq_fields_->cache_reallocate(this->element_cache_map_, multidim_assembly_[1_d]->used_fields_);
102  // DebugOut() << "Order of evaluated fields (" << DimAssembly<1>::name() << "):" << multidim_assembly_[1_d]->eq_fields_->print_dependency();
103  }
104 
106  std::shared_ptr<Observe> observe_; ///< Shared Observe object.
107  RevertableList<BulkIntegralData> bulk_integral_data_; ///< Holds data for computing bulk integrals.
108 };
109 
110 
111 template <unsigned int dim>
113 {
114 public:
116 
117  static constexpr const char * name() { return "AssemblyObserveOutput"; }
118 
119  /// Constructor.
120  AssemblyObserveOutput(EqFields *eq_fields, const std::unordered_set<string> &observe_fields_list, Observe *observe)
121  : AssemblyBase<dim>(), eq_fields_(eq_fields), observe_(observe) {
123  offsets_.resize(1.1 * CacheMapElementNumber::get());
124 
125  for (auto observe_field : observe_fields_list) {
126  auto found_field = eq_fields_->field(observe_field);
127  used_fields_ += *found_field;
128  }
129  }
130 
131  /// Destructor.
133 
134  /// Initialize auxiliary vectors and other data members
135  void initialize(ElementCacheMap *element_cache_map) {
136  this->element_cache_map_ = element_cache_map;
137  }
138 
139  /// Assembles the cell integrals for the given dimension.
141  unsigned int element_patch_idx, field_value_cache_position, val_idx;
142  this->reset_offsets();
143  for (unsigned int i=0; i<bulk_integral_data.permanent_size(); ++i) {
144  if (bulk_integral_data[i].cell.dim() != dim) continue;
145  element_patch_idx = this->element_cache_map_->position_in_cache(bulk_integral_data[i].cell.elm_idx());
146  auto p = *( this->bulk_points(element_patch_idx).begin()); // evaluation point
147  field_value_cache_position = this->element_cache_map_->element_eval_point(element_patch_idx, p.eval_point_idx() + bulk_integral_data[i].subset_index);
149  this->offsets_[field_value_cache_position] = val_idx;
150  }
151  for (FieldListAccessor f_acc : this->used_fields_.fields_range()) {
152  f_acc->fill_observe_value(observe_->get_output_cache(f_acc->name()), this->offsets_);
153  }
154  }
155 
156 
157  /// Create bulk integral according to dim
158  void create_observe_integrals(std::shared_ptr<EvalPoints> eval_points, AssemblyIntegrals &integrals) {
159  std::vector<arma::vec> reg_points;
160 
161  auto &patch_point_data = observe_->patch_point_data();
162  for(auto & p_data : patch_point_data) {
163  auto el_acc = eq_fields_->mesh()->element_accessor(p_data.elem_idx);
164  if (el_acc.dim()!=dim) continue;
165  p_data.i_reg = el_acc.region_idx().idx();
166  p_data.i_quad = el_acc.dim() - 1;
167  p_data.i_quad_point = reg_points.size();
168  reg_points.push_back(p_data.local_coords);
169  }
170 
171  if (reg_points.size() > 0) {
172  this->quad_ = new Quadrature(dim, reg_points.size());
173  for (uint j=0; j<reg_points.size(); j++) {
174  arma::vec::fixed<dim> fix_p = reg_points[j].subvec(0, dim-1);
175  this->quad_->weight(j) = 1.0;
176  this->quad_->set(j) = fix_p;
177  }
178  this->integrals_.bulk_ = eval_points->add_bulk<dim>(*this->quad_);
179  integrals.bulk_[dim-1] = this->integrals_.bulk_;
180  }
181  }
182 
183 private:
184  void reset_offsets() {
185  std::fill(offsets_.begin(), offsets_.end(), -1);
186  }
187 
188  /// Data objects shared with EquationOutput
191 
192  FieldSet used_fields_; ///< Sub field set contains fields performed to output
193  std::vector<int> offsets_; ///< Holds indices (offsets) of cached data to output data vector
194 
195 
196  template < template<IntDim...> class DimAssembly>
198 };
199 
200 
201 
202 
203 #endif /* ASSEMBLY_OBSERVE_HH_ */
Range< BulkPoint > bulk_points(unsigned int element_patch_idx) const
Return BulkPoint range of appropriate dimension.
Quadrature * quad_
Quadrature used in assembling methods.
DimIntegrals integrals_
Set of used integrals.
int active_integrals_
Holds mask of active integrals.
ElementCacheMap * element_cache_map_
ElementCacheMap shared with GenericAssembly object.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
void assemble_cell_integrals(const RevertableList< GenericAssemblyBase::BulkIntegralData > &bulk_integral_data)
Assembles the cell integrals for the given dimension.
EqFields * eq_fields_
Data objects shared with EquationOutput.
AssemblyObserveOutput(EqFields *eq_fields, const std::unordered_set< string > &observe_fields_list, Observe *observe)
Constructor.
void create_observe_integrals(std::shared_ptr< EvalPoints > eval_points, AssemblyIntegrals &integrals)
Create bulk integral according to dim.
FieldSet used_fields_
Sub field set contains fields performed to output.
~AssemblyObserveOutput()
Destructor.
std::vector< int > offsets_
Holds indices (offsets) of cached data to output data vector.
static constexpr const char * name()
static unsigned int get()
Return number of stored elements.
Cell accessor allow iterate over DOF handler cells.
RegionIdx region_idx() const
Definition: accessors.hh:201
Directing class of FieldValueCache.
void clear_element_eval_points_map()
Reset all items of elements_eval_points_map.
void create_patch()
Create patch of cached elements before reading data to cache.
RevertableList< EvalPointData > eval_point_data_
void init(std::shared_ptr< EvalPoints > eval_points)
Init cache.
int element_eval_point(unsigned int i_elem_in_cache, unsigned int i_eval_point) const
unsigned int position_in_cache(unsigned mesh_elm_idx, bool bdr=false) const
Return position of element stored in ElementCacheMap.
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
Range< FieldListAccessor > fields_range() const
Returns range of Fields held in field_list.
Definition: field_set.cc:275
const Mesh * mesh() const
Returns pointer to mesh.
Definition: field_set.hh:393
FieldCommon * field(const std::string &field_name) const
Definition: field_set.cc:173
ElementCacheMap element_cache_map_
ElementCacheMap according to EvalPoints.
std::shared_ptr< EvalPoints > eval_points_
EvalPoints object shared by all integrals.
AssemblyIntegrals integrals_
Holds integral objects.
Generic class of observe output assemblation.
std::shared_ptr< Observe > observe_
Shared Observe object.
void reallocate_cache()
Calls cache_reallocate method on set of used fields.
GenericAssemblyObserve(typename DimAssembly< 1 >::EqFields *eq_fields, const std::unordered_set< string > &observe_fields_list, std::shared_ptr< Observe > observe)
Constructor.
RevertableList< BulkIntegralData > bulk_integral_data_
Holds data for computing bulk integrals.
MixedPtr< DimAssembly, 1 > multidim_assembly_
Assembly object.
MixedPtr< DimAssembly, 1 > multidim_assembly() const
Getter to set of assembly objects.
void assemble(std::shared_ptr< DOFHandlerMultiDim > dh) override
General assemble methods.
ElementAccessor< 3 > element_accessor(unsigned int idx) const
Create and return ElementAccessor to element of given idx.
Definition: mesh.cc:866
Point accessor allow iterate over local Observe points.
Definition: observe.hh:372
unsigned int loc_point_time_index() const
Return local index in data cache (combination of local point index and index of stored time)
Definition: observe.hh:402
OutputDataPtr get_output_cache(std::string field_name)
Definition: observe.hh:300
PatchPointVec & patch_point_data()
Getter of patch_point_data.
Definition: observe.hh:307
Base class for quadrature rules on simplices in arbitrary dimensions.
Definition: quadrature.hh:48
Armor::Array< double >::ArrayMatSet set(uint i)
Definition: quadrature.hh:93
double weight(unsigned int i) const
Returns the ith weight.
Definition: quadrature.hh:103
unsigned int idx() const
Returns a global index of the region.
Definition: region.hh:81
Declaration of class which handles the ordering of degrees of freedom (dof) and mappings between loca...
@ bulk
unsigned int uint
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
std::shared_ptr< BulkIntegral > bulk_
Bulk integrals of elements.
Set of all used integral necessary in assemblation.
std::array< std::shared_ptr< BulkIntegral >, 3 > bulk_
Bulk integrals of elements of dimensions 1, 2, 3.
std::size_t make_permanent()
Finalize temporary part of data.
void reset()
Clear the list.
std::size_t emplace_back(Args &&... args)
std::size_t permanent_size() const
Return permanent size of list.
#define END_TIMER(tag)
Ends a timer with specified tag.
#define START_TIMER(tag)
Starts a timer with specified tag.