GetFEM  5.4.2
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:
556  case GA_NODE_SECONDARY_DOMAIN_X:
557  case GA_NODE_INTERPOLATE_ELT_K: case GA_NODE_INTERPOLATE_ELT_B:
558  case GA_NODE_INTERPOLATE_NORMAL:
559  case GA_NODE_SECONDARY_DOMAIN_NORMAL:
560  if (pnode1->interpolate_name.compare(pnode2->interpolate_name))
561  return false;
562  break;
563  case GA_NODE_INTERPOLATE_DERIVATIVE:
564  if (pnode1->interpolate_name_der.compare(pnode2->interpolate_name_der))
565  return false;
566  if (pnode1->interpolate_name.compare(pnode2->interpolate_name) ||
567  pnode1->elementary_name.compare(pnode2->elementary_name))
568  return false;
569  {
570  const mesh_fem *mf1 = workspace.associated_mf(pnode1->name);
571  const mesh_fem *mf2 = workspace.associated_mf(pnode2->name);
572  switch (version) {
573  case 0:
574  if (pnode1->name.compare(pnode2->name) ||
575  pnode1->test_function_type != pnode2->test_function_type)
576  return false;
577  break;
578  case 1:
579  if (mf1 != mf2 ||
580  workspace.qdim(pnode1->name) != workspace.qdim(pnode2->name) ||
581  pnode1->test_function_type != pnode2->test_function_type)
582  return false;
583  break;
584  case 2:
585  if (mf1 != mf2 ||
586  workspace.qdim(pnode1->name) != workspace.qdim(pnode2->name) ||
587  pnode1->test_function_type == pnode2->test_function_type)
588  return false;
589  break;
590  }
591  }
592  break;
593  case GA_NODE_INTERPOLATE_VAL_TEST: case GA_NODE_INTERPOLATE_GRAD_TEST:
594  case GA_NODE_INTERPOLATE_HESS_TEST: case GA_NODE_INTERPOLATE_DIVERG_TEST:
595  case GA_NODE_ELEMENTARY_VAL_TEST: case GA_NODE_ELEMENTARY_GRAD_TEST:
596  case GA_NODE_ELEMENTARY_HESS_TEST: case GA_NODE_ELEMENTARY_DIVERG_TEST:
597  case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
598  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
599  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
600  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
601  case GA_NODE_XFEM_PLUS_VAL_TEST: case GA_NODE_XFEM_PLUS_GRAD_TEST:
602  case GA_NODE_XFEM_PLUS_HESS_TEST: case GA_NODE_XFEM_PLUS_DIVERG_TEST:
603  case GA_NODE_XFEM_MINUS_VAL_TEST: case GA_NODE_XFEM_MINUS_GRAD_TEST:
604  case GA_NODE_XFEM_MINUS_HESS_TEST: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
605  if (pnode1->interpolate_name.compare(pnode2->interpolate_name) ||
606  pnode1->elementary_name.compare(pnode2->elementary_name))
607  return false;
608  {
609  const mesh_fem *mf1 = workspace.associated_mf(pnode1->name);
610  const mesh_fem *mf2 = workspace.associated_mf(pnode2->name);
611  switch (version) {
612  case 0:
613  if (pnode1->name.compare(pnode2->name) ||
614  pnode1->test_function_type != pnode2->test_function_type)
615  return false;
616  break;
617  case 1:
618  if (mf1 != mf2 ||
619  workspace.qdim(pnode1->name) != workspace.qdim(pnode2->name) ||
620  pnode1->test_function_type != pnode2->test_function_type)
621  return false;
622  break;
623  case 2:
624  if (mf1 != mf2 ||
625  workspace.qdim(pnode1->name) != workspace.qdim(pnode2->name) ||
626  pnode1->test_function_type == pnode2->test_function_type)
627  return false;
628  break;
629  }
630  }
631  break;
632  case GA_NODE_VAL_TEST: case GA_NODE_GRAD_TEST:
633  case GA_NODE_HESS_TEST: case GA_NODE_DIVERG_TEST:
634  {
635  const mesh_fem *mf1 = workspace.associated_mf(pnode1->name);
636  const mesh_fem *mf2 = workspace.associated_mf(pnode2->name);
637  switch (version) {
638  case 0:
639  if (pnode1->name.compare(pnode2->name) ||
640  pnode1->test_function_type != pnode2->test_function_type)
641  return false;
642  break;
643  case 1:
644  if (mf1 != mf2 ||
645  workspace.qdim(pnode1->name) != workspace.qdim(pnode2->name) ||
646  pnode1->test_function_type != pnode2->test_function_type)
647  return false;
648  break;
649  case 2:
650  if (mf1 != mf2 ||
651  workspace.qdim(pnode1->name) != workspace.qdim(pnode2->name) ||
652  pnode1->test_function_type == pnode2->test_function_type)
653  return false;
654  break;
655  }
656  }
657  break;
658  case GA_NODE_VAL: case GA_NODE_GRAD:
659  case GA_NODE_HESS: case GA_NODE_DIVERG:
660  if (pnode1->name.compare(pnode2->name)) return false;
661  break;
662  case GA_NODE_INTERPOLATE_VAL: case GA_NODE_INTERPOLATE_GRAD:
663  case GA_NODE_INTERPOLATE_HESS: case GA_NODE_INTERPOLATE_DIVERG:
664  case GA_NODE_ELEMENTARY_VAL: case GA_NODE_ELEMENTARY_GRAD:
665  case GA_NODE_ELEMENTARY_HESS: case GA_NODE_ELEMENTARY_DIVERG:
666  case GA_NODE_SECONDARY_DOMAIN_VAL: case GA_NODE_SECONDARY_DOMAIN_GRAD:
667  case GA_NODE_SECONDARY_DOMAIN_HESS: case GA_NODE_SECONDARY_DOMAIN_DIVERG:
668  case GA_NODE_XFEM_PLUS_VAL: case GA_NODE_XFEM_PLUS_GRAD:
669  case GA_NODE_XFEM_PLUS_HESS: case GA_NODE_XFEM_PLUS_DIVERG:
670  case GA_NODE_XFEM_MINUS_VAL: case GA_NODE_XFEM_MINUS_GRAD:
671  case GA_NODE_XFEM_MINUS_HESS: case GA_NODE_XFEM_MINUS_DIVERG:
672  if (pnode1->interpolate_name.compare(pnode2->interpolate_name) ||
673  pnode1->elementary_name.compare(pnode2->elementary_name) ||
674  pnode1->name.compare(pnode2->name))
675  return false;
676  break;
677  case GA_NODE_X:
678  if (pnode1->nbc1 != pnode2->nbc1) return false;
679  break;
680 
681  default:break;
682  }
683 
684  if (version && ntype1 == GA_NODE_OP && pnode1->symmetric_op) {
685  if (sub_tree_are_equal(pnode1->children[0], pnode2->children[0],
686  workspace, version) &&
687  sub_tree_are_equal(pnode1->children[1], pnode2->children[1],
688  workspace, version))
689  return true;
690  if (sub_tree_are_equal(pnode1->children[1], pnode2->children[0],
691  workspace, version) &&
692  sub_tree_are_equal(pnode1->children[0], pnode2->children[1],
693  workspace, version) )
694  return true;
695  return false;
696  } else {
697  for (size_type i = 0; i < pnode1->children.size(); ++i)
698  if (!(sub_tree_are_equal(pnode1->children[i], pnode2->children[i],
699  workspace, version)))
700  return false;
701  }
702  return true;
703  }
704 
705  static void verify_tree(const pga_tree_node pnode,
706  const pga_tree_node parent) {
707  GMM_ASSERT1(pnode->parent == parent,
708  "Invalid tree node " << pnode->node_type);
709  for (pga_tree_node &child : pnode->children)
710  verify_tree(child, pnode);
711  }
712 
713 
714  static void ga_print_constant_tensor(const pga_tree_node pnode,
715  std::ostream &str) {
716  size_type nt = pnode->nb_test_functions(); // for printing zero tensors
717  switch (pnode->tensor_order()) {
718  case 0:
719  str << (nt ? scalar_type(0) : pnode->tensor()[0]);
720  break;
721 
722  case 1:
723  str << "[";
724  for (size_type i = 0; i < pnode->tensor_proper_size(0); ++i) {
725  if (i != 0) str << ",";
726  str << (nt ? scalar_type(0) : pnode->tensor()[i]);
727  }
728  str << "]";
729  break;
730 
731  case 2: case 3: case 4:
732  {
733  size_type ii(0);
734  size_type n0 = pnode->tensor_proper_size(0);
735  size_type n1 = pnode->tensor_proper_size(1);
736  size_type n2 = ((pnode->tensor_order() > 2) ?
737  pnode->tensor_proper_size(2) : 1);
738  size_type n3 = ((pnode->tensor_order() > 3) ?
739  pnode->tensor_proper_size(3) : 1);
740  if (n3 > 1) str << "[";
741  for (size_type l = 0; l < n3; ++l) {
742  if (l != 0) str << ",";
743  if (n2 > 1) str << "[";
744  for (size_type k = 0; k < n2; ++k) {
745  if (k != 0) str << ",";
746  if (n1 > 1) str << "[";
747  for (size_type j = 0; j < n1; ++j) {
748  if (j != 0) str << ",";
749  if (n0 > 1) str << "[";
750  for (size_type i = 0; i < n0; ++i) {
751  if (i != 0) str << ",";
752  str << (nt ? scalar_type(0) : pnode->tensor()[ii++]);
753  }
754  if (n0 > 1) str << "]";
755  }
756  if (n1 > 1) str << "]";
757  }
758  if (n2 > 1) str << "]";
759  }
760  if (n3 > 1) str << "]";
761  }
762  break;
763 
764  case 5: case 6: case 7: case 8:
765  str << "Reshape([";
766  for (size_type i = 0; i < pnode->tensor_proper_size(); ++i) {
767  if (i != 0) str << ";";
768  str << (nt ? scalar_type(0) : pnode->tensor()[i]);
769  }
770  str << "]";
771  for (size_type i = 0; i < pnode->tensor_order(); ++i) {
772  if (i != 0) str << ", ";
773  str << pnode->tensor_proper_size(i);
774  }
775  str << ")";
776  break;
777 
778  default: GMM_ASSERT1(false, "Invalid tensor dimension");
779  }
780  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
781  }
782 
783  void ga_print_node(const pga_tree_node pnode,
784  std::ostream &str) {
785  if (!pnode) return;
786  long prec = str.precision(16);
787 
788  bool is_interpolate(false), is_elementary(false), is_secondary(false);
789  bool is_xfem_plus(false), is_xfem_minus(false);
790  switch(pnode->node_type) {
791  case GA_NODE_INTERPOLATE:
792  case GA_NODE_INTERPOLATE_X:
793  case GA_NODE_INTERPOLATE_ELT_K: case GA_NODE_INTERPOLATE_ELT_B:
794  case GA_NODE_INTERPOLATE_NORMAL:
795  case GA_NODE_INTERPOLATE_VAL:
796  case GA_NODE_INTERPOLATE_GRAD:
797  case GA_NODE_INTERPOLATE_HESS:
798  case GA_NODE_INTERPOLATE_DIVERG:
799  case GA_NODE_INTERPOLATE_VAL_TEST:
800  case GA_NODE_INTERPOLATE_GRAD_TEST:
801  case GA_NODE_INTERPOLATE_HESS_TEST:
802  case GA_NODE_INTERPOLATE_DIVERG_TEST:
803  str << "Interpolate(";
804  is_interpolate = true;
805  break;
806  case GA_NODE_ELEMENTARY:
807  case GA_NODE_ELEMENTARY_VAL:
808  case GA_NODE_ELEMENTARY_GRAD:
809  case GA_NODE_ELEMENTARY_HESS:
810  case GA_NODE_ELEMENTARY_DIVERG:
811  case GA_NODE_ELEMENTARY_VAL_TEST:
812  case GA_NODE_ELEMENTARY_GRAD_TEST:
813  case GA_NODE_ELEMENTARY_HESS_TEST:
814  case GA_NODE_ELEMENTARY_DIVERG_TEST:
815  is_elementary = true;
816  str << "Elementary_transformation(";
817  break;
818  case GA_NODE_SECONDARY_DOMAIN:
819  case GA_NODE_SECONDARY_DOMAIN_X:
820  case GA_NODE_SECONDARY_DOMAIN_NORMAL:
821  case GA_NODE_SECONDARY_DOMAIN_VAL:
822  case GA_NODE_SECONDARY_DOMAIN_GRAD:
823  case GA_NODE_SECONDARY_DOMAIN_HESS:
824  case GA_NODE_SECONDARY_DOMAIN_DIVERG:
825  case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
826  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
827  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
828  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
829  str << "Secondary_domain(";
830  is_secondary = true;
831  break;
832 
833  case GA_NODE_XFEM_PLUS:
834  case GA_NODE_XFEM_PLUS_VAL:
835  case GA_NODE_XFEM_PLUS_GRAD:
836  case GA_NODE_XFEM_PLUS_HESS:
837  case GA_NODE_XFEM_PLUS_DIVERG:
838  case GA_NODE_XFEM_PLUS_VAL_TEST:
839  case GA_NODE_XFEM_PLUS_GRAD_TEST:
840  case GA_NODE_XFEM_PLUS_HESS_TEST:
841  case GA_NODE_XFEM_PLUS_DIVERG_TEST:
842  is_xfem_plus = true;
843  str << "Xfem_plus(";
844  break;
845  case GA_NODE_XFEM_MINUS:
846  case GA_NODE_XFEM_MINUS_VAL:
847  case GA_NODE_XFEM_MINUS_GRAD:
848  case GA_NODE_XFEM_MINUS_HESS:
849  case GA_NODE_XFEM_MINUS_DIVERG:
850  case GA_NODE_XFEM_MINUS_VAL_TEST:
851  case GA_NODE_XFEM_MINUS_GRAD_TEST:
852  case GA_NODE_XFEM_MINUS_HESS_TEST:
853  case GA_NODE_XFEM_MINUS_DIVERG_TEST:
854  is_xfem_minus = true;
855  str << "Xfem_minus(";
856  break;
857  default:
858  break;
859  }
860 
861  switch(pnode->node_type) {
862  case GA_NODE_GRAD:
863  case GA_NODE_INTERPOLATE_GRAD:
864  case GA_NODE_ELEMENTARY_GRAD:
865  case GA_NODE_SECONDARY_DOMAIN_GRAD:
866  case GA_NODE_XFEM_PLUS_GRAD:
867  case GA_NODE_XFEM_MINUS_GRAD:
868  case GA_NODE_GRAD_TEST:
869  case GA_NODE_INTERPOLATE_GRAD_TEST:
870  case GA_NODE_ELEMENTARY_GRAD_TEST:
871  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
872  case GA_NODE_XFEM_PLUS_GRAD_TEST:
873  case GA_NODE_XFEM_MINUS_GRAD_TEST:
874  str << "Grad_";
875  break;
876  case GA_NODE_HESS:
877  case GA_NODE_INTERPOLATE_HESS:
878  case GA_NODE_ELEMENTARY_HESS:
879  case GA_NODE_SECONDARY_DOMAIN_HESS:
880  case GA_NODE_XFEM_PLUS_HESS:
881  case GA_NODE_XFEM_MINUS_HESS:
882  case GA_NODE_HESS_TEST:
883  case GA_NODE_INTERPOLATE_HESS_TEST:
884  case GA_NODE_ELEMENTARY_HESS_TEST:
885  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
886  case GA_NODE_XFEM_PLUS_HESS_TEST:
887  case GA_NODE_XFEM_MINUS_HESS_TEST:
888  str << "Hess_";
889  break;
890  case GA_NODE_DIVERG:
891  case GA_NODE_INTERPOLATE_DIVERG:
892  case GA_NODE_SECONDARY_DOMAIN_DIVERG:
893  case GA_NODE_ELEMENTARY_DIVERG:
894  case GA_NODE_XFEM_PLUS_DIVERG:
895  case GA_NODE_XFEM_MINUS_DIVERG:
896  case GA_NODE_DIVERG_TEST:
897  case GA_NODE_INTERPOLATE_DIVERG_TEST:
898  case GA_NODE_ELEMENTARY_DIVERG_TEST:
899  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
900  case GA_NODE_XFEM_PLUS_DIVERG_TEST:
901  case GA_NODE_XFEM_MINUS_DIVERG_TEST:
902  str << "Div_";
903  break;
904  default:
905  break;
906  }
907 
908  switch(pnode->node_type) {
909  case GA_NODE_OP:
910  {
911  bool par = false;
912  if (pnode->parent) {
913  if (pnode->parent->node_type == GA_NODE_OP &&
914  (ga_operator_priorities[pnode->op_type] >= 2 ||
915  ga_operator_priorities[pnode->op_type]
916  < ga_operator_priorities[pnode->parent->op_type]))
917  par = true;
918  if (pnode->parent->node_type == GA_NODE_PARAMS) par = true;
919  }
920 
921 
922  if (par) str << "(";
923  if (pnode->op_type == GA_UNARY_MINUS) {
924  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
925  str << "-"; ga_print_node(pnode->children[0], str);
926  } else if (pnode->op_type == GA_QUOTE) {
927  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
928  ga_print_node(pnode->children[0], str); str << "'";
929  } else if (pnode->op_type == GA_SYM) {
930  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
931  str << "Sym("; ga_print_node(pnode->children[0], str); str << ")";
932  } else if (pnode->op_type == GA_SKEW) {
933  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
934  str << "Skew("; ga_print_node(pnode->children[0], str); str << ")";
935  } else if (pnode->op_type == GA_TRACE) {
936  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
937  str << "Trace("; ga_print_node(pnode->children[0], str); str << ")";
938  } else if (pnode->op_type == GA_DEVIATOR) {
939  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree with "
940  << pnode->children.size() << " children instead of 1");
941  str << "Deviator("; ga_print_node(pnode->children[0], str); str<<")";
942  } else if (pnode->op_type == GA_PRINT) {
943  GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
944  str << "Print("; ga_print_node(pnode->children[0], str); str << ")";
945  } else {
946  if (!par && pnode->op_type == GA_MULT &&
947  (pnode->children.size() == 1 ||
948  pnode->test_function_type == size_type(-1) ||
949  (pnode->children[0]->tensor_order() == 4 &&
950  pnode->children[1]->tensor_order() == 2)))
951  { par = true; str << "("; }
952  ga_print_node(pnode->children[0], str);
953  switch (pnode->op_type) {
954  case GA_PLUS: str << "+"; break;
955  case GA_MINUS: str << "-"; break;
956  case GA_MULT: str << "*"; break;
957  case GA_DIV: str << "/"; break;
958  case GA_COLON: str << ":"; break;
959  case GA_DOT: str << "."; break;
960  case GA_DOTMULT: str << ".*"; break;
961  case GA_DOTDIV: str << "./"; break;
962  case GA_TMULT: str << "@"; break;
963  default: GMM_ASSERT1(false, "Invalid or not taken into account "
964  "operation");
965  }
966  if (pnode->children.size() >= 2)
967  ga_print_node(pnode->children[1], str);
968  else
969  str << "(unknown second argument)";
970  }
971  if (par) str << ")";
972  }
973  break;
974 
975  case GA_NODE_X:
976  if (pnode->nbc1) str << "X(" << pnode->nbc1 << ")"; else str << "X";
977  break;
978  case GA_NODE_ELT_SIZE: str << "element_size"; break;
979  case GA_NODE_ELT_K: str << "element_K"; break;
980  case GA_NODE_ELT_B: str << "element_B"; break;
981  case GA_NODE_NORMAL: str << "Normal"; break;
982  case GA_NODE_INTERPOLATE_FILTER:
983  str << "Interpolate_filter(" << pnode->interpolate_name << ",";
984  ga_print_node(pnode->children[0], str);
985  if (pnode->children.size() == 2)
986  { str << ","; ga_print_node(pnode->children[1], str); }
987  else if (pnode->nbc1 != size_type(-1)) str << "," << pnode->nbc1;
988  str << ")";
989  break;
990  case GA_NODE_INTERPOLATE_X: case GA_NODE_SECONDARY_DOMAIN_X:
991  str << "X";
992  break;
993  case GA_NODE_INTERPOLATE_NORMAL: case GA_NODE_SECONDARY_DOMAIN_NORMAL:
994  str << "Normal";
995  break;
996  case GA_NODE_INTERPOLATE_ELT_K:
997  str << "element_K";
998  break;
999  case GA_NODE_INTERPOLATE_ELT_B:
1000  str << "element_B";
1001  break;
1002  case GA_NODE_INTERPOLATE_DERIVATIVE:
1003  str << (pnode->test_function_type == 1 ? "Test_" : "Test2_")
1004  << "Interpolate_derivative(" << pnode->interpolate_name_der << ","
1005  << pnode->name;
1006  if (pnode->interpolate_name.size())
1007  str << "," << pnode->interpolate_name;
1008  str << ")";
1009  break;
1010  case GA_NODE_INTERPOLATE:
1011  case GA_NODE_ELEMENTARY:
1012  case GA_NODE_SECONDARY_DOMAIN:
1013  case GA_NODE_XFEM_PLUS:
1014  case GA_NODE_XFEM_MINUS:
1015  case GA_NODE_VAL:
1016  case GA_NODE_INTERPOLATE_VAL:
1017  case GA_NODE_ELEMENTARY_VAL:
1018  case GA_NODE_SECONDARY_DOMAIN_VAL:
1019  case GA_NODE_XFEM_PLUS_VAL:
1020  case GA_NODE_XFEM_MINUS_VAL:
1021  case GA_NODE_GRAD:
1022  case GA_NODE_INTERPOLATE_GRAD:
1023  case GA_NODE_SECONDARY_DOMAIN_GRAD:
1024  case GA_NODE_ELEMENTARY_GRAD:
1025  case GA_NODE_XFEM_PLUS_GRAD:
1026  case GA_NODE_XFEM_MINUS_GRAD:
1027  case GA_NODE_HESS:
1028  case GA_NODE_INTERPOLATE_HESS:
1029  case GA_NODE_SECONDARY_DOMAIN_HESS:
1030  case GA_NODE_ELEMENTARY_HESS:
1031  case GA_NODE_XFEM_PLUS_HESS:
1032  case GA_NODE_XFEM_MINUS_HESS:
1033  case GA_NODE_DIVERG:
1034  case GA_NODE_INTERPOLATE_DIVERG:
1035  case GA_NODE_ELEMENTARY_DIVERG:
1036  case GA_NODE_SECONDARY_DOMAIN_DIVERG:
1037  case GA_NODE_XFEM_PLUS_DIVERG:
1038  case GA_NODE_XFEM_MINUS_DIVERG:
1039  str << pnode->name;
1040  break;
1041  case GA_NODE_VAL_TEST:
1042  case GA_NODE_INTERPOLATE_VAL_TEST:
1043  case GA_NODE_ELEMENTARY_VAL_TEST:
1044  case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
1045  case GA_NODE_XFEM_PLUS_VAL_TEST:
1046  case GA_NODE_XFEM_MINUS_VAL_TEST:
1047  case GA_NODE_GRAD_TEST:
1048  case GA_NODE_INTERPOLATE_GRAD_TEST:
1049  case GA_NODE_ELEMENTARY_GRAD_TEST:
1050  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
1051  case GA_NODE_XFEM_PLUS_GRAD_TEST:
1052  case GA_NODE_XFEM_MINUS_GRAD_TEST:
1053  case GA_NODE_HESS_TEST:
1054  case GA_NODE_INTERPOLATE_HESS_TEST:
1055  case GA_NODE_ELEMENTARY_HESS_TEST:
1056  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
1057  case GA_NODE_XFEM_PLUS_HESS_TEST:
1058  case GA_NODE_XFEM_MINUS_HESS_TEST:
1059  case GA_NODE_DIVERG_TEST:
1060  case GA_NODE_INTERPOLATE_DIVERG_TEST:
1061  case GA_NODE_ELEMENTARY_DIVERG_TEST:
1062  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
1063  case GA_NODE_XFEM_PLUS_DIVERG_TEST:
1064  case GA_NODE_XFEM_MINUS_DIVERG_TEST:
1065  str << (pnode->test_function_type == 1 ? "Test_" : "Test2_")
1066  << pnode->name;
1067  break;
1068  case GA_NODE_SPEC_FUNC: str << pnode->name; break;
1069  case GA_NODE_OPERATOR:
1070  case GA_NODE_PREDEF_FUNC:
1071  if (pnode->der1) {
1072  str << "Derivative_" << pnode->der1 << "_";
1073  if (pnode->der2) str << pnode->der2 << "_";
1074  }
1075  str << pnode->name; break;
1076  case GA_NODE_ZERO:
1077  GMM_ASSERT1(pnode->test_function_type != size_type(-1),
1078  "Internal error");
1079  if (pnode->test_function_type) str << "(";
1080  ga_print_constant_tensor(pnode, str);
1081  if (pnode->name_test1.size()) {
1082  GMM_ASSERT1(pnode->qdim1 > 0, "Internal error");
1083  if (pnode->qdim1 == 1)
1084  str << "*Test_" << pnode->name_test1;
1085  else {
1086  str << "*(Reshape(Test_" << pnode->name_test1 << ","
1087  << pnode->qdim1<< ")(1))";
1088  }
1089  }
1090  if (pnode->name_test2.size()) {
1091  GMM_ASSERT1(pnode->qdim2 > 0, "Internal error");
1092  if (pnode->qdim2 == 1)
1093  str << "*Test2_" << pnode->name_test2;
1094  else {
1095  str << "*(Reshape(Test2_" << pnode->name_test2 << ","
1096  << pnode->qdim2<< ")(1))";
1097  }
1098  }
1099  if (pnode->test_function_type) str << ")";
1100  break;
1101 
1102  case GA_NODE_CONSTANT:
1103  ga_print_constant_tensor(pnode, str);
1104  break;
1105 
1106  case GA_NODE_ALLINDICES:
1107  str << ":";
1108  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1109  break;
1110 
1111  case GA_NODE_PARAMS:
1112  GMM_ASSERT1(pnode->children.size(), "Invalid tree");
1113  ga_print_node(pnode->children[0], str);
1114  str << "(";
1115  for (size_type i = 1; i < pnode->children.size(); ++i) {
1116  if (i > 1) str << ", ";
1117  ga_print_node(pnode->children[i], str);
1118  }
1119  str << ")";
1120  break;
1121 
1122  case GA_NODE_NAME:
1123  str << pnode->name;
1124  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1125  break;
1126 
1127  case GA_NODE_MACRO_PARAM:
1128  if (pnode->nbc2 == 1) str << "Grad_";
1129  if (pnode->nbc2 == 2) str << "Hess_";
1130  if (pnode->nbc2 == 3) str << "Div_";
1131  if (pnode->nbc3 == 1) str << "Test_";
1132  if (pnode->nbc3 == 2) str << "Test2_";
1133  str << "P" << pnode->nbc1;
1134  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1135  break;
1136 
1137  case GA_NODE_RESHAPE:
1138  str << "Reshape";
1139  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1140  break;
1141 
1142  case GA_NODE_CROSS_PRODUCT:
1143  str << "Cross_product";
1144  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1145  break;
1146 
1147  case GA_NODE_SWAP_IND:
1148  str << "Swap_indices";
1149  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1150  break;
1151 
1152  case GA_NODE_IND_MOVE_LAST:
1153  str << "Index_move_last";
1154  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1155  break;
1156 
1157  case GA_NODE_CONTRACT:
1158  str << "Contract";
1159  GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
1160  break;
1161 
1162  case GA_NODE_C_MATRIX:
1163  GMM_ASSERT1(pnode->children.size(), "Invalid tree");
1164  GMM_ASSERT1(pnode->nbc1 == pnode->tensor_order(), "Invalid C_MATRIX");
1165  switch (pnode->tensor_order()) {
1166  case 0:
1167  ga_print_node(pnode->children[0], str);
1168  break;
1169 
1170  case 1:
1171  str << "[";
1172  for (size_type i = 0; i < pnode->tensor_proper_size(0); ++i) {
1173  if (i != 0) str << ",";
1174  ga_print_node(pnode->children[i], str);
1175  }
1176  str << "]";
1177  break;
1178 
1179  case 2: case 3: case 4:
1180  {
1181  size_type ii(0);
1182  size_type n0 = pnode->tensor_proper_size(0);
1183  size_type n1 = pnode->tensor_proper_size(1);
1184  size_type n2 = ((pnode->tensor_order() > 2) ?
1185  pnode->tensor_proper_size(2) : 1);
1186  size_type n3 = ((pnode->tensor_order() > 3) ?
1187  pnode->tensor_proper_size(3) : 1);
1188  if (n3 > 1) str << "[";
1189  for (size_type l = 0; l < n3; ++l) {
1190  if (l != 0) str << ",";
1191  if (n2 > 1) str << "[";
1192  for (size_type k = 0; k < n2; ++k) {
1193  if (k != 0) str << ",";
1194  if (n1 > 1) str << "[";
1195  for (size_type j = 0; j < n1; ++j) {
1196  if (j != 0) str << ",";
1197  if (n0 > 1) str << "[";
1198  for (size_type i = 0; i < n0; ++i) {
1199  if (i != 0) str << ",";
1200  ga_print_node(pnode->children[ii++], str);
1201  }
1202  if (n0 > 1) str << "]";
1203  }
1204  if (n1 > 1) str << "]";
1205  }
1206  if (n2 > 1) str << "]";
1207  }
1208  if (n3 > 1) str << "]";
1209  }
1210  break;
1211 
1212  case 5: case 6: case 7: case 8:
1213  str << "Reshape([";
1214  for (size_type i = 0; i < pnode->tensor_proper_size(); ++i) {
1215  if (i != 0) str << ";";
1216  ga_print_node(pnode->children[i], str);
1217  }
1218  str << "]";
1219  for (size_type i = 0; i < pnode->tensor_order(); ++i) {
1220  if (i != 0) str << ", ";
1221  str << pnode->tensor_proper_size(i);
1222  }
1223  str << ")";
1224  break;
1225 
1226  default: GMM_ASSERT1(false, "Invalid tensor dimension");
1227  }
1228  break;
1229 
1230  default:
1231  str << "Invalid or not taken into account node type "
1232  << pnode->node_type;
1233  break;
1234  }
1235 
1236  if (is_interpolate)
1237  str << "," << pnode->interpolate_name << ")";
1238  else if (is_elementary) {
1239  str << "," << pnode->elementary_name;
1240  if (pnode->name.compare(pnode->elementary_target) != 0)
1241  str << "," << pnode->elementary_target;
1242  str << ")";
1243  } else if (is_secondary)
1244  str << ")";
1245  else if (is_xfem_plus || is_xfem_minus)
1246  str << ")";
1247 
1248  str.precision(prec);
1249  }
1250 
1251  std::string ga_tree_to_string(const ga_tree &tree) {
1252  std::stringstream str;
1253  str.precision(16);
1254  if (tree.root) verify_tree(tree.root, 0);
1255  if (tree.root) ga_print_node(tree.root, str); else str << "0";
1256  return str.str();
1257  }
1258 
1259  size_type ga_parse_prefix_operator(std::string &name) {
1260  if (name.size() >= 5 && name.compare(0, 5, "Grad_") == 0)
1261  { name = name.substr(5); return 1; }
1262  else if (name.size() >= 5 && name.compare(0, 5, "Hess_") == 0)
1263  { name = name.substr(5); return 2; }
1264  else if (name.size() >= 4 && name.compare(0, 4, "Div_") == 0)
1265  { name = name.substr(4); return 3; }
1266  return 0;
1267  }
1268 
1269  size_type ga_parse_prefix_test(std::string &name) {
1270  if (name.size() >= 5 && name.compare(0, 5, "Test_") == 0)
1271  { name = name.substr(5); return 1; }
1272  else if (name.size() >= 6 && name.compare(0, 6, "Test2_") == 0)
1273  { name = name.substr(6); return 2; }
1274  return 0;
1275  }
1276 
1277  // 0 : ok
1278  // 1 : function or operator name or "X"
1279  // 2 : reserved prefix Grad, Hess, Div, Derivative_ Test and Test2
1280  // 3 : reserved prefix Dot and Previous
1281  int ga_check_name_validity(const std::string &name) {
1282  if (name.compare(0, 11, "Derivative_") == 0)
1283  return 2;
1284 
1285  const ga_predef_operator_tab &PREDEF_OPERATORS
1286  = dal::singleton<ga_predef_operator_tab>::instance(0);
1287  const ga_spec_function_tab &SPEC_FUNCTIONS
1288  = dal::singleton<ga_spec_function_tab>::instance(0);
1289  const ga_spec_op_tab &SPEC_OP
1290  = dal::singleton<ga_spec_op_tab>::instance(0);
1291  const ga_predef_function_tab &PREDEF_FUNCTIONS
1292  = dal::singleton<ga_predef_function_tab>::instance(0);
1293 
1294  if (SPEC_OP.find(name) != SPEC_OP.end())
1295  return 1;
1296 
1297  if (PREDEF_FUNCTIONS.find(name) != PREDEF_FUNCTIONS.end())
1298  return 1;
1299 
1300  if (SPEC_FUNCTIONS.find(name) != SPEC_FUNCTIONS.end())
1301  return 1;
1302 
1303  if (PREDEF_OPERATORS.tab.find(name) != PREDEF_OPERATORS.tab.end())
1304  return 1;
1305 
1306  if (name.size() >= 5 && name.compare(0, 5, "Grad_") == 0)
1307  return 2;
1308 
1309  if (name.size() >= 5 && name.compare(0, 5, "Hess_") == 0)
1310  return 2;
1311 
1312  if (name.size() >= 4 && name.compare(0, 4, "Div_") == 0)
1313  return 2;
1314 
1315  if (name.size() >= 6 && name.compare(0, 6, "Test2_") == 0)
1316  return 2;
1317 
1318  if (name.size() >= 5 && name.compare(0, 5, "Test_") == 0)
1319  return 2;
1320 
1321 // if (name.size() >= 4 && name.compare(0, 4, "Dot_") == 0)
1322 // return 3;
1323 // if (name.size() >= 5 && name.compare(0, 5, "Dot2_") == 0)
1324 // return 3;
1325 
1326 // if (name.size() >= 9 && name.compare(0, 9, "Previous_") == 0)
1327 // return 3;
1328 // if (name.size() >= 10 && name.compare(0, 10, "Previous2_") == 0)
1329 // return 3;
1330 // if (name.size() >= 12 && name.compare(0, 12, "Previous1_2_") == 0)
1331 // return 3;
1332 
1333  return 0;
1334  }
1335 
1336  //=========================================================================
1337  // Structure dealing with macros.
1338  //=========================================================================
1339 
1340  ga_macro::ga_macro() : ptree(new ga_tree), nbp(0) {}
1341  ga_macro::~ga_macro() { delete ptree; }
1342  ga_macro::ga_macro(const std::string &name, const ga_tree &t, size_type nbp_)
1343  : ptree(new ga_tree(t)), macro_name_(name), nbp(nbp_) {}
1344  ga_macro::ga_macro(const ga_macro &gam)
1345  : ptree(new ga_tree(gam.tree())), macro_name_(gam.name()),
1346  nbp(gam.nb_params()) {}
1347  ga_macro &ga_macro::operator =(const ga_macro &gam) {
1348  delete ptree; ptree = new ga_tree(gam.tree());
1349  macro_name_ = gam.name();
1350  nbp = gam.nb_params();
1351  return *this;
1352  }
1353 
1354  static void ga_replace_macro_params
1355  (ga_tree &tree, pga_tree_node pnode,
1356  const std::vector<pga_tree_node> &children) {
1357  if (!pnode) return;
1358  for (size_type i = 0; i < pnode->children.size(); ++i)
1359  ga_replace_macro_params(tree, pnode->children[i], children);
1360 
1361  if (pnode->node_type == GA_NODE_MACRO_PARAM) {
1362  size_type po = pnode->nbc2;
1363  size_type pt = pnode->nbc3;
1364  GMM_ASSERT1(pnode->nbc1+1 < children.size(), "Internal error");
1365  pga_tree_node pchild = children[pnode->nbc1+1];
1366 
1367  if (po || pt || pnode->op_type != GA_NAME) {
1368  if (!(pchild->children.empty()) || pchild->node_type != GA_NODE_NAME)
1369  ga_throw_error(pchild->expr, pchild->pos, "Error in macro "
1370  "expansion. Only variable name are allowed for macro "
1371  "parameter preceded by Grad_ Hess_ Test_ or Test2_ "
1372  "prefixes.");
1373  switch(pnode->op_type) {
1374  case GA_NAME : pnode->node_type = GA_NODE_NAME; break;
1375  case GA_INTERPOLATE : pnode->node_type = GA_NODE_INTERPOLATE; break;
1376  case GA_INTERPOLATE_DERIVATIVE :
1377  pnode->node_type = GA_NODE_INTERPOLATE_DERIVATIVE; break;
1378  case GA_ELEMENTARY : pnode->node_type = GA_NODE_ELEMENTARY; break;
1379  case GA_SECONDARY_DOMAIN :
1380  pnode->node_type = GA_NODE_SECONDARY_DOMAIN; break;
1381  case GA_XFEM_PLUS : pnode->node_type = GA_NODE_XFEM_PLUS; break;
1382  case GA_XFEM_MINUS: pnode->node_type = GA_NODE_XFEM_MINUS; break;
1383  default:break;
1384  }
1385  pnode->name = pchild->name;
1386  if (pt == 1) pnode->name = "Test_" + pnode->name;
1387  if (pt == 2) pnode->name = "Test2_" + pnode->name;
1388  if (po == 1) pnode->name = "Grad_" + pnode->name;
1389  if (po == 2) pnode->name = "Hess_" + pnode->name;
1390  if (po == 3) pnode->name = "Div_" + pnode->name;
1391  } else {
1392  pga_tree_node pnode_old = pnode;
1393  pnode = nullptr;
1394  tree.copy_node(pchild, pnode_old->parent, pnode);
1395  if (pnode_old->parent)
1396  pnode_old->parent->replace_child(pnode_old, pnode);
1397  else
1398  tree.root = pnode;
1399  GMM_ASSERT1(pnode_old->children.empty(), "Internal error");
1400  delete pnode_old;
1401  }
1402  }
1403  }
1404 
1405  static void ga_expand_macro(ga_tree &tree, pga_tree_node pnode,
1406  const ga_macro_dictionary &macro_dict) {
1407  if (!pnode) return;
1408 
1409  if (pnode->node_type == GA_NODE_PARAMS) {
1410 
1411  for (size_type i = 1; i < pnode->children.size(); ++i)
1412  ga_expand_macro(tree, pnode->children[i], macro_dict);
1413 
1414  if (pnode->children[0]->node_type != GA_NODE_NAME) {
1415  ga_expand_macro(tree, pnode->children[0], macro_dict);
1416  } else {
1417 
1418  if (macro_dict.macro_exists(pnode->children[0]->name)) {
1419 
1420  const ga_macro &gam = macro_dict.get_macro(pnode->children[0]->name);
1421 
1422  if (gam.nb_params()==0) { // Macro without parameters
1423  pga_tree_node pnode_old = pnode->children[0];
1424  pnode->children[0] = nullptr;
1425  tree.copy_node(gam.tree().root,
1426  pnode_old->parent,pnode->children[0]);
1427  GMM_ASSERT1(pnode_old->children.empty(), "Internal error");
1428  delete pnode_old;
1429  } else { // Macro with parameters
1430  if (gam.nb_params()+1 != pnode->children.size())
1431  ga_throw_error(pnode->expr, pnode->pos,
1432  "Bad number of parameters in the use of macro '"
1433  << gam.name() << "'. Expected " << gam.nb_params()
1434  << " found " << pnode->children.size()-1 << ".");
1435 
1436  pga_tree_node pnode_old = pnode;
1437  pnode = nullptr;
1438  tree.copy_node(gam.tree().root, pnode_old->parent, pnode);
1439  if (pnode_old->parent)
1440  pnode_old->parent->replace_child(pnode_old, pnode);
1441  else
1442  tree.root = pnode;
1443  ga_replace_macro_params(tree, pnode, pnode_old->children);
1444  tree.clear_node_rec(pnode_old);
1445  }
1446  }
1447  }
1448 
1449  } else if (pnode->node_type == GA_NODE_NAME &&
1450  macro_dict.macro_exists(pnode->name)) {
1451  // Macro without parameters
1452  const ga_macro &gam = macro_dict.get_macro(pnode->name);
1453  if (gam.nb_params() != 0)
1454  ga_throw_error(pnode->expr, pnode->pos,
1455  "Bad number of parameters in the use of macro '"
1456  << gam.name() << "'. Expected " << gam.nb_params()
1457  << " none found.");
1458 
1459  pga_tree_node pnode_old = pnode;
1460  pnode = nullptr;
1461  tree.copy_node(gam.tree().root, pnode_old->parent, pnode);
1462  if (pnode_old->parent)
1463  pnode_old->parent->replace_child(pnode_old, pnode);
1464  else
1465  tree.root = pnode;
1466  GMM_ASSERT1(pnode_old->children.empty(), "Internal error");
1467  delete pnode_old;
1468  } else {
1469  for (size_type i = 0; i < pnode->children.size(); ++i)
1470  ga_expand_macro(tree, pnode->children[i], macro_dict);
1471  }
1472  }
1473 
1474  static void ga_mark_macro_params_rec(const pga_tree_node pnode,
1475  const std::vector<std::string> &params) {
1476  if (!pnode) return;
1477  for (size_type i = 0; i < pnode->children.size(); ++i)
1478  ga_mark_macro_params_rec(pnode->children[i], params);
1479 
1480  if (pnode->node_type == GA_NODE_NAME ||
1481  pnode->node_type == GA_NODE_INTERPOLATE ||
1482  pnode->node_type == GA_NODE_ELEMENTARY ||
1483  pnode->node_type == GA_NODE_SECONDARY_DOMAIN ||
1484  pnode->node_type == GA_NODE_XFEM_PLUS ||
1485  pnode->node_type == GA_NODE_XFEM_MINUS) {
1486  std::string name = pnode->name;
1487  size_type po = ga_parse_prefix_operator(name);
1488  size_type pt = ga_parse_prefix_test(name);
1489 
1490  for (size_type i = 0; i < params.size(); ++i)
1491  if (name.compare(params[i]) == 0) {
1492  pnode->name = name;
1493  switch(pnode->node_type) {
1494  case GA_NODE_NAME : pnode->op_type = GA_NAME; break;
1495  case GA_NODE_INTERPOLATE : pnode->op_type = GA_INTERPOLATE; break;
1496  case GA_NODE_INTERPOLATE_DERIVATIVE :
1497  pnode->op_type = GA_INTERPOLATE_DERIVATIVE; break;
1498  case GA_NODE_ELEMENTARY : pnode->op_type = GA_ELEMENTARY; break;
1499  case GA_NODE_SECONDARY_DOMAIN :
1500  pnode->op_type = GA_SECONDARY_DOMAIN; break;
1501  case GA_NODE_XFEM_PLUS : pnode->op_type = GA_XFEM_PLUS; break;
1502  case GA_NODE_XFEM_MINUS: pnode->op_type = GA_XFEM_MINUS; break;
1503  default:break;
1504  }
1505  pnode->node_type = GA_NODE_MACRO_PARAM;
1506  pnode->nbc1 = i; pnode->nbc2 = po; pnode->nbc3 = pt;
1507  }
1508  }
1509  }
1510 
1511  static void ga_mark_macro_params(ga_macro &gam,
1512  const std::vector<std::string> &params,
1513  const ga_macro_dictionary &macro_dict) {
1514  if (gam.tree().root) {
1515  ga_mark_macro_params_rec(gam.tree().root, params);
1516  ga_expand_macro(gam.tree(), gam.tree().root, macro_dict);
1517  }
1518  }
1519 
1520  bool ga_macro_dictionary::macro_exists(const std::string &name) const {
1521  if (macros.find(name) != macros.end()) return true;
1522  if (parent && parent->macro_exists(name)) return true;
1523  return false;
1524  }
1525 
1526  const ga_macro &
1527  ga_macro_dictionary::get_macro(const std::string &name) const {
1528  auto it = macros.find(name);
1529  if (it != macros.end()) return it->second;
1530  if (parent) return parent->get_macro(name);
1531  GMM_ASSERT1(false, "Undefined macro");
1532  }
1533 
1534  void ga_macro_dictionary::add_macro(const ga_macro &gam)
1535  { macros[gam.name()] = gam; }
1536 
1537  void ga_macro_dictionary::add_macro(const std::string &name,
1538  const std::string &expr)
1539  { ga_tree tree; ga_read_string_reg("Def "+name+":="+expr, tree, *this); }
1540 
1541  void ga_macro_dictionary::del_macro(const std::string &name) {
1542  auto it = macros.find(name);
1543  GMM_ASSERT1(it != macros.end(), "Undefined macro (at this level)");
1544  macros.erase(it);
1545  }
1546 
1547 
1548  //=========================================================================
1549  // Syntax analysis for the generic assembly language
1550  //=========================================================================
1551 
1552  // Read a term with an (implicit) pushdown automaton.
1553  static GA_TOKEN_TYPE ga_read_term(pstring expr, size_type &pos,
1554  ga_tree &tree,
1555  ga_macro_dictionary &macro_dict) {
1556  size_type token_pos, token_length;
1557  GA_TOKEN_TYPE t_type;
1558  int state = 1; // 1 = reading term, 2 = reading after term
1559 
1560  for (;;) {
1561 
1562  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1563 
1564  switch (state) {
1565 
1566  case 1:
1567  switch (t_type) {
1568  case GA_SCALAR:
1569  {
1570  char *endptr; const char *nptr = &((*expr)[token_pos]);
1571  scalar_type s_read = ::strtod(nptr, &endptr);
1572  if (endptr == nptr)
1573  ga_throw_error(expr, token_pos, "Bad numeric format.");
1574  tree.add_scalar(s_read, token_pos, expr);
1575  }
1576  state = 2; break;
1577 
1578  case GA_COLON:
1579  tree.add_allindices(token_pos, expr);
1580  state = 2; break;
1581 
1582  case GA_NAME:
1583  tree.add_name(&((*expr)[token_pos]), token_length, token_pos, expr);
1584  state = 2; break;
1585 
1586  case GA_MINUS: // unary -
1587  tree.add_op(GA_UNARY_MINUS, token_pos, expr);
1588  state = 1; break;
1589 
1590  case GA_PLUS: // unary +
1591  state = 1; break;
1592 
1593  case GA_SYM:
1594  tree.add_op(GA_SYM, token_pos, expr);
1595  state = 1; break;
1596 
1597  case GA_SKEW:
1598  tree.add_op(GA_SKEW, token_pos, expr);
1599  state = 1; break;
1600 
1601  case GA_TRACE:
1602  tree.add_op(GA_TRACE, token_pos, expr);
1603  state = 1; break;
1604 
1605  case GA_DEVIATOR:
1606  tree.add_op(GA_DEVIATOR, token_pos, expr);
1607  state = 1; break;
1608 
1609  case GA_DEF:
1610  {
1611  ga_macro gam;
1612  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1613  if (t_type != GA_NAME)
1614  ga_throw_error(expr, pos,
1615  "Macro definition should begin with macro name");
1616  gam.name() = std::string(&((*expr)[token_pos]), token_length);
1617  if (ga_check_name_validity(gam.name()))
1618  ga_throw_error(expr, pos-1, "Invalid macro name.")
1619  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1620  std::vector<std::string> params;
1621  if (t_type == GA_LPAR) {
1622  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1623  while (t_type == GA_NAME) {
1624  params.push_back(std::string(&((*expr)[token_pos]),
1625  token_length));
1626  if (ga_check_name_validity(params.back()))
1627  ga_throw_error(expr, pos-1, "Invalid macro parameter name.");
1628  for (size_type i = 0; i+1 < params.size(); ++i)
1629  if (params.back().compare(params[i]) == 0)
1630  ga_throw_error(expr, pos-1,
1631  "Invalid repeated macro parameter name.");
1632  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1633  if (t_type == GA_COMMA)
1634  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1635  }
1636  if (t_type != GA_RPAR)
1637  ga_throw_error(expr, pos-1,
1638  "Missing right parenthesis in macro definition.");
1639  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1640  }
1641  if (t_type != GA_COLON_EQ)
1642  ga_throw_error(expr, pos-1, "Missing := for macro definition.");
1643 
1644  t_type = ga_read_term(expr, pos, gam.tree(), macro_dict);
1645  if (gam.tree().root)
1646  ga_expand_macro(gam.tree(), gam.tree().root, macro_dict);
1647  gam.nb_params() = params.size();
1648  if (params.size())
1649  ga_mark_macro_params(gam, params, macro_dict);
1650  macro_dict.add_macro(gam);
1651 
1652  // cout << "macro \"" << gam.name() << "\" registered with "
1653  // << gam.nb_params() << " params := "
1654  // << ga_tree_to_string(gam.tree()) << endl;
1655 
1656  if (t_type == GA_END) return t_type;
1657  else if (t_type != GA_SEMICOLON)
1658  ga_throw_error(expr, pos-1,
1659  "Syntax error at the end of macro definition.");
1660  state = 1;
1661  }
1662  break;
1663 
1664  case GA_INTERPOLATE:
1665  {
1666  tree.add_scalar(scalar_type(0), token_pos, expr);
1667  tree.current_node->node_type = GA_NODE_INTERPOLATE;
1668  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1669  if (t_type != GA_LPAR)
1670  ga_throw_error(expr, pos-1, "Missing interpolate arguments.");
1671  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1672  if (t_type != GA_NAME)
1673  ga_throw_error(expr, pos,
1674  "First argument of Interpolate should be a "
1675  "variable, test function, X or Normal.");
1676  tree.current_node->name = std::string(&((*expr)[token_pos]),
1677  token_length);
1678 
1679  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1680  if (t_type != GA_COMMA)
1681  ga_throw_error(expr, pos, "Bad format for Interpolate "
1682  "arguments.");
1683  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1684  if (t_type != GA_NAME)
1685  ga_throw_error(expr, pos,
1686  "Second argument of Interpolate should be a "
1687  "transformation name.");
1688  tree.current_node->interpolate_name
1689  = std::string(&((*expr)[token_pos]), token_length);
1690  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1691  if (t_type != GA_RPAR)
1692  ga_throw_error(expr, pos-1, "Missing a parenthesis after "
1693  "interpolate arguments.");
1694  state = 2;
1695  }
1696  break;
1697 
1698  case GA_INTERPOLATE_DERIVATIVE:
1699  {
1700  tree.add_scalar(scalar_type(0), token_pos, expr);
1701  tree.current_node->node_type = GA_NODE_INTERPOLATE_DERIVATIVE;
1702  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1703  if (t_type != GA_LPAR)
1704  ga_throw_error(expr, pos-1,
1705  "Missing Interpolate_derivative arguments.");
1706  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1707  if (t_type != GA_NAME)
1708  ga_throw_error(expr, pos,
1709  "First argument of Interpolate should the "
1710  "interpolate transformtion name ");
1711  tree.current_node->interpolate_name_der
1712  = std::string(&((*expr)[token_pos]), token_length);
1713  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1714  if (t_type != GA_COMMA)
1715  ga_throw_error(expr, pos, "Bad format for Interpolate_derivative "
1716  "arguments.");
1717  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1718  if (t_type != GA_NAME)
1719  ga_throw_error(expr, pos,
1720  "Second argument of Interpolate should be a "
1721  "variable name.");
1722  tree.current_node->name
1723  = std::string(&((*expr)[token_pos]), token_length);
1724 
1725  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1726  tree.current_node->interpolate_name = "";
1727  if (t_type == GA_COMMA) {
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  "Third argument of Interpolate should be a "
1732  "interpolate 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  }
1737  if (t_type != GA_RPAR)
1738  ga_throw_error(expr, pos-1, "Missing a parenthesis after "
1739  "Interpolate_derivative arguments.");
1740  state = 2;
1741  }
1742  break;
1743 
1744  case GA_ELEMENTARY:
1745  {
1746  tree.add_scalar(scalar_type(0), token_pos, expr);
1747  tree.current_node->node_type = GA_NODE_ELEMENTARY;
1748  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1749  if (t_type != GA_LPAR)
1750  ga_throw_error(expr, pos-1,
1751  "Missing Elementary_transformation arguments.");
1752  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1753  if (t_type != GA_NAME)
1754  ga_throw_error(expr, pos,
1755  "First argument of Elementary_transformation "
1756  "should be a variable or a test function.");
1757  tree.current_node->name = std::string(&((*expr)[token_pos]),
1758  token_length);
1759  tree.current_node->elementary_target = tree.current_node->name;
1760 
1761  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1762  if (t_type != GA_COMMA)
1763  ga_throw_error(expr, pos, "Bad format for "
1764  "Elementary_transformation arguments.");
1765  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1766  if (t_type != GA_NAME)
1767  ga_throw_error(expr, pos,
1768  "Second argument of Elementary_transformation "
1769  "should be a transformation name.");
1770  tree.current_node->elementary_name
1771  = std::string(&((*expr)[token_pos]), token_length);
1772  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1773 
1774  if (t_type == GA_COMMA) {
1775  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1776  if (t_type != GA_NAME)
1777  ga_throw_error(expr, pos,
1778  "Third argument of Elementary_transformation "
1779  "should be a variable or data name.");
1780 
1781  tree.current_node->elementary_target =
1782  std::string(&((*expr)[token_pos]), token_length);
1783  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1784  }
1785 
1786  if (t_type != GA_RPAR)
1787  ga_throw_error(expr, pos-1, "Missing a parenthesis after "
1788  "Elementary_transformation arguments.");
1789  state = 2;
1790  }
1791  break;
1792 
1793  case GA_SECONDARY_DOMAIN:
1794  {
1795  tree.add_scalar(scalar_type(0), token_pos, expr);
1796  tree.current_node->node_type = GA_NODE_SECONDARY_DOMAIN;
1797  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1798  if (t_type != GA_LPAR)
1799  ga_throw_error(expr, pos-1,"Missing Secondary_domain arguments.");
1800  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1801  if (t_type != GA_NAME)
1802  ga_throw_error(expr, pos,
1803  "First argument of Secondary_domain should be a "
1804  "variable, test function, X or Normal.");
1805  tree.current_node->name = std::string(&((*expr)[token_pos]),
1806  token_length);
1807  tree.current_node->interpolate_name = tree.secondary_domain;
1808  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1809  if (t_type != GA_RPAR)
1810  ga_throw_error(expr, pos-1, "Missing a parenthesis after "
1811  "Secondary_domain arguments.");
1812  state = 2;
1813  }
1814  break;
1815 
1816  case GA_XFEM_PLUS:
1817  {
1818  tree.add_scalar(scalar_type(0), token_pos, expr);
1819  tree.current_node->node_type = GA_NODE_XFEM_PLUS;
1820  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1821  if (t_type != GA_LPAR)
1822  ga_throw_error(expr, pos-1,
1823  "Missing Xfem_plus arguments.");
1824  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1825  if (t_type != GA_NAME)
1826  ga_throw_error(expr, pos,
1827  "The argument of Xfem_plus should be a "
1828  "variable or a test function.");
1829  tree.current_node->name = std::string(&((*expr)[token_pos]),
1830  token_length);
1831  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1832  if (t_type != GA_RPAR)
1833  ga_throw_error(expr, pos-1, "Missing a parenthesis after "
1834  "Xfem_plus argument.");
1835  state = 2;
1836  }
1837  break;
1838 
1839  case GA_XFEM_MINUS:
1840  {
1841  tree.add_scalar(scalar_type(0), token_pos, expr);
1842  tree.current_node->node_type = GA_NODE_XFEM_MINUS;
1843  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1844  if (t_type != GA_LPAR)
1845  ga_throw_error(expr, pos-1,
1846  "Missing Xfem_minus arguments.");
1847  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1848  if (t_type != GA_NAME)
1849  ga_throw_error(expr, pos,
1850  "The argument of Xfem_minus should be a "
1851  "variable or a test function.");
1852  tree.current_node->name = std::string(&((*expr)[token_pos]),
1853  token_length);
1854  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1855  if (t_type != GA_RPAR)
1856  ga_throw_error(expr, pos-1, "Missing a parenthesis after "
1857  "Xfem_minus argument.");
1858  state = 2;
1859  }
1860  break;
1861 
1862  case GA_INTERPOLATE_FILTER:
1863  {
1864  tree.add_scalar(scalar_type(0), token_pos, expr);
1865  tree.current_node->node_type = GA_NODE_INTERPOLATE_FILTER;
1866  tree.current_node->nbc1 = size_type(-1);
1867  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1868  if (t_type != GA_LPAR)
1869  ga_throw_error(expr, pos-1, "Missing interpolate arguments.");
1870  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1871  if (t_type != GA_NAME)
1872  ga_throw_error(expr, pos, "First argument of Interpolate_filter "
1873  "should be a transformation name.");
1874  tree.current_node->interpolate_name
1875  = std::string(&((*expr)[token_pos]), token_length);
1876  t_type = ga_get_token(*expr, pos, token_pos, token_length);
1877  if (t_type != GA_COMMA)
1878  ga_throw_error(expr, pos,
1879  "Bad format for Interpolate_filter arguments.");
1880  ga_tree sub_tree;
1881  t_type = ga_read_term(expr, pos, sub_tree, macro_dict);
1882  if (t_type != GA_RPAR && t_type != GA_COMMA)
1883  ga_throw_error(expr, pos-1,
1884  "Bad format for Interpolate_filter arguments.");
1885  tree.add_sub_tree(sub_tree);
1886  if (t_type == GA_COMMA) {
1887  ga_tree sub_tree2;
1888  t_type = ga_read_term(expr, pos, sub_tree2, macro_dict);
1889  tree.add_sub_tree(sub_tree2);
1890  }
1891  if (t_type != GA_RPAR)
1892  ga_throw_error(expr, pos-1, "Unbalanced parenthesis.");
1893  state = 2;
1894  }
1895  break;
1896 
1897  case GA_PRINT:
1898  tree.add_op(GA_PRINT, token_pos, expr);
1899  state = 1; break;
1900 
1901  case GA_LPAR: // Parenthesed expression
1902  {
1903  ga_tree sub_tree;
1904  GA_TOKEN_TYPE r_type;
1905  r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
1906  if (r_type != GA_RPAR)
1907  ga_throw_error(expr, pos-1, "Unbalanced parenthesis.");
1908  tree.add_sub_tree(sub_tree);
1909  state = 2;
1910  }
1911  break;
1912 
1913  case GA_LBRACKET: // Explicit vector/matrix or tensor
1914  {
1915  ga_tree sub_tree;
1916  GA_TOKEN_TYPE r_type;
1917  size_type nbc1(0), nbc2(0), nbc3(0), n1(0), n2(0), n3(0);
1918  size_type tensor_order(1);
1919  bool foundcomma(false), foundsemi(false);
1920 
1921  r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
1922  size_type nb_comp = 0;
1923  tree.add_matrix(token_pos, expr);
1924 
1925  if (sub_tree.root->node_type == GA_NODE_C_MATRIX) { // nested format
1926  bgeot::multi_index mii;
1927  do {
1928  if (nb_comp) {
1929  sub_tree.clear();
1930  r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
1931  }
1932  // in the nested format only "," and "]" are expected
1933  if (sub_tree.root->node_type != GA_NODE_C_MATRIX ||
1934  (r_type != GA_COMMA && r_type != GA_RBRACKET))
1935  ga_throw_error(expr, pos-1, "Bad explicit "
1936  "vector/matrix/tensor format.");
1937 
1938  // convert a row vector [a,b] to a column vector [a;b]
1939  if (sub_tree.root->marked &&
1940  sub_tree.root->tensor().sizes()[0] == 1 &&
1941  sub_tree.root->tensor().size() != 1) {
1942  bgeot::multi_index mi = sub_tree.root->tensor().sizes();
1943  for (size_type i = mi.size()-1; i > 0; i--)
1944  mi[i-1] = mi[i];
1945  mi.pop_back();
1946  sub_tree.root->tensor().adjust_sizes(mi);
1947  }
1948  if (!nb_comp) mii = sub_tree.root->tensor().sizes();
1949  else {
1950  const bgeot::multi_index &mi=sub_tree.root->tensor().sizes();
1951  bool cmp = true;
1952  if (mii.size() == mi.size()) {
1953  for (size_type i = 0; i < mi.size(); ++i)
1954  if (mi[i] != mii[i]) cmp = false;
1955  } else cmp = false;
1956  if (!cmp)
1957  ga_throw_error(expr, pos-1, "Bad explicit "
1958  "vector/matrix/tensor format.");
1959  }
1960  for (size_type i = 0; i < sub_tree.root->children.size(); ++i) {
1961  sub_tree.root->children[i]->parent = tree.current_node;
1962  tree.current_node->children.push_back
1963  (sub_tree.root->children[i]);
1964  }
1965  sub_tree.root->children.resize(0);
1966  nb_comp++;
1967  } while (r_type != GA_RBRACKET);
1968  tree.current_node->marked = false;
1969  mii.push_back(nb_comp);
1970  tree.current_node->tensor().adjust_sizes(mii);
1971  } else { // non nested format
1972  do {
1973  if (nb_comp) {
1974  sub_tree.clear();
1975  r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
1976  }
1977  nb_comp++;
1978 
1979  tree.add_sub_tree(sub_tree);
1980 
1981  ++n1; ++n2; ++n3;
1982  if (tensor_order < 2) ++nbc1;
1983  if (tensor_order < 3) ++nbc2;
1984  if (tensor_order < 4) ++nbc3;
1985 
1986  if (r_type == GA_COMMA) {
1987  if (!foundcomma && tensor_order > 1)
1988  ga_throw_error(expr, pos-1, "Bad explicit "
1989  "vector/matrix/tensor format.");
1990  foundcomma = true;
1991  } else if (r_type == GA_SEMICOLON) {
1992  if (n1 != nbc1)
1993  ga_throw_error(expr, pos-1, "Bad explicit "
1994  "vector/matrix/tensor format.");
1995  n1 = 0;
1996  tensor_order = std::max(tensor_order, size_type(2));
1997  } else if (r_type == GA_DCOMMA) {
1998  if (n1 != nbc1 || n2 != nbc2)
1999  ga_throw_error(expr, pos-1, "Bad explicit "
2000  "vector/matrix/tensor format.");
2001  foundsemi = true;
2002  n2 = n1 = 0;
2003  tensor_order = std::max(tensor_order, size_type(3));
2004  } else if (r_type == GA_DSEMICOLON) {
2005  if (n1 != nbc1 || n2 != nbc2 || n3 != nbc3 ||
2006  tensor_order < 3)
2007  ga_throw_error(expr, pos-1, "Bad explicit "
2008  "vector/matrix/tensor format.");
2009  n3 = n2 = n1 = 0;
2010  tensor_order = std::max(tensor_order, size_type(4));
2011  } else if (r_type == GA_RBRACKET) {
2012  if (n1 != nbc1 || n2 != nbc2 || n3 != nbc3 ||
2013  tensor_order == 3)
2014  ga_throw_error(expr, pos-1, "Bad explicit "
2015  "vector/matrix/tensor format.");
2016  tree.current_node->nbc1 = nbc1;
2017  if (tensor_order == 4) {
2018  tree.current_node->nbc2 = nbc2/nbc1;
2019  tree.current_node->nbc3 = nbc3/nbc2;
2020  } else {
2021  tree.current_node->nbc2 = tree.current_node->nbc3 = 1;
2022  }
2023  } else {
2024  ga_throw_error(expr, pos-1, "The explicit "
2025  "vector/matrix/tensor components should be "
2026  "separated by ',', ';', ',,' and ';;' and "
2027  "be ended by ']'.");
2028  }
2029 
2030  } while (r_type != GA_RBRACKET);
2031  bgeot::multi_index mi;
2032  nbc1 = tree.current_node->nbc1; nbc2 = tree.current_node->nbc2;
2033  nbc3 = tree.current_node->nbc3;
2034 
2035  size_type nbl = tree.current_node->children.size()
2036  / (nbc2 * nbc1 * nbc3);
2037  switch(tensor_order) {
2038  case 1:
2039  /* mi.push_back(1); */ mi.push_back(nbc1); break;
2040  case 2:
2041  mi.push_back(nbl); if (nbc1 > 1) mi.push_back(nbc1); break;
2042  case 3:
2043  mi.push_back(nbl); mi.push_back(nbc2);
2044  mi.push_back(nbc1);
2045  break;
2046  case 4:
2047  mi.push_back(nbl); mi.push_back(nbc3);
2048  mi.push_back(nbc2); mi.push_back(nbc1);
2049  break;
2050  default: GMM_ASSERT1(false, "Internal error");
2051  }
2052  tree.current_node->tensor().adjust_sizes(mi);
2053  std::vector<pga_tree_node> children = tree.current_node->children;
2054  auto it = tree.current_node->children.begin();
2055  for (size_type i = 0; i < nbc1; ++i)
2056  for (size_type j = 0; j < nbc2; ++j)
2057  for (size_type k = 0; k < nbc3; ++k)
2058  for (size_type l = 0; l < nbl; ++l, ++it)
2059  *it = children[i+nbc1*(j+nbc2*(k+nbc3*l))];
2060  tree.current_node->marked = true;
2061  }
2062  }
2063  tree.current_node->nbc1 = tree.current_node->tensor().sizes().size();
2064  state = 2;
2065  break;
2066 
2067  default:
2068  ga_throw_error(expr, token_pos, "Unexpected token.");
2069  }
2070  break;
2071 
2072  case 2:
2073  switch (t_type) {
2074  case GA_PLUS: case GA_MINUS: case GA_MULT: case GA_DIV:
2075  case GA_COLON: case GA_DOT: case GA_DOTMULT: case GA_DOTDIV:
2076  case GA_TMULT:
2077  tree.add_op(t_type, token_pos, expr);
2078  state = 1; break;
2079  case GA_QUOTE:
2080  tree.add_op(t_type, token_pos, expr);
2081  state = 2; break;
2082  case GA_END: case GA_RPAR: case GA_COMMA: case GA_DCOMMA:
2083  case GA_RBRACKET: case GA_SEMICOLON: case GA_DSEMICOLON:
2084  return t_type;
2085  case GA_LPAR: // Parameter list
2086  {
2087  ga_tree sub_tree;
2088  GA_TOKEN_TYPE r_type;
2089  tree.add_params(token_pos, expr);
2090  do {
2091  r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
2092  if (r_type != GA_RPAR && r_type != GA_COMMA)
2093  ga_throw_error(expr, pos-((r_type != GA_END)?1:0),
2094  "Parameters should be separated "
2095  "by ',' and parameter list ended by ')'.");
2096  tree.add_sub_tree(sub_tree);
2097  } while (r_type != GA_RPAR);
2098  state = 2;
2099  }
2100  break;
2101 
2102  default:
2103  ga_throw_error(expr, token_pos, "Unexpected token.");
2104  }
2105  break;
2106  }
2107  }
2108 
2109  return GA_INVALID;
2110  }
2111 
2112  // Syntax analysis of a string. Conversion to a tree. register the macros.
2113  void ga_read_string_reg(const std::string &expr, ga_tree &tree,
2114  ga_macro_dictionary &macro_dict) {
2115  size_type pos = 0, token_pos, token_length;
2116  tree.clear();
2117  GA_TOKEN_TYPE t = ga_get_token(expr, pos, token_pos, token_length);
2118  if (t == GA_END) return;
2119  pos = 0;
2120  pstring nexpr(new std::string(expr));
2121 
2122  t = ga_read_term(nexpr, pos, tree, macro_dict);
2123  if (tree.root) ga_expand_macro(tree, tree.root, macro_dict);
2124 
2125  switch (t) {
2126  case GA_RPAR: ga_throw_error(nexpr, pos-1, "Unbalanced parenthesis.");
2127  case GA_RBRACKET: ga_throw_error(nexpr, pos-1, "Unbalanced braket.");
2128  case GA_END: break;
2129  default: ga_throw_error(nexpr, pos-1, "Unexpected token.");
2130  }
2131  }
2132 
2133  // Syntax analysis of a string. Conversion to a tree.
2134  // Do not register the macros (but expand them).
2135  void ga_read_string(const std::string &expr, ga_tree &tree,
2136  const ga_macro_dictionary &macro_dict) {
2137  ga_macro_dictionary macro_dict_loc(true, macro_dict);
2138  ga_read_string_reg(expr, tree, macro_dict_loc);
2139  }
2140 
2141  // Small tool to make basic substitutions into an assembly string
2142  // Should be replaced by macros now.
2143  std::string ga_substitute(const std::string &expr,
2144  const std::map<std::string, std::string> &dict) {
2145  if (dict.size()) {
2146  size_type pos = 0, token_pos, token_length;
2147  std::stringstream exprs;
2148 
2149  while (true) {
2150  GA_TOKEN_TYPE t_type = ga_get_token(expr, pos, token_pos, token_length);
2151  if (t_type == GA_END) return exprs.str();
2152  std::string name(&(expr[token_pos]), token_length);
2153  if (t_type == GA_NAME && dict.find(name) != dict.end())
2154  exprs << dict.at(name); else exprs << name;
2155  }
2156  }
2157  return expr;
2158  }
2159 
2160 } /* 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