GetFEM  5.4.3
getfem_generic_assembly_semantic.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 // Semantic analysis of assembly trees and semantic manipulations.
23 
24 
25 #include <getfem/getfem_generic_assembly_functions_and_operators.h>
26 #include <getfem/getfem_generic_assembly_semantic.h>
27 #include <getfem/getfem_generic_assembly_compile_and_exec.h>
28 
29 namespace getfem {
30 
31  extern bool predef_operators_nonlinear_elasticity_initialized;
32  extern bool predef_operators_plasticity_initialized;
33  extern bool predef_operators_contact_initialized;
34 
35  static void ga_node_derivation
36  (ga_tree &tree, const ga_workspace &workspace, const mesh &m,
37  pga_tree_node pnode, const std::string &varname,
38  const std::string &interpolatename, size_type order, bool any_trans = false);
39 
40  static void ga_node_grad(ga_tree &tree, const ga_workspace &workspace,
41  const mesh &m, pga_tree_node pnode);
42  static bool ga_node_mark_tree_for_grad(pga_tree_node pnode,
43  const ga_workspace &workspace);
44  static void ga_node_analysis(ga_tree &tree,
45  const ga_workspace &workspace,
46  pga_tree_node pnode, const mesh &me,
47  size_type ref_elt_dim, bool eval_fixed_size,
48  bool ignore_X, int option);
49 
50 
51  bool ga_extract_variables(const pga_tree_node pnode,
52  const ga_workspace &workspace,
53  const mesh &m,
54  std::set<var_trans_pair> &vars,
55  bool ignore_data) {
56  bool expand_groups = !ignore_data;
57  bool found_var = false;
58  if (pnode->node_type == GA_NODE_VAL ||
59  pnode->node_type == GA_NODE_GRAD ||
60  pnode->node_type == GA_NODE_HESS ||
61  pnode->node_type == GA_NODE_DIVERG ||
62  pnode->node_type == GA_NODE_INTERPOLATE_VAL ||
63  pnode->node_type == GA_NODE_INTERPOLATE_GRAD ||
64  pnode->node_type == GA_NODE_INTERPOLATE_HESS ||
65  pnode->node_type == GA_NODE_INTERPOLATE_DIVERG ||
66  pnode->node_type == GA_NODE_ELEMENTARY_VAL ||
67  pnode->node_type == GA_NODE_ELEMENTARY_GRAD ||
68  pnode->node_type == GA_NODE_ELEMENTARY_HESS ||
69  pnode->node_type == GA_NODE_ELEMENTARY_DIVERG ||
70  pnode->node_type == GA_NODE_SECONDARY_DOMAIN_VAL ||
71  pnode->node_type == GA_NODE_SECONDARY_DOMAIN_GRAD ||
72  pnode->node_type == GA_NODE_SECONDARY_DOMAIN_HESS ||
73  pnode->node_type == GA_NODE_SECONDARY_DOMAIN_DIVERG ||
74  pnode->node_type == GA_NODE_XFEM_PLUS_VAL ||
75  pnode->node_type == GA_NODE_XFEM_PLUS_GRAD ||
76  pnode->node_type == GA_NODE_XFEM_PLUS_HESS ||
77  pnode->node_type == GA_NODE_XFEM_PLUS_DIVERG ||
78  pnode->node_type == GA_NODE_XFEM_MINUS_VAL ||
79  pnode->node_type == GA_NODE_XFEM_MINUS_GRAD ||
80  pnode->node_type == GA_NODE_XFEM_MINUS_HESS ||
81  pnode->node_type == GA_NODE_XFEM_MINUS_DIVERG) {
82  bool group = workspace.variable_group_exists(pnode->name);
83  bool iscte = (!group) && workspace.is_constant(pnode->name);
84  if (!iscte) found_var = true;
85  if (!ignore_data || !iscte) {
86  if (group && expand_groups) {
87  for (const std::string &t : workspace.variable_group(pnode->name))
88  vars.insert(var_trans_pair(t, pnode->interpolate_name));
89  } else
90  vars.insert(var_trans_pair(pnode->name, pnode->interpolate_name));
91  }
92  }
93  if (pnode->node_type == GA_NODE_INTERPOLATE_VAL ||
94  pnode->node_type == GA_NODE_INTERPOLATE_GRAD ||
95  pnode->node_type == GA_NODE_INTERPOLATE_HESS ||
96  pnode->node_type == GA_NODE_INTERPOLATE_DIVERG ||
97  pnode->node_type == GA_NODE_INTERPOLATE_VAL_TEST ||
98  pnode->node_type == GA_NODE_INTERPOLATE_GRAD_TEST ||
99  pnode->node_type == GA_NODE_INTERPOLATE_HESS_TEST ||
100  pnode->node_type == GA_NODE_INTERPOLATE_DIVERG_TEST ||
101  pnode->node_type == GA_NODE_INTERPOLATE_X ||
102  pnode->node_type == GA_NODE_INTERPOLATE_ELT_K ||
103  pnode->node_type == GA_NODE_INTERPOLATE_ELT_B ||
104  pnode->node_type == GA_NODE_INTERPOLATE_NORMAL) {
105  workspace.interpolate_transformation(pnode->interpolate_name)
106  ->extract_variables(workspace, vars, ignore_data, m,
107  pnode->interpolate_name);
108  }
109  for (auto&& child : pnode->children)
110  found_var = ga_extract_variables(child, workspace, m,
111  vars, ignore_data)
112  || found_var;
113  return found_var;
114  }
115 
116 
117  static bool ga_node_mark_tree_for_variable
118  (pga_tree_node pnode, const ga_workspace &workspace, const mesh &m,
119  const std::string &varname,
120  const std::string &interpolatename, bool any_trans = false) {
121  bool marked = false;
122  for (pga_tree_node &child : pnode->children)
123  if (ga_node_mark_tree_for_variable(child, workspace, m,
124  varname, interpolatename, any_trans))
125  marked = true;
126 
127  bool plain_node(pnode->node_type == GA_NODE_VAL ||
128  pnode->node_type == GA_NODE_GRAD ||
129  pnode->node_type == GA_NODE_HESS ||
130  pnode->node_type == GA_NODE_DIVERG);
131  bool interpolate_node(pnode->node_type == GA_NODE_INTERPOLATE_VAL ||
132  pnode->node_type == GA_NODE_INTERPOLATE_GRAD ||
133  pnode->node_type == GA_NODE_INTERPOLATE_HESS ||
134  pnode->node_type == GA_NODE_INTERPOLATE_DIVERG);
135  bool elementary_node(pnode->node_type == GA_NODE_ELEMENTARY_VAL ||
136  pnode->node_type == GA_NODE_ELEMENTARY_GRAD ||
137  pnode->node_type == GA_NODE_ELEMENTARY_HESS ||
138  pnode->node_type == GA_NODE_ELEMENTARY_DIVERG);
139  bool secondary_node(pnode->node_type == GA_NODE_SECONDARY_DOMAIN_VAL ||
140  pnode->node_type == GA_NODE_SECONDARY_DOMAIN_GRAD ||
141  pnode->node_type == GA_NODE_SECONDARY_DOMAIN_HESS ||
142  pnode->node_type == GA_NODE_SECONDARY_DOMAIN_DIVERG);
143  bool xfem_node(pnode->node_type == GA_NODE_XFEM_PLUS_VAL ||
144  pnode->node_type == GA_NODE_XFEM_PLUS_GRAD ||
145  pnode->node_type == GA_NODE_XFEM_PLUS_HESS ||
146  pnode->node_type == GA_NODE_XFEM_PLUS_DIVERG ||
147  pnode->node_type == GA_NODE_XFEM_MINUS_VAL ||
148  pnode->node_type == GA_NODE_XFEM_MINUS_GRAD ||
149  pnode->node_type == GA_NODE_XFEM_MINUS_HESS ||
150  pnode->node_type == GA_NODE_XFEM_MINUS_DIVERG);
151  bool interpolate_test_node
152  (pnode->node_type == GA_NODE_INTERPOLATE_VAL_TEST ||
153  pnode->node_type == GA_NODE_INTERPOLATE_GRAD_TEST ||
154  pnode->node_type == GA_NODE_INTERPOLATE_HESS_TEST ||
155  pnode->node_type == GA_NODE_INTERPOLATE_DIVERG_TEST);
156 
157  if ((plain_node || interpolate_node || secondary_node ||
158  elementary_node || xfem_node) &&
159  (pnode->name.compare(varname) == 0 &&
160  (any_trans || pnode->interpolate_name.compare(interpolatename) == 0)))
161  marked = true;
162 
163  if (interpolate_node || interpolate_test_node ||
164  pnode->node_type == GA_NODE_INTERPOLATE_X ||
165  pnode->node_type == GA_NODE_INTERPOLATE_ELT_K ||
166  pnode->node_type == GA_NODE_INTERPOLATE_ELT_B ||
167  pnode->node_type == GA_NODE_INTERPOLATE_NORMAL) {
168  std::set<var_trans_pair> vars;
169  workspace.interpolate_transformation(pnode->interpolate_name)
170  ->extract_variables(workspace, vars, true,
171  m, pnode->interpolate_name);
172  for (const var_trans_pair &pair : vars)
173  if (pair.varname.compare(varname) == 0 &&
174  (any_trans || pair.transname.compare(interpolatename) == 0))
175  marked = true;
176  }
177  pnode->marked = marked;
178  return marked;
179  }
180 
181  static void ga_node_expand_expression_in_place_of_test
182  (ga_tree &tree, const ga_workspace &workspace,
183  pga_tree_node pnode, const std::string &varname,
184  pga_tree_node pexpr, ga_tree &grad_expr, ga_tree &hess_expr,
185  size_type order, const mesh &me, size_type ref_elt_dim, bool eval_fixed_size,
186  bool ignore_X, int option) {
187  pga_tree_node parent = pnode->parent;
188  for (pga_tree_node &child : pnode->children)
189  ga_node_expand_expression_in_place_of_test
190  (tree, workspace, child, varname, pexpr, grad_expr,
191  hess_expr, order, me, ref_elt_dim, eval_fixed_size, ignore_X, option);
192  const std::string &name = pnode->name;
193  size_type loc_order = pnode->test_function_type;
194 
195  if (loc_order == order && !(name.compare(varname))) {
196  bool need_grad = (pnode->node_type == GA_NODE_GRAD_TEST ||
197  pnode->node_type == GA_NODE_DIVERG_TEST ||
198  pnode->node_type == GA_NODE_HESS_TEST);
199  bool need_hess = (pnode->node_type == GA_NODE_HESS_TEST);
200 
201  if (need_grad && grad_expr.root == nullptr) {
202  tree.copy_node(pexpr, grad_expr.root);
203  if (ga_node_mark_tree_for_grad(grad_expr.root, workspace)) {
204  ga_node_grad(grad_expr, workspace, me, grad_expr.root);
205  ga_node_analysis(grad_expr, workspace, grad_expr.root, me,
206  ref_elt_dim, eval_fixed_size, ignore_X, option);
207  } else {
208  bgeot::multi_index mi = grad_expr.root->t.sizes();
209  mi.push_back(me.dim());
210  grad_expr.root->t.adjust_sizes(mi);
211  grad_expr.root->node_type = GA_NODE_ZERO;
212  gmm::clear(grad_expr.root->tensor().as_vector());
213  grad_expr.clear_children(grad_expr.root);
214  }
215  }
216 
217  if (need_hess && hess_expr.root == nullptr) {
218  tree.copy_node(grad_expr.root, hess_expr.root);
219  if (ga_node_mark_tree_for_grad(hess_expr.root, workspace)) {
220  ga_node_grad(hess_expr, workspace, me, hess_expr.root);
221  ga_node_analysis(hess_expr, workspace, hess_expr.root, me,
222  ref_elt_dim, eval_fixed_size, ignore_X, option);
223  } else {
224  bgeot::multi_index mi = hess_expr.root->t.sizes();
225  mi.push_back(me.dim());
226  hess_expr.root->t.adjust_sizes(mi);
227  hess_expr.root->node_type = GA_NODE_ZERO;
228  gmm::clear(hess_expr.root->tensor().as_vector());
229  hess_expr.clear_children(grad_expr.root);
230  }
231  }
232  switch(pnode->node_type) {
233  case GA_NODE_VAL_TEST:
234  case GA_NODE_GRAD_TEST:
235  case GA_NODE_HESS_TEST:
236  case GA_NODE_DIVERG_TEST:
237  {
238  pga_tree_node pnode_new = nullptr;
239  if (pnode->node_type == GA_NODE_VAL_TEST)
240  tree.copy_node(pexpr, pnode_new); // allocates new
241  else if (pnode->node_type == GA_NODE_GRAD_TEST ||
242  pnode->node_type == GA_NODE_DIVERG_TEST)
243  tree.copy_node(grad_expr.root, pnode_new); // allocates new
244  else if (pnode->node_type == GA_NODE_HESS_TEST)
245  tree.copy_node(hess_expr.root, pnode_new); // allocates new
246  pnode_new->parent = parent;
247  parent->replace_child(pnode, pnode_new);
248  if (pnode->node_type == GA_NODE_DIVERG_TEST) {
249  tree.insert_node(pnode_new, GA_NODE_OP);
250  pnode_new->parent->op_type = GA_TRACE;
251  }
252  delete pnode; // deallocates old
253  pnode = nullptr;
254  }
255  break;
256  case GA_NODE_INTERPOLATE_VAL_TEST: case GA_NODE_INTERPOLATE_GRAD_TEST:
257  case GA_NODE_INTERPOLATE_HESS_TEST: case GA_NODE_INTERPOLATE_DIVERG_TEST:
258  GMM_ASSERT1(pexpr->node_type == GA_NODE_VAL_TEST,
259  "Sorry, directional derivative does not work for the "
260  "moment with interpolate transformations. Future work.");
261  pnode->name = pexpr->name;
262  break;
263  case GA_NODE_ELEMENTARY_VAL_TEST: case GA_NODE_ELEMENTARY_GRAD_TEST:
264  case GA_NODE_ELEMENTARY_HESS_TEST: case GA_NODE_ELEMENTARY_DIVERG_TEST:
265  GMM_ASSERT1(pexpr->node_type == GA_NODE_VAL_TEST,
266  "Sorry, directional derivative does not work for the "
267  "moment with elementary transformations. Future work.");
268  pnode->name = pexpr->name;
269  break;
270  case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
271  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
272  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
273  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
274  GMM_ASSERT1(pexpr->node_type == GA_NODE_VAL_TEST,
275  "Sorry, directional derivative does not work for the "
276  "moment with secondary domains. Future work.");
277  pnode->name = pexpr->name;
278  break;
279  case GA_NODE_XFEM_PLUS_VAL_TEST: case GA_NODE_XFEM_PLUS_GRAD_TEST:
280  case GA_NODE_XFEM_PLUS_HESS_TEST: case GA_NODE_XFEM_PLUS_DIVERG_TEST:
281  case GA_NODE_XFEM_MINUS_VAL_TEST: case GA_NODE_XFEM_MINUS_GRAD_TEST:
282  case GA_NODE_XFEM_MINUS_HESS_TEST: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
283  GMM_ASSERT1(pexpr->node_type == GA_NODE_VAL_TEST,
284  "Sorry, directional derivative does not work for the "
285  "moment with Xfem_plus and Xfem_minus operations. "
286  "Future work.");
287  pnode->name = pexpr->name;
288  break;
289  default:
290  break;
291  }
292  }
293  }
294 
295  //=========================================================================
296  // Some hash code functions for node identification
297  //=========================================================================
298 
299  static scalar_type ga_hash_code(const std::string &s) {
300  scalar_type c(0);
301  for (size_type i = 0; i < s.size(); ++i)
302  c += sin(M_E+scalar_type(s[i]))+M_PI*M_E*scalar_type(i+1);
303  return c;
304  }
305 
306  static scalar_type ga_hash_code(const base_tensor &t) {
307  scalar_type c(0);
308  for (size_type i = 0; i < t.size(); ++i)
309  c += sin((1.+M_E)*t[i]+M_E*M_E*scalar_type(i+1))+scalar_type(i+1)*M_PI;
310  return c;
311  }
312 
313  static scalar_type ga_hash_code(GA_NODE_TYPE e) {
314  return cos(M_E + scalar_type((e == GA_NODE_ZERO) ? GA_NODE_CONSTANT : e));
315  }
316 
317  static scalar_type ga_hash_code(const pga_tree_node pnode) {
318  scalar_type c = ga_hash_code(pnode->node_type);
319 
320  switch (pnode->node_type) {
321  case GA_NODE_CONSTANT: case GA_NODE_ZERO:
322  c += ga_hash_code(pnode->tensor());
323  if (pnode->test_function_type & 1)
324  c += 34.731 * ga_hash_code(pnode->name_test1);
325  if (pnode->test_function_type & 2)
326  c += 34.731 * ga_hash_code(pnode->name_test2);
327  break;
328 
329  case GA_NODE_OP: c += scalar_type(pnode->op_type)*M_E*M_PI*M_PI; break;
330  case GA_NODE_X: c += scalar_type(pnode->nbc1) + M_E*M_PI; break;
331  case GA_NODE_VAL: case GA_NODE_GRAD:
332  case GA_NODE_HESS: case GA_NODE_DIVERG:
333  case GA_NODE_VAL_TEST: case GA_NODE_GRAD_TEST:
334  case GA_NODE_HESS_TEST: case GA_NODE_DIVERG_TEST:
335  c += ga_hash_code(pnode->name); break;
336 
337  case GA_NODE_INTERPOLATE_FILTER:
338  c += 1.73*ga_hash_code(pnode->interpolate_name)
339  + 2.486*double(pnode->nbc1 + 1);
340  break;
341  case GA_NODE_INTERPOLATE_DERIVATIVE:
342  c += 2.321*ga_hash_code(pnode->interpolate_name_der);
343  [[fallthrough]]; // The hash code is completed with the next item
344  case GA_NODE_INTERPOLATE_VAL: case GA_NODE_INTERPOLATE_GRAD:
345  case GA_NODE_INTERPOLATE_HESS: case GA_NODE_INTERPOLATE_DIVERG:
346  case GA_NODE_INTERPOLATE_VAL_TEST: case GA_NODE_INTERPOLATE_GRAD_TEST:
347  case GA_NODE_INTERPOLATE_HESS_TEST: case GA_NODE_INTERPOLATE_DIVERG_TEST:
348  case GA_NODE_SECONDARY_DOMAIN_VAL: case GA_NODE_SECONDARY_DOMAIN_GRAD:
349  case GA_NODE_SECONDARY_DOMAIN_HESS: case GA_NODE_SECONDARY_DOMAIN_DIVERG:
350  case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
351  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
352  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
353  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
354  c += 1.33*(1.22+ga_hash_code(pnode->name))
355  + 1.66*ga_hash_code(pnode->interpolate_name);
356  break;
357  case GA_NODE_ELEMENTARY_VAL: case GA_NODE_ELEMENTARY_GRAD:
358  case GA_NODE_ELEMENTARY_HESS: case GA_NODE_ELEMENTARY_DIVERG:
359  case GA_NODE_ELEMENTARY_VAL_TEST: case GA_NODE_ELEMENTARY_GRAD_TEST:
360  case GA_NODE_ELEMENTARY_HESS_TEST: case GA_NODE_ELEMENTARY_DIVERG_TEST:
361  c += 1.33*(1.22+ga_hash_code(pnode->name))
362  + 2.63*ga_hash_code(pnode->elementary_name)
363  + 3.47*ga_hash_code(pnode->elementary_target);
364  break;
365  case GA_NODE_XFEM_PLUS_VAL: case GA_NODE_XFEM_PLUS_GRAD:
366  case GA_NODE_XFEM_PLUS_HESS: case GA_NODE_XFEM_PLUS_DIVERG:
367  case GA_NODE_XFEM_PLUS_VAL_TEST: case GA_NODE_XFEM_PLUS_GRAD_TEST:
368  case GA_NODE_XFEM_PLUS_HESS_TEST: case GA_NODE_XFEM_PLUS_DIVERG_TEST:
369  case GA_NODE_XFEM_MINUS_VAL: case GA_NODE_XFEM_MINUS_GRAD:
370  case GA_NODE_XFEM_MINUS_HESS: case GA_NODE_XFEM_MINUS_DIVERG:
371  case GA_NODE_XFEM_MINUS_VAL_TEST: case GA_NODE_XFEM_MINUS_GRAD_TEST:
372  case GA_NODE_XFEM_MINUS_HESS_TEST: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
373  c += 1.33*(1.22+ga_hash_code(pnode->name));
374  break;
375  case GA_NODE_INTERPOLATE_X:
376  case GA_NODE_INTERPOLATE_ELT_K: case GA_NODE_INTERPOLATE_ELT_B:
377  case GA_NODE_INTERPOLATE_NORMAL:
378  case GA_NODE_SECONDARY_DOMAIN_X: case GA_NODE_SECONDARY_DOMAIN_NORMAL:
379  c += M_PI*1.33*ga_hash_code(pnode->interpolate_name);
380  break;
381  case GA_NODE_PREDEF_FUNC: case GA_NODE_SPEC_FUNC: case GA_NODE_OPERATOR:
382  c += ga_hash_code(pnode->name)
383  + tanh(scalar_type(pnode->der1)/M_PI + scalar_type(pnode->der2)*M_PI);
384  break;
385  default: break;
386  }
387  return c;
388  }
389 
390 # define ga_valid_operand(pnode) \
391  { \
392  if (pnode && (pnode->node_type == GA_NODE_PREDEF_FUNC || \
393  pnode->node_type == GA_NODE_SPEC_FUNC || \
394  pnode->node_type == GA_NODE_NAME || \
395  pnode->node_type == GA_NODE_OPERATOR || \
396  pnode->node_type == GA_NODE_ALLINDICES)) \
397  ga_throw_error(pnode->expr, pnode->pos, "Invalid term"); \
398  }
399 
400  static void ga_node_analysis(ga_tree &tree,
401  const ga_workspace &workspace,
402  pga_tree_node pnode, const mesh &me,
403  size_type ref_elt_dim, bool eval_fixed_size,
404  bool ignore_X, int option) {
405  // cout << "Analysis of "; ga_print_node(pnode, cout); cout << endl;
406  bool all_cte = true, all_sc = true;
407  size_type meshdim = (&me == &dummy_mesh()) ? 1 : me.dim();
408  pnode->symmetric_op = false;
409 
410  for (size_type i = 0; i < pnode->children.size(); ++i) {
411  ga_node_analysis(tree, workspace, pnode->children[i], me,
412  ref_elt_dim, eval_fixed_size, ignore_X, option);
413  all_cte = all_cte && (pnode->children[i]->is_constant());
414  all_sc = all_sc && (pnode->children[i]->tensor_proper_size() == 1);
415  if (pnode->children[i]->test_function_type == size_type(-1)) {
416  cerr << "node : "; ga_print_node(pnode, cerr); cerr << endl;
417  GMM_ASSERT1(false, "internal error on child " << i);
418  }
419  if (pnode->node_type != GA_NODE_PARAMS)
420  ga_valid_operand(pnode->children[i]);
421  }
422 
423  size_type nbch = pnode->children.size();
424  pga_tree_node child0 = (nbch > 0) ? pnode->children[0] : 0;
425  pga_tree_node child1 = (nbch > 1) ? pnode->children[1] : 0;
426  bgeot::multi_index mi;
427  const bgeot::multi_index &size0 = child0 ? child0->t.sizes() : mi;
428  const bgeot::multi_index &size1 = child1 ? child1->t.sizes() : mi;
429  size_type dim0 = child0 ? child0->tensor_order() : 0;
430  size_type dim1 = child1 ? child1->tensor_order() : 0;
431 
432  const ga_predef_function_tab &PREDEF_FUNCTIONS
433  = dal::singleton<ga_predef_function_tab>::instance(0);
434  const ga_predef_operator_tab &PREDEF_OPERATORS
435  = dal::singleton<ga_predef_operator_tab>::instance(0);
436  const ga_spec_function_tab &SPEC_FUNCTIONS
437  = dal::singleton<ga_spec_function_tab>::instance(0);
438 
439  switch (pnode->node_type) {
440  case GA_NODE_PREDEF_FUNC: case GA_NODE_OPERATOR: case GA_NODE_SPEC_FUNC:
441  case GA_NODE_CONSTANT: case GA_NODE_X: case GA_NODE_ELT_SIZE:
442  case GA_NODE_ELT_K: case GA_NODE_ELT_B: case GA_NODE_NORMAL:
443  case GA_NODE_RESHAPE: case GA_NODE_CROSS_PRODUCT:
444  case GA_NODE_IND_MOVE_LAST: case GA_NODE_SWAP_IND:
445  case GA_NODE_CONTRACT: case GA_NODE_INTERPOLATE_X:
446  case GA_NODE_INTERPOLATE_ELT_K: case GA_NODE_INTERPOLATE_ELT_B:
447  case GA_NODE_INTERPOLATE_NORMAL: case GA_NODE_SECONDARY_DOMAIN_X:
448  case GA_NODE_SECONDARY_DOMAIN_NORMAL:
449  pnode->test_function_type = 0;
450  break;
451 
452  case GA_NODE_ALLINDICES:
453  pnode->test_function_type = 0;
454  break;
455 
456  case GA_NODE_VAL:
457  if (eval_fixed_size && !(workspace.associated_mf(pnode->name))
458  && !(workspace.associated_im_data(pnode->name))) {
459  gmm::copy(workspace.value(pnode->name), pnode->tensor().as_vector());
460  pnode->node_type = GA_NODE_CONSTANT;
461  }
462  break;
463 
464  case GA_NODE_ZERO: case GA_NODE_GRAD:
465  case GA_NODE_HESS: case GA_NODE_DIVERG:
466  case GA_NODE_INTERPOLATE_VAL: case GA_NODE_INTERPOLATE_GRAD:
467  case GA_NODE_INTERPOLATE_HESS: case GA_NODE_INTERPOLATE_DIVERG:
468  case GA_NODE_ELEMENTARY_VAL: case GA_NODE_ELEMENTARY_GRAD:
469  case GA_NODE_ELEMENTARY_HESS: case GA_NODE_ELEMENTARY_DIVERG:
470  case GA_NODE_SECONDARY_DOMAIN_VAL: case GA_NODE_SECONDARY_DOMAIN_GRAD:
471  case GA_NODE_SECONDARY_DOMAIN_HESS: case GA_NODE_SECONDARY_DOMAIN_DIVERG:
472  case GA_NODE_XFEM_PLUS_VAL: case GA_NODE_XFEM_PLUS_GRAD:
473  case GA_NODE_XFEM_PLUS_HESS: case GA_NODE_XFEM_PLUS_DIVERG:
474  case GA_NODE_XFEM_MINUS_VAL: case GA_NODE_XFEM_MINUS_GRAD:
475  case GA_NODE_XFEM_MINUS_HESS: case GA_NODE_XFEM_MINUS_DIVERG:
476  break;
477 
478  case GA_NODE_VAL_TEST: case GA_NODE_GRAD_TEST:
479  case GA_NODE_HESS_TEST: case GA_NODE_DIVERG_TEST:
480  case GA_NODE_INTERPOLATE_VAL_TEST: case GA_NODE_INTERPOLATE_GRAD_TEST:
481  case GA_NODE_INTERPOLATE_HESS_TEST: case GA_NODE_INTERPOLATE_DIVERG_TEST:
482  case GA_NODE_INTERPOLATE_DERIVATIVE:
483  case GA_NODE_ELEMENTARY_VAL_TEST: case GA_NODE_ELEMENTARY_GRAD_TEST:
484  case GA_NODE_ELEMENTARY_HESS_TEST: case GA_NODE_ELEMENTARY_DIVERG_TEST:
485  case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
486  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
487  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
488  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
489  case GA_NODE_XFEM_PLUS_VAL_TEST: case GA_NODE_XFEM_PLUS_GRAD_TEST:
490  case GA_NODE_XFEM_PLUS_HESS_TEST: case GA_NODE_XFEM_PLUS_DIVERG_TEST:
491  case GA_NODE_XFEM_MINUS_VAL_TEST: case GA_NODE_XFEM_MINUS_GRAD_TEST:
492  case GA_NODE_XFEM_MINUS_HESS_TEST: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
493  {
494  const mesh_fem *mf = workspace.associated_mf(pnode->name);
495  const im_data *imd = workspace.associated_im_data(pnode->name);
496  size_type t_type = pnode->test_function_type;
497  if (t_type == 1) {
498  pnode->name_test1 = pnode->name;
499  pnode->interpolate_name_test1 = pnode->interpolate_name;
500  pnode->interpolate_name_test2 = pnode->name_test2 = "";
501  pnode->qdim1 = (mf || imd)
502  ? workspace.qdim(pnode->name)
503  : gmm::vect_size(workspace.value(pnode->name));
504  if (option == 1)
505  workspace.test1.insert
506  (var_trans_pair(pnode->name_test1,
507  pnode->interpolate_name_test1));
508  if (!(pnode->qdim1))
509  ga_throw_error(pnode->expr, pnode->pos,
510  "Invalid null size of variable");
511  } else {
512  pnode->interpolate_name_test1 = pnode->name_test1 = "";
513  pnode->name_test2 = pnode->name;
514  pnode->interpolate_name_test2 = pnode->interpolate_name;
515  pnode->qdim2 = (mf || imd)
516  ? workspace.qdim(pnode->name)
517  : gmm::vect_size(workspace.value(pnode->name));
518  if (option == 1)
519  workspace.test2.insert
520  (var_trans_pair(pnode->name_test2,
521  pnode->interpolate_name_test2));
522  if (!(pnode->qdim2))
523  ga_throw_error(pnode->expr, pnode->pos,
524  "Invalid null size of variable");
525  }
526  size_type q = workspace.qdim(pnode->name);
527  if (!q)
528  ga_throw_error(pnode->expr, pnode->pos,
529  "Invalid null size of variable");
530  if (!mf & !imd) { // global variable
531  if (q == 1) {
532  pnode->init_vector_tensor(1);
533  pnode->tensor()[0] = scalar_type(1);
534  } else
535  pnode->init_identity_matrix_tensor(q);
536  pnode->test_function_type = t_type;
537  } else if (imd) {
538  bgeot::multi_index mii = workspace.qdims(pnode->name);
539  if (q == 1 && mii.size() <= 1) {
540  pnode->init_vector_tensor(1);
541  pnode->tensor()[0] = scalar_type(1);
542  } else {
543  mii.insert(mii.begin(), q);
544  pnode->t.set_to_original();
545  pnode->t.adjust_sizes(mii);
546  auto itw = pnode->tensor().begin();
547  for (size_type i = 0; i < q; ++i) // set identity matrix
548  for (size_type j = 0; j < q; ++j)
549  *itw++ = (i == j) ? scalar_type(1) : scalar_type(0);
550  }
551  pnode->test_function_type = t_type;
552  }
553  }
554  break;
555 
556  case GA_NODE_SECONDARY_DOMAIN:
557  pnode->interpolate_name = tree.secondary_domain;
558  if (tree.secondary_domain.size() == 0)
559  ga_throw_error(pnode->expr, pnode->pos, "Secondary domain used "
560  "in a single domain term.");
561  [[fallthrough]];
562  case GA_NODE_INTERPOLATE:
563  if (pnode->name.compare("X") == 0) {
564  if (pnode->node_type == GA_NODE_INTERPOLATE) {
565  pnode->node_type = GA_NODE_INTERPOLATE_X;
566  pnode->init_vector_tensor(meshdim);
567  } else {
568  auto psd = workspace.secondary_domain(tree.secondary_domain);
569  pnode->node_type = GA_NODE_SECONDARY_DOMAIN_X;
570  pnode->init_vector_tensor(psd->mim().linked_mesh().dim());
571  }
572  break;
573  }
574  if (pnode->name.compare("Normal") == 0) {
575  if (pnode->node_type == GA_NODE_INTERPOLATE) {
576  pnode->node_type = GA_NODE_INTERPOLATE_NORMAL;
577  pnode->init_vector_tensor(meshdim);
578  } else {
579  auto psd = workspace.secondary_domain(tree.secondary_domain);
580  pnode->node_type = GA_NODE_SECONDARY_DOMAIN_NORMAL;
581  pnode->init_vector_tensor(psd->mim().linked_mesh().dim());
582  }
583  break;
584  }
585  if (pnode->name.compare("element_K") == 0) {
586  if (pnode->node_type == GA_NODE_INTERPOLATE) {
587  pnode->node_type = GA_NODE_INTERPOLATE_ELT_K;
588  if (ref_elt_dim == 1)
589  pnode->init_vector_tensor(meshdim);
590  else
591  pnode->init_matrix_tensor(meshdim, ref_elt_dim);
592  }
593  break;
594  }
595  if (pnode->name.compare("element_B") == 0) {
596  if (pnode->node_type == GA_NODE_INTERPOLATE) {
597  pnode->node_type = GA_NODE_INTERPOLATE_ELT_B;
598  pnode->init_matrix_tensor(ref_elt_dim, meshdim);
599  }
600  break;
601  }
602  [[fallthrough]];
603  case GA_NODE_ELEMENTARY: // +GA_NODE_SECONDARY_DOMAIN, GA_NODE_INTERPOLATE:
604  case GA_NODE_XFEM_PLUS:
605  case GA_NODE_XFEM_MINUS:
606  {
607  int ndt = ((pnode->node_type == GA_NODE_INTERPOLATE) ? 1 : 0)
608  + ((pnode->node_type == GA_NODE_ELEMENTARY) ? 2 : 0)
609  + ((pnode->node_type == GA_NODE_SECONDARY_DOMAIN) ? 3 : 0)
610  + ((pnode->node_type == GA_NODE_XFEM_PLUS) ? 4 : 0)
611  + ((pnode->node_type == GA_NODE_XFEM_MINUS) ? 5 : 0);
612  GMM_ASSERT1(ndt > 0 && ndt < 6, "internal error");
613  constexpr std::array<const char*,5>
614  op_name_selector{"Interpolation",
615  "Elementary_transformation",
616  "Secondary_domain",
617  "Xfem_plus",
618  "Xfem_minus"};
619  std::string op__name = op_name_selector.at(ndt-1);
620 
621  std::string name = pnode->name;
622  size_type prefix_id = ga_parse_prefix_operator(name);
623  size_type test = ga_parse_prefix_test(name);
624  pnode->name = name;
625 
626  if (ndt == 2) {
627  name = pnode->elementary_target;
628  ga_parse_prefix_operator(name);
629  ga_parse_prefix_test(name);
630  pnode->elementary_target = name;
631  }
632 
633  // Group must be tested and it should be a fem variable
634  if (!(workspace.variable_or_group_exists(name)))
635  ga_throw_error(pnode->expr, pnode->pos,
636  "Unknown variable or group of variables \""
637  + name + "\"");
638 
639  const mesh_fem *mf = workspace.associated_mf(name);
640  if (!mf)
641  ga_throw_error(pnode->expr, pnode->pos, op__name
642  << " can only apply to finite element variables/data");
643 
644  size_type q = workspace.qdim(name), n = mf->linked_mesh().dim();
645  if (!q) ga_throw_error(pnode->expr, pnode->pos,
646  "Invalid null size of variable");
647 
648  bgeot::multi_index mii = workspace.qdims(name);
649  if (mii.size() > 6)
650  ga_throw_error(pnode->expr, pnode->pos,
651  "Tensor with too many dimensions. Limited to 6");
652 
653  if (test == 1) {
654  pnode->name_test1 = pnode->name;
655  pnode->interpolate_name_test1 = pnode->interpolate_name;
656  if (option == 1)
657  workspace.test1.insert
658  (var_trans_pair(pnode->name_test1,
659  pnode->interpolate_name_test1));
660  pnode->qdim1 = workspace.qdim(name);
661  if (!(pnode->qdim1))
662  ga_throw_error(pnode->expr, pnode->pos,
663  "Invalid null size of variable");
664  } else if (test == 2) {
665  pnode->name_test2 = pnode->name;
666  pnode->interpolate_name_test2 = pnode->interpolate_name;
667  if (option == 1)
668  workspace.test2.insert
669  (var_trans_pair(pnode->name_test2,
670  pnode->interpolate_name_test2));
671  pnode->qdim2 = workspace.qdim(name);
672  if (!(pnode->qdim2))
673  ga_throw_error(pnode->expr, pnode->pos,
674  "Invalid null size of variable");
675  }
676 
677  constexpr std::array<GA_NODE_TYPE,5>
678  node_type_selector_val{GA_NODE_INTERPOLATE_VAL,
679  GA_NODE_ELEMENTARY_VAL,
680  GA_NODE_SECONDARY_DOMAIN_VAL,
681  GA_NODE_XFEM_PLUS_VAL,
682  GA_NODE_XFEM_MINUS_VAL},
683  node_type_selector_val_test{GA_NODE_INTERPOLATE_VAL_TEST,
684  GA_NODE_ELEMENTARY_VAL_TEST,
685  GA_NODE_SECONDARY_DOMAIN_VAL_TEST,
686  GA_NODE_XFEM_PLUS_VAL_TEST,
687  GA_NODE_XFEM_MINUS_VAL_TEST},
688  node_type_selector_grad{GA_NODE_INTERPOLATE_GRAD,
689  GA_NODE_ELEMENTARY_GRAD,
690  GA_NODE_SECONDARY_DOMAIN_GRAD,
691  GA_NODE_XFEM_PLUS_GRAD,
692  GA_NODE_XFEM_MINUS_GRAD},
693  node_type_selector_grad_test{GA_NODE_INTERPOLATE_GRAD_TEST,
694  GA_NODE_ELEMENTARY_GRAD_TEST,
695  GA_NODE_SECONDARY_DOMAIN_GRAD_TEST,
696  GA_NODE_XFEM_PLUS_GRAD_TEST,
697  GA_NODE_XFEM_MINUS_GRAD_TEST},
698  node_type_selector_hess{GA_NODE_INTERPOLATE_HESS,
699  GA_NODE_ELEMENTARY_HESS,
700  GA_NODE_SECONDARY_DOMAIN_HESS,
701  GA_NODE_XFEM_PLUS_HESS,
702  GA_NODE_XFEM_MINUS_HESS},
703  node_type_selector_hess_test{GA_NODE_INTERPOLATE_HESS_TEST,
704  GA_NODE_ELEMENTARY_HESS_TEST,
705  GA_NODE_SECONDARY_DOMAIN_HESS_TEST,
706  GA_NODE_XFEM_PLUS_HESS_TEST,
707  GA_NODE_XFEM_MINUS_HESS_TEST},
708  node_type_selector_div{GA_NODE_INTERPOLATE_DIVERG,
709  GA_NODE_ELEMENTARY_DIVERG,
710  GA_NODE_SECONDARY_DOMAIN_DIVERG,
711  GA_NODE_XFEM_PLUS_DIVERG,
712  GA_NODE_XFEM_MINUS_DIVERG},
713  node_type_selector_div_test{GA_NODE_INTERPOLATE_DIVERG_TEST,
714  GA_NODE_ELEMENTARY_DIVERG_TEST,
715  GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST,
716  GA_NODE_XFEM_PLUS_DIVERG_TEST,
717  GA_NODE_XFEM_MINUS_DIVERG_TEST};
718 
719  switch (prefix_id) {
720  case 0: // value
721  pnode->node_type = test ? node_type_selector_val_test[ndt-1]
722  : node_type_selector_val[ndt-1];
723  if (test) {
724  if (q == 1 && mii.size() <= 1) {
725  mii.resize(1);
726  mii[0] = 2;
727  } else
728  mii.insert(mii.begin(), 2);
729  }
730  break;
731  case 1: // grad
732  pnode->node_type = test ? node_type_selector_grad_test[ndt-1]
733  : node_type_selector_grad[ndt-1];
734  if (test) {
735  if (q == 1 && mii.size() <= 1) {
736  mii.resize(1);
737  mii[0] = 2;
738  } else
739  mii.insert(mii.begin(), 2);
740  if (n > 1) mii.push_back(n);
741  } else if (n > 1) {
742  if (q == 1 && mii.size() == 1)
743  mii[0] = n;
744  else
745  mii.push_back(n);
746  }
747  break;
748  case 2: // Hessian
749  pnode->node_type = test ? node_type_selector_hess_test[ndt-1]
750  : node_type_selector_hess[ndt-1];
751  if (test) {
752  if (q == 1 && mii.size() <= 1) {
753  mii.resize(1);
754  mii[0] = 2;
755  } else
756  mii.insert(mii.begin(), 2);
757  if (n > 1) {
758  mii.push_back(n);
759  mii.push_back(n);
760  }
761  } else if (n > 1) {
762  if (q == 1 && mii.size() == 1)
763  mii[0] = n;
764  else
765  mii.push_back(n);
766  mii.push_back(n);
767  }
768  break;
769  case 3: // divergence
770  if (q != n)
771  ga_throw_error(pnode->expr, pnode->pos,
772  "Divergence operator requires fem qdim ("
773  << q << ") to be equal to dim (" << n << ")");
774  pnode->node_type = test ? node_type_selector_div_test[ndt-1]
775  : node_type_selector_div[ndt-1];
776  mii.resize(1);
777  mii[0] = test ? 2 : 1;
778  break;
779  }
780  pnode->t.adjust_sizes(mii);
781  pnode->test_function_type = test;
782 
783  if (ndt == 1) {
784  if (!(workspace.interpolate_transformation_exists
785  (pnode->interpolate_name)))
786  ga_throw_error(pnode->expr, pnode->pos,
787  "Unknown interpolate transformation");
788  } else if (ndt == 2) {
789  if (!(workspace.elementary_transformation_exists
790  (pnode->elementary_name)))
791  ga_throw_error(pnode->expr, pnode->pos,
792  "Unknown elementary transformation");
793  if (!(workspace.variable_or_group_exists(pnode->elementary_target)))
794  ga_throw_error(pnode->expr, pnode->pos, "Unknown data or variable "
795  << pnode->elementary_target);
796  const mesh_fem *mft = workspace.associated_mf(name);
797  if (!mft)
798  ga_throw_error(pnode->expr, pnode->pos,
799  "Thir argument of the elementary transformation "
800  "should be a finite element variables/data");
801  } else if (ndt == 3) {
802  if (!(workspace.secondary_domain_exists
803  (pnode->interpolate_name)))
804  ga_throw_error(pnode->expr, pnode->pos,
805  "Unknown secondary domain");
806  }
807  }
808  break;
809 
810  case GA_NODE_INTERPOLATE_FILTER:
811  {
812  if (pnode->children.size() == 2) {
813  bool valid = (child1->node_type == GA_NODE_CONSTANT);
814  int n = valid ? int(round(child1->tensor()[0])) : -1;
815  if (n < 0 || n > 100 || child1->tensor_order() > 0)
816  ga_throw_error(pnode->expr, pnode->pos, "The third argument of "
817  "Interpolate_filter should be a (small) "
818  "non-negative integer.");
819  pnode->nbc1 = size_type(n);
820  tree.clear_node(child1);
821  }
822  if (!(workspace.interpolate_transformation_exists
823  (pnode->interpolate_name)))
824  ga_throw_error(pnode->expr, pnode->pos,
825  "Unknown interpolate transformation");
826  pnode->t = child0->t;
827  pnode->test_function_type = child0->test_function_type;
828  pnode->name_test1 = child0->name_test1;
829  pnode->name_test2 = child0->name_test2;
830  pnode->interpolate_name_test1 = child0->interpolate_name_test1;
831  pnode->interpolate_name_test2 = child0->interpolate_name_test2;
832  pnode->qdim1 = child0->qdim1;
833  pnode->qdim2 = child0->qdim2;
834  }
835  break;
836 
837 
838  case GA_NODE_OP:
839  switch(pnode->op_type) {
840 
841  case GA_PLUS: case GA_MINUS:
842  {
843  if (pnode->op_type == GA_PLUS) pnode->symmetric_op = true;
844  size_type c_size = std::min(size0.size(), size1.size());
845  bool compatible = true;
846 
847  size_type f_ind = 0;
848  if (child0->test_function_type &&
849  child1->test_function_type == child0->test_function_type)
850  f_ind = (child0->test_function_type == 3) ? 2:1;
851 
852  for (size_type i = f_ind; i < c_size; ++i)
853  if (size0[i] != size1[i]) compatible = false;
854  for (size_type i = c_size; i < size0.size(); ++i)
855  if (size0[i] != 1) compatible = false;
856  for (size_type i = c_size; i < size1.size(); ++i)
857  if (size1[i] != 1) compatible = false;
858 
859  if (!compatible)
860  ga_throw_error(pnode->expr, pnode->pos,
861  "Addition or substraction of expressions of "
862  "different sizes: " << size0 << " != " << size1);
863  if (child0->test_function_type || child1->test_function_type) {
864  switch (option) {
865  case 0: case 2:
866  if (child0->name_test1.compare(child1->name_test1) ||
867  child0->name_test2.compare(child1->name_test2) ||
868  child0->interpolate_name_test1.compare
869  (child1->interpolate_name_test1) ||
870  child0->interpolate_name_test2.compare
871  (child1->interpolate_name_test2))
872  compatible = false;
873  break;
874  case 1: case 3: break;
875  default: GMM_ASSERT1(false, "Unknown option");
876  }
877  }
878 
879  if (child0->test_function_type != child1->test_function_type ||
880  (!compatible && option != 2))
881  ga_throw_error(pnode->expr, pnode->pos, "Addition or substraction "
882  "of incompatible test functions");
883  if (all_cte) {
884  pnode->node_type = GA_NODE_CONSTANT;
885  pnode->test_function_type = 0;
886  pnode->tensor() = pnode->children[0]->tensor();
887  if (pnode->op_type == GA_MINUS)
888  pnode->tensor() -= pnode->children[1]->tensor();
889  else
890  pnode->tensor() += pnode->children[1]->tensor();
891  tree.clear_children(pnode);
892  } else {
893  pnode->t = child0->t;
894  pnode->test_function_type = child0->test_function_type;
895  pnode->name_test1 = child0->name_test1;
896  pnode->name_test2 = child0->name_test2;
897  pnode->interpolate_name_test1 = child0->interpolate_name_test1;
898  pnode->interpolate_name_test2 = child0->interpolate_name_test2;
899  pnode->qdim1 = child0->qdim1;
900  pnode->qdim2 = child0->qdim2;
901 
902  // simplification if one of the two operands is constant and zero
903  if (child0->tensor_is_zero()) {
904  if (pnode->op_type == GA_MINUS) {
905  pnode->op_type = GA_UNARY_MINUS;
906  tree.clear_node(child0);
907  } else {
908  tree.replace_node_by_child(pnode, 1);
909  pnode = child1;
910  }
911  } else if (child1->tensor_is_zero()) {
912  tree.replace_node_by_child(pnode, 0);
913  pnode = child0;
914  } else if (option == 2 && !compatible) {
915  bool child0_compatible = true, child1_compatible = true;
916  if (pnode->test_function_type & 1) {
917  if (child0->name_test1.compare(workspace.selected_test1.varname)
918  || child0->interpolate_name_test1.compare
919  (workspace.selected_test1.transname))
920  child0_compatible = false;
921  if (child1->name_test1.compare(workspace.selected_test1.varname)
922  || child1->interpolate_name_test1.compare
923  (workspace.selected_test1.transname))
924  child1_compatible = false;
925  }
926  if (pnode->test_function_type & 2) {
927  if (child0->name_test2.compare(workspace.selected_test2.varname)
928  || child0->interpolate_name_test2.compare
929  (workspace.selected_test2.transname))
930  child0_compatible = false;
931  if (child1->name_test2.compare(workspace.selected_test2.varname)
932  || child1->interpolate_name_test1.compare
933  (workspace.selected_test2.transname))
934  child1_compatible = false;
935  }
936  if (child0_compatible) {
937  tree.replace_node_by_child(pnode, 0);
938  pnode = child0;
939  } else if (child1_compatible) {
940  if (pnode->op_type == GA_MINUS) {
941  pnode->op_type = GA_UNARY_MINUS;
942  pnode->t = child1->t;
943  pnode->test_function_type = child1->test_function_type;
944  pnode->name_test1 = child1->name_test1;
945  pnode->name_test2 = child1->name_test2;
946  pnode->interpolate_name_test1=child1->interpolate_name_test1;
947  pnode->interpolate_name_test2=child1->interpolate_name_test2;
948  pnode->qdim1 = child1->qdim1;
949  pnode->qdim2 = child1->qdim2;
950  tree.clear_node(child0);
951  } else {
952  tree.replace_node_by_child(pnode, 1);
953  pnode = child1;
954  }
955  }
956  }
957  }
958  }
959  break;
960 
961  case GA_DOTMULT: case GA_DOTDIV:
962  {
963  if (pnode->op_type == GA_DOTMULT) pnode->symmetric_op = true;
964  bool compatible = true;
965  if (child0->tensor_proper_size() != child1->tensor_proper_size())
966  compatible = false;
967 
968  if (child0->tensor_proper_size() != 1) {
969  if (child0->tensor_order() != child1->tensor_order())
970  compatible = false;
971 
972  for (size_type i = 0; i < child0->tensor_order(); ++i)
973  if (child0->tensor_proper_size(i)!=child1->tensor_proper_size(i))
974  compatible = false;
975  }
976 
977  if (!compatible)
978  ga_throw_error(pnode->expr, pnode->pos,
979  "Arguments of different sizes for .* or ./");
980 
981  if (pnode->op_type == GA_DOTDIV && child1->test_function_type)
982  ga_throw_error(pnode->expr, pnode->pos,
983  "Division by test functions is not allowed");
984 
985  pnode->mult_test(child0, child1);
986  mi = pnode->t.sizes();
987  for (size_type i = 0; i < child0->tensor_order(); ++i)
988  mi.push_back(child0->tensor_proper_size(i));
989  pnode->t.adjust_sizes(mi);
990 
991  if (all_cte) {
992  pnode->node_type = GA_NODE_CONSTANT;
993  pnode->test_function_type = 0;
994  if (pnode->op_type == GA_DOTMULT) {
995  for (size_type i = 0; i < child0->tensor().size(); ++i)
996  pnode->tensor()[i] = child0->tensor()[i] * child1->tensor()[i];
997  } else {
998  for (size_type i = 0; i < child0->tensor().size(); ++i) {
999  if (child1->tensor()[i] == scalar_type(0))
1000  ga_throw_error(pnode->expr, pnode->pos, "Division by zero.");
1001  pnode->tensor()[i] = child0->tensor()[i] / child1->tensor()[i];
1002  }
1003  }
1004  tree.clear_children(pnode);
1005  } else {
1006  if (child0->tensor_is_zero() || child1->tensor_is_zero()) {
1007  gmm::clear(pnode->tensor().as_vector());
1008  pnode->node_type = GA_NODE_ZERO;
1009  tree.clear_children(pnode);
1010  }
1011  if (child1->tensor_is_zero() && pnode->op_type == GA_DOTDIV)
1012  ga_throw_error(pnode->expr, pnode->pos, "Division by zero.");
1013  }
1014  }
1015  break;
1016 
1017  case GA_UNARY_MINUS:
1018  pnode->t = child0->t;
1019  pnode->test_function_type = child0->test_function_type;
1020  pnode->name_test1 = child0->name_test1;
1021  pnode->name_test2 = child0->name_test2;
1022  pnode->interpolate_name_test1 = child0->interpolate_name_test1;
1023  pnode->interpolate_name_test2 = child0->interpolate_name_test2;
1024  pnode->qdim1 = child0->qdim1;
1025  pnode->qdim2 = child0->qdim2;
1026  if (all_cte) {
1027  pnode->node_type = GA_NODE_CONSTANT;
1028  pnode->test_function_type = 0;
1029  gmm::scale(pnode->tensor().as_vector(), scalar_type(-1));
1030  tree.clear_children(pnode);
1031  } else if (child0->node_type == GA_NODE_ZERO) {
1032  tree.replace_node_by_child(pnode, 0);
1033  pnode = child0;
1034  }
1035  break;
1036 
1037  case GA_QUOTE:
1038  mi = size0;
1039  if (child0->tensor_proper_size() == 1)
1040  { tree.replace_node_by_child(pnode, 0); pnode = child0; break; }
1041  else if (dim0 == 1)
1042  { size_type N = mi.back(); mi.back() = 1; mi.push_back(N); }
1043  else std::swap(mi[child0->nb_test_functions()],
1044  mi[child0->nb_test_functions()+1]);
1045 
1046 
1047  pnode->t.adjust_sizes(mi);
1048  pnode->test_function_type = child0->test_function_type;
1049  pnode->name_test1 = child0->name_test1;
1050  pnode->name_test2 = child0->name_test2;
1051  pnode->interpolate_name_test1 = child0->interpolate_name_test1;
1052  pnode->interpolate_name_test2 = child0->interpolate_name_test2;
1053  pnode->qdim1 = child0->qdim1;
1054  pnode->qdim2 = child0->qdim2;
1055  if (child0->node_type == GA_NODE_ZERO) {
1056  pnode->node_type = GA_NODE_ZERO;
1057  gmm::clear(pnode->tensor().as_vector());
1058  tree.clear_children(pnode);
1059  } else if (all_cte) {
1060  pnode->node_type = GA_NODE_CONSTANT;
1061  pnode->test_function_type = 0;
1062 
1063  if (dim0 == 1) {
1064  for (size_type i = 0; i < mi.back(); ++i)
1065  pnode->tensor()(0, i) = child0->tensor()[i];
1066  } else {
1067  size_type n1 = child0->tensor_proper_size(0);
1068  size_type n2 = child0->tensor_proper_size(1);
1069  size_type nn = child0->tensor().size()/(n1*n2);
1070  auto it = pnode->tensor().begin();
1071  for (size_type i = 0; i < nn; ++i)
1072  for (size_type j = 0; j < n1; ++j)
1073  for (size_type k = 0; k < n2; ++k, ++it)
1074  *it = child0->tensor()[j+k*n1+i*n1*n2];
1075  GA_DEBUG_ASSERT(it == pnode->tensor().end(), "Wrong sizes");
1076  }
1077  tree.clear_children(pnode);
1078  }
1079  break;
1080 
1081  case GA_SYM: case GA_SKEW:
1082  if (child0->tensor_proper_size() != 1 &&
1083  (dim0 != 2 || size0.back() != size0[size0.size()-2]))
1084  ga_throw_error(pnode->expr, pnode->pos, "Sym and Skew operators are "
1085  "for square matrices only.");
1086  mi = size0;
1087  if (child0->tensor_proper_size() == 1) {
1088  if (pnode->op_type == GA_SYM)
1089  { tree.replace_node_by_child(pnode, 0); pnode = child0; break; }
1090  else {
1091  pnode->node_type = GA_NODE_ZERO;
1092  gmm::clear(pnode->tensor().as_vector());
1093  tree.clear_children(pnode);
1094  break;
1095  }
1096  }
1097 
1098  pnode->t.adjust_sizes(mi);
1099  pnode->test_function_type = child0->test_function_type;
1100  pnode->name_test1 = child0->name_test1;
1101  pnode->name_test2 = child0->name_test2;
1102  pnode->interpolate_name_test1 = child0->interpolate_name_test1;
1103  pnode->interpolate_name_test2 = child0->interpolate_name_test2;
1104  pnode->qdim1 = child0->qdim1;
1105  pnode->qdim2 = child0->qdim2;
1106  if (all_cte) {
1107  pnode->node_type = GA_NODE_CONSTANT;
1108  pnode->test_function_type = 0;
1109  for (size_type i = 0; i < mi.back(); ++i)
1110  for (size_type j = 0; j < mi.back(); ++j)
1111  if (pnode->op_type == GA_SYM)
1112  pnode->tensor()(j, i) = 0.5*(child0->tensor()(j,i)
1113  + child0->tensor()(i,j));
1114  else
1115  pnode->tensor()(j, i) = 0.5*(child0->tensor()(j,i)
1116  - child0->tensor()(i,j));
1117  tree.clear_children(pnode);
1118  } else if (child0->node_type == GA_NODE_ZERO) {
1119  pnode->node_type = GA_NODE_ZERO;
1120  gmm::clear(pnode->tensor().as_vector());
1121  tree.clear_children(pnode);
1122  }
1123  break;
1124 
1125  case GA_TRACE:
1126  {
1127  mi = size0;
1128  size_type N = (child0->tensor_proper_size() == 1) ? 1 : mi.back();
1129 
1130  if (child0->tensor_proper_size() == 1)
1131  { tree.replace_node_by_child(pnode, 0); pnode = child0; break; }
1132 
1133  if ((dim0 != 2 && child0->tensor_proper_size() != 1) ||
1134  (dim0 == 2 && mi[mi.size()-2] != N))
1135  ga_throw_error(pnode->expr, pnode->pos,
1136  "Trace operator is for square matrices only.");
1137 
1138  if (dim0 == 2) { mi.pop_back(); mi.pop_back(); }
1139  pnode->t.adjust_sizes(mi);
1140  pnode->test_function_type = child0->test_function_type;
1141  pnode->name_test1 = child0->name_test1;
1142  pnode->name_test2 = child0->name_test2;
1143  pnode->interpolate_name_test1 = child0->interpolate_name_test1;
1144  pnode->interpolate_name_test2 = child0->interpolate_name_test2;
1145  pnode->qdim1 = child0->qdim1;
1146  pnode->qdim2 = child0->qdim2;
1147  if (all_cte) {
1148  pnode->node_type = GA_NODE_CONSTANT;
1149  pnode->test_function_type = 0;
1150  if (dim0 == 2) {
1151  pnode->tensor()[0] = scalar_type(0);
1152  for (size_type i = 0; i < N; ++i)
1153  pnode->tensor()[0] += child0->tensor()(i,i);
1154  } else {
1155  pnode->tensor()[0] += child0->tensor()[0];
1156  }
1157  tree.clear_children(pnode);
1158  } else if (child0->node_type == GA_NODE_ZERO) {
1159  pnode->node_type = GA_NODE_ZERO;
1160  gmm::clear(pnode->tensor().as_vector());
1161  tree.clear_children(pnode);
1162  }
1163  }
1164  break;
1165 
1166  case GA_DEVIATOR:
1167  {
1168  mi = size0;
1169  size_type N = (child0->tensor_proper_size() == 1) ? 1 : mi.back();
1170 
1171  if ((dim0 != 2 && child0->tensor_proper_size() != 1) ||
1172  (dim0 == 2 && mi[mi.size()-2] != N))
1173  ga_throw_error(pnode->expr, pnode->pos,
1174  "Deviator operator is for square matrices only.");
1175 
1176  if (child0->tensor_proper_size() == 1) {
1177  pnode->node_type = GA_NODE_ZERO;
1178  gmm::clear(pnode->tensor().as_vector());
1179  tree.clear_children(pnode);
1180  break;
1181  }
1182 
1183  pnode->t.adjust_sizes(mi);
1184  pnode->test_function_type = child0->test_function_type;
1185  pnode->name_test1 = child0->name_test1;
1186  pnode->name_test2 = child0->name_test2;
1187  pnode->interpolate_name_test1 = child0->interpolate_name_test1;
1188  pnode->interpolate_name_test2 = child0->interpolate_name_test2;
1189  pnode->qdim1 = child0->qdim1;
1190  pnode->qdim2 = child0->qdim2;
1191  if (all_cte) {
1192  pnode->node_type = GA_NODE_CONSTANT;
1193  pnode->test_function_type = 0;
1194  if (dim0 == 2) {
1195  scalar_type tr(0);
1196  gmm::copy(child0->tensor().as_vector(),
1197  pnode->tensor().as_vector());
1198  for (size_type i = 0; i < N; ++i)
1199  tr += child0->tensor()(i,i);
1200  for (size_type i = 0; i < N; ++i)
1201  pnode->tensor()(i,i) -= tr / scalar_type(N);
1202  } else {
1203  pnode->tensor()[0] = scalar_type(0);
1204  }
1205  tree.clear_children(pnode);
1206  } else if (child0->node_type == GA_NODE_ZERO) {
1207  pnode->node_type = GA_NODE_ZERO;
1208  gmm::clear(pnode->tensor().as_vector());
1209  tree.clear_children(pnode);
1210  }
1211  }
1212  break;
1213 
1214  case GA_PRINT:
1215  {
1216  pnode->t = child0->t;
1217  pnode->test_function_type = child0->test_function_type;
1218  pnode->name_test1 = child0->name_test1;
1219  pnode->name_test2 = child0->name_test2;
1220  pnode->interpolate_name_test1 = child0->interpolate_name_test1;
1221  pnode->interpolate_name_test2 = child0->interpolate_name_test2;
1222  pnode->qdim1 = child0->qdim1;
1223  pnode->qdim2 = child0->qdim2;
1224  if (all_cte) {
1225  pnode->node_type = GA_NODE_CONSTANT;
1226  cout << "Print constant term "; ga_print_node(child0, cout);
1227  cout << ": " << pnode->tensor() << endl;
1228  tree.clear_children(pnode);
1229  } else if (child0->node_type == GA_NODE_ZERO) {
1230  pnode->node_type = GA_NODE_ZERO;
1231  gmm::clear(pnode->tensor().as_vector());
1232  cout << "Print zero term "; ga_print_node(child0, cout);
1233  cout << ": " << pnode->tensor() << endl;
1234  tree.clear_children(pnode);
1235  }
1236  }
1237  break;
1238 
1239  case GA_DOT:
1240  {
1241  size_type s0 = dim0 == 0 ? 1 : size0.back();
1242  size_type s1 = dim1 == 0 ? 1 : size1[size1.size()-dim1];
1243 
1244  if (s0 != s1) ga_throw_error(pnode->expr, pnode->pos, "Dot product "
1245  "of expressions of different sizes ("
1246  << s0 << " != " << s1 << ").");
1247  if (dim0 <= 1 && dim1 <= 1) pnode->symmetric_op = true;
1248  pnode->mult_test(child0, child1);
1249  if (dim0 > 1 || dim1 > 1) {
1250  mi = pnode->t.sizes();
1251  for (size_type i = 1; i < dim0; ++i)
1252  mi.push_back(child0->tensor_proper_size(i-1));
1253  for (size_type i = 1; i < dim1; ++i)
1254  mi.push_back(child1->tensor_proper_size(i));
1255  pnode->t.adjust_sizes(mi);
1256  }
1257 
1258  if (child0->tensor_is_zero() || child1->tensor_is_zero()) {
1259  gmm::clear(pnode->tensor().as_vector());
1260  pnode->node_type = GA_NODE_ZERO;
1261  tree.clear_children(pnode);
1262  } else if (all_cte) {
1263  gmm::clear(pnode->tensor().as_vector());
1264  size_type m0 = child0->tensor().size() / s0;
1265  size_type m1 = child1->tensor().size() / s1;
1266  for (size_type i = 0; i < m0; ++i)
1267  for (size_type j = 0; j < m1; ++j)
1268  for (size_type k = 0; k < s0; ++k)
1269  pnode->tensor()[i*m1+j]
1270  += child0->tensor()[i*s0+k] * child1->tensor()[k*m1+j];
1271  pnode->node_type = GA_NODE_CONSTANT;
1272  tree.clear_children(pnode);
1273  }
1274  }
1275  break;
1276 
1277  case GA_COLON:
1278  {
1279  size_type s00 = (dim0 == 0) ? 1
1280  : (dim0 == 1 ? size0.back() : size0[size0.size()-2]);
1281  size_type s01 = (dim0 >= 2) ? size0.back() : 1;
1282  size_type s10 = (dim1 == 0) ? 1 : child1->tensor_proper_size(0);
1283  size_type s11 = (dim1 < 2) ? 1 : child1->tensor_proper_size(1);
1284  if (s00 != s10 || s01 != s11)
1285  ga_throw_error(pnode->expr, pnode->pos, "Frobenius product "
1286  "of expressions of different sizes ("
1287  << s00 << "," << s01 << " != " << s10 << ","
1288  << s11 << ").");
1289  if (child0->tensor_order() <= 2 && child1->tensor_order() <= 2)
1290  pnode->symmetric_op = true;
1291  pnode->mult_test(child0, child1);
1292  if (dim0 > 2 || dim1 > 2) {
1293  mi = pnode->t.sizes();
1294  for (size_type i = 0; i < dim0-2; ++i)
1295  mi.push_back(child0->tensor_proper_size(i));
1296  for (size_type i = 2; i < dim1; ++i)
1297  mi.push_back(child1->tensor_proper_size(i));
1298  pnode->t.adjust_sizes(mi);
1299  }
1300 
1301  if (child0->tensor_is_zero() || child1->tensor_is_zero()) {
1302  gmm::clear(pnode->tensor().as_vector());
1303  pnode->node_type = GA_NODE_ZERO;
1304  tree.clear_children(pnode);
1305  } else if (all_cte) {
1306  pnode->node_type = GA_NODE_CONSTANT;
1307  gmm::clear(pnode->tensor().as_vector());
1308  size_type k = 0;
1309  for (size_type i = 0, j = 0; i < child0->tensor().size(); ++i) {
1310  pnode->tensor()[j] += child0->tensor()[i] * child1->tensor()[k];
1311  ++j; if (j == pnode->tensor().size()) { j = 0; ++k; }
1312  }
1313  GMM_ASSERT1(k == child1->tensor().size(), "Internal error");
1314  tree.clear_children(pnode);
1315  }
1316  }
1317  break;
1318 
1319  case GA_TMULT:
1320  if (all_cte) {
1321  pnode->node_type = GA_NODE_CONSTANT;
1322  pnode->test_function_type = 0;
1323  if (child0->tensor().size() == 1 && child1->tensor().size() == 1) {
1324  pnode->init_scalar_tensor
1325  (child0->tensor()[0] * child1->tensor()[0]);
1326  } else if (child0->tensor().size() == 1) {
1327  pnode->t = child1->t;
1328  gmm::scale(pnode->tensor().as_vector(),
1329  scalar_type(child0->tensor()[0]));
1330  } else if (child1->tensor().size() == 1) {
1331  pnode->t = child0->t;
1332  gmm::scale(pnode->tensor().as_vector(),
1333  scalar_type(child1->tensor()[0]));
1334  } else {
1335  if (dim0+dim1 > 6)
1336  ga_throw_error(pnode->expr, pnode->pos, "Unauthorized "
1337  "tensor multiplication.");
1338  for (size_type i = 0; i < dim0; ++i)
1339  mi.push_back(child0->tensor().size(i));
1340  for (size_type i = 0; i < dim1; ++i)
1341  mi.push_back(child1->tensor().size(i));
1342  pnode->t.adjust_sizes(mi);
1343  size_type n0 = child0->tensor().size();
1344  size_type n1 = child1->tensor().size();
1345  for (size_type i = 0; i < n0; ++i)
1346  for (size_type j = 0; j < n1; ++j)
1347  pnode->tensor()[i+j*n0]=child0->tensor()[i]*child1->tensor()[j];
1348  }
1349  tree.clear_children(pnode);
1350  } else {
1351  pnode->mult_test(child0, child1);
1352  mi = pnode->t.sizes();
1353  if (child0->tensor_proper_size() != 1
1354  || child1->tensor_proper_size() != 1) {
1355  if (child0->tensor_proper_size() == 1) {
1356  for (size_type i = 0; i < dim1; ++i)
1357  mi.push_back(child1->tensor_proper_size(i));
1358  } else if (child1->tensor().size() == 1) {
1359  for (size_type i = 0; i < dim0; ++i)
1360  mi.push_back(child0->tensor_proper_size(i));
1361  } else {
1362  if (dim0+dim1 > 6)
1363  ga_throw_error(pnode->expr, pnode->pos, "Unauthorized "
1364  "tensor multiplication.");
1365  for (size_type i = 0; i < dim0; ++i)
1366  mi.push_back(child0->tensor_proper_size(i));
1367  for (size_type i = 0; i < dim1; ++i)
1368  mi.push_back(child1->tensor_proper_size(i));
1369  }
1370  pnode->t.adjust_sizes(mi);
1371  }
1372  if (child0->tensor_is_zero() || child1->tensor_is_zero()) {
1373  gmm::clear(pnode->tensor().as_vector());
1374  pnode->node_type = GA_NODE_ZERO;
1375  tree.clear_children(pnode);
1376  }
1377  }
1378  break;
1379 
1380  case GA_MULT:
1381  if (all_cte) {
1382  pnode->node_type = GA_NODE_CONSTANT;
1383  pnode->test_function_type = 0;
1384  if (child0->tensor_proper_size() == 1 &&
1385  child1->tensor_proper_size() == 1) {
1386  pnode->init_scalar_tensor(child0->tensor()[0]*child1->tensor()[0]);
1387  } else if (child0->tensor_proper_size() == 1) {
1388  pnode->t = child1->t;
1389  gmm::scale(pnode->tensor().as_vector(), child0->tensor()[0]);
1390  } else if (child1->tensor_proper_size() == 1) {
1391  pnode->t = child0->t;
1392  gmm::scale(pnode->tensor().as_vector(), child1->tensor()[0]);
1393  } else if (dim0 == 2 && dim1 == 1) {
1394  size_type m=child0->tensor().size(0), n=child0->tensor().size(1);
1395  if (n != child1->tensor().size(0))
1396  ga_throw_error(pnode->expr, pnode->pos,
1397  "Incompatible sizes in matrix-vector "
1398  "multiplication (" << n << " != "
1399  << child1->tensor().size(0) << ").");
1400  pnode->init_vector_tensor(m);
1401  gmm::clear(pnode->tensor().as_vector());
1402  for (size_type i = 0; i < m; ++i)
1403  for (size_type j = 0; j < n; ++j)
1404  pnode->tensor()[i] += child0->tensor()(i,j)*child1->tensor()[j];
1405  } else if (dim0 == 2 && dim1 == 2) {
1406  size_type m = child0->tensor().size(0);
1407  size_type n = child0->tensor().size(1);
1408  size_type p = child1->tensor().size(1);
1409  if (n != child1->tensor().size(0))
1410  ga_throw_error(pnode->expr, pnode->pos,
1411  "Incompatible sizes in matrix-matrix "
1412  "multiplication (" << n << " != "
1413  << child1->tensor().size(0) << ").");
1414  pnode->init_matrix_tensor(m,p);
1415  gmm::clear(pnode->tensor().as_vector());
1416  for (size_type i = 0; i < m; ++i)
1417  for (size_type j = 0; j < n; ++j)
1418  for (size_type k = 0; k < p; ++k)
1419  pnode->tensor()(i,k) += child0->tensor()(i,j)
1420  * child1->tensor()(j,k);
1421  }
1422  else if (dim0 == 4 && dim1 == 2) {
1423  size_type m=child0->tensor().size(0), n=child0->tensor().size(1);
1424  size_type o=child0->tensor().size(2), p=child0->tensor().size(3);
1425  if (o != child1->tensor().size(0) || p != child1->tensor().size(1))
1426  ga_throw_error(pnode->expr, pnode->pos,
1427  "Incompatible sizes in tensor-matrix "
1428  "multiplication (" << o << "," << p << " != "
1429  << child1->tensor().size(0) << ","
1430  << child1->tensor().size(1) << ").");
1431  pnode->init_matrix_tensor(m,n);
1432  gmm::clear(pnode->tensor().as_vector());
1433  for (size_type i = 0; i < m; ++i)
1434  for (size_type j = 0; j < n; ++j)
1435  for (size_type k = 0; k < o; ++k)
1436  for (size_type l = 0; l < p; ++l)
1437  pnode->tensor()(i,j) += child0->tensor()(i,j,k,l)
1438  * child1->tensor()(k,l);
1439  } else ga_throw_error(pnode->expr, pnode->pos,
1440  "Unauthorized multiplication.");
1441  tree.clear_children(pnode);
1442  } else {
1443  pnode->mult_test(child0, child1);
1444  mi = pnode->t.sizes();
1445 
1446  if (child0->tensor_proper_size() == 1 &&
1447  child1->tensor_proper_size() == 1) {
1448  pnode->symmetric_op = true;
1449  } else if (child0->tensor_proper_size() == 1) {
1450  pnode->symmetric_op = true;
1451  for (size_type i = 0; i < dim1; ++i)
1452  mi.push_back(child1->tensor_proper_size(i));
1453  } else if (child1->tensor_proper_size() == 1) {
1454  pnode->symmetric_op = true;
1455  for (size_type i = 0; i < dim0; ++i)
1456  mi.push_back(child0->tensor_proper_size(i));
1457  } else if (child0->tensor_order() == 2 &&
1458  child1->tensor_order() == 1) {
1459  size_type m = child0->tensor_proper_size(0);
1460  size_type n = child0->tensor_proper_size(1);
1461  mi.push_back(m);
1462  if (n != child1->tensor_proper_size(0))
1463  ga_throw_error(pnode->expr, pnode->pos,
1464  "Incompatible sizes in matrix-vector "
1465  "multiplication (" << n << " != "
1466  << child1->tensor_proper_size(0) << ").");
1467  } else if (child0->tensor_order() == 2 &&
1468  child1->tensor_order() == 2) {
1469  size_type m = child0->tensor_proper_size(0);
1470  size_type n = child0->tensor_proper_size(1);
1471  size_type p = child1->tensor_proper_size(1);
1472  mi.push_back(m); mi.push_back(p);
1473  if (n != child1->tensor_proper_size(0))
1474  ga_throw_error(pnode->expr, pnode->pos,
1475  "Incompatible sizes in matrix-matrix "
1476  "multiplication (" << n << " != "
1477  << child1->tensor_proper_size(0) << ").");
1478  }
1479  else if (pnode->children[0]->tensor_order() == 4 &&
1480  pnode->children[1]->tensor_order() == 2) {
1481  size_type m = child0->tensor_proper_size(0);
1482  size_type n = child0->tensor_proper_size(1);
1483  size_type o = child0->tensor_proper_size(2);
1484  size_type p = child0->tensor_proper_size(3);
1485  mi.push_back(m); mi.push_back(n);
1486  if (o != child1->tensor_proper_size(0) ||
1487  p != child1->tensor_proper_size(1))
1488  ga_throw_error(pnode->expr, pnode->pos,
1489  "Incompatible sizes in tensor-matrix "
1490  "multiplication (" << o << "," << p << " != "
1491  << child1->tensor_proper_size(0) << ","
1492  << child1->tensor_proper_size(1) << ").");
1493  } else ga_throw_error(pnode->expr, pnode->pos,
1494  "Unauthorized multiplication.");
1495  pnode->t.adjust_sizes(mi);
1496  // Simplifications
1497  if (child0->tensor_is_zero() || child1->tensor_is_zero()) {
1498  gmm::clear(pnode->tensor().as_vector());
1499  pnode->node_type = GA_NODE_ZERO;
1500  tree.clear_children(pnode);
1501  } else if (child0->node_type == GA_NODE_CONSTANT &&
1502  child0->tensor().size() == 1 &&
1503  child0->tensor()[0] == scalar_type(1)) {
1504  tree.replace_node_by_child(pnode, 1);
1505  pnode = child1;
1506  } else if (child1->node_type == GA_NODE_CONSTANT &&
1507  child1->tensor().size() == 1 &&
1508  child1->tensor()[0] == scalar_type(1)) {
1509  tree.replace_node_by_child(pnode, 0);
1510  pnode = child0;
1511  }
1512  }
1513  break;
1514 
1515  case GA_DIV:
1516  if (child1->tensor_proper_size() > 1)
1517  ga_throw_error(pnode->expr, pnode->pos,
1518  "Only the division by a scalar is allowed. "
1519  "Got a size of " << child1->tensor_proper_size());
1520  if (child1->test_function_type)
1521  ga_throw_error(pnode->expr, pnode->pos,
1522  "Division by test functions is not allowed.");
1523  if (child1->node_type == GA_NODE_CONSTANT &&
1524  child1->tensor()[0] == scalar_type(0))
1525  ga_throw_error(pnode->expr, pnode->children[1]->pos,
1526  "Division by zero");
1527 
1528  pnode->t = child0->t;
1529  pnode->test_function_type = child0->test_function_type;
1530  pnode->name_test1 = child0->name_test1;
1531  pnode->name_test2 = child0->name_test2;
1532  pnode->interpolate_name_test1 = child0->interpolate_name_test1;
1533  pnode->interpolate_name_test2 = child0->interpolate_name_test2;
1534  pnode->qdim1 = child0->qdim1;
1535  pnode->qdim2 = child0->qdim2;
1536 
1537  if (all_cte) {
1538  pnode->node_type = GA_NODE_CONSTANT;
1539  pnode->t = pnode->children[0]->t;
1540  pnode->test_function_type = 0;
1541  gmm::scale(pnode->tensor().as_vector(),
1542  scalar_type(1) / pnode->children[1]->tensor()[0]);
1543  tree.clear_children(pnode);
1544  } else if (child0->tensor_is_zero()) {
1545  gmm::clear(pnode->tensor().as_vector());
1546  pnode->node_type = GA_NODE_ZERO;
1547  tree.clear_children(pnode);
1548  } else if (child1->node_type == GA_NODE_CONSTANT &&
1549  child1->tensor().size() == 1 &&
1550  child1->tensor()[0] == scalar_type(1)) {
1551  tree.replace_node_by_child(pnode, 0);
1552  pnode = child0;
1553  }
1554  break;
1555 
1556  default:GMM_ASSERT1(false, "Unexpected operation. Internal error.");
1557  }
1558  break;
1559 
1560  case GA_NODE_C_MATRIX:
1561  {
1562  if (!all_sc) ga_throw_error(pnode->expr, pnode->pos,
1563  "Constant vector/matrix/tensor "
1564  "components should be scalar valued.");
1565  GMM_ASSERT1(pnode->children.size() == pnode->tensor_proper_size(),
1566  "Internal error");
1567 
1568  pnode->test_function_type = 0;
1569  for (size_type i = 0; i < pnode->children.size(); ++i) {
1570  if (pnode->children[i]->test_function_type) {
1571  if (pnode->test_function_type == 0) {
1572  pnode->test_function_type=pnode->children[i]->test_function_type;
1573  pnode->name_test1 = pnode->children[i]->name_test1;
1574  pnode->name_test2 = pnode->children[i]->name_test2;
1575  pnode->interpolate_name_test1
1576  = pnode->children[i]->interpolate_name_test1;
1577  pnode->interpolate_name_test2
1578  = pnode->children[i]->interpolate_name_test2;
1579  pnode->qdim1 = pnode->children[i]->qdim1;
1580  pnode->qdim2 = pnode->children[i]->qdim2;
1581  } else {
1582  if (pnode->test_function_type !=
1583  pnode->children[i]->test_function_type ||
1584  pnode->name_test1.compare(pnode->children[i]->name_test1) ||
1585  pnode->name_test2.compare(pnode->children[i]->name_test2) ||
1586  pnode->interpolate_name_test1.compare
1587  (pnode->children[i]->interpolate_name_test1) ||
1588  pnode->interpolate_name_test2.compare
1589  (pnode->children[i]->interpolate_name_test2))
1590  ga_throw_error(pnode->expr, pnode->pos, "Inconsistent use of "
1591  "test function in constant matrix.");
1592  }
1593  }
1594  }
1595  int to_add = int(pnode->nb_test_functions() + pnode->nbc1)
1596  - int(pnode->tensor().sizes().size());
1597  GMM_ASSERT1(to_add >= 0 && to_add <= 2, "Internal error");
1598  if (to_add) {
1599  mi = pnode->tensor().sizes();
1600  mi.resize(pnode->nbc1+pnode->nb_test_functions());
1601  for (int i = int(mi.size()-1); i >= to_add; --i)
1602  mi[i] = mi[i-to_add];
1603  for (int i = 0; i < to_add; ++i)
1604  mi[i] = 2;
1605 
1606  if (pnode->test_function_type & 1) {
1607  const mesh_fem *mf1 = workspace.associated_mf(pnode->name_test1);
1608  const im_data *imd1 = workspace.associated_im_data(pnode->name_test1);
1609  if (mf1 == 0 && imd1 == 0) // global variable
1610  mi[0] = gmm::vect_size(workspace.value(pnode->name_test1));
1611  else if (imd1) // im_data variable
1612  mi[0] = workspace.qdim(pnode->name_test1); // == 1 because of all_sc = true
1613  }
1614  if (pnode->test_function_type & 2) {
1615  const mesh_fem *mf2 = workspace.associated_mf(pnode->name_test2);
1616  const im_data *imd2 = workspace.associated_im_data(pnode->name_test2);
1617  if (mf2 == 0 && imd2 == 0) // global variable
1618  mi[(pnode->test_function_type & 1) ? 1 : 0]
1619  = gmm::vect_size(workspace.value(pnode->name_test2)); // == 1 because of all_sc = true
1620  else if (imd2) // im_data variable
1621  mi[(pnode->test_function_type & 1) ? 1 : 0]
1622  = workspace.qdim(pnode->name_test2);
1623  }
1624  pnode->tensor().adjust_sizes(mi);
1625  }
1626 
1627  if (all_cte) {
1628  bool all_zero = true;
1629  for (size_type i = 0; i < pnode->children.size(); ++i) {
1630  pnode->tensor()[i] = pnode->children[i]->tensor()[0];
1631  if (pnode->tensor()[i] != scalar_type(0)) all_zero = false;
1632  }
1633  if (all_zero)
1634  pnode->node_type = GA_NODE_ZERO;
1635  else
1636  pnode->node_type = GA_NODE_CONSTANT;
1637  tree.clear_children(pnode);
1638  }
1639  }
1640  break;
1641 
1642 
1643  case GA_NODE_NAME:
1644  {
1645  std::string name = pnode->name;
1646 
1647  if (!ignore_X && !(name.compare("X"))) {
1648  pnode->node_type = GA_NODE_X;
1649  pnode->nbc1 = 0;
1650  pnode->init_vector_tensor(meshdim);
1651  break;
1652  }
1653  if (!(name.compare("Diff"))) {
1654  pnode->test_function_type = 0;
1655  break;
1656  }
1657  if (!(name.compare("Grad"))) {
1658  pnode->test_function_type = 0;
1659  break;
1660  }
1661  if (!(name.compare("element_size"))) {
1662  pnode->node_type = GA_NODE_ELT_SIZE;
1663  pnode->init_scalar_tensor(0);
1664  break;
1665  }
1666  if (!(name.compare("Cross_product"))) {
1667  pnode->node_type = GA_NODE_CROSS_PRODUCT;
1668  pnode->test_function_type = 0;
1669  break;
1670  }
1671  if (!(name.compare("element_K"))) {
1672  pnode->node_type = GA_NODE_ELT_K;
1673  if (ref_elt_dim == 1)
1674  pnode->init_vector_tensor(meshdim);
1675  else
1676  pnode->init_matrix_tensor(meshdim, ref_elt_dim);
1677  break;
1678  }
1679  if (!(name.compare("element_B"))) {
1680  pnode->node_type = GA_NODE_ELT_B;
1681  pnode->init_matrix_tensor(ref_elt_dim, meshdim);
1682  break;
1683  }
1684  if (!(name.compare("Normal"))) {
1685  pnode->node_type = GA_NODE_NORMAL;
1686  pnode->init_vector_tensor(meshdim);
1687  break;
1688  }
1689  if (!(name.compare("Reshape"))) {
1690  pnode->node_type = GA_NODE_RESHAPE;
1691  pnode->init_scalar_tensor(0);
1692  break;
1693  }
1694  if (!(name.compare("Swap_indices"))) {
1695  pnode->node_type = GA_NODE_SWAP_IND;
1696  pnode->init_scalar_tensor(0);
1697  break;
1698  }
1699  if (!(name.compare("Index_move_last"))) {
1700  pnode->node_type = GA_NODE_IND_MOVE_LAST;
1701  pnode->init_scalar_tensor(0);
1702  break;
1703  }
1704  if (!(name.compare("Contract"))) {
1705  pnode->node_type = GA_NODE_CONTRACT;
1706  pnode->init_scalar_tensor(0);
1707  break;
1708  }
1709 
1710  if (name.compare(0, 11, "Derivative_") == 0) {
1711  name = name.substr(11);
1712  bool valid = true;
1713  pnode->der1 = 1; pnode->der2 = 0;
1714  char *p;
1715  size_type d = strtol(name.c_str(), &p, 10);
1716  size_type s = p - name.c_str();
1717  if (s > 0) {
1718  pnode->der1 = d;
1719  if (name[s] != '_') valid = false; else
1720  name = name.substr(s+1);
1721  }
1722  d = strtol(name.c_str(), &p, 10);
1723  s = p - name.c_str();
1724  if (s > 0) {
1725  pnode->der2 = d;
1726  if (name[s] != '_') valid = false; else
1727  name = name.substr(s+1);
1728  }
1729  if (!valid || pnode->der1 == 0)
1730  ga_throw_error(pnode->expr, pnode->pos,"Invalid derivative format");
1731  }
1732 
1733  ga_predef_function_tab::const_iterator it=PREDEF_FUNCTIONS.find(name);
1734  if (it != PREDEF_FUNCTIONS.end()) {
1735  // Predefined function found
1736  pnode->node_type = GA_NODE_PREDEF_FUNC;
1737  pnode->name = name;
1738  pnode->test_function_type = 0;
1739  if (pnode->der1) {
1740  if (pnode->der1 > it->second.nbargs()
1741  || pnode->der2 > it->second.nbargs())
1742  ga_throw_error(pnode->expr, pnode->pos, "Invalid derivative.");
1743  const ga_predef_function &F = it->second;
1744  if ((F.ftype() == 0 || F.dtype() == 2) && !(pnode->der2)) {
1745  pnode->name = ((pnode->der1 == 1) ?
1746  F.derivative1() : F.derivative2());
1747  pnode->der1 = pnode->der2 = 0;
1748  }
1749  }
1750  } else if (SPEC_FUNCTIONS.find(name) != SPEC_FUNCTIONS.end()) {
1751  // Special function found
1752  if (pnode->der1)
1753  ga_throw_error(pnode->expr, pnode->pos, "Special functions do not "
1754  "support derivatives.");
1755  pnode->node_type = GA_NODE_SPEC_FUNC;
1756  pnode->name = name;
1757  pnode->test_function_type = 0;
1758  if (!name.compare("pi")) {
1759  pnode->node_type = GA_NODE_CONSTANT;
1760  pnode->init_scalar_tensor(M_PI);
1761  } else if (!name.compare("meshdim")) {
1762  pnode->node_type = GA_NODE_CONSTANT;
1763  pnode->init_scalar_tensor(scalar_type(meshdim));
1764  } else if (!name.compare("timestep")) {
1765  pnode->node_type = GA_NODE_CONSTANT;
1766  pnode->init_scalar_tensor(scalar_type(workspace.get_time_step()));
1767  }
1768  } else if (PREDEF_OPERATORS.tab.find(name)
1769  != PREDEF_OPERATORS.tab.end()) {
1770  // Nonlinear operator found
1771  pnode->node_type = GA_NODE_OPERATOR;
1772  pnode->name = name;
1773  pnode->test_function_type = 0;
1774  } else {
1775  // Search for a variable name with optional gradient, Hessian,
1776  // divergence or test functions
1777 
1778  size_type prefix_id = ga_parse_prefix_operator(name);
1779  size_type test = ga_parse_prefix_test(name);
1780 
1781  if (!(workspace.variable_exists(name)))
1782  ga_throw_error(pnode->expr, pnode->pos, "Unknown variable, "
1783  "function, operator or data \"" + name + "\"");
1784 
1785  if (pnode->der1)
1786  ga_throw_error(pnode->expr, pnode->pos, "Derivative is for "
1787  "functions or operators, not for variables. "
1788  "Use Grad instead.");
1789  pnode->name = name;
1790 
1791  const mesh_fem *mf = workspace.associated_mf(name);
1792  const im_data *imd = workspace.associated_im_data(name);
1793 
1794  if (test && workspace.is_constant(name) &&
1795  !(workspace.is_disabled_variable(name)))
1796  ga_throw_error(pnode->expr, pnode->pos, "Test functions of "
1797  "constants are not allowed.");
1798  if (test == 1) {
1799  pnode->name_test1 = name;
1800  pnode->interpolate_name_test1 = "";
1801  if (option == 1)
1802  workspace.test1.insert
1803  (var_trans_pair(pnode->name_test1,
1804  pnode->interpolate_name_test1));
1805  pnode->qdim1 = (mf || imd) ? workspace.qdim(name)
1806  : gmm::vect_size(workspace.value(name));
1807  if (!(pnode->qdim1))
1808  ga_throw_error(pnode->expr, pnode->pos,
1809  "Invalid null size of variable");
1810  } else if (test == 2) {
1811  pnode->name_test2 = name;
1812  pnode->interpolate_name_test2 = "";
1813  if (option == 1)
1814  workspace.test2.insert
1815  (var_trans_pair(pnode->name_test2,
1816  pnode->interpolate_name_test2));
1817  pnode->qdim2 = (mf || imd) ? workspace.qdim(name)
1818  : gmm::vect_size(workspace.value(name));
1819  if (!(pnode->qdim2))
1820  ga_throw_error(pnode->expr, pnode->pos,
1821  "Invalid null size of variable");
1822  }
1823 
1824  if (!mf && !imd) { // global variable
1825  if (prefix_id)
1826  ga_throw_error(pnode->expr, pnode->pos, "Gradient, Hessian or "
1827  "Divergence cannot be evaluated for fixed size data.");
1828  if (test)
1829  pnode->node_type = GA_NODE_VAL_TEST;
1830  else if (eval_fixed_size)
1831  pnode->node_type = GA_NODE_CONSTANT;
1832  else
1833  pnode->node_type = GA_NODE_VAL;
1834 
1835  size_type q = gmm::vect_size(workspace.value(name));
1836  if (q == 1) {
1837  if (test) {
1838  pnode->init_vector_tensor(1);
1839  pnode->tensor()[0] = scalar_type(1);
1840  }
1841  else
1842  pnode->init_scalar_tensor(workspace.value(name)[0]);
1843  } else {
1844  if (test) {
1845  pnode->init_identity_matrix_tensor(q);
1846  } else {
1847  pnode->t.adjust_sizes(workspace.qdims(name));
1848  gmm::copy(workspace.value(name), pnode->tensor().as_vector());
1849  }
1850  }
1851  } else if (imd) { // im_data variable
1852  size_type q = workspace.qdim(name);
1853  bgeot::multi_index mii = workspace.qdims(name);
1854 
1855  if (!q) ga_throw_error(pnode->expr, pnode->pos,
1856  "Invalid null size of variable " << name);
1857  if (mii.size() > 6)
1858  ga_throw_error(pnode->expr, pnode->pos,
1859  "Tensor with too many dimensions. Limited to 6");
1860  if (prefix_id)
1861  ga_throw_error(pnode->expr, pnode->pos, "Gradient, Hessian or "
1862  "Divergence cannot be evaluated for im data.");
1863 
1864  pnode->node_type = test ? GA_NODE_VAL_TEST : GA_NODE_VAL;
1865 
1866  if (test) {
1867  if (q == 1 && mii.size() <= 1) {
1868  pnode->init_vector_tensor(1);
1869  pnode->tensor()[0] = scalar_type(1);
1870  } else {
1871  mii.insert(mii.begin(), q);
1872  pnode->t.set_to_original();
1873  pnode->t.adjust_sizes(mii);
1874  auto itw = pnode->tensor().begin();
1875  for (size_type i = 0; i < q; ++i) // set identity matrix
1876  for (size_type j = 0; j < q; ++j)
1877  *itw++ = (i == j) ? scalar_type(1) : scalar_type(0);
1878  }
1879  } else
1880  pnode->t.adjust_sizes(mii);
1881  } else { // mesh_fem variable
1882  size_type q = workspace.qdim(name);
1883  size_type n = mf->linked_mesh().dim();
1884  bgeot::multi_index mii = workspace.qdims(name);
1885 
1886  if (!q) ga_throw_error(pnode->expr, pnode->pos,
1887  "Invalid null size of variable " << name);
1888  if (mii.size() > 6)
1889  ga_throw_error(pnode->expr, pnode->pos,
1890  "Tensor with too many dimensions. Limited to 6");
1891 
1892  switch (prefix_id) {
1893  case 0: // value
1894  pnode->node_type = test ? GA_NODE_VAL_TEST : GA_NODE_VAL;
1895  // For Test nodes a first dimension of size equal to 2 has to be
1896  // prepended by convention (to be adapted later)
1897  if (test && q == 1 && mii.size() <= 1) {
1898  mii.resize(1);
1899  mii[0] = 2;
1900  } else if (test)
1901  mii.insert(mii.begin(), 2);
1902  break;
1903  case 1: // grad
1904  pnode->node_type = test ? GA_NODE_GRAD_TEST : GA_NODE_GRAD;
1905  if (test) {
1906  if (q == 1 && mii.size() <= 1) {
1907  mii.resize(1);
1908  mii[0] = 2;
1909  } else
1910  mii.insert(mii.begin(), 2);
1911  }
1912  if (n > 1) {
1913  if (mii.size() == 1 && mii[0] == 1) mii[0] = n;
1914  else mii.push_back(n);
1915  }
1916  break;
1917  case 2: // Hessian
1918  pnode->node_type = test ? GA_NODE_HESS_TEST : GA_NODE_HESS;
1919  if (test) {
1920  if (q == 1 && mii.size() <= 1) {
1921  mii.resize(1);
1922  mii[0] = 2;
1923  } else
1924  mii.insert(mii.begin(), 2);
1925  }
1926  if (n > 1) {
1927  if (mii.size() == 1 && mii[0] == 1) mii[0] = n;
1928  else mii.push_back(n);
1929  mii.push_back(n);
1930  }
1931  break;
1932  case 3: // divergence
1933  pnode->node_type = test ? GA_NODE_DIVERG_TEST : GA_NODE_DIVERG;
1934  if (q != n)
1935  ga_throw_error(pnode->expr, pnode->pos,
1936  "Divergence operator can only be applied to"
1937  "Fields with qdim (" << q << ") equal to dim ("
1938  << n << ")");
1939  mii.resize(1);
1940  mii[0] = test ? 2 : 1;
1941  break;
1942  }
1943  pnode->t.adjust_sizes(mii);
1944  }
1945  pnode->test_function_type = test;
1946  }
1947  }
1948  break;
1949 
1950  case GA_NODE_PARAMS:
1951 
1952  // Grad and Diff operators
1953  if (child0->node_type == GA_NODE_NAME) {
1954  if (child0->name.compare("Grad") == 0) {
1955  // cout<<"Compute gradient of ";ga_print_node(child1,cout);cout<<endl;
1956  if (pnode->children.size() != 2)
1957  ga_throw_error(pnode->expr, child0->pos,
1958  "Bad number of parameters for Grad operator");
1959  if (ga_node_mark_tree_for_grad(child1, workspace)) {
1960  ga_node_grad(tree, workspace, me, child1);
1961  ga_node_analysis(tree, workspace, pnode->children[1], me,
1962  ref_elt_dim, eval_fixed_size, ignore_X, option);
1963  child1 = pnode->children[1];
1964  } else {
1965  mi = child1->t.sizes(); mi.push_back(meshdim);
1966  child1->t.adjust_sizes(mi);
1967  child1->node_type = GA_NODE_ZERO;
1968  gmm::clear(child1->tensor().as_vector());
1969  tree.clear_children(child1);
1970  }
1971  tree.replace_node_by_child(pnode, 1);
1972  pnode = child1;
1973  } else if (child0->name.compare("Diff") == 0) {
1974 
1975  if (pnode->children.size() != 3 && pnode->children.size() != 4)
1976  ga_throw_error(pnode->expr, child0->pos,
1977  "Bad number of parameters for Diff operator");
1978  pga_tree_node child2 = pnode->children[2];
1979  if (child2->node_type != GA_NODE_VAL)
1980  ga_throw_error(pnode->expr, child2->pos, "Second parameter of "
1981  "Diff operator has to be a variable name");
1982  std::string vardiff = child2->name;
1983  size_type order = child1->test_function_type;
1984  if (order > 1)
1985  ga_throw_error(pnode->expr, child2->pos, "Cannot derive further "
1986  "this order two expression");
1987 
1988  if (ga_node_mark_tree_for_variable(child1, workspace, me,
1989  vardiff, "", true)) {
1990  ga_node_derivation(tree, workspace, me, child1,
1991  vardiff, "", order+1, true);
1992  child1 = pnode->children[1];
1993  ga_node_analysis(tree, workspace, child1, me, ref_elt_dim,
1994  eval_fixed_size, ignore_X, option);
1995  child1 = pnode->children[1];
1996  } else {
1997  mi = child1->t.sizes(); mi.insert(mi.begin(), 2);
1998  child1->t.adjust_sizes(mi);
1999  child1->node_type = GA_NODE_ZERO;
2000  child1->test_function_type = order ? 3 : 1;
2001  gmm::clear(child1->tensor().as_vector());
2002  tree.clear_children(child1);
2003  }
2004  if (pnode->children.size() == 4) {
2005  ga_tree grad_expr, hess_expr;
2006  ga_node_expand_expression_in_place_of_test
2007  (tree, workspace, pnode->children[1], vardiff, pnode->children[3],
2008  grad_expr, hess_expr, order+1, me, ref_elt_dim, eval_fixed_size,
2009  ignore_X, option);
2010  ga_node_analysis(tree, workspace, pnode->children[1], me,
2011  ref_elt_dim, eval_fixed_size, ignore_X, option);
2012  }
2013  child1 = pnode->children[1];
2014  tree.replace_node_by_child(pnode, 1);
2015  pnode = child1;
2016  } else
2017  ga_throw_error(pnode->expr, child0->pos, "Unknown special operator");
2018  }
2019 
2020  // X (current coordinates on the mesh)
2021  else if (child0->node_type == GA_NODE_X) {
2022  child0->init_scalar_tensor(0);
2023  if (pnode->children.size() != 2)
2024  ga_throw_error(pnode->expr, child1->pos,
2025  "X stands for the coordinates on "
2026  "the real elements. It accepts only one index.");
2027  if (!(child1->node_type == GA_NODE_CONSTANT) ||
2028  child1->tensor().size() != 1)
2029  ga_throw_error(pnode->expr, child1->pos,
2030  "Index for X has to be constant and of size 1.");
2031  child0->nbc1 = size_type(round(child1->tensor()[0]));
2032  if (child0->nbc1 == 0 || child0->nbc1 > meshdim)
2033  ga_throw_error(pnode->expr, child1->pos,"Index for X not convenient. "
2034  "Found " << child0->nbc1 << " with meshdim = "
2035  << meshdim);
2036  tree.replace_node_by_child(pnode, 0);
2037  pnode = child0;
2038  }
2039 
2040  // Reshape operation
2041  else if (child0->node_type == GA_NODE_RESHAPE) {
2042  if (pnode->children.size() < 3)
2043  ga_throw_error(pnode->expr, child1->pos,
2044  "Not enough parameters for Reshape");
2045  if (pnode->children.size() > 12)
2046  ga_throw_error(pnode->expr, child1->pos,
2047  "Too many parameters for Reshape");
2048  pnode->t = child1->t;
2049  pnode->test_function_type = child1->test_function_type;
2050  pnode->name_test1 = child1->name_test1;
2051  pnode->name_test2 = child1->name_test2;
2052  pnode->interpolate_name_test1 = child1->interpolate_name_test1;
2053  pnode->interpolate_name_test2 = child1->interpolate_name_test2;
2054  pnode->qdim1 = child1->qdim1;
2055  pnode->qdim2 = child1->qdim2;
2056  mi.resize(0);
2057  for (size_type i = 0; i < pnode->nb_test_functions(); ++i)
2058  mi.push_back(size1[i]);
2059 
2060  for (size_type i = 2; i < pnode->children.size(); ++i) {
2061  if (pnode->children[i]->node_type != GA_NODE_CONSTANT)
2062  ga_throw_error(pnode->expr, pnode->children[i]->pos,"Reshape sizes "
2063  "should be constant positive integers.");
2064  mi.push_back(size_type(round(pnode->children[i]->tensor()[0])));
2065  if (mi.back() == 0)
2066  ga_throw_error(pnode->expr, pnode->children[i]->pos,
2067  "Wrong zero size for Reshape.");
2068  }
2069  size_type total_size = 1;
2070  for (size_type i = pnode->nb_test_functions(); i < mi.size(); ++i)
2071  total_size *= mi[i];
2072  if (total_size != pnode->tensor_proper_size())
2073  ga_throw_error(pnode->expr, pnode->pos,"Invalid sizes for reshape, "
2074  "found a total of " << total_size << " should be " <<
2075  pnode->tensor_proper_size() << ".");
2076  pnode->t.adjust_sizes(mi);
2077 
2078  if (child1->node_type == GA_NODE_CONSTANT) {
2079  pnode->node_type = GA_NODE_CONSTANT;
2080  tree.clear_children(pnode);
2081  } else if (child1->node_type == GA_NODE_ZERO) {
2082  pnode->node_type = GA_NODE_ZERO;
2083  tree.clear_children(pnode);
2084  }
2085  }
2086 
2087  // Cross product of two vectors
2088  else if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
2089  if (pnode->children.size() != 3)
2090  ga_throw_error(pnode->expr, child1->pos,
2091  "Wrong number of parameters for Cross_product");
2092  pga_tree_node child2 = pnode->children[2];
2093 
2094  if (false && child1->is_constant() && child2->is_constant()) {
2095  pnode->node_type = GA_NODE_CONSTANT;
2096  pnode->test_function_type = 0;
2097  if (child1->tensor_proper_size() != 3 ||
2098  child2->tensor_proper_size() != 3)
2099  ga_throw_error(pnode->expr, child1->pos, "Cross_product is only "
2100  "defined on three-dimensional vectors");
2101  pnode->t = child1->t;
2102  base_tensor &t0 = pnode->tensor();
2103  base_tensor &t1 = child1->tensor(), &t2 = child2->tensor();
2104  t0[0] = t1[1]*t2[2] - t1[2]*t2[1];
2105  t0[1] = t1[2]*t2[0] - t1[0]*t2[2];
2106  t0[2] = t1[0]*t2[1] - t1[1]*t2[0];
2107  if (pnode->tensor_is_zero())
2108  pnode->node_type = GA_NODE_ZERO;
2109  tree.clear_children(pnode);
2110  } else {
2111  pnode->mult_test(child1, child2);
2112  mi = pnode->t.sizes();
2113  mi.push_back(3);
2114  pnode->t.adjust_sizes(mi);
2115  }
2116  }
2117 
2118  // Swap_indices operation
2119  else if (child0->node_type == GA_NODE_SWAP_IND) {
2120  if (pnode->children.size() != 4)
2121  ga_throw_error(pnode->expr, child1->pos,
2122  "Wrong number of parameters for Swap_indices");
2123  pnode->t = child1->t;
2124  pnode->test_function_type = child1->test_function_type;
2125  pnode->name_test1 = child1->name_test1;
2126  pnode->name_test2 = child1->name_test2;
2127  pnode->interpolate_name_test1 = child1->interpolate_name_test1;
2128  pnode->interpolate_name_test2 = child1->interpolate_name_test2;
2129  pnode->qdim1 = child1->qdim1;
2130  pnode->qdim2 = child1->qdim2;
2131  size_type ind[4];
2132 
2133  for (size_type i = 2; i < 4; ++i) {
2134  if (pnode->children[i]->node_type != GA_NODE_CONSTANT ||
2135  pnode->children[i]->tensor().size() != 1)
2136  ga_throw_error(pnode->expr, pnode->children[i]->pos, "Indices "
2137  "for swap should be constant positive integers.");
2138  ind[i] = size_type(round(pnode->children[i]->tensor()[0]));
2139  if (ind[i] < 1 || ind[i] > child1->tensor_proper_size())
2140  ga_throw_error(pnode->expr, pnode->children[i]->pos, "Index "
2141  "out of range for Swap_indices.");
2142  ind[i]--;
2143  }
2144  if (ind[2] == ind[3]) {
2145  tree.replace_node_by_child(pnode, 1);
2146  pnode = child1;
2147  } else {
2148  mi = pnode->tensor().sizes();
2149  size_type nbtf = child1->nb_test_functions();
2150  std::swap(mi[ind[2]+nbtf], mi[ind[3]+nbtf]);
2151  pnode->t.adjust_sizes(mi);
2152 
2153  if (child1->node_type == GA_NODE_CONSTANT) {
2154  pnode->node_type = GA_NODE_CONSTANT;
2155  if (ind[2] > ind[3]) std::swap(ind[2], ind[3]);
2156  size_type ii1 = 1, ii2 = 1, ii3 = 1;
2157  for (size_type i = 0; i < child1->tensor_order(); ++i) {
2158  if (i<ind[2]) ii1 *= child1->tensor_proper_size(i);
2159  if (i>ind[2] && i<ind[3]) ii2 *= child1->tensor_proper_size(i);
2160  if (i>ind[3]) ii3 *= child1->tensor_proper_size(i);
2161  }
2162  size_type nn1 = child1->tensor_proper_size(ind[2]);
2163  size_type nn2 = child1->tensor_proper_size(ind[3]);
2164  auto it = pnode->tensor().begin();
2165  for (size_type i = 0; i < ii3; ++i)
2166  for (size_type j = 0; j < nn1; ++j)
2167  for (size_type k = 0; k < ii2; ++k)
2168  for (size_type l = 0; l < nn2; ++l)
2169  for (size_type m = 0; m < ii1; ++m, ++it)
2170  *it = child0->tensor()[m+j*ii1+k*ii1*nn1+l*ii1*nn1*ii2+
2171  i*ii1*nn1*ii2*nn2];
2172  tree.clear_children(pnode);
2173  } else if (child1->node_type == GA_NODE_ZERO) {
2174  pnode->node_type = GA_NODE_ZERO;
2175  tree.clear_children(pnode);
2176  }
2177  }
2178  }
2179 
2180  // Index_move_last operation
2181  else if (child0->node_type == GA_NODE_IND_MOVE_LAST) {
2182  if (pnode->children.size() != 3)
2183  ga_throw_error(pnode->expr, child1->pos,
2184  "Wrong number of parameters for Index_move_last");
2185  pnode->t = child1->t;
2186  pnode->test_function_type = child1->test_function_type;
2187  pnode->name_test1 = child1->name_test1;
2188  pnode->name_test2 = child1->name_test2;
2189  pnode->interpolate_name_test1 = child1->interpolate_name_test1;
2190  pnode->interpolate_name_test2 = child1->interpolate_name_test2;
2191  pnode->qdim1 = child1->qdim1;
2192  pnode->qdim2 = child1->qdim2;
2193  size_type ind;
2194 
2195  if (pnode->children[2]->node_type != GA_NODE_CONSTANT ||
2196  pnode->children[2]->tensor().size() != 1)
2197  ga_throw_error(pnode->expr, pnode->children[2]->pos, "Index for "
2198  "Index_move_last should be constant positive integers.");
2199  ind = size_type(round(pnode->children[2]->tensor()[0]));
2200  if (ind < 1 || ind > child1->tensor_proper_size())
2201  ga_throw_error(pnode->expr, pnode->children[2]->pos, "Index "
2202  "out of range for Index_move_last.");
2203 
2204  if (ind == child1->tensor_order()) {
2205  tree.replace_node_by_child(pnode, 1);
2206  pnode = child1;
2207  } else {
2208  mi = pnode->tensor().sizes();
2209  size_type nbtf = child1->nb_test_functions();
2210  for (size_type i = ind; i < child1->tensor_order(); ++i)
2211  std::swap(mi[i-1+nbtf], mi[i+nbtf]);
2212  pnode->t.adjust_sizes(mi);
2213 
2214  if (child1->node_type == GA_NODE_CONSTANT) {
2215  pnode->node_type = GA_NODE_CONSTANT;
2216  ind--;
2217  size_type ii1 = 1, ii2 = 1;
2218  for (size_type i = 0; i < child1->tensor_order(); ++i) {
2219  if (i<ind) ii1 *= child1->tensor_proper_size(i);
2220  if (i>ind) ii2 *= child1->tensor_proper_size(i);
2221  }
2222  size_type nn = child1->tensor_proper_size(ind);
2223  auto it = pnode->tensor().begin();
2224  for (size_type i = 0; i < nn; ++i)
2225  for (size_type j = 0; j < ii2; ++j)
2226  for (size_type k = 0; k < ii1; ++k, ++it)
2227  *it = child0->tensor()[k+i*ii1+j*ii1*nn];
2228  tree.clear_children(pnode);
2229  } else if (child1->node_type == GA_NODE_ZERO) {
2230  pnode->node_type = GA_NODE_ZERO;
2231  tree.clear_children(pnode);
2232  }
2233  }
2234  }
2235 
2236  // Tensor contraction operator
2237  else if (child0->node_type == GA_NODE_CONTRACT) {
2238  std::vector<size_type> ind(2), indsize(2);
2239  if (pnode->children.size() == 4)
2240  { ind[0] = 2; ind[1] = 3; }
2241  else if (pnode->children.size() == 5)
2242  { ind[0] = 2; ind[1] = 4; }
2243  else if (pnode->children.size() == 7) {
2244  ind.resize(4);
2245  indsize.resize(4);
2246  ind[0] = 2; ind[1] = 3; ind[2] = 5; ind[3] = 6;
2247  }
2248 
2249  size_type order = 0, kk = 0, ll = 1; // Indices extraction and controls
2250  for (size_type i = 1; i < pnode->children.size(); ++i) {
2251  if (i == ind[kk]) {
2252  if (pnode->children[i]->node_type != GA_NODE_CONSTANT ||
2253  pnode->children[i]->tensor().size() != 1)
2254  ga_throw_error(pnode->children[i]->expr, pnode->children[i]->pos,
2255  "Invalid parameter for Contract. "
2256  "Should be an index number.");
2257  ind[kk] = size_type(round(pnode->children[i]->tensor()[0]));
2258  order = pnode->children[ll]->tensor_order();
2259  if (ind[kk] < 1 || ind[kk] > order)
2260  ga_throw_error(pnode->children[i]->expr, pnode->children[i]->pos,
2261  "Parameter out of range for Contract (should be "
2262  "between 1 and " << order << ")");
2263  ind[kk]--;
2264  indsize[kk] = pnode->children[ll]->tensor_proper_size(ind[kk]);
2265  if (kk >= ind.size()/2 && indsize[kk] != indsize[kk-ind.size()/2])
2266  ga_throw_error(child0->expr, child0->pos,
2267  "Invalid parameters for Contract. Cannot "
2268  "contract on indices of different sizes");
2269  ++kk;
2270  } else ll = i;
2271  }
2272 
2273  if (pnode->children.size() == 4) {
2274  // Contraction of a single tensor on a single pair of indices
2275  pnode->test_function_type = child1->test_function_type;
2276  pnode->name_test1 = child1->name_test1;
2277  pnode->name_test2 = child1->name_test2;
2278  pnode->interpolate_name_test1 = child1->interpolate_name_test1;
2279  pnode->interpolate_name_test2 = child1->interpolate_name_test2;
2280  pnode->qdim1 = child1->qdim1;
2281  pnode->qdim2 = child1->qdim2;
2282 
2283  size_type i1 = ind[0], i2 = ind[1];
2284  if (i1 == i2)
2285  ga_throw_error(child0->expr, child0->pos,
2286  "Invalid parameters for Contract. Repeated index.");
2287 
2288  mi.resize(0);
2289  for (size_type i = 0; i < pnode->nb_test_functions(); ++i)
2290  mi.push_back(size1[i]);
2291  for (size_type i = 0; i < order; ++i)
2292  if (i != i1 && i != i2)
2293  mi.push_back(child1->tensor_proper_size(i));
2294  pnode->t.adjust_sizes(mi);
2295 
2296  if (child1->node_type == GA_NODE_CONSTANT) {
2297 
2298  if (i1 > i2) std::swap(i1, i2);
2299  size_type ii1 = 1, ii2 = 1, ii3 = 1;
2300  size_type nn = indsize[0];
2301  for (size_type i = 0; i < order; ++i) {
2302  if (i < i1) ii1 *= child1->tensor_proper_size(i);
2303  if (i > i1 && i < i2) ii2 *= child1->tensor_proper_size(i);
2304  if (i > i2) ii3 *= child1->tensor_proper_size(i);
2305  }
2306 
2307  auto it = pnode->tensor().begin();
2308  for (size_type i = 0; i < ii3; ++i)
2309  for (size_type j = 0; j < ii2; ++j)
2310  for (size_type k = 0; k < ii1; ++k, ++it) {
2311  *it = scalar_type(0);
2312  size_type pre_ind = k+j*ii1*nn+i*ii1*nn*ii2*nn;
2313  for (size_type n = 0; n < nn; ++n)
2314  *it += child1->tensor()[pre_ind+n*ii1+n*ii1*nn*ii2];
2315  }
2316 
2317  pnode->node_type = GA_NODE_CONSTANT;
2318  tree.clear_children(pnode);
2319  } else if (child1->node_type == GA_NODE_ZERO) {
2320  pnode->node_type = GA_NODE_ZERO;
2321  tree.clear_children(pnode);
2322  }
2323 
2324  } else if (pnode->children.size() == 5) {
2325  // Contraction of two tensors on a single pair of indices
2326  pga_tree_node child2 = pnode->children[3];
2327  pnode->mult_test(child1, child2);
2328 
2329  size_type i1 = ind[0], i2 = ind[1];
2330  mi = pnode->tensor().sizes();
2331  for (size_type i = 0; i < dim1; ++i)
2332  if (i != i1) mi.push_back(child1->tensor_proper_size(i));
2333  for (size_type i = 0; i < order; ++i)
2334  if (i != i2) mi.push_back(child2->tensor_proper_size(i));
2335  pnode->t.adjust_sizes(mi);
2336 
2337  if (child1->node_type == GA_NODE_CONSTANT &&
2338  child2->node_type == GA_NODE_CONSTANT) {
2339  size_type ii1 = 1, ii2 = 1, ii3 = 1, ii4 = 1;
2340  size_type nn = indsize[0];
2341  for (size_type i = 0; i < dim1; ++i) {
2342  if (i < i1) ii1 *= child1->tensor_proper_size(i);
2343  if (i > i1) ii2 *= child1->tensor_proper_size(i);
2344  }
2345  for (size_type i = 0; i < order; ++i) {
2346  if (i < i2) ii3 *= child2->tensor_proper_size(i);
2347  if (i > i2) ii4 *= child2->tensor_proper_size(i);
2348  }
2349 
2350  auto it = pnode->tensor().begin();
2351  for (size_type i = 0; i < ii4; ++i)
2352  for (size_type j = 0; j < ii3; ++j)
2353  for (size_type k = 0; k < ii2; ++k)
2354  for (size_type l = 0; l < ii1; ++l, ++it) {
2355  *it = scalar_type(0);
2356  for (size_type n = 0; n < nn; ++n)
2357  *it += child1->tensor()[l+n*ii1+k*ii1*nn]
2358  * child2->tensor()[j+n*ii3+i*ii3*nn];
2359  }
2360 
2361  } else if (child1->tensor_is_zero() || child2->tensor_is_zero()) {
2362  pnode->node_type = GA_NODE_ZERO;
2363  tree.clear_children(pnode);
2364  }
2365 
2366  } else if (pnode->children.size() == 7) {
2367  // Contraction of two tensors on two pairs of indices
2368  pga_tree_node child2 = pnode->children[4];
2369  pnode->mult_test(child1, child2);
2370  if (ind[0] == ind[1] || ind[2] == ind[3])
2371  ga_throw_error(child0->expr, child0->pos,
2372  "Invalid parameters for Contract. Repeated index.");
2373 
2374  size_type i1 = ind[0], i2 = ind[1], i3 = ind[2], i4 = ind[3];
2375  mi = pnode->tensor().sizes();
2376  for (size_type i = 0; i < dim1; ++i)
2377  if (i != i1 && i != i2) mi.push_back(child1->tensor_proper_size(i));
2378  for (size_type i = 0; i < order; ++i)
2379  if (i != i3 && i != i4) mi.push_back(child2->tensor_proper_size(i));
2380  pnode->t.adjust_sizes(mi);
2381 
2382  if (child1->tensor_is_zero() || child2->tensor_is_zero()) {
2383  pnode->node_type = GA_NODE_ZERO;
2384  tree.clear_children(pnode);
2385  } else if (child1->node_type == GA_NODE_CONSTANT &&
2386  child2->node_type == GA_NODE_CONSTANT) {
2387  size_type nn1 = indsize[0], nn2 = indsize[1];
2388  if (i1 > i2)
2389  { std::swap(i1, i2); std::swap(i3, i4); std::swap(nn1, nn2); }
2390 
2391  size_type ii1 = 1, ii2 = 1, ii3 = 1, ii4 = 1, ii5 = 1, ii6 = 1;
2392  for (size_type i = 0; i < dim1; ++i) {
2393  if (i < i1) ii1 *= child1->tensor_proper_size(i);
2394  if (i > i1 && i < i2) ii2 *= child1->tensor_proper_size(i);
2395  if (i > i2) ii3 *= child1->tensor_proper_size(i);
2396  }
2397  for (size_type i = 0; i < order; ++i) {
2398  if (i < i3 && i < i4) ii4 *= child2->tensor_proper_size(i);
2399  if ((i > i3 && i < i4) || (i > i4 && i < i3))
2400  ii5 *= child2->tensor_proper_size(i);
2401  if (i > i3 && i > i4) ii6 *= child2->tensor_proper_size(i);
2402  }
2403 
2404  auto it = pnode->tensor().begin();
2405  size_type m1 = ii4, m2 = ii4*nn1*ii5;
2406  if (i3 < i4) std::swap(m1, m2);
2407  for (size_type i = 0; i < ii6; ++i)
2408  for (size_type j = 0; j < ii5; ++j)
2409  for (size_type k = 0; k < ii4; ++k)
2410  for (size_type l = 0; l < ii3; ++l)
2411  for (size_type m = 0; m < ii2; ++m)
2412  for (size_type p = 0; p < ii1; ++p, ++it) {
2413  *it = scalar_type(0);
2414  size_type q1 = p+m*ii1*nn1+l*ii1*nn1*ii2*nn2;
2415  size_type q2 = k+j*ii4*nn1+i*ii4*nn1*ii5*nn2;
2416  for (size_type n1 = 0; n1 < nn1; ++n1)
2417  for (size_type n2 = 0; n2 < nn2; ++n2)
2418  *it += child1->tensor()[q1+n1*ii1+n2*ii1*nn1*ii2]
2419  * child2->tensor()[q2+n1*m1+n2*m2];
2420  }
2421  }
2422  } else
2423  ga_throw_error(pnode->expr, child1->pos,
2424  "Wrong number of parameters for Contract");
2425  }
2426 
2427  // Evaluation of a predefined function
2428  else if (child0->node_type == GA_NODE_PREDEF_FUNC) {
2429 
2430  for (size_type i = 1; i < pnode->children.size(); ++i)
2431  ga_valid_operand(pnode->children[i]);
2432  std::string name = child0->name;
2433  ga_predef_function_tab::const_iterator it = PREDEF_FUNCTIONS.find(name);
2434  const ga_predef_function &F = it->second;
2435  size_type nbargs = F.nbargs();
2436  if (nbargs+1 != pnode->children.size()) {
2437  ga_throw_error(pnode->expr, pnode->pos, "Bad number of arguments "
2438  "for predefined function " << name << ". Found "
2439  << pnode->children.size()-1 << ", should be "<<nbargs << ".");
2440  }
2441  pnode->test_function_type = 0;
2442  pga_tree_node child2 = (nbargs == 2) ? pnode->children[2] : child1;
2443  all_cte = child1->node_type == GA_NODE_CONSTANT;
2444  if (nbargs == 2)
2445  all_cte = all_cte && (child2->node_type == GA_NODE_CONSTANT);
2446  if (child1->test_function_type || child2->test_function_type)
2447  ga_throw_error(pnode->expr, pnode->pos, "Test functions cannot be "
2448  "passed as argument of a predefined function.");
2449  // if (child1->tensor_order() > 2 || child2->tensor_order() > 2)
2450  // ga_throw_error(pnode->expr, pnode->pos, "Sorry, function can be "
2451  // "applied to scalar, vector and matrices only.");
2452  size_type s1 = child1->tensor().size();
2453  size_type s2 = (nbargs == 2) ? child2->tensor().size() : s1;
2454  if (s1 != s2 && (s1 != 1 || s2 != 1))
2455  ga_throw_error(pnode->expr, pnode->pos,
2456  "Invalid argument size for a scalar function. "
2457  "Size of first argument: " << s1 <<
2458  ". Size of second argument: " << s2 << ".");
2459 
2460  if (nbargs == 1) {
2461  pnode->t = child1->t;
2462  } else {
2463  if (s1 == s2) {
2464  pnode->t = child1->t;
2465  } else if (s1 == 1) {
2466  pnode->t = child2->t;
2467  } else {
2468  pnode->t = child1->t;
2469  }
2470  }
2471 
2472  if (all_cte) {
2473  if (pnode->der1)
2474  GMM_ASSERT1(false, "Sorry, to be done");
2475  pnode->node_type = GA_NODE_CONSTANT;
2476  if (nbargs == 1) {
2477  for (size_type i = 0; i < s1; ++i)
2478  pnode->tensor()[i] = F(child1->tensor()[i]);
2479  } else {
2480  if (s1 == s2) {
2481  for (size_type i = 0; i < s1; ++i)
2482  pnode->tensor()[i] = F(child1->tensor()[i],
2483  child2->tensor()[i]);
2484  } else if (s1 == 1) {
2485  for (size_type i = 0; i < s2; ++i)
2486  pnode->tensor()[i] = F(child1->tensor()[0],
2487  child2->tensor()[i]);
2488  } else {
2489  for (size_type i = 0; i < s1; ++i)
2490  pnode->tensor()[i] = F(child1->tensor()[i],
2491  child2->tensor()[0]);
2492  }
2493  }
2494  tree.clear_children(pnode);
2495  }
2496  }
2497 
2498  // Special constant functions: meshdim, qdim(u) ...
2499  else if (child0->node_type == GA_NODE_SPEC_FUNC) {
2500 
2501  for (size_type i = 1; i < pnode->children.size(); ++i)
2502  ga_valid_operand(pnode->children[i]);
2503  if (pnode->children.size() != 2)
2504  ga_throw_error(pnode->expr, pnode->pos,
2505  "One and only one argument is allowed for function "
2506  +child0->name+".");
2507 
2508  if (!(child0->name.compare("qdim"))) {
2509  if (child1->node_type != GA_NODE_VAL)
2510  ga_throw_error(pnode->expr, pnode->pos, "The argument of qdim "
2511  "function can only be a variable name.");
2512  pnode->node_type = GA_NODE_CONSTANT;
2513  pnode->init_scalar_tensor(scalar_type(workspace.qdim(child1->name)));
2514  if (pnode->tensor()[0] <= 0)
2515  ga_throw_error(pnode->expr, pnode->pos,
2516  "Invalid null size of variable");
2517  } else if (!(child0->name.compare("qdims"))) {
2518  if (child1->node_type != GA_NODE_VAL)
2519  ga_throw_error(pnode->expr, pnode->pos, "The argument of qdim "
2520  "function can only be a variable name.");
2521  pnode->node_type = GA_NODE_CONSTANT;
2522  bgeot::multi_index mii = workspace.qdims(child1->name);
2523  if (mii.size() > 6)
2524  ga_throw_error(pnode->expr, pnode->pos,
2525  "Tensor with too much dimensions. Limited to 6");
2526  if (mii.size() == 0 || scalar_type(mii[0]) <= 0)
2527  ga_throw_error(pnode->expr, pnode->pos,
2528  "Invalid null size of variable");
2529  if (mii.size() == 1)
2530  pnode->init_scalar_tensor(scalar_type(mii[0]));
2531  if (mii.size() >= 1) {
2532  pnode->init_vector_tensor(mii.size());
2533  for (size_type i = 0; i < mii.size(); ++i)
2534  pnode->tensor()[i] = scalar_type(mii[i]);
2535  }
2536  } else if (!(child0->name.compare("Id"))) {
2537  bool valid = (child1->node_type == GA_NODE_CONSTANT);
2538  int n = valid ? int(round(child1->tensor()[0])) : -1;
2539  if (n <= 0 || n > 100 || child1->tensor_order() > 0)
2540  ga_throw_error(pnode->expr, pnode->pos, "The argument of Id "
2541  "should be a (small) positive integer.");
2542  pnode->node_type = GA_NODE_CONSTANT;
2543  if (n == 1)
2544  pnode->init_scalar_tensor(scalar_type(1));
2545  else {
2546  pnode->init_matrix_tensor(n,n);
2547  gmm::clear(pnode->tensor().as_vector());
2548  for (int i = 0; i < n; ++i) pnode->tensor()(i,i) = scalar_type(1);
2549  }
2550  } else ga_throw_error(pnode->expr, pnode->children[0]->pos,
2551  "Unknown special function.");
2552  tree.clear_children(pnode);
2553  }
2554 
2555  // Nonlinear operator call
2556  else if (child0->node_type == GA_NODE_OPERATOR) {
2557 
2558  for (size_type i = 1; i < pnode->children.size(); ++i)
2559  ga_valid_operand(pnode->children[i]);
2560  all_cte = true;
2561  ga_nonlinear_operator::arg_list args;
2562  for (size_type i = 1; i < pnode->children.size(); ++i) {
2563  all_cte = all_cte
2564  && (pnode->children[i]->node_type == GA_NODE_CONSTANT);
2565  args.push_back(&(pnode->children[i]->tensor()));
2566  if (pnode->children[i]->node_type == GA_NODE_ALLINDICES)
2567  ga_throw_error(pnode->expr, pnode->children[i]->pos,
2568  "Colon operator is not allowed in nonlinear "
2569  "operator call.");
2570  if (pnode->children[i]->test_function_type)
2571  ga_throw_error(pnode->expr, pnode->pos, "Test functions cannot be "
2572  "passed as argument of a nonlinear operator.");
2573  if (pnode->children[i]->tensor_order() > 2)
2574  ga_throw_error(pnode->expr, pnode->pos,
2575  "Sorry, arguments to nonlinear operators should "
2576  "only be scalar, vector or matrices");
2577  }
2578  ga_predef_operator_tab::T::const_iterator it
2579  = PREDEF_OPERATORS.tab.find(child0->name);
2580  const ga_nonlinear_operator &OP = *(it->second);
2581  mi.resize(0);
2582  if (!(OP.result_size(args, mi)))
2583  ga_throw_error(pnode->expr, pnode->pos,
2584  "Wrong number or wrong size of arguments for the "
2585  "call of nonlinear operator " + child0->name);
2586 
2587  pnode->test_function_type = 0;
2588 
2589  if (child0->der1 > args.size() || child0->der2 > args.size())
2590  ga_throw_error(pnode->expr, child0->pos,
2591  "Invalid derivative number for nonlinear operator "
2592  + child0->name);
2593 
2594  if (child0->der1 && child0->der2 == 0) {
2595  for (size_type i = 0; i < args[child0->der1-1]->sizes().size(); ++i)
2596  mi.push_back(args[child0->der1-1]->sizes()[i]);
2597  pnode->t.adjust_sizes(mi);
2598  if (all_cte) {
2599  pnode->node_type = GA_NODE_CONSTANT;
2600  OP.derivative(args, child0->der1, pnode->tensor());
2601  tree.clear_children(pnode);
2602  }
2603  } else if (child0->der1 && child0->der2) {
2604  for (size_type i = 0; i < args[child0->der1-1]->sizes().size(); ++i)
2605  mi.push_back(args[child0->der1-1]->sizes()[i]);
2606  for (size_type i = 0; i < args[child0->der2-1]->sizes().size(); ++i)
2607  mi.push_back(args[child0->der2-1]->sizes()[i]);
2608  pnode->t.adjust_sizes(mi);
2609  if (all_cte) {
2610  pnode->node_type = GA_NODE_CONSTANT;
2611  OP.second_derivative(args, child0->der1, child0->der2,
2612  pnode->tensor());
2613  tree.clear_children(pnode);
2614  }
2615  } else {
2616  pnode->t.adjust_sizes(mi);
2617  if (all_cte) {
2618  pnode->node_type = GA_NODE_CONSTANT;
2619  OP.value(args, pnode->tensor());
2620  tree.clear_children(pnode);
2621  }
2622  }
2623  }
2624 
2625  // Access to components of a tensor
2626  else {
2627  all_cte = (child0->node_type == GA_NODE_CONSTANT);
2628  if (pnode->children.size() != child0->tensor_order() + 1)
2629  ga_throw_error(pnode->expr, pnode->pos, "Bad number of indices.");
2630  for (size_type i = 1; i < pnode->children.size(); ++i)
2631  if (pnode->children[i]->node_type != GA_NODE_ALLINDICES &&
2632  (pnode->children[i]->node_type != GA_NODE_CONSTANT ||
2633  pnode->children[i]->tensor().size() != 1))
2634  ga_throw_error(pnode->expr, pnode->children[i]->pos,
2635  "Indices should be constant integers or colon.");
2636 
2637  bgeot::multi_index mi1(size0.size()), mi2, indices;
2638  for (size_type i = 0; i < child0->tensor_order(); ++i) {
2639  if (pnode->children[i+1]->node_type == GA_NODE_ALLINDICES) {
2640  mi2.push_back(child0->tensor_proper_size(i));
2641  indices.push_back(i);
2642  mi1[i] = 0;
2643  } else {
2644  mi1[i] = size_type(round(pnode->children[i+1]->tensor()[0])-1);
2645  if (mi1[i] >= child0->tensor_proper_size(i))
2646  ga_throw_error(pnode->expr, pnode->children[i+1]->pos,
2647  "Index out of range, " << mi1[i]+1
2648  << ". Should be between 1 and "
2649  << child0->tensor_proper_size(i) << ".");
2650  }
2651  }
2652  mi.resize(0);
2653  for (size_type i = 0; i < child0->nb_test_functions(); ++i)
2654  mi.push_back(child0->t.sizes()[i]);
2655  for (size_type i = 0; i < mi2.size(); ++i) mi.push_back(mi2[i]);
2656  pnode->t.adjust_sizes(mi);
2657  pnode->test_function_type = child0->test_function_type;
2658  pnode->name_test1 = child0->name_test1;
2659  pnode->name_test2 = child0->name_test2;
2660  pnode->interpolate_name_test1 = child0->interpolate_name_test1;
2661  pnode->interpolate_name_test2 = child0->interpolate_name_test2;
2662  pnode->qdim1 = child0->qdim1;
2663  pnode->qdim2 = child0->qdim2;
2664 
2665  if (child0->tensor_is_zero()) {
2666  gmm::clear(pnode->tensor().as_vector());
2667  pnode->node_type = GA_NODE_ZERO;
2668  tree.clear_children(pnode);
2669  } else if (all_cte) {
2670  pnode->node_type = GA_NODE_CONSTANT;
2671  for (bgeot::multi_index mi3(mi2.size()); !mi3.finished(mi2);
2672  mi3.incrementation(mi2)) {
2673  for (size_type j = 0; j < mi2.size(); ++j) {
2674  mi1[indices[j]] = mi3[j];
2675  }
2676  pnode->tensor()(mi3) = pnode->children[0]->tensor()(mi1);
2677  }
2678  tree.clear_children(pnode);
2679  }
2680  }
2681  break;
2682 
2683  default:GMM_ASSERT1(false, "Unexpected node type " << pnode->node_type
2684  << " in semantic analysis. Internal error.");
2685  }
2686  pnode->hash_value = ga_hash_code(pnode);
2687  for (size_type i = 0; i < pnode->children.size(); ++i) {
2688  pnode->hash_value += (pnode->children[i]->hash_value)
2689  * 1.0101*(pnode->symmetric_op ? scalar_type(1) : scalar_type(1+i));
2690  }
2691  pnode->hash_value = sin(1.2003*pnode->hash_value);
2692  }
2693 
2694 
2695  void ga_semantic_analysis(ga_tree &tree,
2696  const ga_workspace &workspace,
2697  const mesh &m,
2698  size_type ref_elt_dim,
2699  bool eval_fixed_size,
2700  bool ignore_X, int option) {
2701  GMM_ASSERT1(predef_operators_nonlinear_elasticity_initialized &&
2702  predef_operators_plasticity_initialized &&
2703  predef_operators_contact_initialized, "Internal error");
2704  if (!(tree.root))
2705  return;
2706  // cout << "Begin semantic analysis with ";
2707  // ga_print_node(tree.root, cout); cout << endl;
2708 
2709  if (option == 1) { workspace.test1.clear(); workspace.test2.clear(); }
2710  ga_node_analysis(tree, workspace, tree.root, m, ref_elt_dim,
2711  eval_fixed_size, ignore_X, option);
2712  if (tree.root && option == 2) {
2713  if (((tree.root->test_function_type & 1) &&
2714  (tree.root->name_test1.compare(workspace.selected_test1.varname)
2715  || tree.root->interpolate_name_test1.compare
2716  (workspace.selected_test1.transname)))
2717  ||
2718  ((tree.root->test_function_type & 2) &&
2719  (tree.root->name_test2.compare(workspace.selected_test2.varname)
2720  || tree.root->interpolate_name_test2.compare
2721  (workspace.selected_test2.transname))))
2722  tree.clear();
2723  }
2724  ga_valid_operand(tree.root);
2725  // cout << "End of semantic analysis";
2726  // if (tree.root)
2727  // ga_print_node(tree.root, cout);
2728  // cout << endl;
2729  }
2730 
2731 
2732  void ga_extract_factor(ga_tree &result_tree, pga_tree_node pnode,
2733  pga_tree_node &new_pnode) {
2734  result_tree.clear();
2735  result_tree.copy_node(pnode, result_tree.root);
2736  new_pnode = result_tree.root;
2737 
2738  bool minus_sign = false;
2739 
2740  pga_tree_node pnode_child = pnode;
2741  pnode = pnode->parent;
2742 
2743  while (pnode) {
2744 
2745  size_type nbch = pnode->children.size();
2746  pga_tree_node child0 = (nbch > 0) ? pnode->children[0] : 0;
2747  pga_tree_node child1 = (nbch > 1) ? pnode->children[1] : 0;
2748 
2749  switch (pnode->node_type) {
2750 
2751  case GA_NODE_OP:
2752  switch(pnode->op_type) {
2753  case GA_PLUS:
2754  // Nothing to do
2755  break;
2756  case GA_MINUS:
2757  if (child1 == pnode_child) minus_sign = !(minus_sign);
2758  // A remaining minus sign is added at the end if necessary.
2759  break;
2760  case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM: case GA_SKEW:
2761  case GA_TRACE: case GA_DEVIATOR: case GA_PRINT:
2762  // Copy of the term
2763  result_tree.insert_node(result_tree.root, pnode->node_type);
2764  result_tree.root->op_type = pnode->op_type;
2765  result_tree.root->pos = pnode->pos;
2766  result_tree.root->expr = pnode->expr;
2767  break;
2768  case GA_DOT: case GA_MULT: case GA_COLON: case GA_TMULT:
2769  case GA_DOTMULT: case GA_DIV: case GA_DOTDIV:
2770  // Copy of the term and other child
2771  result_tree.insert_node(result_tree.root, pnode->node_type);
2772  result_tree.root->op_type = pnode->op_type;
2773  result_tree.root->pos = pnode->pos;
2774  result_tree.root->expr = pnode->expr;
2775  result_tree.root->children.resize(2, nullptr);
2776  if (child0 == pnode_child) {
2777  result_tree.copy_node(child1, result_tree.root->children[1]);
2778  result_tree.root->accept_child(1);
2779  } else if (child1 == pnode_child) {
2780  std::swap(result_tree.root->children[1],
2781  result_tree.root->children[0]);
2782  result_tree.copy_node(child0, result_tree.root->children[0]);
2783  result_tree.root->accept_child(0);
2784  } else GMM_ASSERT1(false, "Corrupted tree");
2785  break;
2786  default: GMM_ASSERT1(false, "Unexpected operation. Internal error.");
2787  }
2788  break;
2789 
2790  case GA_NODE_PARAMS:
2791  if (child0->node_type == GA_NODE_RESHAPE) {
2792  GMM_ASSERT1(child1 == pnode_child, "Cannot extract a factor of a "
2793  "Reshape size parameter");
2794  // Copy of the term and other children
2795  result_tree.insert_node(result_tree.root, pnode->node_type);
2796  result_tree.root->pos = pnode->pos;
2797  result_tree.root->expr = pnode->expr;
2798  result_tree.root->children.resize(pnode->children.size(), nullptr);
2799  std::swap(result_tree.root->children[1],
2800  result_tree.root->children[0]);
2801  for (size_type i = 0; i < pnode->children.size(); ++i)
2802  if (i != 1) {
2803  result_tree.copy_node(pnode->children[i],
2804  result_tree.root->children[i]);
2805  result_tree.root->accept_child(i);
2806  }
2807  } else if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
2808  pga_tree_node child2 = pnode->children[2];
2809  result_tree.insert_node(result_tree.root, pnode->node_type);
2810  result_tree.root->pos = pnode->pos;
2811  result_tree.root->expr = pnode->expr;
2812  result_tree.root->children.resize(3, nullptr);
2813  if (child1 == pnode_child) {
2814  std::swap(result_tree.root->children[1],
2815  result_tree.root->children[0]);
2816  result_tree.copy_node(pnode->children[0],
2817  result_tree.root->children[0]);
2818  result_tree.root->accept_child(0);
2819  result_tree.copy_node(pnode->children[2],
2820  result_tree.root->children[2]);
2821  result_tree.root->accept_child(2);
2822  } else if (child2 == pnode_child) {
2823  std::swap(result_tree.root->children[2],
2824  result_tree.root->children[0]);
2825  result_tree.copy_node(pnode->children[0],
2826  result_tree.root->children[0]);
2827  result_tree.root->accept_child(0);
2828  result_tree.copy_node(pnode->children[1],
2829  result_tree.root->children[1]);
2830  result_tree.root->accept_child(1);
2831  } else GMM_ASSERT1(false, "Corrupted tree");
2832  } else if (child0->node_type == GA_NODE_SWAP_IND) {
2833  GMM_ASSERT1(child1 == pnode_child, "Cannot extract a factor of a "
2834  "Swap_indices size parameter");
2835  // Copy of the term and other children
2836  result_tree.insert_node(result_tree.root, pnode->node_type);
2837  result_tree.root->pos = pnode->pos;
2838  result_tree.root->expr = pnode->expr;
2839  result_tree.root->children.resize(pnode->children.size(), nullptr);
2840  for (size_type i = 0; i < pnode->children.size(); ++i)
2841  if (pnode->children[i] == pnode_child)
2842  std::swap(result_tree.root->children[i],
2843  result_tree.root->children[0]);
2844  for (size_type i = 0; i < pnode->children.size(); ++i)
2845  if (pnode->children[i] != pnode_child) {
2846  result_tree.copy_node(pnode->children[i],
2847  result_tree.root->children[i]);
2848  result_tree.root->accept_child(i);
2849  }
2850  } else if (child0->node_type == GA_NODE_IND_MOVE_LAST) {
2851  GMM_ASSERT1(child1 == pnode_child, "Cannot extract a factor of a "
2852  "Index_move_last size parameter");
2853  // Copy of the term and other children
2854  result_tree.insert_node(result_tree.root, pnode->node_type);
2855  result_tree.root->pos = pnode->pos;
2856  result_tree.root->expr = pnode->expr;
2857  result_tree.root->children.resize(pnode->children.size(), nullptr);
2858  for (size_type i = 0; i < pnode->children.size(); ++i)
2859  if (pnode->children[i] == pnode_child)
2860  std::swap(result_tree.root->children[i],
2861  result_tree.root->children[0]);
2862  for (size_type i = 0; i < pnode->children.size(); ++i)
2863  if (pnode->children[i] != pnode_child) {
2864  result_tree.copy_node(pnode->children[i],
2865  result_tree.root->children[i]);
2866  result_tree.root->accept_child(i);
2867  }
2868  } else if (child0->node_type == GA_NODE_CONTRACT) {
2869  // Copy of the term and other children
2870  result_tree.insert_node(result_tree.root, pnode->node_type);
2871  result_tree.root->pos = pnode->pos;
2872  result_tree.root->expr = pnode->expr;
2873  result_tree.root->children.resize(pnode->children.size(), nullptr);
2874  for (size_type i = 0; i < pnode->children.size(); ++i)
2875  if (pnode->children[i] == pnode_child)
2876  std::swap(result_tree.root->children[i],
2877  result_tree.root->children[0]);
2878  for (size_type i = 0; i < pnode->children.size(); ++i)
2879  if (pnode->children[i] != pnode_child) {
2880  result_tree.copy_node(pnode->children[i],
2881  result_tree.root->children[i]);
2882  result_tree.root->accept_child(i);
2883  }
2884  } else
2885  GMM_ASSERT1(false, "Cannot extract a factor which is a parameter "
2886  "of a nonlinear operator/function");
2887  break;
2888 
2889  case GA_NODE_C_MATRIX:
2890  result_tree.insert_node(result_tree.root, pnode->node_type);
2891  result_tree.root->pos = pnode->pos;
2892  result_tree.root->expr = pnode->expr;
2893  result_tree.root->children.resize(pnode->children.size());
2894  for (size_type i = 0; i < pnode->children.size(); ++i)
2895  if (pnode_child == pnode->children[i]) {
2896  result_tree.root->children[i] = result_tree.root->children[0];
2897  result_tree.root->children[0] = 0;
2898  }
2899 
2900  for (size_type i = 0; i < pnode->children.size(); ++i) {
2901  if (pnode_child == pnode->children[i]) {
2902  pnode->children[i]
2903  = new ga_tree_node(GA_NODE_ZERO, pnode->pos, pnode->expr);
2904  pnode->children[i]->init_scalar_tensor(scalar_type(0));
2905  pnode->accept_child(i);
2906  }
2907  }
2908  break;
2909 
2910  default: GMM_ASSERT1(false, "Unexpected node type " << pnode->node_type
2911  << " in extract constant term. Internal error.");
2912  }
2913 
2914  pnode_child = pnode;
2915  pnode = pnode->parent;
2916  }
2917 
2918  if (minus_sign) {
2919  result_tree.insert_node(result_tree.root, GA_NODE_OP);
2920  result_tree.root->op_type = GA_UNARY_MINUS;
2921  result_tree.root->pos = pnode->children[0]->pos;
2922  result_tree.root->expr = pnode->children[0]->expr;
2923  }
2924  }
2925 
2926  bool ga_node_extract_constant_term
2927  (ga_tree &tree, pga_tree_node pnode, const ga_workspace &workspace,
2928  const mesh &m) {
2929  bool is_constant = true;
2930  size_type nbch = pnode->children.size();
2931  pga_tree_node child0 = (nbch > 0) ? pnode->children[0] : 0;
2932  // pga_tree_node child1 = (nbch > 1) ? pnode->children[1] : 0;
2933  bool child_0_is_constant = (nbch <= 0) ? true :
2934  ga_node_extract_constant_term(tree, pnode->children[0], workspace, m);
2935  bool child_1_is_constant = (nbch <= 1) ? true :
2936  ga_node_extract_constant_term(tree, pnode->children[1], workspace, m);
2937 
2938  switch (pnode->node_type) {
2939  case GA_NODE_ZERO: is_constant = false; break;
2940 
2941  case GA_NODE_ELEMENTARY_VAL_TEST: case GA_NODE_ELEMENTARY_GRAD_TEST:
2942  case GA_NODE_ELEMENTARY_HESS_TEST: case GA_NODE_ELEMENTARY_DIVERG_TEST:
2943  case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
2944  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
2945  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
2946  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
2947  case GA_NODE_XFEM_PLUS_VAL_TEST: case GA_NODE_XFEM_PLUS_GRAD_TEST:
2948  case GA_NODE_XFEM_PLUS_HESS_TEST: case GA_NODE_XFEM_PLUS_DIVERG_TEST:
2949  case GA_NODE_XFEM_MINUS_VAL_TEST: case GA_NODE_XFEM_MINUS_GRAD_TEST:
2950  case GA_NODE_XFEM_MINUS_HESS_TEST: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
2951  case GA_NODE_VAL_TEST: case GA_NODE_GRAD_TEST: case GA_NODE_PREDEF_FUNC:
2952  case GA_NODE_HESS_TEST: case GA_NODE_DIVERG_TEST: case GA_NODE_RESHAPE:
2953  case GA_NODE_CROSS_PRODUCT:
2954  case GA_NODE_SWAP_IND: case GA_NODE_IND_MOVE_LAST: case GA_NODE_CONTRACT:
2955  case GA_NODE_ELT_SIZE: case GA_NODE_ELT_K: case GA_NODE_ELT_B:
2956  case GA_NODE_CONSTANT: case GA_NODE_X: case GA_NODE_NORMAL:
2957  case GA_NODE_SECONDARY_DOMAIN_X: case GA_NODE_SECONDARY_DOMAIN_NORMAL:
2958  case GA_NODE_OPERATOR:
2959  is_constant = true; break;
2960 
2961  case GA_NODE_ELEMENTARY_VAL: case GA_NODE_ELEMENTARY_GRAD:
2962  case GA_NODE_ELEMENTARY_HESS: case GA_NODE_ELEMENTARY_DIVERG:
2963  case GA_NODE_SECONDARY_DOMAIN_VAL: case GA_NODE_SECONDARY_DOMAIN_GRAD:
2964  case GA_NODE_SECONDARY_DOMAIN_HESS: case GA_NODE_SECONDARY_DOMAIN_DIVERG:
2965  case GA_NODE_XFEM_PLUS_VAL: case GA_NODE_XFEM_PLUS_GRAD:
2966  case GA_NODE_XFEM_PLUS_HESS: case GA_NODE_XFEM_PLUS_DIVERG:
2967  case GA_NODE_XFEM_MINUS_VAL: case GA_NODE_XFEM_MINUS_GRAD:
2968  case GA_NODE_XFEM_MINUS_HESS: case GA_NODE_XFEM_MINUS_DIVERG:
2969  case GA_NODE_VAL: case GA_NODE_GRAD: case GA_NODE_HESS:
2970  case GA_NODE_DIVERG:
2971  is_constant = workspace.is_constant(pnode->name);
2972  break;
2973 
2974  case GA_NODE_INTERPOLATE_VAL:
2975  case GA_NODE_INTERPOLATE_GRAD:
2976  case GA_NODE_INTERPOLATE_HESS:
2977  case GA_NODE_INTERPOLATE_DIVERG:
2978  {
2979  if (!(workspace.is_constant(pnode->name))) {
2980  is_constant = false; break;
2981  }
2982  std::set<var_trans_pair> vars;
2983  workspace.interpolate_transformation(pnode->interpolate_name)
2984  ->extract_variables(workspace, vars, true, m,
2985  pnode->interpolate_name);
2986  for (const var_trans_pair &var : vars) {
2987  if (!(workspace.is_constant(var.varname)))
2988  { is_constant = false; break; }
2989  }
2990  }
2991  break;
2992 
2993  case GA_NODE_INTERPOLATE_FILTER:
2994  if (!child_0_is_constant) {
2995  is_constant = false;
2996  break;
2997  }
2998  [[fallthrough]];
2999  case GA_NODE_INTERPOLATE_VAL_TEST:
3000  case GA_NODE_INTERPOLATE_GRAD_TEST:
3001  case GA_NODE_INTERPOLATE_DIVERG_TEST:
3002  case GA_NODE_INTERPOLATE_HESS_TEST:
3003  case GA_NODE_INTERPOLATE_X:
3004  case GA_NODE_INTERPOLATE_ELT_K: case GA_NODE_INTERPOLATE_ELT_B:
3005  case GA_NODE_INTERPOLATE_NORMAL:
3006  case GA_NODE_INTERPOLATE_DERIVATIVE:
3007  {
3008  std::set<var_trans_pair> vars;
3009  workspace.interpolate_transformation(pnode->interpolate_name)
3010  ->extract_variables(workspace, vars, true, m,
3011  pnode->interpolate_name);
3012  for (const var_trans_pair &var : vars) {
3013  if (!(workspace.is_constant(var.varname)))
3014  { is_constant = false; break; }
3015  }
3016  }
3017  break;
3018 
3019  case GA_NODE_OP:
3020  switch(pnode->op_type) {
3021  case GA_PLUS: case GA_MINUS:
3022  if (!child_0_is_constant && !child_1_is_constant)
3023  { is_constant = false; break; }
3024  break;
3025 
3026  case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM: case GA_SKEW:
3027  case GA_TRACE: case GA_DEVIATOR: case GA_PRINT:
3028  is_constant = child_0_is_constant;
3029  break;
3030 
3031  case GA_DOT: case GA_MULT: case GA_COLON: case GA_TMULT:
3032  case GA_DOTMULT: case GA_DIV: case GA_DOTDIV:
3033  is_constant = (child_0_is_constant && child_1_is_constant);
3034  break;
3035 
3036  default: GMM_ASSERT1(false, "Unexpected operation. Internal error.");
3037  }
3038  break;
3039 
3040  case GA_NODE_C_MATRIX:
3041  for (size_type i = 0; i < pnode->children.size(); ++i) {
3042  if (!(ga_node_extract_constant_term(tree, pnode->children[i],
3043  workspace, m)))
3044  { is_constant = false; break; }
3045  }
3046  break;
3047 
3048  case GA_NODE_PARAMS:
3049  if (child0->node_type == GA_NODE_RESHAPE ||
3050  child0->node_type == GA_NODE_SWAP_IND ||
3051  child0->node_type == GA_NODE_IND_MOVE_LAST) {
3052  is_constant = child_1_is_constant;
3053  } else if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
3054  bool child_2_is_constant
3055  = ga_node_extract_constant_term(tree,pnode->children[2],workspace,m);
3056  is_constant = (child_1_is_constant && child_2_is_constant);
3057  } else if (child0->node_type == GA_NODE_CONTRACT) {
3058  for (size_type i = 1; i < pnode->children.size(); ++i) {
3059  if (!(ga_node_extract_constant_term(tree, pnode->children[i],
3060  workspace, m)))
3061  { is_constant = false; break; }
3062  }
3063  } else if (child0->node_type == GA_NODE_PREDEF_FUNC) {
3064  for (size_type i = 1; i < pnode->children.size(); ++i) {
3065  if (!(ga_node_extract_constant_term(tree, pnode->children[i],
3066  workspace, m)))
3067  { is_constant = false; break; }
3068  }
3069  } else if (child0->node_type == GA_NODE_SPEC_FUNC) {
3070  GMM_ASSERT1(false, "internal error");
3071  } else if (child0->node_type == GA_NODE_OPERATOR) {
3072  for (size_type i = 1; i < pnode->children.size(); ++i) {
3073  if (!(ga_node_extract_constant_term(tree, pnode->children[i],
3074  workspace, m)))
3075  { is_constant = false; break; }
3076  }
3077  } else {
3078  is_constant = child_0_is_constant;
3079  }
3080  break;
3081 
3082  default: GMM_ASSERT1(false, "Unexpected node type " << pnode->node_type
3083  << " in extract constant term. Internal error.");
3084  }
3085 
3086  if (!is_constant) {
3087  pnode->node_type = GA_NODE_ZERO;
3088  tree.clear_children(pnode);
3089  }
3090  return is_constant;
3091  }
3092 
3093 
3094  //=========================================================================
3095  // Extract Neumann terms
3096  //=========================================================================
3097 
3098  static std::string ga_extract_one_Neumann_term
3099  (const std::string &varname,
3100  ga_workspace &workspace, pga_tree_node pnode) {
3101  size_type N = workspace.qdim(varname);
3102  const mesh_fem *mf = workspace.associated_mf(varname);
3103  GMM_ASSERT1(mf, "Works only with fem variables.");
3104  size_type meshdim = mf->linked_mesh().dim();
3105  ga_tree factor;
3106  pga_tree_node new_pnode = nullptr;
3107  ga_extract_factor(factor, pnode, new_pnode);
3108  std::string result;
3109  pga_tree_node nnew_pnode = new_pnode;
3110 
3111  int cas = new_pnode->node_type == GA_NODE_GRAD_TEST ? 0 : 1;
3112  // Allow to detect Trace(Grad_Test_u)
3113  if (cas == 0 && new_pnode->parent &&
3114  new_pnode->parent->node_type == GA_NODE_OP &&
3115  new_pnode->parent->op_type == GA_TRACE) {
3116  cas = 2;
3117  nnew_pnode = new_pnode->parent;
3118  }
3119  bool ok = true;
3120  pga_tree_node colon_pnode = nullptr;
3121  bool quote_before_colon = false;
3122 
3123  // A:Grad_Test_u --> A*Normal if A is a matrix
3124  // Grad_Test_u:A --> A*Normal if A is a matrix
3125  // A*Div_Test_u --> A*Normal if A is a scalar
3126  // Div_Test_u*A --> Normal*A if A is a scalar
3127  // A*(Grad_Test_u)' --> (A)'*Normal if A is a matrix
3128  // intercaled scalar multiplications and divisions are taken into account
3129  while (nnew_pnode->parent) {
3130  pga_tree_node previous_node = nnew_pnode;
3131  nnew_pnode = nnew_pnode->parent;
3132 
3133  if (nnew_pnode->node_type == GA_NODE_OP &&
3134  nnew_pnode->op_type == GA_MULT &&
3135  nnew_pnode->children[0] == previous_node &&
3136  nnew_pnode->children[1]->tensor_proper_size() == 1) {
3137  } else if (nnew_pnode->node_type == GA_NODE_OP &&
3138  nnew_pnode->op_type == GA_MULT &&
3139  nnew_pnode->children[1] == previous_node &&
3140  nnew_pnode->children[0]->tensor_proper_size() == 1) {
3141 
3142  } else if (nnew_pnode->node_type == GA_NODE_OP &&
3143  nnew_pnode->op_type == GA_DIV &&
3144  nnew_pnode->children[0] == previous_node &&
3145  nnew_pnode->children[1]->tensor_proper_size() == 1) {
3146 
3147  } else if (nnew_pnode->node_type == GA_NODE_OP &&
3148  nnew_pnode->op_type == GA_COLON &&
3149  nnew_pnode->children[0] == previous_node &&
3150  nnew_pnode->children[1]->tensor_order() == 2 &&
3151  colon_pnode == 0 && cas == 0) {
3152  std::swap(nnew_pnode->children[0], nnew_pnode->children[1]);
3153  colon_pnode = nnew_pnode;
3154  } else if (nnew_pnode->node_type == GA_NODE_OP &&
3155  nnew_pnode->op_type == GA_COLON &&
3156  nnew_pnode->children[1] == previous_node &&
3157  nnew_pnode->children[0]->tensor_order() == 2 &&
3158  colon_pnode == 0 && cas == 0) {
3159  colon_pnode = nnew_pnode;
3160  } else if (nnew_pnode->node_type == GA_NODE_OP &&
3161  nnew_pnode->op_type == GA_QUOTE &&
3162  colon_pnode == 0 && cas == 0 && !quote_before_colon) {
3163  quote_before_colon = true;
3164  } else ok = false;
3165  }
3166 
3167  if (ok && cas == 0 && !colon_pnode) ok = false;
3168 
3169  if (N == 1) {
3170  new_pnode->node_type = GA_NODE_NORMAL;
3171  result = "(" + ga_tree_to_string(factor) + ")";
3172  } else if (ok) {
3173  switch (cas) {
3174  case 0:
3175  new_pnode->node_type = GA_NODE_NORMAL;
3176  colon_pnode->op_type = GA_MULT;
3177  if (quote_before_colon) {
3178  factor.insert_node(colon_pnode->children[0], GA_NODE_OP);
3179  colon_pnode->children[0]->op_type = GA_QUOTE;
3180  nnew_pnode = new_pnode->parent;
3181  while(nnew_pnode->node_type != GA_NODE_OP ||
3182  nnew_pnode->op_type != GA_QUOTE)
3183  nnew_pnode = nnew_pnode->parent;
3184  factor.replace_node_by_child(nnew_pnode,0);
3185  }
3186  break;
3187  case 1:
3188  new_pnode->node_type = GA_NODE_NORMAL;
3189  break;
3190  case 2:
3191  new_pnode->parent->node_type = GA_NODE_NORMAL;
3192  factor.clear_children(new_pnode->parent);
3193  break;
3194  }
3195  result = "(" + ga_tree_to_string(factor) + ")";
3196 
3197  } else {
3198  // General case
3199 
3200  result = "[";
3201  bgeot::multi_index mi(2); mi[0] = N; mi[1] = meshdim;
3202 
3203  for (size_type i = 0; i < N; ++i) {
3204  factor.clear_children(new_pnode);
3205  new_pnode->node_type = GA_NODE_C_MATRIX;
3206  new_pnode->t.adjust_sizes(mi);
3207  new_pnode->children.resize(N*meshdim);
3208  for (size_type j = 0; j < N; ++j) {
3209  for (size_type k = 0; k < meshdim; ++k) {
3210  if (j == i) {
3211  pga_tree_node param_node = new_pnode->children[k*N+j]
3212  = new ga_tree_node(GA_NODE_PARAMS, pnode->pos, pnode->expr);
3213  new_pnode->accept_child(k+j*meshdim);
3214  param_node->children.resize(2);
3215  param_node->children[0]
3216  = new ga_tree_node(GA_NODE_NORMAL, pnode->pos, pnode->expr);
3217  param_node->accept_child(0);
3218  param_node->children[1]
3219  = new ga_tree_node(GA_NODE_CONSTANT, pnode->pos, pnode->expr);
3220  param_node->accept_child(1);
3221  param_node->children[1]->init_scalar_tensor(scalar_type(k));
3222 
3223  } else {
3224  new_pnode->children[k+j*meshdim]
3225  = new ga_tree_node(GA_NODE_ZERO, pnode->pos, pnode->expr);
3226  new_pnode->children[k+j*meshdim]
3227  ->init_scalar_tensor(scalar_type(0));
3228  new_pnode->accept_child(k+j*meshdim);
3229  }
3230  }
3231  }
3232  result += "(" + ga_tree_to_string(factor) + ")";
3233  if (i < N-1) result += ";";
3234  }
3235  result += "]";
3236  GMM_TRACE2("Warning, generic Neumann term used: " << result);
3237  }
3238 
3239  return result;
3240  }
3241 
3242 
3243  void ga_extract_Neumann_term
3244  (ga_tree &tree, const std::string &varname,
3245  ga_workspace &workspace, pga_tree_node pnode, std::string &result) {
3246 
3247  for (size_type i = 0; i < pnode->children.size(); ++i)
3248  ga_extract_Neumann_term(tree, varname, workspace,
3249  pnode->children[i], result);
3250 
3251  switch (pnode->node_type) {
3252  case GA_NODE_DIVERG_TEST: case GA_NODE_GRAD_TEST:
3253  case GA_NODE_ELEMENTARY_GRAD_TEST: case GA_NODE_ELEMENTARY_DIVERG_TEST:
3254  if (pnode->name.compare(varname) == 0) {
3255  if (result.size()) result += " + ";
3256  result += ga_extract_one_Neumann_term(varname, workspace, pnode);
3257  }
3258  break;
3259  case GA_NODE_INTERPOLATE_GRAD_TEST: case GA_NODE_INTERPOLATE_DIVERG_TEST:
3260  if (pnode->name.compare(varname) == 0)
3261  GMM_ASSERT1(false, "Do not know how to extract a "
3262  "Neumann term with an interpolate transformation");
3263  break;
3264  case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
3265  case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
3266  if (pnode->name.compare(varname) == 0)
3267  GMM_ASSERT1(false, "Do not know how to extract a "
3268  "Neumann term with an two-domain term");
3269  break;
3270  default:
3271  break;
3272  }
3273  }
3274 
3275  //=========================================================================
3276  // Derivation algorithm: derivation of a tree with respect to a variable
3277  // The result tree is not ready to use. It has to be passed again in
3278  // ga_semantic_analysis for enrichment.
3279  //=========================================================================
3280 
3281  static void ga_node_derivation(ga_tree &tree, const ga_workspace &workspace,
3282  const mesh &m,
3283  pga_tree_node pnode,
3284  const std::string &varname,
3285  const std::string &interpolatename,
3286  size_type order, bool any_trans) {
3287 
3288  size_type nbch = pnode->children.size();
3289  pga_tree_node child0 = (nbch > 0) ? pnode->children[0] : 0;
3290  pga_tree_node child1 = (nbch > 1) ? pnode->children[1] : 0;
3291  bool mark0 = ((nbch > 0) ? child0->marked : false);
3292  bool mark1 = ((nbch > 1) ? child1->marked : false);
3293  bgeot::multi_index mi;
3294 
3295  const ga_predef_function_tab &PREDEF_FUNCTIONS
3296  = dal::singleton<ga_predef_function_tab>::instance(0);
3297 
3298  switch (pnode->node_type) {
3299  case GA_NODE_VAL: case GA_NODE_GRAD:
3300  case GA_NODE_HESS: case GA_NODE_DIVERG:
3301  mi.resize(1); mi[0] = 2;
3302  for (size_type i = 0; i < pnode->tensor_order(); ++i)
3303  mi.push_back(pnode->tensor_proper_size(i));
3304  pnode->t.adjust_sizes(mi);
3305  if (pnode->node_type == GA_NODE_VAL)
3306  pnode->node_type = GA_NODE_VAL_TEST;
3307  else if (pnode->node_type == GA_NODE_GRAD)
3308  pnode->node_type = GA_NODE_GRAD_TEST;
3309  else if (pnode->node_type == GA_NODE_HESS)
3310  pnode->node_type = GA_NODE_HESS_TEST;
3311  else if (pnode->node_type == GA_NODE_DIVERG)
3312  pnode->node_type = GA_NODE_DIVERG_TEST;
3313  pnode->test_function_type = order;
3314  break;
3315 
3316  case GA_NODE_INTERPOLATE_VAL:
3317  case GA_NODE_INTERPOLATE_GRAD:
3318  case GA_NODE_INTERPOLATE_HESS:
3319  case GA_NODE_INTERPOLATE_DIVERG:
3320  {
3321  bool is_val(pnode->node_type == GA_NODE_INTERPOLATE_VAL);
3322  bool is_grad(pnode->node_type == GA_NODE_INTERPOLATE_GRAD);
3323  bool is_hess(pnode->node_type == GA_NODE_INTERPOLATE_HESS);
3324  bool is_diverg(pnode->node_type == GA_NODE_INTERPOLATE_DIVERG);
3325 
3326  bool ivar = (pnode->name.compare(varname) == 0 &&
3327  (any_trans ||
3328  pnode->interpolate_name.compare(interpolatename) == 0));
3329  bool itrans = !ivar;
3330  if (!itrans) {
3331  std::set<var_trans_pair> vars;
3332  workspace.interpolate_transformation(pnode->interpolate_name)
3333  ->extract_variables(workspace, vars, true, m,
3334  pnode->interpolate_name);
3335  for (const var_trans_pair &var : vars) {
3336  if (var.varname.compare(varname) == 0 &&
3337  var.transname.compare(interpolatename) == 0)
3338  itrans = true;
3339  }
3340  }
3341 
3342  pga_tree_node pnode_trans = pnode;
3343  if (is_hess) {
3344  GMM_ASSERT1(!itrans, "Sorry, cannot derive a hessian once more");
3345  } else if (itrans && ivar) {
3346  tree.duplicate_with_addition(pnode);
3347  pnode_trans = pnode->parent->children[1];
3348  }
3349  if (ivar) { // Derivative wrt the interpolated variable
3350  mi.resize(1); mi[0] = 2;
3351  for (size_type i = 0; i < pnode->tensor_order(); ++i)
3352  mi.push_back(pnode->tensor_proper_size(i));
3353  pnode->t.adjust_sizes(mi);
3354  if (is_val) // --> t(Qmult*ndof,Qmult*target_dim)
3355  pnode->node_type = GA_NODE_INTERPOLATE_VAL_TEST;
3356  else if (is_grad) // --> t(Qmult*ndof,Qmult*target_dim,N)
3357  pnode->node_type = GA_NODE_INTERPOLATE_GRAD_TEST;
3358  else if (is_hess) // --> t(Qmult*ndof,Qmult*target_dim,N,N)
3359  pnode->node_type = GA_NODE_INTERPOLATE_HESS_TEST;
3360  else if (is_diverg) // --> t(Qmult*ndof)
3361  pnode->node_type = GA_NODE_INTERPOLATE_DIVERG_TEST;
3362  pnode->test_function_type = order;
3363  }
3364 
3365  if (itrans) { // Derivative with respect to a variable that the
3366  // interpolate transformation depends on
3367  const mesh_fem *mf = workspace.associated_mf(pnode_trans->name);
3368  size_type q = workspace.qdim(pnode_trans->name);
3369  size_type n = mf->linked_mesh().dim();
3370  size_type qv = workspace.qdim(varname);
3371  bgeot::multi_index mii = workspace.qdims(pnode_trans->name);
3372 
3373  if (is_val) // --> t(target_dim*Qmult,N)
3374  pnode_trans->node_type = GA_NODE_INTERPOLATE_GRAD;
3375  else if (is_grad || is_diverg) // --> t(target_dim*Qmult,N,N)
3376  pnode_trans->node_type = GA_NODE_INTERPOLATE_HESS;
3377 
3378  if (n > 1) {
3379  if (q == 1 && mii.size() <= 1) { mii.resize(1); mii[0] = n; }
3380  else mii.push_back(qv);
3381  if (is_grad || is_diverg) mii.push_back(n);
3382  }
3383  pnode_trans->t.adjust_sizes(mii);
3384  tree.duplicate_with_operation(pnode_trans,
3385  (n > 1) ? GA_DOT : GA_MULT);
3386  pga_tree_node pnode_der = pnode_trans->parent->children[1];
3387  pnode_der->node_type = GA_NODE_INTERPOLATE_DERIVATIVE;
3388  if (n == 1)
3389  pnode_der->init_vector_tensor(2);
3390  else
3391  pnode_der->init_matrix_tensor(2, n);
3392  pnode_der->test_function_type = order;
3393  pnode_der->name = varname;
3394  pnode_der->interpolate_name_der = pnode_der->interpolate_name;
3395  pnode_der->interpolate_name = interpolatename;
3396 
3397  if (is_diverg) { // --> t(Qmult*ndof)
3398  tree.insert_node(pnode_trans->parent, GA_NODE_OP);
3399  pga_tree_node pnode_tr = pnode_trans->parent->parent;
3400  pnode_tr->op_type = GA_TRACE;
3401  pnode_tr->init_vector_tensor(2);
3402 // pnode_tr->test_function_type = order;
3403 // pnode_tr->name_test1 = pnode_trans->name_test1;
3404 // pnode_tr->name_test2 = pnode_trans->name_test2;
3405  }
3406  }
3407  }
3408  break;
3409 
3410  case GA_NODE_INTERPOLATE_VAL_TEST:
3411  case GA_NODE_INTERPOLATE_GRAD_TEST:
3412  case GA_NODE_INTERPOLATE_DIVERG_TEST:
3413  {
3414  bool is_val(pnode->node_type == GA_NODE_INTERPOLATE_VAL_TEST);
3415  bool is_grad(pnode->node_type == GA_NODE_INTERPOLATE_GRAD_TEST);
3416  bool is_diverg(pnode->node_type == GA_NODE_INTERPOLATE_DIVERG_TEST);
3417 
3418  pga_tree_node pnode_trans = pnode;
3419  const mesh_fem *mf = workspace.associated_mf(pnode_trans->name);
3420  size_type q = workspace.qdim(pnode_trans->name);
3421  size_type n = mf->linked_mesh().dim();
3422  size_type qv = workspace.qdim(varname);
3423  bgeot::multi_index mii = workspace.qdims(pnode_trans->name);
3424  if (is_val) // --> t(Qmult*ndof,Qmult*target_dim,N)
3425  pnode_trans->node_type = GA_NODE_INTERPOLATE_GRAD_TEST;
3426  else if (is_grad || is_diverg) // --> t(Qmult*ndof,Qmult*target_dim,N,N)
3427  pnode_trans->node_type = GA_NODE_INTERPOLATE_HESS_TEST;
3428 
3429  if (q == 1 && mii.size() <= 1) { mii.resize(1); mii[0] = 2; }
3430  else mii.insert(mii.begin(), 2);
3431 
3432  if (n > 1) {
3433  mii.push_back(qv);
3434  if (is_grad || is_diverg) mii.push_back(n);
3435  }
3436  pnode_trans->t.adjust_sizes(mii);
3437  // pnode_trans->test_function_type = order;
3438  tree.duplicate_with_operation(pnode_trans,
3439  (n > 1 ? GA_DOT : GA_MULT));
3440  pga_tree_node pnode_der = pnode_trans->parent->children[1];
3441  pnode_der->node_type = GA_NODE_INTERPOLATE_DERIVATIVE;
3442  if (n == 1)
3443  pnode_der->init_vector_tensor(2);
3444  else
3445  pnode_der->init_matrix_tensor(2, n);
3446  pnode_der->test_function_type = order;
3447  pnode_der->name = varname;
3448  pnode_der->interpolate_name_der = pnode_der->interpolate_name;
3449  pnode_der->interpolate_name = interpolatename;
3450 
3451  if (is_diverg) { // --> t(Qmult*ndof)
3452  tree.insert_node(pnode_trans->parent, GA_NODE_OP);
3453  pga_tree_node pnode_tr = pnode_trans->parent->parent;
3454  pnode_tr->op_type = GA_TRACE;
3455  pnode_tr->init_vector_tensor(2);
3456 // pnode_tr->test_function_type = order;
3457 // pnode_tr->name_test1 = pnode_trans->name_test1;
3458 // pnode_tr->name_test2 = pnode_trans->name_test2;
3459  }
3460  }
3461  break;
3462 
3463  case GA_NODE_INTERPOLATE_HESS_TEST:
3464  GMM_ASSERT1(false, "Sorry, cannot derive a hessian once more");
3465  break;
3466 
3467  case GA_NODE_INTERPOLATE_X:
3468  {
3469  size_type n = m.dim();
3470  pga_tree_node pnode_der = pnode;
3471  pnode_der->node_type = GA_NODE_INTERPOLATE_DERIVATIVE;
3472  if (n == 1)
3473  pnode_der->init_vector_tensor(2);
3474  else
3475  pnode_der->init_matrix_tensor(2, n);
3476  pnode_der->test_function_type = order;
3477  pnode_der->name = varname;
3478  pnode_der->interpolate_name_der = pnode_der->interpolate_name;
3479  pnode_der->interpolate_name = interpolatename;
3480  }
3481  break;
3482 
3483  case GA_NODE_INTERPOLATE_ELT_K:
3484  GMM_ASSERT1(false, "Sorry, cannot derive the interpolated element_K");
3485  break;
3486 
3487  case GA_NODE_INTERPOLATE_ELT_B:
3488  GMM_ASSERT1(false, "Sorry, cannot derive the interpolated element_B");
3489  break;
3490 
3491  case GA_NODE_INTERPOLATE_NORMAL:
3492  GMM_ASSERT1(false, "Sorry, cannot derive the interpolated Normal");
3493  break;
3494 
3495  case GA_NODE_INTERPOLATE_DERIVATIVE:
3496  GMM_ASSERT1(false, "Sorry, second order transformation derivative "
3497  "not taken into account");
3498  break;
3499 
3500  case GA_NODE_INTERPOLATE_FILTER:
3501  ga_node_derivation(tree, workspace, m, child0, varname,
3502  interpolatename, order, any_trans);
3503  break;
3504 
3505  case GA_NODE_SECONDARY_DOMAIN_VAL:
3506  case GA_NODE_SECONDARY_DOMAIN_GRAD:
3507  case GA_NODE_SECONDARY_DOMAIN_HESS:
3508  case GA_NODE_SECONDARY_DOMAIN_DIVERG:
3509  case GA_NODE_ELEMENTARY_VAL:
3510  case GA_NODE_ELEMENTARY_GRAD:
3511  case GA_NODE_ELEMENTARY_HESS:
3512  case GA_NODE_ELEMENTARY_DIVERG:
3513  case GA_NODE_XFEM_PLUS_VAL:
3514  case GA_NODE_XFEM_PLUS_GRAD:
3515  case GA_NODE_XFEM_PLUS_HESS:
3516  case GA_NODE_XFEM_PLUS_DIVERG:
3517  case GA_NODE_XFEM_MINUS_VAL:
3518  case GA_NODE_XFEM_MINUS_GRAD:
3519  case GA_NODE_XFEM_MINUS_HESS:
3520  case GA_NODE_XFEM_MINUS_DIVERG:
3521  mi.resize(1); mi[0] = 2;
3522  for (size_type i = 0; i < pnode->tensor_order(); ++i)
3523  mi.push_back(pnode->tensor_proper_size(i));
3524  pnode->t.adjust_sizes(mi);
3525  switch(pnode->node_type) {
3526  case GA_NODE_SECONDARY_DOMAIN_VAL:
3527  pnode->node_type = GA_NODE_SECONDARY_DOMAIN_VAL_TEST; break;
3528  case GA_NODE_SECONDARY_DOMAIN_GRAD:
3529  pnode->node_type = GA_NODE_SECONDARY_DOMAIN_GRAD_TEST; break;
3530  case GA_NODE_SECONDARY_DOMAIN_HESS:
3531  pnode->node_type = GA_NODE_SECONDARY_DOMAIN_HESS_TEST; break;
3532  case GA_NODE_SECONDARY_DOMAIN_DIVERG:
3533  pnode->node_type = GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST; break;
3534  case GA_NODE_ELEMENTARY_VAL:
3535  pnode->node_type = GA_NODE_ELEMENTARY_VAL_TEST; break;
3536  case GA_NODE_ELEMENTARY_GRAD:
3537  pnode->node_type = GA_NODE_ELEMENTARY_GRAD_TEST; break;
3538  case GA_NODE_ELEMENTARY_HESS:
3539  pnode->node_type = GA_NODE_ELEMENTARY_HESS_TEST; break;
3540  case GA_NODE_ELEMENTARY_DIVERG:
3541  pnode->node_type = GA_NODE_ELEMENTARY_DIVERG_TEST; break;
3542  case GA_NODE_XFEM_PLUS_VAL:
3543  pnode->node_type = GA_NODE_XFEM_PLUS_VAL_TEST; break;
3544  case GA_NODE_XFEM_PLUS_GRAD:
3545  pnode->node_type = GA_NODE_XFEM_PLUS_GRAD_TEST; break;
3546  case GA_NODE_XFEM_PLUS_HESS:
3547  pnode->node_type = GA_NODE_XFEM_PLUS_HESS_TEST; break;
3548  case GA_NODE_XFEM_PLUS_DIVERG:
3549  pnode->node_type = GA_NODE_XFEM_PLUS_DIVERG_TEST; break;
3550  case GA_NODE_XFEM_MINUS_VAL:
3551  pnode->node_type = GA_NODE_XFEM_MINUS_VAL_TEST; break;
3552  case GA_NODE_XFEM_MINUS_GRAD:
3553  pnode->node_type = GA_NODE_XFEM_MINUS_GRAD_TEST; break;
3554  case GA_NODE_XFEM_MINUS_HESS:
3555  pnode->node_type = GA_NODE_XFEM_MINUS_HESS_TEST; break;
3556  case GA_NODE_XFEM_MINUS_DIVERG:
3557  pnode->node_type = GA_NODE_XFEM_MINUS_DIVERG_TEST; break;
3558  default : GMM_ASSERT1(false, "internal error");
3559  }
3560  pnode->test_function_type = order;
3561  break;
3562 
3563  case GA_NODE_OP:
3564  switch(pnode->op_type) {
3565  case GA_PLUS: case GA_MINUS:
3566  if (mark0 && mark1) {
3567  ga_node_derivation(tree, workspace, m, child0, varname,
3568  interpolatename, order, any_trans);
3569  ga_node_derivation(tree, workspace, m, child1, varname,
3570  interpolatename, order, any_trans);
3571  } else if (mark0) {
3572  ga_node_derivation(tree, workspace, m, child0, varname,
3573  interpolatename, order, any_trans);
3574  tree.replace_node_by_child(pnode, 0);
3575  } else {
3576  ga_node_derivation(tree, workspace, m, child1, varname,
3577  interpolatename, order, any_trans);
3578  if (pnode->op_type == GA_MINUS) {
3579  pnode->op_type = GA_UNARY_MINUS;
3580  tree.clear_node(child0);
3581  }
3582  else
3583  tree.replace_node_by_child(pnode, 1);
3584  }
3585  break;
3586 
3587  case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM: case GA_SKEW:
3588  case GA_TRACE: case GA_DEVIATOR: case GA_PRINT:
3589  ga_node_derivation(tree, workspace, m, child0, varname,
3590  interpolatename, order, any_trans);
3591  break;
3592 
3593  case GA_DOT: case GA_MULT: case GA_COLON: case GA_TMULT:
3594  case GA_DOTMULT:
3595  if (mark0 && mark1) {
3596  if (sub_tree_are_equal(child0, child1, workspace, 0) &&
3597  (pnode->op_type != GA_MULT || child0->tensor_order() < 2) &&
3598  (pnode->op_type != GA_DOT || child0->tensor_order() < 2)) {
3599  ga_node_derivation(tree, workspace, m, child1, varname,
3600  interpolatename, order, any_trans);
3601  tree.insert_node(pnode, GA_NODE_OP);
3602  pnode->parent->op_type = GA_MULT;
3603  tree.add_child(pnode->parent);
3604  pnode->parent->children[1]->node_type = GA_NODE_CONSTANT;
3605  pnode->parent->children[1]->init_scalar_tensor(scalar_type(2));
3606  } else {
3607  tree.duplicate_with_addition(pnode);
3608  if ((pnode->op_type == GA_COLON && child0->tensor_order() == 2) ||
3609  (pnode->op_type == GA_DOT && child0->tensor_order() == 1 &&
3610  child1->tensor_order() == 1) ||
3611  pnode->op_type == GA_DOTMULT ||
3612  (child0->tensor_proper_size()== 1 &&
3613  child1->tensor_proper_size()== 1))
3614  std::swap(pnode->children[0], pnode->children[1]);
3615  ga_node_derivation(tree, workspace, m, child0, varname,
3616  interpolatename, order, any_trans);
3617  ga_node_derivation(tree, workspace, m,
3618  pnode->parent->children[1]->children[1],
3619  varname, interpolatename, order, any_trans);
3620  }
3621  } else if (mark0) {
3622  ga_node_derivation(tree, workspace, m, child0, varname,
3623  interpolatename, order, any_trans);
3624  } else
3625  ga_node_derivation(tree, workspace, m, child1, varname,
3626  interpolatename, order, any_trans);
3627  break;
3628 
3629  case GA_DIV: case GA_DOTDIV:
3630  if (mark1) {
3631  if (pnode->children[0]->node_type == GA_NODE_CONSTANT)
3632  gmm::scale(pnode->children[0]->tensor().as_vector(),
3633  scalar_type(-1));
3634  else {
3635  if (mark0) {
3636  tree.duplicate_with_substraction(pnode);
3637  ga_node_derivation(tree, workspace, m, child0, varname,
3638  interpolatename, order, any_trans);
3639  pnode = pnode->parent->children[1];
3640  } else {
3641  tree.insert_node(pnode, GA_NODE_OP);
3642  pnode->parent->op_type = GA_UNARY_MINUS;
3643  }
3644  }
3645  tree.insert_node(pnode->children[1], GA_NODE_PARAMS);
3646  pga_tree_node pnode_param = pnode->children[1];
3647  tree.add_child(pnode_param);
3648  std::swap(pnode_param->children[0], pnode_param->children[1]);
3649  pnode_param->children[0]->node_type = GA_NODE_PREDEF_FUNC;
3650  pnode_param->children[0]->name = "sqr";
3651  tree.insert_node(pnode, GA_NODE_OP);
3652  pga_tree_node pnode_mult = pnode->parent;
3653  if (pnode->op_type == GA_DOTDIV)
3654  pnode_mult->op_type = GA_DOTMULT;
3655  else
3656  pnode_mult->op_type = GA_MULT;
3657  pnode_mult->children.resize(2, nullptr);
3658  tree.copy_node(pnode_param->children[1], pnode_mult->children[1]);
3659  pnode_mult->accept_child(1);
3660  ga_node_derivation(tree, workspace, m, pnode_mult->children[1],
3661  varname, interpolatename, order, any_trans);
3662  } else {
3663  ga_node_derivation(tree, workspace, m, child0, varname,
3664  interpolatename, order, any_trans);
3665  }
3666  break;
3667 
3668  default: GMM_ASSERT1(false, "Unexpected operation. Internal error.");
3669  }
3670  break;
3671 
3672  case GA_NODE_C_MATRIX:
3673  for (size_type i = 0; i < pnode->children.size(); ++i) {
3674  if (pnode->children[i]->marked)
3675  ga_node_derivation(tree, workspace, m, pnode->children[i],
3676  varname, interpolatename, order, any_trans);
3677  else {
3678  pnode->children[i]->init_scalar_tensor(scalar_type(0));
3679  pnode->children[i]->node_type = GA_NODE_ZERO;
3680  tree.clear_children(pnode->children[i]);
3681  }
3682  }
3683  break;
3684 
3685  case GA_NODE_PARAMS:
3686  if (child0->node_type == GA_NODE_RESHAPE ||
3687  child0->node_type == GA_NODE_SWAP_IND||
3688  child0->node_type == GA_NODE_IND_MOVE_LAST) {
3689  ga_node_derivation(tree, workspace, m, pnode->children[1],
3690  varname, interpolatename, order, any_trans);
3691  } else if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
3692  pga_tree_node child2 = pnode->children[2];
3693  bool mark2 = child2->marked;
3694  if (mark1 && mark2) {
3695  tree.duplicate_with_addition(pnode);
3696  ga_node_derivation(tree, workspace, m, child1, varname,
3697  interpolatename, order, any_trans);
3698  ga_node_derivation(tree, workspace, m,
3699  pnode->parent->children[1]->children[2],
3700  varname, interpolatename, order, any_trans);
3701  } else if (mark1) {
3702  ga_node_derivation(tree, workspace, m, child1, varname,
3703  interpolatename, order, any_trans);
3704  } else
3705  ga_node_derivation(tree, workspace, m, child2, varname,
3706  interpolatename, order, any_trans);
3707  } else if (child0->node_type == GA_NODE_CONTRACT) {
3708 
3709  if (pnode->children.size() == 4) {
3710  ga_node_derivation(tree, workspace, m, pnode->children[1],
3711  varname, interpolatename, order, any_trans);
3712  } else if (pnode->children.size() == 5 || pnode->children.size() == 7) {
3713  size_type n2 = (pnode->children.size()==5) ? 3 : 4;
3714  pga_tree_node child2 = pnode->children[n2];
3715 
3716  if (mark1 && child2->marked) {
3717  tree.duplicate_with_addition(pnode);
3718  ga_node_derivation(tree, workspace, m, child1, varname,
3719  interpolatename, order, any_trans);
3720  ga_node_derivation(tree, workspace, m,
3721  pnode->parent->children[1]->children[n2],
3722  varname, interpolatename, order, any_trans);
3723  } else if (mark1) {
3724  ga_node_derivation(tree, workspace, m, child1, varname,
3725  interpolatename, order, any_trans);
3726  } else
3727  ga_node_derivation(tree, workspace, m, child2, varname,
3728  interpolatename, order, any_trans);
3729 
3730  } else GMM_ASSERT1(false, "internal error");
3731  } else if (child0->node_type == GA_NODE_PREDEF_FUNC) {
3732  std::string name = child0->name;
3733  ga_predef_function_tab::const_iterator it = PREDEF_FUNCTIONS.find(name);
3734  const ga_predef_function &F = it->second;
3735 
3736  if (F.nbargs() == 1) {
3737  switch (F.dtype()) {
3738  case 0:
3739  GMM_ASSERT1(false, "Cannot derive function " << child0->name
3740  << ". No derivative provided or not derivable function.");
3741  break;
3742  case 1:
3743  child0->name = F.derivative1();
3744  break;
3745  case 2: case 3:
3746  {
3747  child0->name = "DER_PDFUNC_" + child0->name;
3748  if (F.dtype() == 2)
3749  ga_define_function(child0->name, 1, F.derivative1());
3750  else {
3751  std::string expr = ga_derivative_scalar_function(F.expr(),"t");
3752  ga_define_function(child0->name, 1, expr);
3753  }
3754  // Inline extension if the derivative is affine (for instance
3755  // for sqr)
3756  ga_predef_function_tab::const_iterator
3757  itp = PREDEF_FUNCTIONS.find(child0->name);
3758  const ga_predef_function &Fp = itp->second;
3759  if (Fp.is_affine("t")) {
3760  scalar_type b = Fp(scalar_type(0));
3761  scalar_type a = Fp(scalar_type(1)) - b;
3762  pnode->node_type = GA_NODE_OP;
3763  pnode->op_type = GA_MULT;
3764  child0->init_scalar_tensor(a);
3765  child0->node_type = ((a == scalar_type(0)) ?
3766  GA_NODE_ZERO : GA_NODE_CONSTANT);
3767  if (b != scalar_type(0)) {
3768  tree.insert_node(pnode, GA_NODE_OP);
3769  pnode->parent->op_type = (b > 0) ? GA_PLUS : GA_MINUS;
3770  tree.add_child(pnode->parent);
3771  pga_tree_node pnode_cte = pnode->parent->children[1];
3772  pnode_cte->node_type = GA_NODE_CONSTANT;
3773  pnode_cte->t = pnode->t;
3774  std::fill(pnode_cte->tensor().begin(),
3775  pnode_cte->tensor().end(), gmm::abs(b));
3776  pnode = pnode->parent;
3777  }
3778  }
3779  }
3780  break;
3781  }
3782  if (pnode->children.size() >= 2) {
3783  tree.insert_node(pnode, GA_NODE_OP);
3784  pga_tree_node pnode_op = pnode->parent;
3785  if (child1->tensor_order() == 0)
3786  pnode_op->op_type = GA_MULT;
3787  else
3788  pnode_op->op_type = GA_DOTMULT;
3789  pnode_op->children.resize(2, nullptr);
3790  tree.copy_node(child1, pnode_op->children[1]);
3791  pnode_op->accept_child(1);
3792  ga_node_derivation(tree, workspace, m, pnode_op->children[1],
3793  varname, interpolatename, order, any_trans);
3794  }
3795  } else {
3796  pga_tree_node child2 = pnode->children[2];
3797  pga_tree_node pg2 = pnode;
3798 
3799  if (child1->marked && child2->marked) {
3800  tree.duplicate_with_addition(pnode);
3801  pg2 = pnode->parent->children[1];
3802  }
3803 
3804  if (child1->marked) {
3805  switch (F.dtype()) {
3806  case 0:
3807  GMM_ASSERT1(false, "Cannot derive function " << child0->name
3808  << ". No derivative provided");
3809  break;
3810  case 1:
3811  child0->name = F.derivative1();
3812  break;
3813  case 2:
3814  child0->name = "DER_PDFUNC1_" + child0->name;
3815  ga_define_function(child0->name, 2, F.derivative1());
3816  break;
3817  case 3:
3818  child0->name = "DER_PDFUNC1_" + child0->name;
3819  std::string expr = ga_derivative_scalar_function(F.expr(), "t");
3820  ga_define_function(child0->name, 2, expr);
3821  break;
3822  }
3823  tree.insert_node(pnode, GA_NODE_OP);
3824  pga_tree_node pnode_op = pnode->parent;
3825  if (child1->tensor_order() == 0)
3826  pnode_op->op_type = GA_MULT;
3827  else
3828  pnode_op->op_type = GA_DOTMULT;
3829  pnode_op->children.resize(2, nullptr);
3830  tree.copy_node(child1, pnode_op->children[1]);
3831  pnode_op->accept_child(1);
3832  ga_node_derivation(tree, workspace, m, pnode_op->children[1],
3833  varname, interpolatename, order, any_trans);
3834  }
3835  if (child2->marked) {
3836  pnode = pg2;
3837  child0 = pnode->children[0]; child1 = pnode->children[1];
3838  child2 = pnode->children[2];
3839 
3840  switch (F.dtype()) {
3841  case 0:
3842  GMM_ASSERT1(false, "Cannot derive function " << child0->name
3843  << ". No derivative provided");
3844  break;
3845  case 1:
3846  child0->name = F.derivative2();
3847  break;
3848  case 2:
3849  child0->name = "DER_PDFUNC2_" + child0->name;
3850  ga_define_function(child0->name, 2, F.derivative2());
3851  break;
3852  case 3:
3853  child0->name = "DER_PDFUNC2_" + child0->name;
3854  std::string expr = ga_derivative_scalar_function(F.expr(), "u");
3855  ga_define_function(child0->name, 2, expr);
3856  break;
3857  }
3858  tree.insert_node(pnode, GA_NODE_OP);
3859  pga_tree_node pnode_op = pnode->parent;
3860  if (child2->tensor_order() == 0)
3861  pnode_op->op_type = GA_MULT;
3862  else
3863  pnode_op->op_type = GA_DOTMULT;
3864  pnode_op->children.resize(2, nullptr);
3865  tree.copy_node(child2, pnode_op->children[1]);
3866  pnode_op->accept_child(1);
3867  ga_node_derivation(tree, workspace, m, pnode_op->children[1],
3868  varname, interpolatename, order, any_trans);
3869  }
3870  }
3871  } else if (child0->node_type == GA_NODE_SPEC_FUNC) {
3872  GMM_ASSERT1(false, "internal error");
3873  } else if (child0->node_type == GA_NODE_OPERATOR) {
3874  if (child0->der2)
3875  GMM_ASSERT1(false, "Error in derivation of the assembly string. "
3876  "Cannot derive again operator " << child0->name);
3877 
3878  size_type nbargs_der = 0;
3879  for (size_type i = 1; i < pnode->children.size(); ++i)
3880  if (pnode->children[i]->marked) ++nbargs_der;
3881  pga_tree_node pnode2 = 0;
3882 
3883  size_type j = 0;
3884  for (size_type i = 1; i < pnode->children.size(); ++i) {
3885  if (pnode->children[i]->marked) {
3886  ++j;
3887  if (j != nbargs_der) {
3888  tree.insert_node(pnode, GA_NODE_OP);
3889  pga_tree_node pnode_op = pnode->parent;
3890  pnode_op->node_type = GA_NODE_OP;
3891  pnode_op->op_type = GA_PLUS;
3892  pnode_op->children.resize(2, nullptr);
3893  tree.copy_node(pnode, pnode_op->children[1]);
3894  pnode_op->accept_child(1);
3895  pnode2 = pnode_op->children[1];
3896  }
3897  else pnode2 = pnode;
3898 
3899  if (child0->der1)
3900  pnode2->children[0]->der2 = i;
3901  else
3902  pnode2->children[0]->der1 = i;
3903  tree.insert_node(pnode2, GA_NODE_OP);
3904  pga_tree_node pnode_op = pnode2->parent;
3905  // Computation of contraction order
3906  size_type red = pnode->children[i]->tensor_order();
3907  switch (red) {
3908  case 0 : pnode_op->op_type = GA_MULT; break;
3909  case 1 : pnode_op->op_type = GA_DOT; break;
3910  case 2 : pnode_op->op_type = GA_COLON; break;
3911  default: GMM_ASSERT1(false, "Error in derivation of the assembly "
3912  "string. Bad reduction order.")
3913  }
3914  pnode_op->children.resize(2, nullptr);
3915  tree.copy_node(pnode->children[i], pnode_op->children[1]);
3916  pnode_op->accept_child(1);
3917  ga_node_derivation(tree, workspace, m, pnode_op->children[1],
3918  varname, interpolatename, order, any_trans);
3919 
3920  if (pnode2->children[0]->name.compare("Norm_sqr") == 0
3921  && pnode2->children[0]->der1 == 1) {
3922  pnode2->node_type = GA_NODE_OP;
3923  pnode2->op_type = GA_MULT;
3924  pnode2->children[0]->node_type = GA_NODE_CONSTANT;
3925  pnode2->children[0]->init_scalar_tensor(scalar_type(2));
3926  }
3927 
3928 
3929  }
3930  }
3931 
3932  } else {
3933  ga_node_derivation(tree, workspace, m, child0, varname,
3934  interpolatename, order, any_trans);
3935  }
3936  break;
3937 
3938  default: GMM_ASSERT1(false, "Unexpected node type " << pnode->node_type
3939  << " in derivation. Internal error.");
3940  }
3941  }
3942 
3943  // The tree is modified. Should be copied first and passed to
3944  // ga_semantic_analysis after for enrichment
3945  void ga_derivative(ga_tree &tree, const ga_workspace &workspace,
3946  const mesh &m, const std::string &varname,
3947  const std::string &interpolatename, size_type order) {
3948  // cout << "Will derive : " << ga_tree_to_string(tree) << endl;
3949  if (!(tree.root)) return;
3950  if (ga_node_mark_tree_for_variable(tree.root, workspace, m, varname,
3951  interpolatename))
3952  ga_node_derivation(tree, workspace, m, tree.root, varname,
3953  interpolatename, order);
3954  else
3955  tree.clear();
3956  // cout << "Derivation done : " << ga_tree_to_string(tree) << endl;
3957  }
3958 
3959  //=========================================================================
3960  // Gradient algorithm: gradient of a tree.
3961  // The result tree is not ready to use. It has to be passed again in
3962  // ga_semantic_analysis for enrichment.
3963  //=========================================================================
3964 
3965  static bool ga_node_mark_tree_for_grad(pga_tree_node pnode,
3966  const ga_workspace &workspace) {
3967  bool marked = false;
3968  for (size_type i = 0; i < pnode->children.size(); ++i)
3969  if (ga_node_mark_tree_for_grad(pnode->children[i], workspace))
3970  marked = true;
3971 
3972  bool plain_node(pnode->node_type == GA_NODE_VAL ||
3973  pnode->node_type == GA_NODE_GRAD ||
3974  pnode->node_type == GA_NODE_HESS ||
3975  pnode->node_type == GA_NODE_DIVERG);
3976  bool test_node(pnode->node_type == GA_NODE_VAL_TEST ||
3977  pnode->node_type == GA_NODE_GRAD_TEST ||
3978  pnode->node_type == GA_NODE_HESS_TEST ||
3979  pnode->node_type == GA_NODE_DIVERG_TEST);
3980  bool interpolate_node(pnode->node_type == GA_NODE_INTERPOLATE_VAL ||
3981  pnode->node_type == GA_NODE_INTERPOLATE_GRAD ||
3982  pnode->node_type == GA_NODE_INTERPOLATE_HESS ||
3983  pnode->node_type == GA_NODE_INTERPOLATE_DIVERG);
3984  bool elementary_node(pnode->node_type == GA_NODE_ELEMENTARY_VAL ||
3985  pnode->node_type == GA_NODE_ELEMENTARY_GRAD ||
3986  pnode->node_type == GA_NODE_ELEMENTARY_HESS ||
3987  pnode->node_type == GA_NODE_ELEMENTARY_DIVERG);
3988  bool xfem_node(pnode->node_type == GA_NODE_XFEM_PLUS_VAL ||
3989  pnode->node_type == GA_NODE_XFEM_PLUS_GRAD ||
3990  pnode->node_type == GA_NODE_XFEM_PLUS_HESS ||
3991  pnode->node_type == GA_NODE_XFEM_PLUS_DIVERG ||
3992  pnode->node_type == GA_NODE_XFEM_MINUS_VAL ||
3993  pnode->node_type == GA_NODE_XFEM_MINUS_GRAD ||
3994  pnode->node_type == GA_NODE_XFEM_MINUS_HESS ||
3995  pnode->node_type == GA_NODE_XFEM_MINUS_DIVERG);
3996  bool interpolate_test_node
3997  (pnode->node_type == GA_NODE_INTERPOLATE_VAL_TEST ||
3998  pnode->node_type == GA_NODE_INTERPOLATE_GRAD_TEST ||
3999  pnode->node_type == GA_NODE_INTERPOLATE_HESS_TEST ||
4000  pnode->node_type == GA_NODE_INTERPOLATE_DIVERG_TEST);
4001 
4002  if ((plain_node || test_node || interpolate_node ||
4003  elementary_node || xfem_node) &&
4004  (workspace.associated_mf(pnode->name) != 0)) marked = true;
4005 
4006  if (pnode->node_type == GA_NODE_X ||
4007  pnode->node_type == GA_NODE_NORMAL) marked = true;
4008 
4009  if ((interpolate_node || interpolate_test_node) &&
4010  (workspace.associated_mf(pnode->name) != 0)) marked = true;
4011 
4012  if (pnode->node_type == GA_NODE_INTERPOLATE_X ||
4013  pnode->node_type == GA_NODE_INTERPOLATE_ELT_K ||
4014  pnode->node_type == GA_NODE_INTERPOLATE_ELT_B ||
4015  pnode->node_type == GA_NODE_INTERPOLATE_NORMAL) marked = true;
4016 
4017  pnode->marked = marked;
4018  return marked;
4019  }
4020 
4021  static void ga_node_grad(ga_tree &tree, const ga_workspace &workspace,
4022  const mesh &m, pga_tree_node pnode) {
4023  // cout << "Gradient of "; ga_print_node(pnode, cout); cout << endl;
4024  size_type meshdim = (&m == &dummy_mesh()) ? 1 : m.dim();
4025  size_type nbch = pnode->children.size();
4026  pga_tree_node child0 = (nbch > 0) ? pnode->children[0] : 0;
4027  pga_tree_node child1 = (nbch > 1) ? pnode->children[1] : 0;
4028  bool mark0 = ((nbch > 0) ? child0->marked : false);
4029  bool mark1 = ((nbch > 1) ? child1->marked : false);
4030  bgeot::multi_index mi;
4031 
4032  const ga_predef_function_tab &PREDEF_FUNCTIONS
4033  = dal::singleton<ga_predef_function_tab>::instance(0);
4034 
4035  switch (pnode->node_type) {
4036  case GA_NODE_VAL: case GA_NODE_VAL_TEST:
4037  if (pnode->node_type == GA_NODE_VAL)
4038  pnode->node_type = GA_NODE_GRAD;
4039  else
4040  pnode->node_type = GA_NODE_GRAD_TEST;
4041  mi = pnode->tensor().sizes();
4042  if (mi.back() != 1) mi.push_back(m.dim()); else mi.back() = m.dim();
4043  pnode->t.adjust_sizes(mi);
4044  break;
4045  case GA_NODE_GRAD: case GA_NODE_GRAD_TEST:
4046  if (pnode->node_type == GA_NODE_GRAD)
4047  pnode->node_type = GA_NODE_HESS;
4048  else
4049  pnode->node_type = GA_NODE_HESS_TEST;
4050  mi = pnode->tensor().sizes();
4051  if (mi.back() != 1) mi.push_back(m.dim()); else mi.back() = m.dim();
4052  pnode->t.adjust_sizes(mi);
4053  break;
4054  case GA_NODE_HESS: case GA_NODE_HESS_TEST:
4055  GMM_ASSERT1(false, "Sorry, cannot derive an Hessian once more");
4056  break;
4057  case GA_NODE_DIVERG: case GA_NODE_DIVERG_TEST: // Hess_u : Id(meshdim)
4058  if (pnode->node_type == GA_NODE_DIVERG)
4059  pnode->node_type = GA_NODE_HESS;
4060  else
4061  pnode->node_type = GA_NODE_HESS_TEST;
4062  mi = pnode->tensor().sizes();
4063  mi.pop_back(), mi.push_back(m.dim());
4064  if (m.dim() > 1) mi.push_back(m.dim());
4065  pnode->t.adjust_sizes(mi);
4066  tree.duplicate_with_operation(pnode, GA_COLON);
4067  child0 = pnode;
4068  pnode = pnode->parent;
4069  child1 = pnode->children[1];
4070  child1->init_matrix_tensor(meshdim, meshdim);
4071  gmm::clear(pnode->tensor().as_vector());
4072  for (size_type i = 0; i < meshdim; ++i)
4073  child1->tensor()(i,i) = scalar_type(1);
4074  child1->node_type = GA_NODE_CONSTANT;
4075  break;
4076 
4077  case GA_NODE_INTERPOLATE_HESS_TEST:
4078  case GA_NODE_INTERPOLATE_HESS:
4079  case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
4080  case GA_NODE_SECONDARY_DOMAIN_HESS:
4081  GMM_ASSERT1(false, "Sorry, cannot derive a hessian once more");
4082  break;
4083 
4084  case GA_NODE_INTERPOLATE_VAL:
4085  case GA_NODE_INTERPOLATE_GRAD:
4086  case GA_NODE_INTERPOLATE_DIVERG:
4087  case GA_NODE_INTERPOLATE_VAL_TEST:
4088  case GA_NODE_INTERPOLATE_GRAD_TEST:
4089  case GA_NODE_INTERPOLATE_DIVERG_TEST:
4090  case GA_NODE_INTERPOLATE_X:
4091  {
4092  std::string &tname = pnode->interpolate_name;
4093  auto ptrans = workspace.interpolate_transformation(tname);
4094  std::string expr_trans = ptrans->expression();
4095  if (expr_trans.size() == 0)
4096  GMM_ASSERT1(false, "Sorry, the gradient of tranformation "
4097  << tname << " cannot be calculated. "
4098  "The gradient computation is available only for "
4099  "transformations having an explicit expression");
4100 
4101  ga_tree trans_tree;
4102  ga_read_string(expr_trans, trans_tree, workspace.macro_dictionary());
4103  ga_semantic_analysis(trans_tree, workspace, m,
4104  ref_elt_dim_of_mesh(m, -1), false, false, 1);
4105  if (trans_tree.root) {
4106  ga_node_grad(trans_tree, workspace, m, trans_tree.root);
4107  ga_semantic_analysis(trans_tree, workspace, m,
4108  ref_elt_dim_of_mesh(m, -1), false, false, 1);
4109 
4110  GMM_ASSERT1(trans_tree.root->tensor().sizes().size() == 2,
4111  "Problem in transformation" << tname);
4112  size_type trans_dim = trans_tree.root->tensor().sizes()[1];
4113 
4114  tree.duplicate_with_operation(pnode, GA_DOT);
4115 
4116  if (pnode->node_type == GA_NODE_INTERPOLATE_VAL) {
4117  pnode->node_type = GA_NODE_INTERPOLATE_GRAD;
4118  mi = pnode->tensor().sizes();
4119  if (mi.back() != 1) mi.push_back(trans_dim);
4120  else mi.back() = trans_dim;
4121  pnode->t.adjust_sizes(mi);
4122  } else if (pnode->node_type == GA_NODE_INTERPOLATE_GRAD) {
4123  pnode->node_type = GA_NODE_INTERPOLATE_HESS;
4124  mi = pnode->tensor().sizes();
4125  if (mi.back() != 1) mi.push_back(trans_dim);
4126  else mi.back() = trans_dim;
4127  pnode->t.adjust_sizes(mi);
4128  } else if (pnode->node_type == GA_NODE_INTERPOLATE_VAL_TEST) {
4129  pnode->node_type = GA_NODE_INTERPOLATE_GRAD_TEST;
4130  mi = pnode->tensor().sizes();
4131  if (mi.back() != 1) mi.push_back(trans_dim);
4132  else mi.back() = trans_dim;
4133  pnode->t.adjust_sizes(mi);
4134  } else if (pnode->node_type == GA_NODE_INTERPOLATE_GRAD_TEST) {
4135  pnode->node_type = GA_NODE_INTERPOLATE_HESS_TEST;
4136  mi = pnode->tensor().sizes();
4137  if (mi.back() != 1) mi.push_back(trans_dim);
4138  else mi.back() = trans_dim;
4139  pnode->t.adjust_sizes(mi);
4140  } else if (pnode->node_type == GA_NODE_INTERPOLATE_DIVERG ||
4141  pnode->node_type == GA_NODE_INTERPOLATE_DIVERG_TEST) {
4142  if (pnode->node_type == GA_NODE_INTERPOLATE_DIVERG)
4143  pnode->node_type = GA_NODE_INTERPOLATE_HESS;
4144  else
4145  pnode->node_type = GA_NODE_INTERPOLATE_HESS_TEST;
4146  mi = pnode->tensor().sizes();
4147  mi.pop_back(), mi.push_back(trans_dim);
4148  if (trans_dim > 1) mi.push_back(trans_dim);
4149  pnode->t.adjust_sizes(mi);
4150  tree.duplicate_with_operation(pnode, GA_COLON);
4151  child0 = pnode;
4152  pnode = pnode->parent;
4153  child1 = pnode->children[1];
4154  child1->init_matrix_tensor(trans_dim, trans_dim);
4155  gmm::clear(pnode->tensor().as_vector());
4156  for (size_type i = 0; i < trans_dim; ++i)
4157  child1->tensor()(i,i) = scalar_type(1);
4158  child1->node_type = GA_NODE_CONSTANT;
4159  } else if (pnode->node_type == GA_NODE_INTERPOLATE_X) {
4160  pnode->node_type = GA_NODE_CONSTANT;
4161  if (pnode->nbc1) {
4162  pnode->init_vector_tensor(trans_dim);
4163  gmm::clear(pnode->tensor().as_vector());
4164  pnode->tensor()[pnode->nbc1-1] = scalar_type(1);
4165  } else {
4166  pnode->init_matrix_tensor(trans_dim, trans_dim);
4167  gmm::clear(pnode->tensor().as_vector());
4168  for (size_type i = 0; i < trans_dim; ++i)
4169  pnode->tensor()(i,i) = scalar_type(1);
4170  }
4171  }
4172 
4173  tree.clear_node_rec(pnode->parent->children[1]);
4174  pnode->parent->children[1] = nullptr;
4175  tree.copy_node(trans_tree.root, pnode->parent->children[1]);
4176  pnode->parent->accept_child(1);
4177  } else {
4178  pnode->node_type = GA_NODE_ZERO;
4179  mi = pnode->tensor().sizes();
4180  mi.push_back(m.dim());
4181  gmm::clear(pnode->tensor().as_vector());
4182  }
4183  }
4184  break;
4185 
4186  case GA_NODE_X:
4187  pnode->node_type = GA_NODE_CONSTANT;
4188  if (pnode->nbc1) {
4189  pnode->init_vector_tensor(meshdim);
4190  gmm::clear(pnode->tensor().as_vector());
4191  pnode->tensor()[pnode->nbc1-1] = scalar_type(1);
4192  } else {
4193  pnode->init_matrix_tensor(meshdim, meshdim);
4194  gmm::clear(pnode->tensor().as_vector());
4195  for (size_type i = 0; i < meshdim; ++i)
4196  pnode->tensor()(i,i) = scalar_type(1);
4197  }
4198  break;
4199 
4200  case GA_NODE_NORMAL:
4201  case GA_NODE_INTERPOLATE_NORMAL:
4202  GMM_ASSERT1(false, "Sorry, Gradient of Normal vector not implemented");
4203  break;
4204 
4205  case GA_NODE_ELT_K: case GA_NODE_ELT_B:
4206  case GA_NODE_INTERPOLATE_ELT_K: case GA_NODE_INTERPOLATE_ELT_B:
4207  GMM_ASSERT1(false, "Sorry, Gradient of element_K or element_B "
4208  "not implemented");
4209  break;
4210 
4211  case GA_NODE_INTERPOLATE_DERIVATIVE:
4212  GMM_ASSERT1(false, "Sorry, gradient of the derivative of a "
4213  "tranformation not implemented");
4214  break;
4215 
4216  case GA_NODE_INTERPOLATE_FILTER:
4217  ga_node_grad(tree, workspace, m, child0);
4218  break;
4219 
4220  case GA_NODE_ELEMENTARY_VAL: case GA_NODE_ELEMENTARY_VAL_TEST:
4221  case GA_NODE_XFEM_PLUS_VAL: case GA_NODE_XFEM_PLUS_VAL_TEST:
4222  case GA_NODE_XFEM_MINUS_VAL: case GA_NODE_XFEM_MINUS_VAL_TEST:
4223  case GA_NODE_ELEMENTARY_GRAD: case GA_NODE_ELEMENTARY_GRAD_TEST:
4224  case GA_NODE_XFEM_PLUS_GRAD: case GA_NODE_XFEM_PLUS_GRAD_TEST:
4225  case GA_NODE_XFEM_MINUS_GRAD: case GA_NODE_XFEM_MINUS_GRAD_TEST:
4226  {
4227  static const std::map<GA_NODE_TYPE, GA_NODE_TYPE>
4228  replacement_table =
4229  { {GA_NODE_ELEMENTARY_VAL, GA_NODE_ELEMENTARY_GRAD},
4230  {GA_NODE_ELEMENTARY_VAL_TEST, GA_NODE_ELEMENTARY_GRAD_TEST},
4231  {GA_NODE_ELEMENTARY_GRAD, GA_NODE_ELEMENTARY_HESS},
4232  {GA_NODE_ELEMENTARY_GRAD_TEST, GA_NODE_ELEMENTARY_HESS_TEST},
4233  {GA_NODE_XFEM_PLUS_VAL, GA_NODE_XFEM_PLUS_GRAD},
4234  {GA_NODE_XFEM_PLUS_VAL_TEST, GA_NODE_XFEM_PLUS_GRAD_TEST},
4235  {GA_NODE_XFEM_PLUS_GRAD, GA_NODE_XFEM_PLUS_HESS},
4236  {GA_NODE_XFEM_PLUS_GRAD_TEST, GA_NODE_XFEM_PLUS_HESS_TEST},
4237  {GA_NODE_XFEM_MINUS_VAL, GA_NODE_XFEM_MINUS_GRAD},
4238  {GA_NODE_XFEM_MINUS_VAL_TEST, GA_NODE_XFEM_MINUS_GRAD_TEST},
4239  {GA_NODE_XFEM_MINUS_GRAD, GA_NODE_XFEM_MINUS_HESS},
4240  {GA_NODE_XFEM_MINUS_GRAD_TEST, GA_NODE_XFEM_MINUS_HESS_TEST}
4241  };
4242  pnode->node_type = replacement_table.at(pnode->node_type);
4243  }
4244  mi = pnode->tensor().sizes();
4245  if (mi.back() == 1)
4246  mi.back() = m.dim();
4247  else
4248  mi.push_back(m.dim());
4249  pnode->t.adjust_sizes(mi);
4250  break;
4251  case GA_NODE_ELEMENTARY_HESS: case GA_NODE_ELEMENTARY_HESS_TEST:
4252  case GA_NODE_XFEM_PLUS_HESS: case GA_NODE_XFEM_PLUS_HESS_TEST:
4253  case GA_NODE_XFEM_MINUS_HESS: case GA_NODE_XFEM_MINUS_HESS_TEST:
4254  GMM_ASSERT1(false, "Sorry, cannot derive an Hessian once more");
4255  break;
4256  case GA_NODE_ELEMENTARY_DIVERG: case GA_NODE_ELEMENTARY_DIVERG_TEST:
4257  case GA_NODE_XFEM_PLUS_DIVERG: case GA_NODE_XFEM_PLUS_DIVERG_TEST:
4258  case GA_NODE_XFEM_MINUS_DIVERG: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
4259  {
4260  static const std::map<GA_NODE_TYPE, GA_NODE_TYPE>
4261  replacement_table =
4262  { {GA_NODE_ELEMENTARY_DIVERG, GA_NODE_ELEMENTARY_HESS},
4263  {GA_NODE_ELEMENTARY_DIVERG_TEST, GA_NODE_ELEMENTARY_HESS_TEST},
4264  {GA_NODE_XFEM_PLUS_DIVERG, GA_NODE_XFEM_PLUS_HESS},
4265  {GA_NODE_XFEM_PLUS_DIVERG_TEST, GA_NODE_XFEM_PLUS_HESS_TEST},
4266  {GA_NODE_XFEM_MINUS_DIVERG, GA_NODE_XFEM_MINUS_HESS},
4267  {GA_NODE_XFEM_MINUS_DIVERG_TEST, GA_NODE_XFEM_MINUS_HESS_TEST}
4268  };
4269  pnode->node_type = replacement_table.at(pnode->node_type);
4270  }
4271  mi = pnode->tensor().sizes();
4272  mi.pop_back();
4273  mi.push_back(m.dim());
4274  if (m.dim() > 1)
4275  mi.push_back(m.dim());
4276  pnode->t.adjust_sizes(mi);
4277  tree.duplicate_with_operation(pnode, GA_COLON);
4278  child0 = pnode;
4279  pnode = pnode->parent;
4280  child1 = pnode->children[1];
4281  child1->init_matrix_tensor(meshdim, meshdim);
4282  gmm::clear(pnode->tensor().as_vector());
4283  for (size_type i = 0; i < meshdim; ++i)
4284  child1->tensor()(i,i) = scalar_type(1);
4285  child1->node_type = GA_NODE_CONSTANT;
4286  break;
4287 
4288  case GA_NODE_OP:
4289  switch(pnode->op_type) {
4290  case GA_PLUS: case GA_MINUS:
4291  if (mark0 && mark1) {
4292  ga_node_grad(tree, workspace, m, child0);
4293  ga_node_grad(tree, workspace, m, child1);
4294  } else if (mark0) {
4295  ga_node_grad(tree, workspace, m, child0);
4296  tree.replace_node_by_child(pnode, 0);
4297  } else {
4298  ga_node_grad(tree, workspace, m, child1);
4299  if (pnode->op_type == GA_MINUS) {
4300  pnode->op_type = GA_UNARY_MINUS;
4301  tree.clear_node(child0);
4302  }
4303  else
4304  tree.replace_node_by_child(pnode, 1);
4305  }
4306  break;
4307 
4308  case GA_UNARY_MINUS: case GA_PRINT:
4309  ga_node_grad(tree, workspace, m, child0);
4310  break;
4311 
4312  case GA_QUOTE:
4313  if (child0->tensor_order() == 1) {
4314  size_type nn = child0->tensor_proper_size(0);
4315  ga_node_grad(tree, workspace, m, child0);
4316  pnode->node_type = GA_NODE_PARAMS;
4317  tree.add_child(pnode);
4318  std::swap(pnode->children[0], pnode->children[1]);
4319  pnode->children[0]->node_type = GA_NODE_RESHAPE;
4320  tree.add_child(pnode); tree.add_child(pnode); tree.add_child(pnode);
4321  pnode->children[2]->node_type = GA_NODE_CONSTANT;
4322  pnode->children[3]->node_type = GA_NODE_CONSTANT;
4323  pnode->children[4]->node_type = GA_NODE_CONSTANT;
4324  pnode->children[2]->init_scalar_tensor(scalar_type(1));
4325  pnode->children[3]->init_scalar_tensor(scalar_type(nn));
4326  pnode->children[4]->init_scalar_tensor(scalar_type(m.dim()));
4327  } else {
4328  ga_node_grad(tree, workspace, m, child0);
4329  }
4330  break;
4331 
4332  case GA_SYM: // Replace Sym(T) by (T+T')*0.5
4333  case GA_SKEW: // Replace Skew(T) by (T-T')*0.5
4334  if (pnode->op_type == GA_SYM) {
4335  tree.replace_node_by_child(pnode, 0); // cannot be before the if
4336  tree.duplicate_with_addition(child0);
4337  } else { // if (pnode->node_type == GA_SKEW)
4338  tree.replace_node_by_child(pnode, 0); // cannot be before the if
4339  tree.duplicate_with_substraction(child0);
4340  }
4341  tree.insert_node(child0->parent, GA_NODE_OP);
4342  tree.add_child(child0->parent->parent);
4343  child0->parent->parent->op_type = GA_MULT;
4344  child0->parent->parent->children[1]->node_type = GA_NODE_CONSTANT;
4345  child0->parent->parent->children[1]->init_scalar_tensor(0.5);
4346  tree.insert_node(child0->parent->children[1], GA_NODE_OP);
4347  child0->parent->children[1]->op_type = GA_QUOTE;
4348  ga_node_grad(tree, workspace, m, child0);
4349  ga_node_grad(tree, workspace, m,
4350  child0->parent->children[1]->children[0]);
4351  break;
4352 
4353  case GA_TRACE:
4354  ga_node_grad(tree, workspace, m, child0);
4355  pnode->node_type = GA_NODE_PARAMS;
4356  tree.add_child(pnode);
4357  std::swap(pnode->children[0], pnode->children[1]);
4358  pnode->children[0]->node_type = GA_NODE_NAME;
4359  pnode->children[0]->name = "Contract";
4360  tree.add_child(pnode); tree.add_child(pnode);
4361  pnode->children[2]->node_type = GA_NODE_CONSTANT;
4362  pnode->children[3]->node_type = GA_NODE_CONSTANT;
4363  pnode->children[2]->init_scalar_tensor(scalar_type(1));
4364  pnode->children[3]->init_scalar_tensor(scalar_type(2));
4365  break;
4366 
4367  case GA_DEVIATOR: // Replace Deviator(T) by (T-Trace(T)*Id(mdim)/mdim)
4368  {
4369  tree.duplicate_with_substraction(child0);
4370  child1 = child0->parent->children[1];
4371  tree.insert_node(child1, GA_NODE_OP);
4372  child1->parent->op_type = GA_TRACE;
4373  tree.insert_node(child1->parent, GA_NODE_OP);
4374  child1->parent->parent->op_type = GA_TMULT;
4375  tree.add_child(child1->parent->parent);
4376  std::swap(child1->parent->parent->children[0],
4377  child1->parent->parent->children[1]);
4378  pga_tree_node pid = child1->parent->parent->children[0];
4379  tree.duplicate_with_operation(pid, GA_DIV);
4380  pid->parent->children[1]->node_type = GA_NODE_CONSTANT;
4381  pid->parent->children[1]->init_scalar_tensor(m.dim());
4382  pid->node_type = GA_NODE_PARAMS;
4383  tree.add_child(pid); tree.add_child(pid);
4384  pid->children[0]->node_type = GA_NODE_NAME;
4385  pid->children[0]->name = "Id";
4386  pid->children[1]->node_type = GA_NODE_CONSTANT;
4387  pid->children[1]->init_scalar_tensor(m.dim());
4388  ga_node_grad(tree, workspace, m, child0);
4389  child1->parent->marked = true;
4390  ga_node_grad(tree, workspace, m, child1->parent);
4391  tree.replace_node_by_child(pnode, 0);
4392  }
4393  break;
4394 
4395 
4396  case GA_DOT: case GA_MULT: case GA_DOTMULT: case GA_COLON: case GA_TMULT:
4397  {
4398  pga_tree_node pg1(0), pg2(0);
4399  if (mark0 && mark1) {
4400  if (sub_tree_are_equal(child0, child1, workspace, 0) &&
4401  (pnode->op_type != GA_MULT || child0->tensor_order() < 2) &&
4402  (pnode->op_type != GA_DOT || child0->tensor_order() < 2)) {
4403  pg2 = pnode;
4404  tree.insert_node(pnode, GA_NODE_OP);
4405  pnode->parent->op_type = GA_MULT;
4406  tree.add_child(pnode->parent);
4407  pnode->parent->children[1]->node_type = GA_NODE_CONSTANT;
4408  pnode->parent->children[1]->init_scalar_tensor(scalar_type(2));
4409  } else {
4410  tree.duplicate_with_addition(pnode);
4411  pg1 = pnode;
4412  pg2 = pnode->parent->children[1];
4413  }
4414  } else if (mark0) pg1 = pnode; else pg2 = pnode;
4415 
4416  if (pg1) {
4417  if ((pg1->op_type == GA_COLON && child0->tensor_order() == 2) ||
4418  (pg1->op_type == GA_DOT && child0->tensor_order() == 1 &&
4419  child1->tensor_order() == 1) ||
4420  pg1->op_type == GA_DOTMULT ||
4421  (child0->tensor_proper_size()== 1 ||
4422  child1->tensor_proper_size()== 1)) {
4423  std::swap(pg1->children[0], pg1->children[1]);
4424  } else {
4425  size_type nred = 0;
4426  if (pnode->op_type == GA_MULT || pnode->op_type == GA_COLON ||
4427  pnode->op_type == GA_DOT) {
4428  if ((pg1->children[0]->tensor_order() <= 2 &&
4429  pnode->op_type == GA_MULT) ||
4430  pnode->op_type == GA_DOT) {
4431  nred = pg1->children[0]->tensor_order();
4432  pg1->node_type = GA_NODE_PARAMS;
4433  tree.add_child(pg1);tree.add_child(pg1);tree.add_child(pg1);
4434  std::swap(pg1->children[1], pg1->children[3]);
4435  std::swap(pg1->children[0], pg1->children[1]);
4436  pg1->children[0]->node_type = GA_NODE_NAME;
4437  pg1->children[0]->name = "Contract";
4438  pg1->children[2]->node_type = GA_NODE_CONSTANT;
4439  pg1->children[2]->init_scalar_tensor
4440  (scalar_type(pg1->children[1]->tensor_order()));
4441  pg1->children[4]->node_type = GA_NODE_CONSTANT;
4442  pg1->children[4]->init_scalar_tensor(scalar_type(1));
4443  ga_node_grad(tree, workspace, m, pg1->children[1]);
4444  } else {
4445  nred = pg1->children[0]->tensor_order()-1;
4446  pg1->node_type = GA_NODE_PARAMS;
4447  tree.add_child(pg1);tree.add_child(pg1);tree.add_child(pg1);
4448  tree.add_child(pg1);tree.add_child(pg1);
4449  std::swap(pg1->children[1], pg1->children[4]);
4450  std::swap(pg1->children[0], pg1->children[1]);
4451  pg1->children[0]->node_type = GA_NODE_NAME;
4452  pg1->children[0]->name = "Contract";
4453  pg1->children[2]->node_type = GA_NODE_CONSTANT;
4454  pg1->children[2]->init_scalar_tensor
4455  (scalar_type(pg1->children[1]->tensor_order()-1));
4456  pg1->children[3]->node_type = GA_NODE_CONSTANT;
4457  pg1->children[3]->init_scalar_tensor
4458  (scalar_type(pg1->children[1]->tensor_order()));
4459  pg1->children[5]->node_type = GA_NODE_CONSTANT;
4460  pg1->children[5]->init_scalar_tensor(scalar_type(1));
4461  pg1->children[6]->node_type = GA_NODE_CONSTANT;
4462  pg1->children[6]->init_scalar_tensor(scalar_type(2));
4463  ga_node_grad(tree, workspace, m, pg1->children[1]);
4464  }
4465  } else if (pnode->op_type == GA_TMULT) {
4466  nred = pg1->children[0]->tensor_order()+1;
4467  ga_node_grad(tree, workspace, m, pg1->children[0]);
4468  } else GMM_ASSERT1(false, "internal error");
4469  tree.insert_node(pg1, GA_NODE_PARAMS);
4470  tree.add_child(pg1->parent);
4471  tree.add_child(pg1->parent);
4472  std::swap(pg1->parent->children[0], pg1->parent->children[1]);
4473  pg1->parent->children[0]->node_type = GA_NODE_IND_MOVE_LAST;
4474  pg1->parent->children[2]->node_type = GA_NODE_CONSTANT;
4475  pg1->parent->children[2]->init_scalar_tensor(scalar_type(nred));
4476  pg1 = 0;
4477  }
4478  }
4479 
4480  for (; pg2||pg1;pg2=pg1, pg1=0) {
4481  if (pg2) {
4482  if (pnode->op_type == GA_MULT || pnode->op_type == GA_COLON ||
4483  pnode->op_type == GA_DOT) {
4484  if (pg2->children[1]->tensor_proper_size() == 1) {
4485  if (pg2->children[0]->tensor_proper_size() == 1)
4486  pg2->op_type = GA_MULT;
4487  else
4488  pg2->op_type = GA_TMULT;
4489  ga_node_grad(tree, workspace, m, pg2->children[1]);
4490  } else if (pg2->children[0]->tensor_proper_size() == 1) {
4491  pg2->op_type = GA_MULT;
4492  ga_node_grad(tree, workspace, m, pg2->children[1]);
4493  } else if ((pg2->children[0]->tensor_order() <= 2 &&
4494  pnode->op_type == GA_MULT) ||
4495  pnode->op_type == GA_DOT) {
4496  pg2->op_type = GA_DOT;
4497  ga_node_grad(tree, workspace, m, pg2->children[1]);
4498  } else {
4499  pg2->node_type = GA_NODE_PARAMS;
4500  tree.add_child(pg2);tree.add_child(pg2);tree.add_child(pg2);
4501  tree.add_child(pg2);tree.add_child(pg2);
4502  std::swap(pg2->children[1], pg2->children[4]);
4503  std::swap(pg2->children[0], pg2->children[1]);
4504  pg2->children[0]->node_type = GA_NODE_NAME;
4505  pg2->children[0]->name = "Contract";
4506  pg2->children[2]->node_type = GA_NODE_CONSTANT;
4507  pg2->children[2]->init_scalar_tensor
4508  (scalar_type(pg2->children[4]->tensor_order()-1));
4509  pg2->children[3]->node_type = GA_NODE_CONSTANT;
4510  pg2->children[3]->init_scalar_tensor
4511  (scalar_type(pg2->children[4]->tensor_order()));
4512  pg2->children[5]->node_type = GA_NODE_CONSTANT;
4513  pg2->children[5]->init_scalar_tensor(scalar_type(1));
4514  pg2->children[6]->node_type = GA_NODE_CONSTANT;
4515  pg2->children[6]->init_scalar_tensor(scalar_type(2));
4516  ga_node_grad(tree, workspace, m, pg2->children[4]);
4517  }
4518  } else if (pnode->op_type == GA_TMULT) {
4519  ga_node_grad(tree, workspace, m, pg2->children[1]);
4520  } else if (pnode->op_type == GA_DOTMULT) {
4521  if (pg2->children[0]->tensor_proper_size() == 1) {
4522  pg2->op_type = GA_MULT;
4523  ga_node_grad(tree, workspace, m, pg2->children[1]);
4524  } else {
4525  tree.insert_node(pg2->children[0], GA_NODE_OP);
4526  tree.add_child(pg2->children[0], GA_NODE_CONSTANT);
4527  pg2->children[0]->op_type = GA_TMULT;
4528  pg2->children[0]->children[1]->init_vector_tensor(m.dim());
4529  gmm::fill(pg2->children[0]->children[1]->tensor().as_vector(),
4530  scalar_type(1));
4531  ga_node_grad(tree, workspace, m, pg2->children[1]);
4532  }
4533  }
4534  }
4535  }
4536  }
4537  break;
4538 
4539  case GA_DIV: case GA_DOTDIV:
4540  {
4541  pga_tree_node pg1 = child0;
4542  if (mark1) {
4543  if (pnode->children[0]->node_type == GA_NODE_CONSTANT) {
4544  gmm::scale(pnode->children[0]->tensor().as_vector(),
4545  scalar_type(-1));
4546  pg1=0;
4547  } else {
4548  if (mark0) {
4549  tree.duplicate_with_substraction(pnode);
4550  pnode = pnode->parent->children[1];
4551  pg1 = child0;
4552  } else {
4553  tree.insert_node(pnode, GA_NODE_OP);
4554  pnode->parent->op_type = GA_UNARY_MINUS;
4555  pg1 = 0;
4556  }
4557  }
4558  tree.insert_node(pnode->children[1], GA_NODE_PARAMS);
4559  pga_tree_node pnode_param = pnode->children[1];
4560  tree.add_child(pnode_param);
4561  std::swap(pnode_param->children[0], pnode_param->children[1]);
4562  pnode_param->children[0]->node_type = GA_NODE_PREDEF_FUNC;
4563  pnode_param->children[0]->name = "sqr";
4564  tree.insert_node(pnode, GA_NODE_OP);
4565  pga_tree_node pnode_mult = pnode->parent;
4566  if (pnode->op_type == GA_DOTDIV) {
4567  pnode_mult->op_type = GA_DOTMULT;
4568  tree.insert_node(pnode_mult->children[0], GA_NODE_OP);
4569  pga_tree_node ptmult = pnode_mult->children[0];
4570  ptmult->op_type = GA_TMULT;
4571  tree.add_child(ptmult, GA_NODE_CONSTANT);
4572  ptmult->children[1]->init_vector_tensor(m.dim());
4573  gmm::fill(ptmult->children[1]->tensor().as_vector(),
4574  scalar_type(1));
4575  } else {
4576  pnode_mult->op_type = GA_TMULT;
4577  }
4578  pnode_mult->children.resize(2, nullptr);
4579  tree.copy_node(pnode_param->children[1],
4580  pnode_mult->children[1]);
4581  pnode_mult->accept_child(1);
4582  ga_node_grad(tree, workspace, m, pnode_mult->children[1]);
4583  }
4584 
4585  if (pg1) {
4586  ga_node_grad(tree, workspace, m, pg1);
4587  if (pnode->op_type == GA_DOTDIV) {
4588  tree.insert_node(pg1->parent->children[1], GA_NODE_OP);
4589  pga_tree_node ptmult = pg1->parent->children[1];
4590  ptmult->op_type = GA_TMULT;
4591  tree.add_child(ptmult, GA_NODE_CONSTANT);
4592  ptmult->children[1]->init_vector_tensor(m.dim());
4593  gmm::fill(ptmult->children[1]->tensor().as_vector(),
4594  scalar_type(1));
4595  }
4596  }
4597  }
4598  break;
4599  default: GMM_ASSERT1(false, "Unexpected operation. Internal error.");
4600  }
4601  break;
4602 
4603  case GA_NODE_C_MATRIX:
4604  {
4605  for (size_type i = 0; i < pnode->children.size(); ++i) {
4606  if (pnode->children[i]->marked)
4607  ga_node_grad(tree, workspace, m, pnode->children[i]);
4608  else {
4609  pnode->children[i]->init_scalar_tensor(scalar_type(0));
4610  pnode->children[i]->node_type = GA_NODE_ZERO;
4611  tree.clear_children(pnode->children[i]);
4612  }
4613  }
4614  if (m.dim() > 1) {
4615  size_type orgsize = pnode->children.size();
4616  mi = pnode->tensor().sizes();
4617  mi.push_back(m.dim());
4618  pnode->nbc1 += 1;
4619  pnode->tensor().adjust_sizes(mi);
4620 
4621  pnode->children.resize(orgsize*m.dim(), nullptr);
4622  for (size_type i = orgsize; i < pnode->children.size(); ++i) {
4623  tree.copy_node(pnode->children[i % orgsize],
4624  pnode->children[i]);
4625  pnode->accept_child(i);
4626  }
4627  for (size_type i = 0; i < pnode->children.size(); ++i) {
4628  pga_tree_node child = pnode->children[i];
4629  if (child->node_type != GA_NODE_ZERO) {
4630  tree.insert_node(child, GA_NODE_PARAMS);
4631  tree.add_child(child->parent, GA_NODE_CONSTANT);
4632  child->parent->children[1]
4633  ->init_scalar_tensor(scalar_type(1+i/orgsize));
4634  }
4635  }
4636  }
4637  }
4638  break;
4639 
4640  case GA_NODE_PARAMS:
4641  if (child0->node_type == GA_NODE_RESHAPE) {
4642  ga_node_grad(tree, workspace, m, pnode->children[1]);
4643  tree.add_child(pnode, GA_NODE_CONSTANT);
4644  pnode->children.back()->init_scalar_tensor(scalar_type(m.dim()));
4645  } else if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
4646  GMM_ASSERT1(false, "Sorry, the gradient of a cross product"
4647  "has not been implemented. To be done.");
4648  } else if (child0->node_type == GA_NODE_IND_MOVE_LAST) {
4649  size_type order = pnode->tensor_order();
4650  ga_node_grad(tree, workspace, m, pnode->children[1]);
4651  tree.insert_node(pnode, GA_NODE_PARAMS);
4652  tree.add_child(pnode->parent);
4653  tree.add_child(pnode->parent);
4654  tree.add_child(pnode->parent);
4655  std::swap(pnode->parent->children[0], pnode->parent->children[1]);
4656  pnode->parent->children[0]->node_type = GA_NODE_SWAP_IND;
4657  pnode->parent->children[2]->node_type = GA_NODE_CONSTANT;
4658  pnode->parent->children[3]->node_type = GA_NODE_CONSTANT;
4659  pnode->parent->children[2]->init_scalar_tensor(scalar_type(order));
4660  pnode->parent->children[3]->init_scalar_tensor(scalar_type(order+1));
4661  } else if (child0->node_type == GA_NODE_SWAP_IND) {
4662  ga_node_grad(tree, workspace, m, pnode->children[1]);
4663  } else if (child0->node_type == GA_NODE_CONTRACT) {
4664  mark0 = mark1;
4665  size_type ch2 = 0;
4666  if (pnode->children.size() == 5) ch2 = 3;
4667  if (pnode->children.size() == 7) ch2 = 4;
4668  mark1 = pnode->children[ch2]->marked;
4669 
4670  if (pnode->children.size() == 4) {
4671  ga_node_grad(tree, workspace, m, pnode->children[1]);
4672  } else {
4673  pga_tree_node pg1(pnode), pg2(pnode);
4674  if (mark0 && mark1) {
4675  tree.duplicate_with_addition(pnode);
4676  pg2 = pnode->parent->children[1];
4677  }
4678  if (mark0) {
4679  size_type nred = pg1->children[1]->tensor_order();
4680  if (pnode->children.size() == 7) nred--;
4681  ga_node_grad(tree, workspace, m, pg1->children[1]);
4682  tree.insert_node(pg1, GA_NODE_PARAMS);
4683  tree.add_child(pg1->parent);
4684  tree.add_child(pg1->parent);
4685  std::swap(pg1->parent->children[0], pg1->parent->children[1]);
4686  pg1->parent->children[0]->node_type = GA_NODE_IND_MOVE_LAST;
4687  pg1->parent->children[2]->node_type = GA_NODE_CONSTANT;
4688  pg1->parent->children[2]->init_scalar_tensor(scalar_type(nred));
4689  }
4690  if (mark1) {
4691  ga_node_grad(tree, workspace, m, pg2->children[ch2]);
4692  }
4693  }
4694 
4695  } else if (child0->node_type == GA_NODE_PREDEF_FUNC) {
4696  std::string name = child0->name;
4697  ga_predef_function_tab::const_iterator it = PREDEF_FUNCTIONS.find(name);
4698  const ga_predef_function &F = it->second;
4699 
4700  if (F.nbargs() == 1) {
4701  switch (F.dtype()) {
4702  case 0:
4703  GMM_ASSERT1(false, "Cannot derive function " << child0->name
4704  << ". No derivative provided or not derivable function.");
4705  break;
4706  case 1:
4707  child0->name = F.derivative1();
4708  break;
4709  case 2: case 3:
4710  {
4711  child0->name = "DER_PDFUNC_" + child0->name;
4712  if (F.dtype() == 2)
4713  ga_define_function(child0->name, 1, F.derivative1());
4714  else {
4715  std::string expr=ga_derivative_scalar_function(F.expr(),"t");
4716  ga_define_function(child0->name, 1, expr);
4717  }
4718  // Inline extension if the derivative is affine (for instance
4719  // for sqr)
4720  ga_predef_function_tab::const_iterator
4721  itp = PREDEF_FUNCTIONS.find(child0->name);
4722  const ga_predef_function &Fp = itp->second;
4723  if (Fp.is_affine("t")) {
4724  scalar_type b = Fp(scalar_type(0));
4725  scalar_type a = Fp(scalar_type(1)) - b;
4726  pnode->node_type = GA_NODE_OP;
4727  pnode->op_type = GA_MULT;
4728  child0->init_scalar_tensor(a);
4729  child0->node_type = ((a == scalar_type(0)) ?
4730  GA_NODE_ZERO : GA_NODE_CONSTANT);
4731  if (b != scalar_type(0)) {
4732  tree.insert_node(pnode, GA_NODE_OP);
4733  pnode->parent->op_type = (b > 0) ? GA_PLUS : GA_MINUS;
4734  tree.add_child(pnode->parent);
4735  pga_tree_node pnode_cte = pnode->parent->children[1];
4736  pnode_cte->node_type = GA_NODE_CONSTANT;
4737  pnode_cte->t = pnode->t;
4738  std::fill(pnode_cte->tensor().begin(),
4739  pnode_cte->tensor().end(), gmm::abs(b));
4740  pnode = pnode->parent;
4741  }
4742  }
4743  }
4744  break;
4745  }
4746  if (pnode->children.size() >= 2) {
4747  tree.insert_node(pnode, GA_NODE_OP);
4748  pga_tree_node pnode_op = pnode->parent;
4749  if (child1->tensor_order() == 0)
4750  pnode_op->op_type = GA_MULT;
4751  else {
4752  pnode_op->op_type = GA_DOTMULT;
4753  tree.insert_node(pnode, GA_NODE_OP);
4754  pnode->parent->op_type = GA_TMULT;
4755  tree.add_child(pnode->parent, GA_NODE_CONSTANT);
4756  pnode->parent->children[1]->init_vector_tensor(m.dim());
4757  gmm::fill(pnode->parent->children[1]->tensor().as_vector(),
4758  scalar_type(1));
4759  }
4760  pnode_op->children.resize(2, nullptr);
4761  tree.copy_node(child1, pnode_op->children[1]);
4762  pnode_op->accept_child(1);
4763  ga_node_grad(tree, workspace, m, pnode_op->children[1]);
4764  }
4765  } else {
4766  pga_tree_node child2 = pnode->children[2];
4767  pga_tree_node pg2 = pnode;
4768 
4769  if (child1->marked && child2->marked) {
4770  tree.duplicate_with_addition(pnode);
4771  pg2 = pnode->parent->children[1];
4772  }
4773 
4774  if (child1->marked) {
4775  switch (F.dtype()) {
4776  case 0:
4777  GMM_ASSERT1(false, "Cannot derive function " << child0->name
4778  << ". No derivative provided");
4779  break;
4780  case 1:
4781  child0->name = F.derivative1();
4782  break;
4783  case 2:
4784  child0->name = "DER_PDFUNC1_" + child0->name;
4785  ga_define_function(child0->name, 2, F.derivative1());
4786  break;
4787  case 3:
4788  child0->name = "DER_PDFUNC1_" + child0->name;
4789  std::string expr = ga_derivative_scalar_function(F.expr(), "t");
4790  ga_define_function(child0->name, 2, expr);
4791  break;
4792  }
4793  tree.insert_node(pnode, GA_NODE_OP);
4794  pga_tree_node pnode_op = pnode->parent;
4795  if (child1->tensor_order() == 0)
4796  pnode_op->op_type = GA_MULT;
4797  else {
4798  pnode_op->op_type = GA_DOTMULT;
4799  tree.insert_node(pnode, GA_NODE_OP);
4800  pnode->parent->op_type = GA_TMULT;
4801  tree.add_child(pnode->parent, GA_NODE_CONSTANT);
4802  pnode->parent->children[1]->init_vector_tensor(m.dim());
4803  gmm::fill(pnode->parent->children[1]->tensor().as_vector(),
4804  scalar_type(1));
4805  }
4806  pnode_op->children.resize(2, nullptr);
4807  tree.copy_node(child1, pnode_op->children[1]);
4808  pnode_op->accept_child(1);
4809  ga_node_grad(tree, workspace, m, pnode_op->children[1]);
4810  }
4811  if (child2->marked) {
4812  pnode = pg2;
4813  child0 = pnode->children[0]; child1 = pnode->children[1];
4814  child2 = pnode->children[2];
4815 
4816  switch (F.dtype()) {
4817  case 0:
4818  GMM_ASSERT1(false, "Cannot derive function " << child0->name
4819  << ". No derivative provided");
4820  break;
4821  case 1:
4822  child0->name = F.derivative2();
4823  break;
4824  case 2:
4825  child0->name = "DER_PDFUNC2_" + child0->name;
4826  ga_define_function(child0->name, 2, F.derivative2());
4827  break;
4828  case 3:
4829  child0->name = "DER_PDFUNC2_" + child0->name;
4830  std::string expr = ga_derivative_scalar_function(F.expr(), "u");
4831  ga_define_function(child0->name, 2, expr);
4832  break;
4833  }
4834  tree.insert_node(pnode, GA_NODE_OP);
4835  pga_tree_node pnode_op = pnode->parent;
4836  if (child2->tensor_order() == 0)
4837  pnode_op->op_type = GA_MULT;
4838  else {
4839  pnode_op->op_type = GA_DOTMULT;
4840  tree.insert_node(pnode, GA_NODE_OP);
4841  pnode->parent->op_type = GA_TMULT;
4842  tree.add_child(pnode->parent, GA_NODE_CONSTANT);
4843  pnode->parent->children[1]->init_vector_tensor(m.dim());
4844  gmm::fill(pnode->parent->children[1]->tensor().as_vector(),
4845  scalar_type(1));
4846  }
4847  pnode_op->children.resize(2, nullptr);
4848  tree.copy_node(child2, pnode_op->children[1]);
4849  pnode_op->accept_child(1);
4850  ga_node_grad(tree, workspace, m, pnode_op->children[1]);
4851  }
4852  }
4853  } else if (child0->node_type == GA_NODE_SPEC_FUNC) {
4854  GMM_ASSERT1(false, "internal error");
4855  } else if (child0->node_type == GA_NODE_OPERATOR) {
4856  if (child0->der2)
4857  GMM_ASSERT1(false, "Error in derivation of the assembly string. "
4858  "Cannot derive again operator " << child0->name);
4859 
4860  size_type nbargs_der = 0;
4861  for (size_type i = 1; i < pnode->children.size(); ++i)
4862  if (pnode->children[i]->marked) ++nbargs_der;
4863  pga_tree_node pnode2 = 0;
4864 
4865  size_type j = 0;
4866  for (size_type i = 1; i < pnode->children.size(); ++i) {
4867  if (pnode->children[i]->marked) {
4868  ++j;
4869  if (j != nbargs_der) {
4870  tree.insert_node(pnode, GA_NODE_OP);
4871  pga_tree_node pnode_op = pnode->parent;
4872  pnode_op->node_type = GA_NODE_OP;
4873  pnode_op->op_type = GA_PLUS;
4874  pnode_op->children.resize(2, nullptr);
4875  tree.copy_node(pnode, pnode_op->children[1]);
4876  pnode_op->accept_child(1);
4877  pnode2 = pnode_op->children[1];
4878  }
4879  else pnode2 = pnode;
4880 
4881  if (child0->der1)
4882  pnode2->children[0]->der2 = i;
4883  else
4884  pnode2->children[0]->der1 = i;
4885  tree.insert_node(pnode2, GA_NODE_OP);
4886  pga_tree_node pnode_op = pnode2->parent;
4887  // calcul de l'ordre de reduction
4888  size_type red = pnode->children[i]->tensor_order();
4889  switch (red) {
4890  case 0 : pnode_op->op_type = GA_MULT; break;
4891  case 1 : pnode_op->op_type = GA_DOT; break;
4892  case 2 : pnode_op->op_type = GA_COLON; break;
4893  default: GMM_ASSERT1(false, "Error in derivation of the assembly "
4894  "string. Bad reduction order.")
4895  }
4896  pnode_op->children.resize(2, nullptr);
4897  tree.copy_node(pnode->children[i], pnode_op->children[1]);
4898  pnode_op->accept_child(1);
4899  ga_node_grad(tree, workspace, m, pnode_op->children[1]);
4900 
4901  if (pnode2->children[0]->name.compare("Norm_sqr") == 0
4902  && pnode2->children[0]->der1 == 1) {
4903  pnode2->node_type = GA_NODE_OP;
4904  pnode2->op_type = GA_MULT;
4905  pnode2->children[0]->node_type = GA_NODE_CONSTANT;
4906  pnode2->children[0]->init_scalar_tensor(scalar_type(2));
4907  }
4908  }
4909  }
4910  } else {
4911  ga_node_grad(tree, workspace, m, child0);
4912  tree.add_child(pnode, GA_NODE_ALLINDICES);
4913  }
4914  break;
4915 
4916 
4917  default: GMM_ASSERT1(false, "Unexpected node type " << pnode->node_type
4918  << " in gradient. Internal error.");
4919  }
4920  // cout << "End of gradient of "; ga_print_node(pnode, cout); cout << endl;
4921  }
4922 
4923 
4924  // The tree is modified. Should be copied first and passed to
4925  // ga_semantic_analysis after for enrichment
4926  void ga_grad(ga_tree &tree, const ga_workspace &workspace, const mesh &m) {
4927  if (!(tree.root)) return;
4928  if (ga_node_mark_tree_for_grad(tree.root, workspace))
4929  ga_node_grad(tree, workspace, m, tree.root);
4930  else
4931  tree.clear();
4932  }
4933 
4934 
4935 
4936  static void ga_replace_test_by_cte(pga_tree_node pnode, bool full_replace) {
4937  for (size_type i = 0; i < pnode->children.size(); ++i)
4938  ga_replace_test_by_cte(pnode->children[i], full_replace);
4939  GMM_ASSERT1(pnode->node_type != GA_NODE_GRAD_TEST, "Invalid tree");
4940  GMM_ASSERT1(pnode->node_type != GA_NODE_HESS_TEST, "Invalid tree");
4941  GMM_ASSERT1(pnode->node_type != GA_NODE_DIVERG_TEST, "Invalid tree");
4942  if (pnode->node_type == GA_NODE_VAL_TEST) {
4943  pnode->node_type = GA_NODE_CONSTANT;
4944  if (full_replace) pnode->init_scalar_tensor(scalar_type(1));
4945  }
4946  }
4947 
4948  std::string ga_derivative_scalar_function(const std::string &expr,
4949  const std::string &var) {
4950  base_vector t(1), u(1);
4951  ga_workspace workspace;
4952  workspace.add_fixed_size_variable("t", gmm::sub_interval(0,1), t);
4953  workspace.add_fixed_size_variable("u", gmm::sub_interval(0,1), u);
4954  workspace.add_function_expression(expr);
4955  GMM_ASSERT1(workspace.nb_trees() <= 1, "Internal error");
4956  if (workspace.nb_trees()) {
4957  ga_tree tree = *(workspace.tree_info(0).ptree);
4958  ga_derivative(tree, workspace, dummy_mesh(), var, "", 1);
4959  if (tree.root) {
4960  ga_replace_test_by_cte(tree.root, true);
4961  ga_semantic_analysis(tree, workspace, dummy_mesh(), 1,
4962  false, true);
4963  }
4964  return ga_tree_to_string(tree);
4965  } else return "0";
4966  }
4967 
4968  static bool ga_node_is_affine(const pga_tree_node pnode) {
4969 
4970  size_type nbch = pnode->children.size();
4971  pga_tree_node child0 = (nbch > 0) ? pnode->children[0] : 0;
4972  pga_tree_node child1 = (nbch > 1) ? pnode->children[1] : 0;
4973  bool mark0 = ((nbch > 0) ? child0->marked : false);
4974  bool mark1 = ((nbch > 1) ? child1->marked : false);
4975 
4976  switch (pnode->node_type) {
4977  case GA_NODE_VAL: case GA_NODE_GRAD:
4978  case GA_NODE_HESS: case GA_NODE_DIVERG:
4979  case GA_NODE_INTERPOLATE_VAL: case GA_NODE_INTERPOLATE_GRAD:
4980  case GA_NODE_INTERPOLATE_HESS: case GA_NODE_INTERPOLATE_DIVERG:
4981  case GA_NODE_INTERPOLATE_DERIVATIVE:
4982  case GA_NODE_ELEMENTARY_VAL: case GA_NODE_ELEMENTARY_GRAD:
4983  case GA_NODE_ELEMENTARY_HESS: case GA_NODE_ELEMENTARY_DIVERG:
4984  case GA_NODE_SECONDARY_DOMAIN_VAL: case GA_NODE_SECONDARY_DOMAIN_GRAD:
4985  case GA_NODE_SECONDARY_DOMAIN_HESS: case GA_NODE_SECONDARY_DOMAIN_DIVERG:
4986  case GA_NODE_XFEM_PLUS_VAL: case GA_NODE_XFEM_PLUS_GRAD:
4987  case GA_NODE_XFEM_PLUS_HESS: case GA_NODE_XFEM_PLUS_DIVERG:
4988  case GA_NODE_XFEM_MINUS_VAL: case GA_NODE_XFEM_MINUS_GRAD:
4989  case GA_NODE_XFEM_MINUS_HESS: case GA_NODE_XFEM_MINUS_DIVERG:
4990  return true;
4991  case GA_NODE_INTERPOLATE_FILTER:
4992  return ga_node_is_affine(child0);
4993  case GA_NODE_OP:
4994  switch(pnode->op_type) {
4995  case GA_PLUS: case GA_MINUS:
4996  if (mark0 && mark1)
4997  return ga_node_is_affine(child0) &&
4998  ga_node_is_affine(child1);
4999  if (mark0) return ga_node_is_affine(child0);
5000  return ga_node_is_affine(child1);
5001 
5002  case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM: case GA_SKEW:
5003  case GA_TRACE: case GA_DEVIATOR: case GA_PRINT:
5004  return ga_node_is_affine(child0);
5005 
5006  case GA_DOT: case GA_MULT: case GA_COLON: case GA_TMULT:
5007  case GA_DOTMULT:
5008  if (mark0 && mark1) return false;
5009  if (mark0) return ga_node_is_affine(child0);
5010  return ga_node_is_affine(child1);
5011 
5012  case GA_DIV: case GA_DOTDIV:
5013  if (mark1) return false;
5014  if (mark0) return ga_node_is_affine(child0);
5015  return true;
5016 
5017  default: GMM_ASSERT1(false, "Unexpected operation. Internal error.");
5018  }
5019  break;
5020 
5021  case GA_NODE_C_MATRIX:
5022  for (size_type i = 0; i < pnode->children.size(); ++i)
5023  if (pnode->children[i]->marked &&
5024  !(ga_node_is_affine(pnode->children[i])))
5025  return false;
5026  return true;
5027 
5028  case GA_NODE_PARAMS:
5029  if (child0->node_type == GA_NODE_RESHAPE ||
5030  child0->node_type == GA_NODE_SWAP_IND ||
5031  child0->node_type == GA_NODE_IND_MOVE_LAST)
5032  return ga_node_is_affine(child1);
5033  if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
5034  pga_tree_node child2 = pnode->children[2];
5035  bool mark2 = child2->marked;
5036  if (mark1 && mark2) return false;
5037  if (mark1) return ga_node_is_affine(child1);
5038  return ga_node_is_affine(child2);
5039  }
5040  if (child0->node_type == GA_NODE_CONTRACT) {
5041  if (pnode->children.size() == 4) {
5042  return ga_node_is_affine(child1);
5043  } else if (pnode->children.size() == 5) {
5044  if (mark1 && pnode->children[3]->marked) return false;
5045  if (mark1) return ga_node_is_affine(child1);
5046  return ga_node_is_affine(pnode->children[3]);
5047  } else if (pnode->children.size() == 7) {
5048  if (mark1 && pnode->children[4]->marked) return false;
5049  if (mark1) return ga_node_is_affine(child1);
5050  return ga_node_is_affine(pnode->children[4]);
5051  }
5052  }
5053  if (child0->node_type == GA_NODE_PREDEF_FUNC)
5054  return false;
5055  if (child0->node_type == GA_NODE_OPERATOR)
5056  return false;
5057  return ga_node_is_affine(child0);
5058 
5059  default: GMM_ASSERT1(false, "Unexpected node type " << pnode->node_type
5060  << " in derivation. Internal error.");
5061  }
5062  }
5063 
5064  bool ga_is_affine(const ga_tree &tree, const ga_workspace &workspace,
5065  const std::string &varname,
5066  const std::string &interpolatename) {
5067  const mesh &m = dummy_mesh();
5068  if (tree.root && ga_node_mark_tree_for_variable(tree.root, workspace, m,
5069  varname, interpolatename))
5070  return ga_node_is_affine(tree.root);
5071  return true;
5072  }
5073 
5074 } /* end of namespace */
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:165
getfem_generic_assembly_semantic.h
Semantic analysis of assembly trees and semantic manipulations.
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:978
gmm::fill
void fill(L &l, typename gmm::linalg_traits< L >::value_type x)
*‍/
Definition: gmm_blas.h:104
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
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_accumulated_distro.h:47