GetFEM  5.4.1
getfem_generic_assembly_tree.cc
1 /*===========================================================================
2 
3  Copyright (C) 2013-2020 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 #include "getfem/getfem_generic_assembly_tree.h"
23 #include "getfem/getfem_generic_assembly_functions_and_operators.h"
24 #include "getfem/getfem_generic_assembly_compile_and_exec.h"
25 
26 
27 namespace getfem {
28 
29  //=========================================================================
30  // Lexical analysis for the generic assembly language
31  //=========================================================================
32 
33  static GA_TOKEN_TYPE ga_char_type[256];
34  static int ga_operator_priorities[GA_NB_TOKEN_TYPE];
35 
36  // Initialize ga_char_type and ga_operator_priorities arrays
37  static bool init_ga_char_type() {
38  for (int i = 0; i < 256; ++i) ga_char_type[i] = GA_INVALID;
39  ga_char_type['+'] = GA_PLUS; ga_char_type['-'] = GA_MINUS;
40  ga_char_type['*'] = GA_MULT; ga_char_type['/'] = GA_DIV;
41  ga_char_type[':'] = GA_COLON; ga_char_type['\''] = GA_QUOTE;
42  ga_char_type['.'] = GA_DOT; ga_char_type['@'] = GA_TMULT;
43  ga_char_type[','] = GA_COMMA; ga_char_type[';'] = GA_SEMICOLON;
44  ga_char_type['('] = GA_LPAR; ga_char_type[')'] = GA_RPAR;
45  ga_char_type['['] = GA_LBRACKET; ga_char_type[']'] = GA_RBRACKET;
46  ga_char_type['_'] = GA_NAME; ga_char_type['='] = GA_COLON_EQ;
47  for (unsigned i = 'a'; i <= 'z'; ++i) ga_char_type[i] = GA_NAME;
48  for (unsigned i = 'A'; i <= 'Z'; ++i) ga_char_type[i] = GA_NAME;
49  for (unsigned i = '0'; i <= '9'; ++i) ga_char_type[i] = GA_SCALAR;
50 
51  for (unsigned i=0; i < GA_NB_TOKEN_TYPE; ++i) ga_operator_priorities[i] = 0;
52  ga_operator_priorities[GA_PLUS] = 1;
53  ga_operator_priorities[GA_MINUS] = 1;
54  ga_operator_priorities[GA_MULT] = 2;
55  ga_operator_priorities[GA_DIV] = 2;
56  ga_operator_priorities[GA_COLON] = 2;
57  ga_operator_priorities[GA_DOT] = 2;
58  ga_operator_priorities[GA_DOTMULT] = 2;
59  ga_operator_priorities[GA_DOTDIV] = 2;
60  ga_operator_priorities[GA_TMULT] = 2;
61  ga_operator_priorities[GA_QUOTE] = 3;
62  ga_operator_priorities[GA_UNARY_MINUS] = 3;
63  ga_operator_priorities[GA_SYM] = 4;
64  ga_operator_priorities[GA_SKEW] = 4;
65  ga_operator_priorities[GA_TRACE] = 4;
66  ga_operator_priorities[GA_DEVIATOR] = 4;
67  ga_operator_priorities[GA_PRINT] = 4;
68 
69  return true;
70  }
71 
72  static bool ga_initialized = init_ga_char_type();
73 
74  // Get the next token in the string at position 'pos' end return its type
75  static GA_TOKEN_TYPE ga_get_token(const std::string &expr,
76  size_type &pos,
77  size_type &token_pos,
78  size_type &token_length) {
79  bool fdot = false, fE = false;
80  GMM_ASSERT1(ga_initialized, "Internal error");
81 
82  // Ignore white spaces
83  while (expr[pos] == ' ' && pos < expr.size()) ++pos;
84  token_pos = pos;
85  token_length = 0;
86 
87  // Detecting end of expression
88  if (pos >= expr.size()) return GA_END;
89 
90  // Treating the different cases (Operation, name or number)
91  GA_TOKEN_TYPE type = ga_char_type[unsigned(expr[pos])];
92  ++pos; ++token_length;
93  if (type == GA_DOT) {
94  if (pos >= expr.size()) return type;
95  if (expr[pos] == '*') { ++pos; ++token_length; return GA_DOTMULT; }
96  if (expr[pos] == '/') { ++pos; ++token_length; return GA_DOTDIV; }
97  if (ga_char_type[unsigned(expr[pos])] != GA_SCALAR) return type;
98  fdot = true; type = GA_SCALAR;
99  }
100  switch (type) {
101  case GA_SCALAR:
102  while (pos < expr.size()) {
103  GA_TOKEN_TYPE ctype = ga_char_type[unsigned(expr[pos])];
104  switch (ctype) {
105  case GA_DOT:
106  if (fdot) return type;
107  fdot = true; ++pos; ++token_length;
108  break;
109  case GA_NAME:
110  if (fE || (expr[pos] != 'E' && expr[pos] != 'e')) return type;
111  fE = true; fdot = true; ++pos; ++token_length;
112  if (pos < expr.size()) {
113  if (expr[pos] == '+' || expr[pos] == '-')
114  { ++pos; ++token_length; }
115  }
116  if (pos >= expr.size()
117  || ga_char_type[unsigned(expr[pos])] != GA_SCALAR)
118  return GA_INVALID;
119  break;
120  case GA_SCALAR:
121  ++pos; ++token_length; break;
122  default:
123  return type;
124  }
125  }
126  return type;
127  case GA_NAME:
128  while (pos < expr.size()) {
129  GA_TOKEN_TYPE ctype = ga_char_type[unsigned(expr[pos])];
130  if (ctype != GA_SCALAR && ctype != GA_NAME) break;
131  ++pos; ++token_length;
132  }
133  if (expr.compare(token_pos, token_length, "Sym") == 0)
134  return GA_SYM;
135  if (expr.compare(token_pos, token_length, "Def") == 0)
136  return GA_DEF;
137  if (expr.compare(token_pos, token_length, "Skew") == 0)
138  return GA_SKEW;
139  if (expr.compare(token_pos, token_length, "Trace") == 0)
140  return GA_TRACE;
141  if (expr.compare(token_pos, token_length, "Deviator") == 0)
142  return GA_DEVIATOR;
143  if (expr.compare(token_pos, token_length, "Interpolate") == 0)
144  return GA_INTERPOLATE;
145  if (expr.compare(token_pos, token_length, "Interpolate_derivative") == 0)
146  return GA_INTERPOLATE_DERIVATIVE;
147  if (expr.compare(token_pos, token_length, "Interpolate_filter") == 0)
148  return GA_INTERPOLATE_FILTER;
149  if (expr.compare(token_pos, token_length,
150  "Elementary_transformation") == 0)
151  return GA_ELEMENTARY;
152  if (expr.compare(token_pos, token_length, "Secondary_domain") == 0 ||
153  expr.compare(token_pos, token_length, "Secondary_Domain") == 0)
154  return GA_SECONDARY_DOMAIN;
155  if (expr.compare(token_pos, token_length, "Xfem_plus") == 0)
156  return GA_XFEM_PLUS;
157  if (expr.compare(token_pos, token_length, "Xfem_minus") == 0)
158  return GA_XFEM_MINUS;
159  if (expr.compare(token_pos, token_length, "Print") == 0)
160  return GA_PRINT;
161  return type;
162  case GA_COMMA:
163  if (pos < expr.size() &&
164  ga_char_type[unsigned(expr[pos])] == GA_COMMA) {
165  ++pos; return GA_DCOMMA;
166  }
167  return type;
168  case GA_SEMICOLON:
169  if (pos < expr.size() &&
170  ga_char_type[unsigned(expr[pos])] == GA_SEMICOLON) {
171  ++pos; return GA_DSEMICOLON;
172  }
173  return type;
174  case GA_COLON:
175  if (pos < expr.size() &&
176  ga_char_type[unsigned(expr[pos])] == GA_COLON_EQ) {
177  ++pos; return GA_COLON_EQ;
178  }
179  return type;
180  case GA_COLON_EQ:
181  return GA_INVALID;
182  default: return type;
183  }
184  }
185 
186  //=========================================================================
187  // Error handling
188  //=========================================================================
189 
190  void ga_throw_error_msg(pstring expr, size_type pos,
191  const std::string &msg) {
192  int length_before = 70, length_after = 70;
193  if (expr && expr->size()) {
194  int first = std::max(0, int(pos)-length_before);
195  int last = std::min(int(pos)+length_after, int(expr->size()));
196  if (last - first < length_before+length_after)
197  first = std::max(0, int(pos)-length_before
198  -(length_before+length_after-last+first));
199  if (last - first < length_before+length_after)
200  last = std::min(int(pos)+length_after
201  +(length_before+length_after-last+first),
202  int(expr->size()));
203  if (first > 0) cerr << "...";
204  cerr << expr->substr(first, last-first);
205  if (last < int(expr->size())) cerr << "...";
206  cerr << endl;
207  if (first > 0) cerr << " ";
208  if (int(pos) > first)
209  cerr << std::setfill ('-') << std::setw(int(pos)-first) << '-'
210  << std::setfill (' ');
211  cerr << "^" << endl;
212  }
213  cerr << msg << endl;
214  }
215 
216  //=========================================================================
217  // Tree structure
218  //=========================================================================
219 
220  void ga_tree_node::mult_test(const pga_tree_node n0, const pga_tree_node n1) {
221 
222  size_type test0 = n0->test_function_type, test1 = n1->test_function_type;
223  if (test0 && test1 && (test0 == test1 || test0 >= 3 || test1 >= 3))
224  ga_throw_error(expr, pos,
225  "Incompatibility of test functions in product.");
226  GMM_ASSERT1(test0 != size_type(-1) && test1 != size_type(-1),
227  "internal error");
228 
229  test_function_type = test0 + test1;
230 
231  size_type st = nb_test_functions();
232  bgeot::multi_index mi(st);
233 
234  switch (test0) {
235  case 1: mi[0] = n0->t.sizes()[0]; break;
236  case 2: mi[st-1] = n0->t.sizes()[0]; break;
237  case 3: mi[0] = n0->t.sizes()[0]; mi[1] = n0->t.sizes()[1]; break;
238  }
239  switch (test1) {
240  case 1: mi[0] = n1->t.sizes()[0]; break;
241  case 2: mi[st-1] = n1->t.sizes()[0]; break;
242  case 3: mi[0] = n1->t.sizes()[0]; mi[1] = n1->t.sizes()[1]; break;
243  }
244 
245  if (n0->name_test1.size()) {
246  name_test1 = n0->name_test1; qdim1 = n0->qdim1;
247  interpolate_name_test1 = n0->interpolate_name_test1;
248  } else {
249  name_test1 = n1->name_test1; qdim1 = n1->qdim1;
250  interpolate_name_test1 = n1->interpolate_name_test1;
251  }
252 
253  if (n0->name_test2.size()) {
254  name_test2 = n0->name_test2; qdim2 = n0->qdim2;
255  interpolate_name_test2 = n0->interpolate_name_test2;
256  } else {
257  name_test2 = n1->name_test2; qdim2 = n1->qdim2;
258  interpolate_name_test2 = n1->interpolate_name_test2;
259  }
260  t.adjust_sizes(mi);
261  }
262 
263  void ga_tree::add_scalar(scalar_type val, size_type pos, pstring expr) {
264  while (current_node && current_node->node_type != GA_NODE_OP)
265  current_node = current_node->parent;
266  if (current_node) {
267  current_node->adopt_child(new ga_tree_node(val, pos, expr));
268  current_node = current_node->children.back();
269  }
270  else {
271  GMM_ASSERT1(root == nullptr, "Invalid tree operation");
272  current_node = root = new ga_tree_node(val, pos, expr);
273  root->parent = nullptr;
274  }
275  }
276 
277  void ga_tree::add_allindices(size_type pos, pstring expr) {
278  while (current_node && current_node->node_type != GA_NODE_OP)
279  current_node = current_node->parent;
280  if (current_node) {
281  current_node->adopt_child(new ga_tree_node(GA_NODE_ALLINDICES, pos,expr));
282  current_node = current_node->children.back();
283  }
284  else {
285  GMM_ASSERT1(root == nullptr, "Invalid tree operation");
286  current_node = root = new ga_tree_node(GA_NODE_ALLINDICES, pos, expr);
287  root->parent = nullptr;
288  }
289  }
290 
291  void ga_tree::add_name(const char *name, size_type length, size_type pos,
292  pstring expr) {
293  while (current_node && current_node->node_type != GA_NODE_OP)
294  current_node = current_node->parent;
295  if (current_node) {
296  current_node->adopt_child(new ga_tree_node(name, length, pos, expr));
297  current_node = current_node->children.back();
298  }
299  else {
300  GMM_ASSERT1(root == nullptr, "Invalid tree operation");
301  current_node = root = new ga_tree_node(name, length, pos, expr);
302  root->parent = nullptr;
303  }
304  }
305 
306  void ga_tree::add_sub_tree(ga_tree &sub_tree) {
307  if (current_node &&
308  (current_node->node_type == GA_NODE_PARAMS ||
309  current_node->node_type == GA_NODE_INTERPOLATE_FILTER ||
310  current_node->node_type == GA_NODE_C_MATRIX)) {
311  GMM_ASSERT1(sub_tree.root, "Invalid tree operation");
312  current_node->adopt_child(sub_tree.root);
313  } else {
314  GMM_ASSERT1(sub_tree.root, "Invalid tree operation");
315  while (current_node && current_node->node_type != GA_NODE_OP)
316  current_node = current_node->parent;
317  if (current_node) {
318  current_node->adopt_child(sub_tree.root);
319  current_node = sub_tree.root;
320  }
321  else {
322  GMM_ASSERT1(root == nullptr, "Invalid tree operation");
323  current_node = root = sub_tree.root;
324  root->parent = nullptr;
325  }
326  }
327  sub_tree.root = sub_tree.current_node = nullptr;
328  }
329 
330  void ga_tree::add_params(size_type pos, pstring expr) {
331  GMM_ASSERT1(current_node, "internal error");
332  while (current_node && current_node->parent &&
333  current_node->parent->node_type == GA_NODE_OP &&
334  ga_operator_priorities[current_node->parent->op_type] >= 4)
335  current_node = current_node->parent;
336  pga_tree_node new_node = new ga_tree_node(GA_NODE_PARAMS, pos, expr);
337  new_node->parent = current_node->parent;
338  if (current_node->parent)
339  current_node->parent->replace_child(current_node, new_node);
340  else
341  root = new_node;
342  new_node->adopt_child(current_node);
343  current_node = new_node;
344  }
345 
346  void ga_tree::add_matrix(size_type pos, pstring expr) {
347  while (current_node && current_node->node_type != GA_NODE_OP)
348  current_node = current_node->parent;
349  if (current_node) {
350  current_node->adopt_child(new ga_tree_node(GA_NODE_C_MATRIX, pos, expr));
351  current_node = current_node->children.back();
352  }
353  else {
354  GMM_ASSERT1(root == nullptr, "Invalid tree operation");
355  current_node = root = new ga_tree_node(GA_NODE_C_MATRIX, pos, expr);
356  root->parent = nullptr;
357  }
358  current_node->nbc1 = current_node->nbc2 = current_node->nbc3 = 0;
359  }
360 
361  void ga_tree::add_op(GA_TOKEN_TYPE op_type, size_type pos,
362  pstring expr) {
363  while (current_node && current_node->parent &&
364  current_node->parent->node_type == GA_NODE_OP &&
365  ga_operator_priorities[current_node->parent->op_type]
366  >= ga_operator_priorities[op_type])
367  current_node = current_node->parent;
368  pga_tree_node new_node = new ga_tree_node(op_type, pos, expr);
369  if (current_node) {
370  if (op_type == GA_UNARY_MINUS
371  || op_type == GA_SYM || op_type == GA_SKEW
372  || op_type == GA_TRACE || op_type == GA_DEVIATOR
373  || op_type == GA_PRINT) {
374  current_node->adopt_child(new_node);
375  } else {
376  new_node->parent = current_node->parent;
377  if (current_node->parent)
378  current_node->parent->replace_child(current_node, new_node);
379  else
380  root = new_node;
381  new_node->adopt_child(current_node);
382  }
383  } else {
384  if (root) new_node->adopt_child(root);
385  root = new_node;
386  root->parent = nullptr;
387  }
388  current_node = new_node;
389  }
390 
391  void ga_tree::clear_node_rec(pga_tree_node pnode) {
392  if (pnode) {
393  for (pga_tree_node &child : pnode->children)
394  clear_node_rec(child);
395  delete pnode;
396  current_node = nullptr;
397  }
398  }
399 
400  void ga_tree::clear_node(pga_tree_node pnode) {
401  if (pnode) {
402  pga_tree_node parent = pnode->parent;
403  if (parent) { // keep all siblings of pnode
404  size_type j = 0;
405  for (pga_tree_node &sibling : parent->children)
406  if (sibling != pnode)
407  parent->children[j++] = sibling;
408  parent->children.resize(j);
409  } else
410  root = nullptr;
411  }
412  clear_node_rec(pnode);
413  }
414 
415  void ga_tree::clear_children(pga_tree_node pnode) {
416  for (pga_tree_node &child : pnode->children)
417  clear_node_rec(child);
418  pnode->children.resize(0);
419  }
420 
421  void ga_tree::replace_node_by_child(pga_tree_node pnode, size_type i) {
422  GMM_ASSERT1(i < pnode->children.size(), "Internal error");
423  pga_tree_node child = pnode->children[i];
424  child->parent = pnode->parent;
425  if (pnode->parent)
426  pnode->parent->replace_child(pnode, child);
427  else
428  root = child;
429  current_node = nullptr;
430  for (pga_tree_node &sibling : pnode->children)
431  if (sibling != child) clear_node_rec(sibling);
432  delete pnode;
433  }
434 
435  void ga_tree::copy_node(pga_tree_node pnode, pga_tree_node parent,
436  pga_tree_node &child) {
437  GMM_ASSERT1(child == nullptr, "Internal error");
438  child = new ga_tree_node();
439  *child = *pnode;
440  child->parent = parent;
441  for (pga_tree_node &grandchild : child->children)
442  grandchild = nullptr;
443  for (size_type j = 0; j < child->children.size(); ++j)
444  copy_node(pnode->children[j], child, child->children[j]);
445  }
446 
447  void ga_tree::duplicate_with_operation(pga_tree_node pnode,
448  GA_TOKEN_TYPE op_type) {
449  pga_tree_node newop = new ga_tree_node(op_type, pnode->pos, pnode->expr);
450  newop->children.resize(2, nullptr);
451  newop->children[0] = pnode;
452  newop->pos = pnode->pos; newop->expr = pnode->expr;
453  newop->parent = pnode->parent;
454  if (pnode->parent)
455  pnode->parent->replace_child(pnode, newop);
456  else
457  root = newop;
458  pnode->parent = newop;
459  copy_node(pnode, newop, newop->children[1]);
460  }
461 
462  void ga_tree::add_child(pga_tree_node pnode, GA_NODE_TYPE node_type) {
463  pga_tree_node newnode=new ga_tree_node();
464  newnode->pos = pnode->pos; newnode->expr = pnode->expr;
465  newnode->node_type = node_type; pnode->adopt_child(newnode);
466  }
467 
468  void ga_tree::insert_node(pga_tree_node pnode, GA_NODE_TYPE node_type) {
469  pga_tree_node newnode = new ga_tree_node();
470  newnode->node_type = node_type;
471  newnode->parent = pnode->parent;
472  newnode->pos = pnode->pos; newnode->expr = pnode->expr;
473  if (pnode->parent)
474  pnode->parent->replace_child(pnode, newnode);
475  else
476  root = newnode;
477  newnode->adopt_child(pnode);
478  }
479 
480  bool sub_tree_are_equal
481  (const pga_tree_node pnode1, const pga_tree_node pnode2,
482  const ga_workspace &workspace, int version) {
483 
484  size_type ntype1 = pnode1->node_type;
485  if (ntype1 == GA_NODE_ZERO) ntype1 = GA_NODE_CONSTANT;
486  size_type ntype2 = pnode2->node_type;
487  if (ntype2 == GA_NODE_ZERO) ntype2 = GA_NODE_CONSTANT;
488  if (ntype1 != ntype2) return false;
489  if (pnode1->children.size() != pnode2->children.size()) return false;
490 
491  switch(ntype1) {
492  case GA_NODE_OP:
493  if (pnode1->op_type != pnode2->op_type) return false;
494  if (pnode1->symmetric_op != pnode2->symmetric_op) return false;
495  break;
496  case GA_NODE_OPERATOR:
497  if (pnode1->der1 != pnode2->der1 || pnode1->der2 != pnode2->der2)
498  return false;
499  if (pnode1->name.compare(pnode2->name)) return false;
500  break;
501  case GA_NODE_PREDEF_FUNC: case GA_NODE_SPEC_FUNC:
502  if (pnode1->name.compare(pnode2->name)) return false;
503  break;
504  case GA_NODE_CONSTANT: case GA_NODE_ZERO:
505  if (pnode1->tensor().size() != pnode2->tensor().size()) return false;
506 
507  switch(version) {
508  case 0: case 1:
509  if (pnode1->test_function_type != pnode2->test_function_type)
510  return false;
511  if ((pnode1->test_function_type & 1) &&
512  pnode1->name_test1.compare(pnode2->name_test1) != 0)
513  return false;
514  if ((pnode1->test_function_type & 2) &&
515  pnode1->name_test2.compare(pnode2->name_test2) != 0)
516  return false;
517  break;
518  case 2:
519  if ((pnode1->test_function_type == 1 &&
520  pnode2->test_function_type == 1) ||
521  (pnode1->test_function_type == 2 &&
522  pnode2->test_function_type == 2))
523  return false;
524  if ((pnode1->test_function_type & 1) &&
525  pnode1->name_test1.compare(pnode2->name_test2) != 0)
526  return false;
527  if ((pnode1->test_function_type & 2) &&
528  pnode1->name_test2.compare(pnode2->name_test1) != 0)
529  return false;
530  break;
531  }
532  if (pnode1->tensor().size() != 1 &&
533  pnode1->t.sizes().size() != pnode2->t.sizes().size()) return false;
534  for (size_type i = 0; i < pnode1->t.sizes().size(); ++i)
535  if (pnode1->t.sizes()[i] != pnode2->t.sizes()[i]) return false;
536  for (size_type i = 0; i < pnode1->tensor().size(); ++i)
537  if (gmm::abs(pnode1->tensor()[i] - pnode2->tensor()[i]) > 1E-25)
538  return false;
539  break;
540  case GA_NODE_C_MATRIX:
541  if (pnode1->children.size() != pnode2->children.size()) return false;
542  if (pnode1->nb_test_functions() != pnode2->nb_test_functions())
543  return false;
544  if (pnode1->t.sizes().size() != pnode2->t.sizes().size()) return false;
545  for (size_type i=pnode1->nb_test_functions();
546  i<pnode1->t.sizes().size(); ++i)
547  if (pnode1->t.sizes()[i] != pnode2->t.sizes()[i]) return false;
548  if (pnode1->nbc1 != pnode2->nbc1) return false;
549  break;
550  case GA_NODE_INTERPOLATE_FILTER:
551  if (pnode1->interpolate_name.compare(pnode2->interpolate_name) ||
552  pnode1->nbc1 != pnode2->nbc1)
553  return false;
554  break;
555  case GA_NODE_INTERPOLATE_X: case GA_NODE_SECONDARY_DOMAIN_X:
556  case GA_NODE_INTERPOLATE_NORMAL: case GA_NODE_SECONDARY_DOMAIN_NORMAL:
557  if (pnode1->interpolate_name.compare(pnode2->interpolate_name))
558  return false;
559  break;
560  case GA_NODE_INTERPOLATE_DERIVATIVE:
561  if (pnode1->interpolate_name_der.compare(pnode2->interpolate_name_der))
562  return false;
563  if (pnode1->interpolate_name.compare(pnode2->interpolate_name) ||
564  pnode1->elementary_name.compare(pnode2->elementary_name))
565  return false;
566  {
567  const mesh_fem *mf1 = workspace.associated_mf(pnode1->name);
568  const mesh_fem *mf2 = workspace.associated_mf(pnode2->name);
569  switch (version) {
570  case 0:
571  if (pnode1->name.compare(pnode2->name) ||
572  pnode1->test_function_type != pnode2->test_function_type)
573  return false;
574  break;
575  case 1:
576  if (mf1 != mf2 ||
577  workspace.qdim(pnode1->name) != workspace.qdim(pnode2->name) ||
578  pnode1->test_function_type != pnode2->test_function_type)
579  return false;
580  break;
581  case 2:
582  if (mf1 != mf2 ||
583  workspace.qdim(pnode1->name) != workspace.qdim(pnode2->name) ||
584  pnode1->test_function_type == pnode2->test_function_type)
585  return false;
586  break;
587  }
588  }
589  break;
590  case GA_NODE_INTERPOLATE_VAL_TEST: case GA_NODE_INTERPOLATE_GRAD_TEST:
591  case GA_NODE_INTERPOLATE_HESS_TEST: case GA_NODE_INTERPOLATE_DIVERG_TEST:
592  case GA_NODE_ELEMENTARY_VAL_TEST: case GA_NODE_ELEMENTARY_GRAD_TEST:
593  case GA_NODE_ELEMENTARY_HESS_TEST: case GA_NODE_ELEMENTARY_DIVERG_TEST:
594  case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
595  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
596  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
597  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
598  case GA_NODE_XFEM_PLUS_VAL_TEST: case GA_NODE_XFEM_PLUS_GRAD_TEST:
599  case GA_NODE_XFEM_PLUS_HESS_TEST: case GA_NODE_XFEM_PLUS_DIVERG_TEST:
600  case GA_NODE_XFEM_MINUS_VAL_TEST: case GA_NODE_XFEM_MINUS_GRAD_TEST:
601  case GA_NODE_XFEM_MINUS_HESS_TEST: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
602  if (pnode1->interpolate_name.compare(pnode2->interpolate_name) ||
603  pnode1->elementary_name.compare(pnode2->elementary_name))
604  return false;
605  {
606  const mesh_fem *mf1 = workspace.associated_mf(pnode1->name);
607  const mesh_fem *mf2 = workspace.associated_mf(pnode2->name);
608  switch (version) {
609  case 0:
610  if (pnode1->name.compare(pnode2->name) ||
611  pnode1->test_function_type != pnode2->test_function_type)
612  return false;
613  break;
614  case 1:
615  if (mf1 != mf2 ||
616  workspace.qdim(pnode1->name) != workspace.qdim(pnode2->name) ||
617  pnode1->test_function_type != pnode2->test_function_type)
618  return false;
619  break;
620  case 2:
621  if (mf1 != mf2 ||
622  workspace.qdim(pnode1->name) != workspace.qdim(pnode2->name) ||
623  pnode1->test_function_type == pnode2->test_function_type)
624  return false;
625  break;
626  }
627  }
628  break;
629  case GA_NODE_VAL_TEST: case GA_NODE_GRAD_TEST:
630  case GA_NODE_HESS_TEST: case GA_NODE_DIVERG_TEST:
631  {
632  const mesh_fem *mf1 = workspace.associated_mf(pnode1->name);
633  const mesh_fem *mf2 = workspace.associated_mf(pnode2->name);
634  switch (version) {
635  case 0:
636  if (pnode1->name.compare(pnode2->name) ||
637  pnode1->test_function_type != pnode2->test_function_type)
638  return false;
639  break;
640  case 1:
641  if (mf1 != mf2 ||
642  workspace.qdim(pnode1->name) != workspace.qdim(pnode2->name) ||
643  pnode1->test_function_type != pnode2->test_function_type)
644  return false;
645  break;
646  case 2:
647  if (mf1 != mf2 ||
648  workspace.qdim(pnode1->name) != workspace.qdim(pnode2->name) ||
649  pnode1->test_function_type == pnode2->test_function_type)
650  return false;
651  break;
652  }
653  }
654  break;
655  case GA_NODE_VAL: case GA_NODE_GRAD:
656  case GA_NODE_HESS: case GA_NODE_DIVERG:
657  if (pnode1->name.compare(pnode2->name)) return false;
658  break;
659  case GA_NODE_INTERPOLATE_VAL: case GA_NODE_INTERPOLATE_GRAD:
660  case GA_NODE_INTERPOLATE_HESS: case GA_NODE_INTERPOLATE_DIVERG:
661  case GA_NODE_ELEMENTARY_VAL: case GA_NODE_ELEMENTARY_GRAD:
662  case GA_NODE_ELEMENTARY_HESS: case GA_NODE_ELEMENTARY_DIVERG:
663  case GA_NODE_SECONDARY_DOMAIN_VAL: case GA_NODE_SECONDARY_DOMAIN_GRAD:
664  case GA_NODE_SECONDARY_DOMAIN_HESS: case GA_NODE_SECONDARY_DOMAIN_DIVERG:
665  case GA_NODE_XFEM_PLUS_VAL: case GA_NODE_XFEM_PLUS_GRAD:
666  case GA_NODE_XFEM_PLUS_HESS: case GA_NODE_XFEM_PLUS_DIVERG:
667  case GA_NODE_XFEM_MINUS_VAL: case GA_NODE_XFEM_MINUS_GRAD:
668  case GA_NODE_XFEM_MINUS_HESS: case GA_NODE_XFEM_MINUS_DIVERG:
669  if (pnode1->interpolate_name.compare(pnode2->interpolate_name) ||
670  pnode1->elementary_name.compare(pnode2->elementary_name) ||
671  pnode1->name.compare(pnode2->name))
672  return false;
673  break;
674  case GA_NODE_X:
675  if (pnode1->nbc1 != pnode2->nbc1) return false;
676  break;
677 
678  default:break;
679  }
680 
681  if (version && ntype1 == GA_NODE_OP && pnode1->symmetric_op) {
682  if (sub_tree_are_equal(pnode1->children[0], pnode2->children[0],
683  workspace, version) &&
684  sub_tree_are_equal(pnode1->children[1], pnode2->children[1],
685  workspace, version))
686  return true;
687  if (sub_tree_are_equal(pnode1->children[1], pnode2->children[0],
688  workspace, version) &&
689  sub_tree_are_equal(pnode1->children[0], pnode2->children[1],
690  workspace, version) )
691  return true;
692  return false;
693  } else {
694  for (size_type i = 0; i < pnode1->children.size(); ++i)
695  if (!(sub_tree_are_equal(pnode1->children[i], pnode2->children[i],
696  workspace, version)))
697  return false;
698  }
699  return true;
700  }
701 
702  static void verify_tree(const pga_tree_node pnode,
703  const pga_tree_node parent) {
704  GMM_ASSERT1(pnode->parent == parent,
705  "Invalid tree node " << pnode->node_type);
706  for (pga_tree_node &child : pnode->children)
707  verify_tree(child, pnode);
708  }
709 
710 
711  static void ga_print_constant_tensor(const pga_tree_node pnode,
712  std::ostream &str) {
713  size_type nt = pnode->nb_test_functions(); // for printing zero tensors
714  switch (pnode->tensor_order()) {
715  case 0:
716  str << (nt ? scalar_type(0) : pnode->tensor()[0]);
717  break;
718 
719  case 1:
720  str << "[";
721  for (size_type i = 0; i < pnode->tensor_proper_size(0); ++i) {
722  if (i != 0) str << ",";
723  str << (nt ? scalar_type(0) : pnode->tensor()[i]);
724  }
725  str << "]";
726  break;
727 
728  case 2: case 3: case 4:
729  {
730  size_type ii(0);
731  size_type n0 = pnode->tensor_proper_size(0);
732  size_type n1 = pnode->tensor_proper_size(1);
733  size_type n2 = ((pnode->tensor_order() > 2) ?
734  pnode->tensor_proper_size(2) : 1);
735  size_type n3 = ((pnode->tensor_order() > 3) ?
736  pnode->tensor_proper_size(3) : 1);
737  if (n3 > 1) str << "[";
738  for (size_type l = 0; l < n3; ++l) {
739  if (l != 0) str << ",";
740  if (n2 > 1) str << "[";
741  for (size_type k = 0; k < n2; ++k) {
742  if (k != 0) str << ",";
743  if (n1 > 1) str << "[";
744  for (size_type j = 0; j < n1; ++j) {
745  if (j != 0) str << ",";
746  if (n0 > 1) str << "[";
747  for (size_type i = 0; i < n0; ++i) {
748  if (i != 0) str << ",";
749  str << (nt ? scalar_type(0) : pnode->tensor()[ii++]);
750  }
751  if (n0 > 1) str << "]";
752  }
753  if (n1 > 1) str << "]";
754  }
755  if (n2 > 1) str << "]";
756  }
757  if (n3 > 1) str << "]";
758  }
759  break;
760 
761  case 5: case 6: case 7: case 8:
762  str << "Reshape([";
763  for (size_type i = 0; i < pnode->tensor_proper_size(); ++i) {
764  if (i != 0) str << ";";
765  str << (nt ? scalar_type(0) : pnode->tensor()[i]);
766  }
767  str << "]";
768  for (size_type i = 0; i < pnode->tensor_order(); ++i) {
769  if (i != 0) str << ", ";
770  str << pnode->tensor_proper_size(i);
771  }
772  str << ")";
773  break;
774 
775  default: GMM_ASSERT1(false, "Invalid tensor dimension");
776  }
777  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
778  }
779 
780  void ga_print_node(const pga_tree_node pnode,
781  std::ostream &str) {
782  if (!pnode) return;
783  long prec = str.precision(16);
784 
785  bool is_interpolate(false), is_elementary(false), is_secondary(false);
786  bool is_xfem_plus(false), is_xfem_minus(false);
787  switch(pnode->node_type) {
788  case GA_NODE_INTERPOLATE:
789  case GA_NODE_INTERPOLATE_X:
790  case GA_NODE_INTERPOLATE_NORMAL:
791  case GA_NODE_INTERPOLATE_VAL:
792  case GA_NODE_INTERPOLATE_GRAD:
793  case GA_NODE_INTERPOLATE_HESS:
794  case GA_NODE_INTERPOLATE_DIVERG:
795  case GA_NODE_INTERPOLATE_VAL_TEST:
796  case GA_NODE_INTERPOLATE_GRAD_TEST:
797  case GA_NODE_INTERPOLATE_HESS_TEST:
798  case GA_NODE_INTERPOLATE_DIVERG_TEST:
799  str << "Interpolate(";
800  is_interpolate = true;
801  break;
802  case GA_NODE_ELEMENTARY:
803  case GA_NODE_ELEMENTARY_VAL:
804  case GA_NODE_ELEMENTARY_GRAD:
805  case GA_NODE_ELEMENTARY_HESS:
806  case GA_NODE_ELEMENTARY_DIVERG:
807  case GA_NODE_ELEMENTARY_VAL_TEST:
808  case GA_NODE_ELEMENTARY_GRAD_TEST:
809  case GA_NODE_ELEMENTARY_HESS_TEST:
810  case GA_NODE_ELEMENTARY_DIVERG_TEST:
811  is_elementary = true;
812  str << "Elementary_transformation(";
813  break;
814  case GA_NODE_SECONDARY_DOMAIN:
815  case GA_NODE_SECONDARY_DOMAIN_X:
816  case GA_NODE_SECONDARY_DOMAIN_NORMAL:
817  case GA_NODE_SECONDARY_DOMAIN_VAL:
818  case GA_NODE_SECONDARY_DOMAIN_GRAD:
819  case GA_NODE_SECONDARY_DOMAIN_HESS:
820  case GA_NODE_SECONDARY_DOMAIN_DIVERG:
821  case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
822  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
823  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
824  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
825  str << "Secondary_domain(";
826  is_secondary = true;
827  break;
828 
829  case GA_NODE_XFEM_PLUS:
830  case GA_NODE_XFEM_PLUS_VAL:
831  case GA_NODE_XFEM_PLUS_GRAD:
832  case GA_NODE_XFEM_PLUS_HESS:
833  case GA_NODE_XFEM_PLUS_DIVERG:
834  case GA_NODE_XFEM_PLUS_VAL_TEST:
835  case GA_NODE_XFEM_PLUS_GRAD_TEST:
836  case GA_NODE_XFEM_PLUS_HESS_TEST:
837  case GA_NODE_XFEM_PLUS_DIVERG_TEST:
838  is_xfem_plus = true;
839  str << "Xfem_plus(";
840  break;
841  case GA_NODE_XFEM_MINUS:
842  case GA_NODE_XFEM_MINUS_VAL:
843  case GA_NODE_XFEM_MINUS_GRAD:
844  case GA_NODE_XFEM_MINUS_HESS:
845  case GA_NODE_XFEM_MINUS_DIVERG:
846  case GA_NODE_XFEM_MINUS_VAL_TEST:
847  case GA_NODE_XFEM_MINUS_GRAD_TEST:
848  case GA_NODE_XFEM_MINUS_HESS_TEST:
849  case GA_NODE_XFEM_MINUS_DIVERG_TEST:
850  is_xfem_minus = true;
851  str << "Xfem_minus(";
852  break;
853  default:
854  break;
855  }
856 
857  switch(pnode->node_type) {
858  case GA_NODE_GRAD:
859  case GA_NODE_INTERPOLATE_GRAD:
860  case GA_NODE_ELEMENTARY_GRAD:
861  case GA_NODE_SECONDARY_DOMAIN_GRAD:
862  case GA_NODE_XFEM_PLUS_GRAD:
863  case GA_NODE_XFEM_MINUS_GRAD:
864  case GA_NODE_GRAD_TEST:
865  case GA_NODE_INTERPOLATE_GRAD_TEST:
866  case GA_NODE_ELEMENTARY_GRAD_TEST:
867  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
868  case GA_NODE_XFEM_PLUS_GRAD_TEST:
869  case GA_NODE_XFEM_MINUS_GRAD_TEST:
870  str << "Grad_";
871  break;
872  case GA_NODE_HESS:
873  case GA_NODE_INTERPOLATE_HESS:
874  case GA_NODE_ELEMENTARY_HESS:
875  case GA_NODE_SECONDARY_DOMAIN_HESS:
876  case GA_NODE_XFEM_PLUS_HESS:
877  case GA_NODE_XFEM_MINUS_HESS:
878  case GA_NODE_HESS_TEST:
879  case GA_NODE_INTERPOLATE_HESS_TEST:
880  case GA_NODE_ELEMENTARY_HESS_TEST:
881  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
882  case GA_NODE_XFEM_PLUS_HESS_TEST:
883  case GA_NODE_XFEM_MINUS_HESS_TEST:
884  str << "Hess_";
885  break;
886  case GA_NODE_DIVERG:
887  case GA_NODE_INTERPOLATE_DIVERG:
888  case GA_NODE_SECONDARY_DOMAIN_DIVERG:
889  case GA_NODE_ELEMENTARY_DIVERG:
890  case GA_NODE_XFEM_PLUS_DIVERG:
891  case GA_NODE_XFEM_MINUS_DIVERG:
892  case GA_NODE_DIVERG_TEST:
893  case GA_NODE_INTERPOLATE_DIVERG_TEST:
894  case GA_NODE_ELEMENTARY_DIVERG_TEST:
895  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
896  case GA_NODE_XFEM_PLUS_DIVERG_TEST:
897  case GA_NODE_XFEM_MINUS_DIVERG_TEST:
898  str << "Div_";
899  break;
900  default:
901  break;
902  }
903 
904  switch(pnode->node_type) {
905  case GA_NODE_OP:
906  {
907  bool par = false;
908  if (pnode->parent) {
909  if (pnode->parent->node_type == GA_NODE_OP &&
910  (ga_operator_priorities[pnode->op_type] >= 2 ||
911  ga_operator_priorities[pnode->op_type]
912  < ga_operator_priorities[pnode->parent->op_type]))
913  par = true;
914  if (pnode->parent->node_type == GA_NODE_PARAMS) par = true;
915  }
916 
917 
918  if (par) str << "(";
919  if (pnode->op_type == GA_UNARY_MINUS) {
920  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
921  str << "-"; ga_print_node(pnode->children[0], str);
922  } else if (pnode->op_type == GA_QUOTE) {
923  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
924  ga_print_node(pnode->children[0], str); str << "'";
925  } else if (pnode->op_type == GA_SYM) {
926  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
927  str << "Sym("; ga_print_node(pnode->children[0], str); str << ")";
928  } else if (pnode->op_type == GA_SKEW) {
929  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
930  str << "Skew("; ga_print_node(pnode->children[0], str); str << ")";
931  } else if (pnode->op_type == GA_TRACE) {
932  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
933  str << "Trace("; ga_print_node(pnode->children[0], str); str << ")";
934  } else if (pnode->op_type == GA_DEVIATOR) {
935  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree with "
936  << pnode->children.size() << " children instead of 1");
937  str << "Deviator("; ga_print_node(pnode->children[0], str); str<<")";
938  } else if (pnode->op_type == GA_PRINT) {
939  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
940  str << "Print("; ga_print_node(pnode->children[0], str); str << ")";
941  } else {
942  if (!par && pnode->op_type == GA_MULT &&
943  (pnode->children.size() == 1 ||
944  pnode->test_function_type == size_type(-1) ||
945  (pnode->children[0]->tensor_order() == 4 &&
946  pnode->children[1]->tensor_order() == 2)))
947  { par = true; str << "("; }
948  ga_print_node(pnode->children[0], str);
949  switch (pnode->op_type) {
950  case GA_PLUS: str << "+"; break;
951  case GA_MINUS: str << "-"; break;
952  case GA_MULT: str << "*"; break;
953  case GA_DIV: str << "/"; break;
954  case GA_COLON: str << ":"; break;
955  case GA_DOT: str << "."; break;
956  case GA_DOTMULT: str << ".*"; break;
957  case GA_DOTDIV: str << "./"; break;
958  case GA_TMULT: str << "@"; break;
959  default: GMM_ASSERT1(false, "Invalid or not taken into account "
960  "operation");
961  }
962  if (pnode->children.size() >= 2)
963  ga_print_node(pnode->children[1], str);
964  else
965  str << "(unknown second argument)";
966  }
967  if (par) str << ")";
968  }
969  break;
970 
971  case GA_NODE_X:
972  if (pnode->nbc1) str << "X(" << pnode->nbc1 << ")"; else str << "X";
973  break;
974  case GA_NODE_ELT_SIZE: str << "element_size"; break;
975  case GA_NODE_ELT_K: str << "element_K"; break;
976  case GA_NODE_ELT_B: str << "element_B"; break;
977  case GA_NODE_NORMAL: str << "Normal"; break;
978  case GA_NODE_INTERPOLATE_FILTER:
979  str << "Interpolate_filter(" << pnode->interpolate_name << ",";
980  ga_print_node(pnode->children[0], str);
981  if (pnode->children.size() == 2)
982  { str << ","; ga_print_node(pnode->children[1], str); }
983  else if (pnode->nbc1 != size_type(-1)) str << "," << pnode->nbc1;
984  str << ")";
985  break;
986  case GA_NODE_INTERPOLATE_X: case GA_NODE_SECONDARY_DOMAIN_X:
987  str << "X";
988  break;
989  case GA_NODE_INTERPOLATE_NORMAL: case GA_NODE_SECONDARY_DOMAIN_NORMAL:
990  str << "Normal";
991  break;
992  case GA_NODE_INTERPOLATE_DERIVATIVE:
993  str << (pnode->test_function_type == 1 ? "Test_" : "Test2_")
994  << "Interpolate_derivative(" << pnode->interpolate_name_der << ","
995  << pnode->name;
996  if (pnode->interpolate_name.size())
997  str << "," << pnode->interpolate_name;
998  str << ")";
999  break;
1000  case GA_NODE_INTERPOLATE:
1001  case GA_NODE_ELEMENTARY:
1002  case GA_NODE_SECONDARY_DOMAIN:
1003  case GA_NODE_XFEM_PLUS:
1004  case GA_NODE_XFEM_MINUS:
1005  case GA_NODE_VAL:
1006  case GA_NODE_INTERPOLATE_VAL:
1007  case GA_NODE_ELEMENTARY_VAL:
1008  case GA_NODE_SECONDARY_DOMAIN_VAL:
1009  case GA_NODE_XFEM_PLUS_VAL:
1010  case GA_NODE_XFEM_MINUS_VAL:
1011  case GA_NODE_GRAD:
1012  case GA_NODE_INTERPOLATE_GRAD:
1013  case GA_NODE_SECONDARY_DOMAIN_GRAD:
1014  case GA_NODE_ELEMENTARY_GRAD:
1015  case GA_NODE_XFEM_PLUS_GRAD:
1016  case GA_NODE_XFEM_MINUS_GRAD:
1017  case GA_NODE_HESS:
1018  case GA_NODE_INTERPOLATE_HESS:
1019  case GA_NODE_SECONDARY_DOMAIN_HESS:
1020  case GA_NODE_ELEMENTARY_HESS:
1021  case GA_NODE_XFEM_PLUS_HESS:
1022  case GA_NODE_XFEM_MINUS_HESS:
1023  case GA_NODE_DIVERG:
1024  case GA_NODE_INTERPOLATE_DIVERG:
1025  case GA_NODE_ELEMENTARY_DIVERG:
1026  case GA_NODE_SECONDARY_DOMAIN_DIVERG:
1027  case GA_NODE_XFEM_PLUS_DIVERG:
1028  case GA_NODE_XFEM_MINUS_DIVERG:
1029  str << pnode->name;
1030  break;
1031  case GA_NODE_VAL_TEST:
1032  case GA_NODE_INTERPOLATE_VAL_TEST:
1033  case GA_NODE_ELEMENTARY_VAL_TEST:
1034  case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
1035  case GA_NODE_XFEM_PLUS_VAL_TEST:
1036  case GA_NODE_XFEM_MINUS_VAL_TEST:
1037  case GA_NODE_GRAD_TEST:
1038  case GA_NODE_INTERPOLATE_GRAD_TEST:
1039  case GA_NODE_ELEMENTARY_GRAD_TEST:
1040  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
1041  case GA_NODE_XFEM_PLUS_GRAD_TEST:
1042  case GA_NODE_XFEM_MINUS_GRAD_TEST:
1043  case GA_NODE_HESS_TEST:
1044  case GA_NODE_INTERPOLATE_HESS_TEST:
1045  case GA_NODE_ELEMENTARY_HESS_TEST:
1046  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
1047  case GA_NODE_XFEM_PLUS_HESS_TEST:
1048  case GA_NODE_XFEM_MINUS_HESS_TEST:
1049  case GA_NODE_DIVERG_TEST:
1050  case GA_NODE_INTERPOLATE_DIVERG_TEST:
1051  case GA_NODE_ELEMENTARY_DIVERG_TEST:
1052  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
1053  case GA_NODE_XFEM_PLUS_DIVERG_TEST:
1054  case GA_NODE_XFEM_MINUS_DIVERG_TEST:
1055  str << (pnode->test_function_type == 1 ? "Test_" : "Test2_")
1056  << pnode->name;
1057  break;
1058  case GA_NODE_SPEC_FUNC: str << pnode->name; break;
1059  case GA_NODE_OPERATOR:
1060  case GA_NODE_PREDEF_FUNC:
1061  if (pnode->der1) {
1062  str << "Derivative_" << pnode->der1 << "_";
1063  if (pnode->der2) str << pnode->der2 << "_";
1064  }
1065  str << pnode->name; break;
1066  case GA_NODE_ZERO:
1067  GMM_ASSERT1(pnode->test_function_type != size_type(-1),
1068  "Internal error");
1069  if (pnode->test_function_type) str << "(";
1070  ga_print_constant_tensor(pnode, str);
1071  if (pnode->name_test1.size()) {
1072  GMM_ASSERT1(pnode->qdim1 > 0, "Internal error");
1073  if (pnode->qdim1 == 1)
1074  str << "*Test_" << pnode->name_test1;
1075  else {
1076  str << "*(Reshape(Test_" << pnode->name_test1 << ","
1077  << pnode->qdim1<< ")(1))";
1078  }
1079  }
1080  if (pnode->name_test2.size()) {
1081  GMM_ASSERT1(pnode->qdim2 > 0, "Internal error");
1082  if (pnode->qdim2 == 1)
1083  str << "*Test2_" << pnode->name_test2;
1084  else {
1085  str << "*(Reshape(Test2_" << pnode->name_test2 << ","
1086  << pnode->qdim2<< ")(1))";
1087  }
1088  }
1089  if (pnode->test_function_type) str << ")";
1090  break;
1091 
1092  case GA_NODE_CONSTANT:
1093  ga_print_constant_tensor(pnode, str);
1094  break;
1095 
1096  case GA_NODE_ALLINDICES:
1097  str << ":";
1098  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1099  break;
1100 
1101  case GA_NODE_PARAMS:
1102  GMM_ASSERT1(pnode->children.size(), "Invalid tree");
1103  ga_print_node(pnode->children[0], str);
1104  str << "(";
1105  for (size_type i = 1; i < pnode->children.size(); ++i) {
1106  if (i > 1) str << ", ";
1107  ga_print_node(pnode->children[i], str);
1108  }
1109  str << ")";
1110  break;
1111 
1112  case GA_NODE_NAME:
1113  str << pnode->name;
1114  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1115  break;
1116 
1117  case GA_NODE_MACRO_PARAM:
1118  if (pnode->nbc2 == 1) str << "Grad_";
1119  if (pnode->nbc2 == 2) str << "Hess_";
1120  if (pnode->nbc2 == 3) str << "Div_";
1121  if (pnode->nbc3 == 1) str << "Test_";
1122  if (pnode->nbc3 == 2) str << "Test2_";
1123  str << "P" << pnode->nbc1;
1124  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1125  break;
1126 
1127  case GA_NODE_RESHAPE:
1128  str << "Reshape";
1129  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1130  break;
1131 
1132  case GA_NODE_CROSS_PRODUCT:
1133  str << "Cross_product";
1134  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1135  break;
1136 
1137  case GA_NODE_SWAP_IND:
1138  str << "Swap_indices";
1139  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1140  break;
1141 
1142  case GA_NODE_IND_MOVE_LAST:
1143  str << "Index_move_last";
1144  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1145  break;
1146 
1147  case GA_NODE_CONTRACT:
1148  str << "Contract";
1149  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1150  break;
1151 
1152  case GA_NODE_C_MATRIX:
1153  GMM_ASSERT1(pnode->children.size(), "Invalid tree");
1154  GMM_ASSERT1(pnode->nbc1 == pnode->tensor_order(), "Invalid C_MATRIX");
1155  switch (pnode->tensor_order()) {
1156  case 0:
1157  ga_print_node(pnode->children[0], str);
1158  break;
1159 
1160  case 1:
1161  str << "[";
1162  for (size_type i = 0; i < pnode->tensor_proper_size(0); ++i) {
1163  if (i != 0) str << ",";
1164  ga_print_node(pnode->children[i], str);
1165  }
1166  str << "]";
1167  break;
1168 
1169  case 2: case 3: case 4:
1170  {
1171  size_type ii(0);
1172  size_type n0 = pnode->tensor_proper_size(0);
1173  size_type n1 = pnode->tensor_proper_size(1);
1174  size_type n2 = ((pnode->tensor_order() > 2) ?
1175  pnode->tensor_proper_size(2) : 1);
1176  size_type n3 = ((pnode->tensor_order() > 3) ?
1177  pnode->tensor_proper_size(3) : 1);
1178  if (n3 > 1) str << "[";
1179  for (size_type l = 0; l < n3; ++l) {
1180  if (l != 0) str << ",";
1181  if (n2 > 1) str << "[";
1182  for (size_type k = 0; k < n2; ++k) {
1183  if (k != 0) str << ",";
1184  if (n1 > 1) str << "[";
1185  for (size_type j = 0; j < n1; ++j) {
1186  if (j != 0) str << ",";
1187  if (n0 > 1) str << "[";
1188  for (size_type i = 0; i < n0; ++i) {
1189  if (i != 0) str << ",";
1190  ga_print_node(pnode->children[ii++], str);
1191  }
1192  if (n0 > 1) str << "]";
1193  }
1194  if (n1 > 1) str << "]";
1195  }
1196  if (n2 > 1) str << "]";
1197  }
1198  if (n3 > 1) str << "]";
1199  }
1200  break;
1201 
1202  case 5: case 6: case 7: case 8:
1203  str << "Reshape([";
1204  for (size_type i = 0; i < pnode->tensor_proper_size(); ++i) {
1205  if (i != 0) str << ";";
1206  ga_print_node(pnode->children[i], str);
1207  }
1208  str << "]";
1209  for (size_type i = 0; i < pnode->tensor_order(); ++i) {
1210  if (i != 0) str << ", ";
1211  str << pnode->tensor_proper_size(i);
1212  }
1213  str << ")";
1214  break;
1215 
1216  default: GMM_ASSERT1(false, "Invalid tensor dimension");
1217  }
1218  break;
1219 
1220  default:
1221  str << "Invalid or not taken into account node type "
1222  << pnode->node_type;
1223  break;
1224  }
1225 
1226  if (is_interpolate)
1227  str << "," << pnode->interpolate_name << ")";
1228  else if (is_elementary) {
1229  str << "," << pnode->elementary_name;
1230  if (pnode->name.compare(pnode->elementary_target) != 0)
1231  str << "," << pnode->elementary_target;
1232  str << ")";
1233  } else if (is_secondary)
1234  str << ")";
1235  else if (is_xfem_plus || is_xfem_minus)
1236  str << ")";
1237 
1238  str.precision(prec);
1239  }
1240 
1241  std::string ga_tree_to_string(const ga_tree &tree) {
1242  std::stringstream str;
1243  str.precision(16);
1244  if (tree.root) verify_tree(tree.root, 0);
1245  if (tree.root) ga_print_node(tree.root, str); else str << "0";
1246  return str.str();
1247  }
1248 
1249  size_type ga_parse_prefix_operator(std::string &name) {
1250  if (name.size() >= 5 && name.compare(0, 5, "Grad_") == 0)
1251  { name = name.substr(5); return 1; }
1252  else if (name.size() >= 5 && name.compare(0, 5, "Hess_") == 0)
1253  { name = name.substr(5); return 2; }
1254  else if (name.size() >= 4 && name.compare(0, 4, "Div_") == 0)
1255  { name = name.substr(4); return 3; }
1256  return 0;
1257  }
1258 
1259  size_type ga_parse_prefix_test(std::string &name) {
1260  if (name.size() >= 5 && name.compare(0, 5, "Test_") == 0)
1261  { name = name.substr(5); return 1; }
1262  else if (name.size() >= 6 && name.compare(0, 6, "Test2_") == 0)
1263  { name = name.substr(6); return 2; }
1264  return 0;
1265  }
1266 
1267  // 0 : ok
1268  // 1 : function or operator name or "X"
1269  // 2 : reserved prefix Grad, Hess, Div, Derivative_ Test and Test2
1270  // 3 : reserved prefix Dot and Previous
1271  int ga_check_name_validity(const std::string &name) {
1272  if (name.compare(0, 11, "Derivative_") == 0)
1273  return 2;
1274 
1275  const ga_predef_operator_tab &PREDEF_OPERATORS
1276  = dal::singleton<ga_predef_operator_tab>::instance(0);
1277  const ga_spec_function_tab &SPEC_FUNCTIONS
1278  = dal::singleton<ga_spec_function_tab>::instance(0);
1279  const ga_spec_op_tab &SPEC_OP
1280  = dal::singleton<ga_spec_op_tab>::instance(0);
1281  const ga_predef_function_tab &PREDEF_FUNCTIONS
1282  = dal::singleton<ga_predef_function_tab>::instance(0);
1283 
1284  if (SPEC_OP.find(name) != SPEC_OP.end())
1285  return 1;
1286 
1287  if (PREDEF_FUNCTIONS.find(name) != PREDEF_FUNCTIONS.end())
1288  return 1;
1289 
1290  if (SPEC_FUNCTIONS.find(name) != SPEC_FUNCTIONS.end())
1291  return 1;
1292 
1293  if (PREDEF_OPERATORS.tab.find(name) != PREDEF_OPERATORS.tab.end())
1294  return 1;
1295 
1296  if (name.size() >= 5 && name.compare(0, 5, "Grad_") == 0)
1297  return 2;
1298 
1299  if (name.size() >= 5 && name.compare(0, 5, "Hess_") == 0)
1300  return 2;
1301 
1302  if (name.size() >= 4 && name.compare(0, 4, "Div_") == 0)
1303  return 2;
1304 
1305  if (name.size() >= 6 && name.compare(0, 6, "Test2_") == 0)
1306  return 2;
1307 
1308  if (name.size() >= 5 && name.compare(0, 5, "Test_") == 0)
1309  return 2;
1310 
1311 // if (name.size() >= 4 && name.compare(0, 4, "Dot_") == 0)
1312 // return 3;
1313 // if (name.size() >= 5 && name.compare(0, 5, "Dot2_") == 0)
1314 // return 3;
1315 
1316 // if (name.size() >= 9 && name.compare(0, 9, "Previous_") == 0)
1317 // return 3;
1318 // if (name.size() >= 10 && name.compare(0, 10, "Previous2_") == 0)
1319 // return 3;
1320 // if (name.size() >= 12 && name.compare(0, 12, "Previous1_2_") == 0)
1321 // return 3;
1322 
1323  return 0;
1324  }
1325 
1326  //=========================================================================
1327  // Structure dealing with macros.
1328  //=========================================================================
1329 
1330  ga_macro::ga_macro() : ptree(new ga_tree), nbp(0) {}
1331  ga_macro::~ga_macro() { delete ptree; }
1332  ga_macro::ga_macro(const std::string &name, const ga_tree &t, size_type nbp_)
1333  : ptree(new ga_tree(t)), macro_name_(name), nbp(nbp_) {}
1334  ga_macro::ga_macro(const ga_macro &gam)
1335  : ptree(new ga_tree(gam.tree())), macro_name_(gam.name()),
1336  nbp(gam.nb_params()) {}
1337  ga_macro &ga_macro::operator =(const ga_macro &gam) {
1338  delete ptree; ptree = new ga_tree(gam.tree());
1339  macro_name_ = gam.name();
1340  nbp = gam.nb_params();
1341  return *this;
1342  }
1343 
1344  static void ga_replace_macro_params
1345  (ga_tree &tree, pga_tree_node pnode,
1346  const std::vector<pga_tree_node> &children) {
1347  if (!pnode) return;
1348  for (size_type i = 0; i < pnode->children.size(); ++i)
1349  ga_replace_macro_params(tree, pnode->children[i], children);
1350 
1351  if (pnode->node_type == GA_NODE_MACRO_PARAM) {
1352  size_type po = pnode->nbc2;
1353  size_type pt = pnode->nbc3;
1354  GMM_ASSERT1(pnode->nbc1+1 < children.size(), "Internal error");
1355  pga_tree_node pchild = children[pnode->nbc1+1];
1356 
1357  if (po || pt || pnode->op_type != GA_NAME) {
1358  if (!(pchild->children.empty()) || pchild->node_type != GA_NODE_NAME)
1359  ga_throw_error(pchild->expr, pchild->pos, "Error in macro "
1360  "expansion. Only variable name are allowed for macro "
1361  "parameter preceded by Grad_ Hess_ Test_ or Test2_ "
1362  "prefixes.");
1363  switch(pnode->op_type) {
1364  case GA_NAME : pnode->node_type = GA_NODE_NAME; break;
1365  case GA_INTERPOLATE : pnode->node_type = GA_NODE_INTERPOLATE; break;
1366  case GA_INTERPOLATE_DERIVATIVE :
1367  pnode->node_type = GA_NODE_INTERPOLATE_DERIVATIVE; break;
1368  case GA_ELEMENTARY : pnode->node_type = GA_NODE_ELEMENTARY; break;
1369  case GA_SECONDARY_DOMAIN :
1370  pnode->node_type = GA_NODE_SECONDARY_DOMAIN; break;
1371  case GA_XFEM_PLUS : pnode->node_type = GA_NODE_XFEM_PLUS; break;
1372  case GA_XFEM_MINUS: pnode->node_type = GA_NODE_XFEM_MINUS; break;
1373  default:break;
1374  }
1375  pnode->name = pchild->name;
1376  if (pt == 1) pnode->name = "Test_" + pnode->name;
1377  if (pt == 2) pnode->name = "Test2_" + pnode->name;
1378  if (po == 1) pnode->name = "Grad_" + pnode->name;
1379  if (po == 2) pnode->name = "Hess_" + pnode->name;
1380  if (po == 3) pnode->name = "Div_" + pnode->name;
1381  } else {
1382  pga_tree_node pnode_old = pnode;
1383  pnode = nullptr;
1384  tree.copy_node(pchild, pnode_old->parent, pnode);
1385  if (pnode_old->parent)
1386  pnode_old->parent->replace_child(pnode_old, pnode);
1387  else
1388  tree.root = pnode;
1389  GMM_ASSERT1(pnode_old->children.empty(), "Internal error");
1390  delete pnode_old;
1391  }
1392  }
1393  }
1394 
1395  static void ga_expand_macro(ga_tree &tree, pga_tree_node pnode,
1396  const ga_macro_dictionary &macro_dict) {
1397  if (!pnode) return;
1398 
1399  if (pnode->node_type == GA_NODE_PARAMS) {
1400 
1401  for (size_type i = 1; i < pnode->children.size(); ++i)
1402  ga_expand_macro(tree, pnode->children[i], macro_dict);
1403 
1404  if (pnode->children[0]->node_type != GA_NODE_NAME) {
1405  ga_expand_macro(tree, pnode->children[0], macro_dict);
1406  } else {
1407 
1408  if (macro_dict.macro_exists(pnode->children[0]->name)) {
1409 
1410  const ga_macro &gam = macro_dict.get_macro(pnode->children[0]->name);
1411 
1412  if (gam.nb_params()==0) { // Macro without parameters
1413  pga_tree_node pnode_old = pnode->children[0];
1414  pnode->children[0] = nullptr;
1415  tree.copy_node(gam.tree().root,
1416  pnode_old->parent,pnode->children[0]);
1417  GMM_ASSERT1(pnode_old->children.empty(), "Internal error");
1418  delete pnode_old;
1419 
1420  } else { // Macro with parameters
1421 
1422  if (gam.nb_params()+1 != pnode->children.size())
1423  ga_throw_error(pnode->expr, pnode->pos,
1424  "Bad number of parameters in the use of macro '"
1425  << gam.name() << "'. Expected " << gam.nb_params()
1426  << " found " << pnode->children.size()-1 << ".");
1427 
1428  pga_tree_node pnode_old = pnode;
1429  pnode = nullptr;
1430  tree.copy_node(gam.tree().root, pnode_old->parent, pnode);
1431  if (pnode_old->parent)
1432  pnode_old->parent->replace_child(pnode_old, pnode);
1433  else
1434  tree.root = pnode;
1435  ga_replace_macro_params(tree, pnode, pnode_old->children);
1436  }
1437  }
1438  }
1439 
1440  } else if (pnode->node_type == GA_NODE_NAME &&
1441  macro_dict.macro_exists(pnode->name)) {
1442  // Macro without parameters
1443  const ga_macro &gam = macro_dict.get_macro(pnode->name);
1444  if (gam.nb_params() != 0)
1445  ga_throw_error(pnode->expr, pnode->pos,
1446  "Bad number of parameters in the use of macro '"
1447  << gam.name() << "'. Expected " << gam.nb_params()
1448  << " none found.");
1449 
1450  pga_tree_node pnode_old = pnode;
1451  pnode = nullptr;
1452  tree.copy_node(gam.tree().root, pnode_old->parent, pnode);
1453  if (pnode_old->parent)
1454  pnode_old->parent->replace_child(pnode_old, pnode);
1455  else
1456  tree.root = pnode;
1457  GMM_ASSERT1(pnode_old->children.empty(), "Internal error");
1458  delete pnode_old;
1459  } else {
1460  for (size_type i = 0; i < pnode->children.size(); ++i)
1461  ga_expand_macro(tree, pnode->children[i], macro_dict);
1462  }
1463  }
1464 
1465  static void ga_mark_macro_params_rec(const pga_tree_node pnode,
1466  const std::vector<std::string> &params) {
1467  if (!pnode) return;
1468  for (size_type i = 0; i < pnode->children.size(); ++i)
1469  ga_mark_macro_params_rec(pnode->children[i], params);
1470 
1471  if (pnode->node_type == GA_NODE_NAME ||
1472  pnode->node_type == GA_NODE_INTERPOLATE ||
1473  pnode->node_type == GA_NODE_ELEMENTARY ||
1474  pnode->node_type == GA_NODE_SECONDARY_DOMAIN ||
1475  pnode->node_type == GA_NODE_XFEM_PLUS ||
1476  pnode->node_type == GA_NODE_XFEM_MINUS) {
1477  std::string name = pnode->name;
1478  size_type po = ga_parse_prefix_operator(name);
1479  size_type pt = ga_parse_prefix_test(name);
1480 
1481  for (size_type i = 0; i < params.size(); ++i)
1482  if (name.compare(params[i]) == 0) {
1483  pnode->name = name;
1484  switch(pnode->node_type) {
1485  case GA_NODE_NAME : pnode->op_type = GA_NAME; break;
1486  case GA_NODE_INTERPOLATE : pnode->op_type = GA_INTERPOLATE; break;
1487  case GA_NODE_INTERPOLATE_DERIVATIVE :
1488  pnode->op_type = GA_INTERPOLATE_DERIVATIVE; break;
1489  case GA_NODE_ELEMENTARY : pnode->op_type = GA_ELEMENTARY; break;
1490  case GA_NODE_SECONDARY_DOMAIN :
1491  pnode->op_type = GA_SECONDARY_DOMAIN; break;
1492  case GA_NODE_XFEM_PLUS : pnode->op_type = GA_XFEM_PLUS; break;
1493  case GA_NODE_XFEM_MINUS: pnode->op_type = GA_XFEM_MINUS; break;
1494  default:break;
1495  }
1496  pnode->node_type = GA_NODE_MACRO_PARAM;
1497  pnode->nbc1 = i; pnode->nbc2 = po; pnode->nbc3 = pt;
1498  }
1499  }
1500  }
1501 
1502  static void ga_mark_macro_params(ga_macro &gam,
1503  const std::vector<std::string> &params,
1504  const ga_macro_dictionary &macro_dict) {
1505  if (gam.tree().root) {
1506  ga_mark_macro_params_rec(gam.tree().root, params);
1507  ga_expand_macro(gam.tree(), gam.tree().root, macro_dict);
1508  }
1509  }
1510 
1511  bool ga_macro_dictionary::macro_exists(const std::string &name) const {
1512  if (macros.find(name) != macros.end()) return true;
1513  if (parent && parent->macro_exists(name)) return true;
1514  return false;
1515  }
1516 
1517  const ga_macro &
1518  ga_macro_dictionary::get_macro(const std::string &name) const {
1519  auto it = macros.find(name);
1520  if (it != macros.end()) return it->second;
1521  if (parent) return parent->get_macro(name);
1522  GMM_ASSERT1(false, "Undefined macro");
1523  }
1524 
1525  void ga_macro_dictionary::add_macro(const ga_macro &gam)
1526  { macros[gam.name()] = gam; }
1527 
1528  void ga_macro_dictionary::add_macro(const std::string &name,
1529  const std::string &expr)
1530  { ga_tree tree; ga_read_string_reg("Def "+name+":="+expr, tree, *this); }
1531 
1532  void ga_macro_dictionary::del_macro(const std::string &name) {
1533  auto it = macros.find(name);
1534  GMM_ASSERT1(it != macros.end(), "Undefined macro (at this level)");
1535  macros.erase(it);
1536  }
1537 
1538 
1539  //=========================================================================
1540  // Syntax analysis for the generic assembly language
1541  //=========================================================================
1542 
1543  // Read a term with an (implicit) pushdown automaton.
1544  static GA_TOKEN_TYPE ga_read_term(pstring expr, size_type &pos,
1545  ga_tree &tree,
1546  ga_macro_dictionary &macro_dict) {
1547  size_type token_pos, token_length;
1548  GA_TOKEN_TYPE t_type;
1549  int state = 1; // 1 = reading term, 2 = reading after term
1550 
1551  for (;;) {
1552 
1553  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1554 
1555  switch (state) {
1556 
1557  case 1:
1558  switch (t_type) {
1559  case GA_SCALAR:
1560  {
1561  char *endptr; const char *nptr = &((*expr)[token_pos]);
1562  scalar_type s_read = ::strtod(nptr, &endptr);
1563  if (endptr == nptr)
1564  ga_throw_error(expr, token_pos, "Bad numeric format.");
1565  tree.add_scalar(s_read, token_pos, expr);
1566  }
1567  state = 2; break;
1568 
1569  case GA_COLON:
1570  tree.add_allindices(token_pos, expr);
1571  state = 2; break;
1572 
1573  case GA_NAME:
1574  tree.add_name(&((*expr)[token_pos]), token_length, token_pos, expr);
1575  state = 2; break;
1576 
1577  case GA_MINUS: // unary -
1578  tree.add_op(GA_UNARY_MINUS, token_pos, expr);
1579  state = 1; break;
1580 
1581  case GA_PLUS: // unary +
1582  state = 1; break;
1583 
1584  case GA_SYM:
1585  tree.add_op(GA_SYM, token_pos, expr);
1586  state = 1; break;
1587 
1588  case GA_SKEW:
1589  tree.add_op(GA_SKEW, token_pos, expr);
1590  state = 1; break;
1591 
1592  case GA_TRACE:
1593  tree.add_op(GA_TRACE, token_pos, expr);
1594  state = 1; break;
1595 
1596  case GA_DEVIATOR:
1597  tree.add_op(GA_DEVIATOR, token_pos, expr);
1598  state = 1; break;
1599 
1600  case GA_DEF:
1601  {
1602  ga_macro gam;
1603  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1604  if (t_type != GA_NAME)
1605  ga_throw_error(expr, pos,
1606  "Macro definition should begin with macro name");
1607  gam.name() = std::string(&((*expr)[token_pos]), token_length);
1608  if (ga_check_name_validity(gam.name()))
1609  ga_throw_error(expr, pos-1, "Invalid macro name.")
1610  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1611  std::vector<std::string> params;
1612  if (t_type == GA_LPAR) {
1613  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1614  while (t_type == GA_NAME) {
1615  params.push_back(std::string(&((*expr)[token_pos]),
1616  token_length));
1617  if (ga_check_name_validity(params.back()))
1618  ga_throw_error(expr, pos-1, "Invalid macro parameter name.");
1619  for (size_type i = 0; i+1 < params.size(); ++i)
1620  if (params.back().compare(params[i]) == 0)
1621  ga_throw_error(expr, pos-1,
1622  "Invalid repeated macro parameter name.");
1623  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1624  if (t_type == GA_COMMA)
1625  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1626  }
1627  if (t_type != GA_RPAR)
1628  ga_throw_error(expr, pos-1,
1629  "Missing right parenthesis in macro definition.");
1630  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1631  }
1632  if (t_type != GA_COLON_EQ)
1633  ga_throw_error(expr, pos-1, "Missing := for macro definition.");
1634 
1635  t_type = ga_read_term(expr, pos, gam.tree(), macro_dict);
1636  if (gam.tree().root)
1637  ga_expand_macro(gam.tree(), gam.tree().root, macro_dict);
1638  gam.nb_params() = params.size();
1639  if (params.size())
1640  ga_mark_macro_params(gam, params, macro_dict);
1641  macro_dict.add_macro(gam);
1642 
1643  // cout << "macro \"" << gam.name() << "\" registered with "
1644  // << gam.nb_params() << " params := "
1645  // << ga_tree_to_string(gam.tree()) << endl;
1646 
1647  if (t_type == GA_END) return t_type;
1648  else if (t_type != GA_SEMICOLON)
1649  ga_throw_error(expr, pos-1,
1650  "Syntax error at the end of macro definition.");
1651  state = 1;
1652  }
1653  break;
1654 
1655  case GA_INTERPOLATE:
1656  {
1657  tree.add_scalar(scalar_type(0), token_pos, expr);
1658  tree.current_node->node_type = GA_NODE_INTERPOLATE;
1659  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1660  if (t_type != GA_LPAR)
1661  ga_throw_error(expr, pos-1, "Missing interpolate arguments.");
1662  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1663  if (t_type != GA_NAME)
1664  ga_throw_error(expr, pos,
1665  "First argument of Interpolate should be a "
1666  "variable, test function, X or Normal.");
1667  tree.current_node->name = std::string(&((*expr)[token_pos]),
1668  token_length);
1669 
1670  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1671  if (t_type != GA_COMMA)
1672  ga_throw_error(expr, pos, "Bad format for Interpolate "
1673  "arguments.");
1674  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1675  if (t_type != GA_NAME)
1676  ga_throw_error(expr, pos,
1677  "Second argument of Interpolate should be a "
1678  "transformation name.");
1679  tree.current_node->interpolate_name
1680  = std::string(&((*expr)[token_pos]), token_length);
1681  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1682  if (t_type != GA_RPAR)
1683  ga_throw_error(expr, pos-1, "Missing a parenthesis after "
1684  "interpolate arguments.");
1685  state = 2;
1686  }
1687  break;
1688 
1689  case GA_INTERPOLATE_DERIVATIVE:
1690  {
1691  tree.add_scalar(scalar_type(0), token_pos, expr);
1692  tree.current_node->node_type = GA_NODE_INTERPOLATE_DERIVATIVE;
1693  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1694  if (t_type != GA_LPAR)
1695  ga_throw_error(expr, pos-1,
1696  "Missing Interpolate_derivative arguments.");
1697  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1698  if (t_type != GA_NAME)
1699  ga_throw_error(expr, pos,
1700  "First argument of Interpolate should the "
1701  "interpolate transformtion name ");
1702  tree.current_node->interpolate_name_der
1703  = std::string(&((*expr)[token_pos]), token_length);
1704  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1705  if (t_type != GA_COMMA)
1706  ga_throw_error(expr, pos, "Bad format for Interpolate_derivative "
1707  "arguments.");
1708  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1709  if (t_type != GA_NAME)
1710  ga_throw_error(expr, pos,
1711  "Second argument of Interpolate should be a "
1712  "variable name.");
1713  tree.current_node->name
1714  = std::string(&((*expr)[token_pos]), token_length);
1715 
1716  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1717  tree.current_node->interpolate_name = "";
1718  if (t_type == GA_COMMA) {
1719  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1720  if (t_type != GA_NAME)
1721  ga_throw_error(expr, pos,
1722  "Third argument of Interpolate should be a "
1723  "interpolate transformation name.");
1724  tree.current_node->interpolate_name
1725  = std::string(&((*expr)[token_pos]), token_length);
1726  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1727  }
1728  if (t_type != GA_RPAR)
1729  ga_throw_error(expr, pos-1, "Missing a parenthesis after "
1730  "Interpolate_derivative arguments.");
1731  state = 2;
1732  }
1733  break;
1734 
1735  case GA_ELEMENTARY:
1736  {
1737  tree.add_scalar(scalar_type(0), token_pos, expr);
1738  tree.current_node->node_type = GA_NODE_ELEMENTARY;
1739  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1740  if (t_type != GA_LPAR)
1741  ga_throw_error(expr, pos-1,
1742  "Missing Elementary_transformation arguments.");
1743  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1744  if (t_type != GA_NAME)
1745  ga_throw_error(expr, pos,
1746  "First argument of Elementary_transformation "
1747  "should be a variable or a test function.");
1748  tree.current_node->name = std::string(&((*expr)[token_pos]),
1749  token_length);
1750  tree.current_node->elementary_target = tree.current_node->name;
1751 
1752  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1753  if (t_type != GA_COMMA)
1754  ga_throw_error(expr, pos, "Bad format for "
1755  "Elementary_transformation arguments.");
1756  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1757  if (t_type != GA_NAME)
1758  ga_throw_error(expr, pos,
1759  "Second argument of Elementary_transformation "
1760  "should be a transformation name.");
1761  tree.current_node->elementary_name
1762  = std::string(&((*expr)[token_pos]), token_length);
1763  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1764 
1765  if (t_type == GA_COMMA) {
1766  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1767  if (t_type != GA_NAME)
1768  ga_throw_error(expr, pos,
1769  "Third argument of Elementary_transformation "
1770  "should be a variable or data name.");
1771 
1772  tree.current_node->elementary_target =
1773  std::string(&((*expr)[token_pos]), token_length);
1774  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1775  }
1776 
1777  if (t_type != GA_RPAR)
1778  ga_throw_error(expr, pos-1, "Missing a parenthesis after "
1779  "Elementary_transformation arguments.");
1780  state = 2;
1781  }
1782  break;
1783 
1784  case GA_SECONDARY_DOMAIN:
1785  {
1786  tree.add_scalar(scalar_type(0), token_pos, expr);
1787  tree.current_node->node_type = GA_NODE_SECONDARY_DOMAIN;
1788  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1789  if (t_type != GA_LPAR)
1790  ga_throw_error(expr, pos-1,"Missing Secondary_domain arguments.");
1791  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1792  if (t_type != GA_NAME)
1793  ga_throw_error(expr, pos,
1794  "First argument of Secondary_domain should be a "
1795  "variable, test function, X or Normal.");
1796  tree.current_node->name = std::string(&((*expr)[token_pos]),
1797  token_length);
1798  tree.current_node->interpolate_name = tree.secondary_domain;
1799  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1800  if (t_type != GA_RPAR)
1801  ga_throw_error(expr, pos-1, "Missing a parenthesis after "
1802  "Secondary_domain arguments.");
1803  state = 2;
1804  }
1805  break;
1806 
1807  case GA_XFEM_PLUS:
1808  {
1809  tree.add_scalar(scalar_type(0), token_pos, expr);
1810  tree.current_node->node_type = GA_NODE_XFEM_PLUS;
1811  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1812  if (t_type != GA_LPAR)
1813  ga_throw_error(expr, pos-1,
1814  "Missing Xfem_plus arguments.");
1815  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1816  if (t_type != GA_NAME)
1817  ga_throw_error(expr, pos,
1818  "The argument of Xfem_plus should be a "
1819  "variable or a test function.");
1820  tree.current_node->name = std::string(&((*expr)[token_pos]),
1821  token_length);
1822  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1823  if (t_type != GA_RPAR)
1824  ga_throw_error(expr, pos-1, "Missing a parenthesis after "
1825  "Xfem_plus argument.");
1826  state = 2;
1827  }
1828  break;
1829 
1830  case GA_XFEM_MINUS:
1831  {
1832  tree.add_scalar(scalar_type(0), token_pos, expr);
1833  tree.current_node->node_type = GA_NODE_XFEM_MINUS;
1834  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1835  if (t_type != GA_LPAR)
1836  ga_throw_error(expr, pos-1,
1837  "Missing Xfem_minus arguments.");
1838  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1839  if (t_type != GA_NAME)
1840  ga_throw_error(expr, pos,
1841  "The argument of Xfem_minus should be a "
1842  "variable or a test function.");
1843  tree.current_node->name = std::string(&((*expr)[token_pos]),
1844  token_length);
1845  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1846  if (t_type != GA_RPAR)
1847  ga_throw_error(expr, pos-1, "Missing a parenthesis after "
1848  "Xfem_minus argument.");
1849  state = 2;
1850  }
1851  break;
1852 
1853  case GA_INTERPOLATE_FILTER:
1854  {
1855  tree.add_scalar(scalar_type(0), token_pos, expr);
1856  tree.current_node->node_type = GA_NODE_INTERPOLATE_FILTER;
1857  tree.current_node->nbc1 = size_type(-1);
1858  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1859  if (t_type != GA_LPAR)
1860  ga_throw_error(expr, pos-1, "Missing interpolate arguments.");
1861  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1862  if (t_type != GA_NAME)
1863  ga_throw_error(expr, pos, "First argument of Interpolate_filter "
1864  "should be a transformation name.");
1865  tree.current_node->interpolate_name
1866  = std::string(&((*expr)[token_pos]), token_length);
1867  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1868  if (t_type != GA_COMMA)
1869  ga_throw_error(expr, pos,
1870  "Bad format for Interpolate_filter arguments.");
1871  ga_tree sub_tree;
1872  t_type = ga_read_term(expr, pos, sub_tree, macro_dict);
1873  if (t_type != GA_RPAR && t_type != GA_COMMA)
1874  ga_throw_error(expr, pos-1,
1875  "Bad format for Interpolate_filter arguments.");
1876  tree.add_sub_tree(sub_tree);
1877  if (t_type == GA_COMMA) {
1878  ga_tree sub_tree2;
1879  t_type = ga_read_term(expr, pos, sub_tree2, macro_dict);
1880  tree.add_sub_tree(sub_tree2);
1881  }
1882  if (t_type != GA_RPAR)
1883  ga_throw_error(expr, pos-1, "Unbalanced parenthesis.");
1884  state = 2;
1885  }
1886  break;
1887 
1888  case GA_PRINT:
1889  tree.add_op(GA_PRINT, token_pos, expr);
1890  state = 1; break;
1891 
1892  case GA_LPAR: // Parenthesed expression
1893  {
1894  ga_tree sub_tree;
1895  GA_TOKEN_TYPE r_type;
1896  r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
1897  if (r_type != GA_RPAR)
1898  ga_throw_error(expr, pos-1, "Unbalanced parenthesis.");
1899  tree.add_sub_tree(sub_tree);
1900  state = 2;
1901  }
1902  break;
1903 
1904  case GA_LBRACKET: // Explicit vector/matrix or tensor
1905  {
1906  ga_tree sub_tree;
1907  GA_TOKEN_TYPE r_type;
1908  size_type nbc1(0), nbc2(0), nbc3(0), n1(0), n2(0), n3(0);
1909  size_type tensor_order(1);
1910  bool foundcomma(false), foundsemi(false);
1911 
1912  r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
1913  size_type nb_comp = 0;
1914  tree.add_matrix(token_pos, expr);
1915 
1916  if (sub_tree.root->node_type == GA_NODE_C_MATRIX) { // nested format
1917  bgeot::multi_index mii;
1918  do {
1919  if (nb_comp) {
1920  sub_tree.clear();
1921  r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
1922  }
1923  // in the nested format only "," and "]" are expected
1924  if (sub_tree.root->node_type != GA_NODE_C_MATRIX ||
1925  (r_type != GA_COMMA && r_type != GA_RBRACKET))
1926  ga_throw_error(expr, pos-1, "Bad explicit "
1927  "vector/matrix/tensor format.");
1928 
1929  // convert a row vector [a,b] to a column vector [a;b]
1930  if (sub_tree.root->marked &&
1931  sub_tree.root->tensor().sizes()[0] == 1 &&
1932  sub_tree.root->tensor().size() != 1) {
1933  bgeot::multi_index mi = sub_tree.root->tensor().sizes();
1934  for (size_type i = mi.size()-1; i > 0; i--)
1935  mi[i-1] = mi[i];
1936  mi.pop_back();
1937  sub_tree.root->tensor().adjust_sizes(mi);
1938  }
1939  if (!nb_comp) mii = sub_tree.root->tensor().sizes();
1940  else {
1941  const bgeot::multi_index &mi=sub_tree.root->tensor().sizes();
1942  bool cmp = true;
1943  if (mii.size() == mi.size()) {
1944  for (size_type i = 0; i < mi.size(); ++i)
1945  if (mi[i] != mii[i]) cmp = false;
1946  } else cmp = false;
1947  if (!cmp)
1948  ga_throw_error(expr, pos-1, "Bad explicit "
1949  "vector/matrix/tensor format.");
1950  }
1951  for (size_type i = 0; i < sub_tree.root->children.size(); ++i) {
1952  sub_tree.root->children[i]->parent = tree.current_node;
1953  tree.current_node->children.push_back
1954  (sub_tree.root->children[i]);
1955  }
1956  sub_tree.root->children.resize(0);
1957  nb_comp++;
1958  } while (r_type != GA_RBRACKET);
1959  tree.current_node->marked = false;
1960  mii.push_back(nb_comp);
1961  tree.current_node->tensor().adjust_sizes(mii);
1962  } else { // non nested format
1963  do {
1964  if (nb_comp) {
1965  sub_tree.clear();
1966  r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
1967  }
1968  nb_comp++;
1969 
1970  tree.add_sub_tree(sub_tree);
1971 
1972  ++n1; ++n2; ++n3;
1973  if (tensor_order < 2) ++nbc1;
1974  if (tensor_order < 3) ++nbc2;
1975  if (tensor_order < 4) ++nbc3;
1976 
1977  if (r_type == GA_COMMA) {
1978  if (!foundcomma && tensor_order > 1)
1979  ga_throw_error(expr, pos-1, "Bad explicit "
1980  "vector/matrix/tensor format.");
1981  foundcomma = true;
1982  } else if (r_type == GA_SEMICOLON) {
1983  if (n1 != nbc1)
1984  ga_throw_error(expr, pos-1, "Bad explicit "
1985  "vector/matrix/tensor format.");
1986  n1 = 0;
1987  tensor_order = std::max(tensor_order, size_type(2));
1988  } else if (r_type == GA_DCOMMA) {
1989  if (n1 != nbc1 || n2 != nbc2)
1990  ga_throw_error(expr, pos-1, "Bad explicit "
1991  "vector/matrix/tensor format.");
1992  foundsemi = true;
1993  n2 = n1 = 0;
1994  tensor_order = std::max(tensor_order, size_type(3));
1995  } else if (r_type == GA_DSEMICOLON) {
1996  if (n1 != nbc1 || n2 != nbc2 || n3 != nbc3 ||
1997  tensor_order < 3)
1998  ga_throw_error(expr, pos-1, "Bad explicit "
1999  "vector/matrix/tensor format.");
2000  n3 = n2 = n1 = 0;
2001  tensor_order = std::max(tensor_order, size_type(4));
2002  } else if (r_type == GA_RBRACKET) {
2003  if (n1 != nbc1 || n2 != nbc2 || n3 != nbc3 ||
2004  tensor_order == 3)
2005  ga_throw_error(expr, pos-1, "Bad explicit "
2006  "vector/matrix/tensor format.");
2007  tree.current_node->nbc1 = nbc1;
2008  if (tensor_order == 4) {
2009  tree.current_node->nbc2 = nbc2/nbc1;
2010  tree.current_node->nbc3 = nbc3/nbc2;
2011  } else {
2012  tree.current_node->nbc2 = tree.current_node->nbc3 = 1;
2013  }
2014  } else {
2015  ga_throw_error(expr, pos-1, "The explicit "
2016  "vector/matrix/tensor components should be "
2017  "separated by ',', ';', ',,' and ';;' and "
2018  "be ended by ']'.");
2019  }
2020 
2021  } while (r_type != GA_RBRACKET);
2022  bgeot::multi_index mi;
2023  nbc1 = tree.current_node->nbc1; nbc2 = tree.current_node->nbc2;
2024  nbc3 = tree.current_node->nbc3;
2025 
2026  size_type nbl = tree.current_node->children.size()
2027  / (nbc2 * nbc1 * nbc3);
2028  switch(tensor_order) {
2029  case 1:
2030  /* mi.push_back(1); */ mi.push_back(nbc1); break;
2031  case 2:
2032  mi.push_back(nbl); if (nbc1 > 1) mi.push_back(nbc1); break;
2033  case 3:
2034  mi.push_back(nbl); mi.push_back(nbc2);
2035  mi.push_back(nbc1);
2036  break;
2037  case 4:
2038  mi.push_back(nbl); mi.push_back(nbc3);
2039  mi.push_back(nbc2); mi.push_back(nbc1);
2040  break;
2041  default: GMM_ASSERT1(false, "Internal error");
2042  }
2043  tree.current_node->tensor().adjust_sizes(mi);
2044  std::vector<pga_tree_node> children = tree.current_node->children;
2045  auto it = tree.current_node->children.begin();
2046  for (size_type i = 0; i < nbc1; ++i)
2047  for (size_type j = 0; j < nbc2; ++j)
2048  for (size_type k = 0; k < nbc3; ++k)
2049  for (size_type l = 0; l < nbl; ++l, ++it)
2050  *it = children[i+nbc1*(j+nbc2*(k+nbc3*l))];
2051  tree.current_node->marked = true;
2052  }
2053  }
2054  tree.current_node->nbc1 = tree.current_node->tensor().sizes().size();
2055  state = 2;
2056  break;
2057 
2058  default:
2059  ga_throw_error(expr, token_pos, "Unexpected token.");
2060  }
2061  break;
2062 
2063  case 2:
2064  switch (t_type) {
2065  case GA_PLUS: case GA_MINUS: case GA_MULT: case GA_DIV:
2066  case GA_COLON: case GA_DOT: case GA_DOTMULT: case GA_DOTDIV:
2067  case GA_TMULT:
2068  tree.add_op(t_type, token_pos, expr);
2069  state = 1; break;
2070  case GA_QUOTE:
2071  tree.add_op(t_type, token_pos, expr);
2072  state = 2; break;
2073  case GA_END: case GA_RPAR: case GA_COMMA: case GA_DCOMMA:
2074  case GA_RBRACKET: case GA_SEMICOLON: case GA_DSEMICOLON:
2075  return t_type;
2076  case GA_LPAR: // Parameter list
2077  {
2078  ga_tree sub_tree;
2079  GA_TOKEN_TYPE r_type;
2080  tree.add_params(token_pos, expr);
2081  do {
2082  r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
2083  if (r_type != GA_RPAR && r_type != GA_COMMA)
2084  ga_throw_error(expr, pos-((r_type != GA_END)?1:0),
2085  "Parameters should be separated "
2086  "by ',' and parameter list ended by ')'.");
2087  tree.add_sub_tree(sub_tree);
2088  } while (r_type != GA_RPAR);
2089  state = 2;
2090  }
2091  break;
2092 
2093  default:
2094  ga_throw_error(expr, token_pos, "Unexpected token.");
2095  }
2096  break;
2097  }
2098  }
2099 
2100  return GA_INVALID;
2101  }
2102 
2103  // Syntax analysis of a string. Conversion to a tree. register the macros.
2104  void ga_read_string_reg(const std::string &expr, ga_tree &tree,
2105  ga_macro_dictionary &macro_dict) {
2106  size_type pos = 0, token_pos, token_length;
2107  tree.clear();
2108  GA_TOKEN_TYPE t = ga_get_token(expr, pos, token_pos, token_length);
2109  if (t == GA_END) return;
2110  pos = 0;
2111  pstring nexpr(new std::string(expr));
2112 
2113  t = ga_read_term(nexpr, pos, tree, macro_dict);
2114  if (tree.root) ga_expand_macro(tree, tree.root, macro_dict);
2115 
2116  switch (t) {
2117  case GA_RPAR: ga_throw_error(nexpr, pos-1, "Unbalanced parenthesis.");
2118  case GA_RBRACKET: ga_throw_error(nexpr, pos-1, "Unbalanced braket.");
2119  case GA_END: break;
2120  default: ga_throw_error(nexpr, pos-1, "Unexpected token.");
2121  }
2122  }
2123 
2124  // Syntax analysis of a string. Conversion to a tree.
2125  // Do not register the macros (but expand them).
2126  void ga_read_string(const std::string &expr, ga_tree &tree,
2127  const ga_macro_dictionary &macro_dict) {
2128  ga_macro_dictionary macro_dict_loc(true, macro_dict);
2129  ga_read_string_reg(expr, tree, macro_dict_loc);
2130  }
2131 
2132  // Small tool to make basic substitutions into an assembly string
2133  // Should be replaced by macros now.
2134  std::string ga_substitute(const std::string &expr,
2135  const std::map<std::string, std::string> &dict) {
2136  if (dict.size()) {
2137  size_type pos = 0, token_pos, token_length;
2138  std::stringstream exprs;
2139 
2140  while (true) {
2141  GA_TOKEN_TYPE t_type = ga_get_token(expr, pos, token_pos, token_length);
2142  if (t_type == GA_END) return exprs.str();
2143  std::string name(&(expr[token_pos]), token_length);
2144  if (t_type == GA_NAME && dict.find(name) != dict.end())
2145  exprs << dict.at(name); else exprs << name;
2146  }
2147  }
2148  return expr;
2149  }
2150 
2151 } /* end of namespace */
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:165
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem_generic_assembly_tree.h
Compilation and execution operations.
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46