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