Flow123d  master-f44eb46
field_fe.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 field_fe.cc
15  * @brief
16  */
17 
18 
19 #include <limits>
20 
21 #include "fields/field_fe.hh"
22 #include "la/vector_mpi.hh"
23 #include "fields/field_instances.hh" // for instantiation macros
24 #include "input/input_type.hh"
25 #include "fem/fe_p.hh"
26 #include "fem/fe_system.hh"
27 #include "fem/dh_cell_accessor.hh"
28 #include "fem/mapping_p1.hh"
29 #include "io/reader_cache.hh"
30 #include "io/msh_gmshreader.h"
31 #include "mesh/accessors.hh"
32 #include "mesh/range_wrapper.hh"
33 #include "mesh/bc_mesh.hh"
35 
36 #include "system/sys_profiler.hh"
37 #include "tools/unit_converter.hh"
41 
42 
43 
44 
45 /// Implementation.
46 
47 namespace it = Input::Type;
48 
49 
50 
51 
53 
54 
55 
56 /************************************************************************************
57  * Implementation of FieldFE methods
58  */
59 template <int spacedim, class Value>
61 {
62  return it::Record("FieldFE", FieldAlgorithmBase<spacedim,Value>::template_name()+" Field given by finite element approximation.")
66  "GMSH mesh with data. Can be different from actual computational mesh.")
68  "Section where to find the field.\n Some sections are specific to file format: "
69  "point_data/node_data, cell_data/element_data, -/element_node_data, native/-.\n"
70  "If not given by a user, we try to find the field in all sections, but we report an error "
71  "if it is found in more than one section.")
73  "The values of the Field are read from the ```$ElementData``` section with field name given by this key.")
74  .declare_key("default_value", IT::Double(), IT::Default::optional(),
75  "Default value is set on elements which values have not been listed in the mesh data file.")
77  "Definition of the unit of all times defined in the mesh data file.")
78  .declare_key("read_time_shift", TimeGovernor::get_input_time_type(), IT::Default("0.0"),
79  "This key allows reading field data from the mesh data file shifted in time. Considering the time 't', field descriptor with time 'T', "
80  "time shift 'S', then if 't > T', we read the time frame 't + S'.")
82  IT::Default("\"equivalent_mesh\""), "Type of interpolation applied to the input spatial data.\n"
83  "The default value 'equivalent_mesh' assumes the data being constant on elements living on the same mesh "
84  "as the computational mesh, but possibly with different numbering. In the case of the same numbering, "
85  "the user can set 'identical_mesh' to omit algorithm for guessing node and element renumbering. "
86  "Alternatively, in case of different input mesh, several interpolation algorithms are available.")
87  .declare_key("is_boundary", IT::Bool(), IT::Default("false"),
88  "Distinguishes bulk / boundary FieldFE.")
89  .close();
90 }
91 
92 template <int spacedim, class Value>
94 {
95  return it::Selection("FE_discretization",
96  "Specify the section in mesh input file where field data is listed.\nSome sections are specific to file format.")
97  //.add_value(OutputTime::DiscreteSpace::NODE_DATA, "node_data", "point_data (VTK) / node_data (GMSH)")
98  .add_value(OutputTime::DiscreteSpace::ELEM_DATA, "element_data", "cell_data (VTK) / element_data (GMSH)")
99  //.add_value(OutputTime::DiscreteSpace::CORNER_DATA, "element_node_data", "element_node_data (only for GMSH)")
100  .add_value(OutputTime::DiscreteSpace::NATIVE_DATA, "native_data", "native_data (only for VTK)")
101  .close();
102 }
103 
104 template <int spacedim, class Value>
106 {
107  return it::Selection("interpolation", "Specify interpolation of the input data from its input mesh to the computational mesh.")
108  .add_value(DataInterpolation::identic_msh, "identic_mesh", "Topology and indices of nodes and elements of"
109  "the input mesh and the computational mesh are identical. "
110  "This interpolation is typically used for GMSH input files containing only the field values without "
111  "explicit mesh specification.")
112  .add_value(DataInterpolation::equivalent_msh, "equivalent_mesh", "Topologies of the input mesh and the computational mesh "
113  "are the same, the node and element numbering may differ. "
114  "This interpolation can be used also for VTK input data.") // default value
115  .add_value(DataInterpolation::gauss_p0, "P0_gauss", "Topologies of the input mesh and the computational mesh may differ. "
116  "Constant values on the elements of the computational mesh are evaluated using the Gaussian quadrature of the fixed order 4, "
117  "where the quadrature points and their values are found in the input mesh and input data using the BIH tree search."
118  )
119  .add_value(DataInterpolation::interp_p0, "P0_intersection", "Topologies of the input mesh and the computational mesh may differ. "
120  "Can be applied only for boundary fields. For every (boundary) element of the computational mesh the "
121  "intersection with the input mesh is computed. Constant values on the elements of the computational mesh "
122  "are evaluated as the weighted average of the (constant) values on the intersecting elements of the input mesh.")
123  .close();
124 }
125 
126 template <int spacedim, class Value>
128  Input::register_class< FieldFE<spacedim, Value>, unsigned int >("FieldFE") +
130 
131 
132 
133 template <int spacedim, class Value>
135 : FieldAlgorithmBase<spacedim, Value>(n_comp),
136  dh_(nullptr), field_name_(""), discretization_(OutputTime::DiscreteSpace::UNDEFINED),
137  boundary_domain_(false), fe_values_(4)
138 {
139  this->is_constant_in_space_ = false;
140 }
141 
142 
143 template<int spacedim, class Value>
145  Input::Array multifield_arr;
146  if (rec.opt_val(field.input_name(), multifield_arr))
147  {
148  unsigned int position = 0;
149  auto it = multifield_arr.begin<Input::AbstractRecord>();
150  if (multifield_arr.size() > 1)
151  while (index_ != position) {
152  ++it; ++position;
153  }
154 
155  Input::Record field_rec = (Input::Record)(*it);
156  if (field_rec.val<std::string>("TYPE") == "FieldFE") {
157  OutputTime::DiscreteSpace discretization;
158  if (field_rec.opt_val<OutputTime::DiscreteSpace>("input_discretization", discretization)) {
159  if (discretization == OutputTime::DiscreteSpace::NATIVE_DATA) {
160  std::shared_ptr< FieldFE<spacedim, Value> > field_fe = std::make_shared< FieldFE<spacedim, Value> >(field.n_comp());
161  FieldAlgoBaseInitData init_data(field.input_name(), field.n_comp(), field.units(), field.limits(), field.get_flags());
162  field_fe->init_from_input( *it, init_data );
163  field_fe->set_fe_data(conc_dof_handler_, dof_vector_);
164  return field_fe;
165  }
166  }
167  }
168  }
169 
170  return NULL;
171 }
172 
173 
174 template <int spacedim, class Value>
175 VectorMPI FieldFE<spacedim, Value>::set_fe_data(std::shared_ptr<DOFHandlerMultiDim> dh, VectorMPI dof_values, unsigned int block_index)
176 {
177  dh_ = dh;
178  if (dof_values.size()==0) { //create data vector according to dof handler - Warning not tested yet
179  data_vec_ = dh_->create_vector();
181  } else {
182  data_vec_ = dof_values;
183  }
184 
185  if ( block_index == FieldFE<spacedim, Value>::undef_uint ) {
186  this->fill_fe_item<0>();
187  this->fill_fe_item<1>();
188  this->fill_fe_item<2>();
189  this->fill_fe_item<3>();
190  this->fe_ = dh_->ds()->fe();
191  } else {
192  this->fill_fe_system_data<0>(block_index);
193  this->fill_fe_system_data<1>(block_index);
194  this->fill_fe_system_data<2>(block_index);
195  this->fill_fe_system_data<3>(block_index);
197  std::dynamic_pointer_cast<FESystem<0>>( dh_->ds()->fe()[Dim<0>{}] )->fe()[block_index],
198  std::dynamic_pointer_cast<FESystem<1>>( dh_->ds()->fe()[Dim<1>{}] )->fe()[block_index],
199  std::dynamic_pointer_cast<FESystem<2>>( dh_->ds()->fe()[Dim<2>{}] )->fe()[block_index],
200  std::dynamic_pointer_cast<FESystem<3>>( dh_->ds()->fe()[Dim<3>{}] )->fe()[block_index]
201  );
202  }
203 
204  // set interpolation
205  interpolation_ = DataInterpolation::equivalent_msh;
206 
207  region_value_err_.resize(dh_->mesh()->region_db().size());
208 
209  return data_vec_;
210 }
211 
212 
213 
214 template <int spacedim, class Value>
216  ElementCacheMap &cache_map, unsigned int region_patch_idx)
217 {
218  auto region_idx = cache_map.region_idx_from_chunk_position(region_patch_idx);
219  if ( (region_idx % 2) == this->boundary_domain_ ) {
220  // Skip evaluation of boundary fields on bulk regions and vice versa
221  return;
222  }
223 
225 
226  unsigned int reg_chunk_begin = cache_map.region_chunk_begin(region_patch_idx);
227  unsigned int reg_chunk_end = cache_map.region_chunk_end(region_patch_idx);
228  unsigned int last_element_idx = -1;
229  DHCellAccessor cell = *( dh_->local_range().begin() ); //needs set variable for correct compiling
230  LocDofVec loc_dofs;
231  unsigned int range_bgn=0, range_end=0;
232 
233  // Throws exception if any element value of processed region is NaN
234  unsigned int r_idx = cache_map.eval_point_data(reg_chunk_begin).i_reg_;
235  if (region_value_err_[r_idx].is_invalid_)
236  THROW( ExcUndefElementValue() << EI_Field(field_name_) << EI_File(reader_file_.filename()) );
237 
238  for (unsigned int i_data = reg_chunk_begin; i_data < reg_chunk_end; ++i_data) { // i_eval_point_data
239  unsigned int elm_idx = cache_map.eval_point_data(i_data).i_element_;
240  if (elm_idx != last_element_idx) {
241  ElementAccessor<spacedim> elm(dh_->mesh(), elm_idx);
242  fe_values_[elm.dim()].reinit( elm );
243  cell = dh_->cell_accessor_from_element( elm_idx );
244  loc_dofs = cell.get_loc_dof_indices();
245  last_element_idx = elm_idx;
246  range_bgn = this->fe_item_[elm.dim()].range_begin_;
247  range_end = this->fe_item_[elm.dim()].range_end_;
248  }
249 
250  unsigned int i_ep=cache_map.eval_point_data(i_data).i_eval_point_;
251  //DHCellAccessor cache_cell = cache_map(cell);
252  mat_value.fill(0.0);
253  for (unsigned int i_dof=range_bgn, i_cdof=0; i_dof<range_end; i_dof++, i_cdof++) {
254  mat_value += data_vec_.get(loc_dofs[i_dof]) * this->handle_fe_shape(cell.dim(), i_cdof, i_ep);
255  }
256  data_cache.set(i_data) = mat_value;
257  }
258 }
259 
260 
261 template <int spacedim, class Value>
263 {
264  std::shared_ptr<EvalPoints> eval_points = cache_map.eval_points();
265  std::array<Quadrature, 4> quads{QGauss(0, 1), this->init_quad<1>(eval_points), this->init_quad<2>(eval_points), this->init_quad<3>(eval_points)};
266  fe_values_[0].initialize(quads[0], *this->fe_[0_d], update_values);
267  fe_values_[1].initialize(quads[1], *this->fe_[1_d], update_values);
268  fe_values_[2].initialize(quads[2], *this->fe_[2_d], update_values);
269  fe_values_[3].initialize(quads[3], *this->fe_[3_d], update_values);
270 }
271 
272 
273 template <int spacedim, class Value>
274 template <unsigned int dim>
275 Quadrature FieldFE<spacedim, Value>::init_quad(std::shared_ptr<EvalPoints> eval_points)
276 {
277  Quadrature quad(dim, eval_points->size(dim));
278  for (unsigned int k=0; k<eval_points->size(dim); k++)
279  quad.set(k) = eval_points->local_point<dim>(k);
280  return quad;
281 }
282 
283 
284 template <int spacedim, class Value>
286  this->init_unit_conversion_coefficient(rec, init_data);
287  this->in_rec_ = rec;
288  flags_ = init_data.flags_;
289 
290 
291  // read data from input record
292  reader_file_ = FilePath( rec.val<FilePath>("mesh_data_file") );
293  field_name_ = rec.val<std::string>("field_name");
294  this->boundary_domain_ = rec.val<bool>("is_boundary");
295  if (! rec.opt_val<OutputTime::DiscreteSpace>("input_discretization", discretization_) ) {
296  discretization_ = OutputTime::DiscreteSpace::UNDEFINED;
297  }
298  if (! rec.opt_val<DataInterpolation>("interpolation", interpolation_) ) {
299  interpolation_ = DataInterpolation::equivalent_msh;
300  }
301  if (! rec.opt_val("default_value", default_value_) ) {
302  default_value_ = numeric_limits<double>::signaling_NaN();
303  }
304 }
305 
306 
307 
308 template <int spacedim, class Value>
310  // Mesh can be set only for field initialized from input.
312  ASSERT(field_name_ != "").error("Uninitialized FieldFE, did you call init_from_input()?\n");
313  if (this->interpolation_ == DataInterpolation::identic_msh) {
314  //DebugOut() << "Identic mesh branch\n";
316  } else {
317  //auto source_mesh = ReaderCache::get_mesh(reader_file_);
318  // TODO: move the call into equivalent_mesh_map, get rd of get_element_ids method.
320  if (this->interpolation_ == DataInterpolation::equivalent_msh) {
321  if (source_target_mesh_elm_map_->empty()) { // incompatible meshes
322  this->interpolation_ = DataInterpolation::gauss_p0;
323  WarningOut().fmt("Source mesh of FieldFE '{}' is not compatible with target mesh.\nInterpolation of input data will be changed to 'P0_gauss'.\n",
324  field_name_);
325  }
326  } else if (this->interpolation_ == DataInterpolation::interp_p0) {
327  if (!this->boundary_domain_) {
328  this->interpolation_ = DataInterpolation::gauss_p0;
329  WarningOut().fmt("Interpolation 'P0_intersection' of FieldFE '{}' can't be used on bulk region.\nIt will be changed to 'P0_gauss'.\n",
330  field_name_);
331  }
332  }
333  }
334  if (dh_ == nullptr) {
335  if (this->boundary_domain_) this->make_dof_handler( mesh->bc_mesh() );
336  else this->make_dof_handler( mesh );
337  }
338  region_value_err_.resize(mesh->region_db().size());
339  }
340 }
341 
342 
343 
344 template <int spacedim, class Value>
346 
347  // temporary solution - these objects will be set through FieldCommon
349  switch (this->value_.n_rows() * this->value_.n_cols()) { // by number of components
350  case 1: { // scalar
351  fe = MixedPtr<FE_P_disc>(0);
352  break;
353  }
354  case 3: { // vector
355  MixedPtr<FE_P_disc> fe_base(0) ;
356  fe = mixed_fe_system(fe_base, FEType::FEVector, 3);
357  break;
358  }
359  case 9: { // tensor
360  MixedPtr<FE_P_disc> fe_base(0) ;
361  fe = mixed_fe_system(fe_base, FEType::FETensor, 9);
362  break;
363  }
364  default:
365  ASSERT_PERMANENT(false).error("Should not happen!\n");
366  }
367 
368  std::shared_ptr<DOFHandlerMultiDim> dh_par = std::make_shared<DOFHandlerMultiDim>( const_cast<MeshBase &>(*mesh) );
369  std::shared_ptr<DiscreteSpace> ds = std::make_shared<EqualOrderDiscreteSpace>( &const_cast<MeshBase &>(*mesh), fe);
370  dh_par->distribute_dofs(ds);
371  dh_ = dh_par;
372 
373  this->fill_fe_item<0>();
374  this->fill_fe_item<1>();
375  this->fill_fe_item<2>();
376  this->fill_fe_item<3>();
377  this->fe_ = dh_->ds()->fe();
378 
379  data_vec_ = VectorMPI::sequential( dh_->lsize() ); // allocate data_vec_
380 }
381 
382 
383 
384 template <int spacedim, class Value>
386  // Time can be set only for field initialized from input.
388  ASSERT(field_name_ != "").error("Uninitialized FieldFE, did you call init_from_input()?\n");
389  ASSERT_PTR(dh_)(field_name_).error("Null target mesh pointer of finite element field, did you call set_mesh()?\n");
390  if ( reader_file_ == FilePath() ) return false;
391 
392  unsigned int n_components = this->value_.n_rows() * this->value_.n_cols();
393  double time_unit_coef = time.read_coef(in_rec_.find<Input::Record>("time_unit"));
394  double time_shift = time.read_time( in_rec_.find<Input::Tuple>("read_time_shift") );
395  double read_time = (time.end()+time_shift) / time_unit_coef;
396 
397  unsigned int n_entities, bdr_shift;
398  bool is_native = (this->discretization_ == OutputTime::DiscreteSpace::NATIVE_DATA);
399  if (is_native) {
400  n_components *= dh_->max_elem_dofs();
401  }
402  if (this->interpolation_==DataInterpolation::identic_msh) {
403  n_entities = source_target_mesh_elm_map_->bulk.size() + source_target_mesh_elm_map_->boundary.size();
404  bdr_shift = source_target_mesh_elm_map_->bulk.size();
405  } else {
406  auto reader_mesh = ReaderCache::get_mesh(reader_file_);
407  n_entities = reader_mesh->n_elements() + reader_mesh->bc_mesh()->n_elements();
408  bdr_shift = reader_mesh->n_elements();
409  }
410 
411  BaseMeshReader::HeaderQuery header_query(field_name_, read_time, this->discretization_, dh_->hash());
412  auto reader = ReaderCache::get_reader(reader_file_);
413  auto header = reader->find_header(header_query);
414  this->input_data_cache_ = reader->template get_element_data<double>(
415  header, n_entities, n_components, bdr_shift);
416 
417  if (is_native) {
418  this->calculate_element_values();
419  } else if (this->interpolation_==DataInterpolation::identic_msh) {
420  this->calculate_element_values();
421  } else if (this->interpolation_==DataInterpolation::equivalent_msh) {
422  this->calculate_element_values();
423  } else if (this->interpolation_==DataInterpolation::gauss_p0) {
424  this->interpolate_gauss();
425  } else { // DataInterpolation::interp_p0
426  this->interpolate_intersection();
427  }
428 
429  return true;
430  } else return false;
431 
432 }
433 
434 
435 template <int spacedim, class Value>
437 {
438  static const unsigned int quadrature_order = 4; // parameter of quadrature
439  std::shared_ptr<Mesh> source_mesh = ReaderCache::get_mesh(reader_file_);
440  std::vector<unsigned int> searched_elements; // stored suspect elements in calculating the intersection
441  std::vector<arma::vec::fixed<3>> q_points; // real coordinates of quadrature points
442  std::vector<double> q_weights; // weights of quadrature points
443  unsigned int quadrature_size=0; // size of quadrature point and weight vector
444  std::vector<double> sum_val(dh_->max_elem_dofs()); // sum of value of one quadrature point
445  unsigned int elem_count; // count of intersect (source) elements of one quadrature point
446  std::vector<double> elem_value(dh_->max_elem_dofs()); // computed value of one (target) element
447  bool contains; // sign if source element contains quadrature point
448 
449  {
450  // set size of vectors to maximal count of quadrature points
451  QGauss quad(3, quadrature_order);
452  q_points.resize(quad.size());
453  q_weights.resize(quad.size());
454  }
455 
456  for (auto cell : dh_->own_range()) {
457  auto ele = cell.elm();
458  std::fill(elem_value.begin(), elem_value.end(), 0.0);
459  switch (cell.dim()) {
460  case 0:
461  quadrature_size = 1;
462  q_points[0] = *ele.node(0);
463  q_weights[0] = 1.0;
464  break;
465  case 1:
466  quadrature_size = compute_fe_quadrature<1>(q_points, q_weights, ele, quadrature_order);
467  break;
468  case 2:
469  quadrature_size = compute_fe_quadrature<2>(q_points, q_weights, ele, quadrature_order);
470  break;
471  case 3:
472  quadrature_size = compute_fe_quadrature<3>(q_points, q_weights, ele, quadrature_order);
473  break;
474  }
475  searched_elements.clear();
476  source_mesh->get_bih_tree().find_bounding_box(ele.bounding_box(), searched_elements);
477 
478  auto r_idx = cell.elm().region_idx().idx();
479  std::string reg_name = cell.elm().region().label();
480  for (unsigned int i=0; i<quadrature_size; ++i) {
481  std::fill(sum_val.begin(), sum_val.end(), 0.0);
482  elem_count = 0;
483  for (std::vector<unsigned int>::iterator it = searched_elements.begin(); it!=searched_elements.end(); it++) {
484  ElementAccessor<3> elm = source_mesh->element_accessor(*it);
485  contains=false;
486  switch (elm->dim()) {
487  case 0:
488  contains = arma::norm(*elm.node(0) - q_points[i], 2) < 4*std::numeric_limits<double>::epsilon();
489  break;
490  case 1:
491  contains = MappingP1<1,3>::contains_point(q_points[i], elm);
492  break;
493  case 2:
494  contains = MappingP1<2,3>::contains_point(q_points[i], elm);
495  break;
496  case 3:
497  contains = MappingP1<3,3>::contains_point(q_points[i], elm);
498  break;
499  default:
500  ASSERT_PERMANENT(false).error("Invalid element dimension!");
501  }
502  if ( contains ) {
503  // projection point in element
504  unsigned int index = sum_val.size() * (*it);
505  for (unsigned int j=0; j < sum_val.size(); j++) {
506  sum_val[j] += get_scaled_value(index+j, dh_->mesh()->elem_index( cell.elm_idx() ), reg_name, region_value_err_[r_idx]);
507  }
508  ++elem_count;
509  }
510  }
511 
512  if (elem_count > 0) {
513  for (unsigned int j=0; j < sum_val.size(); j++) {
514  elem_value[j] += (sum_val[j] / elem_count) * q_weights[i];
515  }
516  }
517  }
518 
519  LocDofVec loc_dofs;
520  loc_dofs = cell.get_loc_dof_indices();
521 
522  ASSERT_LE(loc_dofs.n_elem, elem_value.size());
523  for (unsigned int i=0; i < elem_value.size(); i++) {
524  ASSERT_LT( loc_dofs[i], (int)data_vec_.size());
525  data_vec_.set( loc_dofs[i], elem_value[i] );
526  }
527  }
528 }
529 
530 
531 template <int spacedim, class Value>
533 {
534  std::shared_ptr<Mesh> source_mesh = ReaderCache::get_mesh(reader_file_);
535  std::vector<unsigned int> searched_elements; // stored suspect elements in calculating the intersection
536  std::vector<double> value(dh_->max_elem_dofs());
537  double total_measure;
538  double measure = 0;
539 
540  for (auto elm : dh_->mesh()->elements_range()) {
541  if (elm.dim() == 3) {
542  THROW( ExcInvalidElemeDim() << EI_ElemIdx(elm.idx()) );
543  }
544 
545  double epsilon = 4* numeric_limits<double>::epsilon() * elm.measure();
546  auto r_idx = elm.region_idx().idx();
547  std::string reg_name = elm.region().label();
548 
549  // gets suspect elements
550  if (elm.dim() == 0) {
551  searched_elements.clear();
552  source_mesh->get_bih_tree().find_point(*elm.node(0), searched_elements);
553  } else {
554  BoundingBox bb = elm.bounding_box();
555  searched_elements.clear();
556  source_mesh->get_bih_tree().find_bounding_box(bb, searched_elements);
557  }
558 
559  // set zero values of value object
560  std::fill(value.begin(), value.end(), 0.0);
561  total_measure=0.0;
562 
563  START_TIMER("compute_pressure");
564  ADD_CALLS(searched_elements.size());
565 
566 
567  for (std::vector<unsigned int>::iterator it = searched_elements.begin(); it!=searched_elements.end(); it++)
568  {
569  ElementAccessor<3> source_elm = source_mesh->element_accessor(*it);
570  if (source_elm->dim() == 3) {
571  // get intersection (set measure = 0 if intersection doesn't exist)
572  switch (elm.dim()) {
573  case 0: {
574  arma::vec::fixed<3> real_point = *elm.node(0);
575  arma::mat::fixed<3, 4> elm_map = MappingP1<3,3>::element_map(source_elm);
576  arma::vec::fixed<4> unit_point = MappingP1<3,3>::project_real_to_unit(real_point, elm_map);
577 
578  measure = (std::fabs(arma::sum( unit_point )-1) <= 1e-14
579  && arma::min( unit_point ) >= 0)
580  ? 1.0 : 0.0;
581  break;
582  }
583  case 1: {
584  IntersectionAux<1,3> is(elm.idx(), source_elm.idx());
585  ComputeIntersection<1,3> CI(elm, source_elm);
586  CI.init();
587  CI.compute(is);
588 
589  IntersectionLocal<1,3> ilc(is);
590  measure = ilc.compute_measure() * elm.measure();
591  break;
592  }
593  case 2: {
594  IntersectionAux<2,3> is(elm.idx(), source_elm.idx());
595  ComputeIntersection<2,3> CI(elm, source_elm);
596  CI.init();
597  CI.compute(is);
598 
599  IntersectionLocal<2,3> ilc(is);
600  measure = 2 * ilc.compute_measure() * elm.measure();
601  break;
602  }
603  }
604 
605  //adds values to value_ object if intersection exists
606  if (measure > epsilon) {
607  unsigned int index = value.size() * (*it);
608  for (unsigned int i=0; i < value.size(); i++) {
609  value[i] += get_scaled_value(index+i, dh_->mesh()->elem_index( elm.idx() ), reg_name, region_value_err_[r_idx]) * measure;
610  }
611  total_measure += measure;
612  }
613  }
614  }
615 
616  // computes weighted average, store it to data vector
617  if (total_measure > epsilon) {
618  DHCellAccessor cell = dh_->cell_accessor_from_element(elm.idx());
619  LocDofVec loc_dofs = cell.get_loc_dof_indices();
620 
621  ASSERT_LE(loc_dofs.n_elem, value.size());
622  for (unsigned int i=0; i < value.size(); i++) {
623  data_vec_.set(loc_dofs[i], value[i] / total_measure);
624  }
625  } else {
626  WarningOut().fmt("Processed element with idx {} is out of source mesh!\n", elm.idx());
627  }
628  END_TIMER("compute_pressure");
629 
630  }
631 }
632 
633 
634 template <int spacedim, class Value>
636 {
637  // Same algorithm as in output of Node_data. Possibly code reuse.
638  unsigned int data_vec_i, vec_inc;
639  std::vector<unsigned int> count_vector(data_vec_.size(), 0);
641  std::vector<LongIdx> &source_target_vec = (dynamic_cast<BCMesh*>(dh_->mesh()) != nullptr) ? source_target_mesh_elm_map_->boundary : source_target_mesh_elm_map_->bulk;
642 
643  unsigned int shift;
644  if (this->boundary_domain_) {
645  if (this->interpolation_==DataInterpolation::identic_msh) shift = source_target_mesh_elm_map_->bulk.size();
646  else shift = ReaderCache::get_mesh(reader_file_)->n_elements();
647  } else {
648  shift = 0;
649  }
650 
651  ASSERT_GT(region_value_err_.size(), 0)(field_name_).error("Vector of region isNaN flags is not initialized. Did you call set_mesh or set_fe_data?\n");
652  for (auto r : region_value_err_)
653  r = RegionValueErr();
654 
655  for (auto cell : dh_->own_range()) {
656  LocDofVec loc_dofs = cell.get_loc_dof_indices();
657  //DebugOut() << cell.elm_idx() << " < " << source_target_vec.size() << "\n";
658  int source_idx = source_target_vec[cell.elm_idx()];
659 
660  if (source_idx == (int)(Mesh::undef_idx)) {
661  data_vec_i = source_idx;
662  vec_inc = 0;
663  } else {
664  data_vec_i = (source_idx + shift) * dh_->max_elem_dofs();
665  vec_inc = 1;
666  }
667  auto r_idx = cell.elm().region_idx().idx();
668  std::string reg_name = cell.elm().region().label();
669  for (unsigned int i=0; i<loc_dofs.n_elem; ++i) {
670  ASSERT_LT(loc_dofs[i], (LongIdx)data_vec_.size());
671  data_vec_.add( loc_dofs[i], get_scaled_value(data_vec_i, dh_->mesh()->elem_index( cell.elm_idx() ), reg_name, region_value_err_[r_idx]) );
672  ++count_vector[ loc_dofs[i] ];
673  data_vec_i += vec_inc;
674  }
675  }
676 
677  // compute averages of values
678  for (unsigned int i=0; i<data_vec_.size(); ++i) {
679  if (count_vector[i]>0) data_vec_.normalize(i, count_vector[i]);
680  }
681 }
682 
683 
684 //template <int spacedim, class Value>
685 //void FieldFE<spacedim, Value>::calculate_native_values(ElementDataCache<double>::CacheData data_cache)
686 //{
687 // // Same algorithm as in output of Node_data. Possibly code reuse.
688 // unsigned int data_vec_i;
689 // std::vector<unsigned int> count_vector(data_vec_.size(), 0);
690 // data_vec_.zero_entries();
691 // std::vector<LongIdx> &source_target_vec = (dynamic_cast<BCMesh*>(dh_->mesh()) != nullptr) ? source_target_mesh_elm_map_->boundary : source_target_mesh_elm_map_->bulk;
692 //
693 // // iterate through cells, assembly MPIVector
694 // for (auto cell : dh_->own_range()) {
695 // LocDofVec loc_dofs = cell.get_loc_dof_indices();
696 // data_vec_i = source_target_vec[cell.elm_idx()] * loc_dofs.n_elem;
697 // for (unsigned int i=0; i<loc_dofs.n_elem; ++i, ++data_vec_i) {
698 // data_vec_.add( loc_dofs[i], (*data_cache)[ data_vec_i ] );
699 // ++count_vector[ loc_dofs[i] ];
700 // }
701 // }
702 //
703 // // compute averages of values
704 // for (unsigned int i=0; i<data_vec_.size(); ++i) {
705 // if (count_vector[i]>0) data_vec_.normalize(i, count_vector[i]);
706 // }
707 //}
708 //
709 //
710 //template <int spacedim, class Value>
711 //void FieldFE<spacedim, Value>::calculate_equivalent_values(ElementDataCache<double>::CacheData data_cache)
712 //{
713 // // Same algorithm as in output of Node_data. Possibly code reuse.
714 // unsigned int data_vec_i;
715 // std::vector<unsigned int> count_vector(data_vec_.size(), 0);
716 // data_vec_.zero_entries();
717 // std::vector<LongIdx> &source_target_vec = (dynamic_cast<BCMesh*>(dh_->mesh()) != nullptr) ? source_target_mesh_elm_map_->boundary : source_target_mesh_elm_map_->bulk;
718 // unsigned int shift;
719 // if (this->boundary_domain_) {
720 // if (this->interpolation_==DataInterpolation::identic_msh) shift = this->comp_mesh_->n_elements();
721 // else shift = ReaderCache::get_mesh(reader_file_)->n_elements();
722 // } else {
723 // shift = 0;
724 // }
725 //
726 // // iterate through elements, assembly global vector and count number of writes
727 // for (auto cell : dh_->own_range()) {
728 // LocDofVec loc_dofs = cell.get_loc_dof_indices();
729 // //DebugOut() << cell.elm_idx() << " < " << source_target_vec.size() << "\n";
730 // int source_idx = source_target_vec[cell.elm_idx()];
731 // if (source_idx == (int)(Mesh::undef_idx)) { // undefined value in input data mesh
732 // if ( std::isnan(default_value_) )
733 // THROW( ExcUndefElementValue() << EI_Field(field_name_) );
734 // for (unsigned int i=0; i<loc_dofs.n_elem; ++i) {
735 // ASSERT_LT(loc_dofs[i], (LongIdx)data_vec_.size());
736 // data_vec_.add( loc_dofs[i], default_value_ * this->unit_conversion_coefficient_ );
737 // ++count_vector[ loc_dofs[i] ];
738 // }
739 // } else {
740 // data_vec_i = (source_idx + shift) * dh_->max_elem_dofs();
741 // for (unsigned int i=0; i<loc_dofs.n_elem; ++i, ++data_vec_i) {
742 // ASSERT_LT(loc_dofs[i], (LongIdx)data_vec_.size());
743 // data_vec_.add( loc_dofs[i], (*data_cache)[data_vec_i] );
744 // ++count_vector[ loc_dofs[i] ];
745 // }
746 // }
747 // }
748 //
749 // // compute averages of values
750 // for (unsigned int i=0; i<data_vec_.size(); ++i) {
751 // if (count_vector[i]>0) data_vec_.normalize(i, count_vector[i]);
752 // }
753 //}
754 //
755 //
756 template <int spacedim, class Value>
758  //ASSERT_EQ(output_data_cache.n_values() * output_data_cache.n_comp(), dh_->distr()->lsize()).error();
759  unsigned int n_vals = output_data_cache.n_comp() * output_data_cache.n_dofs_per_element();
760  double loc_values[n_vals];
761  unsigned int i;
762 
763  for (auto dh_cell : dh_->own_range()) {
764  LocDofVec loc_dofs = dh_cell.get_loc_dof_indices();
765  for (i=0; i<loc_dofs.n_elem; ++i) loc_values[i] = data_vec_.get( loc_dofs[i] );
766  for ( ; i<n_vals; ++i) loc_values[i] = numeric_limits<double>::signaling_NaN();
767  output_data_cache.store_value( dh_cell.local_idx(), loc_values );
768  }
769 
770  output_data_cache.set_dof_handler_hash( dh_->hash() );
771 }
772 
773 
774 
775 template <int spacedim, class Value>
776 inline unsigned int FieldFE<spacedim, Value>::data_size() const {
777  return data_vec_.size();
778 }
779 
780 
781 
782 template <int spacedim, class Value>
785 }
786 
787 
788 
789 template <int spacedim, class Value>
792 }
793 
794 
795 
796 template <int spacedim, class Value>
797 double FieldFE<spacedim, Value>::get_scaled_value(int i_cache_el, unsigned int elm_idx, const std::string &region_name,
798  RegionValueErr &actual_compute_region_error) {
799  double return_val;
800  if (i_cache_el == (int)(Mesh::undef_idx))
801  return_val = default_value_;
802  else if ( std::isnan((*input_data_cache_)[i_cache_el]) )
803  return_val = default_value_;
804  else
805  return_val = (*input_data_cache_)[i_cache_el];
806 
807  if ( std::isnan(return_val) ) {
808  actual_compute_region_error = RegionValueErr(region_name, elm_idx, return_val);
809  } else
810  return_val *= this->unit_conversion_coefficient_;
811 
812  return return_val;
813 }
814 
815 
816 
817 /*template <int spacedim, class Value>
818 Armor::ArmaMat<typename Value::element_type, Value::NRows_, Value::NCols_> FieldFE<spacedim, Value>::handle_fe_shape(unsigned int dim,
819  unsigned int i_dof, unsigned int i_qp, unsigned int comp_index)
820 {
821  Armor::ArmaMat<typename Value::element_type, Value::NCols_, Value::NRows_> v;
822  for (unsigned int c=0; c<Value::NRows_*Value::NCols_; ++c)
823  v(c/spacedim,c%spacedim) = fe_values_[dim].shape_value_component(i_dof, i_qp, comp_index+c);
824  if (Value::NRows_ == Value::NCols_)
825  return v;
826  else
827  return v.t();
828 }*/
829 
830 
831 
832 template <int spacedim, class Value>
834 {}
835 
836 
837 // Instantiations of FieldFE
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:301
#define ASSERT_GT(a, b)
Definition of comparative assert macro (Greater Than) only for debug mode.
Definition: asserts.hh:317
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:341
#define ASSERT_LE(a, b)
Definition of comparative assert macro (Less or Equal) only for debug mode.
Definition: asserts.hh:309
ArrayMatSet set(uint index)
Definition: armor.hh:838
Class represents boundary part of mesh.
Definition: bc_mesh.hh:37
Bounding box in 3d ambient space.
Definition: bounding_box.hh:53
unsigned int compute(std::vector< IPAux > &IP13s)
Computes intersection points for 1D-3D intersection. Computes lower dimensional CIs abscissa vs tetra...
void init()
Initializes lower dimensional objects. Sets correctly the pointers to Plucker coordinates and product...
Class for 2D-2D intersections.
void compute(IntersectionAux< 2, 3 > &intersection)
Computes intersection points for 2D-3D intersection. Computes lower dimensional CIs: 1) 3x triangle s...
void init()
Initializes lower dimensional objects. Sets correctly the pointers to Plucker coordinates and product...
Cell accessor allow iterate over DOF handler cells.
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
unsigned int dim() const
Return dimension of element appropriate to cell.
NodeAccessor< 3 > node(unsigned int ni) const
Definition: accessors.hh:230
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
Definition: accessors.hh:215
unsigned int dim() const
Definition: accessors.hh:188
Directing class of FieldValueCache.
unsigned int region_chunk_end(unsigned int region_patch_idx) const
Return end position of region chunk in FieldValueCache.
std::shared_ptr< EvalPoints > eval_points() const
Getter of eval_points object.
unsigned int region_chunk_begin(unsigned int region_patch_idx) const
Return begin position of region chunk in FieldValueCache.
unsigned int region_idx_from_chunk_position(unsigned int chunk_pos) const
Return begin position of region chunk specified by position in map.
const EvalPointData & eval_point_data(unsigned int point_idx) const
Return item of eval_point_data_ specified by its position.
unsigned int n_comp() const
void set_dof_handler_hash(std::size_t hash)
unsigned int n_dofs_per_element() const
void store_value(unsigned int idx, const T *value)
unsigned int dim() const
Definition: elements.h:120
Compound finite element on dim dimensional simplex.
Definition: fe_system.hh:102
const std::vector< std::shared_ptr< FiniteElement< dim > > > & fe() const
Definition: fe_system.hh:147
double unit_conversion_coefficient_
Coeficient of conversion of user-defined unit.
bool is_constant_in_space_
Flag detects that field is only dependent on time.
void init_unit_conversion_coefficient(const Input::Record &rec, const struct FieldAlgoBaseInitData &init_data)
Init value of unit_conversion_coefficient_ from input.
Value value_
Last value, prevents passing large values (vectors) by value.
Common abstract parent of all Field<...> classes.
Definition: field_common.hh:77
const std::string & input_name() const
std::pair< double, double > limits() const
FieldFlag::Flags get_flags() const
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
unsigned int n_comp() const
Field< spacedim, Value >::FieldBasePtr create_field(Input::Record rec, const FieldCommon &field) override
Definition: field_fe.cc:144
std::shared_ptr< DOFHandlerMultiDim > conc_dof_handler_
Definition: field_fe.hh:119
unsigned int index_
Definition: field_fe.hh:118
void cache_reinit(const ElementCacheMap &cache_map) override
Definition: field_fe.cc:262
Armor::ArmaMat< typename Value::element_type, Value::NRows_, Value::NCols_ > handle_fe_shape(unsigned int dim, unsigned int i_dof, unsigned int i_qp)
Definition: field_fe.hh:259
double get_scaled_value(int i_cache_el, unsigned int elm_idx, const std::string &region_name, RegionValueErr &actual_compute_region_error)
Definition: field_fe.cc:797
MixedPtr< FiniteElement > fe_
Definition: field_fe.hh:360
VectorMPI data_vec_
Store data of Field.
Definition: field_fe.hh:326
static const Input::Type::Selection & get_interp_selection_input_type()
Definition: field_fe.cc:105
unsigned int data_size() const
Definition: field_fe.cc:776
std::vector< RegionValueErr > region_value_err_
Set holds data of valid / invalid element values on all regions.
Definition: field_fe.hh:363
void local_to_ghost_data_scatter_begin()
Call begin scatter functions (local to ghost) on data vector.
Definition: field_fe.cc:783
void calculate_element_values()
Calculate data of equivalent_mesh interpolation or native data on input over all elements of target m...
Definition: field_fe.cc:635
virtual void init_from_input(const Input::Record &rec, const struct FieldAlgoBaseInitData &init_data)
Definition: field_fe.cc:285
virtual ~FieldFE()
Destructor.
Definition: field_fe.cc:833
DataInterpolation
Definition: field_fe.hh:61
void local_to_ghost_data_scatter_end()
Call end scatter functions (local to ghost) on data vector.
Definition: field_fe.cc:790
static const Input::Type::Record & get_input_type()
Implementation.
Definition: field_fe.cc:60
std::vector< FEValues< spacedim > > fe_values_
List of FEValues objects of dimensions 0,1,2,3 used for value calculation.
Definition: field_fe.hh:353
void interpolate_gauss()
Interpolate data (use Gaussian distribution) over all elements of target mesh.
Definition: field_fe.cc:436
void make_dof_handler(const MeshBase *mesh)
Create DofHandler object.
Definition: field_fe.cc:345
Quadrature init_quad(std::shared_ptr< EvalPoints > eval_points)
Initialize FEValues object of given dimension.
Definition: field_fe.cc:275
OutputTime::DiscreteSpace discretization_
Specify section where to find the field data in input mesh file.
Definition: field_fe.hh:335
FilePath reader_file_
mesh reader file
Definition: field_fe.hh:329
std::array< FEItem, 4 > fe_item_
Holds specific data of field evaluation over all dimensions.
Definition: field_fe.hh:359
bool set_time(const TimeStep &time) override
Definition: field_fe.cc:385
void interpolate_intersection()
Interpolate data (use intersection library) over all elements of target mesh.
Definition: field_fe.cc:532
VectorMPI set_fe_data(std::shared_ptr< DOFHandlerMultiDim > dh, VectorMPI dof_values=VectorMPI::sequential(0), unsigned int block_index=FieldFE< spacedim, Value >::undef_uint)
Definition: field_fe.cc:175
double default_value_
Default value of element if not set in mesh data file.
Definition: field_fe.hh:344
void cache_update(FieldValueCache< typename Value::element_type > &data_cache, ElementCacheMap &cache_map, unsigned int region_patch_idx) override
Definition: field_fe.cc:215
bool boundary_domain_
Is set in set_mesh method. Value true means, that we accept only boundary element accessors in the va...
Definition: field_fe.hh:350
DataInterpolation interpolation_
Specify type of FE data interpolation.
Definition: field_fe.hh:338
Input::Record in_rec_
Accessor to Input::Record.
Definition: field_fe.hh:347
FieldFlag::Flags flags_
Field flags.
Definition: field_fe.hh:341
void native_data_to_cache(ElementDataCache< double > &output_data_cache)
Definition: field_fe.cc:757
std::string field_name_
field name read from input
Definition: field_fe.hh:332
FieldFE(unsigned int n_comp=0)
Definition: field_fe.cc:134
std::shared_ptr< DOFHandlerMultiDim > dh_
DOF handler object.
Definition: field_fe.hh:324
std::shared_ptr< EquivalentMeshMap > source_target_mesh_elm_map_
Maps element indices from computational mesh to the source (data).
Definition: field_fe.hh:356
static const Input::Type::Selection & get_disc_selection_input_type()
Definition: field_fe.cc:93
ElementDataCache< double >::CacheData input_data_cache_
Input ElementDataCache is stored in set_time and used in all evaluation and interpolation methods.
Definition: field_fe.hh:366
void set_mesh(const Mesh *mesh) override
Definition: field_fe.cc:309
static constexpr Mask declare_input
The field can be set from input. The key in input field descriptor is declared. (default on)
Definition: field_flag.hh:35
static constexpr Mask equation_input
The field is data parameter of the owning equation. (default on)
Definition: field_flag.hh:33
std::shared_ptr< FieldBaseType > FieldBasePtr
Definition: field.hh:96
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
string filename() const
Definition: file_path.cc:188
constexpr bool match(Mask mask) const
Definition: flag_array.hh:163
Accessor to the polymorphic input data of a type given by an AbstracRecord object.
Definition: accessors.hh:458
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
Iterator< ValueType > begin() const
unsigned int size() const
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
bool opt_val(const string &key, Ret &value) const
const Ret val(const string &key) const
Iterator< Ret > find(const string &key) const
Accessor to the data with type Type::Tuple.
Definition: accessors.hh:411
Class for declaration of the input of type Bool.
Definition: type_base.hh:452
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
static Default optional()
The factory function to make an empty default value which is optional.
Definition: type_record.hh:124
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
static FileName input()
The factory function for declaring type FileName for input files.
Definition: type_base.cc:524
Record type proxy class.
Definition: type_record.hh:182
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:196
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Record & copy_keys(const Record &other)
Copy keys from other record.
Definition: type_record.cc:216
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:503
Template for classes storing finite set of named values.
Selection & add_value(const int value, const std::string &key, const std::string &description="", TypeBase::attribute_map attributes=TypeBase::attribute_map())
Adds one new value with name given by key to the Selection.
const Selection & close() const
Close the Selection, no more values can be added.
Class for declaration of the input data that are in string format.
Definition: type_base.hh:582
Internal auxiliary class representing intersection object of simplex<dimA> and simplex<dimB>.
Class represents intersection of two elements.
double compute_measure() const override
Computes the relative measure of intersection object.
static bool contains_point(arma::vec point, ElementAccessor< 3 > elm)
Test if element contains given point.
Definition: mapping_p1.cc:96
static BaryPoint project_real_to_unit(const RealPoint &point, const ElementMap &map)
Definition: mapping_p1.cc:69
static ElementMap element_map(ElementAccessor< 3 > elm)
Definition: mapping_p1.cc:48
Base class for Mesh and BCMesh.
Definition: mesh.h:96
const RegionDB & region_db() const
Definition: mesh.h:175
static const unsigned int undef_idx
Definition: mesh.h:105
Definition: mesh.h:362
BCMesh * bc_mesh() const override
Implement MeshBase::bc_mesh(), getter of boundary mesh.
Definition: mesh.h:567
The class for outputting data during time.
Definition: output_time.hh:51
Symmetric Gauss-Legendre quadrature formulae on simplices.
Base class for quadrature rules on simplices in arbitrary dimensions.
Definition: quadrature.hh:48
unsigned int size() const
Returns number of quadrature points.
Definition: quadrature.hh:86
Armor::Array< double >::ArrayMatSet set(uint i)
Definition: quadrature.hh:93
static std::shared_ptr< EquivalentMeshMap > identic_mesh_map(const FilePath &file_path, Mesh *computational_mesh)
Definition: reader_cache.cc:94
static std::shared_ptr< EquivalentMeshMap > eqivalent_mesh_map(const FilePath &file_path, Mesh *computational_mesh)
static std::shared_ptr< Mesh > get_mesh(const FilePath &file_path)
Definition: reader_cache.cc:41
static std::shared_ptr< BaseMeshReader > get_reader(const FilePath &file_path)
Definition: reader_cache.cc:37
unsigned int size() const
Definition: region.cc:254
static const Input::Type::Tuple & get_input_time_type(double lower_bound=-std::numeric_limits< double >::max(), double upper_bound=std::numeric_limits< double >::max())
Representation of one time step..
double read_time(Input::Iterator< Input::Tuple > time_it, double default_time=std::numeric_limits< double >::quiet_NaN()) const
double end() const
double read_coef(Input::Iterator< Input::Record > unit_it) const
static Input::Type::Default get_input_default()
static const Input::Type::Record & get_input_type()
double get(unsigned int pos) const
Return value on given position.
Definition: vector_mpi.hh:105
void normalize(unsigned int pos, double divisor)
Normalize value on given position.
Definition: vector_mpi.hh:118
void set(unsigned int pos, double val)
Set value on given position.
Definition: vector_mpi.hh:111
void local_to_ghost_end()
local_to_ghost_{begin,end} updates the ghost values on neighbouring processors from local values
Definition: vector_mpi.hh:160
void add(unsigned int pos, double val)
Add value to item on given position.
Definition: vector_mpi.hh:125
static VectorMPI sequential(unsigned int size)
Definition: vector_mpi.cc:44
void local_to_ghost_begin()
local_to_ghost_{begin,end} updates the ghost values on neighbouring processors from local values
Definition: vector_mpi.hh:156
void zero_entries()
Definition: vector_mpi.hh:82
unsigned int size() const
Return size of output data.
Definition: vector_mpi.hh:98
Fundamental simplicial intersections.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Class FESystem for compound finite elements.
MixedPtr< FESystem > mixed_fe_system(MixedPtr< FiniteElement > fe, Args &&... args)
Definition: fe_system.hh:176
#define INSTANCE_ALL(field)
@ FETensor
@ FEVector
@ Value
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:104
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
arma::Col< IntIdx > LocDofVec
Definition: index_types.hh:28
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Definition: index_types.hh:24
Internal class representing intersection object.
Classes with algorithms for computation of intersections of meshes.
static constexpr bool value
Definition: json.hpp:87
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
Class MappingP1 implements the affine transformation of the unit cell onto the actual cell.
typename arma::Mat< Type >::template fixed< nr, nc > ArmaMat
Definition: armor.hh:502
DummyInt isnan(...)
Definition: format.h:306
Definitions of particular quadrature rules on simplices.
Implementation of range helper class.
Definition: mixed.hh:25
unsigned int i_eval_point_
index of point in EvalPoint object
unsigned int i_element_
mesh_idx of ElementAccessor appropriate to element
unsigned int i_reg_
region_idx of element
Helper struct stores data for initizalize descentants of FieldAlgorithmBase.
FieldFlag::Flags flags_
#define END_TIMER(tag)
Ends a timer with specified tag.
#define ADD_CALLS(n_calls)
Increase number of calls in actual timer.
#define START_TIMER(tag)
Starts a timer with specified tag.
@ update_values
Shape function values.
Definition: update_flags.hh:87