GetFEM  5.4.1
getfem_generic_assembly_tree.h
Go to the documentation of this file.
1 /*===========================================================================
2 
3  Copyright (C) 2013-2020 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20  As a special exception, you may use this file as it is a part of a free
21  software library without restriction. Specifically, if other files
22  instantiate templates or use macros or inline functions from this file,
23  or you compile this file and link it with other files to produce an
24  executable, this file does not by itself cause the resulting executable
25  to be covered by the GNU Lesser General Public License. This exception
26  does not however invalidate any other reasons why the executable file
27  might be covered by the GNU Lesser General Public License.
28 
29 ===========================================================================*/
30 
31 /** @file getfem_generic_assembly_tree.h
32  @author Yves Renard <Yves.Renard@insa-lyon.fr>
33  @date November 18, 2013.
34  @brief Definition of the syntax tree and basic operations on it.
35  Internal header for the generic assembly language part.
36  */
37 
38 
39 #ifndef GETFEM_GENERIC_ASSEMBLY_TREE_H__
40 #define GETFEM_GENERIC_ASSEMBLY_TREE_H__
41 
42 #include "getfem/getfem_generic_assembly.h"
43 #include "getfem/getfem_models.h"
44 #include "gmm/gmm_blas.h"
45 #include <iomanip>
46 #include "getfem/getfem_omp.h"
47 #include "getfem/dal_singleton.h"
48 #include "getfem/bgeot_rtree.h"
49 #include "getfem/bgeot_geotrans_inv.h"
50 #include "getfem/getfem_copyable_ptr.h"
51 #ifndef _WIN32
52 extern "C"{
53 #include <unistd.h>
54 }
55 #endif
56 
57 #define GA_DEBUG_ASSERT(a, b) GMM_ASSERT1(a, b)
58 // #define GA_DEBUG_ASSERT(a, b)
59 
60 #if 1
61 # define GA_TIC
62 # define GA_TOC(a)
63 # define GA_TOCTIC(a)
64 #else
65 # define GA_TIC scalar_type _ga_time_ = gmm::uclock_sec();
66 # define GA_TOC(a) { cout <<(a)<<" : "<<gmm::uclock_sec()-_ga_time_<< endl; }
67 # define GA_TOCTIC(a) { GA_TOC(a); _ga_time_ = gmm::uclock_sec(); }
68 #endif
69 
70 namespace getfem {
71 
72  // Basic token types (basic language components)
73  enum GA_TOKEN_TYPE {
74  GA_INVALID = 0, // invalid token
75  GA_END, // string end
76  GA_NAME, // A variable or user defined nonlinear function name
77  GA_SCALAR, // A real number
78  GA_PLUS, // '+'
79  GA_MINUS, // '-'
80  GA_UNARY_MINUS, // '-'
81  GA_MULT, // '*'
82  GA_DIV, // '/'
83  GA_COLON, // ':'
84  GA_QUOTE, // ''' transpose
85  GA_COLON_EQ, // ':=' macro def
86  GA_DEF, // 'Def' macro def
87  GA_SYM, // 'Sym(M)' operator
88  GA_SKEW, // 'Skew(M)' operator
89  GA_TRACE, // 'Trace(M)' operator
90  GA_DEVIATOR, // 'Deviator' operator
91  GA_INTERPOLATE, // 'Interpolate' operation
92  GA_INTERPOLATE_FILTER, // 'Interpolate_filter' operation
93  GA_INTERPOLATE_DERIVATIVE, // 'Interpolate_derivative' operation
94  GA_ELEMENTARY, // 'Elementary' operation (operation at the element level)
95  GA_SECONDARY_DOMAIN, // For the integration on a product of two domains
96  GA_XFEM_PLUS, // Evaluation on the + side of a level-set for fem_level_set
97  GA_XFEM_MINUS, // Evaluation on the - side of a level-set for fem_level_set
98  GA_PRINT, // 'Print' Print the tensor
99  GA_DOT, // '.'
100  GA_DOTMULT, // '.*' componentwise multiplication
101  GA_DOTDIV, // './' componentwise division
102  GA_TMULT, // '@' tensor product
103  GA_COMMA, // ','
104  GA_DCOMMA, // ',,'
105  GA_SEMICOLON, // ';'
106  GA_DSEMICOLON, // ';;'
107  GA_LPAR, // '('
108  GA_RPAR, // ')'
109  GA_LBRACKET, // '['
110  GA_RBRACKET, // ']'
111  GA_NB_TOKEN_TYPE
112  };
113 
114  // Detects Grad_, Hess_ or Div_
115  size_type ga_parse_prefix_operator(std::string &name);
116  // Detects Test_ and Test2_
117  size_type ga_parse_prefix_test(std::string &name);
118 
119  // Types of nodes for the syntax tree
120  enum GA_NODE_TYPE {
121  GA_NODE_VOID = 0,
122  GA_NODE_OP,
123  GA_NODE_PREDEF_FUNC,
124  GA_NODE_SPEC_FUNC,
125  GA_NODE_OPERATOR,
126  GA_NODE_CONSTANT,
127  GA_NODE_NAME,
128  GA_NODE_MACRO_PARAM,
129  GA_NODE_PARAMS,
130  GA_NODE_RESHAPE,
131  GA_NODE_CROSS_PRODUCT,
132  GA_NODE_SWAP_IND,
133  GA_NODE_IND_MOVE_LAST,
134  GA_NODE_CONTRACT,
135  GA_NODE_ALLINDICES,
136  GA_NODE_C_MATRIX,
137  GA_NODE_X,
138  GA_NODE_ELT_SIZE,
139  GA_NODE_ELT_K,
140  GA_NODE_ELT_B,
141  GA_NODE_NORMAL,
142  GA_NODE_VAL,
143  GA_NODE_GRAD,
144  GA_NODE_HESS,
145  GA_NODE_DIVERG,
146  GA_NODE_VAL_TEST,
147  GA_NODE_GRAD_TEST,
148  GA_NODE_HESS_TEST,
149  GA_NODE_DIVERG_TEST,
150  GA_NODE_INTERPOLATE,
151  GA_NODE_INTERPOLATE_FILTER,
152  GA_NODE_INTERPOLATE_VAL,
153  GA_NODE_INTERPOLATE_GRAD,
154  GA_NODE_INTERPOLATE_HESS,
155  GA_NODE_INTERPOLATE_DIVERG,
156  GA_NODE_INTERPOLATE_VAL_TEST,
157  GA_NODE_INTERPOLATE_GRAD_TEST,
158  GA_NODE_INTERPOLATE_HESS_TEST,
159  GA_NODE_INTERPOLATE_DIVERG_TEST,
160  GA_NODE_INTERPOLATE_X,
161  GA_NODE_INTERPOLATE_NORMAL,
162  GA_NODE_INTERPOLATE_DERIVATIVE,
163  GA_NODE_ELEMENTARY,
164  GA_NODE_ELEMENTARY_VAL,
165  GA_NODE_ELEMENTARY_GRAD,
166  GA_NODE_ELEMENTARY_HESS,
167  GA_NODE_ELEMENTARY_DIVERG,
168  GA_NODE_ELEMENTARY_VAL_TEST,
169  GA_NODE_ELEMENTARY_GRAD_TEST,
170  GA_NODE_ELEMENTARY_HESS_TEST,
171  GA_NODE_ELEMENTARY_DIVERG_TEST,
172  GA_NODE_SECONDARY_DOMAIN,
173  GA_NODE_SECONDARY_DOMAIN_VAL,
174  GA_NODE_SECONDARY_DOMAIN_GRAD,
175  GA_NODE_SECONDARY_DOMAIN_HESS,
176  GA_NODE_SECONDARY_DOMAIN_DIVERG,
177  GA_NODE_SECONDARY_DOMAIN_VAL_TEST,
178  GA_NODE_SECONDARY_DOMAIN_GRAD_TEST,
179  GA_NODE_SECONDARY_DOMAIN_HESS_TEST,
180  GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST,
181  GA_NODE_SECONDARY_DOMAIN_X,
182  GA_NODE_SECONDARY_DOMAIN_NORMAL,
183  GA_NODE_XFEM_PLUS,
184  GA_NODE_XFEM_PLUS_VAL,
185  GA_NODE_XFEM_PLUS_GRAD,
186  GA_NODE_XFEM_PLUS_HESS,
187  GA_NODE_XFEM_PLUS_DIVERG,
188  GA_NODE_XFEM_PLUS_VAL_TEST,
189  GA_NODE_XFEM_PLUS_GRAD_TEST,
190  GA_NODE_XFEM_PLUS_HESS_TEST,
191  GA_NODE_XFEM_PLUS_DIVERG_TEST,
192  GA_NODE_XFEM_MINUS,
193  GA_NODE_XFEM_MINUS_VAL,
194  GA_NODE_XFEM_MINUS_GRAD,
195  GA_NODE_XFEM_MINUS_HESS,
196  GA_NODE_XFEM_MINUS_DIVERG,
197  GA_NODE_XFEM_MINUS_VAL_TEST,
198  GA_NODE_XFEM_MINUS_GRAD_TEST,
199  GA_NODE_XFEM_MINUS_HESS_TEST,
200  GA_NODE_XFEM_MINUS_DIVERG_TEST,
201  GA_NODE_ZERO};
202 
203  typedef std::shared_ptr<std::string> pstring;
204  // Print error message indicating the position in the assembly string
205  void ga_throw_error_msg(pstring expr, size_type pos,
206  const std::string &msg);
207 
208 # define ga_throw_error(expr, pos, msg) \
209  { std::stringstream ss; ss << msg; \
210  ga_throw_error_msg(expr, pos, ss.str()); \
211  GMM_ASSERT1(false, "Error in assembly string" ); \
212  }
213 
214  // Structure for the tensor associated with a tree node
215  struct assembly_tensor {
216  bool is_copied;
217  int sparsity_; // 0: plain, 1: vectorized base, 2: vectorised grad, ...
218  size_type qdim_; // Dimension of the vectorization for sparsity tensors
219  base_tensor t;
220  assembly_tensor *tensor_copied;
221 
222  const base_tensor &org_tensor() const
223  { return is_copied ? tensor_copied->org_tensor() : t; }
224  base_tensor &org_tensor()
225  { return is_copied ? tensor_copied->org_tensor() : t; }
226 
227  const base_tensor &tensor() const
228  { return (is_copied ? tensor_copied->tensor() : t); }
229 
230  base_tensor &tensor()
231  { return (is_copied ? tensor_copied->tensor() : t); }
232 
233  void set_sparsity(int sp, size_type q)
234  { sparsity_ = sp; qdim_ = q; }
235 
236  size_type qdim() { return is_copied ? tensor_copied->qdim() : qdim_; }
237 
238  int sparsity() const
239  { return is_copied ? tensor_copied->sparsity() : sparsity_; }
240 
241  inline void set_to_original() { is_copied = false; }
242  inline void set_to_copy(assembly_tensor &t_) {
243  is_copied = true; sparsity_ = t_.sparsity_; qdim_ = t_.qdim_;
244  t = t_.org_tensor(); tensor_copied = &(t_);
245  }
246 
247  inline void adjust_sizes(const bgeot::multi_index &ssizes)
248  { t.adjust_sizes(ssizes); }
249 
250  inline void adjust_sizes()
251  { if (t.sizes().size() || t.size() != 1) t.init(); }
252 
253  inline void adjust_sizes(size_type i)
254  { if (t.sizes().size() != 1 || t.sizes()[0] != i) t.init(i); }
255 
256  inline void adjust_sizes(size_type i, size_type j) {
257  if (t.sizes().size() != 2 || t.sizes()[0] != i || t.sizes()[1] != j)
258  t.init(i, j);
259  }
260 
261  inline void adjust_sizes(size_type i, size_type j, size_type k) {
262  if (t.sizes().size() != 3 || t.sizes()[0] != i || t.sizes()[1] != j
263  || t.sizes()[2] != k)
264  t.init(i, j, k);
265  }
266  inline void adjust_sizes(size_type i, size_type j,
267  size_type k, size_type l) {
268  if (t.sizes().size() != 3 || t.sizes()[0] != i || t.sizes()[1] != j
269  || t.sizes()[2] != k || t.sizes()[3] != l)
270  t.init(i, j, k, l);
271  }
272 
273  void init_scalar_tensor(scalar_type v)
274  { set_to_original(); t.adjust_sizes(); t[0] = v; }
275 
276  void init_vector_tensor(size_type d)
277  { set_to_original(); t.adjust_sizes(d); }
278 
279  void init_matrix_tensor(size_type n, size_type m)
280  { set_to_original(); t.adjust_sizes(n, m); }
281 
282  void init_identity_matrix_tensor(size_type n) {
283  init_matrix_tensor(n, n);
284  auto itw = t.begin();
285  for (size_type i = 0; i < n; ++i)
286  for (size_type j = 0; j < n; ++j)
287  *itw++ = (i == j) ? scalar_type(1) : scalar_type(0);
288  }
289 
290  void init_third_order_tensor(size_type n, size_type m, size_type l)
291  { set_to_original(); t.adjust_sizes(n, m, l); }
292 
293  void init_fourth_order_tensor(size_type n, size_type m,
294  size_type l, size_type k)
295  { set_to_original(); t.adjust_sizes(n, m, l, k); }
296 
297  const bgeot::multi_index &sizes() const { return t.sizes(); }
298 
299  assembly_tensor()
300  : is_copied(false), sparsity_(0), qdim_(0), tensor_copied(0) {}
301  };
302 
303  struct ga_tree_node;
304  typedef ga_tree_node *pga_tree_node;
305 
306 
307  struct ga_tree_node {
308  GA_NODE_TYPE node_type;
309  GA_TOKEN_TYPE op_type;
310  assembly_tensor t;
311  size_type test_function_type; // -1 = undetermined
312  // 0 = no test function,
313  // 1 = first order, 2 = second order,
314  // 3 = both with always first order in first
315  std::string name_test1, name_test2; // variable names corresponding to test
316  // functions when test_function_type > 0.
317  std::string interpolate_name_test1, interpolate_name_test2; // name
318  // of interpolation transformation if any
319  size_type qdim1, qdim2; // Qdims when test_function_type > 0.
320  size_type nbc1, nbc2, nbc3; // For X (nbc1=coordinate number),
321  // macros (nbc1=param number, nbc2,nbc3 type))
322  // and C_MATRIX (nbc1=order).
323  size_type pos; // Position of the first character in string
324  pstring expr; // Original string, for error messages.
325  std::string name; // variable/constant/function/operator name
326  std::string interpolate_name; // For Interpolate : name of transformation
327  std::string interpolate_name_der; // For Interpolate derivative:
328  // name of transformation
329  std::string elementary_name; // For Elementary_transformation :
330  // name of transformation
331  std::string elementary_target;// For Elementary_transformation :
332  // target variable (for its mesh_fem)
333  size_type der1, der2; // For functions and nonlinear operators,
334  // optional derivative or second derivative.
335  bool symmetric_op;
336  pga_tree_node parent; // Parent node
337  std::vector<pga_tree_node> children; // Children nodes
338  scalar_type hash_value; // Hash value to identify nodes.
339  bool marked; // For specific use of some algorithms
340 
341  inline const base_tensor &tensor() const { return t.tensor(); }
342  inline base_tensor &tensor() { return t.tensor(); }
343  int sparsity() const { return t.sparsity(); }
344 
345  inline size_type nb_test_functions() const {
346  if (test_function_type == size_type(-1)) return 0;
347  return test_function_type - (test_function_type >= 2 ? 1 : 0);
348  }
349 
350  inline size_type tensor_order() const
351  { return t.sizes().size() - nb_test_functions(); }
352 
353  inline size_type tensor_test_size() const {
354  size_type st = nb_test_functions();
355  return (st >= 1 ? t.sizes()[0] : 1) * (st == 2 ? t.sizes()[1] : 1);
356  }
357 
358  inline size_type tensor_proper_size() const
359  { return t.org_tensor().size() / tensor_test_size(); }
360 
361  inline size_type tensor_proper_size(size_type i) const
362  { return t.sizes()[nb_test_functions()+i]; }
363 
364 
365  void mult_test(const pga_tree_node n0, const pga_tree_node n1);
366 
367  bool tensor_is_zero() {
368  if (node_type == GA_NODE_ZERO) return true;
369  if (node_type != GA_NODE_CONSTANT) return false;
370  for (size_type i = 0; i < tensor().size(); ++i)
371  if (tensor()[i] != scalar_type(0)) return false;
372  return true;
373  }
374 
375  inline bool is_constant() {
376  return (node_type == GA_NODE_CONSTANT ||
377  (node_type == GA_NODE_ZERO && test_function_type == 0));
378  }
379 
380  inline void init_scalar_tensor(scalar_type v)
381  { t.init_scalar_tensor(v); test_function_type = 0; }
382 
383  inline void init_vector_tensor(size_type d)
384  { t.init_vector_tensor(d); test_function_type = 0; }
385 
386  inline void init_matrix_tensor(size_type n, size_type m)
387  { t.init_matrix_tensor(n, m); test_function_type = 0; }
388 
389  inline void init_identity_matrix_tensor(size_type n)
390  { t.init_identity_matrix_tensor(n); test_function_type = 0; }
391 
392  inline void init_third_order_tensor(size_type n, size_type m, size_type l)
393  { t.init_third_order_tensor(n, m, l); test_function_type = 0; }
394 
395  inline void init_fourth_order_tensor(size_type n, size_type m,
396  size_type l, size_type k)
397  { t.init_fourth_order_tensor(n, m, l, k); test_function_type = 0; }
398 
399  inline void adopt_child(pga_tree_node new_child)
400  { children.push_back(new_child); children.back()->parent = this; }
401 
402  inline void replace_child(pga_tree_node oldchild,
403  pga_tree_node newchild) {
404  bool found = false;
405  for (pga_tree_node &child : children)
406  if (child == oldchild) { child = newchild; found = true; }
407  GMM_ASSERT1(found, "Internal error");
408  }
409 
410  ga_tree_node()
411  : node_type(GA_NODE_VOID), test_function_type(-1), qdim1(0), qdim2(0),
412  nbc1(0), nbc2(0), nbc3(0), pos(0), expr(0), der1(0), der2(0),
413  symmetric_op(false), hash_value(0) {}
414  ga_tree_node(GA_NODE_TYPE ty, size_type p, pstring expr_)
415  : node_type(ty), test_function_type(-1),
416  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
417  pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false),
418  hash_value(0) {}
419  ga_tree_node(scalar_type v, size_type p, pstring expr_)
420  : node_type(GA_NODE_CONSTANT), test_function_type(-1),
421  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
422  pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false),
423  hash_value(0)
424  { init_scalar_tensor(v); }
425  ga_tree_node(const char *n, size_type l, size_type p, pstring expr_)
426  : node_type(GA_NODE_NAME), test_function_type(-1),
427  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
428  pos(p), expr(expr_), name(n, l), der1(0), der2(0), symmetric_op(false),
429  hash_value(0) {}
430  ga_tree_node(GA_TOKEN_TYPE op, size_type p, pstring expr_)
431  : node_type(GA_NODE_OP), op_type(op), test_function_type(-1),
432  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
433  pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false),
434  hash_value(0) {}
435  };
436 
437  struct ga_tree {
438  pga_tree_node root, current_node;
439  std::string secondary_domain;
440 
441  void add_scalar(scalar_type val, size_type pos, pstring expr);
442  void add_allindices(size_type pos, pstring expr);
443  void add_name(const char *name, size_type length, size_type pos,
444  pstring expr);
445  void add_sub_tree(ga_tree &sub_tree);
446  void add_params(size_type pos, pstring expr);
447  void add_matrix(size_type pos, pstring expr);
448  void add_op(GA_TOKEN_TYPE op_type, size_type pos, pstring expr);
449  void clear_node_rec(pga_tree_node pnode);
450  void clear_node(pga_tree_node pnode);
451  void clear() { clear_node_rec(root); root = current_node = nullptr; }
452  void clear_children(pga_tree_node pnode);
453  void replace_node_by_child(pga_tree_node pnode, size_type i);
454  void copy_node(pga_tree_node pnode, pga_tree_node parent,
455  pga_tree_node &child);
456  void duplicate_with_operation(pga_tree_node pnode, GA_TOKEN_TYPE op_type);
457  void duplicate_with_addition(pga_tree_node pnode)
458  { duplicate_with_operation(pnode, GA_PLUS); }
459  void duplicate_with_substraction(pga_tree_node pnode)
460  { duplicate_with_operation(pnode, GA_MINUS); }
461  void insert_node(pga_tree_node pnode, GA_NODE_TYPE node_type);
462  void add_child(pga_tree_node pnode, GA_NODE_TYPE node_type = GA_NODE_VOID);
463  void swap(ga_tree &tree) {
464  std::swap(root, tree.root);
465  std::swap(current_node, tree.current_node);
466  std::swap(secondary_domain, tree.secondary_domain);
467  }
468 
469  ga_tree() : root(nullptr), current_node(nullptr), secondary_domain() {}
470 
471  ga_tree(const ga_tree &tree) : root(nullptr), current_node(nullptr),
472  secondary_domain(tree.secondary_domain)
473  { if (tree.root) copy_node(tree.root, nullptr, root); }
474 
475  ga_tree &operator =(const ga_tree &tree) {
476  clear(); secondary_domain = tree.secondary_domain;
477  if (tree.root)
478  copy_node(tree.root,nullptr,root);
479  return *this;
480  }
481 
482  ~ga_tree() { clear(); }
483  };
484 
485  // Test equality or equivalence of two sub trees.
486  // version = 0 : strict equality
487  // 1 : give the same result
488  // 2 : give the same result with transposition of test functions
489  bool sub_tree_are_equal
490  (const pga_tree_node pnode1, const pga_tree_node pnode2,
491  const ga_workspace &workspace, int version);
492 
493  // Transform the expression of a node and its sub-nodes in the equivalent
494  // assembly string sent to ostream str
495  void ga_print_node(const pga_tree_node pnode,
496  std::ostream &str);
497  // The same for the whole tree, the result is a std::string
498  std::string ga_tree_to_string(const ga_tree &tree);
499 
500  // Syntax analysis of an assembly string. Conversion to a tree.
501  // No semantic analysis is done. The tree can be inconsistent.
502  void ga_read_string(const std::string &expr, ga_tree &tree,
503  const ga_macro_dictionary &macro_dict);
504  void ga_read_string_reg(const std::string &expr, ga_tree &tree,
505  ga_macro_dictionary &macro_dict);
506 
507 
508 } /* end of namespace */
509 
510 
511 #endif /* GETFEM_GENERIC_ASSEMBLY_TREE_H__ */
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
dal_singleton.h
A simple singleton implementation.
getfem_copyable_ptr.h
A smart pointer that copies the value it points to on copy operations.
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
bgeot_rtree.h
region-tree for window/point search on a set of rectangles.
getfem_models.h
Model representation in Getfem.
getfem_omp.h
Tools for multithreaded, OpenMP and Boost based parallelization.
getfem_generic_assembly.h
A language for generic assembly of pde boundary value problems.
bgeot_geotrans_inv.h
Inversion of geometric transformations.
gmm_blas.h
Basic linear algebra functions.