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