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