INMOST
Mathematical Modelling Toolkit
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
inmost_autodiff.h
Go to the documentation of this file.
1 
2 #ifndef INMOST_AUTODIFF_H_INCLUDED
3 #define INMOST_AUTODIFF_H_INCLUDED
4 #include "inmost_common.h"
5 #include "inmost_mesh.h"
6 #include "inmost_solver.h"
7 #include <sstream> //for debug
8 
9 //#define NEW_VERSION
10 
11 #if defined(USE_AUTODIFF) && (!defined(USE_MESH))
12 #warning "USE_AUTODIFF require USE_MESH"
13 #undef USE_AUTODIFF
14 #endif
15 
16 //#define DPRNT
17 
18 
19 
20 #if defined(USE_AUTODIFF)
21 #include <math.h>
22 
23 
24 #define FILTER_EPS 1e-12
25 
26 #define AD_SPACE 16384 //typical number of registered entries
27 
28 #define AD_NONE 0 //this should generate an error
29 //binary operations below
30 #define AD_PLUS 1 //sum x+y
31 #define AD_MINUS 2 //substraction x-y
32 #define AD_MULT 3 //multiplication x*y
33 #define AD_DIV 4 //devision x/y
34 #define AD_INV 5 //inversion 1/x, depricated
35 #define AD_POW 6 //power operation x^y
36 #define AD_SQRT 7 //square root operation
37 
38 
39 #define AD_PRECOMP 15 // this expression points to precomputed replacement, expr * left is position in array of precomputed values, expr * right points to original expression
40 
41 #define AD_VAL 19 // calculate only value of current expression
42 
43 //unary operations below
44 #define AD_COS 20 //cos(x)
45 
46 #define AD_ABS 22 // |x|
47 #define AD_EXP 23 // exp(x)
48 #define AD_LOG 24 // log(x)
49 #define AD_SIN 25 // sin(x)
50 
51 #define AD_CONST 50 // expr * left represents const value
52 #define AD_MES 51 // volume(e) for cell, area(e) for face, length(e) for edge
53 
54 #define AD_COND 60 //condition
55 #define AD_COND_TYPE 61 //condition on element type
56 #define AD_COND_MARK 62 //condition on marker
57 #define AD_ALTR 65 //alternatives for condition
58 
59 
60 
61 #define AD_EXT 99 // pointer to external expression data
62 
63 #define AD_TAG 100 //tag with dynamic data
64 #define AD_CTAG (AD_TAG+AD_SPACE) //tag with constant data
65 #define AD_STNCL (AD_CTAG+AD_SPACE) //stencil that sets current elements and their coefs
66 #define AD_TABLE (AD_STNCL+AD_SPACE) //table of values
67 #define AD_FUNC (AD_TABLE+AD_SPACE) //register function that calculates some value from element
68 
69 
70 //types of operands
71 
72 #define AD_TYPE_INVALID ENUMUNDEF
73 #define AD_TYPE_ENDPOINT 0
74 #define AD_TYPE_UNARY 1
75 #define AD_TYPE_BINARY 2
76 #define AD_TYPE_TERNARY 3
77 #define AD_TYPE_VALUE 4
78 
79 namespace INMOST
80 {
81  class Automatizator; //forward declaration
83  class expr;
84 
85 
86 
88  {
89  public:
90  typedef INMOST_DATA_REAL_TYPE(*func_callback)(const Storage & current_element, void * user_data);
91  typedef std::pair<HandleType, INMOST_DATA_REAL_TYPE> stencil_pair;
93  typedef void(*stencil_callback)(const Storage & current_element, stencil_pairs & out_stencil, void * user_data);
94  private:
95 #if defined(USE_OMP)
96  std::vector<Sparse::RowMerger> merger;
97 #else
98  Sparse::RowMerger merger;
99 #endif
100  typedef struct{ Tag t; MarkerType domain_mask; } tagdomain;
101  typedef dynarray<tagdomain, 128> const_tag_type;
102  typedef struct{ tagdomain d; Tag indices; } tagpair;
103  typedef dynarray<tagpair, 128> tagpairs_type;
104  typedef std::vector<tagpair> index_enum;
105  typedef struct { Tag elements, coefs; } stencil_tag;
106  typedef struct { std::string name; INMOST_DATA_ENUM_TYPE kind; MarkerType domainmask; void * link; } stencil_kind_domain;
107  typedef dynarray<stencil_kind_domain, 128> stencil_type;
108  typedef struct func_name_callback_t { std::string name; func_callback func; } func_name_callback;
109  typedef dynarray<func_name_callback, 128> func_type;
110  public:
111  typedef struct
112  {
113  std::string name;
117  INMOST_DATA_ENUM_TYPE binary_search(INMOST_DATA_REAL_TYPE arg);
119  INMOST_DATA_REAL_TYPE get_derivative(INMOST_DATA_REAL_TYPE arg);
120  std::pair<INMOST_DATA_REAL_TYPE, INMOST_DATA_REAL_TYPE> get_both(INMOST_DATA_REAL_TYPE arg);
121  } table;
122  typedef table * table_ptr;
123  private:
126  const_tag_type reg_ctags;
127  index_enum index_tags;
128  tagpairs_type reg_tags;
129  table_type reg_tables;
130  stencil_type reg_stencils;
131  func_type reg_funcs;
132  INMOST_DATA_ENUM_TYPE first_num;
133  INMOST_DATA_ENUM_TYPE last_num;
134  Mesh * m;
135 #if defined(NEW_VERSION)
136  INMOST_DATA_REAL_TYPE EvaluateSub(expr & var,INMOST_DATA_ENUM_TYPE element,INMOST_DATA_ENUM_TYPE parent, void * user_data);
137  void DerivativeFill(expr & var, INMOST_DATA_ENUM_TYPE element, INMOST_DATA_ENUM_TYPE parent, Sparse::RowMerger & entries, INMOST_DATA_REAL_TYPE multval, void * user_data);
138 #else
139  INMOST_DATA_REAL_TYPE DerivativePrecompute(const expr & var, const Storage & e, precomp_values_t & values, void * user_data);
140  void DerivativeFill(const expr & var, const Storage & e, Sparse::RowMerger & entries, precomp_values_t & values, INMOST_DATA_REAL_TYPE multval, void * user_data);
141 #endif
142  public:
143  Automatizator(Mesh * m);
144  ~Automatizator();
147  INMOST_DATA_ENUM_TYPE RegisterFunc(std::string name, func_callback func);
148  INMOST_DATA_ENUM_TYPE RegisterStencil(std::string name, Tag elements_tag, Tag coefs_tag, MarkerType domain_mask = 0);
149  INMOST_DATA_ENUM_TYPE RegisterStencil(std::string name, stencil_callback func, MarkerType domain_mask = 0);
151  INMOST_DATA_ENUM_TYPE RegisterDynamicTag(Tag t, ElementType typemask, MarkerType domain_mask = 0);
153  void EnumerateDynamicTags();
154  __INLINE Tag GetDynamicValueTag(INMOST_DATA_ENUM_TYPE ind) { return reg_tags[ind-AD_TAG].d.t; }
155  __INLINE Tag GetDynamicIndexTag(INMOST_DATA_ENUM_TYPE ind) { return reg_tags[ind-AD_TAG].indices; }
156  __INLINE MarkerType GetDynamicMask(INMOST_DATA_ENUM_TYPE ind) { return reg_tags[ind-AD_TAG].d.domain_mask; }
157  __INLINE Tag GetStaticValueTag(INMOST_DATA_ENUM_TYPE ind) { return reg_ctags[ind-AD_CTAG].t; }
158  __INLINE MarkerType GetStaticMask(INMOST_DATA_ENUM_TYPE ind) { return reg_ctags[ind-AD_CTAG].domain_mask; }
161  __INLINE bool isDynamicValid(const Storage & e, INMOST_DATA_ENUM_TYPE ind) { MarkerType mask = GetDynamicMask(ind); return mask == 0 || e->GetMarker(mask); }
163  __INLINE bool isStaticValid(const Storage & e, INMOST_DATA_ENUM_TYPE ind) { MarkerType mask = GetStaticMask(ind); return mask == 0 || e->GetMarker(mask); }
164 #if defined(NEW_VERSION)
165  INMOST_DATA_REAL_TYPE Evaluate(expr & var, const Storage & e, void * user_data);
166  INMOST_DATA_REAL_TYPE Derivative(expr & var, const Storage & e, Sparse::Row & out, Storage::real multiply, void * user_data);
167 #else
168  INMOST_DATA_REAL_TYPE Evaluate(const expr & var, const Storage & e, void * user_data);
169  INMOST_DATA_REAL_TYPE Derivative(const expr & var, const Storage & e, Sparse::Row & out, Storage::real multiply, void * user_data);
170 #endif
173  __INLINE Mesh * GetMesh() { return m; }
174  __INLINE INMOST_DATA_REAL_TYPE * GetTableArguments(INMOST_DATA_ENUM_TYPE tableind) {return reg_tables[tableind-AD_TABLE]->args;}
175  __INLINE INMOST_DATA_REAL_TYPE * GetTableValues(INMOST_DATA_ENUM_TYPE tableind) {return reg_tables[tableind-AD_TABLE]->vals;}
177  __INLINE INMOST_DATA_REAL_TYPE GetTableValue(INMOST_DATA_ENUM_TYPE tableind, INMOST_DATA_REAL_TYPE arg) { return reg_tables[tableind-AD_TABLE]->get_value(arg); }
178  __INLINE INMOST_DATA_REAL_TYPE GetTableDerivative(INMOST_DATA_ENUM_TYPE tableind, INMOST_DATA_REAL_TYPE arg) { return reg_tables[tableind-AD_TABLE]->get_derivative(arg); }
179  __INLINE std::pair<INMOST_DATA_REAL_TYPE, INMOST_DATA_REAL_TYPE> GetTableBoth(INMOST_DATA_ENUM_TYPE tableind, INMOST_DATA_REAL_TYPE arg) { return reg_tables[tableind-AD_TABLE]->get_both(arg); }
180  INMOST_DATA_ENUM_TYPE GetStencil(INMOST_DATA_ENUM_TYPE stnclind, const Storage & elem, void * user_data, stencil_pairs & ret);
181  __INLINE INMOST_DATA_REAL_TYPE GetFunction(INMOST_DATA_ENUM_TYPE funcid, const Storage & elem, void * user_data) { return reg_funcs[funcid-AD_FUNC].func(elem, user_data); }
183  {
184 #if defined(USE_OMP)
185  return merger[omp_get_thread_num()];
186 #else
187  return merger;
188 #endif
189  }
190  };
191 
192 #if defined(NEW_VERSION)
193  class expr
194  {
195 
196  typedef small_hash<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE, 1024> replace_type;
197 
198  class expr_data
199  {
200  typedef union { INMOST_DATA_REAL_TYPE r; INMOST_DATA_ENUM_TYPE i; expr * e; expr_data * q; } operand;
201  INMOST_DATA_ENUM_TYPE op, priority;
202  operand left, right;
203  public:
204  expr_data() : op(AD_NONE) { memset(&left, 0, sizeof(operand)); memset(&right, 0, sizeof(operand)); }
205  expr_data(const expr_data & other) : op(other.op)
206  {
207  if (other.op == AD_STNCL)
208  left.e = new expr(*other.left.e);
209  else if (other.op == AD_ALTR)
210  {
211  left.e = new expr(*other.left.e);
212  right.e = new expr(*other.right.e);
213  }
214  else if (other.op == AD_COND)
215  {
216  left = other.left;
217  right.q = new expr_data(*other.right.q);
218  }
219  else
220  {
221  left = other.left;
222  right = other.right;
223  }
224  }
225  expr_data & operator =(expr_data const & other)
226  {
227  op = other.op;
228  if (other.op == AD_STNCL)
229  left.e = new expr(*other.left.e);
230  else if (other.op == AD_ALTR)
231  {
232  left.e = new expr(*other.left.e);
233  right.e = new expr(*other.right.e);
234  }
235  else if (other.op == AD_COND)
236  {
237  left = other.left;
238  right.q = new expr_data(*other.right.q);
239  }
240  else
241  {
242  left = other.left;
243  right = other.right;
244  }
245  return *this;
246  }
247  expr_data(INMOST_DATA_REAL_TYPE val) : op(AD_CONST) { left.r = val; }
248  expr_data(INMOST_DATA_ENUM_TYPE new_op, INMOST_DATA_ENUM_TYPE l, INMOST_DATA_ENUM_TYPE r) : op(new_op) { left.i = l; right.i = r; }
249  expr_data(INMOST_DATA_ENUM_TYPE new_op, INMOST_DATA_ENUM_TYPE l, INMOST_DATA_REAL_TYPE r) : op(new_op) { left.i = l; right.r = r; }
250  expr_data(INMOST_DATA_ENUM_TYPE new_op, INMOST_DATA_ENUM_TYPE l, const expr & e) : op(new_op) { left.i = l; right.e = new expr(e); }
251  expr_data(INMOST_DATA_ENUM_TYPE l, const expr_data & e) : op(AD_COND) { left.i = l; right.q = new expr_data(e); }
252  expr_data(INMOST_DATA_ENUM_TYPE op, const expr * e) : op(op) { left.e = new expr(*e); }
253  expr_data(const expr * a, const expr * b) : op(AD_ALTR) { left.e = new expr(*a); right.e = new expr(*b); }
254  ~expr_data()
255  {
256  if (op == AD_STNCL)
257  {
258  delete left.e;
259  op = AD_NONE;
260  }
261  else if (op == AD_COND)
262  {
263  delete right.q;
264  op = AD_NONE;
265  }
266  else if (op == AD_ALTR)
267  {
268  delete left.e;
269  delete right.e;
270  op = AD_NONE;
271  }
272  }
273  bool operator ==(const expr_data & other) const
274  {
275  assert(op != AD_NONE);
276  if (op == AD_EXT) return left.e->data[right.i] == other;
277  if (other.op == AD_EXT) return *this == other.left.e->data[other.right.i];
278  if (op != other.op) return false;
279  switch (op)
280  {
281  case AD_MES: return true;
282  case AD_CONST: return fabs(left.r - other.left.r) < 1e-9;
283  case AD_PLUS:
284  case AD_MINUS:
285  case AD_MULT:
286  case AD_DIV:
287  case AD_POW:
288  return left.i == other.left.i && right.i == other.right.i;
289  case AD_COS:
290  case AD_ABS:
291  case AD_EXP:
292  case AD_LOG:
293  case AD_SIN:
294  case AD_SQRT:
295  return left.i == other.left.i;
296  case AD_VAL:
297  return left.i == other.right.i && (right.e == other.right.e || *right.e == *other.right.e);
298  case AD_ALTR:
299  if (left.e == other.left.e && right.e == other.right.e) return true;
300  return *left.e == *other.left.e && *right.e == *other.right.e;
301  case AD_COND:
302  return left.i == other.left.i && (right.q == other.right.q || *right.q == *other.right.q);
303  case AD_COND_TYPE:
304  case AD_COND_MARK:
305  return left.i == other.left.i;
306  }
307  if (op >= AD_FUNC) return true;
308  if (op >= AD_TABLE) return left.i == other.left.i;
309  if (op >= AD_STNCL) return (left.e == other.left.e) || (*left.e == *other.left.e);
310  if (op >= AD_CTAG) return left.i == other.left.i;
311  if (op >= AD_TAG) return left.i == other.left.i;
312  assert(false); //cannot reach here
313  return false;
314  }
315  bool deep_cmp(const expr_data & other, const expr & my_set, const expr & other_set ) const
316  {
317  assert(op != AD_NONE);
318 
319  if (op == AD_EXT) return left.e->data[right.i].deep_cmp(other,*left.e,other_set);
320  if (other.op == AD_EXT) return deep_cmp(other.left.e->data[other.right.i],my_set, *other.left.e);
321  if (op != other.op) return false;
322  switch (op)
323  {
324  case AD_MES: return true;
325  case AD_CONST: return fabs(left.r - other.left.r) < 1e-9;
326  case AD_PLUS:
327  case AD_MINUS:
328  case AD_MULT:
329  case AD_DIV:
330  case AD_POW:
331  return my_set.data[left.i].deep_cmp(other_set.data[other.left.i],my_set,other_set)
332  && my_set.data[right.i].deep_cmp(other_set.data[other.right.i],my_set,other_set);
333  case AD_COS:
334  case AD_ABS:
335  case AD_EXP:
336  case AD_LOG:
337  case AD_SIN:
338  case AD_SQRT:
339  return my_set.data[left.i].deep_cmp(other_set.data[other.left.i], my_set, other_set);
340  case AD_VAL:
341  return (right.e == other.right.e || *right.e == *other.right.e) && my_set.data[left.i].deep_cmp(other_set.data[other.left.i], my_set, other_set);
342  case AD_ALTR:
343  if (left.e == other.left.e && right.e == other.right.e) return true;
344  return *left.e == *other.left.e && *right.e == *other.right.e;
345  case AD_COND:
346  return my_set.data[left.i].deep_cmp(other_set.data[other.left.i], my_set, other_set) && (right.q == other.right.q || *right.q == *other.right.q);
347  case AD_COND_TYPE:
348  case AD_COND_MARK:
349  return left.i == other.left.i;
350  }
351  if (op >= AD_FUNC) return true;
352  if (op >= AD_TABLE) return my_set.data[left.i].deep_cmp(other_set.data[other.left.i], my_set, other_set);
353  if (op >= AD_STNCL) return (left.e == other.left.e) || (*left.e == *other.left.e);
354  if (op >= AD_CTAG) return left.i == other.left.i;
355  if (op >= AD_TAG) return left.i == other.left.i;
356  assert(false); //cannot reach here
357  return false;
358  }
359  bool cmp(const expr_data & other,replace_type & r) const
360  {
361  assert(op != AD_NONE);
362  if (op == AD_EXT) return left.e->data[right.i] == other;
363  if (other.op == AD_EXT) return *this == other.left.e->data[other.right.i];
364  if (op != other.op) return false;
365  switch (op)
366  {
367  case AD_MES: return true;
368  case AD_CONST: return fabs(left.r - other.left.r) < 1e-9;
369  case AD_PLUS:
370  case AD_MINUS:
371  case AD_MULT:
372  case AD_DIV:
373  case AD_POW:
374  return left.i == r[other.left.i] && right.i == r[other.right.i];
375  case AD_COS:
376  case AD_ABS:
377  case AD_EXP:
378  case AD_LOG:
379  case AD_SIN:
380  case AD_SQRT:
381  return left.i == r[other.left.i];
382  case AD_VAL:
383  return left.i == r[other.left.i] && (right.e == other.right.e || *right.e == *other.right.e);
384  case AD_ALTR:
385  if (left.e == other.left.e && right.e == other.right.e) return true;
386  return *left.e == *other.left.e && *right.e == *other.right.e;
387  case AD_COND:
388  return left.i == r[other.left.i] && (right.q == other.right.q || *right.q == *other.right.q);
389  case AD_COND_TYPE:
390  case AD_COND_MARK:
391  return left.i == other.left.i;
392  }
393  if (op >= AD_FUNC) return true;
394  if (op >= AD_TABLE) return left.i == r[other.left.i];
395  if (op >= AD_STNCL) return (left.e == other.left.e) || (*left.e == *other.left.e);
396  if (op >= AD_CTAG) return left.i == other.left.i;
397  if (op >= AD_TAG) return left.i == other.left.i;
398  assert(false); //cannot reach here
399  return false;
400  }
401  INMOST_DATA_ENUM_TYPE type() const
402  {
403  assert(op != AD_NONE);
404  if (op == AD_EXT) return left.e->data[right.i].type();
405  switch (op)
406  {
407  case AD_COND_MARK:
408  case AD_COND_TYPE:
409  case AD_CONST:
410  case AD_MES: return AD_TYPE_ENDPOINT;
411  case AD_PLUS:
412  case AD_MINUS:
413  case AD_MULT:
414  case AD_DIV:
415  case AD_POW:
416  case AD_ALTR: return AD_TYPE_BINARY;
417  case AD_COS:
418  case AD_ABS:
419  case AD_EXP:
420  case AD_LOG:
421  case AD_SQRT:
422  case AD_SIN: return AD_TYPE_UNARY;
423  case AD_COND: return AD_TYPE_TERNARY;
424  case AD_VAL: return AD_TYPE_VALUE;
425  }
426  if (op >= AD_FUNC) return AD_TYPE_ENDPOINT;
427  if (op >= AD_TABLE) return AD_TYPE_UNARY;
428  if (op >= AD_STNCL) return AD_TYPE_UNARY;
429  if (op >= AD_CTAG) return AD_TYPE_ENDPOINT;
430  if (op >= AD_TAG) return AD_TYPE_ENDPOINT;
431  assert(false); //cannot reach here
432  return AD_TYPE_INVALID;
433  }
434  bool is_func() { return op >= AD_FUNC; }
435  bool is_table() { return op >= AD_TABLE && op < AD_FUNC; }
436  bool is_stncl() { return op >= AD_STNCL && op < AD_TABLE; }
437  bool is_ctag() { return op >= AD_CTAG && op < AD_STNCL; }
438  bool is_tag() { return op >= AD_TAG && op < AD_CTAG; }
439  friend class expr;
440  friend class Automatizator;
441  };
442  friend class expr_data;
443  //typedef std::vector<expr_data> data_type;
444  typedef std::vector<expr_data> data_type;
445  data_type data;
446  std::vector<INMOST_DATA_REAL_TYPE> values;
447  Automatizator::stencil_pairs current_stencil;
448  __INLINE INMOST_DATA_ENUM_TYPE values_offset(INMOST_DATA_ENUM_TYPE element) {return element * static_cast<INMOST_DATA_ENUM_TYPE>(data.size());}
449  __INLINE INMOST_DATA_ENUM_TYPE derivatives_offset(INMOST_DATA_ENUM_TYPE element) {return element * static_cast<INMOST_DATA_ENUM_TYPE>(data.size()) + static_cast<INMOST_DATA_ENUM_TYPE>(data.size()*current_stencil.size());}
450  void resize_for_stencil()
451  {
452  values.resize(current_stencil.size()*data.size()*2);
453  memset(&values[current_stencil.size()*data.size()],0,sizeof(INMOST_DATA_REAL_TYPE)*current_stencil.size()*data.size());
454  }
455 
456  void relink_data()
457  {
458  for (data_type::iterator it = data.begin(); it != data.end(); ++it)
459  if (it->op == AD_COND )
460  {
461  for (data_type::iterator jt = it->right.q->left.e->data.begin(); jt != it->right.q->left.e->data.end(); ++jt)
462  {
463  if (jt->op == AD_EXT) jt->left.e = this;
464  }
465  it->right.q->left.e->relink_data();
466  for (data_type::iterator jt = it->right.q->right.e->data.begin(); jt != it->right.q->right.e->data.end(); ++jt)
467  {
468  if (jt->op == AD_EXT) jt->left.e = this;
469  }
470  it->right.q->left.e->relink_data();
471  }
472  }
473  void link_data(data_type & other, replace_type & r)
474  {
475  INMOST_DATA_ENUM_TYPE i = 0, j;
476  for (data_type::iterator it = other.begin(); it != other.end(); ++it)
477  {
478  j = 0;
479  if (it->op == AD_EXT)
480  {
481  it->left.e = this;
482  it->right.i = r[it->right.i];
483  }
484  else for (data_type::iterator jt = data.begin(); jt != data.end(); ++jt)
485  {
486  if (*it == *jt)
487  {
488  it->op = AD_EXT;
489  it->left.e = this;
490  it->right.i = j;
491  break;
492  }
493  j++;
494  }
495  i++;
496  }
497  }
498  void replace_operands(expr_data & e, replace_type & r)
499  {
500  if (e.op == AD_EXT) replace_operands(e.left.e->data[e.right.i],r);
501  else if (e.op == AD_ALTR)
502  {
503  link_data(e.left.e->data,r);
504  link_data(e.right.e->data,r);
505  }
506  else if (e.is_stncl()) {}
507  else switch (e.type())
508  {
509  case AD_TYPE_UNARY:
510  assert(r.is_present(e.left.i));
511  e.left.i = r[e.left.i];
512  break;
513  case AD_TYPE_BINARY:
514  assert(r.is_present(e.left.i));
515  assert(r.is_present(e.right.i));
516  e.left.i = r[e.left.i];
517  e.right.i = r[e.right.i];
518  break;
519  case AD_TYPE_TERNARY:
520  assert(r.is_present(e.left.i));
521  e.left.i = r[e.left.i];
522  replace_operands(*e.right.q, r);
523  break;
524  case AD_TYPE_VALUE:
525  e.left.i = r[e.left.i];
526  link_data(e.right.e->data,r);
527  break;
528  }
529  }
530  INMOST_DATA_ENUM_TYPE merge_data(const data_type & other)
531  {
532  replace_type from_to;
533  INMOST_DATA_ENUM_TYPE i = 0, j;
534  for (data_type::const_iterator it = other.begin(); it != other.end(); ++it)
535  {
536  j = 0;
537  for (data_type::iterator jt = data.begin(); jt != data.end(); ++jt)
538  {
539  if (jt->cmp(*it,from_to))
540  {
541  from_to[i] = j;
542  break;
543  }
544  j++;
545  }
546  if (!from_to.is_present(i))
547  {
548  from_to[i] = static_cast<INMOST_DATA_ENUM_TYPE>(data.size());
549  data.push_back(*it);
550  replace_operands(data.back(),from_to);
551  }
552  i++;
553  }
554  return from_to[static_cast<int>(other.size()) - 1];
555  }
556 
557  public:
558  expr() : data() { }
559  expr(const expr & other) : data(other.data) { relink_data(); }
560  expr(INMOST_DATA_ENUM_TYPE op, const expr * sub) :data() { data.push_back(expr_data(op,sub)); }
561  expr(const expr * l, const expr * r) :data() { data.push_back(expr_data(l,r)); }
562  expr(INMOST_DATA_REAL_TYPE val) : data() { data.push_back(expr_data(val)); }
563  expr(INMOST_DATA_ENUM_TYPE op, INMOST_DATA_ENUM_TYPE comp) : data() { data.push_back(expr_data(op, comp, ENUMUNDEF)); }
564  expr(INMOST_DATA_ENUM_TYPE op, const expr & operand) : data(operand.data) { relink_data(); data.push_back(expr_data(op, static_cast<INMOST_DATA_ENUM_TYPE>(data.size() - 1), ENUMUNDEF)); }
565  expr(const expr & operand, const expr & multiplyer) : data(operand.data) { relink_data(); data.push_back(expr_data(AD_VAL, static_cast<INMOST_DATA_ENUM_TYPE>(data.size() - 1), multiplyer)); }
566  expr(const expr & cond, const expr & if_true, const expr & if_false) :data(cond.data)
567  {
568  relink_data();
569  data.push_back(expr_data(static_cast<INMOST_DATA_ENUM_TYPE>(data.size() - 1), expr_data(&if_true, &if_false)));
570  }
571  expr(INMOST_DATA_ENUM_TYPE op, const expr & l, const expr & r) :data(l.data)
572  {
573  relink_data();
574  INMOST_DATA_ENUM_TYPE lp = static_cast<INMOST_DATA_ENUM_TYPE>(data.size() - 1);
575  INMOST_DATA_ENUM_TYPE rp = merge_data(r.data);
576  data.push_back(expr_data(op, lp, rp));
577  }
578  ~expr() {}
579  expr & operator =(expr const & other) {data = other.data; return *this;}
580  expr operator +() const { return expr(*this); }
581  expr operator -() const { expr zero(0.0); return zero - *this; }
582  expr operator +(const expr & other) const {return expr(AD_PLUS,*this,other);}
583  expr operator -(const expr & other) const { return expr(AD_MINUS, *this, other); }
584  expr operator *(const expr & other) const { return expr(AD_MULT, *this, other); }
585  expr operator /(const expr & other) const { return expr(AD_DIV, *this, other); }
586  expr operator +(INMOST_DATA_REAL_TYPE other) const { return expr(AD_PLUS, *this, expr(other)); }
587  expr operator -(INMOST_DATA_REAL_TYPE other) const { return expr(AD_MINUS, *this, expr(other)); }
588  expr operator *(INMOST_DATA_REAL_TYPE other) const { return expr(AD_MULT, *this, expr(other)); }
589  expr operator /(INMOST_DATA_REAL_TYPE other) const { return expr(AD_DIV, *this, expr(other)); }
590  bool operator ==(const expr & other) //this should account for possible reorder! here x*y is not the same as y*x
591  {
592  if (data.size() != other.data.size()) return false;
593  data_type::iterator it = data.begin();
594  data_type::const_iterator jt = other.data.begin();
595  while (it != data.end() && jt != other.data.end())
596  {
597  if (!it->deep_cmp(*jt, *this, other)) return false;
598  ++it; ++jt;
599  }
600  return true;
601  }
602 
603  friend class Automatizator;
604  };
605 }
606 
607 __INLINE INMOST::expr operator +(INMOST_DATA_REAL_TYPE left, const INMOST::expr & right) { return INMOST::expr(left) + right; }
608 __INLINE INMOST::expr operator -(INMOST_DATA_REAL_TYPE left, const INMOST::expr & right) { return INMOST::expr(left) - right; }
609 __INLINE INMOST::expr operator *(INMOST_DATA_REAL_TYPE left, const INMOST::expr & right) { return INMOST::expr(left) * right; }
610 __INLINE INMOST::expr operator /(INMOST_DATA_REAL_TYPE left, const INMOST::expr & right) { return INMOST::expr(left) / right; }
611 __INLINE INMOST::expr ad_pow(const INMOST::expr & v, const INMOST::expr & n) { return INMOST::expr(AD_POW, v, n); }
612 //__INLINE INMOST::expr ad_pow(const INMOST::expr & v, INMOST_DATA_REAL_TYPE n) { return INMOST::expr(AD_POW, v, INMOST::expr(n)); }
613 //__INLINE INMOST::expr ad_pow(const INMOST::expr & v, INMOST_DATA_ENUM_TYPE n) { return INMOST::expr(AD_POW, v, INMOST::expr(static_cast<INMOST_DATA_REAL_TYPE>(n))); }
614 __INLINE INMOST::expr ad_abs(const INMOST::expr & v) { return INMOST::expr(AD_ABS, v); }
615 __INLINE INMOST::expr ad_exp(const INMOST::expr & v) { return INMOST::expr(AD_EXP, v); }
616 __INLINE INMOST::expr ad_log(const INMOST::expr & v) { return INMOST::expr(AD_LOG, v); }
617 __INLINE INMOST::expr ad_sin(const INMOST::expr & v) { return INMOST::expr(AD_SIN, v); }
618 __INLINE INMOST::expr ad_cos(const INMOST::expr & v) { return INMOST::expr(AD_COS, v); }
620 __INLINE INMOST::expr ad_val(const INMOST::expr & v, const INMOST::expr & multiplyer = INMOST::expr(0.0)) {return INMOST::expr(v, multiplyer);}
622 __INLINE INMOST::expr condition(const INMOST::expr & cond, const INMOST::expr & if_gt_zero, const INMOST::expr & if_le_zero) { return INMOST::expr(cond, if_gt_zero, if_le_zero); }
623 __INLINE INMOST::expr condition_etype(INMOST::ElementType etypes, const INMOST::expr & if_true, const INMOST::expr & if_false) { return INMOST::expr(INMOST::expr(AD_COND_TYPE, etypes), if_true, if_false); }
624 __INLINE INMOST::expr condition_marker(INMOST::MarkerType marker, const INMOST::expr & if_true, const INMOST::expr & if_false) { return INMOST::expr(INMOST::expr(AD_COND_MARK, marker), if_true, if_false); }
625 __INLINE INMOST::expr stencil(INMOST_DATA_ENUM_TYPE stncl, const INMOST::expr & v) { assert(stncl >= AD_STNCL && stncl < AD_TABLE); return INMOST::expr(stncl, &v); }
626 __INLINE INMOST::expr tabval(INMOST_DATA_ENUM_TYPE tabl, const INMOST::expr & v) { assert(tabl >= AD_TABLE && tabl < AD_FUNC); return INMOST::expr(tabl, v); }
627 __INLINE INMOST::expr tagval(INMOST_DATA_ENUM_TYPE reg_tag, INMOST_DATA_ENUM_TYPE comp = 0) { assert(reg_tag >= AD_TAG && reg_tag < AD_STNCL); return INMOST::expr(reg_tag, comp); }
628 __INLINE INMOST::expr funcval(INMOST_DATA_ENUM_TYPE reg_func) { assert(reg_func >= AD_FUNC); return INMOST::expr(reg_func, ENUMUNDEF); }
629 namespace INMOST
630 {
631 #else
632 }
633 
638 
639 namespace INMOST
640 {
641  class expr
642  {
645  expr * left;
646  expr * right;
647  public:
648  expr() : op(AD_NONE), coef(1), left(NULL), right(NULL) { }
649  expr(const expr & other) : op(other.op), coef(other.coef)
650  {
651  if ((other.op >= AD_TAG && other.op < AD_STNCL) || other.op == AD_COND_TYPE || other.op == AD_COND_MARK) left = other.left; //copy component number
652  else if (other.left != NULL) left = new expr(*other.left); else left = NULL;
653  if (other.right != NULL) right = new expr(*other.right); else right = NULL;
654  }
655  expr(INMOST_DATA_REAL_TYPE val) : op(AD_CONST), coef(val), left(NULL), right(NULL) {}
656  //expr(INMOST_DATA_ENUM_TYPE val) : op(AD_CONST), coef(val), left(NULL), right(NULL) {}
657  expr(INMOST_DATA_ENUM_TYPE new_op, expr * l, expr * r) : op(new_op), coef(1), left(l), right(r) {}
658  expr(INMOST_DATA_ENUM_TYPE tag_op, INMOST_DATA_ENUM_TYPE comp) :op(tag_op), coef(1), right(NULL)
659  {
660  *(INMOST_DATA_ENUM_TYPE *)(&left) = comp;
661  }
663  {
664  if (op < AD_COS)
665  {
666  delete left;
667  delete right;
668  }
669  else if (op < AD_CONST)
670  delete left;
671  else if (op >= AD_STNCL && op < AD_TABLE)
672  delete left;
673  }
674  expr & operator =(expr const & other)
675  {
676  op = other.op; coef = other.coef;
677  if (other.op >= AD_TAG && other.op < AD_STNCL)
678  {
679  left = other.left; //copy component number
680  right = other.right;
681  }
682  else if (other.left != NULL) left = new expr(*other.left); else left = NULL;
683  if (other.right != NULL) right = new expr(*other.right); else right = NULL;
684  return *this;
685  }
686  expr operator +() { return expr(*this); }
687  expr operator -() { expr var(*this); var.coef *= -1.0; return var; }
688  expr operator +(const expr & other) const { return expr(AD_PLUS, new expr(*this), new expr(other)); }
689  expr operator -(const expr & other) const { return expr(AD_MINUS, new expr(*this), new expr(other)); }
690  expr operator *(const expr & other) const { return expr(AD_MULT, new expr(*this), new expr(other)); }
691  expr operator /(const expr & other) const { return expr(AD_DIV, new expr(*this), new expr(other)); }
692  expr operator +(const INMOST_DATA_REAL_TYPE & other) const { return expr(AD_PLUS, new expr(*this), new expr(other)); }
693  expr operator -(const INMOST_DATA_REAL_TYPE & other) const { return expr(AD_MINUS, new expr(*this), new expr(other)); }
694  expr operator *(const INMOST_DATA_REAL_TYPE & other) const { expr var(*this); var.coef *= other; return var; }
695  expr operator /(const INMOST_DATA_REAL_TYPE & other) const { expr var(*this); var.coef /= other; return var; }
696  __INLINE friend expr (::operator +)(const INMOST_DATA_REAL_TYPE & left, const expr & right);
697  __INLINE friend expr (::operator -)(const INMOST_DATA_REAL_TYPE & left, const expr & right);
698  __INLINE friend expr (::operator *)(const INMOST_DATA_REAL_TYPE & left, const expr & right);
699  __INLINE friend expr (::operator /)(const INMOST_DATA_REAL_TYPE & left, const expr & right);
700  bool operator ==(const expr & other) {return this == &other || (op == other.op && left == other.left && right == other.right);}
701  bool is_endpoint() { return op >= AD_TAG && op < AD_STNCL; }
702  friend class Automatizator;
703  };
704 
705 }
706 
707 __INLINE INMOST::expr operator +(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right) { return INMOST::expr(AD_PLUS, new INMOST::expr(left), new INMOST::expr(right)); }
708 __INLINE INMOST::expr operator -(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right) { return INMOST::expr(AD_MINUS, new INMOST::expr(left), new INMOST::expr(right)); }
709 __INLINE INMOST::expr operator *(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right) { INMOST::expr var(right); var.coef *= left; return var; }
710 __INLINE INMOST::expr operator /(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right) { INMOST::expr var; var.op = AD_INV; var.coef = left; var.left = new INMOST::expr(right); var.right = NULL; return var; }
718 __INLINE INMOST::expr ad_val(const INMOST::expr & v, const INMOST::expr & multiplyer = INMOST::expr(0.0)) {return INMOST::expr(AD_VAL,new INMOST::expr(v), new INMOST::expr(multiplyer));}
719 __INLINE INMOST::expr measure() { return INMOST::expr(AD_MES, NULL, NULL); }
720 __INLINE INMOST::expr condition_etype(INMOST::ElementType etype, const INMOST::expr & if_true, const INMOST::expr & if_false) { return INMOST::expr(AD_COND, new INMOST::expr(AD_COND_TYPE,etype), new INMOST::expr(AD_ALTR, new INMOST::expr(if_true), new INMOST::expr(if_false))); }
721 __INLINE INMOST::expr condition_marker(INMOST::MarkerType marker, const INMOST::expr & if_true, const INMOST::expr & if_false) { return INMOST::expr(AD_COND, new INMOST::expr(AD_COND_MARK,marker), new INMOST::expr(AD_ALTR, new INMOST::expr(if_true), new INMOST::expr(if_false))); }
722 __INLINE INMOST::expr condition(const INMOST::expr & cond, const INMOST::expr & if_gt_zero, const INMOST::expr & if_le_zero) { return INMOST::expr(AD_COND, new INMOST::expr(cond), new INMOST::expr(AD_ALTR, new INMOST::expr(if_gt_zero), new INMOST::expr(if_le_zero))); }
723 __INLINE INMOST::expr stencil(INMOST_DATA_ENUM_TYPE stncl, const INMOST::expr & v) { assert(stncl >= AD_STNCL && stncl < AD_TABLE); return INMOST::expr(stncl, new INMOST::expr(v), NULL); }
724 __INLINE INMOST::expr tabval(INMOST_DATA_ENUM_TYPE tabl, const INMOST::expr & v) { assert(tabl >= AD_TABLE && tabl < AD_FUNC); return INMOST::expr(tabl, new INMOST::expr(v), NULL); }
725 __INLINE INMOST::expr tagval(INMOST_DATA_ENUM_TYPE reg_tag, INMOST_DATA_ENUM_TYPE comp = 0) { assert(reg_tag >= AD_TAG && reg_tag < AD_STNCL); return INMOST::expr(reg_tag, comp); }
726 __INLINE INMOST::expr funcval(INMOST_DATA_ENUM_TYPE reg_func) { assert(reg_func >= AD_FUNC); return INMOST::expr(reg_func, NULL, NULL); }
727 namespace INMOST
728 {
729 #endif
730 }
731 
732 #endif
733 #endif
__INLINE INMOST::expr operator-(const INMOST_DATA_REAL_TYPE &left, const INMOST::expr &right)
__INLINE INMOST::expr ad_exp(const INMOST::expr &v)
#define AD_EXT
__INLINE INMOST_DATA_REAL_TYPE * GetTableArguments(INMOST_DATA_ENUM_TYPE tableind)
__INLINE INMOST_DATA_REAL_TYPE GetFunction(INMOST_DATA_ENUM_TYPE funcid, const Storage &elem, void *user_data)
__INLINE INMOST::unary_custom_variable< INMOST::abs_expression< typename A::Var >, A > fabs(INMOST::shell_dynamic_variable< typename A::Var, A > const &Arg)
#define AD_ABS
__INLINE INMOST_DATA_REAL_TYPE get_value(INMOST::shell_expression< A > const &Arg)
__INLINE INMOST_DATA_REAL_TYPE GetDynamicValue(const Storage &e, INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_ENUM_TYPE comp=0)
#define AD_TYPE_UNARY
#define AD_INV
#define AD_VAL
__INLINE INMOST_DATA_REAL_TYPE GetStaticValue(const Storage &e, INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_ENUM_TYPE comp=0)
__INLINE MarkerType GetStaticMask(INMOST_DATA_ENUM_TYPE ind)
__INLINE integer_array IntegerArray(const Tag &tag) const
Definition: inmost_mesh.h:3266
__INLINE INMOST::expr tabval(INMOST_DATA_ENUM_TYPE tabl, const INMOST::expr &v)
#define AD_POW
#define AD_TYPE_TERNARY
Class to store the sparse matrix row.
Definition: inmost_sparse.h:90
__INLINE MarkerType GetDynamicMask(INMOST_DATA_ENUM_TYPE ind)
dynarray< INMOST_DATA_REAL_TYPE, 2048 > precomp_values_t
#define AD_FUNC
#define AD_MES
__INLINE std::pair< INMOST_DATA_REAL_TYPE, INMOST_DATA_REAL_TYPE > GetTableBoth(INMOST_DATA_ENUM_TYPE tableind, INMOST_DATA_REAL_TYPE arg)
__INLINE INMOST_DATA_ENUM_TYPE GetComponents(const Storage &e, INMOST_DATA_ENUM_TYPE tagind)
__INLINE INMOST::expr operator+(const INMOST_DATA_REAL_TYPE &left, const INMOST::expr &right)
__INLINE INMOST_DATA_ENUM_TYPE GetLastIndex()
INMOST_DATA_ENUM_TYPE GetStencil(INMOST_DATA_ENUM_TYPE stnclind, const Storage &elem, void *user_data, stencil_pairs &ret)
bool operator==(const expr &other)
__INLINE real_array RealArray(const Tag &tag) const
Definition: inmost_mesh.h:3262
#define AD_TABLE
__INLINE friend const expr & right
#define AD_COND
__INLINE INMOST::expr ad_sin(const INMOST::expr &v)
#define INMOST_DATA_REAL_TYPE
expr(INMOST_DATA_ENUM_TYPE new_op, expr *l, expr *r)
#define AD_LOG
expr(INMOST_DATA_ENUM_TYPE tag_op, INMOST_DATA_ENUM_TYPE comp)
__INLINE INMOST::expr ad_pow(const INMOST::expr &v, const INMOST::expr n)
__INLINE INMOST_DATA_REAL_TYPE GetIndex(const Storage &e, INMOST_DATA_ENUM_TYPE tagind, INMOST_DATA_ENUM_TYPE comp=0)
#define AD_TYPE_VALUE
__INLINE INMOST::expr ad_val(const INMOST::expr &v, const INMOST::expr &multiplyer=INMOST::expr(0.0))
__INLINE INMOST_DATA_ENUM_TYPE GetFirstIndex()
#define AD_COND_MARK
#define AD_PLUS
dynarray< stencil_pair, 64 > stencil_pairs
INMOST_DATA_BULK_TYPE ElementType
Definition: inmost_data.h:23
__INLINE INMOST_DATA_REAL_TYPE * GetTableValues(INMOST_DATA_ENUM_TYPE tableind)
__INLINE INMOST::expr operator/(const INMOST_DATA_REAL_TYPE &left, const INMOST::expr &right)
__INLINE Tag GetDynamicValueTag(INMOST_DATA_ENUM_TYPE ind)
#define AD_DIV
__INLINE INMOST::expr stencil(INMOST_DATA_ENUM_TYPE stncl, const INMOST::expr &v)
#define AD_COS
__INLINE bool isDynamicValid(const Storage &e, INMOST_DATA_ENUM_TYPE ind)
__INLINE INMOST_DATA_REAL_TYPE GetTableValue(INMOST_DATA_ENUM_TYPE tableind, INMOST_DATA_REAL_TYPE arg)
__INLINE INMOST::expr funcval(INMOST_DATA_ENUM_TYPE reg_func)
__INLINE INMOST::expr ad_abs(const INMOST::expr &v)
INMOST_DATA_ENUM_TYPE RegisterDynamicTag(Tag t, ElementType typemask, MarkerType domain_mask=0)
__INLINE size_type size() const
Definition: container.hpp:765
#define AD_TAG
__INLINE Tag GetDynamicIndexTag(INMOST_DATA_ENUM_TYPE ind)
#define AD_CTAG
Sparse::RowMerger & GetMerger()
#define AD_ALTR
__INLINE INMOST_DATA_ENUM_TYPE GetTableSize(INMOST_DATA_ENUM_TYPE tableind)
__INLINE size_type size() const
Definition: container.hpp:1780
INMOST_DATA_REAL_TYPE Derivative(const expr &var, const Storage &e, Sparse::Row &out, Storage::real multiply, void *user_data)
#define AD_TYPE_INVALID
#define __INLINE
Definition: inmost_common.h:75
__INLINE INMOST::expr ad_cos(const INMOST::expr &v)
__INLINE bool isStaticValid(const Storage &e, INMOST_DATA_ENUM_TYPE ind)
__INLINE bool GetMarker(MarkerType n) const
Definition: inmost_mesh.h:3438
expr operator/(const expr &other) const
expr(INMOST_DATA_REAL_TYPE val)
#define AD_TYPE_BINARY
#define ENUMUNDEF
__INLINE INMOST::expr condition(const INMOST::expr &cond, const INMOST::expr &if_gt_zero, const INMOST::expr &if_le_zero)
#define AD_TYPE_ENDPOINT
__INLINE INMOST::expr ad_sqrt(const INMOST::expr &v)
#define INMOST_DATA_ENUM_TYPE
__INLINE INMOST::expr tagval(INMOST_DATA_ENUM_TYPE reg_tag, INMOST_DATA_ENUM_TYPE comp=0)
INMOST_DATA_REAL_TYPE real
Storage type for representing real values.
Definition: inmost_data.h:278
#define AD_COND_TYPE
#define AD_STNCL
__INLINE INMOST::expr operator*(const INMOST_DATA_REAL_TYPE &left, const INMOST::expr &right)
std::pair< HandleType, INMOST_DATA_REAL_TYPE > stencil_pair
INMOST_DATA_REAL_TYPE * vals
__INLINE Mesh * GetMesh()
INMOST_DATA_ENUM_TYPE RegisterStaticTag(Tag t, MarkerType domain_mask=0)
#define AD_NONE
__INLINE INMOST_DATA_REAL_TYPE GetTableDerivative(INMOST_DATA_ENUM_TYPE tableind, INMOST_DATA_REAL_TYPE arg)
INMOST_DATA_ENUM_TYPE RegisterStencil(std::string name, Tag elements_tag, Tag coefs_tag, MarkerType domain_mask=0)
#define AD_CONST
void(* stencil_callback)(const Storage &current_element, stencil_pairs &out_stencil, void *user_data)
__INLINE INMOST::expr ad_log(const INMOST::expr &v)
#define AD_MULT
__INLINE INMOST_DATA_ENUM_TYPE GetDynamicIndex(const Storage &e, INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_ENUM_TYPE comp=0)
INMOST_DATA_REAL_TYPE Evaluate(const expr &var, const Storage &e, void *user_data)
INMOST_DATA_ENUM_TYPE size
INMOST_DATA_REAL_TYPE(* func_callback)(const Storage &current_element, void *user_data)
This class provides the access to the individual mesh datum and general information about it...
Definition: inmost_data.h:185
INMOST_DATA_REAL_TYPE * args
__INLINE INMOST::expr condition_etype(INMOST::ElementType etype, const INMOST::expr &if_true, const INMOST::expr &if_false)
expr operator*(const expr &other) const
__INLINE Tag GetStaticValueTag(INMOST_DATA_ENUM_TYPE ind)
#define AD_MINUS
friend class Automatizator
expr(const expr &other)
INMOST_DATA_ENUM_TYPE MarkerType
Low 8 bits - marker mask, rest high bits - position of marker.
Definition: inmost_data.h:62
__INLINE INMOST::expr condition_marker(INMOST::MarkerType marker, const INMOST::expr &if_true, const INMOST::expr &if_false)
INMOST_DATA_ENUM_TYPE RegisterTable(std::string name, INMOST_DATA_REAL_TYPE *Arguments, INMOST_DATA_REAL_TYPE *Values, INMOST_DATA_ENUM_TYPE size)
expr & operator=(expr const &other)
__INLINE INMOST::expr measure()
INMOST_DATA_ENUM_TYPE RegisterFunc(std::string name, func_callback func)
#define AD_EXP
#define AD_SQRT
#define AD_SIN