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