Flow123d  master-f44eb46
generic_assembly.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 generic_assembly.hh
15  * @brief
16  */
17 
18 #ifndef GENERIC_ASSEMBLY_HH_
19 #define GENERIC_ASSEMBLY_HH_
20 
22 #include "fields/eval_subset.hh"
23 #include "fields/eval_points.hh"
25 #include "tools/revertable_list.hh"
26 #include "system/sys_profiler.hh"
27 
28 
29 
30 /// Allow set mask of active integrals.
32  no_intg = 0,
33  bulk = 0x0001,
34  edge = 0x0002,
35  coupling = 0x0004,
36  boundary = 0x0008
37 };
38 
39 
40 /// Set of all used integral necessary in assemblation
42  std::array<std::shared_ptr<BulkIntegral>, 3> bulk_; ///< Bulk integrals of elements of dimensions 1, 2, 3
43  std::array<std::shared_ptr<EdgeIntegral>, 3> edge_; ///< Edge integrals between elements of dimensions 1, 2, 3
44  std::array<std::shared_ptr<CouplingIntegral>, 2> coupling_; ///< Coupling integrals between elements of dimensions 1-2, 2-3
45  std::array<std::shared_ptr<BoundaryIntegral>, 3> boundary_; ///< Boundary integrals betwwen elements of dimensions 1, 2, 3 and boundaries
46 };
47 
48 
49 /**
50  * Common interface class for all Assembly classes.
51  */
53 {
54 public:
55  /**
56  * Helper structzre holds data of cell (bulk) integral
57  *
58  * Data is specified by cell and subset index in EvalPoint object
59  */
61  /// Default constructor
63 
64  /// Constructor with data mebers initialization
65  BulkIntegralData(DHCellAccessor dhcell, unsigned int subset_idx)
66  : cell(dhcell), subset_index(subset_idx) {}
67 
68  /// Copy constructor
70  : cell(other.cell), subset_index(other.subset_index) {}
71 
72  DHCellAccessor cell; ///< Specified cell (element)
73  unsigned int subset_index; ///< Index (order) of subset in EvalPoints object
74  };
75 
76  /**
77  * Helper structzre holds data of edge integral
78  *
79  * Data is specified by side and subset index in EvalPoint object
80  */
82  /// Default constructor
85 
86  /// Copy constructor
89 
90  /// Constructor with data mebers initialization
92  : edge_side_range(range), subset_index(subset_idx) {}
93 
94  RangeConvert<DHEdgeSide, DHCellSide> edge_side_range; ///< Specified cell side (element)
95  unsigned int subset_index; ///< Index (order) of subset in EvalPoints object
96  };
97 
98  /**
99  * Helper structzre holds data of neighbour (coupling) integral
100  *
101  * Data is specified by cell, side and their subset indices in EvalPoint object
102  */
104  /// Default constructor
106 
107  /// Constructor with data mebers initialization
108  CouplingIntegralData(DHCellAccessor dhcell, unsigned int bulk_idx, DHCellSide dhside, unsigned int side_idx)
109  : cell(dhcell), bulk_subset_index(bulk_idx), side(dhside), side_subset_index(side_idx) {}
110 
111  /// Copy constructor
114 
116  unsigned int bulk_subset_index; ///< Index (order) of lower dim subset in EvalPoints object
117  DHCellSide side; ///< Specified cell side (higher dim element)
118  unsigned int side_subset_index; ///< Index (order) of higher dim subset in EvalPoints object
119  };
120 
121  /**
122  * Helper structzre holds data of boundary integral
123  *
124  * Data is specified by side and subset indices of side and appropriate boundary element in EvalPoint object
125  */
127  /// Default constructor
129 
130  /// Constructor with data mebers initialization
131  BoundaryIntegralData(unsigned int bdr_idx, DHCellSide dhside, unsigned int side_idx)
132  : bdr_subset_index(bdr_idx), side(dhside), side_subset_index(side_idx) {}
133 
134  /// Copy constructor
137 
138  // We don't need hold ElementAccessor of boundary element, side.cond().element_accessor() provides it.
139  unsigned int bdr_subset_index; ///< Index (order) of subset on boundary element in EvalPoints object
140  DHCellSide side; ///< Specified cell side (bulk element)
141  unsigned int side_subset_index; ///< Index (order) of subset on side of bulk element in EvalPoints object
142  };
143 
146  virtual void assemble(std::shared_ptr<DOFHandlerMultiDim> dh) = 0;
147 
148 protected:
149  AssemblyIntegrals integrals_; ///< Holds integral objects.
150  std::shared_ptr<EvalPoints> eval_points_; ///< EvalPoints object shared by all integrals
151  ElementCacheMap element_cache_map_; ///< ElementCacheMap according to EvalPoints
152 };
153 
154 
155 /**
156  * @brief Generic class of assemblation.
157  *
158  * Class
159  * - holds assemblation structures (EvalPoints, Integral objects, Integral data table).
160  * - associates assemblation objects specified by dimension
161  * - provides general assemble method
162  * - provides methods that allow construction of element patches
163  */
164 template < template<IntDim...> class DimAssembly>
166 {
167 public:
168  /// Constructor
169  GenericAssembly( typename DimAssembly<1>::EqFields *eq_fields, typename DimAssembly<1>::EqData *eq_data)
170  : multidim_assembly_(eq_fields, eq_data),
171  min_edge_sides_(2),
172  bulk_integral_data_(20, 10),
173  edge_integral_data_(12, 6),
176  {
177  eval_points_ = std::make_shared<EvalPoints>();
178  // first step - create integrals, then - initialize cache and initialize subobject of dimensions
179  multidim_assembly_[1_d]->create_integrals(eval_points_, integrals_);
180  multidim_assembly_[2_d]->create_integrals(eval_points_, integrals_);
181  multidim_assembly_[3_d]->create_integrals(eval_points_, integrals_);
183  multidim_assembly_[1_d]->initialize(&element_cache_map_);
184  multidim_assembly_[2_d]->initialize(&element_cache_map_);
185  multidim_assembly_[3_d]->initialize(&element_cache_map_);
186  active_integrals_ = multidim_assembly_[1_d]->n_active_integrals();
187  }
188 
189  /// Getter to set of assembly objects
191  return multidim_assembly_;
192  }
193 
194  /// Geter to EvalPoints object
195  inline std::shared_ptr<EvalPoints> eval_points() const {
196  return eval_points_;
197  }
198 
199  void set_min_edge_sides(unsigned int val) {
200  min_edge_sides_ = val;
201  }
202 
203  /**
204  * @brief General assemble methods.
205  *
206  * Loops through local cells and calls assemble methods of assembly
207  * object of each cells over space dimension.
208  */
209  void assemble(std::shared_ptr<DOFHandlerMultiDim> dh) override {
210  START_TIMER( DimAssembly<1>::name() );
211  this->reallocate_cache();
212  multidim_assembly_[1_d]->begin();
213 
214  bool add_into_patch = false; // control variable
215  for(auto cell_it = dh->local_range().begin(); cell_it != dh->local_range().end(); )
216  {
217 
218  if (!add_into_patch) {
220  add_into_patch = true;
221  }
222 
223  //START_TIMER("add_integrals_to_patch");
224  this->add_integrals_of_computing_step(*cell_it);
225  //END_TIMER("add_integrals_to_patch");
226 
233  this->assemble_integrals();
234  add_into_patch = false;
235  } else {
242  this->assemble_integrals();
243  add_into_patch = false;
244  }
245  ++cell_it;
246  }
247  }
248  if (add_into_patch) {
249  this->assemble_integrals();
250  }
251 
252  multidim_assembly_[1_d]->end();
253  END_TIMER( DimAssembly<1>::name() );
254  }
255 
256  /// Return ElementCacheMap
257  inline const ElementCacheMap &cache_map() const {
258  return element_cache_map_;
259  }
260 
261 private:
262  /// Call assemblations when patch is filled
264  START_TIMER("create_patch");
266  END_TIMER("create_patch");
267  START_TIMER("cache_update");
268  multidim_assembly_[1_d]->eq_fields_->cache_update(element_cache_map_); // TODO replace with sub FieldSet
269  END_TIMER("cache_update");
271 
272  {
273  START_TIMER("assemble_volume_integrals");
274  multidim_assembly_[1_d]->assemble_cell_integrals(bulk_integral_data_);
275  multidim_assembly_[2_d]->assemble_cell_integrals(bulk_integral_data_);
276  multidim_assembly_[3_d]->assemble_cell_integrals(bulk_integral_data_);
277  END_TIMER("assemble_volume_integrals");
278  }
279 
280  {
281  START_TIMER("assemble_fluxes_boundary");
282  multidim_assembly_[1_d]->assemble_boundary_side_integrals(boundary_integral_data_);
283  multidim_assembly_[2_d]->assemble_boundary_side_integrals(boundary_integral_data_);
284  multidim_assembly_[3_d]->assemble_boundary_side_integrals(boundary_integral_data_);
285  END_TIMER("assemble_fluxes_boundary");
286  }
287 
288  {
289  START_TIMER("assemble_fluxes_elem_elem");
290  multidim_assembly_[1_d]->assemble_edge_integrals(edge_integral_data_);
291  multidim_assembly_[2_d]->assemble_edge_integrals(edge_integral_data_);
292  multidim_assembly_[3_d]->assemble_edge_integrals(edge_integral_data_);
293  END_TIMER("assemble_fluxes_elem_elem");
294  }
295 
296  {
297  START_TIMER("assemble_fluxes_elem_side");
298  multidim_assembly_[2_d]->assemble_neighbour_integrals(coupling_integral_data_);
299  multidim_assembly_[3_d]->assemble_neighbour_integrals(coupling_integral_data_);
300  END_TIMER("assemble_fluxes_elem_side");
301  }
302  // clean integral data
308  }
309 
310  /**
311  * Add data of integrals to appropriate structure and register elements to ElementCacheMap.
312  *
313  * Types of used integrals must be set in data member \p active_integrals_.
314  */
317  if (cell.is_own()) { // Not ghost
318  this->add_volume_integral(cell);
319  }
320 
321  for( DHCellSide cell_side : cell.side_range() ) {
323  if (cell.is_own()) // Not ghost
324  if ( (cell_side.side().edge().n_sides() == 1) && (cell_side.side().is_boundary()) ) {
325  this->add_boundary_integral(cell_side);
326  continue;
327  }
329  if ( (cell_side.n_edge_sides() >= min_edge_sides_) && (cell_side.edge_sides().begin()->element().idx() == cell.elm_idx())) {
330  this->add_edge_integral(cell_side);
331  }
332  }
333 
335  bool add_low = true;
336  for( DHCellSide neighb_side : cell.neighb_sides() ) { // cell -> elm lower dim, neighb_side -> elm higher dim
337  if (cell.dim() != neighb_side.dim()-1) continue;
338  this->add_coupling_integral(cell, neighb_side, add_low);
339  add_low = false;
340  }
341  }
342  }
343 
344  /// Add data of volume integral to appropriate data structure.
345  inline void add_volume_integral(const DHCellAccessor &cell) {
346  uint subset_idx = integrals_.bulk_[cell.dim()-1]->get_subset_idx();
347  bulk_integral_data_.emplace_back(cell, subset_idx);
348 
349  unsigned int reg_idx = cell.elm().region_idx().idx();
350  // Different access than in other integrals: We can't use range method CellIntegral::points
351  // because it passes element_patch_idx as argument that is not known during patch construction.
352  for (uint i=uint( eval_points_->subset_begin(cell.dim(), subset_idx) );
353  i<uint( eval_points_->subset_end(cell.dim(), subset_idx) ); ++i) {
354  element_cache_map_.add_eval_point(reg_idx, cell.elm_idx(), i, cell.local_idx());
355  }
356  }
357 
358  /// Add data of edge integral to appropriate data structure.
359  inline void add_edge_integral(const DHCellSide &cell_side) {
360  auto range = cell_side.edge_sides();
361  edge_integral_data_.emplace_back(range, integrals_.edge_[range.begin()->dim()-1]->get_subset_idx());
362 
363  for( DHCellSide edge_side : range ) {
364  unsigned int reg_idx = edge_side.element().region_idx().idx();
365  for (auto p : integrals_.edge_[range.begin()->dim()-1]->points(edge_side, &element_cache_map_) ) {
366  element_cache_map_.add_eval_point(reg_idx, edge_side.elem_idx(), p.eval_point_idx(), edge_side.cell().local_idx());
367  }
368  }
369  }
370 
371  /// Add data of coupling integral to appropriate data structure.
372  inline void add_coupling_integral(const DHCellAccessor &cell, const DHCellSide &ngh_side, bool add_low) {
373  coupling_integral_data_.emplace_back(cell, integrals_.coupling_[cell.dim()-1]->get_subset_low_idx(), ngh_side,
374  integrals_.coupling_[cell.dim()-1]->get_subset_high_idx());
375 
376  unsigned int reg_idx_low = cell.elm().region_idx().idx();
377  unsigned int reg_idx_high = ngh_side.element().region_idx().idx();
378  for (auto p : integrals_.coupling_[cell.dim()-1]->points(ngh_side, &element_cache_map_) ) {
379  element_cache_map_.add_eval_point(reg_idx_high, ngh_side.elem_idx(), p.eval_point_idx(), ngh_side.cell().local_idx());
380 
381  if (add_low) {
382  auto p_low = p.lower_dim(cell); // equivalent point on low dim cell
383  element_cache_map_.add_eval_point(reg_idx_low, cell.elm_idx(), p_low.eval_point_idx(), cell.local_idx());
384  }
385  }
386  }
387 
388  /// Add data of boundary integral to appropriate data structure.
389  inline void add_boundary_integral(const DHCellSide &bdr_side) {
390  boundary_integral_data_.emplace_back(integrals_.boundary_[bdr_side.dim()-1]->get_subset_low_idx(), bdr_side,
391  integrals_.boundary_[bdr_side.dim()-1]->get_subset_high_idx());
392 
393  unsigned int reg_idx = bdr_side.element().region_idx().idx();
394  for (auto p : integrals_.boundary_[bdr_side.dim()-1]->points(bdr_side, &element_cache_map_) ) {
395  element_cache_map_.add_eval_point(reg_idx, bdr_side.elem_idx(), p.eval_point_idx(), bdr_side.cell().local_idx());
396 
397  BulkPoint p_bdr = p.point_bdr(bdr_side.cond().element_accessor()); // equivalent point on boundary element
398  unsigned int bdr_reg = bdr_side.cond().element_accessor().region_idx().idx();
399  // invalid local_idx value, DHCellAccessor of boundary element doesn't exist
400  element_cache_map_.add_eval_point(bdr_reg, bdr_side.cond().bc_ele_idx(), p_bdr.eval_point_idx(), -1);
401  }
402  }
403 
404  /// Calls cache_reallocate method on
405  inline void reallocate_cache() {
406  multidim_assembly_[1_d]->eq_fields_->cache_reallocate(this->element_cache_map_, multidim_assembly_[1_d]->used_fields_);
407  // DebugOut() << "Order of evaluated fields (" << DimAssembly<1>::name() << "):" << multidim_assembly_[1_d]->eq_fields_->print_dependency();
408  }
409 
410 
411  /// Assembly object
413 
414  /// Holds mask of active integrals.
416 
417  /**
418  * Minimal number of sides on edge.
419  *
420  * Edge integral is created and calculated if number of sides is greater or equal than this value. Default value
421  * is 2 and can be changed
422  */
423  unsigned int min_edge_sides_;
424 
425  // Following variables hold data of all integrals depending of actual computed element.
426  // TODO sizes of arrays should be set dynamically, depend on number of elements in ElementCacheMap,
427  RevertableList<BulkIntegralData> bulk_integral_data_; ///< Holds data for computing bulk integrals.
428  RevertableList<EdgeIntegralData> edge_integral_data_; ///< Holds data for computing edge integrals.
429  RevertableList<CouplingIntegralData> coupling_integral_data_; ///< Holds data for computing couplings integrals.
430  RevertableList<BoundaryIntegralData> boundary_integral_data_; ///< Holds data for computing boundary integrals.
431 };
432 
433 
434 #endif /* GENERIC_ASSEMBLY_HH_ */
uint bc_ele_idx()
Definition: accessors.hh:374
ElementAccessor< 3 > element_accessor()
Base point accessor class.
Definition: eval_subset.hh:55
unsigned int eval_point_idx() const
Return index in EvalPoints object.
Definition: eval_subset.hh:84
static unsigned int get()
Return number of stored elements.
Cell accessor allow iterate over DOF handler cells.
const ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
RangeConvert< DHNeighbSide, DHCellSide > neighb_sides() const
Returns range of neighbour cell of lower dimension corresponding to cell of higher dimension.
bool is_own() const
Return true if accessor represents own element (false for ghost element)
unsigned int dim() const
Return dimension of element appropriate to cell.
Range< DHCellSide > side_range() const
Returns range of cell sides.
unsigned int elm_idx() const
Return serial idx to element of loc_ele_idx_.
unsigned int local_idx() const
Return local index to element (index of DOF handler).
Side accessor allows to iterate over sides of DOF handler cell.
unsigned int elem_idx() const
Boundary cond() const
const DHCellAccessor & cell() const
Return DHCellAccessor appropriate to the side.
RangeConvert< DHEdgeSide, DHCellSide > edge_sides() const
Returns range of all sides looped over common Edge.
unsigned int dim() const
Return dimension of element appropriate to the side.
ElementAccessor< 3 > element() const
Class allows to iterate over sides of edge.
RegionIdx region_idx() const
Definition: accessors.hh:201
Directing class of FieldValueCache.
unsigned int get_simd_rounded_size()
Returns number of eval. points with addition of max simd duplicates due to regions.
void start_elements_update()
Start update of cache.
void finish_elements_update()
Finish update after reading data to cache.
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 add_eval_point(unsigned int i_reg, unsigned int i_ele, unsigned int i_eval_point, unsigned int dh_loc_idx)
void init(std::shared_ptr< EvalPoints > eval_points)
Init cache.
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.
virtual void assemble(std::shared_ptr< DOFHandlerMultiDim > dh)=0
Generic class of assemblation.
MixedPtr< DimAssembly, 1 > multidim_assembly_
Assembly object.
RevertableList< EdgeIntegralData > edge_integral_data_
Holds data for computing edge integrals.
void assemble_integrals()
Call assemblations when patch is filled.
void assemble(std::shared_ptr< DOFHandlerMultiDim > dh) override
General assemble methods.
unsigned int min_edge_sides_
void add_boundary_integral(const DHCellSide &bdr_side)
Add data of boundary integral to appropriate data structure.
const ElementCacheMap & cache_map() const
Return ElementCacheMap.
void add_volume_integral(const DHCellAccessor &cell)
Add data of volume integral to appropriate data structure.
void add_edge_integral(const DHCellSide &cell_side)
Add data of edge integral to appropriate data structure.
RevertableList< BoundaryIntegralData > boundary_integral_data_
Holds data for computing boundary integrals.
int active_integrals_
Holds mask of active integrals.
RevertableList< CouplingIntegralData > coupling_integral_data_
Holds data for computing couplings integrals.
void set_min_edge_sides(unsigned int val)
std::shared_ptr< EvalPoints > eval_points() const
Geter to EvalPoints object.
void reallocate_cache()
Calls cache_reallocate method on.
GenericAssembly(typename DimAssembly< 1 >::EqFields *eq_fields, typename DimAssembly< 1 >::EqData *eq_data)
Constructor.
void add_coupling_integral(const DHCellAccessor &cell, const DHCellSide &ngh_side, bool add_low)
Add data of coupling integral to appropriate data structure.
void add_integrals_of_computing_step(DHCellAccessor cell)
MixedPtr< DimAssembly, 1 > multidim_assembly() const
Getter to set of assembly objects.
RevertableList< BulkIntegralData > bulk_integral_data_
Holds data for computing bulk integrals.
unsigned int idx() const
Returns a global index of the region.
Definition: region.hh:81
Iter< Object > make_iter(Object obj)
ActiveIntegrals
Allow set mask of active integrals.
@ coupling
@ boundary
@ no_intg
@ bulk
@ edge
unsigned int uint
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
Definitions of particular quadrature rules on simplices.
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::array< std::shared_ptr< CouplingIntegral >, 2 > coupling_
Coupling integrals between elements of dimensions 1-2, 2-3.
std::array< std::shared_ptr< BoundaryIntegral >, 3 > boundary_
Boundary integrals betwwen elements of dimensions 1, 2, 3 and boundaries.
std::array< std::shared_ptr< EdgeIntegral >, 3 > edge_
Edge integrals between elements of dimensions 1, 2, 3.
unsigned int side_subset_index
Index (order) of subset on side of bulk element in EvalPoints object.
BoundaryIntegralData(unsigned int bdr_idx, DHCellSide dhside, unsigned int side_idx)
Constructor with data mebers initialization.
BoundaryIntegralData(const BoundaryIntegralData &other)
Copy constructor.
unsigned int bdr_subset_index
Index (order) of subset on boundary element in EvalPoints object.
DHCellSide side
Specified cell side (bulk element)
BulkIntegralData(DHCellAccessor dhcell, unsigned int subset_idx)
Constructor with data mebers initialization.
DHCellAccessor cell
Specified cell (element)
unsigned int subset_index
Index (order) of subset in EvalPoints object.
BulkIntegralData(const BulkIntegralData &other)
Copy constructor.
unsigned int bulk_subset_index
Index (order) of lower dim subset in EvalPoints object.
CouplingIntegralData(DHCellAccessor dhcell, unsigned int bulk_idx, DHCellSide dhside, unsigned int side_idx)
Constructor with data mebers initialization.
unsigned int side_subset_index
Index (order) of higher dim subset in EvalPoints object.
CouplingIntegralData(const CouplingIntegralData &other)
Copy constructor.
DHCellSide side
Specified cell side (higher dim element)
EdgeIntegralData(const EdgeIntegralData &other)
Copy constructor.
RangeConvert< DHEdgeSide, DHCellSide > edge_side_range
Specified cell side (element)
unsigned int subset_index
Index (order) of subset in EvalPoints object.
EdgeIntegralData(RangeConvert< DHEdgeSide, DHCellSide > range, unsigned int subset_idx)
Constructor with data mebers initialization.
std::size_t revert_temporary()
Erase temporary part of data.
std::size_t make_permanent()
Finalize temporary part of data.
void reset()
Clear the list.
std::size_t emplace_back(Args &&... args)
#define END_TIMER(tag)
Ends a timer with specified tag.
#define START_TIMER(tag)
Starts a timer with specified tag.