Flow123d  master-f44eb46
field_values.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file field_values.hh
15  * @brief
16  */
17 
18 #ifndef FIELD_VALUES_HH_
19 #define FIELD_VALUES_HH_
20 
21 #include <string.h> // for memcpy
22 #include <boost/core/explicit_operator_bool.hpp> // for optional::oper...
23 #include <boost/optional/optional.hpp> // for get_pointer
24 #include <cmath> // for abs
25 #include <cstdlib> // for abs
26 #include <limits> // for numeric_limits
27 #include <new> // for operator new[]
28 #include <ostream> // for operator&, ope...
29 #include <string> // for basic_string
30 #include <system/exceptions.hh> // for THROW, ExcStream
31 #include <type_traits> // for is_floating_point
32 #include <vector> // for vector
33 #include <armadillo>
34 
35 #include "system/fmt/format.h"
36 #include "input/accessors.hh" // for Array, Iterator
37 #include "input/accessors_impl.hh" // for Array::size
38 #include "input/input_exception.hh" // for ExcFV_Input::~...
39 #include "input/type_base.hh" // for Array, String
40 #include "input/type_generic.hh" // for Parameter
41 #include "input/type_selection.hh" // for Selection
42 #include "system/armor.hh" // for Armor::Array
43 
44 namespace IT=Input::Type;
45 
46 /**
47  * @file
48  *
49  * This file contains various dispatch classes to simplify implementation of Fields. Essential is class template
50  * @p FieldValues_ which provides unified access and initialization to scalar, vector and matrix type object in Armadillo library
51  */
52 
53 TYPEDEF_ERR_INFO( EI_InputMsg, const string );
54 DECLARE_INPUT_EXCEPTION( ExcFV_Input, << "Wrong field value input: " << EI_InputMsg::val );
55 
56 typedef unsigned int FieldEnum;
57 
58 
59 
60 namespace internal {
61 
62 // Helper functions to get scalar type name
63 /*
64 std::string type_name_(double);
65 std::string type_name_(int);
66 std::string type_name_(FieldEnum);
67 */
68 
69 /**
70  * InputType dispatch from elementary type @p ET of FieldValues_ to elementary Input::Type, i.e. descendant of Input::Type::Scalar.
71  */
72 template <class ET>
73 struct InputType { typedef Input::Type::Double type; };
74 
75 template <>
76 struct InputType<int> { typedef Input::Type::Integer type; };
77 
78 template <>
80 
81 
82 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
83 // resolution of Value::return_type
84 
85 // general element type
86 template<int NRows, int NCols, class ET>
87 struct ReturnType { typedef typename arma::Mat<ET>::template fixed<NRows, NCols> return_type; };
88 
89 template<class ET>
90 struct ReturnType<1,1,ET> { typedef ET return_type; };
91 
92 template <class ET>
93 struct ReturnType<0,1,ET> { typedef arma::Col<ET> return_type; };
94 
95 template <int NRows, class ET>
96 struct ReturnType<NRows,1,ET> { typedef typename arma::Col<ET>::template fixed<NRows> return_type; };
97 
98 
99 // FiledEnum element type - this just returns types with ET=unsigned int, however input should be different
100 template<int NRows, int NCols>
101 struct ReturnType<NRows, NCols, FieldEnum> { typedef typename arma::Mat<unsigned int>::template fixed<NRows, NCols> return_type; };
102 
103 template<>
104 struct ReturnType<1,1, FieldEnum> { typedef unsigned int return_type; };
105 
106 template <>
107 struct ReturnType<0,1, FieldEnum> { typedef arma::Col<unsigned int> return_type; };
108 
109 template <int NRows>
110 struct ReturnType<NRows,1, FieldEnum> { typedef typename arma::Col<unsigned int>::template fixed<NRows> return_type; };
111 
112 
113 // Resolution of helper functions for raw constructor
114 template <class RT> inline RT & set_raw_scalar(RT &, double *raw_data) { return *raw_data;}
115 template <class RT> inline RT & set_raw_scalar(RT &, int *raw_data) { return *raw_data;}
116 template <class RT> inline RT & set_raw_scalar(RT &, FieldEnum *raw_data) { return *raw_data;}
117 
118 template <class RT> inline RT & set_raw_vec(RT &val, double *raw_data) { arma::access::rw(val.mem) = raw_data; return val;}
119 template <class RT> inline RT & set_raw_vec(RT &val, int *raw_data) { arma::access::rw(val.mem) = raw_data; return val;}
120 template <class RT> inline RT & set_raw_vec(RT &val, FieldEnum *raw_data) { arma::access::rw(val.mem) = raw_data; return val;}
121 
122 template <class RT> inline RT & set_raw_fix(RT &val, double *raw_data) { val = RT(raw_data); return val;}
123 template <class RT> inline RT & set_raw_fix(RT &val, int *raw_data) { val = RT(raw_data); return val;}
124 template <class RT> inline RT & set_raw_fix(RT &val, FieldEnum *raw_data) { val = RT(raw_data); return val;}
125 
126 
127 
128 // Resolution of input accessor for vector values of enums,
129 template <class ET>
130 struct AccessTypeDispatch { typedef ET type;};
131 template <>
132 struct AccessTypeDispatch<unsigned int> { typedef Input::Enum type; };
133 
134 
135 
136 
137 /**
138  * Initialize an Armadillo matrix from the input.
139  * Since only limited methods are used we can use this template to initialize StringTensor as well.
140  * Used methodas are: at(row,col), zeros()
141  */
142 template<class MatrixType>
143 void init_matrix_from_input( MatrixType &value, Input::Array rec ) {
144  typedef typename MatrixType::elem_type ET;
145  unsigned int nrows = value.n_rows;
146  unsigned int ncols = value.n_cols;
147 
149  if (it->size() == 1 && nrows == ncols) {
150  // square tensor
151  // input = 3 expands to [ [ 3 ] ]; init to 3 * (identity matrix)
152  // input = [1, 2, 3] expands to [[1], [2], [3]]; init to diag. matrix
153  // input = [1, 2, 3, .. , (N+1)*N/2], .... ; init to symmetric matrix [ [1 ,2 ,3], [2, 4, 5], [ 3, 5, 6] ]
154  if (rec.size() == 1) {// scalar times identity
155  value.zeros();
156  ET scalar=*(it->begin<ET>());
157  for(unsigned int i=0; i< nrows; i++) value.at(i,i)=scalar;
158  } else if (rec.size() == nrows) { // diagonal vector
159  value.zeros();
160  for(unsigned int i=0; i< nrows; i++, ++it) value.at(i,i)=*(it->begin<ET>());
161  } else if (rec.size() == (nrows+1)*nrows/2) { // symmetric part
162  for( unsigned int row=0; row<nrows; row++)
163  for( unsigned int col=0; col<ncols; col++)
164  if (row <= col) {
165  value.at(row,col) = *(it->begin<ET>());
166  ++it;
167  } else value.at(row,col) = value.at(col,row);
168  } else {
169  THROW( ExcFV_Input()
170  << EI_InputMsg(
171  fmt::format("Initializing symmetric matrix {:d}x{:d} by vector of wrong size {:d}, "
172  "should be 1, {:d}, or {:d}.",
173  nrows, ncols, rec.size(), nrows, (nrows+1)*nrows/2))
174  << rec.ei_address()
175 
176  );
177  }
178 
179  } else {
180  // accept only full tensor
181  if (rec.size() == nrows && it->size() == ncols) {
182 
183  for (unsigned int row = 0; row < nrows; row++, ++it) {
184  if (it->size() != ncols)
185  THROW( ExcFV_Input() << EI_InputMsg("Wrong number of columns.")
186  << rec.ei_address());
187  Input::Iterator<ET> col_it = it->begin<ET>();
188  for (unsigned int col = 0; col < ncols; col++, ++col_it)
189  value.at(row, col) = *col_it;
190  }
191  } else {
192  THROW( ExcFV_Input()
193  << EI_InputMsg(
194  fmt::format("Initializing symmetric matrix {:d}x{:d} by vector of wrong size {:d}x{:d}.",
195  nrows, ncols, rec.size(), it->size()))
196  << rec.ei_address()
197  );
198  }
199  }
200 }
201 
202 /**
203  * Initialize an Armadillo vector from the input.
204  * Since only limited methods are used we can use this template to initialize StringTensor as well.
205  * Used methodas are: at(row,col), zeros()
206  */
207 template<class VectorType>
208 void init_vector_from_input( VectorType &value, Input::Array rec ) {
209  typedef typename VectorType::elem_type ET;
210  unsigned int nrows = value.n_rows;
211 
212  typedef typename AccessTypeDispatch<ET>::type InnerType;
213  auto it = rec.begin<Input::Array>();
214 
215  if ( it->size() == 1 && rec.size() == 1 ) {
216  InnerType v = *(it->begin<InnerType>());
217  for(unsigned int i=0; i< nrows; i++)
218  value.at(i) = v;
219  } else if ( it->size() == 1 && rec.size() == nrows ) {
220  for(unsigned int i=0; i< nrows; i++, ++it) {
221  InnerType v = *(it->begin<InnerType>());
222  value.at(i) = v;
223  }
224  } else {
225  THROW( ExcFV_Input()
226  << EI_InputMsg(fmt::format("Initializing vector of size {:d} by vector of other size or by tensor", nrows))
227  << rec.ei_address()
228  );
229  }
230 }
231 
232 
233 
234 } // namespace internal
235 
236 
237 
238 /**
239  * Template for class representing all possible Filed values. It is just common interface to
240  * scalar (double, int) values, vector and tensor values with fixed size and vector/tensor values with variable size.
241  *
242  * ET is type of elements, n_cols and n_rows gives fixed dimensions of the tensor value (nx1 is vector, 1x1 is scalar,
243  * 0x1 is variable size vector, 0x0 is variable size tensor (not implemented yet) )
244  *
245  *
246  */
247 template <int NRows, int NCols, class ET>
248 class FieldValue_ {
249 public:
250  typedef ET element_type;
254  const static int NRows_ = NRows;
255  const static int NCols_ = NCols;
256  const static int rank_ = 2;
257 
258  static std::string type_name() { return fmt::format("R[{:d},{:d}]", NRows, NCols); }
259  static constexpr bool is_scalable() {
261  }
262 
263 
264  inline FieldValue_(return_type &val) : value_(val) {}
265  inline static const return_type &from_raw(return_type &val, ET *raw_data) {return internal::set_raw_fix(val, raw_data);}
266  const ET * mem_ptr() const {
267  return value_.memptr();
268  }
269 
270  /// Casts value stored in Armor::Array to return type.
272  return arr.template mat<NRows, NCols>(idx);
273  }
274 
277  }
278 
279  void set_n_comp(unsigned int) {};
280  inline unsigned int n_cols() const
281  { return NCols; }
282  inline unsigned int n_rows() const
283  { return NRows; }
284  inline ET &operator() ( unsigned int i, unsigned int j)
285  { return value_.at(i,j); }
286  inline ET operator() ( unsigned int i, unsigned int j) const
287  { return value_.at(i,j); }
288  inline operator return_type() const
289  { return value_;}
290 
291  // Set value to matrix of zeros.
292  void zeros() {
293  value_.zeros();
294  }
295  // Set value to identity matrix.
296  void eye() {
297  value_.eye();
298  }
299  // Set value to matrix of ones.
300  void ones() {
301  value_.ones();
302  }
303  // Elementwise equivalence.
304  bool equal_to(const return_type &other) {
305  return arma::max(arma::max(arma::abs(value_ - other))) < 4*std::numeric_limits<ET>::epsilon();
306  }
307  // Multiplied value_ by double coefficient
308  void scale(double scale_coef) {
309  if (is_scalable())
310  for( unsigned int row=0; row<n_rows(); row++)
311  for( unsigned int col=0; col<n_cols(); col++)
312  value_.at(row,col) = scale_coef * value_.at(row,col);
313  }
314 
315 private:
317 };
318 
319 
320 
321 
322 
323 /// **********************************************************************
324 /// Specialization for scalars
325 template <class ET>
326 class FieldValue_<1,1,ET> {
327 public:
328  typedef ET element_type;
332  const static int NRows_ = 1;
333  const static int NCols_ = 1;
334  const static int rank_ = 0;
335 
336  static std::string type_name() { return "R"; }
337  static constexpr bool is_scalable() {
339  }
340 
341  inline FieldValue_(return_type &val) : value_(val) {}
342 
343  /**
344  * Returns reference to the return_type (i.e. double, or arma::vec or arma::mat); with data provided by the parameter @p raw_data.
345  * A reference to a work space @p val has to be provided for efficient work with vector and matrix values.
346  */
347  inline static const return_type &from_raw(return_type &val, ET *raw_data) {return internal::set_raw_scalar(val, raw_data);}
348  const ET * mem_ptr() const { return &(value_); }
349 
350  /// Casts value stored in Armor::Array to return type.
352  return arr.scalar(idx);
353  }
354 
356  auto it = rec.begin<Input::Array>();
357  if (it->size() == 1 && rec.size() == 1) {
358  AccessType scalar = *(it->begin<AccessType>());
359  value_ = return_type(scalar);
360  } else {
361  THROW( ExcFV_Input()
362  << EI_InputMsg("Initializing scalar field value by vector or tensor")
363  << rec.ei_address()
364  );
365 
366  }
367  }
368 
369  void set_n_comp(unsigned int) {};
370  inline unsigned int n_cols() const
371  { return 1; }
372  inline unsigned int n_rows() const
373  { return 1; }
374  inline ET &operator() ( unsigned int, unsigned int )
375  { return value_; }
376  inline ET operator() ( unsigned int, unsigned int) const
377  { return value_; }
378  inline operator return_type() const
379  { return value_;}
380 
381  // Set value to matrix of zeros.
382  void zeros() {
383  value_=0;
384  }
385  // Set value to identity matrix.
386  void eye() {
387  value_=1;
388  }
389  // Set value to matrix of ones.
390  void ones() {
391  value_=1;
392  }
393  bool equal_to(const return_type &other) {
394  return std::abs((double)value_ - other) < 4*std::numeric_limits<ET>::epsilon();
395  }
396  // Multiplied value_ by double coefficient
397  void scale(double scale_coef) {
398  if (is_scalable())
399  value_ = scale_coef * value_;
400  }
401 
402 private:
404 };
405 
406 
407 
408 
409 /// **********************************************************************
410 /// Specialization for variable size vectors
411 template <class ET>
412 class FieldValue_<0,1,ET> {
413 public:
414  typedef ET element_type;
418  const static int NRows_ = 0;
419  const static int NCols_ = 1;
420  const static int rank_ = 10;
421 
422 
423  static std::string type_name() { return "R[n]"; }
424  static constexpr bool is_scalable() {
426  }
427  inline static const return_type &from_raw(return_type &val, ET *raw_data) {return internal::set_raw_vec(val, raw_data);}
428  const ET * mem_ptr() const { return value_.memptr(); }
429 
430  /// Casts value stored in Armor::Array to return type.
432  return arr.template vec<NRows_>(idx);
433  }
434 
435  inline FieldValue_(return_type &val) : value_(val) {}
436 
437 
440  }
441 
442  void set_n_comp(unsigned int n_comp) { value_ = return_type(n_comp,1); };
443  inline unsigned int n_cols() const
444  { return 1; }
445  inline unsigned int n_rows() const
446  { return value_.n_rows; }
447  inline ET &operator() ( unsigned int i, unsigned int )
448  { return value_.at(i); }
449  inline ET operator() ( unsigned int i, unsigned int ) const
450  { return value_.at(i); }
451 
452  inline operator return_type() const
453  { return value_;}
454 
455  // Set value to matrix of zeros.
456  void zeros() {
457  value_.zeros();
458  }
459  // Set value to identity matrix.
460  void eye() {
461  value_.ones();
462  }
463  // Set value to matrix of ones.
464  void ones() {
465  value_.ones();
466  }
467  // Elementwise equivalence.
468  bool equal_to(const return_type &other) {
469  return arma::max(arma::abs(value_ - other)) < 4*std::numeric_limits<ET>::epsilon();
470  }
471  // Multiplied value_ by double coefficient
472  void scale(double scale_coef) {
473  if (is_scalable())
474  for(unsigned int i=0; i< n_rows(); i++) {
475  value_.at(i) = scale_coef * value_.at(i);
476  }
477  }
478 
479 
480 private:
482 };
483 
484 /// **********************************************************************
485 /// Specialization for fixed size vectors
486 template <int NRows, class ET>
487 class FieldValue_<NRows,1,ET> {
488 public:
489  typedef ET element_type;
493  const static int NRows_ = NRows;
494  const static int NCols_ = 1;
495  const static int rank_ = 1;
496 
497 
498  static std::string type_name() { return fmt::format("R[{:d}]", NRows); }
499  static constexpr bool is_scalable() {
501  }
502 
503  inline FieldValue_(return_type &val) : value_(val) {}
504  inline static const return_type &from_raw(return_type &val, ET *raw_data) {return internal::set_raw_fix(val, raw_data);}
505  const ET * mem_ptr() const { return value_.memptr(); }
506 
507  /// Casts value stored in Armor::Array to return type.
509  return arr.template vec<NRows>(idx);
510  }
511 
514  }
515 
516  void set_n_comp(unsigned int) {};
517  inline unsigned int n_cols() const
518  { return 1; }
519  inline unsigned int n_rows() const
520  { return NRows; }
521  inline ET &operator() ( unsigned int i, unsigned int )
522  { return value_.at(i); }
523  inline ET operator() ( unsigned int i, unsigned int ) const
524  { return value_.at(i); }
525 
526  inline operator return_type() const
527  { return value_;}
528 
529  // Set value to matrix of zeros.
530  void zeros() {
531  value_.zeros();
532  }
533  // Set value to identity matrix.
534  void eye() {
535  value_.ones();
536  }
537  // Set value to matrix of ones.
538  void ones() {
539  value_.ones();
540  }
541  // Elementwise equivalence.
542  bool equal_to(const return_type &other) {
543  return arma::max(arma::abs(value_ - other)) < 4*std::numeric_limits<ET>::epsilon();
544  }
545  // Multiplied value_ by double coefficient
546  void scale(double scale_coef) {
547  if (is_scalable())
548  for(unsigned int i=0; i< n_rows(); i++) {
549  value_.at(i) = scale_coef * value_.at(i);
550  }
551  }
552 
553 private:
555 };
556 
557 
558 
559 
560 /**
561  * Class that provides common template-less interface to all templated Field structes.
562  *
563  * Maybe we can delete this since common interface will be given by "Quantity".
564  */
565 
566 template <int spacedim>
567 struct FieldValue {
568  // typedefs for possible field values
577 };
578 
579 
580 
581 
582 
583 
584 
585 
586 
587 
588 #endif /* FIELD_VALUES_HH_ */
Type scalar(uint mat_index) const
Definition: armor.hh:831
unsigned int n_cols() const
void scale(double scale_coef)
internal::ReturnType< 0, 1, ET >::return_type return_type
static const return_type & from_raw(return_type &val, ET *raw_data)
static constexpr bool is_scalable()
void set_n_comp(unsigned int n_comp)
unsigned int n_rows() const
const ET * mem_ptr() const
static std::string type_name()
FieldValue_(return_type &val)
bool equal_to(const return_type &other)
internal::InputType< ET >::type ElementInputType
static return_type get_from_array(const Armor::Array< element_type > &arr, uint idx)
Casts value stored in Armor::Array to return type.
void init_from_input(AccessType rec)
unsigned int n_cols() const
static return_type get_from_array(const Armor::Array< element_type > &arr, uint idx)
Casts value stored in Armor::Array to return type.
static constexpr bool is_scalable()
unsigned int n_rows() const
const ET * mem_ptr() const
FieldValue_(return_type &val)
bool equal_to(const return_type &other)
static std::string type_name()
void set_n_comp(unsigned int)
void init_from_input(Input::Array rec)
internal::ReturnType< 1, 1, ET >::return_type return_type
internal::InputType< ET >::type ElementInputType
void scale(double scale_coef)
internal::AccessTypeDispatch< ET >::type AccessType
static const return_type & from_raw(return_type &val, ET *raw_data)
unsigned int n_cols() const
unsigned int n_rows() const
internal::ReturnType< NRows, 1, ET >::return_type return_type
static std::string type_name()
void set_n_comp(unsigned int)
internal::InputType< ET >::type ElementInputType
bool equal_to(const return_type &other)
static return_type get_from_array(const Armor::Array< element_type > &arr, uint idx)
Casts value stored in Armor::Array to return type.
static constexpr bool is_scalable()
static const return_type & from_raw(return_type &val, ET *raw_data)
void scale(double scale_coef)
FieldValue_(return_type &val)
void init_from_input(AccessType rec)
const ET * mem_ptr() const
void init_from_input(AccessType rec)
Input::Array AccessType
internal::InputType< ET >::type ElementInputType
static return_type get_from_array(const Armor::Array< element_type > &arr, uint idx)
Casts value stored in Armor::Array to return type.
bool equal_to(const return_type &other)
return_type & value_
static std::string type_name()
static const return_type & from_raw(return_type &val, ET *raw_data)
const ET * mem_ptr() const
ET & operator()(unsigned int i, unsigned int j)
internal::ReturnType< NRows, NCols, ET >::return_type return_type
FieldValue_(return_type &val)
static const int NRows_
unsigned int n_cols() const
static constexpr bool is_scalable()
unsigned int n_rows() const
void scale(double scale_coef)
void set_n_comp(unsigned int)
static const int NCols_
static const int rank_
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
EI_Address ei_address() const
Definition: accessors.cc:314
Iterator< ValueType > begin() const
unsigned int size() const
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
Class for declaration of the integral input data.
Definition: type_base.hh:483
Template for classes storing finite set of named values.
TYPEDEF_ERR_INFO(EI_InputMsg, const string)
unsigned int FieldEnum
Definition: field_values.hh:56
DECLARE_INPUT_EXCEPTION(ExcFV_Input,<< "Wrong field value input: "<< EI_InputMsg::val)
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
static constexpr bool value
Definition: json.hpp:87
unsigned int uint
std::string format(CStringRef format_str, ArgList args)
Definition: format.h:3141
void init_matrix_from_input(MatrixType &value, Input::Array rec)
RT & set_raw_fix(RT &val, double *raw_data)
RT & set_raw_vec(RT &val, double *raw_data)
void init_vector_from_input(VectorType &value, Input::Array rec)
RT & set_raw_scalar(RT &, double *raw_data)
FieldValue_< 0, 1, FieldEnum > EnumVector
FieldValue_< 0, 1, int > IntVector
FieldValue_< spacedim, spacedim, double > TensorFixed
FieldValue_< 1, 1, int > Integer
FieldValue_< 0, 1, double > Vector
FieldValue_< 1, 1, FieldEnum > Enum
FieldValue_< 1, 1, double > Scalar
FieldValue_< spacedim, 1, double > VectorFixed
Input::Type::Selection type
Definition: field_values.hh:79
Input::Type::Integer type
Definition: field_values.hh:76
Input::Type::Double type
Definition: field_values.hh:73
arma::Col< unsigned int > return_type
arma::Col< ET >::template fixed< NRows > return_type
Definition: field_values.hh:96
arma::Col< unsigned int >::template fixed< NRows > return_type
arma::Mat< unsigned int >::template fixed< NRows, NCols > return_type
arma::Mat< ET >::template fixed< NRows, NCols > return_type
Definition: field_values.hh:87