Flow123d  master-f44eb46
field_value_cache.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 field_value_cache.hh
15  * @brief
16  * @author David Flanderka
17  */
18 
19 #ifndef FIELD_VALUE_CACHE_HH_
20 #define FIELD_VALUE_CACHE_HH_
21 
22 #include <set>
23 #include <unordered_map>
24 #include <unordered_set>
25 #include <vector>
26 #include "system/armor.hh"
27 #include "fields/eval_points.hh"
28 #include "mesh/accessors.hh"
29 #include "tools/mixed.hh"
30 #include "tools/revertable_list.hh"
31 #include "fem/dofhandler.hh"
32 
33 class EvalPoints;
34 class ElementCacheMap;
35 class DHCellAccessor;
36 class DHCellSide;
37 template < template<IntDim...> class DimAssembly> class GenericAssembly;
38 template < template<IntDim...> class DimAssembly> class GenericAssemblyObserve;
39 
40 
41 /**
42  * @brief Class holds precomputed field values of selected element set.
43  *
44  * Every field in equation use own instance used for elements of all dimensions.
45  */
46 template<class elm_type> using FieldValueCache = Armor::Array<elm_type>;
47 
48 
49 /**
50  * Specifies eval points by idx of region, element and eval point.
51  *
52  * TODO Add better description after finish implementation
53  */
54 struct EvalPointData {
55  EvalPointData() {} ///< Default constructor
56  /// Constructor sets all data members
57  EvalPointData(unsigned int i_reg, unsigned int i_ele, unsigned int i_ep, unsigned int dh_loc_idx)
58  : i_reg_(i_reg), i_element_(i_ele), i_eval_point_(i_ep), dh_loc_idx_(dh_loc_idx) {}
59  /// Copy constructor
62 
63 
64  bool operator < (const EvalPointData &other) {
65  if (i_reg_ == other.i_reg_) {
66  if (i_element_ == other.i_element_)
67  return (i_eval_point_ < other.i_eval_point_);
68  else
69  return (i_element_ < other.i_element_);
70  } else
71  return (i_reg_ < other.i_reg_);
72  }
73 
74  unsigned int i_reg_; ///< region_idx of element
75  unsigned int i_element_; ///< mesh_idx of ElementAccessor appropriate to element
76  unsigned int i_eval_point_; ///< index of point in EvalPoint object
77  unsigned int dh_loc_idx_; ///< local index of cell in DOF handler
78 };
79 
80 
81 
82 
83 /**
84  * @brief Auxiliary data class holds number of elements in cache and allow to set this value
85  * explicitly (e.g. as input parameter).
86  *
87  * Implementation is done as singletone with two access through static methods 'get' and 'set'.
88  *
89  * TODO: This is actually a constant so make it a constant int in element cache map.
90  */
92 public:
93  /// Return number of stored elements
94  static unsigned int get() {
95  return get_instance().n_elem_;
96  }
97 
98  /// Set number of stored elements
99  static void set(unsigned int n_elem) {
100  get_instance().n_elem_ = n_elem;
101  }
102 
103  CacheMapElementNumber(CacheMapElementNumber const&) = delete; ///< We don't need copy constructor.
104  void operator=(CacheMapElementNumber const&) = delete; ///< We don't need assignment operator.
105 
106 private:
107  /// Forbiden default constructor
109 
110 
112  {
113  static CacheMapElementNumber instance;
114  return instance;
115  }
116 
117  /// Maximal number of elements stored in cache.
118  unsigned int n_elem_;
119 };
120 
121 
122 /**
123  * @brief Directing class of FieldValueCache.
124  *
125  * Manage storing and updating element data (elements of same dimension) to cache. We need only
126  * one shared instance of this class for all fields in equation (but typically for dim = 1,2,3).
127  *
128  * IMPORTANT: Because there are combined bulk and boundary elements, we must use mesh_idx value
129  * to correct identification of elements.
130  *
131  * TODO:
132  * 1. Generic assembly pass through the patch collecting needed quadrature points. (PASS ORDER)
133  * 2. Then we sort these points for efficient chae_upadate of the fields (CACHE ORDER)
134  * 3. We pass through the patch again evaluating actual integrals. This second pass is currently inefficient since
135  * we can not map efficiently from the PASS ORDER to the CACHE ORDER), this leads to many complications in
136  * quad point classes.
137  * We should:
138  * 1. have templated patch iteration mechanism, so that we can iterate through the evaluated integrals
139  * twice in consistent way performing:
140  * FIRST: collection of evaluation points
141  * SECOND: evaluation of integrals using the fields and fe_values
142  * Having a consistent implementation allows us to assign unique indices to the integral points on the patch.
143  * 2. Sort collected points, remove duplicities, mark new indices to the original list of points.
144  * 3. SECOND pass, use unique index to find the point in the cache.
145  *
146  * Resulting simplifications:
147  * - no need for associating point operations for Edge, Coupling and Boundary points,
148  * we add assiciated eval point pairs as separate eval points with unique ids. Consistent integral iteration
149  * allows us to simply take two succesive points at the second pass.
150  * - no need for the matrix mapping (element, eval_point) to the cache index
151  */
153 public:
154  /// Index of invalid element in cache.
155  static const unsigned int undef_elem_idx;
156 
157  /// Constructor
158  ElementCacheMap();
159 
160  /// Destructor
162 
163  /// Init cache
164  void init(std::shared_ptr<EvalPoints> eval_points);
165 
166  /// Create patch of cached elements before reading data to cache.
167  void create_patch();
168 
169  /// Reset all items of elements_eval_points_map
172  unsigned int last_element_idx = -1, i_elem_row = -1;
173  for (unsigned int i=0; i<eval_point_data_.permanent_size(); ++i) {
174  if (eval_point_data_[i].i_element_ != last_element_idx) { // new element
175  i_elem_row++;
176  last_element_idx =eval_point_data_[i].i_element_;
177  }
179  }
181  }
182 
183  /// Start update of cache.
184  void start_elements_update();
185 
186  /// Finish update after reading data to cache.
187  void finish_elements_update();
188 
189  /// Getter of eval_points object.
190  inline std::shared_ptr<EvalPoints> eval_points() const {
191  return eval_points_;
192  }
193 
194  /** Adds EvalPointData using emplace_back.
195  * Arguments correspond to constructor of EvalPointData.
196  */
197  inline void add_eval_point(unsigned int i_reg, unsigned int i_ele, unsigned int i_eval_point, unsigned int dh_loc_idx)
198  {
199  eval_point_data_.emplace_back(i_reg, i_ele, i_eval_point, dh_loc_idx);
200  set_of_regions_.insert(i_reg);
201  }
202 
203  /// Returns number of eval. points with addition of max simd duplicates due to regions.
204  inline unsigned int get_simd_rounded_size()
205  {
207  }
208 
209  /*
210  * Access to item of \p element_eval_points_map_ like to two-dimensional array.
211  *
212  * @param i_elem_in_cache idx of ElementAccessor in ElementCacheMap
213  * @param i_eval_point index of local point in EvalPoints
214  * @return index of point in FieldValueCache.
215  */
216  inline int element_eval_point(unsigned int i_elem_in_cache, unsigned int i_eval_point) const {
218  return element_eval_points_map_[i_elem_in_cache*eval_points_->max_size()+i_eval_point];
219  }
220 
221  /// Return mesh_idx of element stored at given position of ElementCacheMap
222  inline unsigned int elm_idx_on_position(unsigned pos) const {
223  return elm_idx_[pos];
224  }
225 
226  /// Return position of element stored in ElementCacheMap
227  inline unsigned int position_in_cache(unsigned mesh_elm_idx, bool bdr=false) const {
228  std::unordered_map<unsigned int, unsigned int>::const_iterator it;
229  if (bdr) {
230  it = element_to_map_bdr_.find(mesh_elm_idx);
231  if ( it != element_to_map_bdr_.end() ) return it->second;
233  } else {
234  it = element_to_map_.find(mesh_elm_idx);
235  if ( it != element_to_map_.end() ) return it->second;
237  }
238  }
239 
240  /// Return number of stored regions.
241  inline unsigned int n_regions() const {
242  return regions_starts_.permanent_size() - 1;
243  }
244 
245  /// Return number of stored elements.
246  inline unsigned int n_elements() const {
247  return element_starts_.permanent_size() - 1;
248  }
249 
250  /// Return begin position of element chunk in FieldValueCache
251  inline unsigned int element_chunk_begin(unsigned int elm_patch_idx) const {
252  ASSERT_LT(elm_patch_idx, n_elements());
253  return element_starts_[elm_patch_idx];
254  }
255 
256  /// Return end position of element chunk in FieldValueCache
257  inline unsigned int element_chunk_end(unsigned int elm_patch_idx) const {
258  ASSERT_LT(elm_patch_idx, n_elements());
259  return element_starts_[elm_patch_idx+1];
260  }
261 
262  /// Return begin position of region chunk in FieldValueCache
263  inline unsigned int region_chunk_begin(unsigned int region_patch_idx) const {
264  ASSERT_LT(region_patch_idx, n_regions());
265  return element_starts_[ regions_starts_[region_patch_idx] ];
266  }
267 
268  /// Return end position of region chunk in FieldValueCache
269  inline unsigned int region_chunk_end(unsigned int region_patch_idx) const {
270  ASSERT_LT(region_patch_idx, n_regions());
271  return element_starts_[ regions_starts_[region_patch_idx+1] ];
272  }
273 
274  /// Return begin position of region chunk specified by position in map
275  inline unsigned int region_chunk_by_map_index(unsigned int r_idx) const {
276  if (r_idx <= n_regions()) return element_starts_[ regions_starts_[r_idx] ];
278  }
279 
280  /// Return begin position of region chunk specified by position in map
281  inline unsigned int region_idx_from_chunk_position(unsigned int chunk_pos) const {
282  return eval_point_data_[ this->region_chunk_by_map_index(chunk_pos) ].i_reg_;
283  }
284 
285  /// Return item of eval_point_data_ specified by its position
286  inline const EvalPointData &eval_point_data(unsigned int point_idx) const {
287  return eval_point_data_[point_idx];
288  }
289 
290  /// Return value of evaluation point given by idx of element in patch and local point idx in EvalPoints from cache.
291  template<class Value>
292  inline typename Value::return_type get_value(const FieldValueCache<typename Value::element_type> &field_cache,
293  unsigned int elem_patch_idx, unsigned int eval_points_idx) const {
294  ASSERT_EQ(Value::NRows_, field_cache.n_rows());
295  ASSERT_EQ(Value::NCols_, field_cache.n_cols());
296  unsigned int value_cache_idx = this->element_eval_point(elem_patch_idx, eval_points_idx);
297  ASSERT(value_cache_idx != ElementCacheMap::undef_elem_idx);
298  return Value::get_from_array(field_cache, value_cache_idx);
299  }
300 
301  /// Size of block (evaluation of FieldFormula) must be multiple of this value.
302  /// TODO We should take this value from BParser and it should be dependent on processor configuration.
303  unsigned int simd_size_double;
304 
305 protected:
306  /// Special constant (@see element_eval_points_map_).
307  static const int unused_point = -1;
308 
309  /// Base number of stored regions in patch
310  static const unsigned int regions_in_chunk = 3;
311 
312  /// Base number of stored elements in patch
313  static const unsigned int elements_in_chunk = 10;
314 
315  /// Set item of \p element_eval_points_map_.
316  inline void set_element_eval_point(unsigned int i_elem_in_cache, unsigned int i_eval_point, int val) const {
318  element_eval_points_map_[i_elem_in_cache*eval_points_->max_size()+i_eval_point] = val;
319  }
320 
321 
322  /// Vector of element indexes stored in cache.
324 
325  /// Pointer to EvalPoints
326  std::shared_ptr<EvalPoints> eval_points_;
327 
328  /// Flag is set down during update of cache when this can't be read
330 
331  /**
332  * This array provides indexes to FieldValueCache.
333  *
334  * This one dimensional array behaves like two dimensional factually.
335  * Size is set to 'n_cached_elements * n_eval_points' and items are
336  * accessible through two indices:
337  *
338  * 1: Over elements holds in ElementCacheMap
339  * 2: Over EvalPoints for each element
340  *
341  * Use always and only methods \p element_eval_point for read and
342  * \p set_element_eval_point (for write) to access to items!
343  *
344  * Array is filled in those three steps:
345  * a. Reset - all items are set to ElementCacheMap::unused_point
346  * b. Used eval points are set to ElementCacheMap::point_in_proggress
347  * c. Eval points marked in previous step are sequentially numbered
348  *
349  * TODO improve description
350  */
352 
353  ///< Holds data of evaluating points in patch.
355 
356  /// @name Holds start positions and orders of region chunks and element chunks
357  // @{
358 
359  RevertableList<unsigned int> regions_starts_; ///< Start positions of elements in regions (size = n_regions+1, last value is end of last region)
360  RevertableList<unsigned int> element_starts_; ///< Start positions of elements in eval_point_data_ (size = n_elements+1)
361  std::unordered_map<unsigned int, unsigned int> element_to_map_; ///< Maps bulk element_idx to element index in patch - TODO remove
362  std::unordered_map<unsigned int, unsigned int> element_to_map_bdr_; ///< Maps boundary element_idx to element index in patch - TODO remove
363 
364  // @}
365 
366  /// Keeps set of unique region indices of added eval. points.
367  std::unordered_set<unsigned int> set_of_regions_;
368 
369  // TODO: remove friend class
370  template < template<IntDim...> class DimAssembly>
371  friend class GenericAssembly;
372  template < template<IntDim...> class DimAssembly>
374 };
375 
376 
377 
378 #endif /* FIELD_VALUE_CACHE_HH_ */
#define ASSERT(expr)
Definition: asserts.hh:351
#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_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:341
uint n_rows() const
Definition: armor.hh:715
uint n_cols() const
Definition: armor.hh:720
Auxiliary data class holds number of elements in cache and allow to set this value explicitly (e....
unsigned int n_elem_
Maximal number of elements stored in cache.
static unsigned int get()
Return number of stored elements.
static CacheMapElementNumber & get_instance()
CacheMapElementNumber()
Forbiden default constructor.
CacheMapElementNumber(CacheMapElementNumber const &)=delete
We don't need copy constructor.
void operator=(CacheMapElementNumber const &)=delete
We don't need assignment operator.
static void set(unsigned int n_elem)
Set number of stored elements.
Cell accessor allow iterate over DOF handler cells.
Side accessor allows to iterate over sides of DOF handler cell.
Directing class of FieldValueCache.
unsigned int region_chunk_end(unsigned int region_patch_idx) const
Return end position of region chunk in FieldValueCache.
void set_element_eval_point(unsigned int i_elem_in_cache, unsigned int i_eval_point, int val) const
Set item of element_eval_points_map_.
unsigned int elm_idx_on_position(unsigned pos) const
Return mesh_idx of element stored at given position of ElementCacheMap.
unsigned int element_chunk_begin(unsigned int elm_patch_idx) const
Return begin position of element chunk in FieldValueCache.
ElementCacheMap()
Constructor.
RevertableList< unsigned int > element_starts_
Start positions of elements in eval_point_data_ (size = n_elements+1)
RevertableList< unsigned int > regions_starts_
Start positions of elements in regions (size = n_regions+1, last value is end of last region)
unsigned int get_simd_rounded_size()
Returns number of eval. points with addition of max simd duplicates due to regions.
std::unordered_map< unsigned int, unsigned int > element_to_map_bdr_
Maps boundary element_idx to element index in patch - TODO remove.
void start_elements_update()
Start update of cache.
void finish_elements_update()
Finish update after reading data to cache.
Value::return_type get_value(const FieldValueCache< typename Value::element_type > &field_cache, unsigned int elem_patch_idx, unsigned int eval_points_idx) const
Return value of evaluation point given by idx of element in patch and local point idx in EvalPoints f...
void clear_element_eval_points_map()
Reset all items of elements_eval_points_map.
unsigned int region_chunk_by_map_index(unsigned int r_idx) const
Return begin position of region chunk specified by position in map.
void create_patch()
Create patch of cached elements before reading data to cache.
std::unordered_set< unsigned int > set_of_regions_
Keeps set of unique region indices of added eval. points.
RevertableList< EvalPointData > eval_point_data_
static const int unused_point
Special constant (.
void add_eval_point(unsigned int i_reg, unsigned int i_ele, unsigned int i_eval_point, unsigned int dh_loc_idx)
unsigned int n_elements() const
Return number of stored elements.
void init(std::shared_ptr< EvalPoints > eval_points)
Init cache.
static const unsigned int undef_elem_idx
Index of invalid element in cache.
static const unsigned int elements_in_chunk
Base number of stored elements in patch.
bool ready_to_reading_
Flag is set down during update of cache when this can't be read.
unsigned int n_regions() const
Return number of stored regions.
static const unsigned int regions_in_chunk
Base number of stored regions in patch.
std::unordered_map< unsigned int, unsigned int > element_to_map_
Maps bulk element_idx to element index in patch - TODO remove.
int element_eval_point(unsigned int i_elem_in_cache, unsigned int i_eval_point) const
std::shared_ptr< EvalPoints > eval_points() const
Getter of eval_points object.
std::shared_ptr< EvalPoints > eval_points_
Pointer to EvalPoints.
unsigned int region_chunk_begin(unsigned int region_patch_idx) const
Return begin position of region chunk in FieldValueCache.
~ElementCacheMap()
Destructor.
unsigned int region_idx_from_chunk_position(unsigned int chunk_pos) const
Return begin position of region chunk specified by position in map.
unsigned int position_in_cache(unsigned mesh_elm_idx, bool bdr=false) const
Return position of element stored in ElementCacheMap.
int * element_eval_points_map_
Holds data of evaluating points in patch.
unsigned int simd_size_double
unsigned int element_chunk_end(unsigned int elm_patch_idx) const
Return end position of element chunk in FieldValueCache.
const EvalPointData & eval_point_data(unsigned int point_idx) const
Return item of eval_point_data_ specified by its position.
std::vector< unsigned int > elm_idx_
Vector of element indexes stored in cache.
Class holds local coordinations of evaluating points (bulk and sides) specified by element dimension.
Definition: eval_points.hh:43
Generic class of observe output assemblation.
Generic class of assemblation.
Declaration of class which handles the ordering of degrees of freedom (dof) and mappings between loca...
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
unsigned int i_eval_point_
index of point in EvalPoint object
unsigned int dh_loc_idx_
local index of cell in DOF handler
unsigned int i_element_
mesh_idx of ElementAccessor appropriate to element
unsigned int i_reg_
region_idx of element
EvalPointData(unsigned int i_reg, unsigned int i_ele, unsigned int i_ep, unsigned int dh_loc_idx)
Constructor sets all data members.
bool operator<(const EvalPointData &other)
EvalPointData(const EvalPointData &other)
Copy constructor.
void reset()
Clear the list.
std::size_t emplace_back(Args &&... args)
std::size_t temporary_size() const
Return temporary size of list (full size of stored data).
std::size_t permanent_size() const
Return permanent size of list.