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