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