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