Flow123d  master-f44eb46
field_formula.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_formula.cc
15  * @brief
16  */
17 
18 
19 
20 #include "fields/field_formula.hh"
21 #include "fields/field_instances.hh" // for instantiation macros
22 #include "fields/surface_depth.hh"
23 #include "input/input_type.hh"
24 //#include "include/arena_alloc.hh" // bparser
25 #include <boost/algorithm/string/replace.hpp>
26 #include <boost/regex.hpp>
27 
28 
29 
30 /// Implementation.
31 
32 namespace it = Input::Type;
33 
34 FLOW123D_FORCE_LINK_IN_CHILD(field_formula)
35 
36 
37 template <int spacedim, class Value>
39 {
40 
41  return it::Record("FieldFormula", FieldAlgorithmBase<spacedim,Value>::template_name()+" Field given by runtime interpreted formula.")
45  "String, array of strings, or matrix of strings with formulas for individual "
46  "entries of scalar, vector, or tensor value respectively.\n"
47  "For vector values, you can use just one string to enter homogeneous vector.\n"
48  "For square (($N\\times N$))-matrix values, you can use:\n\n"
49  " - array of strings of size (($N$)) to enter diagonal matrix\n"
50  " - array of strings of size (($\\frac12N(N+1)$)) to enter symmetric matrix (upper triangle, row by row)\n"
51  " - just one string to enter (spatially variable) multiple of the unit matrix.\n"
52  "Formula can contain variables ```x,y,z,t,d``` and usual operators and functions." )
53  //.declare_key("unit", FieldAlgorithmBase<spacedim, Value>::get_input_type_unit_si(), it::Default::optional(),
54  // "Definition of unit.")
55  .declare_key("surface_direction", it::String(), it::Default("\"0 0 1\""),
56  "The vector used to project evaluation point onto the surface.")
57  .declare_key("surface_region", it::String(), it::Default::optional(),
58  "The name of region set considered as the surface. You have to set surface region if you "
59  "want to use formula variable ```d```.")
60  .allow_auto_conversion("value")
61  .close();
62 }
63 
64 
65 
66 template <int spacedim, class Value>
68  Input::register_class< FieldFormula<spacedim, Value>, unsigned int >("FieldFormula") +
70 
71 
72 
73 template <int spacedim, class Value>
75 : FieldAlgorithmBase<spacedim, Value>(n_comp),
76  b_parser_( CacheMapElementNumber::get() ),
77  arena_alloc_(nullptr)
78 {
79  this->is_constant_in_space_ = false;
80 }
81 
82 
83 
84 template <int spacedim, class Value>
86  this->init_unit_conversion_coefficient(rec, init_data);
87 
88  // read formulas form input
89  this->formula_ = rec.val<std::string>("value");
90  in_rec_ = rec;
91 }
92 
93 
94 template <int spacedim, class Value>
96 
97  this->time_=time;
98  this->is_constant_in_space_ = false;
99  return true;
100 
101 }
102 
103 
104 template <int spacedim, class Value>
106  // create SurfaceDepth object if surface region is set
107  std::string surface_region;
108  if ( in_rec_.opt_val("surface_region", surface_region) ) {
109  surface_depth_ = std::make_shared<SurfaceDepth>(mesh, surface_region, in_rec_.val<std::string>("surface_direction"));
110  }
111 }
112 
113 
114 template <int spacedim, class Value>
116  ElementCacheMap &cache_map, unsigned int region_patch_idx)
117 {
118  unsigned int reg_chunk_begin = cache_map.region_chunk_begin(region_patch_idx);
119  unsigned int reg_chunk_end = cache_map.region_chunk_end(region_patch_idx);
120 
121  for (unsigned int i=reg_chunk_begin; i<reg_chunk_end; ++i) {
122  res_[i] = 0.0;
123  }
124  for (auto it : eval_field_data_) {
125  // Copy data from dependent fields to arena. Temporary solution.
126  // TODO hold field data caches in arena, remove this step
127  auto value_cache = it.first->value_cache();
128  for (unsigned int i=reg_chunk_begin; i<reg_chunk_end; ++i) {
129  if (it.first->name() == "X") {
130  x_[i] = value_cache->template vec<3>(i)(0);
131  y_[i] = value_cache->template vec<3>(i)(1);
132  z_[i] = value_cache->template vec<3>(i)(2);
133  } else
134  it.second[i] = value_cache->data_[i];
135  }
136  }
137 
138  // Get vector of subsets as subarray
139  uint subsets_begin = reg_chunk_begin / cache_map.simd_size_double;
140  uint subsets_end = reg_chunk_end / cache_map.simd_size_double;
141  std::vector<uint> subset_vec;
142  subset_vec.assign(subsets_ + subsets_begin, subsets_ + subsets_end);
143 
144  b_parser_.set_subset(subset_vec);
145  b_parser_.run();
146  uint vec_size = CacheMapElementNumber::get();
147  for(unsigned int row=0; row < this->value_.n_rows(); row++)
148  for(unsigned int col=0; col < this->value_.n_cols(); col++) {
149  uint comp_shift = (row*this->value_.n_cols()+col) * vec_size;
150  for (unsigned int i=reg_chunk_begin; i<reg_chunk_end; ++i) {
151  auto cache_val = data_cache.template mat<Value::NRows_, Value::NCols_>(i);
152  cache_val(row, col) = this->unit_conversion_coefficient_ * res_[i + comp_shift];
153  data_cache.set(i) = cache_val;
154  }
155  }
156 }
157 
158 
159 template <int spacedim, class Value>
161 {
162  if (surface_depth_ && has_depth_var_) {
163  // add value of depth
164  arma::vec p_depth(spacedim+1);
165  p_depth.subvec(0,spacedim-1) = p;
166  try {
167  p_depth(spacedim) = surface_depth_->compute_distance(p);
168  } catch (SurfaceDepth::ExcTooLargeSnapDistance &e) {
169  e << SurfaceDepth::EI_FieldTime(this->time_.end());
170  e << in_rec_.ei_address();
171  throw;
172  }
173  return p_depth;
174  } else {
175  return p;
176  }
177 }
178 
179 template <int spacedim, class Value>
181  required_fields_.clear(); // returned value
182 
183  // set expression and data to BParser
184  try {
185  b_parser_.parse( formula_ );
186  } catch (std::exception const& e) {
187  if (typeid(e) == typeid(bparser::Exception))
188  THROW( ExcParserError() << EI_BParserMsg(e.what()) << EI_Formula(formula_) << Input::EI_Address( in_rec_.address_string() ) );
189  else throw;
190  }
191  std::vector<std::string> variables = b_parser_.free_symbols();
192  std::sort( variables.begin(), variables.end() );
193  variables.erase( std::unique( variables.begin(), variables.end() ), variables.end() );
194 
195  has_time_=false;
196  sum_shape_sizes_=0; // scecifies size of arena
197  for (auto var : variables) {
198  if (var == "X" || var == "x" || var == "y" || var == "z") {
199  required_fields_.push_back( field_set.field("X") );
200  sum_shape_sizes_ += spacedim;
201  }
202  else if (var == "t") has_time_ = true;
203  else {
204  auto field_ptr = field_set.field(var);
205  if (field_ptr != nullptr) required_fields_.push_back( field_ptr );
206  else THROW( FieldSet::ExcUnknownField() << FieldCommon::EI_Field(var) << FieldSet::EI_FieldType("formula") << Input::EI_Address( in_rec_.address_string() ) );
207  // TODO: Test the exception, report input line of the formula.
208  if (field_ptr->value_cache() == nullptr) THROW( ExcNotDoubleField() << EI_Field(var) << Input::EI_Address( in_rec_.address_string() ) );
209  // TODO: Test the exception, report input line of the formula.
210 
211  sum_shape_sizes_ += field_ptr->n_shape();
212  if (var == "d") {
213  field_set.set_surface_depth(this->surface_depth_);
214  }
215  }
216  }
217 
218  return required_fields_;
219 }
220 
221 
222 template <int spacedim, class Value>
224 {
225  // Can not compile expression in set_time as the necessary cache size is not known there yet.
226 
227  if (arena_alloc_!=nullptr) {
228  delete arena_alloc_;
229  }
230  eval_field_data_.clear();
231  uint vec_size = CacheMapElementNumber::get();
232 
233  // number of subset alignment to block size
234  uint n_subsets = vec_size / cache_map.simd_size_double;
235  uint res_comp = Value::NRows_ * Value::NCols_;
236  uint n_vectors = sum_shape_sizes_ + res_comp; // needs add space of result vector
237  arena_alloc_ = new bparser::ArenaAlloc(cache_map.simd_size_double, n_vectors * vec_size * sizeof(double) + n_subsets * sizeof(uint));
238  res_ = arena_alloc_->create_array<double>(vec_size * res_comp);
239  for (auto field : required_fields_) {
240  std::string field_name = field->name();
241  eval_field_data_[field] = arena_alloc_->create_array<double>(field->n_shape() * vec_size);
242  if (field_name == "X") {
243  X_ = eval_field_data_[field] + 0;
244  x_ = eval_field_data_[field] + 0;
245  y_ = eval_field_data_[field] + vec_size;
246  z_ = eval_field_data_[field] + 2*vec_size;
247  }
248  }
249  subsets_ = arena_alloc_->create_array<uint>(n_subsets);
250 
251  // set expression and data to BParser
252  if (has_time_) {
253  b_parser_.set_constant("t", {}, {this->time_.end()});
254  }
255  for (auto field : required_fields_) {
256  std::string field_name = field->name();
257  if (field_name == "X") {
258  b_parser_.set_variable("X", {3}, X_);
259  b_parser_.set_variable("x", {}, x_);
260  b_parser_.set_variable("y", {}, y_);
261  b_parser_.set_variable("z", {}, z_);
262  } else {
263  std::vector<uint> f_shape = {};
264  if (field->n_shape() > 1) f_shape = field->shape_;
265  b_parser_.set_variable(field_name, f_shape, eval_field_data_[field]);
266  }
267  }
268  std::vector<uint> shape = {};
269  if (Value::NRows_ > 1) shape.push_back(Value::NRows_);
270  if (Value::NCols_ > 1) shape.push_back(Value::NCols_);
271  b_parser_.set_variable("_result_", shape, res_);
272  b_parser_.compile();
273  // set subset vector
274  for (uint i=0; i<n_subsets; ++i)
275  subsets_[i] = i;
276 }
277 
278 
279 template <int spacedim, class Value>
281 
282 
283 // Instantiations of FieldFormula
ArrayMatSet set(uint index)
Definition: armor.hh:838
Auxiliary data class holds number of elements in cache and allow to set this value explicitly (e....
static unsigned int get()
Return number of stored elements.
Directing class of FieldValueCache.
unsigned int region_chunk_end(unsigned int region_patch_idx) const
Return end position of region chunk in FieldValueCache.
unsigned int region_chunk_begin(unsigned int region_patch_idx) const
Return begin position of region chunk in FieldValueCache.
unsigned int simd_size_double
bool is_constant_in_space_
Flag detects that field is only dependent on time.
Space< spacedim >::Point Point
void cache_reinit(const ElementCacheMap &cache_map) override
std::vector< const FieldCommon * > set_dependency(FieldSet &field_set) override
static const Input::Type::Record & get_input_type()
Implementation.
FieldFormula(unsigned int n_comp=0)
arma::vec eval_depth_var(const Point &p)
virtual ~FieldFormula()
void cache_update(FieldValueCache< typename Value::element_type > &data_cache, ElementCacheMap &cache_map, unsigned int region_patch_idx) override
virtual void init_from_input(const Input::Record &rec, const struct FieldAlgoBaseInitData &init_data)
void set_mesh(const Mesh *mesh) override
bool set_time(const TimeStep &time) override
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
void set_surface_depth(std::shared_ptr< SurfaceDepth > surface_depth)
Set surface depth object to "d" field.
Definition: field_set.hh:385
FieldCommon * field(const std::string &field_name) const
Definition: field_set.cc:173
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
const Ret val(const string &key) const
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
Record type proxy class.
Definition: type_record.hh:182
virtual Record & allow_auto_conversion(const string &from_key)
Allows shorter input of the Record providing only value of the from_key given as the parameter.
Definition: type_record.cc:133
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
Class for declaration of the input data that are in string format.
Definition: type_base.hh:582
Definition: mesh.h:362
Representation of one time step..
#define INSTANCE_ALL(field)
@ 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
unsigned int uint
ArmaVec< double, N > vec
Definition: armor.hh:885
#define FMT_UNUSED
Definition: posix.h:75
Helper struct stores data for initizalize descentants of FieldAlgorithmBase.